Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Byte your tongue : a computational model of human mandibular-lingual biomechanics for biomedical applications 2010

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2011_spring_stavness_ian.pdf
ubc_2011_spring_stavness_ian.pdf
ubc_2011_spring_stavness_ian.pdf [ 26.33MB ]
Metadata
JSON: 1.0071518.json
JSON-LD: 1.0071518+ld.json
RDF/XML (Pretty): 1.0071518.xml
RDF/JSON: 1.0071518+rdf.json
Turtle: 1.0071518+rdf-turtle.txt
N-Triples: 1.0071518+rdf-ntriples.txt
Citation
1.0071518.ris

Full Text

Byte Your Tongue A Computational Model of Human Mandibular-Lingual Biomechanics for Biomedical Applications by Ian Kent Stavness B.Eng., University of Saskatchewan, 2004 B.Sc., University of Saskatchewan, 2004 M.A.Sc., University of British Columbia, 2006 A DISSERTATION SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2010 c© Ian Kent Stavness, 2010 Abstract Biomechanical models provide a means to analyze movement and forces in highly complex anatomical systems. Models can be used to explain cause and effect in normal body function as well as in abnormal cases where under- lying causes of dysfunction can be clarified. In addition, computer models can be used to simulate surgical changes to bone and muscle structure al- lowing for prediction of functional and aesthetic outcomes. This dissertation proposes a state-of-the-art model of coupled jaw-tongue- hyoid biomechanics for simulating combined jaw and tongue motor tasks, such as chewing, swallowing, and speaking. Simulation results demonstrate that mechanical coupling of tongue muscles acting on the jaw and jaw mus- cles acting on the tongue are significant and should be considered in orofacial modeling studies. Towards validation of the model, simulated tongue veloc- ity and tongue-palate pressure are consistent with published measurements. Inverse simulation methods are also discussed along with the implemen- tation of a technique to automatically compute muscle activations for track- ing a target kinematic trajectory for coupled skeletal and soft-tissue models. Additional target parameters, such as dynamic constraint forces and stiff- ness, are included in the inverse formulation to control muscle activation predictions in redundant models. Simulation results for moving and de- forming muscular-hydrostat models are consistent with published theoreti- cal proposals. Also, muscle activations predicted for lateral jaw movement are consistent with published literature on jaw physiology. As an illustrative case study, models of segmental jaw surgery with and without reconstruction are developed. The models are used to simulate clin- ically observed functional deficits in movement and bite force production. The inverse simulation tools are used to predict muscle forces that could theoretically be used by a patient to compensate for functional deficits fol- lowing jaw surgery. The modeling tools developed and demonstrated in this dissertation provide a foundation for future studies of orofacial function and biomedical applications in oral and maxillofacial surgery and treatment. ii Preface Parts of this dissertation have been published elsewhere. Chapter 3 and Appendix C have been published in Stavness, Lloyd, Payan, and Fels [186] ( c© Wiley (2010), reproduced with permission). For Chapter 3, I wrote the source code for the model, performed the simulations and analysis, and wrote the text. Dr. Payan assisted with analysis and writing. Dr. Fels provided editorial feedback on the manuscript. Dr. Lloyd wrote the material that I have reproduced in Appendix C as background on the simulation techniques used in this dissertation. Versions of Section 4.1, Section 5.1, Section 5.3 have been published in Stavness, Hannam, Lloyd, and Fels [185] ( c© Taylor & Francis (2010), reproduced with permission). I formulated and implemented the inverse solver in consultation with Dr. Lloyd. I developed the jaw surgery models with Dr. Hannam. I performed the inverse simulations and analysis, and wrote the manuscript. Dr. Hannam and Dr. Fels provided feedback on the manuscript. Versions of Section 5.1 and Section 5.2 have been published in Hannam, Stavness, Lloyd, Fels, Miller, and Curtis [72] ( c© Elsevier (2010), reproduced with permission). Dr. Hannam performed the deficit simulations and anal- ysis with my assistance. Dr. Hannam wrote the text for the manuscript that I have adapted and included in Section 5.2. Dr. Fels, Dr. Curtis, and Dr. Miller provided editorial feedback on the manuscript. Section 3.4 describes ongoing collaborative work. I am developing the face model jointly with Dr. Payan and the soft-palate model jointly with Dr. Hui Chen. The integration of hyolaryngeal data and model has been published online in Stavness, Ludlow, Chung, and Fels [184]. Data was collected by Dr. Ludlow and Dr. Chung. I developed the model, formulated the data-model integration plan, and wrote the text for the manuscript under the guidance of Dr. Fels and Dr. Ludlow. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3 Dissertation Outline . . . . . . . . . . . . . . . . . . . . . . . 10 2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1 Functional Data Recording . . . . . . . . . . . . . . . . . . . 13 2.1.1 Movement . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.2 Forces . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.3 Muscle activity . . . . . . . . . . . . . . . . . . . . . 19 2.2 Structural Data Measurement . . . . . . . . . . . . . . . . . 20 2.2.1 Measuring anatomical structure . . . . . . . . . . . . 20 2.2.2 Imaging anatomical structure . . . . . . . . . . . . . 21 2.2.3 Measuring mechanical tissue properties . . . . . . . . 24 2.3 Biomechanical Modeling . . . . . . . . . . . . . . . . . . . . 25 2.3.1 Rigid bone modeling . . . . . . . . . . . . . . . . . . 25 2.3.2 Deformable tissue modeling . . . . . . . . . . . . . . 27 2.3.3 Muscle modeling . . . . . . . . . . . . . . . . . . . . . 29 2.3.4 Isolated orofacial models . . . . . . . . . . . . . . . . 33 iv Table of Contents 2.3.5 Coupled orofacial models . . . . . . . . . . . . . . . . 34 2.4 Inverse Methods . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.5 Biomedical Applications . . . . . . . . . . . . . . . . . . . . . 39 2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3 Jaw-Tongue-Hyoid Model . . . . . . . . . . . . . . . . . . . . 43 3.1 Model Creation . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.1.1 Jaw-hyoid model . . . . . . . . . . . . . . . . . . . . . 45 3.1.2 Tongue model . . . . . . . . . . . . . . . . . . . . . . 46 3.1.3 Registration . . . . . . . . . . . . . . . . . . . . . . . 48 3.1.4 Attachment . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Simulation Descriptions . . . . . . . . . . . . . . . . . . . . . 50 3.2.1 Jaw activated tasks . . . . . . . . . . . . . . . . . . . 51 3.2.2 Tongue activated tasks . . . . . . . . . . . . . . . . . 53 3.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3.1 Jaw-tongue-hyoid coupling . . . . . . . . . . . . . . . 53 3.3.2 Comparison with published data . . . . . . . . . . . . 57 3.3.3 Integration Error . . . . . . . . . . . . . . . . . . . . 60 3.4 Additional Orofacial Sub-Models . . . . . . . . . . . . . . . . 61 3.4.1 Face model . . . . . . . . . . . . . . . . . . . . . . . . 61 3.4.2 Soft-palate model . . . . . . . . . . . . . . . . . . . . 63 3.4.3 Hyoid-larynx model . . . . . . . . . . . . . . . . . . . 63 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.6 Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4 Inverse Simulation Methods . . . . . . . . . . . . . . . . . . . 69 4.1 Inverse Solver Formulation . . . . . . . . . . . . . . . . . . . 71 4.2 Analysis with Canonical Models . . . . . . . . . . . . . . . . 77 4.2.1 Point inverse . . . . . . . . . . . . . . . . . . . . . . . 77 4.2.2 Rigid-body inverse . . . . . . . . . . . . . . . . . . . . 81 4.2.3 Deformable-body inverse . . . . . . . . . . . . . . . . 81 4.3 Analysis with Anatomical Models . . . . . . . . . . . . . . . 84 4.3.1 Jaw inverse . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3.2 Tongue inverse . . . . . . . . . . . . . . . . . . . . . . 88 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.5 Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 v Table of Contents 5 Segmental Jaw Surgery Models . . . . . . . . . . . . . . . . . 95 5.1 Model Creation . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.2 Forward Dynamics: Deficit Simulations . . . . . . . . . . . . 99 5.2.1 Simulation descriptions . . . . . . . . . . . . . . . . . 101 5.2.2 Simulation results . . . . . . . . . . . . . . . . . . . . 101 5.2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3 Inverse Dynamics: Compensatory Simulations . . . . . . . . 107 5.3.1 Simulation descriptions . . . . . . . . . . . . . . . . . 107 5.3.2 Simulation results . . . . . . . . . . . . . . . . . . . . 108 5.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 112 5.4 Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.1 Dissertation Contributions . . . . . . . . . . . . . . . . . . . 118 6.2 Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.3 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . 124 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Appendices A List of Publications . . . . . . . . . . . . . . . . . . . . . . . . . 148 A.1 Journal Publications . . . . . . . . . . . . . . . . . . . . . . . 148 A.2 Conference Publications . . . . . . . . . . . . . . . . . . . . . 149 A.3 Research Visits . . . . . . . . . . . . . . . . . . . . . . . . . . 149 A.4 Research Talks . . . . . . . . . . . . . . . . . . . . . . . . . . 150 A.5 Master’s Publications . . . . . . . . . . . . . . . . . . . . . . 150 A.6 Additional Publications . . . . . . . . . . . . . . . . . . . . . 151 B Head and Neck Anatomy . . . . . . . . . . . . . . . . . . . . . 152 B.1 Bone Structures . . . . . . . . . . . . . . . . . . . . . . . . . 152 B.1.1 Cranium . . . . . . . . . . . . . . . . . . . . . . . . . 152 B.1.2 Mandible . . . . . . . . . . . . . . . . . . . . . . . . . 154 B.1.3 Dentition . . . . . . . . . . . . . . . . . . . . . . . . . 155 B.1.4 Hyoid bone . . . . . . . . . . . . . . . . . . . . . . . . 156 B.1.5 Vertebrae . . . . . . . . . . . . . . . . . . . . . . . . . 156 B.2 Soft-Tissues and Muscles . . . . . . . . . . . . . . . . . . . . 156 B.2.1 Face . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 vi Table of Contents B.2.2 Jaw muscles . . . . . . . . . . . . . . . . . . . . . . . 158 B.2.3 Tongue . . . . . . . . . . . . . . . . . . . . . . . . . . 162 B.2.4 Soft-palate . . . . . . . . . . . . . . . . . . . . . . . . 164 B.2.5 Pharynx . . . . . . . . . . . . . . . . . . . . . . . . . 165 B.2.6 Larynx . . . . . . . . . . . . . . . . . . . . . . . . . . 166 B.2.7 Epiglottis . . . . . . . . . . . . . . . . . . . . . . . . . 169 B.2.8 Temporomandibular joint . . . . . . . . . . . . . . . . 169 C ArtiSynth Simulation Software . . . . . . . . . . . . . . . . . 172 C.1 Physical Simulation Framework . . . . . . . . . . . . . . . . 173 C.1.1 Friction, damping, and stabilization . . . . . . . . . . 175 C.1.2 System solution and complexity . . . . . . . . . . . . 177 C.1.3 Attachments between bodies . . . . . . . . . . . . . . 179 C.1.4 Contact handling . . . . . . . . . . . . . . . . . . . . 181 C.1.5 Simulation engine summary . . . . . . . . . . . . . . 184 C.1.6 Validation using ANSYS . . . . . . . . . . . . . . . . 185 C.2 Graphical Toolset . . . . . . . . . . . . . . . . . . . . . . . . 187 C.2.1 Model component hierarchy . . . . . . . . . . . . . . 188 C.2.2 Viewing, selection, and editing . . . . . . . . . . . . . 188 C.2.3 Properties, control panels, and probes . . . . . . . . . 189 vii List of Tables 3.1 Physiological cross-sectional area of jaw and tongue muscles . 46 3.2 Muscle activation amplitudes for jaw-tongue-hyoid simulations 50 5.1 REST to OPEN to CLOSE inverse simulation results . . . . . 112 5.2 Stable unilateral clenching inverse simulation results . . . . . 114 viii List of Figures 1.1 Upper airway anatomy, CT image and illustration . . . . . . 2 1.2 Jaw-tongue-hyoid model in ArtiSynth . . . . . . . . . . . . . 4 2.1 Early photographic recordings of horse movement by Muybridge 13 2.2 Computed tomography images of the head . . . . . . . . . . . 22 2.3 Magnetic resonance images of the head . . . . . . . . . . . . . 23 2.4 Force characteristics of the Hill-type muscle model . . . . . . 32 3.1 Jaw-hyoid model with labeled muscles . . . . . . . . . . . . . 45 3.2 Tongue model with labeled muscles . . . . . . . . . . . . . . . 47 3.3 Jaw-tongue-hyoid model . . . . . . . . . . . . . . . . . . . . . 49 3.4 Muscle activation patterns for jaw-tongue-hyoid simulations . 51 3.5 Jaw opening simulation results . . . . . . . . . . . . . . . . . 54 3.6 Jaw activated tasks simulation results . . . . . . . . . . . . . 55 3.7 Tongue retraction simulation results . . . . . . . . . . . . . . 56 3.8 Tongue-palate contact simulation results . . . . . . . . . . . . 57 3.9 Tongue velocity during simulated and recorded retraction . . 58 3.10 Assessment of integration error . . . . . . . . . . . . . . . . . 60 3.11 Face model integrated with jaw-tongue-hyoid model . . . . . 61 3.12 Soft-palate model registered to jaw-tongue-hyoid model . . . 62 3.13 Hyoid-larynx model and video fluoroscopy landmarks . . . . . 64 4.1 Point-model inverse simulation: kinematic target . . . . . . . 78 4.2 Point-model inverse simulation: output muscle activations . . 78 4.3 Point-model inverse simulation: output muscle activations . . 79 4.4 Rigid-body model inverse simulation . . . . . . . . . . . . . . 80 4.5 Deformable beam inverse simulation: protrusion . . . . . . . 82 4.6 Deformable beam inverse simulation: output muscle activations 82 4.7 Deformable beam inverse simulation: complex movement . . . 83 4.8 Inverse jaw simulation: lateral incisor movement . . . . . . . 85 4.9 Inverse jaw laterotrusion simulation: output muscle activations 86 4.10 Tongue inverse simulation: tip elevation . . . . . . . . . . . . 89 ix List of Figures 4.11 Tongue tip elevation inverse simulation: output muscle acti- vations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.1 Segmental jaw resection and reconstruction models . . . . . . 97 5.2 Segmental jaw resection model with labeled muscles. . . . . . 98 5.3 Segmental jaw resection model conventions and restraints . . 100 5.4 Incisor displacements during FORCE and OPEN simulations 103 5.5 NOCON model during individual closer muscle activations . . 105 5.6 Target jaw postures to move from REST to OPEN to CLOSE 109 5.7 REST to OPEN inverse simulation results . . . . . . . . . . . 110 5.8 OPEN to CLOSE inverse simulation results . . . . . . . . . . 111 5.9 Target jaw postures for unilateral clenching tasks . . . . . . . 113 B.1 Upper airway anatomy, sagittal cross-section view . . . . . . 153 B.2 Skull, lateral view . . . . . . . . . . . . . . . . . . . . . . . . 153 B.3 Mandible, medial and lateral views . . . . . . . . . . . . . . . 154 B.4 Tooth and alveolar process, lateral view . . . . . . . . . . . . 155 B.5 Lower and upper dentition . . . . . . . . . . . . . . . . . . . . 155 B.6 Facial muscles, lateral view . . . . . . . . . . . . . . . . . . . 157 B.7 Jaw closing muscles . . . . . . . . . . . . . . . . . . . . . . . . 158 B.8 Jaw opening muscles . . . . . . . . . . . . . . . . . . . . . . . 159 B.9 Tongue muscles, lateral view . . . . . . . . . . . . . . . . . . 161 B.10 Tongue muscles, posterior view . . . . . . . . . . . . . . . . . 161 B.11 Soft-palate muscles, posterior view . . . . . . . . . . . . . . . 164 B.12 Larynx and pharynx . . . . . . . . . . . . . . . . . . . . . . . 165 B.13 Oral cavity during swallowing . . . . . . . . . . . . . . . . . . 166 B.14 Neck muscles, frontal view . . . . . . . . . . . . . . . . . . . . 167 B.15 Temporomandibular joint . . . . . . . . . . . . . . . . . . . . 170 B.16 Temporomandibular joint ligaments . . . . . . . . . . . . . . 170 C.1 Errors due to decoupling constraints from the velocity solve . 177 C.2 Factor times for A as a function of system size . . . . . . . . 179 C.3 Contact handling between two deformable models . . . . . . . 182 C.4 ArtiSynth and ANSYS static comparisons . . . . . . . . . . . 185 C.5 ArtiSynth and ANSYS time integration comparisons . . . . . 187 C.6 ArtiSynth model viewer and timeline controller . . . . . . . . 188 C.7 ArtiSynth model hierarchy view and control panel . . . . . . 189 x Glossary 1D One Dimensional 2D Two Dimensional 3D Three Dimensional DOF Degrees-of-Freedom EMG Electromyography EPG Electropalatography EM Electromagnetic EMA Electromagnetic Articulometry VF Video-fluoroscopy US Ultrasound CT Computed Tomography CBCT Cone-Beam Computed Tomography MRI Magnetic Resonance Imaging VHD Visible Human Dataset EPH Equilibrium Point Hypothesis FEM Finite-Element Method OSA Obstructive Sleep Apnea PCA Principal Components Analysis CSA Cross-Sectional Area TMJ Temporomandibular Joint xi Acknowledgements My rich and positive experience in graduate studies is due in large part to my supervisor Sid Fels. Many thanks for his guidance, mentorship, and friendship over the past six years. I have had the pleasure of working closely with a number of fine researchers during this time as well. Thanks to John Lloyd for many fruitful sessions at the white board and at the pub. Thanks to Alan Hannam for many enjoyable discussions and for providing me with a clear model of excellent academic scholarship. Also, thanks to Yohan Payan for a productive collaboration over the past year. Thanks must go to Christy Ludlow for hosting my visit to her lab at the National Institutes of Health as well as to Chris Peck and Alex Wirianski for hosting my visit to their research group at the University of Sydney. Acknowledgement is due to researchers at Gipsa-Lab, Grenoble for kindly providing data to support this project, including Pierre Badin for providing the CT image data as well as Pascal Perrier, Stephanie Buchaillard, and Mohammad Ali Nazari for providing the tongue and face model geometry. Finally, a heart felt thank-you to my family and friends for their love and support during this phase of my life. Special thanks to my parents for providing me the strong foundation upon which my accomplishments are grounded. My deepest thanks go to my dear, Sara, thank you for your daily encouragement, much laughter, and endless love. xii Chapter 1 Introduction Computer simulation is becoming an integral aspect of biomedical research and practice, in applications ranging from basic research in physiology to clinical applications in treatment planning and surgical training. Medical imaging technology has revolutionized our ability to observe internal struc- tures of the body and is one of the most significant biomedical advances of the 20th century. Computer simulation aims to build upon medical imag- ing to further our understanding of how the body functions. Examples of biological simulation include visualizing dynamic movements from static im- ages, predicting unobservable variables such as stresses within tissue, and illuminating mechanisms of sensorimotor control. An important example of biomedical computer simulation is the dynamic modeling of muscle-driven anatomical structures as a means to better understand their function in normal and pathological cases. Modeling such structures is non-trivial, and often involves the combined simulation of hard structures (such as bones), interconnected by constraints (such as joints), attached to various soft tis- sues (such as skin, fat, mucosa, muscle and tendon), and in contact with each other and the environment. This dissertation proposes computational methods to analyze human orofacial biomechanics in order to better under- stand structure-function relationships and motor control strategies in oral movement, including mastication and speech production, towards develop- ing new tools for computer-assisted oral and maxillofacial surgery. 1.1 Motivation The research described in this dissertation is motivated by a desire to create computational tools for analyzing oral and facial biomechanics in the context 1 1.1. Motivation TongueLips Mandible Hyoid bone Epiglottis Trachea Thyroid cartilage Cricoid cartilage Vertebrae Soft palate Hard palate Vocal folds TongueLips Larynx Figure 1.1: Mid-sagittal CT image and illustration of the upper airway anatomy. CT data courtesy of Dr. Pierre Badin, Gipsa-Lab, Grenoble. Illustration c© Elsevier (2010), Drake et al. [47], adapted with permission. of biomedical applications, such as improving our understanding of upper- airway dysfunction and related medical and surgical treatments. Human upper-airway anatomy is structurally complex and critical to a number of life-sustaining functions, making it a good candidate for analysis through biomechanics simulation. Biomedical applications Computational biomechanics has a wide range of viable applications in medicine. Diagnosis of motor system disorders can be improved with biome- chanics simulation to augment and enhance medical images and other ob- servational data. Further, biomechanics simulation can help to determine cause-and-effect in order to establish a deeper understanding of motor sys- tem dysfunction. Biomechanics simulation can also aid in planning of treat- ment and surgical interventions. Computer models can be used to evaluate alternative treatment paths or tailor a particular treatment to a specific patient. Also, enhancing post-operative evaluation with biomechanics sim- 2 1.1. Motivation ulation can help to guide patient rehabilitation. In Chapter 5, we investi- gate segmental jaw resection and reconstruction as a case study of applying biomechanics analysis to post-operative analysis and prediction. Diagnosis Comprehensive dynamic models of the orofacial region will en- hance our understanding of both its normal physiological function and its dysfunctions, such as Obstructive Sleep Apnea (OSA) [86], swallowing dis- orders [110], and speech pathologies. Simulations can predict immeasurable biomechanical quantities that may correlate with dysfunction and there- fore could be used to enhance diagnosis protocols. For example, forces in the Temporomandibular Joint (TMJ) are difficult to measure, but can be determined in a biomechanical analysis of jaw movement and used in the diagnosis of TMJ disorder. Treatment planning Simulating a variety of potential treatments, such as different surgical procedures or prostheses, could inform a treatment plan and complement other factors such as clinician experience and intuition. For example, in surgery planning simulation can be used to predict the consequences and deficits associated with a particular surgical alteration to a musculoskeletal system. Given the mechanical complexity of Three Dimensional (3D) musculoskeletal structures, the small sample size of pa- tients, and the significant risks involved, surgical innovation is necessarily conservative. Computer simulation allows for iterative refinement of surgical procedures with little cost and risk to patients. Section 5.2 reports simula- tions comparing theoretical post-operative deficits for jaw surgery with and without reconstruction. This type of analysis could be used on a patient- specific bases to determine, given the planned extent of tissue resection, whether or not a jaw reconstruction would be beneficial. Post-operative rehabilitation Biomechanics simulation of surgically re- constructed anatomy can also provide post-operative benefit by guiding re- habilitation. Given a model of a specific patient’s reconstruction, simulation of different muscle activation patterns may illuminate new motor strategies 3 1.1. Motivation Figure 1.2: Screen shot of our ArtiSynth modeling toolkit showing the jaw- tongue-hyoid model along with a timeline for controlling input/output data and a control panel for adjusting model properties. to compensate for the altered musculoskeletal structure. Novel motor strate- gies could potentially guide post-operative therapy decisions, for example, which muscles are important to strengthen, or conversely which muscles forces need to be reduced via an intervention such as botulinum toxin in- jection. Section 5.3 reports inverse simulations predicting jaw muscle forces required to compensate for simulated deficits consequent to jaw resection. Importance of orofacial function and dysfunction The human orofacial anatomy, pictured in Figure 1.1, is involved in a num- ber of life sustaining functions including chewing, swallowing, and breathing. Therefore, modeling orofacial biomechanics for understanding and treating orofacial dysfunction has a large potential benefit. Dysphagia, or swallow- ing disorders, is a serious concern for stroke patients as aspiration-related pneumonia leads to 40,000 deaths per year in North America [110]. OSA is a serious respiratory condition involving tissue collapse in the upper airway 4 1.1. Motivation during sleep that afflicts 12 million people in United States. In addition to eating and breathing, the orofacial anatomy is central to speech production. The analysis of speech production mechanisms has the potential to benefit the treatment of speech pathology. Head and neck cancer or trauma is commonly treated with surgical re- section and reconstruction of orofacial tissues, including the jaw and tongue. Head and neck cancer is the sixth most common cancer worldwide and ac- counts for 500,000 new cases per year [24]. Post-surgical deficits can include dysphagia, chewing disorders, OSA [141], speech deficits, as well as signif- icant alteration of facial aesthetics. These surgical procedures are highly complex, highly invasive, and can severely impact patient quality of life. Traditional approaches to observing and recording orofacial function are reviewed in Section 2.1 and we will argue that the inaccessibility of many functionally significant parameters, as well as the structural complexity of the orofacial anatomy, present significant barriers for traditional data collec- tion methods to provide sufficient insight into normal and abnormal orofacial function. Advanced biomechanical toolsets, such as the jaw-tongue-hyoid model developed in Section 3.1 and pictured in Figure 1.2, are required to make new progress in this area. Complexity of orofacial anatomy The high complexity of orofacial anatomy, as shown in Figure 1.1, is an- other motivation for developing new computational techniques to investigate orofacial function; anatomical complexity complicates empirically derived functional hypotheses and prevents complete characterization with empirical techniques. The orofacial anatomy is composed of rigid structures includ- ing the cranium, jaw, and hyoid bone, highly deformable muscle activated tissues such as the tongue, soft palate, and pharynx, and larynx, an intri- cate arrangement of many muscles some of which are capable of exerting very large forces (up to 200 N during tooth clenching), and various con- tact and constraint situations including bite contact and the TMJ. Given the complexity of the upper airway, it is not surprising that its function 5 1.1. Motivation and motor control strategies are not completely understood, especially for complex motor actions such as speech production that involve the coordi- nation of multiple structures in very fast movements up to 20 cm/s [144]. Appendix B provides background information on the orofacial anatomy and muscles. Section 2.2 reviews techniques employed to measure and image anatomical structures and to measure and estimate the mechanical proper- ties of biological tissues. Need for mechanics While movement can be directly recorded with a variety of techniques, other indirect biomechanical variables, such as forces, are also important and use- ful. For example, force predictions are important for understanding tooth loading when designing dental prosthesis to ensure that the prosthesis is suf- ficiently durable, given the expected loading during chewing and clenching. Tooth forces are also important to ensure that reshaping the occlusal sur- faces do not generate abnormal or unbalanced TMJ forces that could lead to articular disc dysfunction. Analyzing force information given kinematics and reaction force record- ings requires information about the mechanics of the system under analysis. As mentioned above, the complexity of human body mechanics is significant, even as compared with modern robotic systems, and in particular biome- chanical systems exhibit unique features: compliant, non-linear actuators as well as kinematic and motor redundancy. Additionally, human movements involve the dynamics of the body and therefore dynamics are important as opposed to simply forces at a quasi-static equilibrium. Section 2.3 reviews previously proposed mechanical models of biological systems, focusing pri- marily on previous models of the human orofacial subsystems, including the jaw, tongue, larynx, and face. Need for dynamics Dynamics are important in speech production, which involves very fast movement of the vocal tract articulators (up to 80 cm/s [149]). The dy- 6 1.1. Motivation namic trajectory of the vocal tract articulations, not just static postures, are important to the acoustics of speech. Also, inertial effects are non- negligible for fast orofacial movements, especially for the tongue. Therefore, we have focused on dynamic simulation techniques, as opposed to quasi- static techniques, both for creating forward dynamic simulations of orofacial movements and in the development of inverse techniques to calculate muscle activations to drive fully dynamic models. Appendix C provides background material on the forward dynamics formulation in ArtiSynth, which was used to create models of the normal jaw-tongue-hyoid (Section 3.1) and abnormal jaws (Section 5.1). Need for holistic models The anatomical structures of the upper airway and face are interconnected. Therefore isolating a subsystem with boundary conditions in a model may be insufficient to fully reproduce orofacial function. This is especially true for the jaw-tongue-hyoid system, where the position of the tongue within the oral cavity is largely determined by the positioning of the jaw and hy- oid, on which the tongue is said to “ride”[78]. Jaw-tongue coupling effects may also play a role in co-articulation effects in speech production [128]. Chapter 5 describes the development and analysis of the first-of-its-kind dy- namically coupled jaw-tongue-hyoid model using dynamic Finite-Element Method (FEM) combined with rigid-body dynamics. Further, an inves- tigation of incorporating the face, soft-palate, and larynx are detailed in Section 3.4. Need for inverse methods While the capability to create representations of musculoskeletal systems and compute the dynamics of the resulting coupled mechanical system has advanced markedly, generating useful, plausible, and accurate motion sim- ulations is still a significant challenge. Generating motion simulations for biomechanical models involves computing an input signal to “drive” a model to perform a desired motor task. Tasks are typically described in terms of 7 1.2. Contributions kinematics (motion) while the drive is thought of as time-varying excitation of motor units or muscles. Controlling dynamics simulations of multi-muscle anatomical systems is challenging due a high-dimensional redundant control space. Manual trial-and-error tuning of muscle activation inputs to a forward dynamics simulation is intractable and difficult to evaluate with respect to recorded subject data. Inverse methods will increase the utility of biome- chanics modeling by systematically predicting activations to achieve target outputs. Further, understanding the motor control strategies underlying human feeding and speech production is an important and active area of research. In speech production, the traditional approach uses statistical analysis of jaw and tongue kinematics; however, this approach cannot determine which aspects of the observed movements arise from central motor commands and which are due to the mechanics of the peripheral musculoskeletal system. Chapter 4 describes the inverse-dynamics techniques developed within Ar- tiSynth for automatic muscle activations prediction for trajectory tracking with the dynamic models. 1.2 Contributions The contributions of this dissertation include creating state-of-the-art mod- els for orofacial biomechanics (Chapter 3), developing new tools for inverse- dynamics simulation (Chapter 4), and applying the models and tools to the analysis of segmental jaw surgery as an example biomedical application (Chapter 5). Publications resulting from this work are listed in Appendix A and the main contributions are summarized here. Modeling of coupled jaw-tongue-hyoid biomechanics i. Created a novel model of the jaw-tongue-hyoid system. We developed a state-of-the-art model of coupled jaw-tongue-hyoid biome- chanics in order to analyze combined jaw and tongue motor tasks, such as chewing, swallowing, and speaking. 8 1.2. Contributions ii. Demonstrated the significance of coupling. We demonstrated that the mechanical coupling of tongue muscles acting on the jaw, and vice versa, are significant and are an important factor for consideration in future orofacial modeling studies. iii. Compared simulations with recorded tongue velocity and pres- sure. We found that simulations of tongue velocity in speech and max- imum voluntary tongue-palate pressure compared well with published measurements. iv. Used as a test case for the ArtiSynth platform. We used the jaw-tongue-hyoid model as a test-case of our simulation platform, Ar- tiSynth, to ensure sufficient capabilities and fidelity for modeling all upper airway structures. Inverse techniques for hard/soft muscle-tissue models i. Formulated trajectory-tracking for muscle-activated dynamic FEM models with constraints. We implemented a technique to automatically compute muscle activations to track a target kinematic trajectory for coupled skeletal and soft-tissue (dynamic FEM) struc- tures. ii. Formulated novel target parameters: constraint forces and stiffness. We extended the inverse formulation to include additional target parameters, including dynamic constraint forces and stiffness, to control muscle activation predictions in redundant models. iii. Predicted beam and tongue muscle activations consistent with muscular-hydrostat theory. We predicted muscle activations needed to move and deform muscular-hydrostat models that are consistent with published theoretical proposals. iv. Predicted plausible muscle activations for lateral jaw move- ment. We predicted muscle activations for lateral jaw movement that are consistent with published literature on jaw physiology. 9 1.3. Dissertation Outline Application to the analysis of segmental jaw surgery i. Created models of segmental jaw surgery with/without re- construction. We developed models of segmental jaw resection and reconstruction through structural alterations to a model of the normal jaw system. ii. Compared mechanical basis of functional deficits between mod- els. We simulated functional deficits in movement and bite force pro- duction that are observed clinically consequent to jaw resection. iii. Applied inverse toolset to predict muscle forces to compen- sate for deficits. We predicted muscle forces that could be used to compensate for functional deficits in a jaw surgery patient. 1.3 Dissertation Outline This dissertation is structured around the three main research contributions. Chapter 2 reviews previous approaches to characterizing orofacial structure and function, including measurement techniques and modeling approaches. Chapter 3 details the jaw-tongue-hyoid model, simulation results, and evalu- ation. Chapter 4 describes the inverse simulation methods developed within the ArtiSynth framework and simulation results for hard and soft tissue models. Chapter 5 details the segmental jaw surgery models used to an- alyze post-operative functional deficits and predict compensatory muscle patterns. Chapter 6 summarizes the dissertation contributions, describes di- rections for future work, and provides concluding remarks. The appendices provide additional background material. Appendix A lists the publications and research talks associated with the dissertation, Appendix B provides a overview of orofacial anatomy, and Appendix C describes the mathematical framework for physics simulation in ArtiSynth. 10 Chapter 2 Related Work This chapter reviews the tools and techniques used to analyze human biome- chanics focusing on studies applied to orofacial structure and function. The traditional approach to biomechanical analysis involves recording observa- tions of human movement and applying statistics to describe relationships within the observations. Such observational studies have generated a wealth of data regarding human motor physiology. However, functional recording techniques are limited by low spatial and temporal resolution, the difficulty of simultaneous recording in multiple modalities, the inability to directly transduce salient variables in humans, and the complexity of the human musculoskeletal and motor systems. These limitations reduce the effective- ness of observational studies for providing a deep understanding of the rela- tionships between peripheral biomechanics and central motor control. The human orofacial region is particularly challenging for observational study be- cause of the 3D nature of soft-tissue deformations and because vocal tract articulators are located within mouth, making them hard to view and ac- cess for movement and muscle recordings. Section 2.1 reviews techniques for functional data recording and discusses related limitations and issues. Recent advances in computation simulation of physical phenomena, such as solid continuum mechanics, have opened new avenues to investigate hu- man biomechanics by developing mathematical representations of the anatom- ical structure. Modeling approaches rely on structural information extracted from medical imaging modalities, such as Computed Tomography (CT) and Magnetic Resonance Imaging (MRI) data, and on mechanical properties of tissues that have been measured in vivo in humans or animals, ex vivo on a test bench, or on cadaver specimens. While based on a range of data sources and founded on a number of assumptions, models can provide explanatory 11 Chapter 2. Related Work power and can fill gaps in functional recordings. Section 2.2 reviews meth- ods for analyzing and quantifying anatomical structure and tissue properties for the purpose of building computational models. Early biomechanical models used greatly simplified representations, such as stick-figure limb models. Increases in computational power and more effi- cient numerical methods have led to models with much richer representations of the 3D structure of bones, muscles, ligaments, joints, and soft-tissues. These advanced modeling techniques are required for complex anatomical structures such as the orofacial region. As models become more accurate, they can be used to analyze casual relationships between anatomical struc- tures and observed functions and ultimately to make predictions extrapo- lating from functional recordings. Section 2.3 describes previously reported orofacial biomechanical models, including different modeling methods for skeletal, soft-tissue, and muscle structures. Biomechanical models require muscle forces as input and simulate re- sulting movements and reaction forces. Integrating models with functional data is challenging because muscle forces are hard to record experimentally. Therefore, inverse simulation methods are important for predicting muscle forces required to match simulated movements with recorded movements. A number of inverse techniques have been proposed and are reviewed in Section 2.4, though little work has been reported for models of the jaw or tongue. As a tool to analyze hypothetical biomechanical situations, models are uniquely suited for biomedical applications, such as analyzing dysfunctional systems and evaluating potential avenues of treatment. In particular, computer- assisted surgery involves integrating computer technologies, including anatom- ical modeling and biomechanical simulation, into surgical planning and ex- ecution. Section 2.5 reviews work on computer-assisted dentistry and oral and maxillofacial surgery. 12 2.1. Functional Data Recording Figure 2.1: Early photographic recordings of horse movement by Muybridge, The Horse in Motion. Subtitle reads, “Sallie Gardner,” owned by Leland Stanford; running at a 1:40 gait over the Palo Alto track, 19th June 1878. 2.1 Functional Data Recording A number of functional data recording techniques have been developed in order to quantitatively observe human movement, forces, and muscle acti- vations. In this section we review how these techniques have been applied to observe jaw, tongue, face, and larynx function. 2.1.1 Movement Jaw and tongue movement has been previously reviewed by Miller [117] and Hiiemae and Palmer [78] respectively. In this section, we organize our review of orofacial movement analysis by recording technique. Optical photographic techniques Quantitative analysis of human movement was enabled by the invention of optical photographic techniques. Muybridge pioneered photographic record- ings of animal and human movement [131]. His studies included horse loco- motion, as pictured in Figure 2.1, illustrating that a horse’s hooves are not 13 2.1. Functional Data Recording in contact with the ground during galloping. Early photographic techniques were modified and refined specifically for the purpose of recording move- ment as described in Chapter 1 of Bernstein [18]. Such techniques included multiple exposures on a single film of a moving subject augmented with high-contrast reflective markers (chronophotography) or incandescent light bulbs (cyclography). Cyclical movements of a stationary subject, such as filing, were recorded on a continuously exposed moving film (kymocyclogra- phy). Mirrors were also used to capture multiple perspectives of a movement of the same exposure in order to reconstruct spatial motion. Quantitative analysis was performed by manual measurements made on the exposed film. Modern video tracking systems use image processing techniques to au- tomatically locate landmarks in each frame of video and multiple cameras to automatically determine the 3D position of landmarks in space. Auto- matic landmark detection is typically aided with active light-emitting diode markers (e.g. NDI Optotrak [134]) or passive reflective markers (e.g. Vicon Motion Systems [207]). Such tracking systems can measure 3D marker po- sitions with sub-millimeter accuracy at sampling rates of 200 Hz or higher. Optical tracking systems have been used to measure jaw movement in six Degrees-of-Freedom (DOF) by rigidly affixing multiple markers to the upper and lower teeth [60, 206]. Our work on the integration of optical jaw tracking data and jaw models is discussed in Section 4.3.1. Facial movement during speech has been recorded with markers distributed over the surface of the face [102, 179]. Section 3.4.1 describes our preliminary work developing a model of face dynamics that could be used to predict muscle forces from optical tracking data of face movements. The main limitation of optical techniques is that line-of-sight is required between markers and cameras. This is problematic for measuring movement of the tongue and other vocal tract articulators, which are internal to the oral cavity. A study by Abd-El-Malek [1] circumvented the problem by us- ing subjects with missing teeth, which provided a window to view inside the mouth during mastication and allowed for visual inspection and illus- tration of tongue shapes. Also, small cameras or optical lenses can be fit through external openings or incisions in the body to view internal body 14 2.1. Functional Data Recording structures. Endoscopy uses a flexible light-transmitting tube that can be inserted through the oral or nasal cavity for a top-down view of the oro- pharyngeal cavity and larynx. Endoscopy is limited to a single top-down view, making 3D laryngeal movements difficult to interpret. However, Selbie et al. [168] used endoscopy to verify a model of 3D cricoarytenoid movement. Other approaches to measuring internal body structures include medical imaging and electromagnetic techniques that do not require line-of-sight. Cineradiography and videofluoroscopy The earliest quantitative recordings of human tongue movement used x- ray cineradiography to analyze kinematics in feeding [13] and speech [149]. Higher resolution, higher frame rate, and lower radiation x-ray recording is achieved with modern Video-fluoroscopy (VF) techniques. VF projects 3D movement to a Two Dimensional (2D) plane and therefore multiple pro- jections, e.g. anteroposterior projection in addition to mediolateral projec- tion, are required to accurately characterize 3D movement [79]. Radiopaque markers have been used to more easily identify landmarks in a similar way as light emitting/reflecting markers are used in optical tracking systems. The position of standardized landmarks for the jaw, tongue, soft-palate, and hy- oid bone, have been digitized from lateral VF recordings used to analyze the kinematic correlation between articulators in both eating and speech move- ments [79, 115]. VF has also been used to observe hyoid position relative to the jaw and tongue during wide jaw opening [130] and vowel postures [23]. A comparative study observed an anterior shift in hyoid position during speech as compared to chewing which was attributed to a need for increased hypopharynx width during speech [79]. VF is also widely used in swallowing analysis and the diagnosis of swal- lowing disorders. The “modified barium swallow” [107] is a standard pro- tocol for assessing risk of aspiration. In the procedure, patients are imaged with lateral projection VF while swallowing a radiopaque liquid (barium), in order to observe whether or not liquid flows into the airway. Section 3.4.3 discusses our preliminary model of hyolaryngeal biomechanics, which we 15 2.1. Functional Data Recording have integrated with a dataset of lateral VF recording of the tongue, hyoid, and larynx movement induced by electrical intramuscular stimulation [31]. Magnetic resonance imaging MRI techniques use the magnetic properties of hydrogen atoms to create images of soft-tissue structures in the body [158]. Static MRI techniques have been used to capture 3D vocal tract shape during held vowel pos- tures [15, 52]. The tissue to air boundary provides high contrast changes in the image, making the airway highly visible. Maintaining a particular tongue posture limits the scan time, which reduces the spatial resolution of the scans. Cricothyroid articulation [191] and tongue muscle structure [190] have also been analyzed during static vowel postures with MRI. Dynamic MRI is limited by long acquisition durations and high speed temporal sam- pling comes at the cost of spatial resolution. Cine-MRI involves acquisition of a time-series of single slice MRI images and requires multiple repetitions of the same speech utterance in order to reconstruct a dynamic MRI image. Tagged cine-MRI has been used to track the position of internal points in the tongue tissue through the cine sequence and provides a measure of local tongue deformation [142]. MRI data of 3D tongue surface shape and internal deformation is a promising modality for integration with 3D biomechanical tongue models, as discussed in Section 4.3.2. Ultrasound Ultrasound (US) imaging uses echoes from pulsed US waves emitted into tis- sue to determine tissue impedance, whereby rapid changes in tissue impedance, such as boundaries between different tissue types or between tissue and air, cause reflections that can be transduced. US has been used to visualize the mid-sagittal contour [182] and the 3D shape [187] of the tongue surface in speech, which is visible due to US reflection at the tongue surface to air boundary. US recordings are fast, but spatial registration of the tongue con- tour to the mandible or palate can be challenging. US has also recently been used to characterize tongue movement in glossectomy patients and reported 16 2.1. Functional Data Recording observations of higher tongue velocity during certain speech movements as compared to normal subjects [159]. Electromagnetic tracking Electromagnetic (EM) tracking systems measure the position of markers within a known EM field by sensing the current induced in small coils within the markers. EM systems do not require line-of-sight and modern systems are self-calibrating (e.g. Polhemus Fastrak [155], Carstens AG500 [34], NDI Wave [135]), however their accuracy degrades if the EM field is distorted by metallic objects in the tracking workspace. EM tracking systems have been proposed specifically for measuring jaw movement, including the ki- nesiograph and the sirognathograph, though optical jaw tracking systems provide higher accuracy and faster sampling rates, as described above. Using EM tracking systems for measuring vocal tract articulations is called Electromagnetic Articulometry (EMA) and has been widely used to analyze tongue and jaw movement in speech [52, 148]. In such studies, markers are glued to the surface of the tongue, lip, and teeth. EMA systems are limited in the number of markers that can be tracked simultaneously, typically up to eight. Therefore, speech studies usually arrange the mark- ers in the mid-sagittal plane since speech movements are mostly bilaterally symmetrical. A sparse sampling of points on the surface of the tongue does not provide detailed spatial information as to the tongue’s 3D shape, but is it effective in recording the kinematics of a few points of interest, e.g. the movement of the tongue tip. The integration of EM tongue tracking data and dynamic tongue models is discussed in Section 4.3.2. Measuring mastication or swallowing kinematics with EMA can be prob- lematic because markers glued to the tongue and teeth tend to fall off during chewing movements. Also, markers placed too far posterior on the tongue may initiate a pharyngeal reflex (gag reflex). Despite these limitations, a recent study has reported EMA recordings for liquid swallowing that show quantitative coordination between the jaw and tongue. 17 2.1. Functional Data Recording Electropalatography Electropalatography (EPG) devices embed an array of contact sensors within a prosthetic palate in order to recorded time-varying patterns of contact be- tween the tongue and palate during tongue movement. In separate studies, EPG was combined with US [188] and MRI [52] to assess tongue movement in speech. 2.1.2 Forces Movement in a biomechanical system is created by muscle forces, which are considered internal forces in the system. External forces arise when a biomechanical system contacts with the environment. Typically only exter- nal forces are available for measurement. Here we discuss mechanisms to record external forces in the orofacial system. Tongue-palate force transducers Pressure between the tongue and palate can be measured with fluid-filled bulbs placed between the tongue and hard-palate. A study of 853 nor- mal subjects reported a maximum tongue-palate pressure of 40.4 ± 9.8 kPa (mean ± standard deviation) for adult males aged 40-49 and showed an age- related decrease in maximum pressure in males, which is a factor in swallow- ing disorders [201]. Recent advances in electronic force sensing technology has led to a new device, which is similar to EPG but capable of measuring force at a number of discrete point on the hard palate [83]. This device has been used to record spatial and temporal patterns of tongue-palate pres- sure force during chewing [82] and swallowing [137]. A preliminary study measured tongue-tip pressure against the anterior hard palate during speech (repeated /ta/ sequences) with a single point force sensor embedded in a full upper denture [91]. Section 3.3 discusses our comparison of maximum vol- untary tongue-palate pressure recordings as a means to validate the tongue muscle force levels in our biomechanical model. 18 2.1. Functional Data Recording Bite force transducers A number of different devices have been proposed for measuring force gen- erated between the upper and lower teeth during jaw function [57]. Bite force during maximal voluntary clench has been measured at different lo- cations around the dental arch, illustrating that tooth forces are largest at the posterior molars and diminish in magnitude at the anterior incisors [76]. Maximal first molar bite force for an intact mandible is within the range of 216-740N [210]. Manometry Manometry uses a flexible tube with embedded pressure sensors to measure pharyngeal pressures [94]. The manometer tube is inserted through the nasal cavity, along the posterior pharyngeal wall and the lowest pressure sensor is placed at the level of the esophageal sphincter in order to calibrate its spatial position. The dataset used to develop our hyolaryngeal model, discussed in Section 3.4.3, includes manometric recordings of pharyngeal pressures during intramuscular stimulation of tongue and laryngeal muscles as well as during liquid swallowing. 2.1.3 Muscle activity Electromyography (EMG) involves transducing electrical signals associated with muscle activation; however EMG can be difficult to record for small, deep muscles in the head and neck, and the relationship between EMG and muscle force is complex for dynamic movements [173]. As discussed in Chapter 4, movement and contact forces are easier to measure directly than muscle forces, making model-based estimation of muscle forces during movement an attractive option in future clinical experiments. EMG has been used in combination with movement recordings in an attempt to better characterize muscle activity. A number of EMG studies have been reported for the jaw, including Moller’s pioneering work on jaw muscle activity during chewing [124] with twenty-six subjects. Other stud- ies have investigated activity of the medial [70] and lateral [214] pterygoid 19 2.2. Structural Data Measurement muscles, which are deep muscles and challenging to access and susceptible to crosstalk from adjacent muscle groups. Fine-wire EMG recordings of regional activation in the upper and lower heads of the lateral pterygoid muscle have also been performed [85, 129] with wire locations within the muscle verified through CT imaging of the subject post-recording before the wires were removed [139]. These studies found differential recruitment in different regions of the muscle and support the hypothesis that the upper and lower heads of the lateral pterygoid muscle functionally co-contract and are both “jaw opener” muscles. In the tongue, EMG recording has proven to be a significant challenge due to the interdigitation of different muscle groups and wire movement during large tongue deformation; however, a few studies have reported EMG recordings for vowel tongue postures [49, 121]. In a recent EMG study, the relative contributions of genioglossus and intrinsic tongue muscles were compared in protrusion tasks and found that both contributed to protrusive tongue movement, but that the intrinsics alone were recruited for generating protrusive force against an external resistance [154]. 2.2 Structural Data Measurement The physical arrangement of bones, soft-tissue, and muscles is important for interpreting the functional observations discussed in the previous section. Orofacial anatomy is reviewed in Appendix B and in this section we discuss three fundamental methods used in anatomical investigation: measurement of cadaver specimens, imaging of living humans, and bench testing of excised tissue samples. 2.2.1 Measuring anatomical structure Traditional methods of study in anatomical science center on the dissection and measurement of cadaver specimens. Skeletal structure is well preserved in cadaver specimens and anatomical features of bone surfaces can be used to determine muscle attachment sites. For example, the mylohyoid line, a 20 2.2. Structural Data Measurement prominent ridge along the inner surface of the mandible, is the attachment site of the mylohyoid muscle, which forms the floor of the mouth. Jaw muscle size and fiber properties have been reported through cadaver mea- surements [203–205]. These measurements form the basis for the muscle properties used in our jaw model described in Section 3.1.1. Tongue muscu- lature has also been examined through cadaveric examination resulting in detailed morphological descriptions of the 3D muscle shapes and fiber struc- ture [120, 192]. These studies form the basis for the spatial and functional muscle definitions used in tongue models, including our model as described in Section 3.1.2. Histology is another analysis technique that involves microscopy of thinly sliced tissue specimens. Staining techniques are used to highlight different tissue structures in the specimen. Histology has been used to determine the percentage of different muscle fiber types (fast-fatigue, fast, or slow type fibers) in tissue samples from different muscle groups. The Visible Human Dataset (VHD) is a macroscale histology of the entire body [132] consisting of photographed slices of frozen cadaver specimens in 1 mm slices for a male specimen and 0.3 mm for a female specimen. The VHD was used to help define the muscle geometry in previous tongue models [29]. The principle limitation of cadaver studies is the fact that the soft-tissue shape and structure degrades shortly after death, reducing the accuracy of muscle size and shape estimates. This is particularly problematic for the tongue, due to its complex arrangement of muscle fibers. Also, anatomi- cal measurements, drawings, and photographs of cadaver dissection do not provide rich information about the 3D nature of the structures under inves- tigation, which is critical for 3D modeling. 2.2.2 Imaging anatomical structure Medical imaging technology has revolutionized our ability to observe internal structures of the body. Modern techniques provide digital, volumetric 3D datasets revealing 3D anatomical structures. Importantly, medical imaging provides access to living subjects, though x-ray imaging modality usage is 21 2.2. Structural Data Measurement CT Image Cone-Beam CT Image (a) (b) Figure 2.2: A comparison of medical computed tomography (a) and cone- beam computed tomography (b) images of different subjects. Medical CT captures a wider range of tissue density at a higher radiation dosage than cone-beam CT. CT image courtesy of Dr. Pierre Badin, Gipsa-Lab. limited due to radiation exposure. Medical image technology is an active area of research and new imaging methods and image processing techniques are improving in terms of spatial and temporal resolution and contrast [158]. Computed tomography CT uses multiple x-ray projections from different angles to compute 3D volumetric data of dense tissue. CT imaging provides high contrast and high spatial resolution (0.5 mm3 voxel) images and is well-suited for bone imaging since dense tissues absorb x-rays. Multi-slice systems use multiple x-ray detectors simultaneously to rapidly image and reconstruct CT data. The principle drawback of CT imaging is radiation exposure. Cone-Beam Computed Tomography (CBCT) reduces the radiation ex- pose to the patient, but images a narrower range of tissue density than medical CT, limiting its use to bone tissues alone. A comparison of CT and CBCT images are shown in Figure 2.2. Our original jaw-hyoid model, described in Section 3.1.1, used skeletal structure derived from CBCT data, while the new jaw-tongue-hyoid model morphology was registered to CT 22 2.2. Structural Data Measurement T1-Weighted MRI T2-Weighted MRI (a) (b) Figure 2.3: A comparison of T1 weighted (a) and T2 weighted (b) magnetic resonance images of the same subject. data (Section 3.1.3). Magnetic resonance imaging MRI uses magnetic resonance in hydrogen atoms to create high-resolution (1 mm3 voxel) images of soft-tissue structures. MRI parameters can be ad- justed to highlight different tissue types. For example, a comparison of T1- weighted versus T2-weighted images is shown in Figure 2.3. High-resolution MRI has been used to study the spatial extent, shape, and path of extrinsic tongue muscles [190] as well as the shape of the laryngeal cartilages [169]. Diffusion tensor MRI can be used to isolate fiber directions within muscle tissue and has been applied to the tongue [64]. Image segmentation Image segmentation is the process of specifying regions within medical im- ages that correspond with tissue boundaries. An example is segmenting the boundary of the mandible from adjacent soft-tissues. The difficulty of image segmentation is related to the degree of intensity variation across the tis- sue boundary. Segmentation of well-defined tissue boundaries, such as the 23 2.2. Structural Data Measurement boundary between bone and soft-tissue or soft-tissue and air in CT data (see Figure 2.2), can be achieved with simple thresholding methods. Segmenta- tion of soft-tissue to soft-tissue boundaries, such as between adjacent muscle groups in MRI data, is more challenging and typically requires semi-manual or manual specification of regions by a trained anatomist. Automatic image segmentation is a active area research, especially with respect to brain and heart imaging domains. Image segmentation in the upper airway domain is complicated by soft-tissue movement artifacts between slices or scans mak- ing automatic segmentation techniques less effective. 2.2.3 Measuring mechanical tissue properties Mechanical properties of human and animal tissues have been examined through ex vivo mechanical testing (see [48] for review). There are only a few reported studies on orofacial tissue properties. Human tongue and cheek tissue properties, which are incorporated in our tongue model (Sec- tion 3.1.2), have been experimentally examined though indentation testing on cadaver specimens [63]. Also, a suction-based device has been devel- oped for in vivo mechanical testing [164]. Soft-palate tissue properties have been measured ex vivo with cadaver specimens [21]. Laryngeal tissue and muscle properties have been experimentally tested on canine tissues [5, 89]. Given the limited literature, experimental measurement of orofacial and up- per airway tissue mechanics is an open research area. In particular, studies examining the anisotropic properties of tongue tissue are needed. New techniques for in vivo measurement are needed to assess living tissue mechanics and mechanical changes during muscle activation. Elastography, with US and more recently MRI, is making progress on the analysis of in vivo tissue mechanics. The majority of elastography studies have focused on assessing tissue elasticity for tumor detection, though recent studies have attempted to estimate elasticity against ground truth data (see Mariappan et al. [112] for a recent review). There are a number of challenges in adapting the techniques to the orofacial region, including the required mechanical vibration of the tissue during the imaging protocol. 24 2.3. Biomechanical Modeling 2.3 Biomechanical Modeling The previous two sections have discussed common methods used by oral bi- ologists and anatomists to gather information about the function and struc- ture of orofacial anatomy. In a general sense, modeling is the process of abstracting general information about an anatomical system from a spe- cific set of physical measurements. Models can take on a variety of forms. Some are created as a means of abstracting higher-level forms of information from specific datasets. For example, Principal Components Analysis (PCA) models reduce a high-dimensional dataset to a small number of principle di- mensions, which can provide insight into those aspects or characteristics of the data that are most important. Other models are created by combining and synthesizing multiple datasets and/or data modalities into a common framework. For example, a geometric model of 3D anatomy may include bone shapes extracted from CT data together with muscle shapes found in MRI data. Synthesis models have the added complication that multi- ple datasets must be transformed into a common format, or in the case of anatomical data, they must be co-registered such that they are spatially con- gruent. The models developed in this dissertation are biomechanical mod- els that combine structural anatomical information with mechanical tissue properties and functional recordings of movement and force. This section describes the process of biomechanical model creation, focusing on skele- tal, soft-tissue, and muscle modeling approaches. In addition, previously proposed biomechanical models of the orofacial structures are reviewed. 2.3.1 Rigid bone modeling Bones are often approximated as rigid bodies in biomechanical models of gross body movement. Multibody techniques [172] can efficiently simulate the dynamics of numerous rigid bodies along with constraints and contact between bodies. Commercial multibody simulation packages include Solid- Works [43] and ADAMS [126]. Several open simulation systems and archi- tectures have been presented to the biomedical community in recent years. OpenSim [46] is a multibody simulator designed for musculoskeletal analysis. 25 2.3. Biomechanical Modeling Our ArtiSynth modeling platform is capable of articulated rigid-body sim- ulation with contact and constraints as well as deformable body simulation as discussed in Section 2.3.2. Geometry Surface meshes are used to define bone boundaries in 3D space and are composed of a set of 2D triangles or other 2D polygons. They are used for visualizing 3D surfaces, as well as to define muscle insertion and origin locations, and to detect collisions between structures in dynamic simulations. Surface mesh generation Generic meshes for bones have been created by anatomists and artists. For example, average-valued dimensions have been reported for the mandible [217] and the laryngeal cartilages [50]. Generic meshes can be co-registered to a subject; however for subject-specific mod- els it is usually more desirable to generate bone mesh surfaces directly from medical imaging data of the subject, if such data are available. Surface mesh generation follows directly from the image segmentation techniques described in Section 2.2.2. Once region boundaries are defined as a closed set of voxels in a 3D image dataset, a surface mesh can be constructed along the boundary with techniques such as marching cubes. Mesh decimation is commonly used to reduce the number of triangles in a mesh while attempt- ing to maintain the same surface shape. However, the decimation process can cause mesh artifacts, such as holes or poorly-conditioned (skinny) tri- angles. A number of open source packages for surface mesh processing and editing are available, including Blender [65] and MeshLab [208]. Mesh registration Models with multiple meshes from different sources or datasets require those sub-meshes to be co-registered into the same spa- tial reference frame. Registering bone meshes from the same subject can be done by finding a rigid transformation (translation and rotation) to mini- mize the sum squared distance between a set of corresponding landmarks defined on each mesh [84]. Registering bone meshes from different subjects, given the wide individual variation in orofacial structure size and shape, 26 2.3. Biomechanical Modeling requires affine or non-linear transformation. Affine transformations can be computed from corresponding landmarks in a similar fashion as rigid trans- formations, but include non-uniform scaling and shearing deformations. An affine transformation will capture the gross differences in size and shape; however, finely detailed shape differences require a non-linear transforma- tion or morphing. Bucki et al. [30] describes a non-elastic mesh-based reg- istration method that we used to adapt generic jaw and skull meshes to CT data of a specific-subject for our jaw-tongue-hyoid model, as discussed in Section 3.1.3. Dynamics Solving for rigid body dynamics requires mass and inertia information for each body. A common approximation computes an inertia based on the sur- face mesh geometry assuming a constant density. Also, the mass of muscle tissue and soft tissue structure is commonly lumped into the mass and in- ertia of the skeletal structures. Experimentally measured masses have been reported for the mandible [217] and the laryngeal cartilages [50]. 2.3.2 Deformable tissue modeling Deformable biological tissues include muscles, connective tissue, fat, and mu- cosa. Cartilage and bone structures are also deformable under sufficiently large loads, though they are often approximated as rigid in models more concerned with their gross movements than their internal deformations. De- formable structures are commonly modeled using FEM approaches [17, 22] that compute approximate solutions to the partial differential equation gov- erning continuum mechanics by spatial and temporal discretization. Commercial software packages for FEM modeling include ANSYS [12] and SIMULIA [42] and are primarily designed for the structural mechan- ics analysis of synthetic materials as opposed to biological tissues. Given the computational complexity of FEM techniques, solution times can be very slow: a one-second simulation of the FEM tongue model described in [29] can require many hours of computing time. A number of open source 27 2.3. Biomechanical Modeling research-oriented toolkits have also been developed, including FEBio [211], a finite element toolkit with special support for tissue modeling and some support for rigid bodies, contact and constraints. Systems geared toward surgical training include Gipsi [35] and Spring [125]. Sofa [6] provides a gen- eral software architecture in which models can be partitioned into different submodels for simulating appearance, behavior, and/or haptic response. To the best of our knowledge none of these open simulation systems yet provide an interactive environment with fully coupled FEM/multibody capabilities. Our ArtiSynth platform combines both FEM and multibody capabilities within an interactive graphical environment and is described in Appendix C. Geometry Volumetric meshes are used to define a volume in 3D space and are composed of a set of 3D tetrahedrons, hexahedrons, or other volumetric elements. The volumetric mesh defines the spatial discretization of the continuum mechanics equations in FEM. For accurate FEM solutions, volumetric mesh elements have quality requirements concerning their size and shape, such as minimum aspect ratios (see Shewchuk [174] for detailed discussion of FEM mesh quality requirements). Volumetric mesh generation Volumetric mesh generation presents a greater challenge than surface meshes as it involves creating a 3D mesh with 3D elements that fill a volume without holes or unconnected nodes. Henshaw [75] provides a current review of the state-of-art in automatic vol- umetric mesh generation. Tetgen [178] is a freely available software package for generating tetrahedral volumetric meshes from surface meshes and com- mercial FEM software packages include similar tools. Hexahedral meshes are more challenging to generate than tetrahedral meshes, but are more desirable for FEM analysis. This is because hexahedral meshes are less susceptible to locking (erroneous increase in stiffness) when simulating incompressible or nearly-incompressible materials. Manual mesh creation is also common, as is the case for our reference tongue model mesh (Section 3.1.2). 28 2.3. Biomechanical Modeling Mesh registration Deformable mesh registration is also a more challeng- ing problem than rigid or affine registration of bone meshes. Deformable tis- sues can require non-linear transformations, or morphing, to provide a good fit between subjects. Bucki et al. [30] report an energy-based mesh mor- phing and registration method that includes constraints on element quality and has been used to adapt face meshes to a specific subject. An initial affine transformation is useful as a good initial guess for more complex non- linear morphing algorithms, which are commonly formulated as optimization problems and therefore can be sensitive to initial conditions. Dynamics Solving for deformable body dynamics requires parameters relating to their mass, stiffness and damping. These parameters depend on the material used to represent the deformable tissue. For example, a linear material uses Young’s Modulus (associated with stiffness) and Poisson’s ratio (associated with compressibility). Non-linear materials require additional parameters as stress varies non-linearly with strain. Elasticity parameters derived from tissue measurement, as described above in Section 2.2.3, are dependent on the chosen material and are fit to the data in order to recreate observed deformations over a range of prescribed loads in the model. 2.3.3 Muscle modeling In addition to passive skeletal and soft-tissue structure, biomechanical mod- els require muscle force in order to simulate active movements. Muscle models include a definition of a muscles spatial extent, its attachment to surrounding structures (origin and insertion locations), and its force gener- ating capabilities and dynamics. Geometry Muscle geometry defines how muscle forces are transmitted to surrounding structures by defining muscle attachments at origin and insertion locations. 29 2.3. Biomechanical Modeling Muscle geometry that accounts for a muscle’s volumetric extent can also define contact forces between the muscle tissue and adjacent structures. Line-based The most common geometric model for muscle forces is one or more connected line segments. Muscle forces are transmitted between the two end-points (insertion and origin sites) and intermediate points are usually treated as ideal pulleys that transmit force without energy dissipa- tion. Line-based muscle models are simple and effective for skeletal muscles that have small origin/insertion areas and direct paths, which is the case for many jaw and laryngeal muscles. Therefore, we use piecewise-linear muscle geometries in our jaw model (Section 3.1.1). Broad or flat muscles, with large origin/insertion areas, can be approximated by a number of muscle “lines” in parallel. Origin/insertion points are typically chosen as the centroid of attachment regions, however a muscle’s effective line of action can be modi- fied by its fiber architecture. Muscles with complex paths, wrapping around bones and/or through joints, are approximated by defining “waypoints” to create a piecewise linear path wrapping around adjacent structures [46]. The “waypoint” formulation can lead to inaccurate transmission of forces to adjacent structures, especially in complex muscle paths, such as with the wrist, shoulder, or knee. Spline-based Spline-based muscle paths are an extension of line-based models to provide additional path DOF and more accurate wrapping of muscles around bones and through joints. This approach has been used to model the complex structure of the forearm and hand, demonstrating force transmission and sliding between musculotendons and bones surfaces [189]. Volumetric A more sophisticated approach for modeling 3D muscle struc- tures is volumetric muscle models that allow for the spatial size of the mus- cle and attachment region to be included in the model. Volumetric muscle models use deformable tissue modeling techniques, as described above, to simulate passive mechanical properties of muscle tissue [193]. Active forces are incorporated in the FEM mesh along the muscle fiber lines-of-action with 30 2.3. Biomechanical Modeling either discrete line segments, as is used in our reference tongue model (see Section 3.1.2 and [29]) or transversely-isotropic FEM materials [212]. Vol- umetric methods are required for muscle-activated soft-tissues such as the face, tongue, and soft-palate, and have also been used to model the spatial extent of the jaw muscles [160]. Pennate Pennate muscles have muscle fibers oriented obliquely to the muscles’ principle line-of-action. The pennate structure is thought to be a mechanism to trade-off a muscle’s capacity to shorten with its capacity to generate force. The jaw closing muscles, which are primarily intended to generate bite forces during mastication, have a pennate structure (see Section B.2.2). For example, the masseter muscle is noted to have multiple muscle sheets of pennate fibers oriented at different angles, which has been suggested as a mechanism for maintaining bite force throughout a range of jaw closing rotation [68]. Pennate-muscle structure is commonly neglected in biomechanics models, making it a interesting area for future investigation. Dynamics A number of parameters are required to characterize the dynamics of mus- cle tissue, including both its passive properties (similar to deformable tissue parameters) and also its force generation under activation. These parame- ters are typically non-linear and are dependent on the chosen muscle model formulation. Sophisticated models of muscle dynamics can also include acti- vation dynamics and fatigue. A common approximation lumps muscle mass into the mass of adjacent bone structure. This approximation does not take into account spatial changes in muscle mass due to muscle lengthen- ing/shortening during movement and can cause significant errors in dynamic simulation [140]. The total mass of jaw muscles [205] is roughly equivalent to the mass of the mandible [217], which would suggest that jaw muscle mass should be explicitly accounted for in dynamic models. However, jaw movements are slow during chewing and small during speech and therefore the inertial effects are less significant. Volumetric muscle models, such as 31 2.3. Biomechanical Modeling Te n si o n Length To ta l  T en sio n Pa ss iv e St re tc h Contractile Re st   L en gt h Muscle TendonCE SE PE CE - Contractile Element SE - Series Elastic Element PE - Parallel Elastic Element (a) (b) Figure 2.4: The standard Hill-type muscle model representation [80] illus- trated with mechanical schematic (a) and force-length characteristics (b). our tongue model, incorporate muscle tissue mass into their dynamics. Hill-type model The Hill-type muscle model [216] is widely used in dy- namic modeling of musculoskeletal systems, including previous jaw models. It includes a non-linear force-length relationship in a three-element model, as shown in Figure 2.4. Muscle properties Another muscle model derived from measurements on feline hindlimb muscles reports detailed measurements of activation dy- namics [27]. Also, canine intrinsic laryngeal muscle properties have been measured with ex vivo experiments [5]. For the jaw muscles, Hill-type model parameters have been derived from measurements of muscle Cross-Sectional Area (CSA) and fiber lengths on cadavers [203–205]. Also, a model-based study reported that low passive tensions in the jaw closer muscles were re- quired to achieve maximum wide gape [146]. We use the Hill-type muscle formulation for the jaw muscles in our model (Section 3.1.1) with CSA pa- rameters given in Table 3.1. 32 2.3. Biomechanical Modeling 2.3.4 Isolated orofacial models Biomechanical models of individual structures of the human upper airway have been developed and used since the 1970s. Model complexity has in- creased due to both the acquisition of new knowledge about anatomical, neurophysiological, and physical characteristics of the articulators, and the vast growth in computational capacity. Early models were based on a 2D mid-sagittal representation of the airway [144, 150, 151, 162], whereas recent models attend to the 3D structure of the articulators. Jaw models 3D models of jaw biomechanics commonly use line-based muscle models to represent the jaw muscles as well as 3D representations of the TMJ and teeth [101, 146]. Recent models have been used to examine TMJ loading during chewing [71], open-close movements [199], and after jaw growth [45]. A few models have used FEM to model volumetric muscles [160] or the articular disc [100]. Tongue models Early tongue models focused on planar 2D representations because speech production is bilaterally symmetric and thought to be primarily a mid- sagittal phenomena. The move to 3D was motivated by the fact that the tongue’s mid-sagittal shape is controlled in part to the mediolateral expan- sion and contraction. Dang and Honda [41] investigated speech production movement with a 2D spring-mass model of the tongue using a mid-sagittal mesh with mediolateral thickness. Recent tongue models employ 3D FEM meshes with hyperelastic constitutive models to better represent non-linear tissue properties [16, 29, 62]. Other oral and facial models A number of 3D biomechanical face models have been proposed mainly for synthesizing visual speech. Early models used a spring-mass system ap- 33 2.3. Biomechanical Modeling proach [194], whereas recent models have employed FEM methods. Sifakis et al. [179] proposed a face model with detail volumetric muscle structure and quasi-static FEM coupled with motion-capture data of face movement during speech utterances. Nazari et al. [133] used line-based muscle models along with a dynamic, hyperelastic and nearly incompressible FEM in or- der to simulate speech movements, including lip rounding and protrusion. Larynx models have primarily focused on intrinsic laryngeal muscles and vocal fold mechanics [90, 122] as opposed to the movement of the hyola- ryngeal complex within the neck. The mechanics of the intrinsic larynx is primarily applicable to speech production, but is also applicable to airway protection during swallowing. Very few biomechanical models of the human soft-palate and pharynx have been reported. Previous soft-plate models used highly simplified mid-sagittal plane 2D FEM meshes and simplified muscles [19, 145]. A more complex pharynx model that used a sagittal 2D FEM model with mediolateral thickness to represent detailed oropharyngeal structures was used to analyze the mechanics of OSA [86]. 2.3.5 Coupled orofacial models While many models of individual orofacial structures have been developed, few models combining and integrating adjacent structures have been re- ported. The focus on individual structures in previous work is likely due to the significant technical challenges with properly modeling the dynamical coupling between soft bodies (tongue, lips, soft palate) and hard structures (jaw, hyoid bone, hard palate) in 3D. Sanguineti et al. [162] have reported a highly simplified 2D jaw-tongue- hyoid model with a 2D FEM to represent the tongue. Fang et al. [54] have reported a 3D model of the jaw-tongue-hyoid using a discrete mass-spring network to represent the tongue tissue. Discrete mass-spring models are known to be less numerically stable and less able to accurately characterize non-linear, incompressible biological tissues as compared to FEM. So, to our knowledge, our approach is unique for providing a 3D modeling framework simulating full dynamical interaction between soft and rigid articulators us- 34 2.4. Inverse Methods ing FEM/rigid-body techniques, as discussed in Chapter 3. Such dynamic coupling is important, especially in speech production where it has been clearly shown that consideration of the dynamic interaction of vocal tract bony and soft structures is needed to correctly account for orofacial dynam- ics [162]. 2.4 Inverse Methods Model-based techniques for estimating muscle forces during movement have been proposed for a variety of musculoskeletal systems. Such inverse meth- ods are useful because muscle force is difficult to measure in vivo. These techniques have been recently reviewed by Erdemir et al. [53] and are briefly described in this section. Some of the techniques are concerned with mod- eling the human motor system and incorporate theories of the human mo- tor control and learning, whereas others have a direct correspondence to recorded data. Humans exhibit tremendous motor skill and dexterity. Therefore, it fol- lows that theories of human motor control may inform inverse methods for biomechanical models. In a complementary manner, computational models provide a means to evaluate theories of human motor control that are based on experimental observation. Two prominent computational models of mo- tor control include the Equilibrium Point Hypothesis (EPH) and supervised learning. Stiffness is also an important mechanism in motor system models. Inverse-dynamics techniques use numerical methods to provide estimates of muscle force with direct relationship to recorded data. Inverse-dynamics techniques, therefore, can be considered “data-driven” approaches as they are formulated with a close association to recorded kinematic and other data describing the motor task under investigation. Equilibrium-point hypothesis The EPH approach controls a model’s individual degrees of freedom by changing the rest-length of spring-like muscle models to drive the model 35 2.4. Inverse Methods between successive quasi-static positions [55]. EPH has been used to study jaw motion in speech, including neuron-motor control [162] and the effects of jaw stiffness [175]. It has also been used with a biomechanical tongue model [29] to predict muscle forces need to achieve static vowel postures. EPH models incorporate sophisticated muscle dynamics and reflexes. The quasi-static target position assumption implicit in the EPH is controversial and has been disputed [66] and defended [56] by a number of researchers. Supervised Learning Machine learning techniques have been used in attempts to recreate prop- erties of the human motor system in computational models. Supervised learning is typically implemented in an artificial neural network structure and uses labeled training data to adjust the parameters in the network in order to reduce prediction error. Learning an inverse model of the motor system requires a motor error training signal, whereas sensory error (mo- tion error) is more readily available to the motor system. The distal-error technique proposed by Jordan and Rumelhart [93] uses a forward model to convert sensory error into motor error for supervised learning of a feedfor- ward controller. Alternatively, the feedback-error learning model proposed by Kawato et al. [95] uses a fixed feedback controller to map sensory-error to motor-error for learning an inverse internal model. Porrill et al. [156] have proposed a recurrent control model that involves supervised learning of a forward model using a sensory-error training signal. Supervised learning techniques for complex inverse mappings typically suffer from slow learning rates and require large amounts of labeled data, which limits their applica- tion to biomechanics simulation. The methods also rely on an approximate inverse solution as a basis (“reference structure” in feedback-error learning model [95] and the “brainstem” in Porrill’s recurrent system [156]). Stiffness Modeling Stiffness is a functionally important mechanism in the biological motor sys- tem [26] that is frequently neglected in computational models. In optimiza- 36 2.4. Inverse Methods tion approaches, the use of a minimum muscle energy constraint is equiv- alently a minimum stiffness solution. Stiffness in a biomechanical system is increased with an increase in muscle activation, through co-activation of antagonist muscle pairs, and by short-latency stretch reflex gains. The control of arm stiffness through co-activation of antagonist muscles has been illustrated in a simple model of a single joint [81]. The joint torques induced by flexor and extensor muscles are given by τF = (T −Ko∆q) aF and τE = (T −Ko∆q) aE , where ∆q is an angular displacement from rest, T is the maximum isometric muscle-induced torque at the joint, and Ko is the intrinsic stiffness of the joint. The net isometric torque at a given joint angle is computed by the sum of the flexor-induced and extensor-induced joint torques (τn = τF + τE), which leads to: τn = T (aF − aE)−Ko(aF + aE)∆q (2.1) and illustrates that torque and stiffness about a joint can be independently controlled by the difference and sum of muscle activations respectively [81]. A generalized formulation of stiffness control through muscle co-activation is described in Section 4.1. Static optimization Static-optimization based inverse dynamics involves estimating net forces in a biomechanical system from recorded kinematics. For example, in simple limb models joint angle trajectories q(t) are measured, numerical differ- entiation is used to estimate joint velocity q̇ and acceleration q̈, and net joint torques are estimated from the kinematic variables using knowledge about the inertia of the limb segments. The net forces in a biomechan- ical system, however, are insufficient to characterize muscle force due to muscle redundancy. Therefore, static optimization is performed at each in- stant of movement to decompose net forces (e.g. joint torques) into muscle forces [39, 215]. The optimization relies on an instantaneous cost function, such as minimum excitation, to resolve muscle redundancy. One drawback of traditional inverse dynamics is that errors in recorded position data are 37 2.4. Inverse Methods magnified in differentiation for velocity and acceleration [180]. A linear pro- gramming formulation of the static-optimization approach was proposed for point-to-point movements of the jaw [98] and is the only previously reported inverse method for jaw movement modeling. Dynamic optimization Dynamic optimization techniques have also been proposed to compute mus- cle forces over the duration of a movement task. These techniques are used to predict muscle activity in motor tasks defined more broadly than specific kinematic trajectories. For example, while trajectory-tracking methods will specify a time-series of positions targets to describe the movement, dynamic optimization methods may only define the start and end positions. Dy- namic optimization allows for more biologically plausible optimality criteria to be introduced into the solution, such as minimizing metabolic energy expenditure per unit distance over the duration of the movement task [8]. Full dynamic optimization, however, remains computationally expensive and current methods are intractable for complex biomechanical models. A com- parison of dynamic and static optimization found little difference in the resulting computed muscle activations, assuming a perfect inverse dynam- ics estimation of joint torques [9], which suggests that static optimization techniques can be sufficient at least for some motor tasks. Optimal control Another approach to dynamic optimization uses optimal control methods to compute a optimal trajectory of muscle activity for a given motor task. A particular formulation, called optimal feedback-control, theorizes that only task-related parameters are optimally controlled by the motor system, and has reported promising results in predicting muscle forces in limb move- ments [197]. The optimal feedback-control method has also been extended to incorporate a hierarchal structure [198], which relates to observations by motor physiologists regarding hierarchical organization of the human motor system [106]. 38 2.5. Biomedical Applications Forward-dynamics assisted tracking To mitigate the inaccuracies of static optimization inverse dynamics, re- cent techniques employ a hybrid approach and compute an inverse solution of muscle activations that accurately drive a forward dynamics simulation through the target kinematics. This approach was termed “forward dynam- ics assisted data tracking” by Erdemir et al. [53] and can improve accuracy as tracking errors in the forward dynamics simulation are used in a feedback control law to adjust the inverse solution. Computationally efficient algo- rithms have been proposed for gait analysis using either static per-timestep optimization [196] or optimal control [171] to decompose joint torques to muscle forces. General formulations of the trajectory-tracking approach have been proposed for quasi-static FEM models [179] and dynamic multi-body simulation with musculotendon elements [189]. Combining and extending these two methods, in Chapter 4 we develop a trajectory-tracking method for dynamic FEM models and combined dynamic rigid-deformable models. We also formulate new target parameters, in addition to movement, includ- ing constraint forces and stiffness. The constraint force targets are used to predict muscle activations needed to generate desired reaction forces in the system, such as a target bite force during jaw clenching simulations. The stiffness target can be used to control co-activation of antagonist muscles. 2.5 Biomedical Applications Biomechanical models can be applied to biomedical situations, such as an- alyzing dysfunctional cases and evaluating potential avenues of treatment. Modeling studies in orofacial biomedicine have been recently reviewed by Han- nam [69]. Here we point to a few studies that are most closely related to our efforts in modeling the functional consequences of jaw surgery described in Chapter 5. Jaw surgery In reconstructive jaw surgery, physical 3D models have been applied to the fabrication of custom-fitting mandible prosthesis with rapid- 39 2.5. Biomedical Applications prototyping techniques [20]. Modeling techniques have also been used to design mechanical templates for reshaping bone grafts for use in jaw recon- struction. The template approach has been shown to reduce the duration of jaw reconstruction procedures from two-years to six weeks [195]. A recent biomechanical jaw simulation of distraction osteogenesis (lengthening one side of the mandible) was used to analyze TMJ forces in a single patient before and after treatment [45]. Dental implants Numerous studies have applied biomechanical modeling to the analysis of tooth forces and dental implants (see [209] for a review). Most studies use FEM in order to assess the force distribution of implanted synthetic teeth in order to optimize implant design and location within the mandible. Force distribution in tooth restorations has also been analyzed through 3D FEM tooth models [92]. Maxillofacial surgery Biomechanical face models with static skull struc- tures have been used to predict aesthetic outcomes of maxillofacial surg- eries [36, 123]. These studies kinematically modify skull structure to mimic a surgical procedure, e.g. lengthening of the mandible, and use passive FEM models of the facial tissue to predict the resulting impact on the face sur- face. To date, such models have not been used to predict post-operative facial motion or function. Glossectomy Biomechanical tongue models have been used to simulate the effect of glossectomy [28, 59]. In a similar methodology as our jaw surgery model, these studies modified the structure of a tongue model to mimic tongue resection and reconstruction with free-flap soft-tissue grafts. The reported simulations deal primarily with the effect of stiffening a sub- region of the tongue, representing a legion or reconstruction, on tongue movements. Scar tissue modeling Modeling scar tissue mechanics is relevant to many biomedical applications as tissue scarring arises from burns, radiation ther- 40 2.6. Summary apy, tissue grafts, and surgical wounds [118]. Few physiological models of scar tissue have been reported in the literature. Along with a better charac- terization of scar tissue mechanics, computational models of scar tissue are needed for use in biomedical modeling applications. 2.6 Summary To summarize, a wide range of techniques have been developed in order to observe and analyze human movement. Earliest observations used photo- graphic techniques to capture body movement and modern tracking systems allow for accurate, fast, and automatic quantification of movement. The oro- facial system presents a challenge as functionally important structures are located within the oral cavity and are therefore less accessible for obser- vation and measurement. Medical imaging techniques are used to capture both anatomical structure of bones, soft-tissues, and muscles and dynamic movements of the face, jaw, tongue, and vocal tract. Modeling is a natu- ral extension of observational analysis, whereby anatomical structure and mechanics are combined within a mathematical representation and used to synthesize hypothetical function. Many models have been proposed in the literature for the jaw, tongue, face, and other orofacial sub-components. However, few models have integrated these to create holistic models of the coupled orofacial system. Also, the utility of biomechanical models is lim- ited by a lack of motor control algorithms to automatically simulate pre- scribed motor tasks. As our capabilities to accurately model biomechanical systems and simulate dynamic movements increase, numerous biomedical applications become feasible, including analyzing dysfunction and planning treatment for patients. This chapter presented a review of biomechanics modeling approaches and demonstrated the utility of modeling approaches. We have identified three areas of investigation for orofacial modeling: firstly, isolated models of individual orofacial structures do not adequately represent the coupled nature of the orofacial system; secondly, forward dynamics models alone are insufficient for analyzing orofacial movements and inverse techniques are 41 2.6. Summary required; and thirdly, models of normal orofacial mechanics can be altered to reflect dysfunctional systems and applied to analyze potential treatments. In the following three chapters we describe our contributions to these open research problems. 42 Chapter 3 Jaw-Tongue-Hyoid Model In the previous chapter, we outlined approaches taken by previous researchers to measure, analyze, and model human orofacial anatomy and function. Many models of the face, jaw, tongue, and larynx have been proposed with a variety of modeling techniques and at a variety of modeling fidelities. These models have almost exclusively focused on individual structures in isolation. The degree to which passive mechanical linkages cause functionally impor- tant coupling between structures has not been analyzed. One reason for the lack of integrated models of orofacial components is the complexity and computational cost of FEM simulations. The ArtiSynth simulation plat- form includes state-of-the-art simulation techniques for simulating coupled rigid-body/FEM systems with efficient computational methods and is de- tailed in Appendix C. The coupled mixed-body approach is well-suited to the analysis of gross movements where the motion of bony structures can be approximated as rigid while also attending to deformation for soft-tissue structures. The structural anatomy of the head and neck includes a complex ar- rangement of bones, muscles, ligaments, and soft-tissues, as discussed in Section 2.2 and reviewed in Appendix B. The jaw muscles connect the mandible above to the cranium and below to the hyoid bone. The tongue fits within the oral cavity and attaches to the jaw (genioglossus, geniohy- oid, mylohyoid), hyoid (hyoglossus), soft-palate (palatoglossus), and cra- nium (styloglossus). Muscles exert different force magnitudes in relation to their size and neural activation. Muscle and soft-tissue interconnections between sub-structures transmit forces throughout the system. A compu- tational model of the interconnected muscles and structures permits the analysis of internal forces that are not readily measurable. 43 3.1. Model Creation The structural and force couplings in the orofacial anatomy is thought to be an important aspect of orofacial function. Chewing is considered to be primarily a jaw muscle action, though coordination of the tongue body and buccinator muscle in the cheek is necessary to form and stabilize the food bolus during the chewing cycle. It is likely that passive stretch of facial tissue and the tongue tissue between the jaw-hyoid has a mechanical effect, but its extent is not clear without analysis of passive tissue stiffnesses and muscle forces. In speech, the coupling of jaw-tongue-hyoid structures is readily apparent. It is known that both the jaw and tongue articulators work together to form the vocal tract shape. Also, the level of muscle force is small as compared with chewing and therefore passive coupling is likely more significant. Co-articulation in speech production, the observation that the relative contributions of jaw and tongue movement to a particular vocal tract articulation are dependent on phonetic context [128], may also be due in part to biomechanical coupling. In this chapter, we develop a new dynamically coupled model of the jaw- tongue-hyoid complex in order to analyze the nature and significance of the interconnected structures. Section 3.1 describes the jaw-tongue-hyoid model implementation. Section 3.2 and Section 3.3 present simulations of simple lingual and mandibular motor tasks, focusing on the dynamic interaction between the tongue soft structure and the jaw rigid body. Section 3.4 re- ports preliminary results on integrating the face, soft-palate, and larynx in order to create a complete model of orofacial and upper airway biomechan- ics. Section 3.5 and Section 3.6 discuss the implications of the model and directions for future refinement. 3.1 Model Creation Building on our experience with the 3D jaw-hyoid [71] and the 3D tongue [29] models, we have developed the first 3D jaw-tongue-hyoid dynamical model taking into account full coupling between the FEM tongue model and the jaw-hyoid bony structures. As described in Appendix C, we depend on the various components of ArtiSynth to provide dynamic simulations of inter- 44 3.1. Model Creation MT PT AT SLP ILP DM SM MP DI PD Front Oblique Sagittal Cut-Away Figure 3.1: Front, oblique, and sagittal cut-away views of the jaw-hyoid model. Muscle groups include the anterior, middle, and posterior temporalis (AT, MT, PT), deep and superficial masseter (DM, SM), medial pterygoid (MP), superior and inferior heads of the lateral pterygoid (SLP, ILP), and anterior and posterior bellies of the digastric muscle (DI, PD). actions due to muscle forces of the jaw-tongue-hyoid complex and contact phenomena such as tongue-palate collisions. 3.1.1 Jaw-hyoid model For the jaw-hyoid model, we started from a previously published model de- veloped in ArtiSynth and used to simulate free jaw movements [183] and chewing [71]. The model included rigid-bodies for the skull, jaw, and hy- oid bone, point-to-point Hill-type actuators for the jaw muscles, constraint surfaces for the TMJs, and planar unilateral constraints for teeth contact. The same Hill-type muscle dynamics were used from the original jaw-hyoid model with force capacity proportional to maximum CSA based on previous studies of jaw [104, 146] and tongue [29] muscles and listed in Table 3.1. The instantaneous force generating capacity of the jaw muscles vary non-linearly with length and linearly with shortening velocity consistent with previous jaw models [104]. The skull is fixed in space. The adapted jaw-hyoid model is shown in Figure 3.1. 45 3.1. Model Creation Jaw Closer Muscles Name AT MT PT SM DM MP MaxForce (N) 158.0 95.6 75.6 190.4 81.6 174.8 CSA (mm2) 395 239 189 476 204 437 Jaw Opener Muscles Name SP IP DI MaxForce (N) 28.7 66.9 40.0 CSA (mm2) 72 167 100 Tongue Intrinsic Muscles Name GGA GGM GGP VERT TRANS IL SL MaxForce (N) 32.8 22.0 67.2 36.4 90.8 16.4 34.4 CSA (mm2) 82 55 168 91 227 41 86 Tongue Extrinsic Muscles Name STY HG MY GH MaxForce (N) 43.6 118 35.4 32.0 CSA (mm2) 109 295 88 80 Table 3.1: Physiological Cross-Sectional Area (CSA) and maximum force generating capability of jaw and tongue muscles. 3.1.2 Tongue model For the dynamic tongue model, we implemented the model by Buchaillard et al. [29] in ArtiSynth, as pictured in Figure 3.2. The original tongue model was based on the anatomy of a single subject using CT data and developed in the ANSYS environment [12] representing the tongue with hexahedral finite elements and hyperelastic properties. Thanks to a collaboration with Buchaillard and colleagues, we were able to obtain data for the 3D tongue mesh and description of the lingual muscular fibers. The mesh and muscle geometry were imported into the ArtiSynth environment, using a large defor- mation FEM framework, hexahedral elements with a density of 1040 kg/m3, and a fifth order incompressible Mooney-Rivlin material with c10 = 1037, c20 = 486, and c01 = c11 = c02 = 0 Pascals. The deviatoric potential energy Ψ̂ of this material is hence Ψ̂ = c10(ĪC − 3) + c20(ĪC − 3)2 (3.1) 46 3.1. Model Creation STY IL MH GGA GH GGP GGM SL STY MH Front Back Sagittal Cut-Away Figure 3.2: Front, back, and sagittal cut-away views of tongue model. At- tachment nodes are also shown for the jaw (front view, red spheres) and the hyoid bone (back view, blue spheres). Muscle groups include the ge- nioglossus (GGA, blue; GGM, green; GGP, red), styloglossus (STY, cyan), geniohyoid (GH, magenta), mylohyoid (MH, orange), hyoglossus, (HG, red), vertical (VERT, green), transverse (TRANS, blue), inferior longitudinal (IL, cyan), and superior longitudinal (SL, magenta) muscles. where ĪC is the first invariant of the deviatoric component on the right Cauchy-Green tensor [22]. The Mooney-Rivlin material was chosen for con- sistency with the reference tongue model and because previous studies have experimentally measured the material parameters from indentation tests on cadaveric tongue and cheek tissue [63]. Incompressibility was implemented using a constraint based approach to maintain the volume of each hexahedral element. Muscles are represented by sets of elements and implemented with node- to-node fiber forces distributed throughout the muscle elements along the principle direction of action. We chose to use a straightforward model for muscle activation with fiber forces directly scaled by input activation, as opposed to the Equilibrium Point Hypothesis (EPH) (λ-model) used in the original tongue model. Our aim was to quantify the coupling between tongue and jaw and not to work on the λ-model motor control assumptions provided by the equilibrium point hypothesis [55]. Tissue stiffening due to muscle activation was also modeled in the same way as Buchaillard and colleagues, i.e. a linear increase of c10 and c20 values ranging between (1037 Pa, 486 Pa) at no activation and (10370 Pa, 4860 Pa) at full muscle activation. As in the 47 3.1. Model Creation original model, each muscle’s force capacity was a function of its CSA (see Table 1 in [29]), with force capacity distributed across fibers weighed by the volume of their surrounding elements. We are currently investigating a more sophisticated muscle modeling approach, as discussed in Section 3.5, using a transversely anisotropic material to more accurately represent muscle tissue mechanics with FEM. 3.1.3 Registration To conform the different morphology of the two models, we adapted skeletal and muscle geometry of jaw-hyoid model to fit CT data (shown in Fig- ure 1.1a) for the subject upon which the Buchaillard tongue model was based. The 3D jaw, skull, and hyoid bone surface meshes were morphed with an energy-based non-linear mesh registration algorithm [30] to a 3D skull surface segmented from CT data. Symmetry was attained by mirroring the left-side of the registered meshes. The inertia of the jaw and hyoid bone were computed from new mesh shapes, assuming uniform density of 3600 kg/m3 and 2000 kg/m3 for the jaw and hyoid bone respectively. Jaw muscle origin and insertions points were adapted with the same non-elastic trans- formation as was applied to the surface meshes and were manually verified to correspond to plausible anatomical landmarks. We removed the point-to- point geniohyoid and mylohyoid muscles from the jaw-hyoid model as these were included in the tongue model. The anterior and posterior digastric muscles were connected to the hyoid bone with the digastric sling modeled as a pulley. 3.1.4 Attachment Hyoid suspension in the reference tongue model was done using eight ver- tical springs with 220 N/m stiffness to connect the hyoid bone to a fixed larynx [29]. Vertically-oriented One Dimensional (1D) springs do not accu- rately represent the 3D stiffness of connective tissue, ligaments, and mus- cles connecting the hyoid bone to the pharynx and larynx within the neck. Therefore, in our coupled jaw-tongue-hyoid model, instead of vertical 1D 48 3.1. Model Creation Front Oblique Sagittal Cut-Away Figure 3.3: Front, oblique and sagittal cut-away views of the dynamically coupled jaw-tongue-hyoid model. springs, we use a linear 6-DOF translational/rotational spring to connect the hyoid bone to a fixed larynx. Lacking sufficient physiological data on the stiffness of the hyoid within the neck, we set the spring stiffness to be consistent with the reference tongue model (8×220 N/m). We are currently developing a dynamic larynx model, which will allow us to more accurately represent the extrinsic laryngeal connective tissue and muscles in order to simulate hyolaryngeal movement within the neck (see Section 3.4.3). Given the limitation of only including passive hyoid stiffness in the model our cur- rent simulation results do not focus in the accuracy of hyoid movement. We couple the dynamics of the jaw, tongue, and hyoid models by defining attachment constraints between the FEM nodes of the tongue and the jaw and hyoid rigid bodies, as described in Section C.1.3. Point to rigid body attachments can be made at arbitrary locations and are not required to be coincident with the rigid-body surface mesh. The attachment points in the model are shown in Figure 3.2. Tongue-jaw attachments include the insertion of genioglossus and geniohyoid onto mandibular geniotubercle and the insertion of mylohyoid along mandibular mylohyoid ridge. Tongue- hyoid attachments include the entire region around the anterior-superior surface of the hyoid bone, including insertions of geniohyoid, mylohyoid, and hyoglossus muscles. The posterior medial surface of tongue is not attached allowing the base of the tongue to move relative to the hyoid bone. The 49 3.2. Simulation Descriptions Jaw Tasks CLR† DI ILP SLP rest — — — — clench 10 — — — open — 15 15 15 hinge-open — 15 — — protrude — — 15 15 right-lateral — — 15? 15? Tongue Tasks CLR† SL TRANS GGP GGM STY retract — — — — — 25 palate 0.5 30 30 60 30 — max-palate 1 100 100 80 30 10 Table 3.2: Percentage muscle activation used in jaw-tongue-hyoid tasks. †Jaw closing muscles (AT, MT, PT, MP, SM, and DM). ?Only left-sided muscles activated. resulting combined model is pictured in Figure 3.3. The soft palate and palatoglossus muscle are not included in the current model as we are not currently interested in soft-palate movement. However, an extended model that includes additional upper-airway structures is under development as described in Section 3.4. 3.2 Simulation Descriptions We chose to simulate a set of tasks similar to those reported for the jaw and tongue models in isolation, including free jaw movements [183], unilateral chewing [71], and tongue movements in speech [29]. All of the tasks, with the exception of unilateral chewing, involve simple piece-wise linear input mus- cle activations so that the passive dynamic coupling effects can be analyzed. Our objective was to analyze the effect of dynamic coupling by using muscle activation to drive the coupled jaw-tongue-hyoid model and observe differ- ences in the movement of the “active” body as well as movement induced on the “passive” body. 50 3.2. Simulation Descriptions 0 100 200 300 400 500 600 700 8000 Time (ms) Ac tiv at ion  (% ) 0 200 4000 Time (ms) Ac tiv at ion  (% ) Figure 3.4: Input muscle activation pattern for jaw tasks (left) and tongue tasks (right). The activation amplitude for each task is given in Table 3.2. 3.2.1 Jaw activated tasks Jaw movement tasks used jaw muscle activation as input and were per- formed both with the jaw-hyoid model alone and with the jaw-tongue-hyoid model. All jaw tasks involved a simple pattern of muscle input (rest, ramp- up, hold, ramp-down, rest) as illustrated in Figure 3.4. The duration of these standardized jaw movements (600 ms) are consistent with an average chewing cycle. The muscle sets and activation amplitudes were chosen to be within physiologically-plausible ranges for jaw movement (10%–15%) and are summarized in Table 3.2. A nominal jaw movement task is rest posture: the equilibrium position of the model with no muscle activation and under downward gravity. In relaxed humans, the jaw typically rests with a 4-6 mm incisal separation. We expect slightly wider gape in jaw-tongue-hyoid model than in the jaw-hyoid model alone, though a majority of the tongue body rests on the hyoid bone, and therefore should result in minimal jaw lowering. Static jaw clenching was simulated with bilateral activation of jaw closing muscles. With no tongue muscle activation we expect the tongue to remain stationary within the mouth during tooth clenching. Opening is simulated by bilateral activation of lateral pterygoids along with the anterior belly of digastric. We chose to simulate a moderate opening gape, with 15% activation in each muscle. Maximum jaw opening in humans is 50 mm on average though both backward head rotation and hyoid posi- tioning become important at wide gape [130] and would complicate the task. We also simulated hinge-like jaw opening with activation of digastric alone. 51 3.2. Simulation Descriptions We expect reduced opening with the jaw-tongue-hyoid model due to passive compression at the floor-of-mouth. The tongue is actively retracted during vowel production, such as /a/, however our opening simulation is with the tongue at rest and therefore we expect the tongue to passively protrude from the mouth during jaw opening [108]. For this reason, the simulation will require contact handling (see Section C.1.4) between the tongue tip and the lower teeth. Protrusion is simulated by bilateral activation of lateral pterygoid mus- cles. The effect of jaw protrusion on the tongue has important implications for OSA as a common therapeutic device is a dental appliance used to ad- vance the jaw and tongue in order to open airway [167]. We expect reduced jaw protrusion in jaw-tongue-hyoid case as compared to the jaw-hyoid model because the lingual elastic connection between the jaw and hyoid should pro- vide some resistance to jaw movement. We also expect forward translation of the base of the tongue. Right laterotrusion is simulated by activation of the left-side lateral pterygoid muscles. We also expect reduced lateral devi- ation in the laterotrusion task with jaw-tongue-hyoid model as compared to the jaw-hyoid model alone. Unilateral chewing involves a complex pattern of jaw muscle activity [124]. We simulated right-sided chewing movement using the same muscle ac- tivation patterns and food bolus that were reported for our original jaw model [71]. An elastic, spherical food bolus (10 mm in diameter) was posi- tioned between the right first molars, which provided resistance during the closing phase of the chewing stroke and collapsed when the applied force exceeded 35 N. Since our adapted jaw-hyoid model has a different bone and muscle geometry, we expect that its chewing movement will be altered, but still plausible, as compared to the original jaw-hyoid model. We also expect that the chewing movement for the jaw-tongue-hyoid model will be signif- icantly altered as the muscle patterns were previously tuned by Hannam et al. [71] to a model without a tongue. 52 3.3. Simulation Results 3.2.2 Tongue activated tasks Tongue movement tasks used tongue muscle forces to move and deform the tongue within the mouth. Tongue tasks involved a ramp-up, hold, and ramp-down pattern of muscle input similar to the jaw tasks, but with faster transitions (50 ms, see Figure 3.4) for consistency with the speed of speech movements. Tongue retraction was simulated by activating styloglossus us- ing the activation trajectory shown in Figure 3.4. We expect that, with the jaw at rest, a retracted tongue posture should induce backward movement of jaw. Tongue-palate contact is an important movement for speech. Tongue tip contact with the anterior hard palate was simulated by activation of superior longitudinal, posterior genioglossus, and transverse muscles as illustrated in Figure 3.8. We stabilized the jaw with low-level (0.5%) activation of jaw closing muscles to maintain a nearly closed jaw posture. We expect that tongue-palate contact will induce a downward movement on the jaw, causing the jaw to open wider. We also performed a maximum tongue-palate pressure simulation by ramping the superior longitudinal and transverse muscles to maximum activation. The ability to generate tongue-to-palate pressure is an important component of healthy swallowing function and we expect that the model’s maximum tongue-palate pressure will be comparable with recorded maximum tongue pressure measurements. The tongue-palate contact simulation required contact handling between the tongue tip and the hard palate surface mesh (see Section C.1.4). 3.3 Simulation Results 3.3.1 Jaw-tongue-hyoid coupling We observed a number of interesting influences of dynamic coupling on the simulated jaw-tongue-hyoid movements. Jaw movements were altered by the presence of the tongue and tongue movements were observed to induce jaw movement. We found a resting jaw posture with 5.6 mm and 6.6 mm incisal gape for the jaw-hyoid and jaw-tongue-hyoid models respectively. 53 3.3. Simulation Results 0 5 10−20 −15 −10 −5 0 Anterior − Posterior (mm) In fe rio r −  S up er io r ( m m ) JTH JH JTH JH Figure 3.5: The jaw-tongue-hyoid model pictured during rest posture (Rest) and at peak jaw opening (Open) with grid spacing of 10 mm. Jaw opening induced passive tongue protrusion such that the tongue tip rested on the lower teeth. Right-most plot shows a lateral view of incisor point movement for opening and hinge-opening with the jaw-hyoid model alone (JH; blue) and the jaw-tongue-hyoid model (JTH; red) with point spacing of 10 ms. The larger incisal gape at rest in the jaw-tongue-hyoid case is due to the added mass of the tongue, but change is small. Also, resting jaw posture is thought to be maintained by low level tonic jaw muscle activity [71] that was not included in our rest posture simulations. Static clenching with the jaw-tongue-hyoid model simulated correctly with the tongue remaining stationary in the mouth. The results of the jaw opening and hinge-opening simulations are shown in Figure 3.5. In both cases the amplitude of jaw opening is reduced in the jaw-tongue-hyoid model due to compression of the lower portion of the tongue between the jaw and hyoid. Figure 3.5 also illustrates the 3D model at the wide gape position showing that jaw opening does indeed cause passive forward protrusion of the tongue such that the tongue tip is resting on the lower teeth. Protrusion and laterotrusion movements are shown in Figure 3.6. The amplitude of protrusion was reduced in the jaw-tongue-hyoid model, likely due to stretching of tongue tissue between the jaw and hyoid, but the am- plitude of lateral deviation was comparable. Interestingly, the tongue also 54 3.3. Simulation Results Protrusion −15 −10 −5 0 −10 −5 0 Anterior − Posterior (mm) In fe rio r −  S up er io r ( m m ) JTH JH Right Lateral Movement −10 −5 0 −10 −5 0 Anterior − Posterior (mm) In fe rio r −  S up er io r ( m m ) 051015 −10 −5 0 Right − Left (mm) In fe rio r −  S up er io r ( m m ) JTH JH JTH JH Right Unilateral Chewing −5 0 5 10−25 −20 −15 −10 −5 0 Anterior − Posterior (mm) In fe rio r −  S up er io r ( m m ) −50510−25 −20 −15 −10 −5 0 Right − Left (mm) In fe rio r −  S up er io r ( m m ) JTH JH JTH JH Figure 3.6: Lateral and frontal views of incisor point movement during jaw- muscle activated simulations for the jaw-hyoid model (JH; blue) and the jaw-tongue-hyoid model (JTH; red) with point spacing of 10 ms. 55 3.3. Simulation Results 0 0.5 1 1.5 2−1.5 −1 −0.5 0 0.5 Anterior − Posterior (mm) In fe rio r −  S up er io r ( m m ) Figure 3.7: The jaw-tongue-hyoid model pictured during rest posture (Rest) and at tongue retraction with styloglossus activation (Retracted) with grid spacing of 10 mm. The right-most panel plots a lateral view of incisor displacement with point spacing of 10 ms. induced significant downward movement of the jaw during protrusion and laterotrusion. The tongue model is pulled forward during jaw protrusion transferring more of its weight from the hyoid bone to the jaw. Right-side chewing movement produced by applying muscle patterns tuned for a different jaw geometry provided a plausible tear-drop shaped incisor movement in the current jaw-hyoid model as shown in Figure 3.6. The movement compares well with one produced by the original jaw-hyoid model (see Figure 2 in [71]). The amplitude of the chewing envelope is re- duced with the jaw-tongue-hyoid model, which is consistent with the jaw opening and laterotrusion simulations. Simulation of tongue retraction in the mouth is pictured in Figure 3.7 (left panels), along with a plot of incisor displacement (right-most panel). Styloglossus activation initially causes an upward and backward movement of the incisor as the condyles move up the articular slope, followed by a back- ward and downward displacement as the tongue retracts farther. The tongue retraction simulation also demonstrates the large range of tongue movement capable with the model and motivates the need for a large-deformation FEM approach. Figure 3.8 shows the mid-sagittal position of jaw, tongue, and hyoid for the tongue-palate contact simulation. Tongue lifting and palate contact 56 3.3. Simulation Results −4 −3 −2 −1 0 3 4 5 6 7 8 Anterior − Posterior (mm) In fe rio r −  S up er io r ( m m ) Figure 3.8: The jaw-tongue-hyoid model pictured during rest posture (Rest) and with tongue lifted into contact with palate (Palate-Contact) with grid spacing of 10 mm. The right-most panel plots a lateral view of incisor point displacement, which starts before tongue-palate contact (× de- notes the beginning of tongue-palate contact) with point spacing of 10 ms. causes a downward jaw movement as expected. The right-most panel plots the incisor point displacement, which starts before tongue-palate contact (as denoted by the × on the plot). The initial downward jaw movement is caused by tongue muscle activation and it increases as force is applied between the tongue and palate. 3.3.2 Comparison with published data The jaw-tongue-hyoid model has been assembled from previously reported reference jaw [71] and tongue models [29]. Qualitative evaluation of the cou- pled model shows similar levels of force and range of movement as exhibited by each individual model. Therefore, the new implementation of these mod- els in the ArtiSynth framework compares well with the previously published versions. As a step toward validation of the jaw-tongue-hyoid model we have made preliminary comparisons to published data on tongue kinematics and forces. We used the tongue retraction simulation as a means to compare tongue velocity generated in the model with measured tongue movement during speech. Payan and Perrier [144] recorded the movement of surface tongue points with EMA during tongue retraction in a [y-o] speech utterance for 57 3.3. Simulation Results 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.90 50 100 150 200 Time (s) Velocity (mm/s) 0 50 100 150 2000 50 100 150 200 Time (ms) Ba ck wa rd  V elo cit y ( m m /s) Data Model Figure 3.9: Backward velocity of a point on the upper tongue surface recorded with an electromagnetic tracking system an during [y-o] speech utterance (Data; reproduced from [144] page 15, figure 9a) and simulated with the tongue model (Model; point shown as cyan sphere in Figure 3.7). the speaker upon which the tongue model was based. Figure 3.9 shows the recorded velocity profile reproduced from Payan and Perrier [144] as well as the simulated anterior-posterior velocity of one node on the tongue’s upper surface. Simulating tongue retraction with a physiological level of styloglossus activation (25%) generates a peak velocity that compares well with the recorded velocity. Both velocity profiles are bell shaped and the peak velocity values are 210 mm/s and 215 mm/s for the recorded and sim- ulated movements respectively. The velocity profile in the simulated case is narrower and also slightly asymmetric (the right side of the bell-shaped curve tapers off more slowly). The asymmetry may be due to the fact that increasing velocity is created by styloglossus activation while decreasing ve- locity is caused only by passive forces slowing down the tongue. In recorded movement antagonist muscles may be recruited to slow down the tongue more quickly. Also, our simulated tongue retraction was generated with a linearly ramped styloglossus activation, whereas the real vowel movement may be caused by a exponential increase in styloglossus activation, which would widen the velocity profile. The damping parameters in the tongue model will affect the tongue velocity profile and could also contribute to the discrepancy in velocity profile shape. We plan to investigate these aspects 58 3.3. Simulation Results of speech-like tongue movement simulations further in future studies. We used the maximum tongue-palate pressure as a metric to evaluate whether or not the parameters for tongue muscle forces in the model are within a plausible range. The maximum force-generating capability of the model’s tongue muscles are proportional to their CSAs (see Table 3.1 and [29, 205]). We also use a CSA-to-maximum-force constant of 40 N/cm2 [146]. Direct measurement of muscle force in vivo is not possible; therefore we rely on external force measurements as an indirect means to evaluate the resultant force generation capability with the jaw-tongue-hyoid model. In particular, Utanohara et al. [201] used a balloon-type disposable oral probe to measure tongue pressure by having subjects compress it onto the palate with maximum voluntary effort. The authors recorded pressures for a large subject pool (850 subjects) and report 40.4 ± 9.8 kPa (mean ± standard deviation) maximum tongue pressures for subjects between forty and forty- nine years of age. Pressure between the tongue-palate contact in the model was calculated by dividing the magnitude of the contact constraints by the area of the contact contours (see Section C.1.4 for discussion of collision detection and handling in ArtiSynth). Maximum tongue-palate pressure simulated with the model was 38.2 kPa, which compares well with the mean value of 40.4 kPa reported by Utanohara et al. [201]. The range of measured pressures (± 9.8 kPa) could also be used to assess the sensitivity of the tongue muscle parameters. For example, applying plausible ranges of tongue muscle CSAs or CSA-to-maximum-force constant in the model could be used to determine if they achieve the measured range of tongue-palate pressures. We plan to move to a more sophisticated model of tongue muscle ac- tuation by incorporating muscle fiber forces into the FEM material with a transversely isotropic constitutive law as discussed below in Section 3.6. This will require additional parameters to describe how the force generating capacity of tongue muscles vary with muscle length and shortening veloc- ity (e.g. a Hill-type muscle formulation, see Section 2.3.3). The maximum voluntary tongue-palate pressure validation will then need to be revisited for the new muscle formulation and additional experimental data may be needed given the additional model parameters. 59 3.3. Simulation Results 0 0.2 0.4 0.6 0.8 10 5 10 15  D isp lac em en t N or m  (m m ) Lower Incisor 0 0.2 0.4 0.6 0.8 10 2 4 6 8  D isp lac em en t N or m  (m m ) Lower Incisor 0 0.2 0.4 0.6 0.8 10 5 10 15  D isp lac em en t N or m  (m m ) Tongue Tip Time (s) 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 14  D isp lac em en t N or m  (m m ) Tongue Tip Time (s) Jaw Opening Tongue-Palate Contact Figure 3.10: Comparison of results for 1 ms (solid lines) and 5 ms (dotted lines) integration steps for the tongue tip (top) and lower incisor (bottom) during jaw opening (left) and tongue-palate contact (right) simulations. 3.3.3 Integration Error The simulations reported above were computed using a second-order New- mark integrator (see Section C.1) with a 5 ms integration step size. In order to assess the numerical error we compared with results computed with a 1 ms integration step. Differences between 5 ms and 1 ms integration steps were found to be small. Figure 3.10 plots lower incisor point and tongue tip displacements for the jaw opening and tongue-palate contact simula- tions computed with 5 ms and 1 ms step sizes. Error was computed as the difference between the displacement trajectories relative to the maximum displacement. The jaw opening simulation showed a very small difference between the 5 ms and 1 ms integration step conditions: the incisor dis- placement error had average and maximum values of 0.5% and 1.4% and the tongue tip displacement error had average and maximum values of 0.7% and 1.6%. Larger discrepancies were found in the tongue-palate contact task, which involves significant contact situations: the incisor displacement error 60 3.4. Additional Orofacial Sub-Models Front Sagittal Cutaway Figure 3.11: Frontal and sagittal cut-away views of an integrated FEM face model integrated with the jaw-tongue-hyoid model showing tongue protru- sion through the lips. had average and maximum values of 1.1% and 3.8% and the tongue tip dis- placement error had average and maximum values of 2.0% and 10.0%. The integration error in the contact simulation is likely due to the discontinuous nature of contact as well as the spatial discretization of the palate/tongue surface meshes on which collisions are detected and responses are generated. Even in the worst case, though, the integration error remains small. 3.4 Additional Orofacial Sub-Models We are also expanding our jaw-tongue-hyoid model to include other anatom- ical structures toward a complete model of orofacial and upper airway biome- chanics. Our preliminary models include the face, soft-palate, and larynx. 3.4.1 Face model We have registered and integrated a muscle-activated FEM model of the face [133], as shown in Figure 3.11. The face model has been adapted to the same subject using CT data in a similar fashion as the skull meshes for the jaw model [30]. The model includes three layers of tissue and line- 61 3.4. Additional Orofacial Sub-Models Front Side Back Figure 3.12: Oblique front, side and oblique back views of the soft-palate model registered to the jaw-tongue-hyoid model. The mandible is cut-away in the side and oblique back views for clarity. based muscle models representing the primary facial muscle groups. We have performed preliminary simulations with the face-jaw-tongue-hyoid model activating the jaw and tongue muscles with the face at rest in order to investigate the effect of passive facial soft-tissue forces on jaw movement. This model has a number of potential high-impact applications. Face and lip motion are important to audio-visual communication. A large amount of face deformation and lip opening is caused by jaw motion and previous face models have not included muscle-activated jaw structures. Also, facial tissues are important in mastication because the cheek (buccinator muscle), together with the tongue, is used to form and manipulate the food bolus during chewing. The integrated model will allow us to simulate chewing with a more realistic, free-floating food bolus. Finally, face modeling has a number of biomedical applications, most notably predicting aesthetic out- comes of alterations to underlying bone structure, such as in orthognathic surgery [36]. Our integrated orofacial model allows for prediction of facial appearance during dynamic jaw and face movements, as opposed to previous models that were limited to static facial aesthetics. 62 3.4. Additional Orofacial Sub-Models 3.4.2 Soft-palate model We have also developed a preliminary FEM model of the soft-palate, as pictured in Figure 3.12. The model geometry was build by manual seg- mentation of MRI data and includes line-based muscle fibers. The model has been approximately registered to the jaw-tongue-hyoid model. On-going work includes adapting and registering the geometry to fit the subject upon which the jaw-tongue-hyoid model was built, which will be informed by existing datasets on that subject’s soft-palate shape [170]. Also, the soft- palate muscles need to be connected to surrounding structures, including the palate elevator muscles to the skull (which may require more detailed bone structure at the base of the skull to determine muscle paths and insertion sites, see Appendix B), the palatoglossus muscle (anterior palatal arch) to the tongue, and the palatopharyngeus muscle (posterior palatal arch) to the pharynx (which has not yet been modeled). The soft-palate is similar in structure to the tongue, but on a smaller scale. It is functionally important in swallowing, as it seals off the nasophar- ynx from the oral cavity, and in respiration, as it is a common site of airway obstruction in OSA. Therefore, the integrated tongue-palate-pharynx model has a number of important potential biomedical applications. 3.4.3 Hyoid-larynx model We have developed a model of the larynx, including the cricoid and thyroid cartilages and extrinsic laryngeal muscles, as pictured in Figure 3.13a. We are interested in the gross movement of the larynx within the neck as it anchors the hyoid bone and tongue through passive connective tissue and the extrinsic laryngeal muscles. Elevation of the hyoid and larynx during swallowing has been associated with airway protection and down-folding of the epiglottis [202]. We are building the hyolarynx model based on a dataset of hyoid and larynx movements, induced by intramuscular stimulation of tongue and la- ryngeal muscles, recorded with VF [31]. One frame of VF data is shown in Figure 3.13b with manually selected landmarks. The mechanical properties 63 3.5. Discussion (a) (b) Figure 3.13: The dynamic hyoid-larynx model (a) and one frame video fluo- roscopic data showing manually selected landmarks (b). Video fluoroscopy is used to record movement of the hyoid bone and larynx induced by intra- muscular stimulation of extrinsic laryngeal and tongue muscles. of the extrinsic laryngeal muscles are largely unknown. Therefore, we are using the muscle stimulation data to estimate the effect of individual muscle forces on hyolaryngeal movement in order to select plausible ranges of force parameters for the model. Clearly, the tongue plays a significant role in hy- oid and larynx movements and therefore the integrated tongue-hyoid-larynx model is required to accurately assess muscle function. We plan to apply the hyolaryngeal model to analyze the mechanics of hyoid elevation and air- way protection. Also, electrical muscle stimulation is being used clinically in dysphagia therapy [109] and biomechanical analysis could help inform strategies for effective muscle stimulation during swallowing. 3.5 Discussion The simulations reported here are a proof of concept demonstrating the ef- fectiveness of our hard-soft tissue simulation and motivating the need to include dynamic coupling in simulations of jaw-tongue-hyoid movements. The reported simulations demonstrate that a wide range of movement, large forces, and large tissue deformations are possible within the current simula- tion framework. We have provided preliminary qualitative comparisons of 64 3.5. Discussion tongue velocity and pressure to illustrate that the model behaves within a plausible range of human movement and force production. Model valida- tion is important and highly dependent on a model’s intended use. Recent reviews have provided high-level guidelines for biomechanical model valida- tion [7], with a focus on modeling domains where detailed experimental data are available and direct validation metrics are applicable. As discussed in Section 2.1 and stated by Hannam [69], acquiring detailed experimental data for upper-airway function with sufficient quality for direct model validation is a significant challenge. Indeed, one of the main motivations for develop- ing upper airway models is the difficulty of experimentally assessing upper- airway biomechanics. We believe that indirect validation metrics, such as the incisor movement, tongue velocity and pressure comparisons made in Section 3.3.2, are the only feasible approaches to validating biomechanical models of the upper airway. Further indirect validation of the jaw-tongue- hyoid model is planned, as additional quantitative data were recorded from the subject upon which the model is based [14]. The passive suspension of the hyoid bone to a fixed larynx in our model is a limitation. The amplitude of hyoid displacement in our simulations is smaller than reported in experimental recordings of eating and speaking movements [79]. This may be attributed to the stiffness of the 6-DOF hy- olaryngeal spring in the model, but is also likely due to the fixed larynx and lack of hyoid depressor muscles. The restricted jaw gape and forward tongue protrusion during jaw opening may be due in part to the reduced hyoid movement, though backward rotation of the head has also been shown to be important to wide jaw opening [99, 130]. We are currently developing a dynamic larynx model, as discussed in Section 3.4.3, including extrinsic la- ryngeal connective tissue and muscles, which will allow us to better analyze hyoid movement in our coupled model. The results of the chewing simulation show that the muscle patterns of [71] are applicable to a different skull morphology, as they produced a very similar chewing pattern, suggesting that they are not overly sensitive to skeletal or muscle geometries. The results also show that the addition of passive tongue tissue has a significant effect on free jaw movements and the 65 3.5. Discussion chewing movement, due to both the passive elastic connection between the jaw and hyoid made by the tongue, especially in compression, e.g. reducing jaw opening, as well as the additional mass of the tongue body, particularly in jaw protrusion. Biomechanical simulations of the type reported here are challenging to create. However, the interactivity afforded by computational performance in ArtiSynth allows for reasonably fast refinement and exploration of the model’s capability. The jaw-tongue-hyoid model described above has two free rigid bodies and 946 FEM nodes for a total of 2505 degrees of freedom. With respect to (Equation C.6), the addition of point-based tongue attach- ments, incompressibility, and jaw joints result in an M̂ that is 2505× 2505 and a G that is typically 2505 × 740 (varying somewhat depending on the number of FEM contacts). In addition, a few unilateral constraints are used to implement bite contact. Solution times for (Equation C.6) using the method described in Section C.1.2 vary from around 130 ms to 200 ms (de- pending on whether unilateral constraints are in play) on a 2.6 GHz Core 2 Duo processor. Overall solution time (including collision detection and all the steps of Section C.1.5) for a 600 ms jaw opening task with a time step size of 5 ms is around 40 seconds and a 400 ms tongue retraction task with a time step of 10 ms is around 20 seconds. Much of this involves Java code that could be significantly optimized. This improves on the computa- tion time reported in [29], where a 100 ms task for the same FEM tongue model (with jaw/palate contact) required 40 minutes of computing time on a similar computer using ANSYS. Complex movements require precise coordination among a large number of muscle input degrees-of-freedom making trial-and-error tuning of muscle inputs to generate simulations tedious and likely over-fitted to a particular model’s geometry. We believe that optimization-based inverse dynamics approaches are a promising direction in this regard, which we discuss further in the following chapter. 66 3.6. Directions 3.6 Directions Our preliminary results also point to a few promising directions. We plan to investigate what changes in muscle patterns are required to improve the chewing stroke in the jaw-tongue-hyoid model, e.g. increasing jaw opener muscle activation in order to attain a larger jaw opening during the simulated chewing cycle. It is noteworthy that the original muscle patterns reported by Hannam et al. were reported as being low in amplitude, therefore muscle activation amplitudes could be increased and remain plausible. We also plan to investigate activating the tongue muscles in concert with the jaw muscles to simulate tongue movements [78] and palate contact pressure patterns [138] during the chewing cycle, which would add a dynamically changing inertia as opposed to the current passive tongue mass. The tongue tissue is currently modeled as isotropic with point-to-point actuators embedded within the material to model anisotropic muscle fiber forces. In the current model we increase the stiffness of elements associated with muscle activation as an approximation of skeletal muscle stiffening. However, we are currently incorporating a transverse-isotropic material [212] into ArtiSynth, which will provide a more realistic representation of skele- tal muscle mechanics. The tongue muscles are particularly challenging to model because multiple muscle groups, with different principal fiber direc- tions, converge and interdigitate within the tongue body. For this reason we are investigating a formulation allowing for the superposition of multiple transverse-isotropic materials with different principal directions. 3.7 Summary To summarize, we have created a 3D model of coupled jaw-tongue-hyoid biomechanics to advance the state-of-art in orofacial modeling. Previously proposed models focused on isolated subcomponents, neglecting the coupling forces from surrounding and adjacent anatomical structures. Our contribu- tion, by creating an integrated model, is the ability to assess the assumption in previously proposed isolated jaw and tongue models that coupling effects 67 3.7. Summary are small. The model was built by adapting and dynamically coupling a reference multi-body jaw/hyoid model and a reference FEM tongue model so that the muscle forces of the tongue imparted forces on the jaw/hyoid and vice- versa. With our coupled model we were able to evaluate the significance of mechanical coupling. Simulations of isolated muscle activations showed that the presence of passive tongue tissue reduced jaw movement and active tongue muscle forces induced jaw movement. Simulations also demonstrated that the model behaved within a plausible range of motion for chewing and speech production. As a step toward validation of the model, we compared simulated tongue velocity and tongue-palate pressure to recorded human measurements and found consistent values for each. Our jaw-tongue-hyoid modeling efforts also served to verify that our underlying simulation plat- form, ArtiSynth, is sufficiently accurate and robust for the challenges of orofacial modeling, namely, large tissue deformations, large muscle forces, and hard/soft tissue coupling and contact. We are currently expanding the jaw-tongue-hyoid model to include deformable face, soft-palate, and lar- ynx models and have created preliminary simulations of integrated face-jaw- tongue-hyoid movements. In order to simulate complex movements we need to specify coordinated, time-varying muscle activations as input to the forward dynamics model. In the next chapter, we describe methods for automatically generating muscle activations to realize specific target outputs. 68 Chapter 4 Inverse Simulation Methods In the previous chapter, we described a model of jaw-tongue-hyoid dynamics and its evaluation with motor tasks generated by simple piece-wise linear muscle activation inputs. A significant challenge in the application of biome- chanical models to the analysis of motor tasks is the specification of muscle inputs. In this chapter, we describe inverse methods for automatically pre- dicting complex patterns of muscle activation input required to drive a model through a prescribed trajectory of movement and/or forces. Predicting muscle activity from kinematic information is an inverse prob- lem because anatomical systems have more muscles than kinematic DOF and are therefore redundant. For example, the human mandibular system has thirty muscle groups (some of which have sub-regions that are independently activated) and only six kinematic DOF. Motor redundancy also exists in the tongue, though it is less obvious than in articulated skeletal systems such as the jaw or limbs. The tongue is a deformable soft-tissue structure and therefore has many kinematic DOF, however kinematic analysis has shown that the tongue changes shape within a low dimensional motion space. For example, in speech movement 3D tongue shape can be characterized by a small number of DOF. The tongue has ten distinct muscle groups, some of which are likely differentially activated, and therefore can be considered re- dundant within the reduced-space of kinematics observed in physiologically relevant tasks. Another challenge with specifying muscle activity as input to forward dynamics models is that muscle activity is difficult to measure. EMG in- volves transducing electrical signals associated with muscle activation; how- ever EMG can be difficult to record for small, deep muscles in the head and neck, and the relationship between EMG and muscle force is complex for dy- 69 Chapter 4. Inverse Simulation Methods namic movements [173]. One approach for validating biomechanical models is to apply recorded EMG patterns as input and compare simulated kine- matics and forces with measurements. However, due to the challenges with EMG recording and the inability to record from all muscle simultaneously, it is desirable to do the opposite: use kinematic and force measurements as input to the model and compare predicted muscle forces to recorded EMG. As described in Section 2.4, recently reported trajectory-tracking tech- niques for inverse-dynamics simulation have shown promising results with limb models [196], complex muscle and tendon models [189], and quasi-static FEM models [179]. These techniques use per-timestep static optimization, as opposed to optimizing over the full time-varying trajectory, but they do incorporate model mechanics and are computationally efficient. In this chapter, we extend the previously proposed trajectory-tracking inverse approach for use with dynamic FEM models, such as our tongue model, as well as with coupled rigid-deformable models, such as our jaw- tongue-hyoid model described in the previous chapter. We also formulate new target parameters, in addition to kinematic targets, that can be used as input to the inverse simulation in order to select desirable muscle activation patterns to overcome motor redundancy. The new target parameters include constraint force magnitudes and stiffness. Constraint force targets can be used to generate desired reaction forces in the system, such as a target bite force during jaw clenching simulations. The stiffness target can be used to control co-activation of antagonist muscles. These new target parameters allow for multiple modalities of target data, such as both motion and contact force recordings, and permit systematic analysis of potential muscle activa- tion patterns. Section 4.1 describes the mathematical formulation of the inverse solver implemented in ArtiSynth. Section 4.2 demonstrates prop- erties of the inverse toolset using simplified canonical models. Section 4.3 reports results of the toolset to more complex jaw and tongue models. Fi- nally, Section 4.4 and Section 4.5 discuss the implications of inverse modeling techniques and directions for further investigation. 70 4.1. Inverse Solver Formulation 4.1 Inverse Solver Formulation The inverse-dynamics tracking algorithm was implemented within the Ar- tiSynth biomechanics simulation toolkit, which is described in more detail in Appendix C. ArtiSynth is designed to simulate the dynamics of hard and soft tissue structures using coupled rigid-body and FEM models. Within ArtiSynth, a mechanical system consists of an assembly of rigid bodies and particles (which include FEM nodes). Let q, u, and f denote the composite position, velocity and force vectors for these components. Following from previous work [179, 189] we let f be partitioned into f = fp+ fa, where fp are the passive forces arising from muscle stretch, ligaments, and scar tissue, and fa are the active forces arising from muscle activation. The system’s dynamic behavior is then determined by Newton’s second law, Mu̇ = f(q,u, t) = fp(q,u) + fa(q,u,a(t)) (4.1) where M is a block diagonal mass matrix. We use Hill-type muscle models that are linear in activation and non- linearly dependent on the length and shortening velocity, so that fa = Λ(q,u)a (4.2) where a is a vector of activation levels bounded between 0− 100% for each muscle. The matrix Λ relates muscle activations to system forces and can be determined either analytically or numerically; we currently use an ana- lytic formulation. We neglect the calcium-dependent activation dynamics of muscle tissue, which is typically modeled as a first-order low-pass filter [216]. The consequence of this assumption is that predicted muscle forces could change faster than physiologically possible, however this was not found to occur for target movements with physiologically plausible velocities. The hu- man tracking measurements place constraints on velocities that get reflected back into physiologically realistic muscle activation values. Situations could arise where muscle activations rapidly switch between agonist muscles, as mentioned by [189] who found that a damping regularizer was needed to 71 4.1. Inverse Solver Formulation reduce such oscillations. The inclusion of muscle activation dynamics in the model would also limit such oscillations in predicted muscle activations. The mechanical system may also contain bilateral and unilateral con- straints; the former include articulating joints between rigid bodies and FEM incompressibility, while the latter include contact conditions and joint limits. Unilateral constraints are not considered here, but we do utilize bi- lateral constraints, which take the form of linear equality constraints on the velocity: G(q)u = 0. (4.3) Differentiating this leads also to acceleration constraints G(q)u̇ = g, g ≡ −Ġu. (4.4) For example, to constrain one point of a rigid-body to a planar surface a constraint is formed to prevent translation of the point normal to the plane, i.e. G = (nx, ny, nz, 0, 0, 0), where (nx, ny, nz) is the normal vector of the plane. Constraints are enforced by forces applied to GT , so that Equation 4.1 becomes Mu̇ = fp(q,u) + Λ(q,u)a + GT (q)λ (4.5) where λ are Lagrange multipliers giving the magnitudes of the constraint reaction forces. Solving the system dynamics involves integrating Equation 4.5 forward in time. At present, this is done using a first order integrator that is semi- implicit with respect to the passive forces fp (which are often stiff). Letting h equal the time step, and using a superscript to denote values at step i, this leads to Muk+1 = Muk + hfk+1p (q,u) + hΛ k(q,u)a + GkT (q)λ (4.6) 72 4.1. Inverse Solver Formulation where λ now denotes constraint impulses. fk+1p is approximated using fk+1p ≈ fkp + ∂fp ∂q ∆q + ∂fp ∂u ∆u = fkp + ∂fp ∂q huk+1 + ∂fp ∂u (uk+1 − uk). Combining this with Equation 4.6 and Equation 4.4 leads to the system( M̂ −GT G 0 )( uk+1 λ ) = ( Muk + hf̂p + hΛa g ) (4.7) where M̂ ≡ ( M− h∂fp ∂u − h2∂fp ∂q ) , f̂p ≡ fkp − ∂fp ∂u uk (4.8) are the mass matrix and force term augmented with Jacobian terms required for the implicit solve. Unlike M, M̂ is neither block diagonal nor symmetric positive definite, but it is sparse and symmetric. Solving for λ in Equation 4.7 we find: λ = (GM̂−1GT )−1g −Q (k + hΛa) (4.9) where Q ≡ (GM̂−1GT )−1GM̂−1 and k ≡ Muk + hf̂p. Back substituting for λ and solving for uk+1 in Equation 4.7 yields uk+1 = QTg + M̂−1Pf (k + hΛa) (4.10) where Pf ≡ (I − GTQ) is a matrix that projects forces into the range compatible with the constraints G, and QT = M̂−1GT (GM̂−1GT )−1 (by the symmetry of M̂). Equations Equation 4.9 and Equation 4.10 relate muscle activations to future constraint forces and velocities and can be used to formulate an optimization over muscle activations with a cost function that includes desired movement and constraint force goals. Movement goal Movement is the traditional goal for inverse dynamics simulation and most recently has been used in rigid-body [189] and quasi-static FEM [179] mod- 73 4.1. Inverse Solver Formulation els. Our contribution includes extending this formulation to dynamic FEM models. The movement goal of the algorithm is given as a target veloc- ity trajectory u∗, and we desire to find muscle activations that minimize 1 2‖u∗ − uk+1‖2. Substituting for uk+1 from Equation 4.10 the optimization term for the movement target can be expressed as a quadratic form in a: φm(a) = 1 2 ‖ū−Hma‖2 (4.11) where ū ≡ u∗ −QTg − M̂−1Pfk and Hm ≡ hM̂−1PfΛ. For a rigid body the target movement can be specified as either: the full 6D position and orientation of the body, the 3D position of a single point on the body (in which case the body’s motion is partly unconstrained), or the 3D position of multiple points on the body (in which case a best least-squares fit to the points is used), as discussed below in Section 4.2.2. Our contribution extends the formulation to work with FEM models, as discussed below in Section 4.2.3. Target velocities can be computed from the target position trajectory at each timestep providing online correction to position errors. Constraint force goal Our contribution also includes extending the inverse dynamics formulation with new goals: constraint forces and stiffnesses. The constraint force goal is given as target values for Lagrange multipliers, λ, which are the magnitudes of the constraint reaction impulses. Given target constraint forces ξ, the corresponding impulses are hξ, and so we desire to find muscle activations that minimize 12‖hξ−λ‖2, which leads to a second term in the optimization cost function which is also a quadratic form in a: φc(a) = 1 2 ‖λ̄−Hca‖2 (4.12) where λ̄ ≡ hξ − (GM̂−1GT )−1g + Qk and Hc ≡ −hQΛ. We can selec- tively include a subset of the constraints into the optimization term. The constraint-force goal is used in Section 5.3 to simulate clenching with the jaw model. 74 4.1. Inverse Solver Formulation Stiffness goal The stiffness goal is given as a desired stiffness in each DOF. It is formulated by extending a model of single joint stiffness control proposed by Hogan [81] (see Section 2.4 for more details), which shows that torque and stiffness about a joint can be independently controlled by the difference and sum of muscle activations. System forces arising from muscle activation is given by fa = Λ(q,u)a (Equation 4.2), whereas stiffness due to co-activation is given by: Ka = Hka (4.13) where Hk ≡ abs(Λ(q,u)), the element-wise absolute value of the Λ matrix. Given the target stiffness values, K∗, we desire to find muscle activations that minimize 12‖K∗−Ka‖2, which leads to a third term in the optimization cost function: φK(a) = 1 2 ‖K∗ −Hka‖2 (4.14) The control of muscle co-activation using the stiffness goal is demonstrated below in Section 4.2.1. It is important to note that intrinsic muscle stiffness increases with mus- cle activation and may be more dominant than co-activation induced stiff- ness in certain biological systems or motor situations. The force-position jacobian, ∂f/∂q, which we calculate for implicit-integration, is a localized system stiffness matrix and could be exploited to control stiffness arising both from individual muscle activation and co-activation of multiple mus- cles. Also, stretch reflexes are currently not included in the biomechanical models we have considered, and reflex feedback gains are used in biological systems to increase effective stiffness. Including intrinsic muscle stiffness and reflex feedback gains in our formulation would result in less co-activation to achieve the same stiffness goal. 75 4.1. Inverse Solver Formulation Per-timestep static optimization Given the above target terms (motion, constraint forces, and/or stiffness) the inverse-dynamics algorithm solves a static optimization problem at each integration timestep of the forward dynamics simulation in order to track the target trajectory. Static optimization allows for more efficient compu- tation because the system dynamics are linearized at each timestep. As discussed in Section 2.4, dynamic optimization or optimal control formula- tions provide an optimization over the full time-varying trajectory. While the per-timestep static optimization may lead to sub-optimal muscle ac- tivations as compared to an optimal control formulation, the technique is significantly more computationally efficient. This is particularly important for coupled rigid-deformable FEM models that have non-trivial forward dy- namics solutions (see Section C.1.2 for a discussion of expected complexity). The per-timestep optimization problem is underdetermined for a biome- chanical system with redundant muscle activations. We include a weighted `2-norm regularization term, 12a TW−1a, where W is a diagonal matrix of muscle CSA, in order to select the most efficient set of activations [4]. Other regularization terms may be used within our formulation, such as an `1-norm term, ‖a‖1, to select the smallest set of non-zero activations as a sparsity constraint or a damping term, 12‖a− ak−1‖2, to enforce smooth activations. Combining the movement and constraint force goals, regularization, and muscle activations bounds, we arrive at the complete optimization problem, which takes the form of a quadratic program: min a wmφm(a) + wcφc(a) + wkφK(a) + wa 2 aTW−1a subject to 0 ≤ a ≤ 1 (4.15) where wm, wc, wk, and wa are weights are used to trade-off between cost terms. For all simulations reported in this chapter the weights were set to wm = 1, wc = 1, and wa = 0.1 so that minimizing tracking error was preferred over small activations. The inverse dynamics optimization is solved at each timestep to pro- 76 4.2. Analysis with Canonical Models vide muscle activations to the forward dynamics simulation. Biomechanics models, especially those that include FEM models, can have many dynamic components resulting in a large, sparse KKT system in Equation 4.7 (see Section C.1 for more details). We use the KKT system solver in ArtiSynth to compute ū and λ̄ as well as Hm and Hc, which are formed by solving the system for each column of Λ. The resulting quadratic program is dense but tends to be small since its dimension is the size of a, i.e. the number of activations being solved for. The quadratic program is also convex, which means it can be solved as a linear complementarity problem, which is done using an implementation of the Cottle-Dantzig algorithm ([38]) contained in ArtiSynth. 4.2 Analysis with Canonical Models Here we evaluate and analyze the inverse-dynamics trajectory tracking algo- rithm proposed in the previous section with a number of canonical models, including a point, rigid-body, and deformable-body models. These simula- tions verify the correctness of the implementation and are used to illustrate properties of the algorithm, which are more readily apparent in these sim- plified cases. The canonical models serve as a basis for the application of the inverse toolset to more complex anatomical models described in Section 4.3. 4.2.1 Point inverse The canonical point model has 2-DOF and includes 16 muscle actuators, arranged radially in the x-y plane, as pictured in Figure 4.1. The cardinal direction muscles (N, S, E, W) have two times larger CSA than the ordinal direction muscles (NE, SE, SW, NW). The kinematic target was specified as a smooth displacement in the north-east direction from the center. Muscle redundancy The point model serves as an illustrative example of motor redundancy and is used to demonstrate the effect of different regu- larization terms and stiffness goals on the output muscle activations. Muscle 77 4.2. Analysis with Canonical Models t = 0s t = 0.25s t = 0.5s Figure 4.1: The canonical 2-DOF point model with 16 muscles shown as red lines. Line thickness denotes muscle CSA. Cardinal direction muscles (N, S, E, W) have two times larger CSA than ordinal direction muscles (NE, SE, SW, NW). The panels picture the inverse simulation of a north-east displacement with the target point shown as a cyan sphere. `2-Norm CSA-`2-Norm `1-Norm 0 1 N 0 1 NNE 0 1 NE Ac tiv at ion  (% ) 0 1 ENE 10 1 E Time (s) 0 1 N 0 1 NNE 0 1 NE Ac tiv at ion  (% ) 0 1 ENE 10 1 E Time (s) 0 1 N 0 1 NNE 0 1 NE Ac tiv at ion  (% ) 0 1 ENE 10 1 E Time (s) Figure 4.2: Muscle activations predicted for the point model simulation with different regularization terms, including an unweighted `2-norm (12a Ta), a muscle CSA-weighted `2-norm (12a TW−1a), and an `1-norm (‖a‖1), where a the muscle activations vector and W is a diagonal matrix of muscle CSA. 78 4.2. Analysis with Canonical Models Low Stiffness High Stiffness 0 1 N 0 1 NE 0 1 E 0 1 SE 0 1 S Ac tiv at ion  (% ) 0 1 SW 0 1 W 10 1 NW Time (s) 0 1 N 0 1 NE 0 1 E 0 1 SE 0 1 S Ac tiv at ion  (% ) 0 1 SW 0 1 W 10 1 NW Time (s) 0 5 100 0.1 0.2 0.3 0.4 0.5 Force Disturbance Magnitude Av er ag e Tr ac kin g Er ro r K*=0 K*=5 K*=10 K*=20 (a) (b) (c) Figure 4.3: Muscle activations predicted for the point model simulation with different stiffness goals, including minimum stiffness (a) and high stiffness (b). Increasing desired stiffness reduces tracking error for a random force disturbance applied to the point (c). activations predicted to drive the point model through the north-east dis- placement trajectory are plotted in Figure 4.2 for different conditions. In all cases the kinematic trajectory was exactly reproduced. The cardinal direc- tion muscles were activated preferentially with CSA-weighted regularization (Figure 4.2b) as compared to unweighted regularization (Figure 4.2a). We also evaluated different stiffness goals. A minimum stiffness goal resulted in zero co-activation (Figure 4.3a), whereas a high stiffness goal produced co-activations that change dynamically throughout the task (Fig- ure 4.3b). Figure 4.3c illustrates that stiffness can be used to regulate force disturbances. We performed a series of simulations with a random force dis- turbance applied to the point with increasing magnitude. The disturbance simulations were performed for different stiffness goals and we found that in- creased stiffness reduced the tracking error induced by the force disturbance, as shown in Figure 4.3c. 79 4.2. Analysis with Canonical Models Unconstrained Point-to-Plane Constraint t = 0s t = 0.25s t = 0.5s Figure 4.4: The canonical 6-DOF rigid block model with 36 muscles shown as red lines. The panels picture inverse simulation of an oblique displacement target (shown as the cyan wire-frame block) for an unconstrained model (upper panels) and a model with the lower back corner (magenta point) constrained to a planar surface (lower panels). 80 4.2. Analysis with Canonical Models 4.2.2 Rigid-body inverse The canonical rigid body model is a block with 6-DOF and includes three sets of orthogonally arranged muscles (36 muscles total) for full controllabil- ity of position and orientation. The rigid body model is used to demonstrate trajectory tracking under dynamic constraints. Constrained dynamics The inverse algorithm is formulated to track a kinematic trajectory in the presence of constraints on the model’s dynamics. Figure 4.4 illustrates an inverse simulation of translation with the canonical rigid-body model. The upper panels show the block following a 6-DOF translation trajectory with no rotation, as shown by the cyan wireframe mesh. The lower panels show the same simulation, but for a model with a point-to-plane constraint on the back-bottom corner of the block. In this case, the simulation computes muscle activations to make the block follow the trajectory as best as possible while maintaining the dynamic constraint, and the block slides along the plane and rotates to minimize the distance to the target position. We can also simulate a desired constraint force that causes increased muscle activation in order to apply force against the point- to-plane constraint during the movement trajectory. 4.2.3 Deformable-body inverse The canonical deformable body model is a FEM beam with orthogonally arranged muscle fibers throughout the FEM mesh. The mesh has 3x3x7 hexahedral elements, 128 nodes (384 DOF), and 276 muscle fibers. Nodes on the back end of the beam are fixed in space. The deformable body model is used to illustrate the effect of muscle grouping and muscular-hydrostat type motor recruitment patterns as proposed by [96]. Muscle grouping Muscle fibers within a model can be grouped based a priori knowledge about the spatial extent of muscle fibers for particular muscle groups. In the beam model we compare the effect of different muscle grouping in an elongation 81 4.2. Analysis with Canonical Models Two Anteroposterior Divisions A P Seven Anteroposterior Divisions A PMMPMMAAM P t = 0s t = 0.25s t = 0.5s Figure 4.5: The canonical 384-DOF deformable body model and inverse simulation of protrusion for a model with two (upper panels) and seven (lower panels) anterior-posterior muscle group divisions. 0 1 A 10 1 P Ac tiv at ion  (% ) Time (s) 0 1 A 0 1 AM 0 1 MA 0 1 M Ac tiv at ion  (% ) 0 1 MP 0 1 PM 10 1 P Time (s) (a) (b) Figure 4.6: Muscle activations predicted to protrude the beam’s tip for a model with two (a) and seven (b) anterior-posterior muscle group divisions. Transverse muscles are activated to compress the posterior section. 82 4.2. Analysis with Canonical Models Figure 4.7: Inverse simulation with the deformable beam model tracking a complex movement trajectory involving combined bending, lengthening, and shortening. motor task. The kinematic target is specified as protrusion of the nodes at the beam’s tip, which is equivalent to elongating the beam while also not compressing the tip laterally or vertically. Figure 4.5 pictures the results of elongation simulation for a model with two anterior-posterior muscle groups (upper panels) and a model with seven anterior-posterior muscle groups (lower panels). The model with more muscle groups deforms more smoothly during elongation. Muscular-hydrostat motor recruitment The muscular-hydrostat theory proposes that movement in non-skeletal mus- cle tissue systems, such tentacles, trunks, and tongues, is generated through the combined effect of tissue incompressibility and the spatial arrangement of muscle fibers [96]. For example, elongation is generated by activation of transversely or helically arranged muscle fibers, while bending is generated by unilateral activation of longitudinal muscle fibers. The theory is sup- 83 4.3. Analysis with Anatomical Models ported by morphological evidence, such as the fact that longitudinal fibers in many tentacles are located along the outer surface of the tentacle in order to provide better leverage for bending [96]. Our simulations of beam elongation support the muscular-hydrostat the- ory as the inverse solution recruits transversely-oriented horizontal and ver- tical muscle groups. Further, we are able to generate combined bending, shortening, and lengthening deformations by specifying a more complex kinematic trajectory for the nodes at the beam’s tip, such as is shown in Fig- ure 4.7. Bending movements were generated by the recruitment longitudinal muscle fibers along the side of the beam toward which bending occurred. 4.3 Analysis with Anatomical Models The canonical models serve to demonstrate specific aspects of the inverse solver implementation within a simplified modeling context. In this sec- tion we discuss the application of inverse modeling tools to more complex anatomical models. 4.3.1 Jaw inverse Inverse simulations with the jaw model predict muscle forces needed to move the jaw in the presence of dynamic constraints at the TMJ. Figure 4.8 pic- tures the results of the inverse simulation of lateral incisor movement. Lat- eral incisor movement is achieved primarily through jaw rotation as straight lateral jaw translation is limited both by ligaments and the shape of the TMJ. We compared two alternative kinematic targets, as shown in Fig- ure 4.8: a full 6-DOF kinematic trajectory (upper panels) and a 3D position trajectory for the lower incisor point (lower panels). The 3D incisor point target provides an under-constrained kinematic target for the jaw, and re- quired an additional lateral TMJ constraint, representing the lateral wall and capsular ligament of the TMJ, in order to induce realistic jaw rota- tion (as opposed to implausible straight lateral translation of the jaw body). Both simulations were found to track the target kinematics with small error. 84 4.3. Analysis with Anatomical Models 6-DOF Jaw Target 3-DOF Incisor Point Target t = 0s t = 0.25s t = 0.5s Figure 4.8: The jaw model during inverse simulation of a right lateral move- ment with a 6-DOF target shown as a cyan wireframe mesh (upper panels) and a 3D incisor point target shown as a cyan point (lower panels). The 3D point target simulation required an additional lateral constraint on the right TMJ to attain rotation from the under-constrained target kinematics. The predicted muscle activation for each condition are plotted in Figure 4.9 and are consistent with the jaw physiology literature that the left-side lateral pterygoid muscle is used in right lateral jaw movement [117]. Interestingly, the 6-DOF-target simulation exhibited larger muscle forces, suggesting in- creased jaw stiffness was required to maintain the specified jaw orientation during the movement, whereas the 3D-point-target simulation required less muscle activity to follow the incisor point trajectory. Application to jaw model validation As stated in the introduction of this chapter, one application of inverse modeling tools is model validation, whereby predicted model muscle acti- vation patterns are compared to recorded EMG for a particular kinematic trajectory. To overcome motor redundancy in the model, the inverse toolset 85 4.3. Analysis with Anatomical Models 6-DOF Jaw 3-DOF Incisor Target Target 0 8 LIP 0 8 LPMH 0 8 RDM Ac tiv at ion  (% ) 0 8 LAMH 0 8 RSM 10 8 RAT Time (s) 0 2 LIP 0 2 LPMH 0 2 RDM Ac tiv at ion  (% ) 0 2 LAMH 0 2 RSM 10 2 RAT Time (s) Figure 4.9: Muscle activations predicted for lateral jaw movement with 6- DOF kinematic target and 3D incisor point kinematic target. Recruited muscle groups included the inferior head of the left lateral pterygoid (LIP), the left mylohyoid (LPMH, LAMH), right masseter (RDM, RSM), and right anterior temporalis (RAT). 86 4.3. Analysis with Anatomical Models provides a number of tunable parameters that choose one particular set of “optimal” muscle activations for a particular movement. However, the “optimality” condition employed by the central nervous system to generate muscle activity is unknown, making comparisons of model predictions to recorded data problematic. Here we propose the application of biomechan- ics modeling to inform the motor recruitment strategies in the context of lateral jaw movement. Lateral jaw movement has been widely analyzed and is known to be primarily driven by contra-lateral (side opposite of the movement) lateral pterygoid muscles [117]. The upper and lower heads of the lateral pterygoid muscle have straight muscle fibers with broad insertion sites on the skull (see Appendix B for detailed description) and therefore have different fiber angles and lengths that change over the course of a lateral jaw movement. The relative activation of different regions of the muscle during lateral movement has been characterized in studies using fine wire EMG [85, 129] with wire locations within the muscle verified through CT imaging of the subject post- recording before the wires were removed [139]. One hypothesis is that the different timing of muscle activation in different regions of the muscle is related to changing biomechanical properties of the muscle fibers throughout the movement, such as the mechanical advantage. Mechanical advantage has been proposed as a means to explain muscle recruitment in both the intercostal muscles [61] and the dorsal interosseous muscles [87]. The inverse simulation tools along with the dynamic jaw model provide a means to perform this analysis. Simulated lateral jaw movement with a biomechanical model that in- cludes detailed muscle fiber information for the lateral pterygoid muscles can make predictions of biomechanical correlates, such as mechanical ad- vantage, during the movement. Mechanical advantage could be formulated with respect to a muscle fiber’s spatial location, e.g. its moment arm for generated the desired torque, or with respect to instantaneous length and shortening velocity which affect its instantaneous force generation capacity (see Section 2.3.3). In a similar way as we have defined a desired stiff- ness goal in the inverse algorithm, one could formulate a optimality term 87 4.3. Analysis with Anatomical Models to maximize mechanical advantage. If the hypothesis of maximizing me- chanical advantage is valid, then inverse simulation of lateral jaw movement should predict similar differential muscle activations in regions of the lateral pterygoid muscle as was observed in the EMG recordings. In addition to the redundancy problem, relating predicted muscle activa- tions to EMG recordings for free jaw movements are challenging because the magnitude of muscle activation in free jaw movements is low, thus providing a low signal-to-noise ratio in recorded EMG signals. We have proposed to use jaw movements under externally applied resistive force in order to elicit larger EMG signal amplitude. Given simultaneously recorded jaw move- ment and external force magnitude and direction, we can use the inverse simulation to predict muscle activity in the model and compare with EMG recordings. We have performed pilot jaw movement and EMG studies with a 1D force sensor to record the subject-applied force resisting the movement direction. We found that the 1D force sensor is insufficient to accurately determine directed applied force and we are currently working to refine the protocol with a 3D force sensor. A further extension would be to apply a controlled force to the jaw from a robotic arm connected to the lower jaw through a custom-fit dental appliance, in a manner similar to that which has been done to measure jaw stiffness experimentally [176], but at higher force magnitudes. 4.3.2 Tongue inverse Similar to the deformable beam model reported in Section 4.2.3, the tongue model includes a larger number of muscle fibers distributed throughout an FEM mesh (see Section 3.1.2 for a detailed description). Controlling the movement and shape deformation of the tongue requires complex patterns of muscle activation and is thus a good candidate for the application of inverse simulation. As an example, we simulated upward and backward movement on the tongue tip, which is an important tongue posture used in the production of an English /r/ sound. We simulated tongue elevation by specifying an arcing upward and backward target trajectory for a set 88 4.3. Analysis with Anatomical Models t = 0s t = 0.25s t = 0.5s t = 0.75s t = 1.0s Figure 4.10: The tongue model during inverse simulation for elevating the tongue tip. Target nodes are shown with semi-transparent cyan points. of forty nodes at the tongue tip. The tongue tip elevation movement is shown in Figure 4.10 and the predicted muscle activations are plotted in Figure 4.11. The superior longitudinal muscle is most active in generating a backward tongue bending, which is consistent with the muscular-hydrostat theory of motor recruitment. In addition, the transverse, mylohyoid, and geniohyoid muscles are activated to prevent lateral and downward expansion of the tongue body in order to maintaining tongue tip protrusion during the tongue tip elevation motion. Application to recorded tongue movement As discussed in Section 2.1.1, Electromagnetic Articulometry (EMA) has been used to record the 3D movement of a small number of points on the tongue surface. This type of data is a natural fit with the inverse trajectory tracking solver that uses the 3D position of FEM nodes as a target kine- matics trajectory. EMA data is limited in spatial resolution as only a small number of markers can be tracked at one time. In speech recording studies, EMA markers are typically located along the mid-sagittal plane in order to best characterize the mid-sagittal tongue contour since speech production is predominantly a bilaterally symmetric motor task. The 3D tongue shape (and consequently the 3D vocal tract shape) is important to the acoustics of vowel sounds and a small number of EMA marker point positions does not adequately describe the tongue’s volumetric extent. In the context of the inverse solver the kinematic information of a small number of surface nodes only provides local movement information for the tongue. Global position information requires the 3D shape of the tongue surface. 89 4.4. Discussion 0 10 MSL 0 10 PSL 0 10 ASL Ac tiv at ion  (% ) 0 10 TRANS 0 10 MH 10 10 GH Time (s) Figure 4.11: Muscle activations predicted for elevating the tongue tip. Re- cruited muscle groups include the middle, posterior, and anterior fibers of the superficial longitudinal muscle (MSL, PSL, ASL), transverse (TRANS), mylohyoid (MH), and geniohyoid (GH) muscles. We are investigating the use of tongue shape information as a means to provide global kinematic targets into the inverse solver. 3D tongue shapes have been segmented from MRI data of static vowel postures [14] for the same subject on whom the tongue model’s morphology is based. Further, the internal deformation of the tongue may be important for determining the relative contribution of antagonist muscle groups. Tagged cine-MRI data can be processed to track the 3D position of tagged points within the tongue tissue [142] and is a promising direction for attaining the necessary kinematic information for inverse tongue model simulations. 4.4 Discussion The inverse-dynamics solver proposed in this chapter is a tool for sys- tematically exploring a biomechanical model’s behavior and for integrating recorded kinematic and other data into biomechanics simulations. Accu- rately measuring all of the kinematic, force, and muscle variables simulta- 90 4.4. Discussion neously for a particular motor task with a single subject is a difficult, if not impossible, task. The inverse-tracking toolset presented in this chapter allows data recorded for a subset of variables to be integrated with a biome- chanical model in order to predict other unmeasured variables. Importantly, inverse solutions predict those variables (muscle forces) that are challenging to record from kinematic variables that are much easier to record. While the inverse tracking method provides a natural fit for integrat- ing recorded data with biomechanical models, it is limited to the analysis of specific motor tasks that are characterized by a full kinematic trajectory. The algorithm requires smooth kinematic trajectories that are differentiated for the velocity trajectory used in Equation 4.11. Kinematic errors arising from numerical differentiation and misalignment of joint centers-of-rotation in limb models have been discussed in the literature [180]. Our method is less sensitive to these errors as the output of the forward-dynamics model is used in a feedback loop, though excessively noisy and/or discontinuous po- sition trajectory are still unacceptable for trajectory-tracking. Interpolation between sparse position measurements will also influence the predicted mus- cle activations. For static position target inputs, a non-linear root-finding optimization, such as the method proposed by Sifakis et al. [179], may be more appropriate and efficient. In addition, low-level trajectory specification is inappropriate for senso- rimotor physiology investigations that are concerned with finding the emer- gent properties of the motor system for high-level motor tasks without pre- supposing the low-level coding (e.g. task-space kinematics versus joint-space kinematics). Such investigations would require higher level sensorimotor components, such as reflexes, the spinal-cord, brainstem, and motor cortex, built “on top of” biomechanics models in order to specify high-level motor tasks as input to inverse simulation. For example, one could specify a start- ing and ending point for a desired movement in a model without specifying the precise trajectory in between. We also currently neglect the calcium-dependent activation dynamics of muscle tissue, which is typically modeled as a first-order low-pass filter [216]. The consequence of this assumption is that predicted muscle forces could 91 4.5. Directions change faster than physiologically possible, however this was not found to occur for target movements with physiologically plausible velocities. Also predicted muscle activations would likely lag recorded EMG signals due to the activation dynamics of real muscle tissue. Thelen and Anderson [196] consider muscle activation dynamics within their Computed Muscle Control framework and incorporating activation dynamics into our model one of our planned directions. Our current trajectory tracking algorithm does not include the unilateral constraints in the forward dynamics that arise from joint limits and contact. The presence of inequality constraints in Equation 4.7 would result in a non quadratic optimization problem for Equation 4.15, which would require more complex non linear programming methods. We have found that contact constraints for deformable bodies can be treated as bilateral for the duration of an integration timestep (see Section C.1.4). This could lead to oscillating or sticking behavior; however, this has not been found to be an issue in our simulations as the constraints arising from deformable contact tend to be reasonably decoupled. Handling contacts in this fashion allows them to be included in our inverse solutions, but the scheme is less appropriate for rigid-body contacts that result in more tightly coupled constraint situations. 4.5 Directions Our short term directions for modeling applications of the inverse toolset focus on integrating recorded human data with models for the purposes of validation. The inverse simulation of lateral jaw movement presented in Section 4.3.1 lays the groundwork for a more comprehensive study comparing the results to EMG recordings. This will require a more sophisticated model of the lateral pterygoid muscle that includes muscle vectors for the sub- region fibers in both the upper and lower heads of the muscle. Mechanical advantage as a goal for motor recruitment would be examined in the context of lateral pterygoid recruitment. We are also interested in investigating motor control strategies in oro- facial motor tasks, including speech production. A long-standing debate in 92 4.6. Summary the speech motor control literature concerns the validity of the EPH, which posits that speech gestures are controlled as static target-to-target posi- tions [66]. EPH has been used in the reference tongue model [29]; however, the technical capability to test alternative control strategies that incorpo- rate dynamics information has been lacking. The inverse simulation toolset described in this chapter fills this technological gap and can be applied to evaluate the significance of controlling dynamics in speech movements. We plan to compare predicted muscle activations from inverse-dynamics versus EPH control for two data-sets from the same speaker: mid-sagittal EMA recordings of speech utterances and 3D tongue shapes segmented from mag- netic resonance images MRI. We expect that dynamics may be important in vowel-vowel tongue transitions that have large movement amplitudes. In future studies, we plan to use the inverse modeling tools with the full coupled jaw-tongue-hyoid model as described in the previous chapter. To date, the inverse methods have been applied to the jaw and tongue models in isolation in order to assess the results within a simpler context. In the previous chapter, we illustrated that coupling effects can be signifi- cant in combined jaw-tongue-hyoid movements, such as speech production, and therefore we expect that accurate prediction of jaw and tongue mus- cle activity using the inverse toolset will require the full model of coupled jaw-tongue-hyoid biomechanics. We expect inverse simulations with the jaw- tongue-hyoid will provide insight into the biomechanical underpinnings of co-articulation effects in speech production. 4.6 Summary Computational models of biomechanical systems have greatly increased in fidelity and complexity. A limiting factor in the application of such mod- els to particular applications in biomedicine and biological science remains the ability to generate movement simulations without trial-and-error speci- fication of input muscle activity. In this chapter, we developed an inverse- dynamics algorithm for automatically predicting muscle activity to drive biomechanics models with combined skeletal and muscle-tissue structures. 93 4.6. Summary We have extended previously reported trajectory tracking-based inverse dynamics approaches to work within the ArtiSynth framework for muscle- activated rigid and deformable body models with dynamic constraints. We have also included additional target parameters that can be used to better control the selection of muscle activation given motor redundancy. Con- straint force magnitude targets can be controlled such as to generate bite forces at the teeth for the jaw model. Co-activation of antagonist muscles can also be modified to control stiffness. In sum, these contributions help to cre- ate a more comprehensive inverse-dynamics toolset than has been previously proposed. Using the inverse toolset, we simulated motion and deformation of muscular-hydrostat models, including a beam model and a physiologically- based tongue model, illustrating muscle activations that are consistent with theoretical proposals. Also, we simulated lateral jaw movement and found muscle activations consistent with published jaw physiology. In the next chapter, we provide a case study analyzing segmental jaw resection and reconstruction using the forward dynamic jaw model described in Chapter 3 and the inverse modeling toolset described in this chapter. 94 Chapter 5 Segmental Jaw Surgery Models In the previous two chapters we have described a new biomechanical model of the jaw-tongue-hyoid complex and new simulation tools for making auto- matic predictions of muscle force patterns required to elicit specific move- ment and force trajectories from the model. Our primary motivation for building these modeling tools is in their application to orofacial biomedicine. In this chapter, we investigate one example biomedical application in detail, segmental jaw surgery, in order to demonstrate the utility of the modeling toolset. Treatment of oral cancer commonly involves surgical resection of cancer- ous tissue, in addition to radiation therapy and chemotherapy. Depending on the size and location of the lesion, tissue resection can involve the por- tions of the mandible (hemimandibulectomy), tongue (glossectomy), floor of mouth, and associated muscles. Vascularized osteocutaneous, osteomy- ocutaneous and alloplastic grafts are commonly used to restore mandibular continuity after hemimandibulectomy [74, 77, 116, 153, 165]. Grafts provide a functional joint on the affected side [114, 152], though articular com- plications can include erosion of the temporal fossa, dental malocclusion, infection, and graft migration [33, 114, 143]. With or without jaw recon- struction, hemimandibulectomy significantly alters jaw biomechanics and deficiencies in mastication, speech and other orofacial functions are often observed [2, 40, 73, 77]. Typically, the mandible deviates to the defective side on opening, and chewing is performed on the normal side [161, 166]. Altered sensation, salivary flow and biomechanics likely affect biting, ma- 95 5.1. Model Creation nipulating and comminuting food [40, 113, 166, 200] and are concerns in oral rehabilitation. Functional jaw recordings, including movement, bite force, and EMG, have been used to analyze the intact masticatory system (see Section 2.1 and [117]), but are rarely performed after hemimandibulectomy. Anecdotal information and clinical observations suggest that jaw-opening is not unduly restricted, and that biting is frequently accompanied by frontal jaw rota- tion. Contributing factors include an asymmetrical musculature, a unilat- eral articulation in mandibles without continuity, and post-operative tissue- scarring. In this chapter, we analyze hemimandibulectomy models with and with- out mandibular continuity. Section 5.1 describes the models’ creation. Sec- tion 5.2 describes forward dynamic simulations performed with the goal of replicating clinically observed post-operative deficits, including atypical resting postures and abnormal movements during opening and clenching. Section 5.3 reports the results of inverse simulations used to determine what is physically possible for the ungrafted hemimandibulectomy model and whether particular patterns of muscle use could be employed to com- pensate for given functional deficiencies. 5.1 Model Creation The hemimandibulectomy models, depicted in Figure 5.1, are based on the intact jaw model described in Section 3.1.1. Model NOCON (no condyle) simulated a left-side, composite jaw resection from the condyle to the left canine, without restoration of mandibular continuity. Model CON (graft- related condyle) simulated a left-side resection with continuity restored by means of an alloplastic graft similar to that described by Marx et al. [114]. The inertial properties of the mandible fragment and graft were computed based on its geometric shape resulting in a new mass of 126 g from 200 g for the normal jaw and a new centre-of-mass shifted from midline in the normal jaw to inside the body of the jaw fragment inferior to the second molar. The hyoid was fixed in both the models, which would be achieved 96 5.1. Model Creation Figure 5.1: The two models used in the study. Model NOCON shows jaw resection without mandibular continuity. Model CON shows continuity re- stored by alloplastic grafting. Both models have left side muscle loss. by activation of infrahyoid muscles to stiffen the hyoid; hyoid movement in hemimandibulectomy patients has not been reported, but likely differs from normal subjects due to scar tissue on the affected side. In both models, left-side muscles were removed as part of the resection. The intact muscles are shown in Figure 5.2 and included the right ante- rior, middle, and posterior temporalis (RAT, RMT, RPT), right deep and superficial masseter (RDM RSM), right medial pterygoid (RMP), right su- perior and inferior lateral pterygoid (RSP, RIP), right mylohyoid (RMH) right and left geniohyoid (RGH, LGH), and right and left digastric (RDI, LDI) muscles. Hill-type muscle models simulated individual muscle CSA and length-tension properties and produced passive forces proportional to muscle stretch and active forces proportional to muscle activation. Muscle CSA are listed in Table 3.1. When the jaw was in the maximal intercuspal position, the right condyles in both models, and the alloplastic condyle in CON, were centered in their articular fossae. The joints were modeled as a bilateral constraint surface, curvilinear in the anteroposterior direction to represent the articular fossa 97 5.1. Model Creation Figure 5.2: Model NOCON showing intact muscle groups including the right anterior, middle and posterior temporalis (RAT, RMT, RPT), right deep and superficial masseter (RDM RSM), right medial pterygoid (RMP), right superior and inferior lateral pterygoid (RSP, RIP), right mylohyoid (RMH), right and left geniohyoid (RGH, LGH), and right and left digastric (RDI, LDI) muscles. The effect of post-operative scarring is represented by the SPRANT spring element. 98 5.2. Forward Dynamics: Deficit Simulations and eminence, with no mediolateral constraint. A posterior constraint was included to represent the posterior aspect of the articular fossa. In all sim- ulations the joint worked in compression, i.e. the forces in the system never worked to pull the joint apart, therefore it was modeled with bilateral con- straints. For clenching simulations we added an additional planar constraint at the right first molar to simulate tooth contact. Passive soft-tissue forces representing tissue scarring were modeled with linear, damped springs. These had stiffnesses of 200 N/m and damping- constants of 10 Ns/m to restrict jaw motion without eliminating it. The springs permitted incisal separations of at least 15 mm, and counteracted the tendency of the jaw to position itself to the right as a result of passive muscle forces on that side. In NOCON, a spring attached to the anterior portion of the jaw fragment (SPRANT) drew this end of the native mandible laterally, posteriorly and inferiorly from its initial position. In CON, a posterior spring attached to the gonial region of the graft (SPRPOST) initially drew the gonial angle of the graft inferiorly and posteriorly, and thereafter created tensile forces proportional to gonial displacement in any direction. The scar springs and spatial coordinate frame are shown in Figure 5.3. The x-y plane was oriented to the Frankfort horizontal plane. Each origin was in the body of the hyoid, with positive x-,y- and z-axes indicating left lateral, posterior and superior respectively. 5.2 Forward Dynamics: Deficit Simulations 1 Forward dynamic simulations were created with ArtiSynth (see Appendix C). The goal of the forward dynamic simulations was to evaluate if the modi- fied jaw model could replicate commonly observed deficits in jaw movement during rest, muscle-driven opening, applied force opening, and unilateral clenching. Additionally, we aimed to compare differences between the NO- CON and CON models. 1A version of the section has been published in Hannam et al. [72]. Dr. Hannam wrote the text for the manuscript, which I have adapted and included in this section. 99 5.2. Forward Dynamics: Deficit Simulations Figure 5.3: Model conventions and restraints. Model illustrated is NOCON with soft-tissue spring SPRANT. In CON, spring SPRPOST was attached to graft’s gonial angle (not shown). For clarity, only left side articular con- straining surfaces are illustrated. Black spheres indicate incisor point, left condylar-center, and molar contact point locations in intercuspal position. 100 5.2. Forward Dynamics: Deficit Simulations 5.2.1 Simulation descriptions Simulated external forces on the jaw, and muscle activation, either singly or in groups, were chosen as a means to provide insight into the movement patterns and jaw instabilities seen clinically. The relaxed rest position of the mandible without postural muscle activity (RRP) was assessed with soft-tissue restraint 1 second after gravitational acceleration from the in- tercuspal position. External force applied downwards to the mandibular incisors perpendicular to the occlusal plane (FORCE) simulated manual jaw-opening from RRP with soft-tissue restraint. This force increased at 10 N/second, and the simulation was terminated when the incisor-point reached approximately 30 mm of inferior displacement. Jaw opening due to muscle activation (OPEN) simulated voluntary jaw opening from RRP with soft-tissue restraint. Here, RSP, RIP, RDI and LDI were driven simul- taneously, reaching 10% of maximum contraction in 0.5 seconds. Unilateral molar contact on the unaffected side on jaw closure from RRP (UNIMOL) was simulated by activating individual jaw-closing muscles in the presence of soft tissue restraint. RAT, RMT, RPT, RDM, RSM and RMP were each driven independently to 10% of maximum contraction in 0.5 seconds, similar to the protocol used by Koolstra and Van Eijden [97]. 5.2.2 Simulation results Relaxed resting jaw posture The displacement of the incisor-point at RRP is plotted in Figure 5.4 for both models. In NOCON, SPRANT caused the jaw to rotate to the defect side, the incisor point displacing markedly left, posteriorly and inferiorly, and there was minimal displacement of the right condyle. In CON, the effect with SPRPOST was less, the jaw remaining closer to the midline at a reduced incisal separation. 101 5.2. Forward Dynamics: Deficit Simulations Force-induced jaw-opening Trajectories of incisor-point displacement during FORCE are compared in Figure 5.4. In both models, inferior displacement of the incisor point closely approximated the opening target of 30 mm with forces less than 5 N (3.29 N and 3.38 N for NOCON and CON). In NOCON, the incisor point ap- proached the midline as the jaw opened, ending 8.38 mm posterior to its maximal intercuspal position (not shown). The path in CON paralleled that in NOCON, but began and ended on the right (unaffected) side, indicating frontal clockwise jaw rotation. Muscle-induced jaw-opening Incisor point displacements during OPEN are compared in Figure 5.4. In both models, the incisor point moved to the left (affected) side. In NOCON, the marked lateral movement of the incisor point to the affected side partly resulted from its initial position in RRP, where it was already displaced to the left. Greater lateral deviation occurred in CON, where the incisor point moved initially to the right, then markedly to the left. Muscle-induced jaw-closing The effects of muscle activation in the UNIMOL simulation were pronounced in the NOCON model and resulted in significant postural instability (large frontal-plane rotations of the jaw fragment), whereas the effects in the CON model were less striking and generally resulted in more stable behavior. Graphic examples of jaw rotation in NOCON are illustrated in Figure 5.5. In NOCON, the actions of RAT, RMT RPT and RDM were similar. The most common effect was movement of the incisor-point laterally to the left and posteriorly, especially for RPT and RDM. This was associated with marked lateral displacement of the right condyle. RPT caused significant superior movement of the incisor point, and excessive movement of the right condyle posteriorly and inferiorly along its posterior planar constraint. These marked translations and rotations are shown collectively in Figure 5.5. RSM activation caused excessive lateral movement of the incisor point to 102 5.2. Forward Dynamics: Deficit Simulations -10 -5 0 5 -30 -25 -20 -15 -10 -5 0 Lateral  Displacement  (mm) Ve rti ca l  D isp lac em en t (m m ) FORCE -10 -5 0 5 10 -30 -25 -20 -15 -10 -5 0 Lateral  Displacement  (mm) Ve rti ca l  D isp lac em en t (m m ) OPEN  10 RRP OPEN NOCON IP CON LEFTRIGHTLEFTRIGHT Figure 5.4: Frontal view of mandibular incisor-point displacements dur- ing FORCE and OPEN. Movements began from jaw’s relaxed rest position (RRP) and are referenced to maximal intercuspal position IP (large crosses). 103 5.2. Forward Dynamics: Deficit Simulations the left (deficient) side, as well as posteriorly and inferiorly. Here, there was minimal displacement of the right condyle, indicating predominant 3- dimensional jaw rotation. RMP activation moved the incisor point exces- sively left, posteriorly and superiorly. Although there was no displacement of the condyle, RMP had a strong rotational action on the jaw, different to that due to RSM activation (Figure 5.4). 5.2.3 Discussion Dynamic modeling can be used to study jaw biomechanics by simulating the effects of mandibular surgery and reconstruction. The approach is physics- based, and suitable for solving complex dynamic interactions amongst mul- tiple components. Anatomical structures can be readily modified, and tis- sue properties can be assigned provided their parameters are known. Pre- dictably, the RRPs occurred at a wider incisal separation than clinical postu- ral rest because low-grade postural muscle activity was not simulated [213]. Soft tissue forces on the defect (operated) side, especially those due to scar- ring, would normally have a restraining effect on this motion, so the RRPs obtained with SPRANT and SPRPOST seem in line with clinical impres- sions. Also, less incisal separation would be anticipated clinically due to postural muscle activity. Jaw deviation in RRP was less for CON than for NOCON, suggesting that scarring in a jaw with a bilateral articulation could result in a relatively normal RRP. The motions caused by FORCE reflected the different components in the two models. Both easily reached their opening targets of 30 mm with a low applied forces of 3-4 N. In NOCON, SPRANT functioned as a simple tether, freely allowing incisal separation. In CON however, SPRPOST acted closer to the jaw’s transverse axis of rotation and limited movement in any direction. High stiffness values assigned here would be expected to restrict jaw motion, and it is significant that the stiffness of wounded porcine skin is higher than the 200 N/m used in the present study [37]. A clinical protocol similar to the FORCE simulations (applying known forces to the mandibular incisor of a patient and tracking jaw motion) might be useful for estimating 104 5.2. Forward Dynamics: Deficit Simulations Figure 5.5: NOCON jaw postures caused by individual closing-muscle acti- vation in right molar contact. Simulations halted when jaw motion became unrealistic at 1.36, 1.45, 1.43 and 1.19 s for RMT, RDM, RMP and RSM respectively. Muscle forces (Fm) and torques (Tm) expressed at jaw centers of mass indicate directions only (scaled for clarity). + denotes incisor-point location at intercuspal position. Grid spacing is 10 mm. 105 5.2. Forward Dynamics: Deficit Simulations scar stiffness in clinical situations. The difference between FORCE and OPEN trajectories can be explained by the primarily inferiorly-directed force in the former, and the primarily oblique muscle forces in the latter. In OPEN, the bilateral articulation in CON reduced this lateral deviation, but did not eliminate it. Deviated jaw- opening is a common clinical observation associated with muscle loss, and has functional implications with respect to mastication, since jaw-closing must begin from the defect side, and the maximal intercuspal position ap- proached mediolaterally. In the present study, wider incisal separations could have been reached with more muscle drive, and using additional mus- cles might have increased it further. More drive in the digastric and ge- niohyoid muscles, however, would have resulted in less lateral jaw deviation because these muscles have poor angles of attack, and their effectiveness diminishes as the jaw opens [99]. Analysis of the biomechanical role of single muscles is a unique feature of computational modeling since living subjects are unable to activate jaw muscles individually. Clinically unrealistic movements resulting from single- muscle activation in both models were therefore not surprising, but were helpful in revealing the actions of muscles likely contributing to mandibu- lar instability. The marked rotation caused by RSM after molar contact in both models, and by RMT in NOCON, substantially explained clinical observations of frontal plane rotation. The tendency of RAT, RMT, RPT and RDM to translate the jaw laterally in a mandible without continuity would normally be resisted by the temporomandibular ligament [97]. Con- tinuous, or perhaps exclusive, use of such muscles may explain mandibular lateral displacement during occlusal contact sometimes observed clinically in mandibular resection patients without reconstruction. The jaw instabilities demonstrated in UNIMOL partly explained the challenges for patients having to find new strategies of muscle contraction. For further investigation in the following section, we intend to determine the extent to which combined muscle use might provide stability by using inverse dynamics simulation. 106 5.3. Inverse Dynamics: Compensatory Simulations 5.3 Inverse Dynamics: Compensatory Simulations Having simulated post-operative deficits in a hemimandibulectomy model, we proposed to use the model to predict possible compensatory muscle pat- terns that hypothetically could be employed by patients and reinforced by post-operative rehabilitation. Without experimental data upon which to base muscle input, however, forward-dynamic simulations require trial-and- error approaches to determine whether a particular jaw posture or movement can be attained by activating the remaining jaw muscles. Here we describe the use of the inverse-dynamics techniques described in Chapter 4 to reveal the muscle forces needed to reduce deficits in the NOCON model. 5.3.1 Simulation descriptions The compensatory simulations were chosen to counteract the two main deficits observed in the forward dynamic simulations: firstly, lateral devia- tion during jaw opening and, secondly, instability during unilateral clench- ing. The ungrafted case, model NOCON, exhibited more pronounced devia- tions and instabilities and was therefore used in compensatory simulations. Hinge movement simulations We were interested in determining the muscle forces needed for midline jaw movements with the model. Therefore we simulated movement between three jaw postures that are illustrated in Figure 5.6: • REST a relaxed rest posture deviated toward the affected side and rotated clockwise in the frontal-plane • OPEN a midline opening posture with 20 mm inter-incisal separation and no frontal-plane rotation, a typically maximal opening gape during chewing • CLOSE a midline posture with the jaw closed to just before first tooth contact 107 5.3. Inverse Dynamics: Compensatory Simulations Two movement simulations were performed: from REST to OPEN mov- ing the jaw fragment from the deviated rest posture to the midline, and from OPEN to CLOSE moving the jaw fragment in centric relation relative to the maxilla, for a “hinge-like” movement with no frontal-plane rotation. Smooth target position trajectories were generated with quintic splines from the start posture to the end posture over a 0.5 s duration. The target veloc- ity trajectories were computed online with finite-differencing which provided on-line correction of position errors as the simulation progressed. We spec- ified the full six degree-of-freedom position trajectory of the jaw in order to control its orientation. The target trajectories referenced to the 3D in- cisor point movement are shown for REST-to-OPEN in Figure 5.7a and for OPEN-to-CLOSE in Figure 5.8a. Unilateral clench simulations We were interested in determining if a stable unilateral clench could be achieved through the recruitment of an ensemble of closing muscle groups with appropriately-balanced activation. We performed clenching simulations with different jaw postures to investigate the effect of jaw position on the ability of the model to generate bite force. 5.3.2 Simulation results Deviated rest to midline retrusive open The muscle patterns predicted by the inverse simulation for the REST-to- OPEN movement are shown in Figure 5.7b and Table 5.1. Incisor point movement to midline along with counter-clockwise rotation of the jaw frag- ment to a neutral orientation was accomplished by co-activation of the right- sided digastric and posterior temporalis muscles. The average position error of the incisor point during the movement was 0.36 mm, 0.22 mm, and 0.39 mm in the left, anterior, and inferior directions respectively. Posterior tem- poralis inserts into the coronoid process and has a force vector best angled to apply the required torque to rotate the fragment to a symmetric midline 108 5.3. Inverse Dynamics: Compensatory Simulations Figure 5.6: Medial and frontal views of the jaw postures for relaxed rest (REST), 20 mm retrusive midline open (OPEN), and retrusive midline close to before tooth contact (CLOSE). + denotes incisor-point location at inter- cuspal position. Grid spacing is 10 mm. 109 5.3. Inverse Dynamics: Compensatory Simulations 0 0.1 0.2 0.3 0.4 0.5 −20 −15 −10 −5 0 5 10 Po sit ion  (m m ) 0 0.1 0.2 0.3 0.4 0.5 −25 −20 −15 −10 −5 0 5 Ve loc ity  (m m /s) Time (s)   0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 M us cle  A cti va tio ns  (% ) 0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 M us cle  F or ce s ( N) Time (s)   RDI RMH RIP RSP RPT RSM TARGET TRAJECTORIES MUSCLE PREDICTIONS Figure 5.7: Movement simulation from REST to OPEN. Target position and velocity trajectories for mandibular incisor point in lateral (solid lines), vertical (dotted lines), and anteroposterior directions (dashed lines) shown in left-side panels. Muscle activations and forces of the primary muscles used to track target trajectory shown in right-side panels. posture. Digastric is activated to open the jaw to 20 mm as is the case in retrusive opening with an intact mandible. Hinge closing from midline open posture The muscle patterns predicted to move the jaw fragment along a midline hinge closing trajectory in the OPEN-to-CLOSE movement are shown in Figure 5.8b and Table 5.1 The average position error of the incisor point during the movement was 0.08 mm, 0.02 mm, and 0.03 mm in the left, pos- terior, and inferior directions respectively. Posterior temporalis was again recruited in order to keep the fragment at midline and its activity increased 110 5.3. Inverse Dynamics: Compensatory Simulations 0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 50 60 70 M us cle  A cti va tio ns  (% ) 0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 50 M us cle  F or ce s ( N) Time (s)   RDI RIP RSP RPT RDM 0 0.1 0.2 0.3 0.4 0.5 −20 −15 −10 −5 0 5 10 Po sit ion  (m m ) 0 0.1 0.2 0.3 0.4 0.5 −40 −20 0 20 40 60 Ve loc ity  (m m /s) Time (s)   TARGET TRAJECTORIES MUSCLE PREDICTIONS Figure 5.8: Movement simulation from OPEN to CLOSE. Target position and velocity trajectories for mandibular mid-incisor point in lateral (solid lines), vertical (dotted lines), and anteroposterior directions (dashed lines) shown in left-side panels. Muscle activations and forces of the primary muscles used to track target trajectory shown in right-side panels. during the closing movement. Lateral pterygoids were also co-contracted, presumably to balance the lateral temporalis force and to maintain a midline movement. Unilateral clench The three simulated clenching postures are illustrated in Figure 5.9 and predicted muscle activation and force magnitudes are provided in Table 5.2. For comparison, maximal first molar bite force for an intact mandible is within the range of 216-740N [210]. Clenching with the jaw fragment de- viated toward the affected side was the most stable act. Moving toward a 111 5.3. Inverse Dynamics: Compensatory Simulations rest posture peak open posture peak close posture rest→open open→close %, [ N ] %, [ N ] %, [ N ] %, [ N ] %, [ N ] RPT 2.3 [ 1.8 ] 11.2 [ 7.9 ] 0.2 [ 0.3 ] 64.9 [ 49.0 ] 40.2 [ 30.4 ] RSM 0.0 [ 0.6 ] 0.0 [ 1.3 ] 0.0 [ 1.3 ] 0.0 [ 1.3 ] 0.0 [ 0.0 ] RDM 0.0 [ 0.2 ] 0.0 [ 0.3 ] 0.0 [ 0.3 ] 7.4 [ 5.8 ] 2.7 [ 2.2 ] RMP 0.1 [ 0.9 ] 1.1 [ 1.9 ] 0.1 [ 1.0 ] 1.1 [ 2.0 ] 0.0 [ 0.0 ] RSP 0.7 [ 0.2 ] 8.6 [ 2.5 ] 4.6 [ 1.3 ] 14.8 [ 4.3 ] 8.7 [ 2.5 ] RIP 1.2 [ 0.8 ] 6.1 [ 4.1 ] 0.0 [ 0.0 ] 38.1 [ 25.5 ] 23.9 [ 16.0 ] RDI 1.5 [ 0.3 ] 20.2 [ 2.5 ] 11.7 [ 1.0 ] 23.1 [ 5.5 ] 6.7 [ 2.7 ] RMH 0.4 [ 0.1 ] 4.2 [ 0.7 ] 3.4 [ 0.5 ] 3.4 [ 0.5 ] 0.0 [ 0.0 ] Table 5.1: Muscle activations (%) and forces (N) for movements between REST, OPEN, and CLOSE postures. Muscle force includes both passive force due to muscle stretch and active force proportional to activation. midline position required muscle recruitment to maintain the posture and reduced bite force magnitudes. At 15 mm lateral and backward deviation of the mandibular mid-incisor point, the simulation was able to generate 124 N of force at the bite constraint; however positioning the jaw medially to a 10 mm deviation reduced the bite force to 111 N and greatly increased the co-activation of opener muscles. In the midline posture, the inverse sim- ulation with a high target bite force was unable to prevent rotation of the jaw fragment. However, a stable midline posture with no rotation of the jaw fragment was possible with a diminished bite force of 25 N. The main closer muscles recruited for clenching were anterior temporalis, medial pterygoid, and deep masseter. 5.3.3 Discussion The simulated recruitment of antagonist muscles for both free jaw move- ments and unilateral clenching suggests that functional deficit caused by unilateral muscle and articular loss may be at least partly overcome by stiffening the system with opposing muscles. The current simulations thus serve as a proof-of-concept of an inverse modeling approach to determine the biomechanical plausibility of hypothesized movements and associated mus- 112 5.3. Inverse Dynamics: Compensatory Simulations Figure 5.9: Superior and frontal views of jaw postures used in clenching tasks. The model was able to generate a 124 N bite force at 15 mm incisor deviation (15DEV); a 111 N bite force at 10 mm incisor deviation (10DEV); and a 25 N bite force at midline intercuspal position (IP). + denotes incisor- point location at intercuspal position. Grid spacing is 10 mm. 113 5.3. Inverse Dynamics: Compensatory Simulations 124 N clench at 111 N clench at 25 N clench at 15DEV 10DEV IP %, [ N ] %, [ N ] %, [ N ] RAT 89.9 [ 140.9 ] 90.1 [ 141.7 ] 16.3 [ 25.8 ] RMT 0.0 [ 0.1 ] 0.0 [ 0.0 ] 4.7 [ 4.5 ] RPT 0.0 [ 0.1 ] 0.0 [ 0.1 ] 2.8 [ 2.1 ] RSM 16.0 [ 29.2 ] 19.3 [ 36.5 ] 6.8 [ 13.0 ] RDM 41.9 [ 34.1 ] 86.8 [ 70.2 ] 1.2 [ 1.0 ] RMP 88.8 [ 154.5 ] 89.0 [ 155.5 ] 11.5 [ 20.1 ] RSP 4.5 [ 0.8 ] 15.6 [ 2.7 ] 0.3 [ 0.0 ] RIP 5.0 [ 2.5 ] 78.0 [ 39.0 ] 5.1 [ 2.5 ] LDI 27.2 [ 4.8 ] 14.9 [ 5.1 ] 0.0 [ 0.0 ] RDI 44.2 [ 20.0 ] 77.5 [ 37.2 ] 0.0 [ 0.0 ] RMH 4.7 [ 0.5 ] 2.3 [ 0.3 ] 0.0 [ 0.0 ] LGH 5.3 [ 0.6 ] 8.3 [ 1.3 ] 0.0 [ 0.0 ] RGH 5.7 [ 0.7 ] 8.5 [ 1.4 ] 0.0 [ 0.0 ] Table 5.2: Muscle activations (%) and forces (N) for unilateral clenching simulations. Muscle force includes both passive force due to muscle stretch and active force proportional to activation. 114 5.3. Inverse Dynamics: Compensatory Simulations cle patterns in an altered jaw system. Comparison of the model predictions with patient data will require significant quantitative clinical measurements on patients. In the reported movement simulations, posterior temporalis was used to provide a torque to move the incisor point to the midline in a right lateral movement to compensate for the left lateral deviation caused by passive scar tissue forces. Such movement would normally be accomplished by the con- tralateral lateral pterygoids which are missing in the hemimandibulectomy case. Ipsilateral temporalis contributes to normal lateral movement [117]. Therefore, it is plausible that right-side posterior temporalis is recruited to compensate for the missing lateral pterygoid muscles in a left-sided deficit. Ipsilateral lateral pterygoids co-activate with posterior temporalis, which is atypical, but necessary in the model to generate medial forces at the condyle to balance the lateral forces generated by posterior temporalis. We expect that the inclusion of lateral constraint at the joint, such as the lateral aspect of the articular fossa or the temporomandibular ligament, would reduce the need for lateral pterygoid co-activation. Predicted muscle activations are affected by the choice of regularization term in the optimization (see Section 4.1 for discussion). Different optimality conditions have been proposed for a variety of physiological and numerical reasons (see Erdemir et al. [53] and Ait-Haddou et al. [4]). We use an `2- norm weighted by the inverse of each muscles CSA. The main difference observed with CSA-weighted regularization was in the relative contribution of synergistic muscles. In the lateral pterygoids, for example, the inferior- head activity was increased and superior-head activity decreased when using the CSA-weighted regularizer as compared to an unweighted `2-norm, be- cause inferior-head is a significantly larger muscle. The relative scaling of antagonist groups was unaffected by regularizer weighting as it is prescribed by the requirements of tracking the target trajectory. Our current tracking algorithm is formulated only for bilateral con- straints. Contact modeling requires unilateral constraints that add inequal- ity conditions on the dynamics equation (see Section C.1). This limitation reduces the number of tasks we can currently model: in hinge jaw opening 115 5.4. Directions the condyle does not move forward off the posterior joint constraint and in clenching with the teeth always in contact. As discussed in Section 4.4, it may be possible to treat contacts as bilateral constraints for the duration of a timestep and break contact if the constraint is found to be pulling apart. The bite constraint, modeled here as a planar contact surface aligned to the occlusion plane, simulates a patient with flat teeth. Here, we found it impossible to hold a midline intercuspal position and generate a bite- force more than 25 N. Moreover, all bite forces were considerably less than those normally generated in the intact jaw during tooth-clenching, i.e. 216- 740N [210]. A reduced capacity to generate bite-force in the present case can be expected given the loss of half the closing muscles. Notwithstanding the additional, atypical contributions of antagonistic muscles made available in our simulations, the postural effect on bite-force generation reflects the extreme biomechanical conditions needed to create any significant unilateral bite force in the hemimandibulectomy patient. 5.4 Directions Our current models were limited by the absence of a tongue, a fixed hyoid apparatus, and arbitrary soft-tissue scar forces. Incorporating a 3D FEM model of tongue tissue is a straight-forward direction, given the coupled jaw- tongue-hyoid model presented in Chapter 3, and would also allow for the investigation of glossectomy cases. Incorporating a physiologically-based 3D FEM model of scar tissue is significantly more challenging. Detailed information on the extent and location of scarring in hemimandibulectomy patients is currently unavailable. The physical properties of soft-tissues altered by wound-healing and radiation therapy also remain undefined. An alternative approach to improve the realism of scarring in the model would involve characterizing the functional effect of scarring on patients. One approach to measure scar-induced jaw stiffness would be to systematically apply external forces to a patient’s lower jaw, while simultaneously tracking jaw motion. Such techniques have been used to characterize jaw stiffness and muscle viscoelasticity in normal subjects [147, 176]. 116 5.5. Summary Our current models use a planar bite constraint, representing contact between flat teeth. It is possible that a modified dental interface, such as inclined tooth surfaces to guide the jaw fragment toward the midline during clenching, has a significant effect on inverse predictions of muscle activity during clenching. In future work, actual occlusal surface shapes used in dental reconstruction could be incorporated in our model, making it possible to study the effects of occlusal configurations on system biomechanics and predicted muscle patterning. 5.5 Summary Computer models of anatomical systems can be used to analyze structurally altered and functionally compromised cases. In this chapter, we have de- scribed a study analyzing jaw movements in compromised mandibles with and without continuity, as would occur in patient after segmental jaw surgery consequent to oral cancer. Through forward dynamic simulations, we illus- trated the functional consequences of missing muscle and bone structure and recreated clinically observed deficits in the model. We also compared the mechanics of deficits in cases with and without jaw grafting. Further, with new inverse modeling techniques, were we able to predict what muscle forces would theoretically be needed to compensate for post-operative deficits. These compensatory muscle predictions could be used to inform patient- specific rehabilitative therapy. Through the case study of jaw surgery, have demonstrated that computer modeling presents a promising approach to un- derstanding the biomechanics of surgically altered musculoskeletal systems. In the final chapter, we summarize the contributions of the dissertation and propose a few broader directions for future investigation. 117 Chapter 6 Conclusions To conclude, this chapter reviews and summarizes the contributions of this dissertation. The research publications and presentations associated with this work are also listed. Finally, longer-term directions drawing on this research are presented along with concluding comments. 6.1 Dissertation Contributions The contributions of this dissertation include an advanced jaw-tongue-hyoid model, inverse simulation toolset, and model-based analysis of jaw surgery. While these contributions are made within the context of orofacial anatomy and function, the framework and methodology is widely applicable to mus- cle activated skeletal and soft-tissue systems and biomedical applications involving biomechanical analysis. Modeling of coupled jaw-tongue-hyoid biomechanics i. Created a novel model of the jaw-tongue-hyoid system. By adapting previously reported reference models to a single-subject within the same simulation platform we have created the most complete phys- iologically based computer model of the oral region to date. The model provides a foundation on which to analyze the biomechanics and neu- romotor control of oral motor tasks, such as mastication, deglutition, and speech production, which involve the coordination of multiple oral articulators. ii. Demonstrated the significance of coupling. Through forward dy- namic simulation with isolated muscle activations we were able to deter- 118 6.1. Dissertation Contributions mine that the mechanical coupling of tongue-muscles acting on the jaw, and vice versa, were significant. In particular, the elastic connection of tongue tissue between the jaw and hyoid were observed to restrict jaw movement. Also, the movement of the tongue within the oral cavity in- duced non-negligible jaw movements with the jaw at rest. Jaw-tongue coupling effects are therefore an important factor for consideration in future orofacial modeling studies. iii. Compared results with recorded tongue velocity and pres- sure. Through forward dynamic simulation of tongue muscle activation we evaluated biomechanical variables in comparison to reported hu- man tongue measurements. The velocity profile of simulated backward tongue movement compared well with recorded EMA data. Tongue to palate pressure simulated with maximum muscle activation in tongue- palate contact agreed with average valued human measurements. These two qualitative comparisons provide a framework for future validation measurements in order to further verify that the model behaves within plausible ranges for physiological and maximum-voluntary tasks. iv. Used as a test case for the ArtiSynth platform. The model was developed concurrently with the underlying simulation platform, ArtiSynth. It provided a sufficiently challenging test case to advance the modeling fidelity of ArtiSynth to be appropriate for the oral, pha- ryngeal, and laryngeal systems; in particular, with respect to hard-soft tissue coupling and contact. Inverse techniques for hard/soft muscle-tissue models i. Formulated trajectory-tracking for muscle-activated dynamic FEM with constraints. We extended inverse trajectory-tracking techniques, previously proposed for articulated rigid-body models of limbs, for a generalized forward dynamics of rigid-body and FEM mod- els within the ArtiSynth platform. The inverse solver is appropriate for complex muscular-hydrostat structures, modeled with dynamic FEM, 119 6.1. Dissertation Contributions as well as coupled skeletal/muscle-tissue structures. The inverse tech- niques provide a significant added-value to the simulation toolset allow- ing for systematic investigation of muscle recruitment, given a partic- ular desired kinematic trajectory, as opposed to trial-and-error based muscle tuning. ii. Formulated novel target parameters: constraint forces and stiffness. We have extended the inverse formulation to include tar- get parameters in additional to kinematics, providing a comprehensive toolset for exploring motor redundancy in motor tasks involving dy- namic constraints or stiffness. Constraint force magnitude targets are used to select muscle activations to increase forces within constrained DOF, such as to generate bite forces in the jaw model. Stiffness tar- gets are used to co-activate antagonist muscle groups to increase system stiffness which is important in many motor tasks. iii. Predicted beam and tongue muscle activations consistent with muscular-hydrostat theory. We used the inverse simulation toolset to automatically compute muscle activations needed to move and de- form muscular-hydrostat models, including a beam model and a phys- iologically based tongue model. The simulation results are consistent with theoretical proposals and provide a significant contribution over previous work by systematically evaluating muscular-hydrostat motor recruitment without a priori selection of muscle patterns. iv. Predicted plausible muscle activations for lateral jaw move- ment. We applied the inverse simulation toolset to automatically com- pute jaw muscle activations in lateral jaw movement. The simulation results are consistent with published jaw physiology. We also proposed a framework to analyze lateral pterygoid recruitment with respect to mechanical advantage in comparison to published EMG studies. 120 6.2. Directions Application to the analysis of segmental jaw surgery i. Created models of segmental jaw surgery with/without re- construction. We developed models of segmental jaw resection and reconstruction through structural alterations to a model of the normal jaw system. The procedure illustrates the application of a biomechan- ical model to biomedical treatment and rehabilitation planning. ii. Compared mechanical basis of functional deficits between mod- els. We simulated movement and bite force deficits that are observed clinically consequent to jaw resection. By comparing a model with and without mandibular grafting we illustrate the biomechanical un- derpinnings of lateral deviation and unstable jaw fragment movements. The simulations provide a foundation for future studies on the effect of scarring and the design of mandible and dental prostheses. iii. Applied inverse toolset to predict muscle forces to compensate for deficits. We applied the inverse simulation toolset to automati- cally predict muscle forces in the model to compensate for functional deficits. Results show that atypical co-activation of antagonist mus- cles was capable of aligning and stabilizing a unilaterally resected jaw model. The technique provides a methodology for analyzing altered and compromised anatomical systems and could inform rehabilitative muscle strengthening or motor retraining. 6.2 Directions Short term directions for the research components of this dissertation are discussed at the end of each research chapter. Here we discuss broader direc- tions for the future potential of integrating computational tools, including digital models and biomechanics simulation, in medical research and prac- tice. Protocols for clinical data collection The modeling work presented in this dissertation has focused on mechanisms to permit the integration 121 6.2. Directions of experimental data with biomechanics models. In particular, the inverse trajectory-tracking methods described in Chapter 4 suggest that simultane- ous recordings of detailed kinematic data and muscle activity are useful to evaluate model-based muscle force predictions. Also, Section 3.3.2 illustrates that maximum voluntary effort experiments, such as maximum tongue- palate pressure, provide a means to calibrate a model’s maximum muscle force parameters. In general, measurements of external contact forces dur- ing motor tasks provide a means to assess simulated internal muscle forces in a model, which are impossible to directly measure in humans. Finally, our efforts to model the mechanical consequences of jaw surgery revealed that the properties of scar tissue play a significant functional role and are not well described in the literature. In addition to detailed measurement of scar volume and mechanical properties, we suggest that clinical assessment of the effect of scar stiffness on jaw mobility would be helpful. Applying known external forces to the patient’s lower jaw in a number of directions and mag- nitudes, while simultaneously measuring jaw movement, would provide data on the effective jaw stiffness induced due to post-operative scarring. These specific experimental measurements are examples of what we see as a promis- ing future promising direction for iterative refinement of both modeling and experimental recording techniques as each process informs the other. Clinical workflow management Clinical and surgical workflows for mul- tistage medical procedures, such as oral and maxillofacial reconstructive surgery, are complex and involve a large volume of patient data from multi- ple sources. Current clinical workflows are organized in ways that can lead to inefficient flow of information and reduced team communication has been associated with higher rates of surgical morbidity and mortality [44]. Digital models present a promising and tangible mechanism to synthesize informa- tion for a single patient across different data sources and imaging modalities as well as to integrate information and statistics across patient populations. Digital modeling tools also provide highly visual artifacts for maintaining communication between treatment stakeholders, especially the patient. An informed patient is more likely to be a satisfied patient and digital mod- 122 6.2. Directions els have the potential to help keep a patient well-informed throughout the treatment process. Treatment alternatives Biomechanical models have the potential to im- pact treatment planning through quantitative analysis of trade-offs between treatment alternatives. For example, radiation therapy in head and neck cancer can cause damage and scarring in healthy tissue, which could po- tentially cause more adverse biomechanical alterations than tissue resection and reconstruction. With accurate, validated, and patient-specific models, these alternative treatment pathways could be evaluated quantitatively in terms of expected biomechanical outcomes. Tissue regeneration Reconstructive medicine can benefit from compu- tational biomechanics because grafting and transplanting tissue has both aesthetic and mechanical motivations. Regenerative medicine, including synthetic tissues and tissue regeneration, has the potential to vastly in- crease our capacity to alter, reshape, and restore anatomical structures and the need for computer tools to plan and guide such activities will increase accordingly [32]. Biological modeling science Advances in biomedical technology, includ- ing models and simulation techniques, are driving an emerging scientific discipline that integrates biological and medical knowledge with engineer- ing expertise. This discipline, the science of biological modeling, requires a broad knowledge base in biological science, empirical techniques, com- puter methods, and mathematics, as demonstrated by the breadth of related work presented in this dissertation. Students undertaking such study will likely be engineers with biological/medical interests, and biologists, doctors, and dentists with a strong technical background. Realizing the potential of this emerging field requires developing both formal curricula and dedi- cated research communities. Paramount to its success will be an open and standards-based research environment that will allow researchers to com- pare and synthesize new models, share common data-sets for both normal 123 6.3. Concluding Remarks subjects and patients, and recreate published modeling and data-processing techniques in order to verify and validate results. 6.3 Concluding Remarks To conclude, this dissertation has presented new models and simulation tech- niques to analyze jaw and tongue movements, muscle forces, and biomechan- ics. We have applied these modeling techniques in the biomedical domain and are working toward the development of new tools for dentists, doctors, and surgeons to better diagnose and treat orofacial and upper airway dys- function. This work constitutes a central component of a larger project developing models of the entire upper airway. The research has been carried out within an interdisciplinary team and has lead to a number of ongoing collaborations and projects. The models and simulation platform are openly available on the web for use by other researchers. 124 Bibliography [1] S. Abd-El-Malek. The part played by the tongue in mastication and deglutition. Journal of Anatomy, 89(2):250–254, 1955. → pages 14 [2] R. Adell, B. Svensson, and T. B̊agenholm. Dental rehabilitation in 101 primarily reconstructed jaws after segmental resections- possibilities and problems. An 18-year study. Journal of Cranio- Maxillofacial Surgery, 36(7):395–402, 2008. → pages 95 [3] M. J. Aftosmis, M. J. Berger, and J. E. Melton. Robust and efficient cartesian mesh generation for component-based geometry. AIAA Journal, 36(6):952–960, 1998. → pages 181 [4] R. Ait-Haddou, A. Jinha, W. Herzog, and P. Binding. Analysis of the force-sharing problem using an optimization model. Mathemati- cal Biosciences, 191(2):111–122, 2004. → pages 76, 115 [5] F. Alipour, I. R. Titze, E. Hunter, and N. Tayama. Active and pas- sive properties of canine abduction/adduction laryngeal muscles. Journal of Voice, 19(3):350–359, 2005. → pages 24, 32 [6] J. Allard, S. Cotin, F. Faure, P.-J. Bensoussan, F. Poyer, C. Duriez, H. Delingette, and L. Grisoni. SOFA – an open source framework for medical simulation. In Medicine Meets Virtual Reality (MMVR ’07), pages 13–18, 2007. → pages 28 [7] A. E. Anderson, B. J. Ellis, and J. A. Weiss. Verification, validation and sensitivity studies in computational biomechanics. Computer Methods in Biomechanics and Biomedical Engineering, 10(3):171– 184, 2007. → pages 65 [8] F. C. Anderson and M. G. Pandy. Dynamic optimization of human walking. Journal of Biomechanical Engineering, 123(5):381–390, 2001. → pages 38 125 Bibliography [9] F. C. Anderson and M. G. Pandy. Static and dynamic optimization solutions for gait are practically equivalent. Journal of Biomechan- ics, 34(2):153–161, 2001. → pages 38 [10] M. Anitescu and G. D. Hart. A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and fric- tion. International Journal for Numerical Methods in Engineering, 60(14):2335–2371, 2004. → pages 176 [11] M. Anitescu and F. A. Potra. A time-stepping method for stiff multibody dynamics with contact and friction. International Jour- nal for Numerical Methods in Engineering, 55(7):753–784, 2002. → pages 175 [12] ANSYS, Canonsburg, PA. ANSYS FEM simulation software. http: //www.ansys.com, 2010. → pages 27, 46 [13] G. Ardran and F. Kemp. The protection of the laryngeal airway during swallowing. The British Journal of Radiology, 25(296):406– 416, 1952. → pages 15 [14] P. Badin and A. Serrurier. Three dimensional modelling of speech organs: articulatory data and models. IEICE Technical Report, 106 (177):29–34, 2006. → pages 65, 90 [15] P. Badin, G. Bailly, L. Reveret, M. Baciu, C. Segebarth, and C. Savariaux. Three-dimensional linear articulatory modeling of tongue, lips and face, based on MRI and video images. Journal of Phonetics, 30(3):533–553, 2002. → pages 16 [16] A. Baker. A biomechanical model of the human tongue for under- standing speech production and other lingual behaviors. PhD thesis, Department of Linguistics, University of Arizona, 2008. → pages 33 [17] T. Belytschko, W. K. Liu, and B. Moran. Nonlinear finite element analysis for continua and structures. Wiley, New York, 2000. → pages 27 [18] N. A. Bernstein. The co-ordination and regulation movements. Perg- amon Press, Oxford, 1967. → pages 14 [19] D. A. Berry, J. B. Moon, and D. P. Kuehn. A finite element model of the soft palate. The Cleft Palate-Craniofacial Journal, 36(3):217– 223, 1999. → pages 34 126 Bibliography [20] R. Bibb, D. Eggbeer, P. Evans, A. Bocca, and A. Sugar. Rapid man- ufacture of custom-fitting surgical guides. Rapid Prototyping Jour- nal, 15(5):346–354, 2009. → pages 40 [21] M. Birch and P. Srodon. Biomechanical properties of the human soft palate. The Cleft Palate-Craniofacial Journal, 46(3):268–274, 2009. → pages 24 [22] J. Bonet and R. D. Wood. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press, Cambridge, UK, 2008. → pages 27, 47 [23] A. Bothorel. Cinéradiographie des voyelles et consonnes du franccais : recueil de documents synchronisés pour quatre sujets: vues latérales du conduit vocal, vues frontales de l’orifice labial, données acous- tiques. Institut de phonétique de Strasbourg, Strasbourg, France, 1986. → pages 15 [24] P. Boyle and B. Levin. World cancer report 2008. IARC Press, online, 2008. → pages 5 [25] B. Brogliato and V. Acary. Numerical methods for nonsmooth dy- namical systems. Springer, Berlin, 2008. → pages 185 [26] I. E. Brown and G. E. Loeb. A reductionist approach to creating and using neuromusculoskeletal models. In J. M. Winters and P. E. Crago, editors, Biomechanics and Neural Control of Posture and Movement, pages 148–163. Springer-Verlag, New York, 2000. → pages 36 [27] I. E. Brown, E. J. Cheng, and G. E. Loeb. Measured and modeled properties of mammalian skeletal muscle. II. The effects of stimulus frequency on force-length and force-velocity relationships. Journal of Muscle Research and Cell Motility, 20(7):627–643, 1999. → pages 32 [28] S. Buchaillard, M. Brix, P. Perrier, and Y. Payan. Simulations of the consequences of tongue surgery on tongue mobility: Implications for speech production in post-surgery conditions. International Journal of Medical Robotics and Computer Assisted Surgery, 3(3):252–261, 2007. → pages 40 [29] S. Buchaillard, P. Perrier, and Y. Payan. A biomechanical model of cardinal vowel production: Muscle activations and the impact of 127 Bibliography gravity on tongue positioning. Journal of the Acoustical Society of America, 126(4):2033–2051, 2009. → pages 21, 27, 31, 33, 36, 44, 45, 46, 48, 50, 57, 59, 66, 93 [30] M. Bucki, M. Nazari, and Y. Payan. Finite element speaker-specific face model generation for the study of speech production. Computer Methods in Biomechanics and Biomedical Engineering., 13(4):459– 467, 2010. → pages 27, 29, 48, 61 [31] T. A. Burnett, E. A. Mann, S. A. Cornell, and C. L. Ludlow. Laryn- geal elevation achieved by neuromuscular stimulation at rest. Jour- nal of Applied Physiology, 94(1):128–134, 2003. → pages 16, 63 [32] D. L. Butler, S. A. Goldstein, R. E. Guldberg, X. E. Guo, R. Kamm, C. T. Laurencin, L. V. McIntire, V. C. Mow, R. M. Nerem, R. L. Sah, L. J. Soslowsky, R. L. Spilker, and R. T. Tran- quillo. The impact of biomechanics in tissue engineering and regen- erative medicine. Tissue Engineering Part B: Reviews, 15(4):477, 2009. → pages 123 [33] E. R. Carlson and R. E. Marx. Mandibular reconstruction using cancellous cellular bone grafts. Journal of Oral and Maxillofacial Surgery, 54(7):889–897, 1996. → pages 95 [34] Cartsens, Munich, Germany. AG500 electromagnetic articulograph. http://www.articulograph.de, 2010. → pages 17 [35] M. C. Cavusoglu, T. Goktekin, and F. Tendick. GiPSi: a frame- work for open source/open architecture software development for organ-level surgical simulation. IEEE Transactions on Information Technology in Biomedicine, 10(2):312–322, 2006. → pages 28 [36] M. Chabanas, V. Luboz, and Y. Payan. Patient specific finite ele- ment model of the face soft tissues for computer-assisted maxillofa- cial surgery. Medical Image Analysis, 7(2):131–151, 2003. → pages 40, 62 [37] D. T. Corr, C. L. Gallant-Behm, N. G. Shrive, and D. A. Hart. Biomechanical behavior of scar tissue and uninjured skin in a porcine model. Wound Repair and Regeneration, 17(2):250–259, 2009. → pages 104 128 Bibliography [38] R. W. Cottle, J.-S. Pang, and R. E. Stone. The Linear Complemen- tarity Problem. Academic Press, Boston, MA, 1992. → pages 77, 177 [39] R. D. Crowninshield and R. A. Brand. A physiologically based crite- rion of muscle force prediction in locomotion. Journal of Biomechan- ics, 14(11):793–801, 1981. → pages 37 [40] D. A. Curtis, O. Plesh, A. J. Miller, T. A. Curtis, A. Sharma, R. Schweitzer, R. L. Hilsinger, L. Schour, and M. Singer. A com- parison of masticatory function in patients with or without recon- struction of the mandible. Head & Neck, 19(4):287–296, 1997. → pages 95, 96 [41] J. Dang and K. Honda. Construction and control of a physiological articulatory model. Journal of the Acoustical Society of America, 115(2):853–870, 2004. → pages 33 [42] Dassault Systèmes, Vélizy-Villacoublay, France. Simulia. http:// www.simulia.com, 2010. → pages 27 [43] Dassault Systèmes, Vélizy-Villacoublay, France. Solidworks. http: //www.solidworks.com/, 2010. → pages 25 [44] D. L. Davenport, W. G. Henderson, C. L. Mosca, S. F. Khuri, and R. M. Mentzer Jr. Risk-adjusted morbidity in teaching hospitals correlates with reported levels of communication and collaboration on surgical teams but not with scale measures of teamwork climate, safety climate, or working conditions. Journal of the American Col- lege of Surgeons, 205(6):778–784, 2007. → pages 122 [45] M. de Zee, P. M. Cattaneo, P. Svensson, T. K. Pedersen, B. Melsene, J. Rasmussen, and M. Dalstra. Prediction of the ar- ticular eminence shape in a patient with unilateral hypoplasia of the right mandibular ramus before and after distraction osteogenesis – a simulation study. Journal of Biomechanics, 42(8):1049–1053, 2009. → pages 33, 40 [46] S. L. Delp, F. C. Anderson, A. S. Arnold, P. Loan, A. Habib, C. T. John, E. Guendelman, and D. G. Thelen. OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE Transactions on Biomedical Engineering, 54(11):1940–1950, 2007. → pages 25, 30 129 Bibliography [47] R. L. Drake, W. Vogl, A. W. M. Mitchell, R. Tibbitts, and P. Richardson. Gray’s anatomy for students. Churchill Livingstone, Elsevier Inc., Philadelphia, PA, 2nd edition, 2010. → pages 2, 152, 153, 154, 155, 157, 158, 159, 161, 164, 165, 167, 170 [48] F. A. Duck. Physical properties of tissue: a comprehensive reference book. Academic Press, London, UK, 1990. → pages 24 [49] P. R. Eastwood, G. T. Allison, K. L. Shepherd, I. Szollosi, and D. R. Hillman. Heterogeneous activity of the human genioglossus muscle assessed by multiple bipolar fine-wire electrodes. Journal of Applied Physiology, 94(5):1849–1858, 2003. → pages 20 [50] H. Eckel, C. Sittel, P. Zorowka, and A. Jerke. Dimensions of the laryngeal framework in adults. Surgical and Radiologic Anatomy, 16 (1):31–36, 1994. → pages 26, 27, 168 [51] H. Edelsbrunner and E. P. Mücke. Simulation of simplicity: a tech- nique to cope with degenerate cases in geometric algorithms. ACM Transactions on Graphics, 9(1):66–104, 1990. → pages 181 [52] O. Engwall. Combining MRI, EMA and EPG measurements in a three-dimensional tongue model. Speech Communication, 41(2-3): 303–329, 2003. → pages 16, 17, 18 [53] A. Erdemir, S. McLean, W. Herzog, and A. J. van den Bogert. Model-based estimation of muscle forces exerted during movements. Clinical Biomechanics, 22(2):131–154, 2007. → pages 35, 39, 115 [54] Q. Fang, S. Fujita, X. Lu, and J. Dang. A model-based investigation of activations of the tongue muscles in vowel production. Acoustical Science and Technology, 30(4):277–287, 2009. ISSN 1346-3969. → pages 34 [55] A. Feldman. Once more on the equilibrium-point hypothesis (lambda model) for motor control. Journal of Motor Behavior, 18 (1):17, 1986. → pages 36, 47 [56] A. Feldman and M. Latash. Testing hypotheses and the advance- ment of science: recent attempts to falsify the equilibrium point hy- pothesis. Experimental Brain Research, 161(1):91–103, 2005. → pages 36 130 Bibliography [57] C. P. Fernandes, P. O. J. Glantz, S. A. Svensson, and A. Bergmark. A novel sensor for bite force determinations. Dental Materials, 19(2): 118–126, 2003. → pages 19 [58] B. R. Fink, R. W. Martin, and C. A. Rohrmann. Biomechanics of the Human Epiglottis. Acta Oto-Laryngologica, 87(3):554–559, 1979. → pages 169 [59] S. Fujita, J. Dang, N. Suzuki, and K. Honda. A computational tongue model and its clinical application. Oral Science International, 4(2):97–109, 2007. → pages 40 [60] L. M. Gallo, K. Fushima, and S. Palla. Mandibular helical axis path- ways during mastication. Journal of Dental Research, 79:1566–1572, 2000. → pages 14 [61] S. C. Gandevia, A. L. Hudson, R. B. Gorman, J. E. Butler, and A. De Troyer. Spatial distribution of inspiratory drive to the parasternal intercostal muscles in humans. The Journal of Physi- ology, 573(1):263, 2006. → pages 87 [62] J. M. Gerard, R. Wilhelms-Tricarico, P. Perrier, and Y. Payan. A 3D dynamical biomechanical tongue model to study speech motor control. Recent Research Developments in Biomechanics, 1:49–64, 2003. → pages 33 [63] J. M. Gerard, J. Ohayon, V. Luboz, P. Perrier, and Y. Payan. Non linear elastic properties of the lingual and facial tissues assessed by indentation technique: application to the biomechanics of speech production. Medical Engineering & Physics, 27(10):884–892, 2005. → pages 24, 47 [64] R. J. Gilbert and V. J. Napadow. Three-dimensional muscular ar- chitecture of the human tongue determined in vivo with diffusion tensor magnetic resonance imaging. Dysphagia, 20(1):1–7, 2005. → pages 23 [65] GNU General Public License. Blender. http://www.blender.org/, 2010. → pages 26 [66] H. Gomi and M. Kawato. Human arm stiffness and equilibrium- point trajectory during multi-joint movement. Biological Cybernet- ics, 76(3):163–171, 1997. → pages 36, 93 131 Bibliography [67] E. Guendelman, R. Bridson, and R. Fedkiw. Nonconvex rigid bodies with stacking. In ACM SIGGRAPH Papers, pages 871–878, 2003. → pages 184 [68] A. Hannam and A. McMillan. Internal organization in the human jaw muscles. Critical Reviews in Oral Biology & Medicine, 5(1):55, 1994. → pages 31, 158, 160 [69] A. G. Hannam. Current computational modelling trends in cran- iomandibular biomechanics and their clinical implications. Journal of Oral Rehabilitation, in press, 2010. → pages 39, 65 [70] A. G. Hannam and W. W. Wood. Medial pterygoid muscle activity during the closing and compressive phases of human mastication. American Journal of Physical Anthropology, 55(3):359–367, 1981. → pages 19 [71] A. G. Hannam, I. Stavness, J. E. Lloyd, and S. Fels. A dynamic model of jaw and hyoid biomechanics during chewing. Journal of Biomechanics, 41(5):1069–1076, 2008. → pages 33, 44, 45, 50, 52, 54, 56, 57, 65, 172 [72] A. G. Hannam, I. Stavness, J. E. Lloyd, S. Fels, A. J. Miller, and D. A. Curtis. A comparison of simulated jaw dynamics in models of segmental mandibular resection versus resection with alloplastic reconstruction. Journal of Prosthetic Dentistry, 104(3):191–198, 2010. → pages iii, 99 [73] M. Haraguchi, H. Mukohyama, D. J. Reisberg, and H. Taniguchi. Electromyographic activity of masticatory muscles and mandibular movement during function in marginal mandibulectomy patients. Journal of Medical and Dental Sciences, 50(4):257–264, 2003. → pages 95 [74] C. Head, D. Alam, J. A. Sercarz, J. T. Lee, J. D. Rawnsley, G. S. Berke, and K. E. Blackwell. Microvascular flap reconstruction of the mandible: a comparison of bone grafts and bridging plates for restoration of mandibular continuity. Otolaryngology-Head and Neck Surgery, 129(1):48–54, 2003. → pages 95 [75] W. D. Henshaw. Automatic grid generation. Acta numerica, 5:121– 148, 2008. → pages 28 132 Bibliography [76] O. Hidaka, M. Iwasaki, M. Saito, and T. Morimoto. Influence of clenching intensity on bite force balance, occlusal contact area, and average bite pressure. Journal of Dental Research, 78(7):1336, 1999. → pages 19 [77] D. A. Hidalgo and A. L. Pusic. Free-flap mandibular reconstruction: a 10-year follow-up study. Plastic and Reconstructive Surgery, 110 (2):438–449, 2002. → pages 95 [78] K. M. Hiiemae and J. B. Palmer. Tongue movements in feeding and speech. Critical Reviews in Oral Biology & Medicine, 14(6): 413, 2003. → pages 7, 13, 67 [79] K. M. Hiiemae, J. B. Palmer, S. W. Medicis, J. Hegener, B. Scott Jackson, and D. E. Lieberman. Hyoid and tongue surface movements in speaking and eating. Archives of Oral Biology, 47(1): 11–27, 2002. → pages 15, 65 [80] A. V. Hill. The heat of shortening and the dynamic constants of muscle. Philosophical Transactions of the Royal Society B: Biological Sciences, 126(843):136–195, 1938. → pages 32 [81] N. Hogan. Adaptive control of mechanical impedance by coactiva- tion of antagonist muscles. IEEE Transactions on Automatic Con- trol, 29(8):681–690, 1984. → pages 37, 75 [82] K. Hori, T. Ono, and T. Nokubi. Coordination of tongue pressure and jaw movement in mastication. Journal of Dental Research, 85 (2):187–191, 2006. → pages 18 [83] K. Hori, T. Ono, K. Tamine, J. Kondo, S. Hamanaka, Y. Maeda, J. Dong, and M. Hatsuda. Newly developed sensor sheet for mea- suring tongue pressure during swallowing. Journal of Prosthodontic Research, 53(1):28–32, 2009. → pages 18 [84] B. K. P. Horn, H. M. Hilden, and S. Negahdaripour. Closed-form solution of absolute orientation using orthonormal matrices. Journal of the Optical Society of America A, 5(7):1127–1135, 1988. → pages 26 [85] B. Y. Huang, T. Whittle, and G. M. Murray. Activity of inferior head of human lateral pterygoid muscle during standardized lateral jaw movements. Archives of Oral Biology, 50(1):49–64, 2005. → pages 20, 87 133 Bibliography [86] Y. Huang, D. P. White, and A. Malhotra. Use of computational modeling to predict responses to upper airway surgery in obstructive sleep apnea. The Laryngoscope, 117(4):648–653, 2007. → pages 3, 34 [87] A. L. Hudson, J. L. Taylor, S. C. Gandevia, and J. E. Butler. Cou- pling between mechanical and neural behaviour in the human first dorsal interosseous muscle. Experimental Physiology, 587(4):917–925, 2009. → pages 87 [88] T. J. R. Hughes. The finite element method: linear static and dy- namic finite element analysis. Dover Publications, New York, 2000. → pages 173, 185, 186 [89] E. J. Hunter and I. R. Titze. Refinements in modeling the passive properties of laryngeal soft tissue. Journal of Applied Physiology, 103(1):206, 2007. → pages 24 [90] E. J. Hunter, I. R. Titze, and F. Alipour. A three-dimensional model of vocal fold abduction/adduction. Journal of the Acoustical Society of America, 115:1747, 2004. → pages 34, 169 [91] C. Jeannin, P. Perrier, Y. Payan, A. Dittmar, and B. Grosgogeat. Tongue pressure recordings during speech using complete denture. Materials Science and Engineering: C, 28(5-6):835–841, 2008. → pages 18 [92] W. Jiang, H. Bo, G. YongChun, and N. LongXing. Stress distri- bution in molars restored with inlays or onlays with or without en- dodontic treatment: A three-dimensional finite element analysis. Journal of Prosthetic Dentistry, 103(1):6–12, 2010. → pages 40 [93] M. I. Jordan and D. E. Rumelhart. Forward models: Supervised learning with a distal teacher. Cognitive Science, 16(3):307–354, 1992. → pages 36 [94] P. Kahrilas, J. Logemann, S. Lin, and G. Ergun. Pharyngeal clear- ance during swallowing: a combined manometric and videofluoro- scopic study. Gastroenterology, 103(1):128, 1992. → pages 19 [95] M. Kawato, K. Furukawa, and R. Suzuki. A hierarchical neural- network model for control and learning of voluntary movement. Bio- logical Cybernetics, 57(3):169–185, 1987. → pages 36 134 Bibliography [96] W. M. Kier and K. K. Smith. Tongues, tentacles and trunks: the biomechanics of movement in muscular-hydrostats. Zoological Jour- nal of the Linnean Society, 83(4):307–324, 1985. → pages 81, 83, 84 [97] J. Koolstra and T. Van Eijden. Three-dimensional dynamical capa- bilities of the human masticatory muscles. Journal of Biomechanics, 32(2):145–152, 1999. → pages 101, 106 [98] J. Koolstra and T. Van Eijden. A method to predict muscle control in the kinematically and mechanically indeterminate human masti- catory system. Journal of Biomechanics, 34(9):1179–1188, 2001. → pages 38 [99] J. Koolstra and T. Van Eijden. Functional significance of the cou- pling between head and jaw movements. Journal of Biomechanics, 37(9):1387–1392, 2004. → pages 65, 106 [100] J. Koolstra and T. Van Eijden. Combined finite-element and rigid- body analysis of human jaw joint dynamics. Journal of Biomechan- ics, 38(12):2431–2439, 2005. → pages 33 [101] J. H. Koolstra and T. Van Eijden. The jaw open-close movements predicted by biomechanical modelling. Journal of Biomechanics, 30 (9):943–950, 1997. → pages 33 [102] T. Kuratate, H. Yehia, and E. Vatikiotis-Bateson. Kinematics-based synthesis of realistic talking faces. In Auditory-Visual Speech Pro- cessing (AVSP ’98), page online, 1998. URL http://marcs.uws.edu. au/links/avisa/avsp98/papers/0028.pdf. → pages 14 [103] C. Lacoursière. Ghosts and machines: regularized variational meth- ods for interactive simulations of multibodies with dry frictional con- tacts. PhD thesis, Computer Science Deptartment, Umea University, Sweden, 2007. → pages 175, 178 [104] G. E. Langenbach and A. G. Hannam. The role of passive muscle tensions in a three-dimensional dynamic model of the human jaw. Archives of Oral Biology, 44(7):557–73, 1999. → pages 45 [105] R. J. Last. Anatomy Regional and Applied. J. & A. Churchill Ltd., London, UK, 4th edition, 1966. → pages 152 135 Bibliography [106] G. E. Loeb, I. E. Brown, and E. J. Cheng. A hierarchical foundation for models of sensorimotor control. Experimental Brain Research, 126(1):1–18, 1999. → pages 38 [107] J. A. Logemann. Role of the modified barium swallow in manage- ment of patients with dysphagia. Otolaryngology-Head and Neck Surgery, 116(3):335–338, 1997. → pages 15 [108] A. A. Lowe and W. D. Johnston. Tongue and jaw muscle activity in response to mandibular rotations in a sample of normal and anterior open-bite subjects. American Journal of Orthodontics, 76(5):565– 576, 1979. → pages 52 [109] C. L. Ludlow. Electrical neuromuscular stimulation in dysphagia: current status. Current Opinion in Otolaryngology & Head and Neck Surgery, 18(3):159—164, 2010. → pages 64 [110] D. S. Lundy, C. Smith, L. Colangelo, P. A. Sullivan, J. A. Loge- mann, C. L. Lazarus, L. A. Newman, T. Murry, L. Lombard, and J. Gaziano. Aspiration: cause and implications. Otolaryngology-Head and Neck Surgery, 120(4):474–478, 1999. → pages 3, 4 [111] C. Lunk and B. Simeon. Solving constrained mechanical systems by the family of newmark and α-methods. Journal of Applied Math- ematics and Mechanics (ZAMM), 86(10):772–784, 2006. → pages 174 [112] Y. K. Mariappan, K. J. Glaser, and R. L. Ehman. Magnetic res- onance elastography: A review. Clinical Anatomy, 23(5):497–511, 2010. → pages 24 [113] M. Marunick, B. Mathes, and B. Klein. Masticatory function in hemimandibulectomy patients. Journal of Oral Rehabilitation, 19(3): 289–295, 1992. → pages 96 [114] R. E. Marx, J. E. Cillo Jr, V. Broumand, and J. J. Ulloa. Outcome analysis of mandibular condylar replacements in tumor and trauma reconstruction: A prospective analysis of 131 cases with long-term follow-up. Journal of Oral and Maxillofacial Surgery, 66(12):2515– 2523, 2008. → pages 95, 96 [115] K. Matsuo and J. B. Palmer. Kinematic linkage of the tongue, jaw, and hyoid during eating and speech. Archives of Oral Biology, 55(4): 325–331, 2010. → pages 15 136 Bibliography [116] B. Mehta and S. Schaal. Forward models in visuomotor control. Journal of Neurophysiology, 88(2):942–953, 2002. → pages 95 [117] A. J. Miller. Craniomandibular muscles: their role in function and form. CRC-Press, Boca Raton, FL, 1991. → pages 13, 85, 87, 96, 115 [118] A. J. Miller and D. J. Van Daele. The pharynx and scar tissue. In OPAL Workshop, 2009. URL http://www.magic.ubc.ca/artisynth/ uploads/OPAL/Arthur Miller.pdf. → pages 41 [119] K. Miller, G. Joldes, D. Lance, and A. Wittek. Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation. Communications in Numerical Methods in Engineering, 23(2):121–134, 2007. → pages 174 [120] K. Miyawaki. A study on the musculature of the human tongue. Annual Bulletin RILP, University of Tokyo, 8:23–49, 1974. URL http://www.umin.ac.jp/memorial/rilp-tokyo/R08/R08 023.pdf. → pages 21 [121] K. Miyawaki, H. Hirose, T. Ushijima, and M. Sawashima. A prelimi- nary report on the electromyographic study of the activity of lingual muscles. Annual Bulletin RILP, University of Tokyo, 9:91–106, 1975. URL http://www.umin.ac.jp/memorial/rilp-tokyo/R09/R09 091.pdf. → pages 20 [122] S. Moisik. A Three-Dimensional Model of the Larynx and the La- ryngeal Constrictor Mechanism: Visually Synthesizing Pharyngeal and Epiglottal Articulations Observed in Laryngoscopy. PhD thesis, Department of Linguistics, University of Victoria, 2008. → pages 34 [123] W. Mollemans, F. Schutyser, N. Nadjmi, F. Maes, and P. Suetens. Predicting soft tissue deformations for a maxillofacial surgery plan- ning system: From computational strategies to a complete clinical validation. Medical Image Analysis, 11(3):282–301, 2007. → pages 40 [124] E. Moller. The chewing apparatus: an electromyographic study of the action of the muscles of mastication and its correlation to facial morphology. Acta Physiologica Scandinavica, Supplementum, 280: 1–229, 1966. → pages 19, 52 137 Bibliography [125] K. Montgomery, C. Bruyns, J. Brown, S. Sorkin, F. Mazzella, G. Thonier, A. Tellier, B. Lerman, and A. Menon. Spring: A gen- eral framework for collaborative, real-time surgical simulation. In Medicine Meets Virtual Reality (MMVR ’02), pages 23–26, 2002. → pages 28 [126] MSC Software Corporation, Santa Ana, California. ADAMS Multi- body Dynamics. http://www.mscsoftware.com, 2010. → pages 25 [127] M. Müller and M. Gross. Interactive virtual materials. Graphics interface (GI ’04), pages 239–246, 2004. → pages 173 [128] K. G. Munhall, M. Kawato, and E. Vatikiotis-Bateson. Coarticula- tion and physical models of speech production. In M. B. Broe and J. B. Pierrehumbert, editors, Papers in Laboratory Phonology V: Ac- quistion and the Lexicon, chapter 1, pages 9–28. Cambridge Univer- sity Press, Cambridge, UK, 2000. → pages 7, 44 [129] G. M. Murray, T. Orfanos, J. Y. Chan, K. Wanigaratne, and I. J. Klineberg. Electromyographic activity of the human lateral ptery- goid muscle during contralateral and protrusive jaw movements. Archives of Oral Biology, 44(3):269–85, 1999. → pages 20, 87 [130] T. Muto and M. Kanazawa. Positional change of the hyoid bone at maximal mouth opening. Oral Surgery, Oral Medicine, Oral Pathol- ogy, 77(5):451–455, 1994. → pages 15, 51, 65 [131] E. Muybridge. Muybridge’s Complete human and animal locomo- tion: all 781 plates from the 1887 Animal locomotion. Courier Dover Publications, 1979. → pages 13 [132] National Library of Medicine. Visible human dataset. http://www. nlm.nih.gov/research/visible/visible human.html, 2010. → pages 21 [133] M. Nazari, P. Perrier, M. Chabanas, and Y. Payan. Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation. Computer Methods in Biomechanics and Biomedi- cal Engineering, 13(4):469–482, 2010. → pages 34, 61 [134] Northern Digital Inc., Waterloo, ON. Optotrak optical position tracking system. http://www.ndigital.com, 2010. → pages 14 [135] Northern Digital Inc., Waterloo, ON. Wave electromagnetric posi- tion tracking system. http://www.ndigital.com, 2010. → pages 17 138 Bibliography [136] J. P. Okeson. Management of Temporomandibular Disorder and Occlusion. Mosby, Saint Louis, MO, 5th edition, 2003. → pages 152, 155 [137] T. Ono, K. Hori, and T. Nokubi. Pattern of tongue pressure on hard palate during swallowing. Dysphagia, 19(4):259–264, 2004. → pages 18 [138] T. Ono, K. Hori, K. Tamine, and Y. Maeda. Evaluation of tongue motor biomechanics during swallowing – from oral feeding models to quantitative sensing methods. Japanese Dental Science Review, 45 (2):65–74, 2009. → pages 67 [139] T. Orfanos, L. Sarinnaphakorn, G. M. Murray, and I. J. Klineberg. Placement and verification of recording electrodes in the superior head of the human lateral pterygoid muscle. Archives of Oral Biol- ogy, 41(5):493–503, 1996. → pages 20, 87 [140] D. K. Pai. Muscle mass in musculoskeletal models. Journal of Biomechanics, 43(11):2093–2098, 2010. → pages 31 [141] W. R. Panje and D. K. Holmes. Mandibulectomy without recon- struction can cause sleep apnea. The Laryngoscope, 94(12):1591– 1594, 2009. → pages 5 [142] V. Parthasarathy, J. L. Prince, M. Stone, E. Z. Murano, and M. NessAiver. Measuring tongue motion from tagged cine-MRI us- ing harmonic phase (HARP) processing. Journal of the Acoustical Society of America, 121:491, 2007. → pages 16, 90 [143] A. Patel and R. Maisel. Condylar prostheses in head and neck cancer reconstruction. Archives of Otolaryngology-Head and Neck Surgery, 127(7):842, 2001. → pages 95 [144] Y. Payan and P. Perrier. Synthesis of VV sequences with a 2D biomechanical tongue model controlled by the Equilibrium Point Hypothesis. Speech Communication, 22(2-3):185–205, 1997. → pages 6, 33, 57, 58 [145] Y. Payan, P. Perrier, C. Vilain, and X. Pelorson. Finite element models of the tongue and velum for a physical understanding of sleep apnea syndrome. In Proceedings of Computer Methods in Biomechanics and Biomedical Engineering (BBE’01), page online, 2001. → pages 34 139 Bibliography [146] C. C. Peck, G. E. J. Langenbach, and A. G. Hannam. Dynamic sim- ulation of muscle and articular properties during human wide jaw opening. Archives of Oral Biology, 45(11):963–982, 2000. → pages 32, 33, 45, 59 [147] C. C. Peck, A. S. Sooch, and A. G. Hannam. Forces resisting jaw displacement in relaxed humans: a predominantly viscous phe- nomenon. Journal of Oral Rehabilitation, 29(2):151–160, 2002. → pages 116 [148] J. Perkell, M. Cohen, M. Svirsky, M. Matthies, I. Garabieta, and M. Jackson. Electromagnetic midsagittal articulometer systems for transducing speech articulatory movements. The Journal of the Acoustical Society of America, 92:3078, 1992. → pages 17 [149] J. S. Perkell. Physiology of speech production: Results and implica- tions of a quantitative cineradiographic study. Research Monograph No. 53. MIT Press, Cambridge, MA, 1969. → pages 6, 15 [150] J. S. Perkell. A physiologically-oriented model of tongue activity in speech production. PhD thesis, Speech Communication Group, Mas- sachusetts Institute of Technology, 1974. → pages 33 [151] P. Perrier, Y. Payan, M. Zandipour, and J. Perkell. Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study. Journal of the Acoustical Society of America, 114(3):1582–1599, 2003. → pages 33 [152] D. H. Perrott, H. Umeda, and L. B. Kaban. Costochondral graft construction/reconstruction of the ramus/condyle unit: long-term follow-up. International Journal of Oral and Maxillofacial Surgery, 23(6):321–328, 1994. → pages 95 [153] G. J. Petruzzelli, K. Cunningham, and D. Vandevender. Impact of mandibular condyle preservation on patterns of failure in head and neck cancer. Otolaryngology-Head and Neck Surgery, 137(5):717–721, 2007. → pages 95 [154] L. J. Pittman and E. F. Bailey. Genioglossus and intrinsic elec- tromyographic activities in impeded and unimpeded protrusion tasks. Journal of Neurophysiology, 101(1):276, 2009. → pages 20 140 Bibliography [155] Polhemus, Colchester, VT. Fastrak magnetic position tracking sys- tem. http://www.polhemus.com, 2010. URL http://www.polhemus. com. → pages 17 [156] J. Porrill, P. Dean, and J. V. Stone. Recurrent cerebellar architec- ture solves the motor-error problem. Proceedings of the Royal Soci- ety B: Biological Sciences, 271(1541):789–796, 2004. → pages 36 [157] F. A. Potra, M. Anitescu, B. Gavrea, and J. Trinkle. A linearly im- plicit trapezoidal method for integrating stiff multibody dynamics with contact, joints, and friction. International Journal for Numeri- cal Methods in Engineering, 66(7), 2006. → pages 175, 184 [158] J. L. Prince and J. M. Links. Medical Imaging Signals and Systems. Pearson Prentice Hall, Upper Saddle River, NJ, 2006. → pages 16, 22 [159] O. Rastadmehr, T. Bressmann, R. Smyth, and J. C. Irish. Increased midsagittal tongue velocity as indication of articulatory compensa- tion in patients with lateral partial glossectomies. Head & Neck, 30 (6):718–726, 2008. → pages 17 [160] O. Röhrle and A. J. Pullan. Three-dimensional finite element mod- elling of muscle forces during mastication. Journal of Biomechanics, 40(15):3363–3372, 2007. → pages 31, 33 [161] E. D. Roumanas, N. Garrett, K. E. Blackwell, E. Freymiller, E. Abe- mayor, W. K. Wong, J. Beumer, K. Fueki, W. Fueki, and K. K. Ka- pur. Masticatory and swallowing threshold performances with con- ventional and implant-supported prostheses after mandibular fibula free-flap reconstruction. Journal of Prosthetic Dentistry, 96(4):289– 297, 2006. → pages 95 [162] V. Sanguineti, R. Laboissière, and D. J. Ostry. A dynamic biome- chanical model for neural control of speech production. Journal of the Acoustical Society of America, 103(3):1615–1627, 1998. → pages 33, 34, 35, 36 [163] O. Schenk and K. Gärtner. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Generation Computer Systems, 20(3):475–487, 2004. → pages 179 141 Bibliography [164] P. Schiavone, E. Promayon, and Y. Payan. LASTIC: A Light Aspi- ration Device for in vivo Soft TIssue Characterization. In Biomedical Simulation (ISBMS ’10), volume Springer LNCS 5958, pages 1–10. Springer, 2010. → pages 24 [165] H. Schliephake, R. Schmelzeisen, R. Schnweiler, T. Schneller, and C. Altenbernd. Speech, deglutition and life quality after intraoral tumour resection: a prospective study. International Journal of Oral and Maxillofacial Surgery, 27(2):99 – 105, 1998. → pages 95 [166] R. Schmelzeisen, F. W. Neukam, T. Shirota, B. Specht, and M. Wichmann. Postoperative function after implant insertion in vascularized bone grafts in maxilla and mandible. Plastic and Re- constructive Surgery, 97(4):719, 1996. → pages 95, 96 [167] W. Schmidt-Nowara, A. Lowe, L. Wiegand, R. Cartwright, F. Perez- Guerra, S. Menn, et al. Oral appliances for the treatment of snoring and obstructive sleep apnea: a review. Sleep, 18(6):501–510, 1995. → pages 52 [168] W. S. Selbie, L. Zhang, W. S. Levine, and C. L. Ludlow. Using joint geometry to determine the motion of the cricoarytenoid joint. Jour- nal of the Acoustical Society of America, 103:1115–1127, 1998. → pages 15 [169] W. S. Selbie, S. L. Gewalt, and C. L. Ludlow. Developing an anatomical model of the human laryngeal cartilages from magnetic resonance imaging. Journal of the Acoustical Society of America, 112:1077, 2002. → pages 23, 168 [170] A. Serrurier and P. Badin. A three-dimensional articulatory model of the velum and nasopharyngeal wall based on MRI and CT data. Journal of the Acoustical Society of America, 123:2335–2355, 2008. → pages 63 [171] A. Seth and M. G. Pandy. A neuromusculoskeletal tracking method for estimating individual muscle forces in human movement. Journal of Biomechanics, 40(2):356–366, 2007. → pages 39 [172] A. A. Shabana. Dynamics of multibody systems. Cambridge Univer- sity Press, Cambridge, UK, 2005. → pages 25 142 Bibliography [173] M. Sherif, R. Gregor, L. Liu, R. Roy, and C. Hager. Correlation of myoelectric activity and muscle force during selected cat treadmill locomotion. Journal of Biomechanics, 16(9):691, 1983. → pages 19, 70 [174] J. R. Shewchuk. What is a good linear element? interpolation, con- ditioning, and quality measures. In Eleventh International Meshing Roundtable, pages 115–126, online, 2002. Sandia National Laborato- ries. → pages 28 [175] D. M. Shiller, R. Laboissiere, and D. J. Ostry. Relationship between jaw stiffness and kinematic variability in speech. Journal of Neuro- physiology, 88(5):2329–2340, 2002. → pages 36 [176] D. M. Shiller, G. Houle, and D. J. Ostry. Voluntary control of hu- man jaw stiffness. Journal of Neurophysiology, 94(3):2207–2217, 2005. → pages 88, 116 [177] T. Shinar, C. Schroeder, and R. Fedkiw. Two-way coupling of rigid and deformable bodies. In ACM SIGGRAPH/Eurographics Sympo- sium on Computer Animation (SCA ’08), 2008. → pages 176, 184 [178] H. Si. Tetgen: A quality tetrahedral mesh generator and a 3d delau- nay triangulator. http://tetgen.berlios.de/, 2009. → pages 28 [179] E. Sifakis, I. Neverov, and R. Fedkiw. Automatic determination of facial muscle activations from sparse motion capture marker data. ACM Transactions on Graphics, 24(3):417–425, July 2005. → pages 14, 34, 39, 70, 71, 73, 91 [180] M. P. T. Silva and J. A. C. Ambrosio. Sensitivity of the results produced by the inverse dynamic analysis of a human stride to per- turbed input data. Gait and Posture, 19(1):35–49, 2004. → pages 38, 91 [181] S. I. Silverman. Oral Physiology. Mosby, Saint Louis, MO, 1961. → pages 166 [182] B. C. Sonies, T. H. Shawker, T. E. Hall, L. H. Gerber, and S. B. Leighton. Ultrasonic visualization of tongue motion during speech. Journal of the Acoustical Society of America, 70:683, 1981. → pages 16 143 Bibliography [183] I. Stavness, A. G. Hannam, J. E. Lloyd, and S. Fels. An integrated, dynamic jaw and laryngeal model constructed from ct data. In Biomedical Simulation (ISBMS ’06), Springer LNCS, volume 4072, pages 169–177, 2006. → pages 45, 50 [184] I. Stavness, C. L. Ludlow, B. Chung, and S. Fels. Hyolaryngeal biomechanics modeling with intramuscular stimulation data. In OPAL Workshop, pages 73–78, 2009. → pages iii [185] I. Stavness, A. G. Hannam, J. Lloyd, and S. Fels. Predicting muscle patterns for hemimandibulectomy models. Computer Methods in Biomechanics and Biomedical Engineering, 13(4):483–491, 2010. → pages iii [186] I. Stavness, J. Lloyd, Y. Payan, and S. Fels. Coupled hard-soft tissue simulation with contact and constraints applied to jaw-tongue-hyoid dynamics. International Journal of Numerical Methods in Biomedi- cal Engineering, in press, 2010. → pages iii, 172 [187] M. Stone and A. Lundberg. Three-dimensional tongue surface shapes of English consonants and vowels. Journal of the Acoustical Society of America, 99(6):3728–3737, 1996. → pages 16 [188] M. Stone and E. Vatikiotis-Bateson. Trade-offs in tongue, jaw, and palatecontributions to speech production. Journal of Phonetics, 23 (1-2):81–100, 1995. → pages 18 [189] S. Sueda, A. Kaufman, and D. K. Pai. Musculotendon simulation for hand animation. In ACM SIGGRAPH Papers, pages 1–8, 2008. → pages 30, 39, 70, 71, 73 [190] S. Takano and K. Honda. An MRI analysis of the extrinsic tongue muscles during vowel production. Speech Communication, 49(1): 49–58, 2007. → pages 16, 23, 162, 163 [191] S. Takano, K. Honda, and K. Kinoshita. Observation of cricothyroid joint motion using 3d high-resolution mri. In Voice Physiology and Biomechanics (ICVPB ’04), page online, 2004. → pages 16 [192] H. Takemoto. Morphological analyses of the human tongue muscu- lature for three-dimensional modeling. Journal of Speech, Language, and Hearing Research, 44(1):95, 2001. → pages 21, 162 144 Bibliography [193] J. Teran, E. Sifakis, S. Blemker, V. Ng-Thow-Hing, C. Lau, and R. Fedkiw. Creating and simulating skeletal muscle from the visible human data set. IEEE Transactions on Visualization and Computer Graphics, 11(3):317–28, 2005. → pages 30 [194] D. Terzopoulos and K. Waters. Analysis and synthesis of facial im- age sequences using physical andanatomical models. IEEE Transac- tions on Pattern Analysis and Machine Intelligence, 15(6):569–579, 1993. → pages 34 [195] E. Teubner, D. Rohner, A. Deak, A. Lorenzon, and C. Marinello. Navigated implant placement using a bone supported CT-guided surgical template: a case report. Schweizer Monatsschrift fur Zahn- medizin, 119(12):1211–1234, 2009. → pages 40 [196] D. G. Thelen and F. C. Anderson. Using computed muscle control to generate forward dynamic simulations of human walking from experimental data. Journal of Biomechanics, 39(6):1107–1115, 2006. → pages 39, 70, 92 [197] E. Todorov and M. I. Jordan. Optimal feedback control as a theory of motor coordination. Nature Neuroscience, 5(11):1226–1235, 2002. → pages 38 [198] E. Todorov, W. Li, and X. Pan. From task parameters to motor synergies: A hierarchical framework for approximately optimal con- trol of redundant manipulators. Journal of Robotic Systems, 22(11): 691–710, 2005. → pages 38 [199] M. Tuijta, J. H. Koolstraa, F. Lobbezoob, and M. Naeijeb. Differ- ences in loading of the temporomandibular joint during opening and closing of the jaw. Journal of Biomechanics, 43(6):1048–1054, 2010. → pages 33 [200] M. L. Urken, D. Buchbinder, H. Weinberg, C. Vickery, A. Sheiner, R. Parker, J. Schaefer, P. Som, A. Shapiro, W. Lawson, et al. Func- tional evaluation following microvascular oromandibular reconstruc- tion of the oral cancer patient: a comparative study of reconstructed and nonreconstructed patients. The Laryngoscope, 101(9), 1991. → pages 96 [201] Y. Utanohara, R. Hayashi, M. Yoshikawa, M. Yoshida, K. Tsuga, and Y. Akagawa. Standard values of maximum tongue pressure 145 Bibliography taken using newly developed disposable tongue pressure measure- ment device. Dysphagia, 23(3):286–290, 2008. → pages 18, 59 [202] D. Van Daele, A. Perlman, and M. Cassell. Intrinsic fibre architec- ture and attachments of the human epiglottis and their contribu- tions to the mechanism of deglutition. Journal of Anatomy, 186(Pt 1):1, 1995. → pages 63, 169 [203] T. Van Eijden, J. Koolstra, and P. Brugman. Architecture of the human pterygoid muscles. Journal of Dental Research, 74(8):1489, 1995. → pages 21, 32 [204] T. Van Eijden, J. H. Koolstra, and P. Brugman. Three-dimensional structure of the human temporalis muscle. The Anatomical Record, 246(4):565–572, 1996. → pages [205] T. Van Eijden, J. Korfage, and P. Brugman. Architecture of the hu- man jaw-closing and jaw-opening muscles. The Anatomical Record, 248(3):464–474, 1997. → pages 21, 31, 32, 59, 160 [206] E. Vatikiotis-Bateson and D. J. Ostry. An analysis of the dimen- sionality of jaw motion in speech. Journal of Phonetics, 23:101–117, 1995. → pages 14 [207] Vicon Motion Systems, Los Angeles, CA. Vicon MX. http://www. vicon.com, 2010. URL http://www.vicon.com. → pages 14 [208] Visual Computing Lab of the Italian National Research Council. Meshlab. http://meshlab.sourceforge.net/, 2010. → pages 26 [209] N. Wakabayashi, M. Ona, T. Suzuki, and Y. Igarashi. Nonlinear finite element analyses: Advances and challenges in dental applica- tions. Journal of dentistry, 36(7):463–471, 2008. → pages 40 [210] A. Waltimo and M. Kononen. A novel bite force recorder and max- imal isometric bite force values for healthy young adults. European Journal of Oral Sciences, 101(3):171–175, 1993. → pages 19, 111, 116 [211] J. Weiss and S. Maas. FEBio - a finite element tool for biomechan- ics, 2007. URL http://imechanica.org/node/1897. → pages 28 [212] J. A. Weiss, B. N. Maker, and S. Govindjee. Finite element imple- mentation of incompressible, transversely isotropic hyperelasticity. 146 Computer Methods in Applied Mechanics and Engineering, 135(1-2): 107–128, 1996. → pages 31, 67 [213] A. Woda, P. Pionchon, and S. Palla. Regulation of mandibular pos- tures: mechanisms and clinical implications. Critical Reviews in Oral Biology & Medicine, 12(2):166, 2001. → pages 104 [214] W. W. Wood, A. G. Hannam, and K. Takada. The electromyo- graphic activity of the inferior part of the human lateral pterygoid muscle during clenching and chewing. Archives of Oral Biology, 31 (4):245–253, 1986. → pages 19 [215] G. T. Yamaguchi, D. W. Moran, and J. Si. A computationally ef- ficient method for solving the redundant problem in biomechanics. Journal of Biomechanics, 28(8):999–1005, 1995. → pages 37 [216] F. Zajac. Muscle and tendon: properties, models, scaling, and ap- plication to biomechanics and motor control. Critical Reviews in Biomedical Engineering, 17(4):359, 1989. → pages 32, 71, 91 [217] F. Zhang, C. C. Peck, and A. G. Hannam. Mass properties of the human mandible. Journal of Biomechanics, 35(7):975–978, 2002. → pages 26, 27, 31 147 Appendix A List of Publications Many of the contributions and ideas described in this dissertation have been previously presented in publications and oral presentations, which are listed here. A.1 Journal Publications Ian Stavness, John Lloyd, Yohan Payan, and Sidney Fels. Coupled hard- soft tissue simulation with contact and constraints applied to jaw- tongue-hyoid dynamics. International Journal of Numerical Methods in Biomedical Engineering, in press, 2010. Ian Stavness, Alan Hannam, John Lloyd, and Sidney Fels. Predicting muscle patterns for hemimandibulectomy models. Computer Methods in Biomechanics and Biomedical Engineering, 13(4):483–491, 2010. Alan G. Hannam, Ian Stavness, John E. Lloyd, Sidney Fels, Don Curtis, and Art Miller. A comparison of simulated jaw dynamics in mod- els of segmental mandibular resection versus resection with alloplastic reconstruction. Journal of Prosthetic Dentistry, 104(3):191–198, 2010. Sandro Palla, W.L. Hylander, A.S. McMillan, E.W.N. Lam, M. Watanabe, G.E.J. Langenbach, I. Stavness, C.C. Peck, and A.G. Hannam. From movement to models: a tribute to Professor Alan Hannam. Journal of Orofacial Pain, 22(4):307–316, 2008. 148 A.2. Conference Publications A.2 Conference Publications Ian Stavness, John Lloyd, Yohan Payan, and Sidney Fels. Dynamic hard-soft tissue models for orofacial biomechanics. In ACM Siggraph Talks, 2010. Ian Stavness, Christy Ludlow, Bethany Chung, and Sidney Fels. Hyola- ryngeal biomechanics modeling with intramuscular stimulation data. In OPAL Workshop, pages 73–78, 2009. Sidney Fels, Ian Stavness, and others. Advanced Tools for Biomechan- ical Modeling of the Oral, Pharyngeal, and Laryngeal Complex. In Proc of 4th International Symposium on Biomechanics, Healthcare and Information Science (ISBHIS ’09), 2009. Ian Stavness, Alan G. Hannam, John E. Lloyd, and Sidney Fels. Towards predicting biomechanical consequences of jaw reconstruction. In Proc of 30th IEEE Engineering in Medicine and Biology Conference (EMBC ’08), pages 4567–4570, 2008. Ian Stavness, Alan G. Hannam, John E. Lloyd, Eric Vatikiotis-Bateson, and Sidney Fels. Tools for predicting biomechanical consequences of alterations to orofacial anatomy. In Proc of 3rd International Sympo- sium on Biomechanics, Healthcare and Information Science (ISBHIS ’08), 2008. Sidney S. Fels, John E. Lloyd, Ian Stavness, Alan Hannam, and Eric Vatikiotis-Bateson. ArtiSynth: A 3d biomechanical simulation toolkit for modeling anatomical structures. In Journal of the Society for Simulation in Healthcare, volume 2, page 148, 2007. A.3 Research Visits October-November 2009 with Dr. Chris Peck and Dr. Greg Murray Jaw Function Research Unit, Westmead Hospital, University of Syd- ney, Australia 149 A.4. Research Talks October 2008 with Dr. Christy Ludlow Laryngeal and Speech Section, The National Institutes of Health Bethesda, MD A.4 Research Talks October 2009 – Upper Airway Biomechanics Modeling • Sleep Research Meeting, Woolcock Institute, Sydney, Australia • Prince of Wales Medical Research Institute, Sydney, Australia • MARCS Auditory Lab, University of Western Sydney, Australia July 2008 – Dynamic Jaw and Tongue Modeling with ArtiSynth • Orofacial Pain and Motor Control: The History of Three Giants in Orofacial Neurosciences Symposium, Toronto June 2008 – Biomechanics and Motor Control Simulation in ArtiSynth • Department of Biomedical Sciences, University of Maryland • Physical Medicine and Rehabilitation, John Hopkins Hospital • Laryngeal and Speech Section, The National Institutes of Health May 2007 – Biomechanical Modeling of Jaw-Hyoid Anatomy with ArtiSynth • Language, Mind and Brain Centre, McGill University, Montreal A.5 Master’s Publications Alan G. Hannam, Ian Stavness, John E. Lloyd, and Sidney Fels. A dy- namic model of jaw and hyoid biomechanics during chewing. Journal of Biomechanics, 41(5):1069–1076, 2008. Ian Stavness, Alan G. Hannam, John E. Lloyd, and Sidney Fels. An inte- grated, dynamic jaw and laryngeal model constructed from ct data. In Proceedings of the International Symposium on Biomedical Simulation (ISBMS 06), Springer LNCS, 4072:169177, 2006. Sidney Fels, John Lloyd, Kees van den Doel, Ian Stavness, and Eric Vatikiotis-Bateson. Developing physically-based, dynamic vocal tract 150 A.6. Additional Publications models using ArtiSynth. In Proceedings of the 7th International Sem- inar on Speech Production (ISSP 06), pages 419426, 2006. A.6 Additional Publications Ian Stavness, Billy Lam, and Sidney Fels. pCubee: A perspective- corrected handheld cubic display. In ACM CHI ’10: Proceedings of the SIGCHI conference on human factors in computing systems, pages 1381–1390, 2010. Billy Lam, Ian Stavness, Ryan Barr, and Sidney Fels. Interacting with a personal cubic 3d display. In ACM Multimedia Technical Demonstra- tions, pages 959–960, 2009. Ian Stavness and Sidney Fels. Cubee: thinking inside the box. In ACM Siggraph Emerging Technologies, page 5, 2006. Ian Stavness and Sidney Fels. Cubee: a cubic 3D display for physics- based interaction. In ACM Siggraph Sketches, page 165, 2006. Ian Stavness, Jen Gluck, Leah Vilhan, and Sidney Fels. The MU- SICtable: a map-based ubiquitous system for social interaction with a digital music collection. In Proceedings of International Conference in Entertainment Computing (ICEC 05), Springer LNCS, (3711):291302, 2005. 151 Appendix B Head and Neck Anatomy This appendix provides an overview of head and neck anatomy focusing on descriptions of jaw, tongue, and laryngeal muscles. More detailed anatomical descriptions can be found elsewhere, including Drake et al. [47], Okeson [136], and Last [105]. Figure B.1 provides a mid-sagittal cut-away view of the upper airway. In this appendix, we describe the bony structures of the head, including the cranium, mandible, teeth, hyoid bone, and vertebrae, as well as soft-tissue structures, including the face, tongue, soft-palate, pharynx, larynx (thyroid, cricoid, arytenoid cartilages), and epiglottis. We review muscles of the jaw, hyoid, tongue, and upper airway that are relevant to the modeling work described in this dissertation. B.1 Bone Structures The head and neck skeleton includes the skull (cranium and mandible), hyoid bone, and vertebrae. B.1.1 Cranium The cranium is composed of a large number of fused bones. The frontal (fa- cial) areas of the cranium includes the frontal, nasal, zygomatic bones, and maxilla. The top and back areas of the skull are formed by the parietal and occipital bones respectively. The temporal bones form the lower side regions of the skull and include the mastoid, styloid, and zygomatic processes. Maxilla The maxilla forms the upper jaw and is fused with the surround- ing bony structure of the cranium. The maxilla extends from the floor of the 152 B.1. Bone Structures TongueLips Mandible Hyoid bone Epiglottis Trachea Thyroid cartilage Cricoid cartilage Vertebrae Soft palateHard palate Vocal folds Figure B.1: Sagittal cross-section of the upper airway showing the skull, lips, tongue, palate, pharynx, and larynx. c© Elsevier (2010), Drake et al. [47], adapted with permission. Parietal Temporal Mandible Maxilla Frontal Occipital Mastoid Process Styloid Process Zygomatic Process Zygomatic Nasal Figure B.2: Lateral view of the of the skull showing the mandible and the fused bones of the cranium. c© Elsevier (2010), Drake et al. [47], adapted with permission. 153 B.1. Bone Structures Mylohyoid line Body Ramus Condylar head Condylar neck Condylar process Coronoid process Figure B.3: Medial and lateral views of the mandible showing and the as- cending ramus forming the coronoid and condylar processes. c© Elsevier (2010), Drake et al. [47], adapted with permission. nasal cavity and lower border of the orbit to the hard palate and alveolar process of the upper dental arch (see Figure B.2). Temporal bone The temporal bone serves as the articular surface for the mandibular condyle. The posterior area of the temporal bone is concave and forms the articular fossa (colloquially referred to as the “jaw joint socket” due to its concave shape). Anterior to the fossa is the articular eminence, a convex bony process that is the pathway for the condyle during jaw opening and protrusion (see Figure B.15). The articular eminence is designed to sustain large joint loads as it consists of thick, dense bone. B.1.2 Mandible The mandible, pictured in Figure B.3, is a u-shaped bone connected to the skull through muscles, ligaments, and contact at the TMJ. The frontal arch- shaped portion, or body of the mandible, forms the alveolar process for the lower teeth. The posterior portion of the mandible extends upward into the ascending ramus that forms two processes. The anterior coronoid process flattens mediolaterally and serves as the insertion site for the temporalis 154 B.1. Bone Structures Tooth crown Alveolar process Tooth root Gingiva Periodontal ligament Figure B.4: A tooth seated within the bony alveolar process and connected by the periodontal ligament. c© Elsevier (2003), Okeson [136] adapted with permission. Canines Premolars MolarsCanines Premolars Molars Incisors Incisors Figure B.5: The lower and upper dentition showing the molor, premolar, canine, and incisal teeth. c© Elsevier (2010), Drake et al. [47], adapted with permission. muscle, extending the mechanical advantage of this muscle for generating closing torque on the jaw. The posterior process of the ascending ramus, which has a convex shape, is termed the “condyle” and articulates with the cranium at the TMJ. B.1.3 Dentition The human dentition consists of an upper and lower arch each containing sixteen teeth, as shown in Figure B.5. The tooth body is divided into an exposed crown and an internal root that anchors it to the bone. The 155 B.2. Soft-Tissues and Muscles tooth root sits within bony sockets called alveolar processes in the mandible and maxilla and are attached by the periodontal ligament, as shown in Figure B.4. The periodontal support fibers provide force distribution and shock absorption for the teeth, as well as sensitive nerve endings to detect tooth loads. The upper dental arch is slightly wider and longer than its lower counterpart in order for tooth crowns to fit together during intercuspation. B.1.4 Hyoid bone The hyoid bone is a floating u-shaped bone located in the neck between the jaw and larynx. It is connected above to the mandible and skull by the submandibular muscles and below to the thyroid cartilage by the hyothyroid membrane along its length, as shown in Figure B.12. The hyoid serves as an anchor for the posterior muscle fibers of the tongue. B.1.5 Vertebrae The vertebrae are the ringed bones that enclose and protect the spinal cord. They extending in a chain from the base of the cranium to the pelvis. The upper section, named the cervical vertebrae (C1 thru C7), form the skele- ton of the neck and provide support to skull. Functionally, each vertebra articulates in sequence to enable movement of the head. The vertebrae also provide convenient anatomical landmarks as their dense structure is easily visible in x-ray and CT images (see Figure 2.2). B.2 Soft-Tissues and Muscles Soft-tissues in the upper airway surround the skeletal structures and gener- ate forces that articulate jaw and vocal tract. The upper airway soft-tissues are shown in Figure B.1 and include cartilage, muscles, ligaments, tendons, fascia (connective tissue), fat, mucous membrane that lines the inner surface of upper airway, and skin that covers the outer surface of the face. 156 B.2. Soft-Tissues and Muscles Figure B.6: Lateral view of the head showing the facial muscles. c© Elsevier (2010), Drake et al. [47], adapted with permission. B.2.1 Face The superficial muscles of the face are shown in Figure B.6. The cheeks and lips are important for shaping the upper airway. During mastication and swallowing the lips seal off the oral cavity and the buccinator muscle in the cheek works with the tongue to form and position the bolus during mastication. Lip shape and movement is also critical to the acoustics and visual appearance of speech production. The face includes a number of thin muscles located within the superficial fascia layer of facial tissue. Facial muscles originate from the skull or fascia and insert into the skin in order to deform the surface of the face and control facial expressions. All facial muscles are innervated by the facial nerve [cranial nerve VII]. 157 B.2. Soft-Tissues and Muscles Temporalis Deep Masseter Superficial Masseter Medial Pterygoid Figure B.7: Jaw closing muscles include the masseter, temporalis, and me- dial pterygoid muscles. c© Elsevier (2010), Drake et al. [47], adapted with permission. B.2.2 Jaw muscles The jaw muscles produce forces that move the jaw and generate tooth forces during chewing and clenching. The submandibular muscles are capable of opening the jaw and lifting the laryngeal complex during swallowing. The jaw muscles also stabilizes the mandible and prevents extreme jaw displace- ments through passive tension generated by muscle stretch. Jaw closing muscles Jaw closing muscles include the masseter, tempo- ralis, and medial pterygoid muscles and are pictured in Figure B.7. Each muscle has a large attachment area and can be further subdivided based on muscle fiber direction [68]. The masseter muscle is typically divided into su- perficial and deep parts; however, the muscle is further compartmentalized into layers of muscle sheets with different fiber angles. The masseter muscle originates along the length of the zygomatic process. The superficial fibers insert into the lateral lower portion of the ramus with a forward angle, while the deep fibers run mostly vertical, inserting into the lateral upper two- 158 B.2. Soft-Tissues and Muscles Anterior Digastric Posterior Digastric Upper Head of Lateral Pterygoid Lower Head of Lateral Pterygoid Stylohyoid Medial Pterygoid Mylohyoid Posterior Digastric Anterior Digastric Figure B.8: Jaw opening muscles include the upper and lower heads of the lateral pterygoid muscles and the anterior and posterior belly of the digastric muscle. c© Elsevier (2010), Drake et al. [47], adapted with permission.styloid thirds of the ramus. The fan-shaped temporalis muscle originates from a large area on the lateral side of the skull and inserts in the coronoid process and the medial side of the ramus. The temporalis fibers run between the zygomatic process and the temporal bone. The muscle is typically grouped into anterior fibers that are directed vertically and posterior fibers that are angled backward. The medial pterygoid muscle is a thick, heavily pennated muscle group originating from the pterygoid fossa and inserting into the medial lower part of the ramus. The jaw closing muscles are innervated by the mandibular branch of the trigeminal nerve [V3]. Jaw opening muscles The jaw opening muscles include the lateral ptery- goid and digastric muscles and are pictured in Figure B.8. The lateral ptery- goid muscle is divided into a superior and inferior heads. The inferior head originates from the outer surface of the lateral pterygoid plate while the upper head originates from the infratemporal surface of the sphenoid bone. Both heads insert onto the anterior neck of the condyle and the capsule of the TMJ (see Figure B.15. The lateral pterygoid is activated during jaw 159 B.2. Soft-Tissues and Muscles opening to cause forward protrusion. The digastric muscle is the primary jaw opener causing the jaw to hinge downward when contracted. The an- terior belly of the digastric originates from the digastric fossa on the lower border of the mandible. The posterior belly originates from the mastoid notch on the temporal bone. These muscle fibers connect and form an in- termediate tendon that is connected to the hyoid bone through a fibrous loop, creating pulley-like mechanism. The fibrous loop is also the insertion site of the stylohyoid muscle, which originates from the styloid process of the temporal bone. The digastric and stylohyoid muscles are also shown in Figure B.14. The hyoid bone is also connected to the mandible by the mylohyoid and geniohyoid muscles. The mylohyoid, lateral pterygoid, and anterior belly of the digastric muscles are innervated by the mandibular branch of the trigeminal nerve [V3]. The posterior belly of the digastric and stylohyoid muscles are innervated by the facial nerve [VII]. Jaw muscle architecture The micro-architecture of the jaw muscles has been described by Van Eijden et al. [205] in detailed dissection studies. The study found the jaw closing muscles to have a number of common archi- tectural characteristics, including short muscle fibers, large percentage of tendon tissue, large pennation angles, large cross sectional sizes, and rela- tively large mass, which are all indicative of physiological design for large force generation. In comparison to the closers, the jaw opening muscles were found to have smaller cross sectional sizes, smaller percentage tendon tis- sue, smaller pennation angles, and longer fiber lengths. These physiological features suggest the openers are designed for larger excursion and higher shortening velocities. The macro-architecture of the jaw muscles has also been described by Hannam and McMillan [68]. The study found the jaw closing muscles to have multiple muscle sheets of pennate fibers oriented at different angles, which was suggested as a mechanism for maintaining bite force throughout a range of jaw closing rotation. 160 B.2. Soft-Tissues and Muscles Mandible Hyoglossus Geniohyoid Hyoid bone Palatoglossus Styloglossus Genioglossus Superior Longitudinal Figure B.9: Lateral and sagittal cross-section views of the tongue muscles. c© Elsevier (2010), Drake et al. [47], adapted with permission. Genioglossus Hyoglossus Superior Longitudinal Mylohyoid Styloglossus Inferior Longitudinal Vertical & Transverse Inferior Longitudinal SeptumHyoid bone Figure B.10: Posterior views of tongue muscles with horizontal and vertical cut-away. c© Elsevier (2010), Drake et al. [47], adapted with permission. 161 B.2. Soft-Tissues and Muscles B.2.3 Tongue The tongue is a large deformable organ and is the main articulator for chang- ing the shape of the oral cavity. It is composed of a number of intrinsic muscles and is connected to surrounding orofacial structures through a set of extrinsic muscles. Tongue model descriptions have been reported from cadaver studies [192] and recently from MRI analysis [190]. The tongue plays a role in mastication by breaking up food as well as forming and po- sitioning the food bolus between chewing strokes. During swallowing, the tongue presses upward against the palate and contracts in a wave-fashion to push the bolus backward to the oropharynx, as illustrated in Figure B.13. Extrinsic tongue muscles The genioglossus muscle accounts for the bulk of the tongue tissue (see Figure B.9). Its fibers originate from the mental spine on the midline medial surface of the mandible, radiate widely in the tongue body, and insert into the mid-sagittal septum of the tongue from tip to base. The genioglossus causes tongue protrusion as the fibers pull the tongue body forward and downward. The mylohyoid is a broad, thin muscle that forms the “floor of the mouth,” as shown in Figure B.10. It originates from the mylohyoid line along the medial surface of the mandible and muscle fibers run mediolater- ally and insert on a midline raphe. Posterior mylohyoid fibers insert on the anterior surface of the hyoid bone. The geniohyoid muscle originates from the inferior mental spine on the midline medial surface of the mandible and runs backward, above the my- lohyoid, inserting on the anterior portion of the hyoid bone, as shown in Figure B.10. Geniohyoid activation pulls the hyoid forward and slightly upward depending on the position of the jaw relative to the hyoid. The hyoglossus muscle originates along the length of each side of the hyoid and insert into the lateral body of the tongue. The hyoglossus fibers run predominantly vertically with the posterior fibers angled slightly for- ward from the back of the hyoid toward the tongue. The fibers interdigitate with the styloglossus muscle at the lateral extent of the tongue body. Hyo- 162 B.2. Soft-Tissues and Muscles glossus activation pulls the tongue downward and slightly backward as well as causing the tongue body to flatten vertically. Also, if the tongue is held forward in the oral cavity, hyoglossus activation will raise the hyoid within the neck. The tongue connects above to the soft-palate with the palatoglossus mus- cle and to the skull with the styloglossus muscle as shown in Figure B.9. The palatoglossus muscle is a thin muscle and primarily used to lower the soft- palate, however it does contribute to tongue retraction. The styloglossus muscle is the main tongue retractor, originating from the styloid process, running downward and forward along the lateral extent of the tongue body, and inserting into the tongue tip. The styloglossus fibers are angled antero- posteriorly and their path changes with tongue deformation [190]. Intrinsic tongue muscles The intrinsic tongue muscles originate and insert within the tongue body and are arranged into muscle groups with roughly orthogonal directions, as shown in Figure B.10. The vertical and transverse muscles include inferior-superior and mediolateral fibers respec- tively and are interdigitated throughout the tongue body. They cause verti- cal and transverse flattening of the tongue and are activated during tongue protrusion in order to prevent lateral expansion of the tongue body, allow- ing the tip to protrude forward with genioglossus activation. The inferior longitudinal muscle has anteroposteriorly angled fibers that run along the inferior side of the tongue tip, through the tongue body, and interdigitate with the styloglossus muscle. The superior longitudinal muscle also has an- teroposterior angled fibers, but is located only along the superior surface of the tongue body. The superior longitudinal muscle’s superficial location provides leverage for backward bending of tongue body and lifting of the tongue tip. Innervation The majority of tongue muscle are innervated by the hy- poglossal nerve [XII]. Exceptions include the palatoglossus muscle (vagus nerve [X], mylohyoid muscle (trigeminal nerve [V3]), and geniohyoid muscle (cervical nerve C1). 163 B.2. Soft-Tissues and Muscles Palatoglossus Palatopharyngeus Superior Constrictor Tensor veli palatini Levator veli palatini Palatine Tonsil Uvula Figure B.11: Posterior view of the soft-palate showing superior pharyngeal constrictor (green) and the soft-palate muscles (pink). c© Elsevier (2010), Drake et al. [47], adapted with permission. B.2.4 Soft-palate The soft-palate is a structure of muscle tissue forming the upper back por- tion of the throat between the oropharynx and nasopharyx. It has numerous extrinsic muscle connections to adjacent structures, as shown in Figure B.11, in order to permit its movement and deformation. It is structurally similar to the tongue, but smaller and with a less complex structure and range of de- formation. The soft-plate is connected above to the skull by the levator veli palatini muscle. The soft-palate is connected below to the tongue and phar- ynx by the palatoglossus and palatopharyngeous muscles. These muscles are visible looking into the mouth, as the palatoglossal and palatopharygeal arches in front of and behind the palatine tonsils (see Figure B.11). Stiff- ening of the soft-palate is achieved by mediolateral forces generated by the tensor veli palatini muscle, which originates above the soft-palate in the skull, runs downward, and bends at a right-angle to insert laterally into the soft-palate body. The uvula is a midline portion of the soft-palate that hangs down in the back of the throat and it is stiffened and elevated by 164 B.2. Soft-Tissues and Muscles Superior Constrictor Middle ConstrictorInferior Constrictor SC MC IC SC MC IC Hyoid bone Cricoid cartilage Thyroid cartilage Hyoglossus Genioglossus Epiglottis Figure B.12: Oblique frontal and lateral view of the pharynx and larynx. Pharyngeal constrictors are shown: superior constrictor (SC), middle con- strictor (MC) and inferior constrictor (IC). The laryngeal cartilages and epiglottis are visible, along with a cut-away of the tongue root revealing the hyoglossus and genioglossus muscles. c© Elsevier (2010), Drake et al. [47], adapted with permission. the small musculus uvulae muscle. The soft-palate muscles are innervated by the vagus nerve [X], except for the tensor veli palatini muscle with in innervated by the mandibular division of the trigeminal nerve [V3]. B.2.5 Pharynx The pharynx is the muscle tissue that forms the back of the throat. The pharyngeal constrictors are illustrated in Figure B.12 and consist of three flat, overlapping cylindrical bands of muscle tissue originating from a mid- line pharyngeal raphe and wrapping laterally around the airway to insert into anterior structures. The superior constrictor inserts into a pterygo- mandibluar raphe (contiguous with the buccinator muscle) and forms the 165 B.2. Soft-Tissues and Muscles Food Bolus Tongue Vertebrae Hard palate Soft palate Hyoid bone Lips Mandible Vocal folds Epiglottis Figure B.13: Mid-sagittal cross-section of the oral cavity during swallowing showing closed lips, upward displaced tongue, clenched jaw, and elevating hyoid. c© Elsevier (1961), Silverman [181], adapted with permission. lateral walls of the orophayrnx. The middle constrictor inserts along the hyoid bone and stylohyoid ligament. The inferior constrictor inserts along the oblique line of the thyroid cartilage. During swallow, the pharyngeal constrictors activate in sequence from top to bottom to transport the food bolus from the oropharynx to the upper esophageal sphincter and into the esophagus. The pharyngeal constrictor muscles are innervated by the vagus nerve [X]. Three longitudinal muscles originate above the phayrnx and insert into the pharyngeal wall and are active in elevation of the pharynx and larynx. They include the stylopharyngeus, palatopharyngeus, and salpingopharyn- geus muscles and originate from the styloid process, soft-palate, and pharyn- gotympanic tubes respectively. The stylopharyngeus muscle is innervated by the glossopharyngeal nerve [IX] while the palatopharyngeus and salpin- gopharyngeus muscles are innervated by the vagus nerve [X]. B.2.6 Larynx The larynx is formed by the thyroid, cricoid and arytenoid cartilages as shown in Figure B.12. The 3D shape of the laryngeal cartilages have been an- 166 B.2. Soft-Tissues and Muscles Hyoid bone Mylohyoid Anterior Digastric Thyrohyoid Sternothyroid Omohyoid Sternohyoid Thyroid cartilage Cricoid cartilage Stylohyoid Posterior Digastric Figure B.14: Frontal view of the submandibular and neck muscles. c© Else- vier (2010), Drake et al. [47], adapted with permission. 167 B.2. Soft-Tissues and Muscles alyzed by measurements on cadaver specimens [50] and MRI analysis [169]. The cricoid cartilage is the posterior border of the larynx and is located directly superior to the trachea. The arytenoid cartilages are small paired pyramidal structures that articulate along the superior posterior surface of the cricoid. Their pyramidal shape allow for the attachment of multiple muscles and the vocal folds. The thyroid cartilage forms the anterior and superior border of the larynx and is open toward the posterior. The two sides of the thyroid cartilage form an anterior prominence, known colloquially as the “Adam’s apple” in men. The vocal and vestibule folds insert into the posterior surface of the thyroid cartilage. During swallowing the arytenoid cartilages approximate with the thyroid cartilage, and along with closure of the vocal and vestibular folds, seal off the airway from the oropharynx. Extrinsic laryngeal muscles Muscles connect the larynx above to the hyoid bone (thyrohyoid muscle) and below to the sternum (sternohyoid, omohyoid, and sternothyroid muscles), as shown in Figure B.14. These muscles function to lower and stabilize the laryngeal complex. The thyro- hyoid muscle tends to approximate the hyoid bone and larynx; it will either raise the larynx or lower the hyoid bone depending on the relative activation of the muscle above the hyoid or below the larynx. The thyrohyoid muscle is innervated by the anterior ramus of cervical nerve C1, and the sternohyoid, omohyoid, and sternothyroid muscles are innervated by the anterior rami of cervical nerves C1 through C3. Intrinsic laryngeal muscles Internal to the larynx, an intricate arrange- ment of small muscles articulate the laryngeal cartilages and manipulate the vocal folds, which are attached between the arytenoid cartilages and the in- ner surface of the thyroid cartilage. The cricothyroid muscles are used to articulate the thyroid cartilage relative to the cricoid cartilage, which is mostly rotation about the cricothyroid joint with a small amount of trans- lation. The cricoarytenoid, transverse, and oblique arytenoid muscles artic- ulate the arytenoid cartilages relative to the cricoid cartilage in a complex motion with combined mediolateral translation and oblique rotation. The 168 B.2. Soft-Tissues and Muscles intrinsic laryngeal muscles are responsible for positioning and stiffening the vocal folds and controlling glottal opening and vocal fold vibration [90]. The intrinsic laryngeal muscles are innervated by the vagus nerve [X]: the cricothyroid muscle from the external branch of superior laryngeal nerve and all others from the recurrent laryngeal branch. B.2.7 Epiglottis The epiglottis is a leaf-shaped cartilage structure located posterior to the thyroid cartilage and hyoid bone and is pictured in Figure B.12. Its broad portion extends upward posterior to the base of the tongue and is easily distinguishable on lateral medical images of the neck (see Figure 1.1). The epiglottis functions in airway protection during swallowing by folding back- ward and downward over the arytenoid cartilages of the larynx. The timing and function of the epiglottis have been examined in VF studies and its in- ternal mechanics has also been investigated [202]; however the mechanisms causing downfolding are not completely understood. There are no explicit muscle attached to the upper portion of the epiglottis that have sufficient leverage to pull down the epiglottis and therefore downfolding has been sug- gested to be a passive mechanism [58] B.2.8 Temporomandibular joint The TMJ is a compound joint composed of the mandibular condyle, tem- poral bone, and a deformable articular disc allowing for combined rotation and translation of the jaw. TMJ is pictured in Figure B.15 with the jaw at rest and during jaw protrusion. The TMJ is formed by the articular fossa and eminence of the temporal bone and the condyle of the mandible bone separated by the articular disc. The disc is composed of fibrous connective tissue and is most dense anteriorly and medially, where the largest joint forces are transmitted. The jaw closing muscles keep the TMJ in compres- sion in normal jaw function and the articular disc is held in place primarily by the concave shape of the articular fossa. A fibrous membrane around the disc and the capsular ligament also help to keep the disc in place. A number 169 B.2. Soft-Tissues and Muscles Articular disc Articular eminence Articular fossa Lateral Pterygoid muscle Jaw at rest Jaw protrusion Figure B.15: Lateral cut-away view of the temporomandibular joint at rest and during jaw protrusion. c© Elsevier (2010), Drake et al. [47], adapted with permission. Sphenomandibular ligament Stylomandibular ligament Lateral ligament Capsule Figure B.16: Lateral view of the TMJ capsule and mandibular ligaments. c© Elsevier (2010), Drake et al. [47], adapted with permission. 170 B.2. Soft-Tissues and Muscles of other mandibular ligaments help to prevent distraction of the TMJ, as shown in Figure B.16. Some muscle fibers of the upper head of the lateral pterygoid muscle insert into the fibrous membrane capsule of the articular disc and pull the disc forward during jaw protrusion. The surfaces of the condyle, articular eminence, and disc are smooth and the joint cavities lining the surfaces produce synovial fluid, which reduce friction during joint motion. The disc deforms to fit the irregular bony contact surfaces in order to distribute forces evenly. Destructive forces or structural joint changes can irreversibly change disc morphology and lead to TMJ disorder. Also, large unbalanced TMJ loading can cause dislocation of the articular disc in some cases. 171 Appendix C ArtiSynth Simulation Software 2 Our physical simulation system is embedded within the ArtiSynth platform (http://www.artisynth.org), an open-source Java-based biomechanical mod- eling toolkit developed at the University of British Columbia under the direction of Dr. Sidney Fels. Originally designed for speech applications, ArtiSynth has evolved into a tool for physiological research (particularly neuromotor control) and clinical treatment planning. I joined the ArtiSynth team in January 2005 during my Master’s thesis during which time I devel- oped the reference jaw-hyoid-larynx model working with Dr. Alan Hannam under the supervision of Dr. Sidney Fels. During my PhD, I have made a technical contribution to the ArtiSynth system by developing the inverse simulation tools described in Chapter 4. The modeling contributions de- scribed in Chapter 3 and Chapter 5 were performed using the ArtiSynth toolkit and have served as a means to develop requirements for ArtiSynth’s features and to evaluate ArtiSynth’s performance. Dr. John Lloyd is the principle designer and developer of the ArtiSynth platform. ArtiSynth models are generally created in Java code, using the packages and classes of the ArtiSynth API. Graphical editing and model creation is also supported. Applications to date have focused on the jaw and oral re- gion [71], but it is broadly applicable to biomechanical modeling in general. Key system features include (1) an architecture that supports extensive in- teractivity, including graphically based model editing and simulation control 2A version of this section has been previously published in Stavness et al. [186]. Dr. Lloyd wrote the material for this section of the manuscript, which I have included in this appendix as background on the relevant simulation techniques used in this dissertation. 172 C.1. Physical Simulation Framework (Figure 1.2), and (2) a physics engine that combines FEM and multibody capabilities, with constraints and contact handling, as described below. C.1 Physical Simulation Framework ArtiSynth models consist of a hierarchy of components, which include dy- namic components such as particles, FEM nodes, or rigid bodies, force effec- tors such as point-to-point muscles (including Hill and other types), linear or nonlinear finite elements, and constraints such as joints or collision spec- ifiers. FEM capabilities include support for tetrahedral, hexahedral, and some higher-order elements, along with linear, corotated linear [127], and some hyperelastic materials including a 5-parameter Mooney-Rivlin mate- rial. We now describe the mathematical framework for the dynamic simu- lation of these models. Let q and u be the generalized positions and velocities of all the dynam- ical components, with q̇ related to u by q̇ = Ωu (Ω generally equals the identity, except for components such as rigid bodies, where it maps angular velocity onto the derivative of a unit quaternion). Let f(q,u, t) be the force produced by all the force effector components (including the finite elements), and let M be the (block-diagonal) composite mass matrix. By representing rigid body velocity and acceleration in body coordinates we can ensure that M is constant. Newton’s second law then gives Mu̇ = f(q,u, t). (C.1) In addition, bilateral and unilateral constraints give rise to locally linear constraints on u of the form G(q)u = 0, N(q)u ≥ 0. (C.2) Bilateral constraints include rigid body joints, FEM incompressibility asso- ciated with the mixed u-P formulation [88], and point-surface constraints, while unilateral constraints include contact and joint limits. Constraints 173 C.1. Physical Simulation Framework give rise to constraint forces (in the directions G(q)T and N(q)T ) which supplement the forces of (Equation C.1) in order to enforce the constraint conditions. In addition, for unilateral constraints, we have a complemen- tarity condition in which Nu > 0 implies no constraint force, and a con- straint force implies Nu = 0. Any given constraint usually involves only a few dynamic components and so G and N are generally sparse. Solving the equations of motion requires integrating (Equation C.1) together with (Equation C.2). The presence of deformable bodies generally makes this system stiff, implying the need for an implicit integrator to obtain efficient performance3. For the work described in this paper, we use a semi-implicit second-order Newmark integrator [111], with γ = 1/2 and β = 1/4 (also known as the trapezoidal rule). Letting k index values at a particular time step, and h denote the time step size, this leads to the update rules uk+1 = uk + h 2 (u̇k + u̇k+1), qk+1 = qk + h 2 (Ωkuk + Ωk+1uk+1), (C.3) subject to Gk+1uk+1 = 0, Nk+1uk+1 ≥ 0. (C.4) Since G and N tend to vary slowly between time steps we can approximate (Equation C.4) using Gkuk+1 = gk, Nkuk+1 ≥ nk, (C.5) where gk ≡ −hĠkuk and nk ≡ −hṄkuk. Likewise, we use the approxima- tion Ωk+1 ≈ Ωk + hΩ̇k. For u̇k+1, recalling that M is constant, an estimate of the (unconstrained) value of u̇k+1 can be obtained from u̇k+1 ≈M−1fk+1, with fk+1 approximated by the first-order Taylor series fk+1 ≈ fk + ∂f k ∂u ∆u + ∂fk ∂q ∆q. 3With very soft tissue, it may sometimes be possible to use explicit methods [119], particularly if stiffness-proportional damping is excluded. 174 C.1. Physical Simulation Framework Placing this into the expression for uk+1 in (Equation C.3), multiplying by M, noting that ∆q = h/2(Ωkuk + Ωk+1uk+1) and ∆u = uk+1 − uk, and incorporating the constraints (Equation C.5), we obtain the systemM̂ k −GkT −NkT Gk 0 0 Nk 0 0  u k+1 λ z + −Mu k − hf̂k −gk −nk  =  00 w  , 0 ≤ z ⊥ w ≥ 0. (C.6) where w is a slack variable, λ and z give the average constraint impulses over the time step, and M̂k ≡M− h 2 ∂fk ∂u − h 2 4 ∂fk ∂q Ωk+1 and f̂k ≡ fk − 1 2 ∂fk ∂u uk + h 4 ∂fk ∂q Ωkuk. The complementarity condition for unilateral constraints is enforced by 0 ≤ z ⊥ w ≥ 0. A more detailed explanation of this formulation can be found in [157]. System (Equation C.6) is a mixed linear complementarity problem, a single solve of which is required to determine uk+1 for each semi-implicit integration step. Other types of integrators give rise to similar systems. A fully implicit integrator (not currently implemented in ArtiSynth) would require (Equation C.6) to be applied iteratively at each time step. For finite element models, the localized stiffness and damping matrices are embedded within ∂fk/∂q and ∂fk/∂u, which means that for models dominated by FEM components M̂ will have an FEM sparsity structure. ArtiSynth FEM models also use a lumped mass model, which ensures that M is block diagonal and makes it easier to interconnect FEM models with mass-spring and rigid body components. C.1.1 Friction, damping, and stabilization Coulomb (dry) friction can be included by extending (Equation C.6) to include either a linearized friction cone [11, 157] or a (more approximate but easier to solve) box friction [103]. ArtiSynth currently implements box 175 C.1. Physical Simulation Framework friction, and since the friction in our system tends to be quite small, we apply this as a post-hoc correction to uk+1 (in the manner of [177]), using a simplified version of (Equation C.6), with M instead of M̂ and extra constraints added in the tangential directions at contact points. Different forms of viscous damping are available, including translational and rotary damping applied directly to particles and rigid bodies, and damp- ing terms embedded in point-to-point springs and muscle actuators. For FEM models, Rayleigh damping is available, which takes the form DF = αMF + βKF where MF is the portion of the (lumped) mass matrix associated with the FEM nodes and KF is the (instantaneous) FEM stiffness matrix. DF is then embedded within the overall system matrix ∂f/∂u. In addition to solving for velocities, it is also necessary to correct posi- tions to account for drift from the constraints, including interpenetrations arising from contact. This can be done at each time step using a modified form of (Equation C.6) which computes an impulse δq that corrects the positions while honoring the constraints:M̂ k −GkT −NkT Gk 0 0 Nk 0 0  δqλ z +  0δg δn  =  00 w  , 0 ≤ z ⊥ w ≥ 0, (C.7) where δg and δn are the constraint displacements that must be corrected. If the corrections are sufficiently small, it is often permissible to use M in place of M̂k, which improves solution efficiency since M is constant and block-diagonal. While such stabilization can sometimes be incorporated directly into (Equation C.6) [10], we prefer to perform the position correction separately as this (a) allows for the possibility of an iterative correction in the case of larger errors, and (b) explicitly separates the computed velocities from the 176 C.1. Physical Simulation Framework (a) (b) (c) (d) Figure C.1: Problems with decoupling constraints from the velocity solve: In (a), a uniform 3x6x3 FEM grid of linear material with a Poisson ratio of 0 is about to be compressed by a block. The decoupled solve causes the top contacting nodes to bunch up on the surface (b), completely squashing the top two element layers, while the lower nodes hardly move at all; the coupled solve (c) causes the correct uniform displacement for all nodes. In (d), a decoupled solve causes a tongue model attached to a jaw to exhibit large vertical errors when the jaw clenches upwards against the bite plane. impulses used to correct errors. C.1.2 System solution and complexity For notational convenience, in this section we will drop the k superscripts from M̂, G, N, g, n, and f̂ in (Equation C.6) and assume that these quan- tities are all evaluated at time step k. System (Equation C.6) is a large, sparse mixed linear complementar- ity problem [38] that is not particularly easy to solve, given the unilateral constraints and the fact that M̂ is not block diagonal. If M̂ is symmet- ric positive definite (SPD), it is equivalent to a convex quadratic program. If there are no unilateral constraints (N = ∅), then it reduces to a linear Karush-Kuhn-Tucker (KKT) system. Generally, M̂ is symmetric (unsymmetric terms sometimes arise from ro- tational effects but these are usually small enough to ignore) and hence will also be SPD for small enough h (since M is SPD). However, the resulting sys- tem is still harder to solve than non-stiff multibody systems where M̂ = M. This is because M̂, while still sparse, is not block-diagonal. Multi-body 177 C.1. Physical Simulation Framework systems are often solved using the projected Gauss-Seidel method [103]. However, this involves a sequence of iterations, each requiring the computa- tion of GiM̂−1GTi or NiM̂ −1NTi , which is easy to do for a block-diagonal M but much more costly for M̂. It is tempting to follow the approach we use for friction and decouple the velocity and constraint solves, by first computing u∗ = M̂−1(Muk +hf̂) and then applying constraints to u∗ in a post-hoc fashion, using a version of (Equation C.6) in which M̂ is replaced with M. This can be done by various methods, including Gauss-Seidel iteration, and is equivalent to projecting u∗ onto the space of legal velocities. Unfortunately, this does not propagate constraint effects properly throughout the system, and can result in very large errors when the constraint forces are large, as illustrated in Figure Figure C.1. At present, ArtiSynth solves (Equation C.6) by using a Schur comple- ment to turn it into a dense regular linear complementarity problem N̄A−1N̄T z + N̄A−1b− n = w 0 ≤ z ⊥ w ≥ 0 (C.8) where A ≡ ( M̂ −GT G 0 ) , N̄ ≡ ( N 0 ) , b ≡ ( Muk + hf̂ g ) which is solved using Keller’s algorithm [103]. uk+1 and λ can then be obtained using back-substitution:( uk+1 λ ) = A−1 ( b + N̄T z ) . (C.9) Keller’s algorithm is a pivoting method with an expected complexity of O(m3), wherem is the number of unilateral constraints. In addition, forming (Equation C.8) and the back solve of (Equation C.9) requires m + 1 solves of a system involving A. This is done using the Pardiso sparse direct solver 178 C.1. Physical Simulation Framework 102 103 104 100 102 104 number of nodes so lv e tim e (m s) Figure C.2: Log-log plot showing factor times for A as a function of the number of nodes (which is proportional to the size of A) for a series of 3D FEM problems with a uniformly increasing node density. The slope of the line indicates a complexity of O(n1.7). [163], and entails a once-per-step factoring of A, plus m+1 solve operations. Experimentally, we have determined that the complexity of factoring A (using Pardiso) for 3D FEM type problems is roughly O(n1.7), where n is the size of A (Figure Figure C.2). Similarly, we have also determined that the complexity of solving a factored A is roughly O(n1.3). Hence we can expect the overall complexity for solving (Equation C.6) to be O(m3) +mO(n1.3) +O(n1.7). This works well provided that the number of unilateral constraints m is small. To help achieve this, we can sometimes treat the unilateral constraints arising from contact as bilateral constraints (i.e., entries in G) on a per-step basis, as described further in Section C.1.4. C.1.3 Attachments between bodies In creating comprehensive anatomical models, it is often necessary to attach various bodies together. Most typically, this is done by connecting points of one body to specific locations on another body. For example, FEM nodes 179 C.1. Physical Simulation Framework may be attached to particular spots on a rigid body, or to other nodes of a different FEM model. To facilitate this, ArtiSynth provides the ability to attach a dynamic component to one or more master components. Let the set of attached components be denoted by β, and the remaining set of unattached active components be denoted by α. In general, the velocity uj of an attached component is related to the velocities uα of the active components by a locally linear velocity constraint of the form uj + Gjαuα = 0. Gjα will be sparse except for entries corresponding to the master compo- nents to which j is attached. Letting Gβα denote the composite matrix formed from Gjα for all attached components, we have Iuβ + Gβαuα = 0 for the constraints that enforce all attachments. We could simply add these constraints to (Equation C.6) and solve the resulting system, but this would increase both the system size and solution time. Instead, we use the attachments to actually reduce the size of (Equa- tion C.6). Consider first the subsystem involving only bilateral constraints. As in Section C.1.2, we drop the k superscripts from M̂, G, N, g, n, and f̂ in (Equation C.6) and assume that these quantities are all evaluated at time step k. Letting b ≡Muk + hf̂ and partitioning the system into active and attached components yields M̂αα M̂αβ GTαα G T βα M̂βα M̂ββ GTαβ I Gαα Gαβ 0 0 Gβα I 0 0   uk+1α uk+1β λα λβ  =  bα bβ gα 0  . 180 C.1. Physical Simulation Framework The identity submatrices make it easy to solve for uk+1β and λβ: uk+1β = −Gβαuk+1α , λβ = bβ − M̂βαuk+1α + M̂ββGβαuk+1α −GTαβλα and hence reduce the system to( M̂′ G′T G′ 0 )( uk+1α λα ) = ( b′ gα ) (C.10) where M̂′ ≡ PM̂PT , G′ ≡ GPT , b′ ≡ Pb, with P ≡ ( I −GTβα ) . Similarly, unilateral constraints can be reduced via N′ = NPT . The reduction operation can be performed in O(n) time and results in a system that is less sparse but generally faster to solve than the original. Most attachments in ArtiSynth are point-based, with the most common kind being the attachment of an FEM node to a rigid body. It is also possible to attach FEM nodes to the faces and edges of an FEM element, allowing us to handle the so-called “T-junction” problem and create FEM models with non-conforming element faces. This is quite useful for creating localized subdivisions of particular elements, particularly hexahedrons. C.1.4 Contact handling Collision detection can be enabled between any combination of rigid or de- formable bodies. It is assumed that the bodies in question contain a triangu- lar surface mesh. A bounding-box hierarchy is used to determine if any two surfaces meshes intersect. If they do, then a tracing algorithm (similar to [3]) is used to compute all the intersection contours between the two meshes as shown in Figure Figure C.3. Such contour tracing can be done relatively quickly but does require the use of robust geometry predicates similar to those in [51]; this is particularly true because collision conditions tend to drive the contacting surfaces into degenerate mesh configurations. Determining the intersection contour allows us to easily create a set of 181 C.1. Physical Simulation Framework Figure C.3: Contact handling between two deformable models (with the topmost rendered as a wireframe), showing the intersection contour (blue) and the contact normals (black lines) of interpenetrating vertices from the upper mesh. constraints for correcting the interpenetration and preventing interpenetrat- ing velocities. It also provides a good estimate of the contact area, which can be used for determining contact pressure. For collisions involving a deformable body, we locate all mesh vertices which are interior to the contour. Each such vertex corresponds to surface FEM node which is interpenetrating the other body. For each interpene- trating node, we then find the nearest point and face on the opposite mesh, and use the face’s normal n as a contact normal. A linear one-dimensional constraint is then created which prohibits relative motion in the negative normal direction between the node and nearest point. If the opposite face is located on a deformable body, this results in a constraint between the ve- locity of the node vn and the velocities v0, v1, v2 of three nodes associated with the nearest face: nTvn − w0nTv0 − w1nTv1 − w2nTv2 ≥ 0 where w0, w1, and w2 are the barycentric coordinates of the nearest point with respect to the face. If the opposite face is located on a rigid body, then the constraint is between vn and the body’s translational and angular 182 C.1. Physical Simulation Framework velocity, vb and ωb (expressed in body coordinates): nTvn − nTRvb + nTR(p× ωb) ≥ 0 where R is the rotational transformation from body to world coordinates and p is the location of the nearest point in body coordinates. Each of these constraints can be expressed in the general form Niu ≥ 0, where Ni is a row of N and u is the vector of all velocities. These constraints serve both to prevent interpenetrating velocities in the contact direction, and to remove interpenetrations during the position correction step (Equation C.7), with the correction distance taken to be the distance d between the interpenetrating node and its nearest point on the other body. Intersection contours are also used to determine contact constraints for rigid-body/rigid-body contact, although we omit the details here for brevity. As mentioned in Section C.1.2, the solution time of (Equation C.6) can be greatly improved if some contact constraints can be temporarily treated as bilateral constraints within a particular time step. By default, ArtiSynth does this for contact involving deformable bodies, since such bodies have many degrees of freedom and their contact constraints tend to be somewhat decoupled. To prevent sticking, each contact’s vertex-face pair is stored between time steps, and if it reappears in the next step, it is used as a contact constraint only if its corresponding λ value computed in (Equation C.6) is ≥ 0, implying that there is no force trying to make it separate. This is in effect an active set method, with the active set used to solve (Equation C.6) being updated between steps. It should be noted that we do not claim that the collision handling scheme described here is optimal for all applications. In particular, we do not currently implement edge-edge type contacts, and so there can be some interpenetration which depends on the coarseness of the surface meshes. However, the collision handling is properly isolated from the rest of the simulation, and other collision handling schemes can be easily used as long as they provide a set of constraints for enforcing the contact and resolving 183 C.1. Physical Simulation Framework interpenetrations. C.1.5 Simulation engine summary The complete ArtiSynth simulation engine is summarized below. It uses the concept of [67, 177] whereby velocities are computed in advance of po- sitions, subject to constraints, to help prevent constraint violations during the subsequent position computation. 1. Compute contacts (as per Section C.1.4) and the bilateral and unilat- eral constraint matrices Gk and Nk. 2. Correct positions qk to remove interpenetration and drift errors, using (Equation C.7). 3. If necessary, adjust Gk and Nk to reflect changes in q. 4. Solve for uk+1 using (Equation C.6). 5. Adjust velocities uk+1 for dry friction, as described in Section C.1.1. 6. Compute new positions: qk+1 = qk + h/2(Ωk+1uk+1 + Ωkuk). This algorithm is generally applicable to any rigid-deformable body dy- namics. In the absence of constraints, the above system turns into a trape- zoidal rule solution of a regular ordinary differential equation, for which global errors can be expected to be proportional to O(h2). The inclusion of constraints, particularly non-smooth unilateral ones, makes formal conver- gence and error analysis more difficult. However, the main velocity update (Equation C.6) is the same as that described in [157], which is shown to be stable and have second-order convergence under certain assumptions. Generally speaking, our method is a time-stepping scheme which uses fixed (or adaptively varying) time steps, as opposed to an event-driven scheme in which the integration time intervals are precisely aligned with contact events but which becomes impractical in the presence of large num- bers of contacts. To the extent to which results exist, time-stepping schemes 184 C.1. Physical Simulation Framework (a) (b) (c) Figure C.4: Examples used for validation, shown in their final positions with a rainbow plot of the resulting Von Mises stresses and the locations of their respective reference nodes. (a) Beam example, ArtiSynth. (b) Cube example, ArtiSynth (maximum stress 2787 Pa). (c) Cube example, ANSYS (maximum stress 2661 Pa). are typically shown to have less accuracy but better convergence properties that event-driven ones [25]. C.1.6 Validation using ANSYS To help assess the performance of our integration scheme, we compared it against the commercial finite element package ANSYS for two test examples: a beam, fixed at one end and allowed to fall under gravity, and a cube, resting on a flat surface and subjected to a downward load applied to several top nodes. It should be noted that ArtiSynth uses several simplifications compared to ANSYS, notably the use of semi-implicit integration and a lumped mass model. The beam example (Figure C.4 (a)) consisted of a beam with dimensions 0.1 x 0.02 x 0.02 m divided uniformly into 8 x 4 x 4 hexahedral elements, with a density of 1040 kg/m3, Rayleigh damping coefficients of α = 20 s−1 and β = 0.015 s, and a five-parameter Mooney Rivlin material with c10 = 10370, c20 = 486 and c01 = c11 = c02 = 0 (Pascals). Incompressibility in both system was enforced using a mixed u-P formulation [88], and time integration was performed for .4 seconds using a one millisecond time step. To assess dynamic performance, the resulting z displacement and velocity of a reference node located in the middle of the free end was compared over time between ArtiSynth and ANSYS (Figure C.5 (left)). The dynamic 185 C.1. Physical Simulation Framework behavior was essentially identical: the resulting displacement error (relative to the maximum displacement) had maximum and average values of 0.3% and 0.08%. Likewise, the resulting velocity error (relative to the maximum velocity) had maximum and average values of 1.3% and 0.4%. We also determined the errors in total displacement and Von Mises stress for all the nodes in the final position: the maximum and average displacement errors (relative to the maximum displacement) were 0.06% and 0.04%, while the maximum and average Von Mises stress errors (relative to the maximum stress value) were 0.9% and 0.13%. The cube example (Figure C.4 (b)) used a cube with a width of 0.1 m in all directions and divided uniformly into 6 x 6 x 6 hexahedral elements, with a density of 1040 kg/m3, Rayleigh damping coefficients of α = 20s−1 and β = 0.015s, and a five-parameter Mooney Rivlin material with c10 = 1037, c20 = 486 and c01 = c11 = c02 = 0 Pascals (identical to the material used for our tongue model described below). In addition to gravity, an immediate external load of -0.8 Newtons was applied in the vertical direction to the nine nodes located in the middle of the top surface, resulting in the deformation shown in Figure C.4. Incompressibility in both system was enforced using the B-bar method [88], and the example was integrated for 0.2 seconds with a one millisecond time step. To assess dynamic performance, a reference node was selected in the middle of the top surface and its z displacements and velocities were compared (Figure C.5 (right)). Displacement errors had maximum and average values of 2.7% and 1.5%, while the velocity errors had maximum and average values of 22% and 1.4%. The large value for the velocity error occurred at the beginning, where ANSYS computed an unexpected initial upward velocity for the node. Compared with ANSYS, the ArtiSynth behavior was slightly more damped. For all nodes in the final position, the maximum and average displacement errors were 3.5% and 0.5%, while the maximum and average Von Mises stress errors were 5.6% and 0.5%. Much of this error was due to differences in the way ArtiSynth and ANSYS compute pressure for the B-bar method, resulting in different dilational displacements: in ArtiSynth the model compressed slightly, while in ANSYS it inflated slightly. 186 C.2. Graphical Toolset 0 100 200 300 400−0.05 −0.04 −0.03 −0.02 −0.01 0 z di sp la ce m en t 0 40 80 120 160 200 −0.02 −0.01 0 z di sp la ce m en t 0 100 200 300 400−0.4 −0.3 −0.2 −0.1 −0 0.1 time (ms) z ve lo ci ty 0 40 80 120 160 200 −0.6 −0.4 −0.2 0 time (ms) z ve lo ci ty Beam example Cube example Figure C.5: Time integration comparisons between ANSYS (solid lines) and ArtiSynth (dotted lines), showing the z displacement (top) and velocity (bottom) of a single reference point in the beam and cube examples. These results help demonstrate that our simulation approach is com- petitive with commercially available codes, in addition to be considerably more efficient: ArtiSynth was 20 and 10 times faster for the beam and cube examples, respectively. C.2 Graphical Toolset ArtiSynth is implemented in Java, both for portability and to take advan- tage of Java’s extensive graphical user interface (GUI) support. Models are generally created in Java code, using the packages and classes of the Ar- tiSynth API. Interactive editing is also possible, as described below. Models can be saved and loaded using a text file format which can also be used for model creation. 187 C.2. Graphical Toolset Figure C.6: ArtiSynth application showing two model views (including a medical image plane), and the timeline for arranging probes. C.2.1 Model component hierarchy An ArtiSynth model is arranged as a hierarchy of components, which include dynamic components such as particles or rigid bodies, force effectors such as point-to-point muscles (including Hill and other types), linear or nonlinear finite elements, and constraints such as joints or collision specifiers. Other models can be included in the hierarchy as submodels. At the top of the hierarchy is the RootModel, which can contain special components such as probes and control panels used for interactive simulation control. C.2.2 Viewing, selection, and editing ArtiSynth provides a graphical interface for model viewing, editing, and simulation control (Figure C.6). Rendering is done using OpenGL, and a Jython console permits scripting and interactive command interaction. Multiple viewers can be created, with aids such as grids, clipping and slicing planes, and image planes overlaid with medical imaging data (Figure C.6, 188 C.2. Graphical Toolset Figure C.7: Navigation panel (left) and muscle control panel (right) for a rat hind limb model (bone meshes courtesy of Dr. Dinesh Pai). right). The component hierarchy can be viewed and components can be selected using a navigation panel (Figure C.7). Components can also be selected in the viewers. Models can be edited using a context-based system in which the current set of selected components determines the available editing functions. Vari- ous components can be added, duplicated, and deleted. Portions of a model can also be written to a file for later use. Dragger fixtures provide transla- tion, rotation and scaling. While it is often easier to create complex models in code, interactive editing provides a way to adjust and modify existing models, and perform certain kinds of specialized tasks, such as inserting muscle fibers into an FEM. C.2.3 Properties, control panels, and probes Model components can export various properties that describe attributes such as mass, damping, force or position. These are implemented using the Java reflection package and can thus be declared easily using one or two lines of code. A collection of properties can be gathered into a composite property; examples of this include render properties that control rendering attributes such as color, shading, and visibility, and material properties, which control 189 C.2. Graphical Toolset the parameters of FEM constitutive laws. Properties can be declared as inheritable, so that their values can be either explicitly set or inherited from an identical property located in a component further up the hierarchy. This provides an inheritance function- ality similar to that found in scene graphs. For example, you can set a default material stiffness value within an FEM model component, and then override this value as necessary within individual FEM elements. Properties provide the main connection handles for various interactive simulation tools. These include control panels, which contain widgets for ad- justing property values. A widget can be automatically created and added to a control panel by simply selecting a specific property of a particular component (or group of components). In a typical application, a user might create a panel containing slider widgets connected to the excitation proper- ties for a set of muscles (Figure C.7). These can then be adjusted during simulation to see the effect of activating specific muscles. A property whose value is a numeric scalar or vector can also be at- tached to a channel of input or output data (known as a probe), which may be used to control (or record) the property’s value over time. Probes can be temporally arranged in a graphical timeline (Figure C.6) which provides control over their start time and duration. In neuromotor control work, input probes are commonly used to drive muscle excitation levels, while output probes are used to record the resulting forces, velocities, and posi- tions. The numeric data within a probe can be displayed graphically, and for input probes can be graphically edited by adding or deleting knot points, changing interpolation, etc. 190

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 31 0
Canada 31 1
Germany 21 5
Unknown 19 0
Taiwan 8 0
Netherlands 8 8
Russia 5 0
Romania 5 0
Pakistan 3 0
United Kingdom 2 0
Finland 2 0
Poland 2 0
Norway 1 0
City Views Downloads
Unknown 39 0
Braunschweig 16 5
Vancouver 14 0
Saskatoon 13 1
Mountain View 11 0
Amsterdam 8 8
Taipei 7 0
Washington 5 0
Ashburn 5 0
Oulu 2 0
Regina 2 0
Karachi 2 0
Herborn 2 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items