A Nonlinear Finite Element Model of the Rat Cervical Spine Validation and Correlation with Histological Measures of Spinal Cord Injury by Colin Macdonald Russell B.A.Sc., The University of British Columbia, 2008 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Biomedical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2012 ➞ Colin Macdonald Russell 2012 Abstract Researchers and clinicians do not currently use the heterogeneity of the primary mechanism of spinal cord injury (SCI) to tailor treatment strategies because the effects of these distinct patterns of acute mechanical damage on long-term neuropathology have not been fully investigated. Computational modelling of SCI enables the analysis of mechanical forces and deformations within the spinal cord tissue that are not visible experimentally. I created a dynamic, hyperviscoelastic three-dimensional finite element (FE) model of the rat cervical spine and simulated contusion and dislocation SCI mechanisms. I investigated the relationship between maximum principal strain and previously published tissue damage patterns, and compared primary injury patterns between mechanisms. My model incorporates the spinal cord white and gray matter, dura mater, cerebrospinal fluid, spinal ligaments, intervertebral discs, a rigid indenter and vertebrae, and failure criteria for ligaments and vertebral endplates. High-speed (1 m/s) contusion and dislocation injuries were simulated between vertebral levels C3 and C6 to match previous animal experiments, and average peak maximum principal strains were calculated for several regions at the injury epicentre and at 1 mm intervals from +5 mm rostral to -5 mm caudal to the lesion. I compared average peak principal strains to tissue damage measured previously via axonal permeability to 10 kD fluorescein-dextran (Choo, 2007). Linear regression of tissue damage against peak maximum principal strain for pooled data within white matter regions yields significant (p < 0.0001) correlations that are similar for both contusion (R2 = 0.86) and dislocation (R2 = 0.54). With additional simulations of cord contusion injuries at lower injury velocities of 3 and 300 mm/s, I found that current material properties used to model the cord are not biofidelic within this velocity range. By fitting existing experimental cord material testing data and plotting alongside the material properties used in several related models, I further demonstrated the remaining divide between experimental data and computational models. My model enhances our understanding of the differences in injury patterns between SCI mechanisms, and provides further evidence for the link between principal strain and tissue damage. Furthermore, my results speak to a continued need to test cord material properties at a range of strains and strain rates to better refine cord hyperviscoelastic properties. ii Preface Dr. Tae-Eun Chung assisted with the choice of finite element modelling software for the project, and created the initial finite element meshes for the model. I conducted the remaining majority of the model development and refinement, and I designed and conducted all of the simulations and analysis presented in this thesis, under the guidance of my supervisor Dr. Thomas Oxland. Portions of this thesis detailing model development and the simulation and results for the injury mechanism experiments have been published in an article titled Maximum Principal Strain Correlates with Spinal Cord Tissue Damage in Contusion and Dislocation Injuries in the Rat Cervical Spine on May 8th , 2012 in the Journal of Neurotrauma [95], and are included here with permission. Dr. Anthony Choo, Dr. Wolfram Tetzlaff, Dr. Tae-Eun Chung, and Dr. Thomas Oxland were co-authors and contributed edits and several paragraphs to the manuscript, of which I wrote the majority. Sections containing content detailed in the article include 1.1, 1.7, 2.1, 2.3.2, 2.4.2, 2.5, 3.2.2, 3.3, 3.4, 3.5, 4.1, 4.3, 5.4. iii Table of Contents Abstract ii Preface iii Table of Contents iv List of Tables vii List of Figures viii List of Abbreviations xi Acknowledgements xii Dedication 1 Introduction 1.1 Background and motivation . . . . . . . . . . . . . . . . . . . 1.1.1 Human FE model . . . . . . . . . . . . . . . . . . . . 1.1.2 Injury mechanism experiments . . . . . . . . . . . . . 1.1.3 Injury velocity experiments . . . . . . . . . . . . . . . 1.2 Project definition . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Anatomy of the rat cervical spine . . . . . . . . . . . . . . . . 1.3.1 Anatomy of the spinal cord . . . . . . . . . . . . . . . 1.4 Experimental mechanisms of spinal cord injury . . . . . . . . 1.5 Modelling of cord injury . . . . . . . . . . . . . . . . . . . . . 1.5.1 Strain theory . . . . . . . . . . . . . . . . . . . . . . . 1.5.2 Material properties of the spinal cord . . . . . . . . . 1.5.3 Material modelling of the spinal cord . . . . . . . . . . 1.5.4 Finite element modelling of the spine, spinal cord, and 1.6 Mechanical indicators of tissue injury . . . . . . . . . . . . . . 1.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . brain . . . . . . . . 2 Methods 2.1 Model development . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Geometry extraction from Magnetic Resonance Imaging 2.1.2 Finite element meshing . . . . . . . . . . . . . . . . . . 2.1.3 Material properties . . . . . . . . . . . . . . . . . . . . . 2.1.4 Fluid-Structure Interaction and the cerebrospinal fluid . 2.2 Material model investigation . . . . . . . . . . . . . . . . . . . 2.2.1 Nonlinear regression of rat cord tensile data . . . . . . . 2.2.2 Tensile coupon simulations . . . . . . . . . . . . . . . . 2.3 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 1 2 3 4 4 4 6 8 9 9 11 12 15 16 17 . . . . . . . . . 19 19 19 21 23 25 26 26 27 28 iv Table of Contents 2.4 2.5 2.3.1 Initial weight-drop validation . . . 2.3.2 Velocity and mechanism validation Injury simulation . . . . . . . . . . . . . . 2.4.1 Injury velocity experiments . . . . 2.4.2 Injury mechanism experiments . . Correlation with histology . . . . . . . . . 2.5.1 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Results 3.1 Material model investigation . . . . . . . . . . . . 3.1.1 Nonlinear regression of rat cord tensile data 3.1.2 Tensile coupon simulations . . . . . . . . . 3.2 Validation . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Initial weight-drop validation . . . . . . . . 3.2.2 Velocity and mechanism simulations . . . . 3.3 Internal strain distributions . . . . . . . . . . . . . 3.4 Correlation with histology . . . . . . . . . . . . . . 3.5 Regional distribution of strain and tissue damage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 28 28 28 29 30 31 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 32 32 35 37 37 39 41 41 43 4 Discussion 4.1 Injury simulation . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Interpretation of cord strain patterns and correlations 4.1.2 Force history delay and subarachnoid space . . . . . . 4.2 Material model investigation . . . . . . . . . . . . . . . . . . 4.3 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Spinal cord material properties . . . . . . . . . . . . . 4.3.2 Strain direction . . . . . . . . . . . . . . . . . . . . . . 4.3.3 Cerebrospinal fluid validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 45 45 47 48 49 49 49 50 5 Conclusion 5.1 Conclusions . . . . . . . . . 5.2 Contributions . . . . . . . . 5.3 Recommendations for future 5.4 Concluding statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 51 52 53 54 . . . . . . . . work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography 55 Appendices A Images of Rat Vertebrae from Literature 65 B Procedure for extracting simulation results for solid elements in each predefined spinal cord zone slice 69 C MATLAB code C.1 Hyperviscoelastic curve fitting - hvstress.m . . . . . . . . . . . . . . C.2 Hyperviscoelastic optimized parameter display - hvlabels.m . . . . . C.3 Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m . . . . C.4 Hyperviscoelastic algorithm validation scripts . . . . . . . . . . . . . C.4.1 Recreation of Greaves hyperelastic plot - greavesogdentest.m C.4.2 Recreation of Miller hyperviscoelastic plots - millertesthv.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 70 71 72 80 80 80 v Table of Contents C.4.3 Recreation of Snedeker hyperviscoelastic plots - snedeckertest.m . . . . . . . 82 C.5 Material model comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 C.5.1 Model parameters - output from hvlabels.m . . . . . . . . . . . . . . . . . . . 89 C.6 FE data extraction and linear regression of tissue damage versus maximum principal strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 C.6.1 Reading of exported FE data from .csv files and saving to .mat - zonesliceread.m 89 C.6.2 Calculation of mean peak values and SD and saving to .mat - mechanismFEdatasave.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 C.6.3 Plotting of regional FE stress and strain as function of slice position - zonesliceplot.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 C.6.4 Pooled and regional correlation plots of FE stress and strain with tissue damage (dislocation epicentre points omitted) - zoneslicecorrplotOMIT.m . . 97 C.6.5 Regional time history plots of FE stress and strain - zonecurveplot.m . . . . 103 D Histological and FE parameters as function of distance from injury epicentre 106 D.1 Histological data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 D.2 FE results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 E Correlation plots 114 E.1 Correlations for pooled white matter regions . . . . . . . . . . . . . . . . . . . . . . . 115 E.2 Regional correlations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 F Embedded 3D Model of Segmented Geometry of the Rat Cervical Spine 121 vi List of Tables Table Table Table Table 2.1 2.2 2.3 2.4 Ratio of rat to human cord cross sectional areas . . . . . . . . . . . . . . . . Cross sectional areas of spinal ligaments . . . . . . . . . . . . . . . . . . . . . Material properties of spinal cord and dura . . . . . . . . . . . . . . . . . . . Material properties of spinal ligaments (Material Type 205 - Nonlinear Tension Only Bar Element) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 2.5 Material properties of intervertebral disc and endplate spotwelds (Material Type 1 - Elastic-Plastic for Solid Elements) . . . . . . . . . . . . . . . . . . . . 23 . 23 . 24 . 25 . 25 Table 3.1 Correlation coefficients for maximum principal strain and tissue damage within cord regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 vii List of Figures Figure 1.1 Mid-sagittal view of von Mises strains for transverse contusion, distraction, and dislocation injuries [Figure and caption text reproduced with permission from Greaves et al. [40].] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Figure 1.2 Illustrations of experimental injury configurations [Figure and caption text reproduced with permission from Choo et al. [16].] . . . . . . . . . . . . . . . 2 Figure 1.3 Parasagittal sections of the spinal cord demonstrating the hemorrhage resulting from contusion for the control, slow, and fast groups. [Figure and caption text reproduced with permission from Sparrey et al. [107].] . . . . . . . . . . . 3 Figure 1.4 Load-displacement curves for the slow and fast contusion groups [Figure and caption text reproduced with permission from Sparrey et al. [107].] . . . . . . 3 Figure 1.5 Spatial terminology with respect to the rat [Illustration reproduced from Wingerd [120].] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Figure 1.6 Diagram of a rat skeleton [Illustration reproduced from Muskopf [78] under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 United States License.] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Figure 1.7 Diagrams of generic rat vertebrae [Illustrations reproduced from Wingerd [120].] 5 Figure 1.8 Diagrammatic transverse section of the medulla spinalis [human spinal cord] and its membranes. [Illustration and caption text reproduced from Gray [39], now in the public domain.] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Figure 1.9 Diagram of a typical neuron. [Reproduced under the GNU Free Documentation and the Creative Commons Attribution-Share Alike 3.0 Unported Licenses (http://commons.wikimedia.org/wiki/File:Neuron.svg).] . . . . . 6 Figure 1.10 Cross section and diagram of the rat spinal cord at C4. [Images and structure abbreviation list reproduced from Watson et al. [117].] . . . . . . . . . . . . . 7 Figure 1.11 2D strain geometry for an infinitesimal material element [Diagram taken from the public domain.] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Figure 1.12 3D framework for finite strain theory [Diagram taken from the public domain.] 12 Figure 1.13 Hyperelastic and viscoelastic properties of the spinal cord [Figures reproduced with permission from a review by Clarke [17].] . . . . . . . . . . . . . . . . . . 13 Figure 1.14 Prony series model schematic [Diagram taken from the public domain.] . . . . 14 Figure 2.1 7T animal MRI scanner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.2 Active contour evolution [Figure and caption text reproduced with permission from [124].] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.3 Extraction of rat cervical spine geometry . . . . . . . . . . . . . . . . . . . . Figure 2.4 Surfaces segmented in ITK-SNAP . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.5 Dura mater and spinal ligaments . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.6 Tensile coupon simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.7 Four-vertebra meshed model . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.8 Contusion and dislocation displacement cross-sections . . . . . . . . . . . . Figure 2.9 Regions of the spinal cord used for quantification of strain and histology. [Illustration adapted with permission from Choo et al. [15].] . . . . . . . . . . 19 . . . . . . . 20 21 22 22 27 29 29 . 30 viii List of Figures Figure 2.10 Sample strain results from spinal cord regions used for quantification . . . . . 31 Figure 2.11 Examples of finite element strain time histories . . . . . . . . . . . . . . . . . 31 Figure Figure Figure Figure Figure 3.1 3.2 3.3 3.4 3.5 Figure 3.6 Figure 3.7 Figure 3.8 Figure 3.9 Figure 3.10 Figure 3.11 Figure 3.12 Figure 3.13 Figure 3.14 Figure 3.15 Recreation of hyperelastic plot from Greaves [41] . . . . . . . . . . . . . . . . Recreation of hyperviscoelastic plots from Miller and Chinzei [74] . . . . . . . Recreation of hyperviscoelastic plots from Snedeker et al. [104] . . . . . . . . Stress-strain plots of experimental data and model predictions . . . . . . . . . Stress-strain plots of experimental fit and model predictions – wide stretch ratio range . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Tensile coupon simulation comparison of PAM-CRASH hyperviscoelastic materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Effect of single versus double precision simulation arithmetic on long-term stress relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Maximum principal strain distribution comparison of Maikos et al. [68] 12.5 mm weight-drop contusion results (a) with UBC model results (b). [Figure in (a) adapted from Maikos et al. [68].] . . . . . . . . . . . . . . . . . . . . . . . . . Results of 12.5 mm weight drop simulations for varying model mesh sizes and other model variations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Force-displacement curves for 12.5 mm and 25 mm weight-drop simulations with α = ±4.7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Force displacement curves for injury velocity simulations . . . . . . . . . . . . Simulated versus experimental applied contusion and dislocation forces . . . . Distribution of maximum principal strain during contusion and dislocation injury simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Correlations of maximum principal strain with tissue damage for contusion and dislocation mechanisms in the white and gray matter. . . . . . . . . . . . Rostrocaudal distributions of experimentally measured tissue damage and computed maximum principal strain . . . . . . . . . . . . . . . . . . . . . . . 32 33 33 34 35 36 36 37 38 39 40 40 41 42 44 Figure 4.1 Pattern of maximum principal strain in cervical dislocation is different from experimental tissue damage in distributed thoracolumbar dislocation [Diagram in (b) adapted from Clarke et al. [18].] . . . . . . . . . . . . . . . . . . . 47 Figure A.1 Cervical vertebrae and nervous system [Illustrations and caption text reproduced from Greene [42].] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure A.2 Cervical and thoracic vertebrae [Illustrations reproduced from Wells [118].] . Figure A.3 Cervical vertebrae [Illustrations reproduced from Rowett [93].] . . . . . . . . Figure A.4 Cervical and thoracic vertebrae [Illustrations reproduced with permission from Johnson et al. [53].] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure Figure Figure Figure Figure Figure Figure D.1 D.2 D.3 D.4 D.5 D.6 D.7 Distance Distance Distance Distance Distance Distance Distance plots plots plots plots plots plots plots – – – – – – – Choo cellular permeability . . . . . . . . . first principal strain (Maikos properties) . . second principal strain (Maikos properties) third principal strain (Maikos properties) . first principal stress (Maikos properties) . . second principal stress (Maikos properties) third principal stress (Maikos properties) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 . 66 . 67 . 68 . . . . . . . 107 108 109 110 111 112 113 Figure E.1 Parameter correlations – pooled data for all white matter regions (Maikos properties) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Figure E.2 Parameter correlations – white matter dorsal column region (Maikos properties)116 ix List of Figures Figure Figure Figure Figure E.3 E.4 E.5 E.6 Parameter correlations – white matter lateral column region (Maikos properties)117 Parameter correlations – white matter ventro-medial region (Maikos properties)118 Parameter correlations – white matter ventro-lateral region (Maikos properties)119 Parameter correlations – gray matter region (Maikos properties) . . . . . . . . 120 x List of Abbreviations 2D - two-dimensional SCI - spinal cord injury MRI - magnetic resonance imaging FE - finite element NURBS - non-uniform rational B-splines AP - anterior-posterior GUI - graphical user interface IV - intervertebral CSF - cerebrospinal fluid ALL - anterior longitudinal ligament PLL - posterior longitudinal ligament JC - joint capsule ligament LF - ligamentum flavum ISL - interspinous ligament DL - denticulate ligaments DA - dural attachments xi Acknowledgements I am deeply grateful to my supervisor, Dr. Thomas Oxland, who gave me the opportunity to begin this project during a co-op work term and went on to encourage me to further pursue my interests in this work through graduate school. Also, many thanks to Dr. Tae-Eun Cheung, Dr. Peter Cripton, Erin Lucas, Tim Bhatnagar, Hannah Gustafson, Dr. Claire Jones, James Boak, Carolyn Van Toen, Cameron Lam, and Maryam Shahrokni for their help and discussions along the way, and to Dr. Anthony Choo and Dr. Wolfram Tetzlaff for their contributions to the Journal of Neurotrauma manuscript. I would also like to thank all of the past and present staff and students at the Orthopaedic and Injury Biomechanics Group for making my time spent there so rewarding and enjoyable. xii To my loving wife, Delphina, and to my family for all of their support along this journey. Either you decide to stay in the shallow end of the pool or you go out in the ocean. - Christopher Reeve xiii Chapter 1 Introduction In the introduction I first outline the background and significance of the project, then define the project and list of primary objectives followed by an overview of strain theory and a review of spinal anatomy and relevant literature. I conclude the chapter with a summary of my project goals and how they address limitations of, and expand on, previous work by others. 1.1 Background and motivation Traumatic spinal cord injury (SCI) often results in a debilitating condition, with estimated incidence rates of 1,800 each year in Canada and 12,000 in the US [1, 25]. A variety of treatments for traumatic SCI have been tested in recent decades, but none have proved widely effective for improving neurological outcomes in humans [49]. The heterogeneity of the SCI population is one possible reason for the lack of effective treatments, in that we do not fully understand the effects of important variables such as age, injury level, severity and mechanism [113]. Further research into these variables is necessary to guide substantial breakthroughs in targeted therapy development. One particular aspect that has received little attention until recently is the possibility of important differences in injury patterns created by the mechanism of primary injury – such as a spinal cord contusion from vertebral burst fracture or a fracture-dislocation, the two most prevalent clinical mechanisms. A second factor that is thought to lead to differences in cord injury pattern is the variation of injury velocity. Differences in cord injury patterns could have implications for differential treatment of patient groups, and can most thoroughly be investigated with a combination of experimental and computational approaches. My thesis continues a line of research at the Orthopaedic and Injury Biomechanics Group (OIBG) at the University of British Columbia that has aimed to investigate the primary response of the spinal cord to mechanical insult through the use of computer and animal models. Specifically, my thesis focuses on the question of how the mechanism and velocity of injury influence the strain distribution during SCI. Several previous projects by OIBG alumni were of immediate influence to the motivation and approach taken in my work, and are briefly discussed next as part of the background to my thesis. 1.1.1 Human FE model Inspiration for the computational aspect of my project came from a finite element (FE) model of the human cervical spine created previously by Greaves et al. [40] using the FE software ANSYS. The model included levels C4-C6 of the human spine and linear elastic material properties for the cord to simulate three distinct injury mechanisms: contusion, dislocation, and distraction 1.1. Greaves’ model showed distinct strain patterns for the three mechanisms, encouraging further research in this area. However, validation of the human model was complicated by the fact that little experimental data is available for the human spine. On the other hand, rats are frequently the subject of SCI research experiments due to the relatively low associated cost and the fact that SCI pathology in the rat closely resembles that in the human [71, 113]. Such experimental data can be used to validate an FE model of the rat spine and allow even more biofidelic simulations. 1 1.1. Background and motivation Figure 1.1: Mid-sagittal view of von Mises strains for (a) transverse contusion, (b) distraction, and (c) dislocation injuries. Transverse view of von Mises strains for (d) transverse contusion injury, at a level adjacent to the cranial edge of the indenter, (e) distraction injury, at a level adjacent to the middle of the C5 vertebral body, and (f) dislocation injury, at a level adjacent to the caudal edge of the C5 lamina. Arrows indicate approximate areas of contact with the dura mater. [Figure and caption text reproduced with permission from Greaves et al. [40].] 1.1.2 Injury mechanism experiments Continuing the line of investigation into SCI differences according to mechanism of injury, Choo et al. [14] developed an experimental rat model of dynamic (100 cm/s) contusion, dislocation, and distraction. This was the first experimental model to systematically compare and contrast multiple injury mechanisms. Cord tissue damage and rostral-caudal extent was assessed by staining histological slices with fluorescein-dextran to identify cord cells – axons in the white matter and neuronal somata in the gray matter – permeable to the marker, indicative of cell membranes ruptured during injury. Figure 1.2: Illustrations of experimental injury configurations. To model dislocation (A), the rostral (left) vertebral clamp was held stationary while the caudal (right) clamp was coupled to the actuator for dorsal translation. For distraction injuries (B), C3 and C4 were held stationary while C5 and C6 were translated caudally. In the contusion model (C), the vertebral clamp holding C4 and C5 was supported while the 2mm spherical head impactor injured the cord through a laminectomy. [Figure and caption text reproduced with permission from Choo et al. [16].] 2 1.1. Background and motivation A notable observation from these experiments was more widespread white matter damage being found in dislocation compared to contusion injuries. Simulation of these same injury mechanisms may shed light on this and other differences. Furthermore, some of the data recorded in these experiments could be used to help validate a FE model of the rat spine. Specifically, the recorded spinal cord indenter displacement and applied force time histories can be used for this purpose. A validated FE model could then be used to compare and possibly correlate simulation computed internal strains of the cord with Choo’s histological measures of tissue injury. 1.1.3 Injury velocity experiments To investigate the influence of injury velocity on cord tissue damage, Sparrey et al. [107] used an experimental model of 1 mm contusion in the thoracic rat spine. Contusions were performed at 3 and 300 mm/s to capture differences over a wide velocity range. Histological results showed increased white matter damage at high velocity compared to low (Figure 1.3). This provided some evidence that there may be an injury velocity threshold for damage to the white matter. Figure 1.3: Parasagittal sections of the spinal cord demonstrating the hemorrhage resulting from contusion for the control (A), slow (B), and fast (C) groups. Tissue was stained with hematoxylin and eosin. Sections are oriented with dorsal surface on the right and caudal aspect towards the top. [Figure and caption text reproduced with permission from Sparrey et al. [107].] Recordings of force and displacement during the experiments also demonstrated stark differences in cord stiffness exhibited during the slow and fast contusions (Figure 1.4). This is indicative of the viscoelastic, rate-dependent material properties of the spinal cord, and is important to capture in simulations within this injury velocity regime. Figure 1.4: The load-displacement curves for the slow (red) and fast (blue) groups. The data shows good repeatability for each contusion and the slope of each line represents the stiffness response of the spinal cord to loading. [Figure and caption text reproduced with permission from Sparrey et al. [107]] 3 1.2. Project definition 1.2 Project definition The previous work described above motivated the development of a finite element model of the rat cervical spine to further investigate the influence of injury velocity and mechanism on spinal cord injury, as well as the relationship between cord tissue injury and strain. Such models may one day aid the design of preventative or emergency treatment devices, but we must first use them alongside experimental methods to gain a better understanding of how strain in the cord is related to observed tissue injury. 1.2.1 Objectives The objectives of my research were to: ❼ Create a biofidelic dynamic, nonlinear finite element model of the rat cervical spine. ❼ Simulate spinal cord contusion experiments at impact velocities of different orders of magnitude and compare to experimental results. ❼ Simulate dynamic spinal cord injury experiments for contusion and dislocation injury mechanisms and compare FE strains to tissue damage. ❼ Validate the FE model by comparing computed injury forces to experimentally measured values. 1.3 Anatomy of the rat cervical spine As vertebrate mammals, the rat cervical spine has much in common with that of the human1 . Accordingly, much of the scientific literature on the anatomy of the human spine is useful in understanding that of the rat, for which there is less published material. However, the rat spine is not simply a scaled down version of the human spine, and in order to accurately model the former some literature specific to the rat is required [32]. In particular, several books and papers on the anatomy of the rat were consulted to aid in the following description of the rat cervical spine [42, 47, 53, 93, 118, 120]. For clarity of discussion, Figure 1.5 demonstrates the spatial terminology used when describing various aspects of the rat anatomy. Figure 1.5: Spatial terminology with respect to the rat [Illustration reproduced from Wingerd [120].] 1 Much of the text in this section first appeared in Russell [94], and is reproduced here for reference. 4 1.3. Anatomy of the rat cervical spine Figure 1.6 shows the vertebrae contained in the rat cervical spine in relation to the rest of the rat skeleton. The seven cervical vertebrae (red) are seen immediately caudal to the skull in the order C1-C7, followed by the first two thoracic vertebrae (blue), T1 and T2 for reference. Figure 1.6: Diagram of a rat skeleton [Illustration reproduced from Muskopf [78] under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 United States License.] Each of these vertebrae has the following basic components (see Figure 1.7): a central body, or centrum; a neural arch which extends from the centrum to create a neural canal axially through the middle of the vertebra; a spinous process which extends dorsally from the neural arch; and articular processes called zygapophyses which form joints between vertebrae [120]. (a) Cervical vertebra (b) Thoracic vertebra Figure 1.7: Diagrams of generic rat vertebrae [Illustrations reproduced from Wingerd [120].] The first two cervical vertebrae are unlike the others in appearance, and have special names. The most cranial is the Atlas, which articulates with the base of the skull. The second is the Axis, with the identifiable odontoid peg about which the Atlas rotates to allow turning of the head, as well as a very pronounced spinous process. Common to all cervical vertebrae are the transverse foramina, or vertebrarterial canals, which house local arteries and veins2 . The sixth cervical vertebra is unique in having two extra ventral processes, the carotid tubercles. The thoracic vertebrae are characterized by long spinous processes, except for T1. Adjacent vertebrae are connected via a network of ligaments and, more distinctly, via an intervertebral (IV) disc. The discs are a mixture of water and fibrous cartilage, with a central nucleus pulposus having a higher water content than the more rigid annulus fibrosus which composes the remainder. Each disc is located between the centra of adjacent vertebrae, and is fused directly to the vertebral bone. Within the protective walls of the vertebral canal lies the spinal cord. The cord itself consists of both gray and white matter, the former making an inner butterfly shape in a cross-section of the cord. Surrounding the cord are the sheathing layers of the meninges – the pia mater, arachnoid, and dura mater (see diagram in Figure 1.8). The subarachnoid cavity between the pia mater and arachnoid layers is occupied by the cerebrospinal fluid (CSF), in which the cord is suspended, and denticulate ligaments that loosely tether the cord to the pia mater. 2 The vertebrarterial canal may be small or absent in the seventh cervical vertebra [118]. 5 1.3. Anatomy of the rat cervical spine Figure 1.8: Diagrammatic transverse section of the medulla spinalis [human spinal cord] and its membranes. [Illustration and caption text reproduced from Gray [39], now in the public domain.] Paired nerve roots branch off from the spinal cord at each spinal level in between vertebra to innervate various parts of the body, depending on the level. In the cervical spine, the nerve roots are named for the lower vertebra of the two-vertebra segment that it runs between. 1.3.1 Anatomy of the spinal cord The distinct white and gray matter of the spinal cord each form a unique and critical part of the central nervous system (CNS). Both parts contain and support neurons, the main functional cells of the CNS (Figure 1.9). The gray matter of the cord contains the neuronal cell bodies and other supportive cells. Neuronal cell bodies send and receive information by passing electrochemical signals, or action potentials, to neighboring neurons via a network of branched projections, called dendrites. Neurons conduct these action potentials far distances throughout the body along their wire-like portions, or axons. The spinal cord white matter is comprised largely of bundles of these axons, surrounding the gray matter and exiting or entering the cord at a nerve root level corresponding to specific function, as well as myelin sheaths wrapped around the axons. Dendrite Axon Terminal Cell body Axon Nucleus Node of Ranvier Schwann cell Myelin sheath Figure 1.9: Diagram of a typical neuron. [Reproduced under the GNU Free Documentation and the Creative Commons Attribution-Share Alike 3.0 Unported Licenses (http://commons.wikimedia. org/wiki/File:Neuron.svg).] The butterfly shape of the gray matter is marked by the dorsal and ventral horns at its tips. The gray matter can be further subdivided into anatomical and physiological regions including areas or strips termed laminae and other smaller zones (Figure 1.10). These regions are marked 6 1.3. Anatomy of the rat cervical spine by cell density and morphological and functional differences between the cells present, in general with larger neuronal cell bodies located within the ventral gray matter. A much higher density of vasculature in the gray matter reflects the higher demands of the cell bodies located there. Figure 1.10: Cross section and diagram of the rat spinal cord at C4. [Images and structure abbreviation list reproduced from Watson et al. [117].] The white matter is not subdivided into as many small zones, but is comprised of several tracts according to function of the axons passing through each area, with axons within a tract all having the same origin, course and termination. Related tracts are referred to as a pathway, more broad groupings of axons within the white matter are termed funiculi, while small groupings of axons with some commonality are referred to as fasciculi. The organization of the spinal cord varies with spinal level, as nerves enter or leave the cord via dorsal nerve roots, which consist largely of afferent sensory fibers, or ventral roots, which include efferent motor fibers. The cervical spinal cord contains the widest range of functional tracts including those serving both the upper and lower body and limbs; injury to the upper cervical cord is thus the most debilitating, with the possibility of quadriplegia and loss of breathing control, depending on injury severity. 7 1.4. Experimental mechanisms of spinal cord injury 1.4 Experimental mechanisms of spinal cord injury A variety of experimental mechanisms have been used to study spinal cord injury using animal models and a brief overview of these is given below. Transection is a common injury mechanism in studies of spinal cord injury recovery as it is relatively simple to implement and allows for precise control of functional deficits depending on the location and size of the cut to the cord [11, 21, 61, 64, 91, 108]. The mechanism involves precise surgical cutting of the cord and can involve either a complete transection – to fully disrupt all axons at a specific spinal level and allow for clearer interpretation of regeneration across the injury – or a partial transection – to allow disruption of specific spinal cord tracts while enabling functional comparison with the uninjured contralateral side. While transection injuries do not reflect the nature of the majority of clinical SCI mechanisms, they are especially useful for precise studies of axonal regeneration [62]. Clip compression is a mechanism that has been used to model sustained cord compression that can result from residual compression following traumatic SCI [26–28, 44, 92, 97]. The clip compression mechanism is created by applying a modified aneurysmal clip to compress the cord, and does not allow for direct measurement of the mechanics of the injury applied to the cord. A contusion injury is used in the majority of animal studies of spinal cord injury [62], and it is an important clinically relevant mechanism. The contusion is defined by the use of a controlled indentation from a rigid impactor, or indenter, to hit the cord surface. This corresponds to burst fracture, a common clinical mechanism, by simulating a vertebral bone fragment impacting the cord. The most widely used implementation of contusion injury is the New York University (NYU) impactor, which uses a 10 g rod dropped from heights of 6.25-50 mm onto the rat thoracic cord to induce graded levels of SCI [6, 43]. This weight-drop contusion mechanism allows for the calculation of velocity at impact and energy delivered to the spinal cord and can allow measurement of the impactor displacement during injury, but does not allow for the measurement of the applied force. Another common contusion device is the Ohio State University (OSU) impactor, an electromagnetic displacement feedback-controlled device to deliver controlled impacts to the thoracic cord, with recording of both the applied displacement and force during injury [9, 52, 79, 80, 109, 110]. A modified version of the OSU impactor was used by Sparrey et al. [107] to create contusions over a wide range of velocity (as mentioned in Section 1.1.3). Other variations of contusion injury mechanisms include those using a pneumatic impactor [2, 85, 99, 121, 126] or using a force-controlled electromagnetic actuator [96]. Fiford et al. [29] were the first to develop a vertebral dislocation model of SCI in the rat, to reflect another clinical SCI mechanism. They conducted lateral dislocations which they found caused greatest axonal injury in the left lateral white matter (where they expected greatest tensile strain in the cord) and vascular injury concentrated within the lateral gray matter, differing from the central injury cavitation typical in contusion injury. As discussed in Section 1.1.2, Choo et al. [14, 15] developed a high speed injury device that they used to investigate primary and secondary cord damage from contusion, anterior dislocation, and distraction of the rat cervical spine. Clarke et al. [18] expanded on this to show that anterior dislocation injury in the rat thoracolumbar spine is more severe than lateral dislocation, in line with clinical observations and emphasizing the utility of studying multiple SCI mechanisms. While critical for examining actual tissue damage during SCI, a key limitation of experimental methods in investigating spinal cord injury mechanics is that they are not able to elucidate internal patterns of stress or strain within the cord during dynamic injuries. Blight and Decrescito [8] began addressing this by using a gelatin surrogate spinal cord and tracing cord deformation during contusion via ink lines injected into the cord – such a model is useful for developing theories of cord strain during contusion that may help explain injury patterns, but is nevertheless a simplified approximation to the material properties of the cord. Recent work by Lucas [65] tracked internal 8 1.5. Modelling of cord injury deformation of the in vivo rat cord during 130 mm/s contusion using high-speed x-ray imaging of fiducial markers, revealing some aspects of internal cord strain but unable to show the full strain distribution. Certainly, this is one area where finite element models of the spine and cord serve as excellent complements to experimental methods. 1.5 1.5.1 Modelling of cord injury Strain theory In order to precisely quantify material deformations, the quantity of strain can be defined at small elements throughout the material. Element strain essentially describes deformation of that element independent of rigid body motion (motion in which the element translates and/or rotates without changing shape). This is precisely the component of motion that is useful for analysis of material failure, since rigid body motion itself has no means to effect failure. With a material modelled as a continuum, a two-dimensional infinitesimal rectangular element of initial dimensions dx by dy can be used to demonstrate the strain resulting from deformation of the material at that point (Figure 1.11). ux y dy d uy y c dy uy(x, y+dy) b C D uy x a y dx dy uy(x, y) A ux x B dx dx ux(x, y) ux(x+dx, y) x Figure 1.11: 2D strain geometry for an infinitesimal material element [Diagram taken from the public domain.] The normal strain of the element in the x-direction is defined as the fractional change in side lengths3 to be ab − AB . x = AB Assuming small strains, as typically encountered by most engineering materials, we have ab ≈ x dx + ∂u ∂x dx, and given AB = dx the strain reduces to x = ∂ux . ∂x , z Similarly for the other directions we have y 3 = ∂uy ∂y = ∂uz ∂z . Note that strain is a unitless quantity as it is expressed as a ratio of lengths. 9 1.5. Modelling of cord injury The change in angle between two of the element’s sides is quantified by the engineering shear strain as γxy = α + β. These angles are defined by the geometry of Figure 1.11 as ∂uy ∂uy dx ∂x ∂x tan α = = ∂ux ∂ux dx + dx 1+ ∂x ∂x , ∂ux ∂ux dy ∂y ∂y tan β = = ∂uy ∂uy dy + 1+ dy ∂y ∂y . ∂ux Again assuming small strains, we know 1 + ≈ 1 and we can further use the small angle ∂x approximations of tan α ≈ α and tan β ≈ β to simplify the angle definition to α= ∂uy ∂x , Therefore γxy = γyx = β= ∂ux ∂y . ∂uy ∂ux + ∂x ∂y and similarly, γyz = γzy = ∂uy ∂uz + ∂z ∂y , γzx = γxz = ∂uz ∂ux + ∂x ∂z . The full strain tensor contains nine components and is expressed in matrix form as γxy /2 γxz /2 xx xy xz xx γyz /2 = yx yy yz = γyx /2 yy γzx /2 γzy /2 zz zx zy zz ∂uy ∂ux ∂uz 1 ∂ux 1 ∂ux + + ∂x 2 ∂y ∂x 2 ∂z ∂x 1 ∂uy ∂uy ∂ux ∂uz 1 ∂uy = 2 ∂x + ∂y . ∂y 2 ∂z + ∂y ∂uy ∂ux ∂uz 1 ∂uz 1 ∂uz 2 ∂x + ∂z 2 ∂y + ∂z ∂z (1.1) (1.2) The strain tensor can alternatively be defined by specifying the component in its i’th row and j’th column as ∂uj 1 ∂ui + , (1.3) ij = 2 ∂xj ∂xi with xi and xj for i or j = 1, 2, 3 corresponding to the axes x, y, z described above. Note that the orthogonal axes x, y, z used in the definition of the strain tensor’s components are arbitrary, and depending on the choice of axes the value of those components is different. There is one particular set of axes which in fact reduces the strain tensor to have zero components everywhere except along the diagonal (1.4). These three diagonal normal strains are unique and referred to as the principal strains, with the corresponding axes called the directions of principal strain. 0 0 1 = 0 2 0 (1.4) 0 0 3 The principal strains are the eigenvalues found by solving the linear algebraic equation ( − i I) ni = 0, 10 1.5. Modelling of cord injury with ni being the eigenvectors corresponding to the directions of principal strain. The principal strains represent the largest magnitude of normal strains, since there are no shear strain components in that set of axes. The maximum principal strain, also called the first principal strain or 1 , is strictly defined as the largest of the three principal strains and is typically tensile (positive valued), or stretching, in nature for most materials and deformation states. The minimum or third principal strain, 3 , is similarly the smallest and is typically compressive (negative valued), or shortening. The second principal strain is in between the other two, and can be either tensile or compressive depending on the deformation state. The relationship between material strain and stress – the distribution of applied forces throughout the material – is referred to as the constitutive equation or material model. Some of these models appropriate for soft tissue such as the spinal cord are discussed in Section 1.5.3. The parameters associated with a material model, which are unique for a specific material being modelled, are referred to as the material properties. For large strains (>5%), which soft biological tissues such as the spinal cord are often subjected to, the finite strain theory framework and a more complex nonlinear strain tensor must be used. This is the Green-Lagrangian strain tensor and is defined as Eij = 1 2 ∂uj ∂ui ∂uk ∂uk + + ∂xj ∂xi ∂xi ∂xj . (1.5) Note that U (components of which appear in the partial derivatives of 1.5) is the tensor mapping the applied stretches to the undeformed configuration, and is followed by a rigid body rotation by R to result in the deformed configuration (Figure 1.12). This strain tensor, Eij , differs from the k ∂uk small strain tensor (1.3) by including the nonlinear product ∂u ; the small strain approximation ∂x ∂x i j is a linearized form of the tensor valid for small strains only, while the Green-Lagrangian strain tensor is exact for any strain value. Another difference is that the Green-Lagrangian strain tensor partial derivative terms are taken with respect to the undeformed configuration (x ), rather than the deformed configuration (x) as in the small strain approximation (1.3), a difference which cannot be neglected at large strains. This finite strain theory formulation leads to a more complicated procedure to solve for a strain field within a material, which explicit finite element solvers are well suited to follow. Nevertheless, the concept of principal strains as unique eigenvalues capturing the essence of the element deformation associated with a given strain tensor still holds in this framework. 1.5.2 Material properties of the spinal cord Characterizing material properties for the spinal cord has long been regarded as important for understanding the biomechanics of spinal cord injury. While various attempts have been made over the years, difficulties associated with accurately and repeatably measuring these properties has prevented reaching a consensus on their values that would cover a wide range of conditions. In particular, only tensile tests have been performed widely enough for comparisons across studies, and these have only been performed at low strain rates (0.001-0.3 s−1 ) and peak strains (<0.1) compared to those typical of traumatic SCI (strain rates >5 s−1 and peak strains >0.2) [12, 17]. Some of the challenges of material testing of the cord include maintenance of testing conditions such as hydration and temperature to mimic the in vivo environment, limitations of standard mechanical testing equipment to strain rates much lower than traumatic loading conditions, characterization of boundary conditions due to the interface of tissue with the testing apparatus, and the effect of specimen characteristics such as species or age on properties. Cheng et al. [13] recently highlighted the effects of preconditioning on mechanical test results of the cord at various strain rates and magnitudes, a factor which had not been addressed previously in the literature and may 11 1.5. Modelling of cord injury Figure 1.12: 3D framework for finite strain theory. U is the tensor mapping the stretches applied to the undeformed configuration, and is followed by a rigid body rotation by R to result in the deformed configuration. The overall deformation is described fully by the deformation gradient, F = RU. [Diagram taken from the public domain.] have contributed considerably to variability between studies with different protocols. Furthermore, accurate quantification of in vivo cord tissue properties from cadaveric specimens is complicated by the fact that cord properties change in the hours to days after death, with the tangent modulus (the stiffness at high strains) increasing by >50% over 72 hours [81]. Another complication is that at high strains localized tissue failure, or damage, may occur, and this damage should be characterized and modelled to fully simulate the cord mechanics during SCI. Some investigation has been done on damage in high strain testing of brain tissue Darvish and Crandall [23], Prange and Margulies [89], Shafieian et al. [101], but no similar work has been done for the spinal cord. Despite these challenges, several measurements of cord tensile properties have been conducted, with the compiled results shown in Figure 1.13. In general, the stress-strain response of the cord is nonlinear over typical ranges of strain and strain-rates. This nonlinearity includes both hyperelastic behaviour – in which the slope of the stress-strain curve, or stiffness, increases with higher levels of strain – and viscoelastic behaviour – in which the stress at a fixed level of strain decays or relaxes over time [7, 10, 17, 19, 31, 48, 81]. Both hyperelasticity and viscoelasticity are characteristics that cause a material to deviate from linear elastic behaviour, for which the stiffness of the material is constant over a wide range of strain and a unique stress value corresponds to a given level of strain. 1.5.3 Material modelling of the spinal cord With simulating deformation of the spinal cord of chief importance to modelling spinal cord injury, an appropriate material model to govern that deformation is required. Several earlier models of the spine used linear elastic material properties for the spinal cord as a first step towards modelling cord deformation and injury [40, 63, 98]. However, material testing data of the spinal cord as well as brain tissue have shown that these tissues exhibit clear hyperelasticity and viscoelasticity [7, 12, 19, 31]. In addition, there is evidence that finite element simulations should model this hyperviscoelasticity as linear elastic cord properties are not sufficient to accurately model stresses and strains within the cord [106]. Hyperviscoelastic properties of the cord are especially important 12 1.5. Modelling of cord injury (a) Hyperelastic stress-strain measurements (b) Viscoelastic relaxation measurements Figure 1.13: Hyperelastic and viscoelastic properties of the spinal cord. Data are compiled from the different studies of cord properties for various peak strains and strain rates. [Figures reproduced with permission from a review by Clarke [17].] when simulating dynamic SCI mechanisms, in which the cord is subjected to large strains at high rates. Some hyperviscoelastic models have been used previously with hyperelasticity based on a polynomial strain energy function and deformation tensor invariants [36, 60, 70, 72]. However, these models were discouraged in favour of a less restrictive – in that it allows different material behaviour in tension versus compression – Ogden hyperelastic model generalized to incorporate Prony series viscoelasticity as proposed by Miller and Chinzei [74] and based on the quasilinear viscoelastic theory of separable hyper- and visco-elastic model components introduced by Fung [34]. The model is practical for finite element simulation as it is currently implemented in most FE solvers, including PAM-CRASH. The mathematical constitutive equations describing each of the Ogden and Prony series material model parts are described in the following subsections. Ogden hyperelasticity Ogden [82] first proposed his theory of hyperelasticity to model incompressible rubberlike solids, and the model has since been used extensively to model biological soft tissues. The Ogden model defines the strain energy density, W , in terms of the principal stretches4 λj , j = 1, 2, 3 as5 : N W (λ1 , λ2 , λ3 ) = 2 i=1 µ i αi (λ + λα2 i + λα3 i − 3) αi 1 (1.6) where N is the order of the model, αi are material constants describing the hyperelastic nonlinearity, and µi define the material shear modulus as µ = N i=1 µi αi . Note that the Ogden model degenerates into the less general Neo-Hookean hyperelastic model for N = 1, α = 2 or the Mooney-Rivlin model for N = 2, α1 = 2, α2 = −2. The principal stresses are then derived from differentiating the strain energy function according to ∂W σj = λ j − p, (1.7) ∂λj where p is a Lagrange multiplier associated with the material incompressibility constraint λ1 λ2 λ3 = 1. Furthermore, to ensure stable behaviour during deformation, the Ogden parameters must satisfy N i=1 µi αi > 0 [24]. 4 Note that the stretch ratio, λ, is defined in terms of the current and initial material sample lengths as λ = 0 and is related to the normal strain by = l−l = λ − 1. l0 5 This definition uses the coefficient convention of PAM-CRASH [24]. l , l0 13 1.5. Modelling of cord injury The model can be further simplified for the case of uniaxial tension, which is useful when fitting −1/2 experimental tissue tests [83]. For uniaxial tension λ2 = λ3 = λ1 and with λ1 = λ the definition of principal stresses reduces to that of a single tensile stress defined by 1 σ(λ) = 2µ λα−1 − λ− 2 α−1 . (1.8) Prony series viscoelasticity The Prony series model of viscoelasticity, also known as the Generalized Maxwell model or the Maxwell-Wiechert model, is the most general form of linear viscoelasticity. Linear viscoelastic models are those that assume separable elastic and viscoelastic responses. Such models yield general equations for stress or strain as a function of time of t F (t − t ) ˙(t )dt σ(t) = Einst,relax (t) + (1.9) 0 and (t) = σ(t) Einst,creep t ˙ )dt K(t − t )σ(t + (1.10) 0 where t is time, σ(t) is stress, (t) is strain, Einst,creep and Einst,relax are instantaneous elastic moduli for creep and relaxation, K(t) is the creep function and F (t) is the relaxation function. Creep is the phenomenon in which tissue strain will increase, or creep, over time when subjected to a constant stress. Relaxation is the inverse process, in which tissue stress decreases, or relaxes, over time when subjected to a constant strain. k1 k2 kj 1 2 j ke Figure 1.14: Prony series model schematic. A series of spring-dashpot Maxwell elements is shown each arranged in parallel with a lone spring. Each dashpot damping element is associated with a relaxation time constant, τi . [Diagram taken from the public domain.] The Prony series models linear viscoelasticity by recognizing that relaxation may not occur at only one time, but at a distribution of times limited only by the order of the series. This is achieved by a sufficiently long series of spring-dashpot Maxwell elements arranged in parallel with a lone spring element to adequately model the viscoelastic behaviour of a material at all relevant timescales (Figure 1.14). As often employed in mathematical and finite element models, the Prony series definition for relaxation of the shear modulus is M Gi e(−t/τi ) , G(t) = G∞ + (1.11) i=1 14 1.5. Modelling of cord injury where G∞ is the long term or steady-state shear modulus, M is the order of the model, Gi is the ith shear modulus component and τi is the associated time constant, or relaxation time, for viscoelastic decay of that component. The equation can alternately be arranged by noting that the instantaneous elastic modulus at time zero is related to the long term modulus by G(t = 0) = G0 = G∞ + Σ M i=1 Gi , leading to M Gi 1 − e(−t/τi ) . G(t) = G0 − (1.12) i=1 A further variation that is sometimes used to express the relaxation parameters is M G(t) = G0 γi 1 − e(−t/τi ) 1− , (1.13) i=1 where γi = 1.5.4 Gi 6 G0 . Finite element modelling of the spine, spinal cord, and brain Thanks to the increasing computational power of personal computers, the last decade has encouraged the development of finite element models in all areas of engineering, and the area of biomechanics is no exception. A range of FE models incorporating either the brain, spinal cord, or both have been developed by groups around the world, and a summary of those relevant to the proposed project will now be presented. Scifert et al. [98] investigated spinal cord mechanics in flexion and extension of the cervical spine using a linear elastic model of the C5-C6 motion segment created in ABAQUS. The model is unique in including nerve roots attached to the spinal cord, though Scifert et al. [98] did not discuss possible influence of this inclusion on results. As mentioned in Section 1.1.1, Greaves et al. [40] also created a linear elastic model of the human cervical spine at levels C5-C6 and simulated contusion, dislocation and distraction injuries. Another linear elastic model by Li and Dai [63] specifically explored hyperextension injury with an isolated FE model of the cord in ANSYS, with model validation performed against static cord compression and axial tension experimental data. Galle et al. [36] created a 2D model of the spinal cord in Matlab/COMSOL formulated to match experimental results of compression of guinea pig cord strips. Their cord model used a hyperelastic Mooney-Rivlin strain energy function with no viscoelastic component. Another 2D model of the spinal cord was used by Ichihara et al. [50] to model their observed differences in white and gray matter properties. They modelled the cord as nonlinearly elastic directly using measured material stress-strain curves with the gray matter generally more stiff than the white matter at higher strains. Applying a 30% quasi-static compression to this cord model, [50] then compared cord deformation to MRI imaging of their experimental results and found better agreement for the non-homogeneous white and gray properties compared to a simulation with homogeneous properties. Viano et al. [115] modelled concussion injury in professional football using a FE model developed in PAM-CRASH, comparing clinical symptoms with simulation results for 28 impact cases. They modelled the brain, brainstem and cerebellum using a Kelvin viscoelastic model, neglecting hyperelastic behaviour. Kleiven [60] developed a hyperviscoelastic brain and head model using LS-DYNA including a Mooney-Rivlin hyperelastic model for the brain combined with a first order Prony series for viscoelasticity. Ho and Kleiven [45, 46] further refined this model to include 3D vasculature within the brain, including nonlinear elastic properties for veins and arteries based on the uniaxial exponential model proposed by Fung [34]. 6 Note that sometimes the symbol gi is used instead of γi . 15 1.6. Mechanical indicators of tissue injury Another LS-DYNA model created by Kimpara et al. [59] included the human head and neck to investigate injury mechanisms during severe frontal impacts. Brain tissue was modelled as viscoelastic while the spinal cord was modelled as hyperelastic, with slightly different material properties for the cord white and gray matter based on direct use of the stress-strain curves obtained from tensile testing of bovine spinal cords by Ichihara et al. [51]. Kimpara et al. [59] additionally performed material testing of the porcine cervical pia mater and included this layer in their model. One of the most recent and relevant contributions is that of Maikos et al. [68], who created the first finite element model of the rat spine. Their explicit hyperviscoelastic model, created in Abaqus, simulated weight-drop contusion experiments in the thoracic spine at 0.49 and 0.69 m/s. Based on an experimental image of tissue damage indicated by albumin extravasation, Maikos et al. [68] then marked elements as injured or uninjured and correlated injury status to maximum principal strain with logistic regression. Correlations with injury for their model were very good in the gray matter and fair in the white matter. Limitations common to many FE models of the spinal cord include validation against quasistatic test conditions that do not reflect dynamic injuries being simulated, oversimplification of geometry or simulation of a very small portion of the spine, and use of cord material properties that do not reflect the in vivo behaviour of the spinal cord during dynamic injury mechanisms. These limitations are, however, beginning to be addressed and overcome. In particular, recent models of the spinal cord or brain typically model these soft tissues as hyperviscoelastic materials. 1.6 Mechanical indicators of tissue injury A variety of different mechanical indicators of tissue injury have been investigated, including kinematic characteristics of injury and finite element simulation results. The goal of such indicators is to be highly correlated with tissue damage so that injury predictions can be made more broadly, or so that computational models can use them to compare different injury mechanisms and severity levels and yield clinically relevant conclusions. Kearney et al. [57] first probed the effects of varying contusion magnitude and velocity, finding that similar functional injury severity could be achieved by either large contusion magnitude or smaller magnitude contusion at high velocity, indicating the importance of both factors. Viano and Lovsund [116] later conducted dynamic cord contusion experiments in a ferret model at velocities of 1.5-6 m/s and displacements of 1.25-3.25 mm (25-65% compression), finding a good correlation between graded SCI and the maximum viscous response, V C, of the injury mechanism defined as the multiplication of velocity (V ) and percent compression (C) of the cord diameter. Other investigators of the relationship of velocity with injury pattern include Jakeman et al. [52] who did not find an association between contusion velocity and behavioural outcome score (assessed with the Basso, Beattie and Bresnahan (BBB) locomotor rating scale [5]) but only varied velocity slightly (14-19 cm/s); Kim et al. [58] who found no significant difference in Basso mouse scale injury severity for 0.8 mm contusions at velocities ranging from 0.1-0.4 m/s; and Maikos and Shreiber [66] who predicted a white matter vasculature injury threshold of 200 mm/s but no gray matter threshold in weight-drop contusions with velocities of 0.5-1 m/s. On a cellular level, Cullen and LaPlaca [22] demonstrated that neural tissue cultured in 3D exhibits more complex loading patterns and different injury thresholds compared to tissue cultured on a 2D plate, underlying the importance of studying spinal cord tissue properties and behaviour in an in vivo setting. Geddes-Klein et al. [37] subjected cultured neurons to uniaxial and biaxial stretches of 0-50% total strain, observing distinct neurophysiological responses to the different injury mechanisms. 16 1.7. Summary Maximum principal strain Several studies have investigated maximum principal strain during brain or spinal cord injury and found it to be a good indicator of tissue damage, making it one of the most widely used indicators and a useful output from finite element models. Shreiber et al. [103] first quantified maximum principal strain in a finite element model of cerebral contusion in the rat, and found it to be a good predictor of damage to the blood-brain barrier (BBB) over a range of loading conditions. Their results also indicated a strain threshold of ∼18.8% below which damage would not be expected. Bain and Meaney [4] used in situ material testing of white matter tissue from the guinea pig optic nerve to quantify maximum principal strain during applied axonal stretch injuries. These results were then compared with assessed morphological tissue injury and electrophysiological impairment from parallel in vivo injuries. This yielded predicted maximum principal strain thresholds of 0.21 for morphological tissue damage and 0.18 for electrophysiological impairment. Zhu et al. [127] investigated non-impact, graded axial rotation injuries in a pediatric pig brain model, finding periods of unconsciousness ranging from 0 to 80 minutes depending on severity. Brain tissue sections were also stained with neurofilament antibody (NF-68) to identify regions of axonal damage. Finite element recreations of the injury grades yielded strain and strain rate throughout the tissue, and Zhu et al. [127] also looked at the product of strain and strain rate. Volume fractions of the tissue showing strains higher than the level they found to predict tissue injury with 90% probability were well correlated with global injury severity assessed by duration of unconsciousness. As mentioned in the previous section, Maikos et al. [68] used their finite element model of thoracic weight-drop injury in the rat to find good correlations of elemental maximum principal strain with injury status of the corresponding locations in the cord tissue. This result shows than maximum principal strain is a practical indicator of tissue damage for use in FE of the spinal cord. McAllister et al. [69] recently compared subject-specific FE results with in vivo diffusion tensor imaging of subjects with diagnosed concussion, finding maximum principal strain and strain rate associated with changes in indicators of white matter integrity. Together these studies demonstrate that maximum principal strain has been widely correlated with neural tissue damage in a range of models. It can yield more specific and localized injury prediction compared to more global quantities such as injury velocity or compression depth. Furthermore, there is an intuitive basis for failure due to high maximum principal strain, as many soft tissue and cellular structures in the cord (such as axons, vasculature, or individual cell membranes) can be imagined to fail under tension, but may be more tolerant of the compressive or shear strains typically encountered – imagine a slightly elastic rope as a simplified model for these structures, for which tension would appear to be the most likely failure mode. 1.7 Summary In vitro cell culture experiments have shown that the biological and mechanical mechanisms of traumatic neuronal injury are influenced by mechanical loading patterns [22, 37]. On a larger scale, dynamic injuries in the in vivo rat model have shown varying patterns of tissue damage for different injury mechanisms [14, 18, 29]. Furthermore, finite element models of the human spine have demonstrated distinct stress and strain patterns within the spinal cord that depend on the biomechanical mechanism of injury [40, 63]. While experimental methods have revealed important details regarding tissue injury thresholds and patterns, they are difficult to fully interpret and apply to clinical injuries as they do not yield information on internal spinal cord deformations during SCI. Finite element models are ideal for investigating these deformations. However, human models to date have been difficult to validate 17 1.7. Summary due to a lack of in vivo loading data. In addition, many models have used linear elastic quasi-static simulations which may not capture the full nature of high speed cord injuries. One group began to address this with an experimentally calibrated dynamic, hyperviscoelastic finite element model of weight-drop contusion in the rat thoracic spine but did not investigate other injury mechanisms or vary injury velocity over several orders of magnitude [68]. The overall goal of my thesis is to develop and validate a dynamic finite element model of the rat cervical spine and to use it to compare internal spinal cord deformations – during contusion injuries at different velocities and during both contusion and dislocation injury mechanisms – with previously observed tissue damage. The indicator of neuronal tissue damage used for this comparison was cellular permeability to fluorescein-dextran as membrane permeability has widely been linked to neuronal pathology [35, 102, 111, 112] and has been used to quantify regional patterns of damage in the spinal cord [14]. I chose maximum principal strain as the primary measure of cord deformation because it has been shown to be a good predictor of neural damage for several animal models [4, 68, 103, 127] and as it is relatively easy to interpret (being generally tensile in nature) and link mechanistically with membrane damage. Improved understanding of the strain distribution in the cord during two distinct mechanisms of SCI and at injury velocities of different orders of magnitude will aid interpretation of tissue damage patterns and may inspire new strategies to treat or prevent injury. 18 Chapter 2 Methods 2.1 2.1.1 Model development Geometry extraction from Magnetic Resonance Imaging The first step in development of the rat cervical spine model was acquisition of spinal geometry to be used in creating a finite element mesh. High resolution magnetic resonance imaging (MRI) was performed on a normal, freshly euthanized Sprague-Dawley rat with a 7 T animal scanner (BioSpec 70/20 USR, Bruker BioSpin Corp., Billerica, MA) at the UBC MRI Research Centre (see 2.1). Two scans, oriented perpendicular to the upper cervical cord (C1–C3) and lower cervical cord (C5–T1), were obtained with 156x156 micron in-plane and 1 mm through-plane resolution. These scans were interpolated to an isotropic 156 micron resolution, zero-padded and registered to each other using Analyze (AnalyzeDirect, Overland Park, KS), yielding a fused image of pixel dimensions 256x256x240 covering the full range of the cervical spine. Figure 2.1: 7T animal MRI scanner Image Segmentation Image segmentation was performed to extract object models of the rat spine components from the MRI data. For this purpose an open source software solution devoted to segmentation, ITKSNAP (http://www.itksnap.org)7 , was chosen for it’s simple yet powerful interface, despite the availability of the commercial Analyze software which has segmentation capabilities. ITK-SNAP was created using the Insight Toolkit (ITK), an image analysis software kit designed to support the images of the Visible Human Project , and the Visualization Toolkit (VTK), a 3D data visualization package. SNAP stands for “SNake Automated Partitioning”, referring to the segmentation algorithms employed by the software which make use of active contour and level set methods to partition elements in an image via snake8 evolution [124]. The method of active contour evolution involves estimation of a target object’s boundaries with a closed surface contour which gradually conforms ➤ 7 8 The version of SNAP used in this project was 1.5.2. The term snake here refers to a closed curve or surface. 19 2.1. Model development to those boundaries9 . This evolution in time is modelled by the following partial differential equation (PDE) in 2D: ∂ C(u, v; t) = F n (2.1) ∂t Where, C = closed surface contour parametrized by spatial variables u,v and time, t n = unit normal to C F = sum of forces acting on C in normal direction Of the two active contour methods provided in SNAP, the Region Competition method (called the ‘Intensity regions’ method within SNAP) was found to achieve the desired segmentations relatively quickly and reliably, and was used throughout the semi-automatic segmentation process. This method, pioneered by [128], uses the following definition of the evolution forces: F = α(Pobj − Pbg ) + βκ (2.2) Where, α,β = weight parameters Pobj = probability of voxel belonging to object Pbg = probability of voxel belonging to background κ = mean curvature of C The respective probabilities are assigned to the image voxels using a fuzzy threshold of image intensity performed in SNAP. As demonstrated in Figure 2.2, the seed contour gradually conforms to the desired object topology through the region competition method. Figure 2.2: Active contour evolution using the feature image based on region competition. The propagation force acts outwards over the ‘foreground’ region (red) and inwards over the ‘background’ region (blue), causing the active contour to reach equilibrium at the boundary of the regions. [Figure and caption text reproduced with permission from [124].] The contour evolution problem is solved by means of the level set method of [84, 100], in which the contour is prescribed as the zeroth level set of some function φ, defined at every voxel in the image. Using the relation n = ∇φ/ ∇φ , (2.1) can be transformed to a PDE in φ: ∂ φ(x; t) = F ∇φ (2.3) ∂t SNAP then efficiently solves (2.3) close to the zeroth level set (the level contour corresponding to φ = 0) using the Extreme Narrow Banding Method proposed by [119]. The weight parameters in (2.2) are left up to the user to define for a given situation in order to achieve the desired segmentation result. One of the great strengths of the SNAP Graphical User Interface (GUI) is that it allows the user to respond to the contour evolution process by altering these parameters in real time, enabling intuitive fine-tuning of the object segmentation. Furthermore, the somewhat abstract parameters are displayed alongside the general effect they have on the contour evolution—α is termed the “balloon force” controlling the magnitude of inward or outward force on the contour, and β is the “curvature force” affecting the smoothness of the contour. 9 The following explanation of the active contour method used in SNAP is adapted from that given by Yushkevich et al. [124]. 20 2.1. Model development Figure 2.3: Extraction of rat cervical spine geometry. Geometry of the rat cervical spine was extracted from 7 T magnetic resonance images using ITK-SNAP, a semi-automated volume segmentation tool using 3D snake evolution [125]. A screenshot of the ITK-SNAP interface is shown after segmentation of the white and gray spinal cord, intervertebral discs, and C1 to T2 vertebrae. (Clockwise from top left) Axial, sagittal, and frontal views of the MRI data are displayed, with a 3D cursor that links all three. A 1 cm long scale bar is shown alongside the spine volumes at bottom left. Figure 2.3 shows a screenshot of the SNAP GUI. Axial (top left), sagittal (top right), and frontal (bottom right) views of the MRI data are displayed, with a 3D cursor that links all three. In the sagittal slice, the fusion of the two MRI scans is evidenced by two intersecting rectangles. The segmentation process for the rat cervical spine began with semi-automated segmentation to achieve rough object boundaries, primarily on the basis of contrast differences. Most objects also required significant “clean-up” work applied manually on a slice-by-slice basis, or via SNAP’s 3D paintbrush or cut-plane tools. This manual work often included removing artifacts located outside an object’s expected boundaries, or creating or enhancing specific attributes of an object such as in the creation of vertebrarterial canals. Manual intervention was especially necessary for creation of the zygapophyses at the boundaries between adjacent vertebrae, as well as for creation of the intervertebral discs; the subtle boundaries of these parts prevented them from being accurately segmented automatically, and they had to be estimated based on images from literature. In the case of the intervertebral discs, the nucleus pulposi were distinguishable by contrast in the MRI data, and as such provided a landmark for the disc locations. Geometry for the white and gray matter of the spinal cord, intervertebral discs, and C1 to T2 vertebrae were exported from ITK-SNAP into Rapidform (INUS Technology, Seoul, Korea) and fit with analytical surfaces using non-uniform rational B-splines (NURBS) [88]. Figure 2.4 shows the geometric surfaces segmented from the MRI data10 . A 3D model of these surfaces is also embedded in the electronic version of this document, in Appendix F. 2.1.2 Finite element meshing The segmented surfaces were meshed initially in HyperMesh (Altair Engineering, Troy, MI) using hexahedral solid elements via the solid map tool for the white and gray cord and tetrahedral elements for the discs and vertebrae. Meshes were imported to PAM-CRASH (ESI Group, Paris, 10 Note that the Atlas and T2 vertebrae are slightly incomplete as they were located at the edge of the image data. 21 2.1. Model development (a) Spinal cord (b) Atlas (c) Axis (d) C3 (e) C4 (f ) C5 (g) C6 (h) C7 (i) T1 (j) T2 (k) IV disc (l) 1-cm scalebar Figure 2.4: Surfaces segmented in ITK-SNAP [not shown to scale] France), an explicit finite element software suitable for impact simulations, for further development in the Visual-Crash pre-processor software. The dura mater could not be reliably identified in the acquired MRI scans and was instead created by expanding the surface of the cord based on MR images outlining the CSF by Franconi et al. [33]. The dura was then assigned a thickness of 90 µm, also based on the images with scale provided by Franconi et al. [33], and meshed with two layers of hexahedral solid elements. Finally, spinal ligaments were created by manually defining two-dimensional bar elements according to anatomical descriptions [40]. Figure 2.5 shows the dura mater and spinal ligaments included in the full cervical rat model. dura mater anterior longitudinal ligament (ALL) posterior longitudinal ligament (PLL) Figure 2.5: The dura mater (shown cut along full cervical rat model. joint capsule ligament (JC) denticulate ligaments (DL) ligamentum flavum (LF) dural attachments (DA) interspinous ligament (ISL) the sagittal plane) and spinal ligaments as modelled in the 22 2.1. Model development Ligament areas The ligament parts consisting of bar elements required assigned cross sectional areas in order for the bar properties to be simulated. Since no data on the cross sectional areas of rat spinal ligaments could be found in the literature and physical examination was difficult due to the extremely small and delicate nature of these tissues, areas were scaled down from human values. Cross sectional areas for each ligament were scaled from those used in a human FE model [40] according to the ratio of the rat to human spinal cord cross sectional areas. This ratio was estimated to be 8.6% using an assumption of elliptical cross section and comparing each cord’s anterior-posterior (AP) and transverse diameter at the mid-cervical C3 level (see Table 2.1). Table 2.1: Ratio of rat to human cord cross sectional areas AP diameter, A (mm) Transverse diameter, B (mm) Area of ellipse = 41 πAB (mm2 ) 8.6 2.55 12.1 3.52 81.7 7.0 Humana Ratb Rat:Human ratio 0.086 a Human cord diameters shown here are taken from the model by Greaves et al. [40], the geometry of which was obtained from transverse cryosection images at 1 mm intervals provided by the Visible Human Project (National Library of Medicine). b Rat cord diameters are those of my model, based on the geometry extracted from MRI as described above. Table 2.2 shows the resulting cross sectional areas assigned to each ligament part in the model. The assigned part cross sectional area represents the area of a single bar element. Ligament areas were distributed equally between the bar elements of each ligament (at one spinal level and side, where appropriate). Table 2.2: Cross sectional areas of spinal ligaments Ligament Human Area (mm2 ) Rat Area (mm2 ) # of Elements Element Area (mm2 ) 20 23 46 47 13 0.25 5 1.72 1.978 3.956 4.042 1.118 0.0215 0.43 5 per level 5 per level 1 per level per side 6 per level 1 per level 1 per region per side 4 per level per side 0.344 0.396 3.956 0.674 1.118 0.0215 0.1075 ALL PLL JC LF ISL DL DA 2.1.3 Material properties Spinal cord and dura Biofidelic spinal cord material properties are crucial to yield reliable tissue deformation during FE simulation. While the difficulty of measuring and modelling soft tissue properties has so far precluded gathering enough data for a consensus on cord material properties, material tests of human and animal spinal cord have all demonstrated hyperelastic and viscoelastic behavior [7, 12, 31, 105]. I chose to use the same hyperviscoelastic Ogden and Prony material properties for the cord and dura presented by Maikos et al. [68] for consistency and to further test these material models in different injury conditions (Table 2.3). Cord hyperelastic properties were based on material tests of 23 2.1. Model development rat spinal cord by Fiford and Bilston (2005), combined with viscoelastic properties of brain tissue from Mendis et al. (1995), and then calibrated by Maikos et al. [68] to fit their own weight-drop experimental behavior at impact velocities of 0.489–0.690 m/s. Properties for the dura were derived by mechanical testing of rat dura mater [67]. Appropriate conversions of material constants were performed because of differing notation conventions between PAM-CRASH and Abaqus (refer to Table 2.3), and Selective Reduced Integration was used for dura and cord elements for the optimal balance between accuracy and computational complexity. Table 2.3: Material properties of spinal cord and dura Tissue Hyperelastic Ogden constants Viscoelastic Prony series constants Spinal cord µ = 40.04 kPaa α = 4.7c ν = 0.45c g1 = 0.5282b τ1 = 8 msc g2 = 0.3018b τ2 = 150 msc Dura µ = 207.41 kPaa α = 16.2c ν = 0.45c g1 = 0.3182b τ1 = 9 msc g2 = 0.1238b τ2 = 81 msc g3 = 0.0997b τ3 = 564 msc g4 = 0.0997b τ4 = 4.69 sc a Adapted from Maikos et al. [68] with conversion µP C = G0,AB /α due to difference in notation convention between PAM-CRASH (PC) and Abaqus (AB). b Adapted from Maikos et al. [68] with conversion gi,P C = Gi,AB /G0,AB . c Taken from Maikos et al. [68]. Spinal ligaments Although spinal ligaments do not play a role in experimental cord contusion injuries, they do play an important role in dislocation injuries, taking on part of the load applied to the vertebrae alongside the intervertebral disc. A Nonlinear Tension Only Bar (Material Type 205) material with linear elastic properties was used for all spinal ligaments [40]. Cross-sectional areas for the ligaments were scaled from the human values used in Greaves et al.’s (2008) model, and maximum strains were assigned based on values reported in the literature (Lee et al., 2006; Quinn and Winkelstein, 2007; Yoganandan et al., 2000) (Table 2.4). Linear elastic material properties were used for the spinal ligaments, using the same elastic moduli as used by Greaves [41]. Material Type 205 - Nonlinear Tension Only Bar was used for all spinal ligaments to prevent them from resisting compression. The density of all ligaments was assumed be similar to water at 0.001 mg/mm3 , since no reported density values could be found. Ligament failure strains were modelled based on the stiffness and maximum principal strain at failure for the rat joint capsule ligament (JC) [90]. Similar values for the ALL, PLL, LF, and ISL were based on published values [123]. Intervertebral discs The annulus fibrosus of the C4/C5 intervertebral disc was modelled as linear elastic (Material Type 1, E = 2.4 MPa) according to the properties used previously by Greaves et al. [40] (Table 2.5). The nucleus pulposus was not modelled for this study to reduce complexity, as it is not expected to play a role in the contusion or dislocation mechanisms. Discs were attached to neighboring vertebrae via 24 2.1. Model development Table 2.4: Material properties of spinal ligaments (Material Type 205 - Nonlinear Tension Only Bar Element) Ligament Cross-sectional area per element, A (mm2 ) Youngs modulus, E (MPa) Stiffness factor, k = EA (N) Mass per unit length, µ (g/mm) Failure strain, EPSLN u 0.344a 0.396a 3.956a 0.674a 1.118a 0.0215a 0.1075a 35.2b 35.7b 4.9a 3.8b 5b 5.8b 35.7b 12.11 14.14 19.38 2.56 5.59 0.1247 3.838 0.000344 0.000396 0.003956 0.000674 0.001118 0.0000215 0.000108 0.308c 0.182c 1.51d 0.77c 0.609c 0.087 0.182 ALL PLL JC LF ISL DL DA a Scaled from Greaves et al. [40]. Taken from Greaves et al. [40]. c Taken from Yoganandan et al. [122]. d Taken from Quinn and Winkelstein [90]. b spot welds, which link element nodes between adjacent parts, to simulate the vertebral endplate connection. Preliminary simulations were performed to calibrate spot-weld rupture criteria to achieve simulated disc endplate failure coinciding with experimental failure predicted by Choo et al. [14] based on force history measured during injury; the resulting criteria were ultimate tensile and ultimate shear strengths of 0.15 MPa. Table 2.5: Material properties of intervertebral disc and endplate spotwelds (Material Type 1 - ElasticPlastic for Solid Elements) Structure Disc annulus fibrosus Endplate spotwelds a Youngs modulus, E (MPa) Ultimate tensile strength (MPa) Ultimate shear strength (MPa) 0.15 0.15 3.4a Taken from Greaves et al. [40]. 2.1.4 Fluid-Structure Interaction and the cerebrospinal fluid Persson et al. [86] recently demonstrated the importance of including the incompressible fluid behavior of the CSF in models of SCI, using an ovine FE model with fluid-structure interaction (FSI). Previously, FE models of SCI have omitted the CSF [40, 63, 98], or modelled it as a quasifluid using solid elements [68]. This study proposes the Smoothed Particle Hydrodynamics (SPH) method [75] as an efficient means to include interaction between the cord, dura and CSF in impact simulations. HyperMesh (Altair Engineering, Troy, MI) was used to define simple cubic elements distributed regularly in the volume between the dura and cord elements (mesh pitch of 0.075 mm and ∼153,000 elements). These were converted to SPH point elements in PAM-CRASH and a Murnaghan Equation of State model (Material Type 28) was used to model the fluid, with pressure defined relative to current and initial density (set to 0.001 g/mm3 ) as: P =B ρ ρ0 7 −1 . This model was proposed previously by Monaghan [76] as an efficient means for modelling fluid 25 2.2. Material model investigation flow when the fluid velocity is much lower than its speed of sound propagation. A parameter, B, is set to artificially reduce the speed of sound in the fluid in order to decrease the minimum solution time step and thus increase computational efficiency. This strategy is shown to have minimal effects on the density variations and fluid behavior provided that the reduced speed of sound is maintained at least ten times the maximum flow velocity [24, 76]. To be conservative, preliminary simulations were run by reducing the value of B until the fluid simulation was no longer the limiting factor in the minimum simulation time step, yielding an optimal value for B in the current simulations of 200 MPa. This value set the speed of sound in the CSF at roughly four hundred times the maximum flow velocity observed in my simulations of ∼3 m/s. 2.2 2.2.1 Material model investigation Nonlinear regression of rat cord tensile data The approach of Goh et al. [38] was used to fit data provided by Bilston, which were previously reported [30]. Assuming separable hyper- and visco-elasticity, the stress for a hyperelastic model combined with Prony series viscoelasticity is expressed as the sum of those two components, N t σ(t) = g∞ σ0 (t) + gi e i=1 N = g∞ σ0 (t) + − t−s τ i 0 dσ0 (s) ds ds (2.4) hi (t) i=1 Briefly, Goh et al. [38] showed that by discretizing Equation 2.4, it is possible to derive an algorithm that allows solving for σ(t) for an arbitrary known applied strain history, (t). Discretizing the integral in, hi (t), Goh et al. [38] derived a recursive equation for stress, N − ∆t τi ∆t 1 − e e − τi hi (tn ) + gi σ(tn+1 ) = g∞ σ0 (tn+1 ) + [σ0 (tn+1 ) − σ0 (tn )] . (2.5) ∆t/τi i=1 Equation 2.5 is flexible in that it allows any arbitrary hyperelastic model to be used. To fit uniaxial tissue testing data to the Ogden hyperelastic model, I used the corresponding uniaxial deformation form of the stress equation11 , 1 σ0 (λ) = 2µ λα−1 − λ− 2 α−1 , (2.6) where µ is the initial shear modulus, α is a parameter describing the hyperelastic nonlinearity, Lf inal and the stretch ratio, λ, is defined as λ = 1 + eng = Linitial . I implemented this algorithm in MATLAB (see Appendix C.1) and used the nonlinear least squares curve fitting function lsqcurvefit to optimize the hyperviscoelastic Ogden and Prony parameters against a known strain history and corresponding stresses. Ogden et al. [83] previously recommended lsqcurvefit for optimizing Ogden hyperelastic model constants, and the function also performed well with my hyperviscoelastic algorithm. Validation of the curve fitting algorithm implementation was performed by plotting and recreating expected results for hyperviscoelastic models by Miller and Chinzei [74], Snedeker et al. [104] and the hyperelastic model used by Greaves [41]. 11 Note that Equation 2.6 was implemented using the PAM-CRASH definition of µP C = values in the alternative convention used by Goh et al. [38] and others. µ α compared to parameter 26 2.2. Material model investigation I also simulated the experiments by Fiford and Bilston [30], using material constants from fitting the hyperviscoelastic model to their 0.2/s strain rate to 5% peak strain test condition to simulate the other conditions. 2.2.2 Tensile coupon simulations In order to assess the hyperviscoelastic material models available in PAM-CRASH in a more controlled fashion, I created a simple tissue coupon model, similar to the shape used experimentally for uniaxial tension tests of engineering materials. I fixed the bottom coupon nodes and applied a velocity of 0.006 mm/ms to the top coupon nodes (corresponding to the 0.2/s strain rate and 30 mm gauge length of experiments by Fiford and Bilston [31]) over 250 ms to a peak of 5% applied total strain. Simulations were performed using both PAM-CRASH Material Type 37 (Viscoelastic Ogden Rubber for Solid Elements, G-Based Viscous Response) and Material Type 38 (Viscoelastic Ogden Rubber for Solid Elements, Ogden-Based Viscous Response12 ), which are both hyperviscoelastic models combining Prony series viscoelasticity with the Ogden hyperelastic model. I assigned the material parameters derived from my fit to the data from Fiford and Bilston [31]. I performed simulations using both single and double precision simulation arithmetic to detect any difference in results. Figure 2.6 shows the coupon simulation model, demonstrating a peak 5.4% local strain along the central element. I plotted the stress time histories of the central coupon element to assess and compare material model behaviour. Figure 2.6: First principal strain distribution for a tensile coupon simulation. 12 The Material 38 model was an undocumented new feature early in my project, but was added as a fully documented feature in later PAM-CRASH revisions. 27 2.3. Validation 2.3 2.3.1 Validation Initial weight-drop validation As a first step of validation and to verify whether the model was behaving similarly to that created by Maikos et al. [68], similar 12.5 mm weight-drop contusion simulations were recreated using the same boundary conditions. For these simulations a 10 g cylindrical impactor with a flat 2.5 mm diameter head was modelled and assigned a 0.489 m/s initial velocity, with a starting position at initial contact with the dura. Vertebrae were modelled as rigid and the simulation was run for 5.5 ms to capture the peak cord compression. Additional simulations were run with nonhomogeneous white and gray matter material properties to gauge the effect of inhomogeneity. The stiffness of the white matter was increased by 20% by increasing G∞ , G1 and G2 . The gross stiffness of the spinal cord was maintained by reducing that of the gray matter to compensate, according to the volume fractions of each cord component [68]. For the full medium mesh model the white matter comprised 96,000 elements for a total volume of 81 mm2 , while the gray matter comprised 73,920 elements and a volume of 53.7 mm2 , yielding volume fractions of 60.1% white and 39.9% gray matter; the appropriate calculated decrease in gray matter stiffness was 30%. Mesh size comparison Simulations were run with three spinal cord element mesh sizes, to assess convergence of the solution for decreasing element size and enable the choice of the optimal element size to balance computational efficiency against solution accuracy. The medium mesh size had element edge length of approximately 0.3 mm, while the coarse and fine mesh size elements were roughly 8 times and 1/8th the volume of the medium mesh elements. 2.3.2 Velocity and mechanism validation Further validation of the model was conducted by comparing force and displacement results to the corresponding measurements reported previously by Sparrey et al. [107] and Choo et al. [14], for each of the injury velocity and injury mechanism simulations. During attempted validation of the injury velocity experiments against force-displacement curves reported by Sparrey et al. [107], close observation revealed that displacement measurements for the 300 mm/s trials appeared to lag behind force measurements by 0.8 ms, with indenter force starting to increase from baseline before displacement had begun. I therefore shifted displacement data for the 300 mm/s experiments forward 0.8 ms and plotted experimental results against both the shifted and unshifted data. Contusion and dislocation mechanism simulation validation was attempted by comparing applied forces to the experimentally reported values [14, 15]. The spinal cord contusion force was validated with the more recent study [15] where improvements in instrumentation provided a more accurate measurement of the smaller forces measured during contusion. 2.4 Injury simulation Prior to injury velocity and mechanism simulations, I reduced the complexity of the model to include only the four vertebrae (C3-C6) located near the injury epicentre at C4/C5 (Figure 2.7). 2.4.1 Injury velocity experiments I performed simulations to model the 300 mm/s and 3 mm/s contusion injuries performed by Sparrey et al. [107] to investigate differences in finite element strain patterns according to injury velocity, which is thought to be an important factor in spinal cord injury [52, 58, 110]. 28 2.4. Injury simulation Figure 2.7: The full rat cervical model was reduced to the four-vertebra (C3-C6) segment shown here to reduce computational complexity. A rigid, flat-headed indenter was modelled after that used experimentally, and the average experimental displacement profile for each of the 300 mm/s and 3 mm/s experiments was applied to the indenter model. Contusion simulations were begun by simulating the 0.015 N dura touch force start position employed by Sparrey et al. [107] with a small 0.3 mm displacement ramp over 30 ms. 2.4.2 Injury mechanism experiments The loading and boundary conditions for the contusion and dislocation simulations were modelled to recreate the experiments by Choo et al. [14, 15] as closely as possible (Figure 2.8). All vertebrae (C3-C6) were modelled as rigid, as well as the intervertebral discs with the exception of the disc directly at injury epicentre (C4/C5). Friction between the contacting vertebrae, dura and cord was not included in the model as it has been found to have negligible influence [68]. Ends of the dura were constrained to prevent axial motion in order to encourage biofidelic membrane behavior and avoid flapping [68]. No boundary conditions were imposed on the CSF particles at either open end of the model as preliminary simulations showed minimal fluid leakage in the short time period up to peak displacement and to avoid non-biofidelic reflections at these locations from confusing results. Figure 2.8: Sagittal cross-sections of the four-vertebra spine model (C3-C6) are shown, demonstrating contusion and dislocation injury mechanism simulations during displacement. Contusion For the contusion simulation, elements from C4 and C5 vertebrae were removed to represent the partial laminectomy performed experimentally. This opening made way for a rigid indenter modelled after the 2 mm spherical headed steel indenter. The indenter was located at the dural surface and aligned normal the surface. Indenter motion was enforced by applying a velocity ramp from 29 2.5. Correlation with histology 110 cm/s (corresponding to the average experimental peak velocity just prior to impact) down to 0 over 3.2725 ms, and continuing to -110 cm/s to return to the starting position. Note that in the FE simulation, strictly defining the velocity profile over time also results in the corresponding displacement profile being strictly enforced, unlike in experimentally controlled devices where some form of feedback is guiding the control. This trajectory results in a peak indenter displacement of 1.8 mm and a peak cord compression of 1.08 mm (defined as the indenter displacement after first contact with the cord). Dislocation Prior to the dislocation simulation, the facet joints and dorsal ligaments between C4/C5 vertebrae were removed as was performed experimentally to increase injury repeatability by eliminating residual facet dislocation. The C3 and C4 vertebrae were constrained in all directions while C5 and C6 were displaced 2.5 mm dorsally from rest by applying velocity ramps up to 95.1 cm/s over 2.629 ms, down to -95.1 cm/s and back to rest. Distraction Distraction simulations were also conducted with the model to see how it would perform under that axial tension mechanism, although the model was not designed specifically with distraction simulation in mind. For these simulations, vertebrae C3 and C4 were rigidly constrained while C5 and C6 were translated caudally by 4.1 mm with a peak velocity of 110 cm/s. 2.5 Correlation with histology Regional zones were outlined in the model and spaced at 1 mm axial slice intervals from the injury epicentre to reflect the dorsal, lateral, ventromedial, and ventrolateral white matter and ventral gray matter regions in which Choo et al. [14] quantified membrane permeability to dextran (Figure 2.9,2.10). For each injury mechanism, mean membrane permeability (axons/mm, or % cells in gray matter) was plotted against mean peak values of maximum principal strain (mm/mm) for each region by matching data according to slice position (see Figure 2.11. A detailed procedure for extracting the maximum principal strain data for elements within the specified regions is given in Appendix B. For reference, all other principal strains and principal stresses were also plotted in this manner, although in general these were found to agree less with membrane permeability than did maximum principal strain (see Appendices D and E). Dorsal Lateral Ventro-medial Ventro-lateral Ventral gray Figure 2.9: Regions of the spinal cord used for quantification of strain and histology. [Illustration adapted with permission from Choo et al. [15].] Data from the white matter regional correlations were pooled to investigate whether an overall relationship existed between strain and tissue damage. Data from the gray matter region were analyzed separately as these data were quantified differently (% cells) from white matter regions [14]. Data corresponding to regions at the injury epicentre for the dislocation mechanism were excluded from all correlations because immediate hemorrhagic necrosis in this region precluded accurate quantification of dextran-positive cells and axons. 30 2.5. Correlation with histology (a) Contusion (b) Dislocation Figure 2.10: Sample strain results from spinal cord regions used for quantification, demonstrating the axial spacing of 1 mm between regional zone slices. WM_DC_DL 0mm 0.5 Maximum principal strain, ε1 Maximum principal strain, ε1 WM_DC_CT 0mm 0.4 0.3 0.2 0.1 0 0 2 4 Time (ms) 6 (a) Contusion strain histories 8 0.5 0.4 0.3 0.2 0.1 0 0 2 4 6 Time (ms) 8 10 12 (b) Dislocation strain histories Figure 2.11: Examples of finite element strain time histories. Strain histories for all elements within the dorsal column of the white matter (WM DC) in the slice at injury epicentre (0 mm) are shown for contusion (CT) and dislocation (DL). The maximum value for each element was found and the average of these maxima across all elements yielded the mean peak strain for that regional slice. 2.5.1 Statistical analysis Linear regression was used to quantify correlations between maximum principal strain in each region and for the pooled white matter data and the gray matter region. R2 correlation coefficients were calculated for each mechanism along with corresponding p-values. 31 Chapter 3 Results 3.1 3.1.1 Material model investigation Nonlinear regression of rat cord tensile data I found my hyper-viscoelastic curve fitting routine (described in Section 2.2.1) to be a useful tool in comparing and understanding the various experimental results and material models found in the literature. Specifically, my code enabled a direct comparison of the stress versus strain behaviour for different choices of spinal cord Ogden and Prony material parameters in both tension and compression, and was used to uniquely fit optimized parameters to describe experimental material testing data. Furthermore, the method allows straightforward demonstration of hyper-viscoelastic model behaviour at various strain rates and ranges. I checked for correct implementation of the hyperviscoelastic material modelling in the routine by successfully recreating stress-strain results from several previous models Greaves [41], Miller and Chinzei [74], Snedeker et al. [104]. Plots generated by using the respective material properties used in each previous study showed the same stress-strain responses shown in those studies (Figures 3.1-3.3). Plot of Bilston (1996) data fit from Greaves thesis 0.12 0.1 Stress (MPa) 0.08 0.06 0.04 0.02 0 (a) Original Ogden hyperelastic plot of spinal cord material properties [Figure reproduced with permission from Greaves [41].] 0 0.02 0.04 0.06 Strain 0.08 0.1 0.12 (b) Recreated Ogden hyperelastic plot Figure 3.1: Recreation of hyperelastic plot from Greaves [41]. I verified correct implementation of my Ogden hyperelastic code by plotting the curve generated from Greaves’ Ogden parameters and observing the same stress-strain response. Figure 3.4 compares the results from a fit to the Fiford and Bilston [31] experimental data to those corresponding to the material model used by Maikos et al. [68]. The blue line and data points show the Fiford and Bilston [31] data corresponding to a peak 5% tensile strain at a strain rate of 0.2/s. A fit to that data is shown by the dotted line, including continuation of the model beyond 5% tensile strain (stretch ratio λ = 1.05) and negative into compressive strain results (λ < 1). The nonlinear least-squares parameter fit to the data yielded a hyperelastic nonlinearity term, α, of 50, and this fit was also used to plot results for higher strain rates of 1/s and 100/s, corresponding 32 3.1. Material model investigation Plot of Miller and Chinzei (2002) data fit 2000 1000 0 Fast Tension α = −4.7 Fast Compression α = −4.7 Slow Tension α = −4.7 Slow Compression α = −4.7 Very Slow Compression α = −4.7 Fast Tension α = +4.7 Fast Compression α = +4.7 Slow Tension α = +4.7 Slow Compression α = +4.7 Very Slow Compression α = +4.7 Stress (Pa) −1000 −2000 −3000 −4000 −5000 −6000 0.6 (a) Original Ogden hyperviscoelastic plots of porcine brain material properties at three strain rates [Figure reproduced with permission from Miller and Chinzei [74]]. Experimental data are denoted with a solid line and model fits with a dashed line. Note that the Lagrange stress is defined as the applied force divided by the original cross-sectional area. 0.8 1 1.2 Stretch ratio, lambda 1.4 1.6 (b) Recreated Ogden hyperviscoelastic plots. Plots for both α = 4.7 (as used by Maikos et al. [68]) and α = −4.7 (as used by Miller and Chinzei [74], matching plots in (a)) are shown, to demonstrate different stress-strain behaviour in compression and tension for the different nonlinear parameter values (see Equation 2.6). Figure 3.2: Recreation of hyperviscoelastic plots from Miller and Chinzei [74] Plot of Snedecker (2005) data fit 30 Fast (250 /s) Slow (0.5 /s) Stress (MPa) 25 20 15 10 5 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Strain (a) Original Ogden hyperviscoelastic plots of kidney capsule material properties [Figure reproduced with permission from Snedeker et al. [104].] (b) Recreated Ogden hyperviscoelastic plots for high (250/s) and low (0.5/s) strain rates Figure 3.3: Recreation of hyperviscoelastic plots from Snedeker et al. [104] 33 3.1. Material model investigation Master stress−strain plot of experimental data and model predictions (zoom) 0.1 0.08 0.06 Stress, σ (MPa) 0.04 0.02 0 100/s − Maikos parameters 100/s − Negative Maikos parameters 100/s − Miller parameters (brain) 100/s − Fit to Fiford and Bilston data 1/s − Maikos parameters 1/s − Negative Maikos parameters 1/s − Miller parameters (brain) 1/s − Fit to Fiford and Bilston data 0.2/s − Fiford and Bilston data 0.2/s − Fit to Fiford and Bilston data −0.02 −0.04 −0.06 −0.08 −0.1 0.9 0.92 0.94 0.96 0.98 1 1.02 Stretch Ratio, λ 1.04 1.06 1.08 1.1 Figure 3.4: Stress-strain plots of experimental data and model predictions to the spinal cord contusion scenarios at impact velocities of 3 mm/s and 300 mm/s, respectively (see dashed and solid black lines). Corresponding results at the same strain rates are also plotted for the parameter set used by Maikos et al. [68] (including a value for α of 4.7), as denoted by the dashed and solid red lines. The first observation from these results is that of the change in stress-strain behaviour when increasing the strain rate – higher rates result in a relative increase in cord stiffness, with higher stress values for the corresponding stretch ratios, indicative of the viscoelastic aspect of the material models. Secondly, the material model exhibits markedly different stress-strain behaviour when using the Maikos et al. [68] parameters compared to the fit to the Fiford and Bilston [31] data, with lower stresses observed at all stretch ratios and a noted absence of the typical J-shaped hyperelastic curve in the tensile region. These differences are even more apparent when seen over a larger range of stretch ratio (Figure 3.5). Such discrepancies indicate that we still have some way to go to reconcile tissue material testing data with finite element material modelling results. Indeed, my own FE simulation attempts using the parameters obtained from fitting the Fiford and Bilston [31] data yielded an unreasonably stiff cord with contact forces orders of magnitude higher than experimental values; certainly the extrapolation of those material properties obtained for low strains (< 0.05) is not valid to the higher strains encountered in traumatic SCI mechanisms. Yet that is not to say we can’t still obtain useful results – that can enhance our understanding of dynamic cord deformation during SCI and be related to risk for tissue injury – using current models for material properties in certain strain and strain rate regimes, as demonstrated by the results of Maikos et al. [68]. 34 3.1. Material model investigation Master stress−strain plot of experimental data and model predictions 0.5 Stress, σ (MPa) 0 −0.5 100/s − Maikos parameters 100/s − Negative Maikos parameters 100/s − Miller parameters (brain) 100/s − Fit to Fiford and Bilston data 1/s − Maikos parameters 1/s − Negative Maikos parameters 1/s − Miller parameters (brain) 1/s − Fit to Fiford and Bilston data −1 −1.5 0.6 0.7 0.8 0.9 1 1.1 Stretch Ratio, λ 1.2 1.3 1.4 Figure 3.5: Stress-strain plots of experimental fit and model predictions – wide stretch ratio range 3.1.2 Tensile coupon simulations Reducing the complexity of a PAM-CRASH finite element model to the case of a simple tensile coupon sample allowed for the focused examination of material model performance. Coupon simulation results demonstrated that the Prony series viscoelasticity as implemented in Material 37 in PAM-CRASH behaved much differently than the standard implementation seen in the literature [38, 67, 68, 73, 104]. The Prony series viscoelastic stress terms for Material 37 are added on top of the underlying hyperelastic stress (and then decay in time following a deformation), while the standard implementation (as implemented in Material 38) instead subtracts the viscoelastic stress terms from the base hyperelastic stress during relaxation; this discrepancy yields different material behaviour despite using equivalent material parameters (Figure 3.6). A longer coupon simulation that was run for 1800s, to include the full viscoelastic relaxation phase of the tensile experiment, demonstrated the importance of using double precision FE simulation arithmetic to properly model the relaxation behaviour – single precision resulted in truncation of the viscoelastic material parameters, leading to different relaxation behaviour than expected (Figure 3.7). Moreover, time savings for simulations run with single versus double precision are very modest (roughly 10-20% in my comparisons), thus double precision was chosen for all further simulations. 35 3.1. Material model investigation 0.05 Stress relaxation comparison of PAM-CRASH Materials 37 and 38 0.05 0.04 Stress-strain comparison of PAM-CRASH Materials 37 and 38 0.04 Material 37 Material 38 0.03 Stress (MPa) Stress (MPa) Material 37 Material 38 0.02 0.01 0 0 0.03 0.02 0.01 400 800 1200 1600 2000 2400 0 0 0.006 Time (ms) 0.012 0.018 0.024 0.030 Strain (a) Stress-relaxation comparison of Materials 37 and 38 (b) Stress-strain behaviour comparison of Materials 37 and 38 Figure 3.6: Tensile coupon simulation comparison of PAM-CRASH hyperviscoelastic materials. Both PAM-CRASH Ogden-Prony models, Materials 37 and 38, used the same Ogden and Prony parameters and were subjected to the same strain history but exhibited markedly different stress responses. Material 38 behaves according to the standard material model implementation in the literature, while Material 37 yields much a much stiffer response that lacks the characteristic stress relaxation. 0.1 Stress relaxation comparison of single and double precision simulation Stress (MPa) 0.08 Single precision Double precision 0.06 0.04 0.02 0 0 400 800 Time (s) 1200 1600 2000 Figure 3.7: Effect of single versus double precision simulation arithmetic on long-term stress relaxation. Single precision simulation resulted in truncated viscoelastic decay terms and a corresponding cutoff of stress relaxation for long term simulation, though single and double precision results matched in the short term (< 50 s of relaxation). 36 3.2. Validation 3.2 3.2.1 Validation Initial weight-drop validation Recreating the weight drop simulations performed previously by Maikos et al. [68] was the first method of validating the UBC finite element model. This rough validation was an important first step to confirm correct gross performance of the model. This step was especially necessary due to the difference in finite element solver software used by the two models, and the corresponding differences in material property notation convention and implementation. (a) (b) Figure 3.8: Maximum principal strain distribution comparison of Maikos et al. [68] 12.5 mm weight-drop contusion results (a) with UBC model results (b). [Figure in (a) adapted from Maikos et al. [68].] Both models show a similar overall pattern with corresponding nominal values in the strain distribution, indicating that both models behave similarly despite using different FE software and model creation methods (Figure 3.8). The results are not expected to be completely the same since the model by Maikos et al. [68] is of the thoracic rat spine, while the UBC model is of the cervical rat spine. One notable difference between the results is that the strain plot from my model appears less smooth than the Maikos results, especially near the dorsal edge of the cord (at the bottom of Figure 3.8). That discontinuity in strain along the dorsal cord appears to be caused by the rough local surface of the vertebra in my model. Figure 3.9 depicts the variation of the simulated 12.5 mm weight-drop impactor trajectory for a variety of mesh sizes and model variations. The course mesh yielded a markedly different trajectory, with a lower peak displacement (1 mm) that occurred earlier (at 3.3 ms), while both the medium and fine mesh results were very similar with a peak of 1.4 mm at 3.75 ms. Furthermore, the nonhomogeneous spinal cord model with gray matter stiffness decreased 30% and white matter stiffness increased 20% had no noticeable effect on the timing or magnitude of the peak displacement, while inclusion of the cerebrospinal fluid resulted in a modest reduction of peak displacement to 1.35 mm. I later ran additional weight-drop simulations for both 12.5 mm and 25 mm drop heights using values for the Ogden hyperelastic nonlinearity parameter, α, of both 4.7 (as used by Maikos et al. [68]) and -4.7 (as proposed by Miller and Chinzei [74]). Force-displacement curves for simulations at both drop heights demonstrated higher peak forces and lower peak displacements for the “stiffer” α = −4.7 model compared to the standard α = 4.7 model (Figure 3.10). As expected, the higher impact velocity 25 mm weight-drop simulations follow the 12.5 mm force-displacement curves closely but are seen to continue on to correspondingly higher peaks. Note that the peak impactor displacement for the standard (α = 4.7) 12.5 mm weight-drop shown here (1.7 mm) was 37 3.2. Validation (a) Weight-drop impactor trajectories for varying mesh sizes (b) Strain distributions for several model variations: (clockwise from top-left) fine mesh cord and dura, fine mesh non-homogeneous spinal cord, medium mesh including cerebrospinal fluid (other variations shown omitted the CSF), and medium mesh. Figure 3.9: Results of 12.5 mm weight drop simulations for varying model mesh sizes and other model variations. higher than that from earlier results (1.4 mm). This difference is due to the switch from using PAM-CRASH Material Type 37 (Viscoelastic Ogden Rubber for Solid Elements, G-Based Viscous Response) in the initial weight-drop simulations to using Material Type 38 (Viscoelastic Ogden Rubber for Solid Elements, Ogden-Based Viscous Response) for later simulations (see previous discussion in Sections 2.2.2 and 3.1.2). 38 3.2. Validation 7 Force-displacement curves for weight-drop simulations 12.5mm negMaikos 12.5mm Maikos 25mm negMaikos 25mm Maikos 6 Force (N) 5 4 3 2 1 0 0 0.3 0.6 0.9 1.2 Displacement (mm) 1.5 1.8 Figure 3.10: Force-displacement curves for 12.5 mm and 25 mm weight-drop simulations with α = ±4.7. “Maikos” denotes α = 4.7, while “negMaikos” denotes α = −4.7. 3.2.2 Velocity and mechanism simulations Injury velocity Figure 3.11 shows the force-displacement curves for the simulations of fast and slow contusions, alongside the experimental data before and after correcting for a time lag in the experimental displacement data. Blue lines depict the results for linear elastic spinal cord material properties, as expected demonstrating no difference between the 3 mm/s (solid) and 300 mm/s (dashed) injury velocities and a constant slope throughout the displacement. In contrast, the black lines corresponding to the hyperviscoelastic properties show a hyperelastic response with increasing slope up to peak displacement, and some evidence of viscoelasticity due to the 300 mm/s (dashed) curve falling above the slower 3 mm/s (solid) result. However, it is obvious that the viscoelastic effect modelled in the simulation is much smaller than that observed experimentally. Furthermore, it is interesting to note that the hyperelastic curves predicted by the model exhibit much different slopes compared to the relatively linear results from the experiments by Sparrey et al. [107]. Due to the clear differences between the experimental and simulation results for this injury velocity study, the model cannot be argued to be valid for these conditions, thus precluding further analysis and comparison of the simulated strain distributions at these two contusion velocities. Injury mechanism I validated the finite element model for contusion and dislocation mechanisms by comparing predicted loads from the enforced displacement profiles to the corresponding experimental corridors from the Choo et al. [14] data (see Figure 3.12). The simulated contusion force, applied directly to the dura and cord, followed a very similar pattern to the mean experimental force but was delayed relative to the experimental traces (Figure 3.12a). The peak force for the contusion simulation was 1.4 N, which was close to the experimental mean of 1.5±0.4 N (± SD). Forces applied during dislocation were much higher than direct cord contusion forces as they were instead applied to the entire vertebral column. The simulated dislocation force demonstrated a 39 3.2. Validation Force−displacement curves (with displacement shift) Force−displacement curves (without displacement shift) 3.5 3 2.5 2.5 2 1.5 2 1.5 1 1 0.5 0.5 0 0 −0.5 0 0.2 Mean force−disp. − 300 mm/s Mean +/− 2*S.D. − 300 mm/s Simulation − 300 mm/s Simulation − 300 mm/s LE Simulation − 3 mm/s Simulation − 3 mm/s LE 3 Force (N) Force (N) 3.5 Mean force−disp. − 300 mm/s Mean +/− 2*S.D. − 300 mm/s Simulation − 300 mm/s Simulation − 300 mm/s LE Simulation − 3 mm/s Simulation − 3 mm/s LE 0.4 0.6 0.8 Displacement (mm) 1 −0.5 0 1.2 0.2 0.4 0.6 0.8 Displacement (mm) (a) 1 1.2 (b) Figure 3.11: Force displacement curves for injury velocity simulations. Experimental data from Sparrey et al. [107] marked with circles and shown (a) before and (b) after shifting to correct for a time lag in the 300 mm/s experimental displacement data. time lag that was not observed as strongly experimentally (Figure 3.12b), but the force-time curve followed a similar path to the experimental corridor. The peak dislocation force of 17.6 N was within the experimental range but below the mean of 24.7±5.7 N. In addition, both experimental and simulated dislocation force traces demonstrated multiple local peaks that indicate sequential failure of local soft tissue components such as the intervertebral disc and spinal ligaments. Dislocation force time histories Contusion force time histories 40 0.5 0 30 Force (N) Force (N) −0.5 −1 20 10 −1.5 Experimental forces Mean experimental force Mean force +/− 2*S.D. Simulated force −2 −2.5 0 5 10 15 20 Time (ms) (a) 25 30 0 −10 0 2 4 6 Time (ms) 8 10 (b) Figure 3.12: Simulated applied forces for contusion and dislocation injuries are plotted alongside corresponding experimental traces and corridors. (a) The simulated contusion force, applied directly to the dura and cord, follows a very similar pattern to the mean experimental force but is seen to lag behind. (b) While the simulated dislocation force deviates substantially from the experimental corridor, the peak force is within the experimental range and both experimental and simulated traces demonstrate multiple local peaks that indicate sequential failure of the intervertebral disc and spinal ligaments. Results from distraction simulation were far off from experimental results, indicating poor biofidelity for simulating this mechanism with the model in its current form. The peak simulated 40 3.3. Internal strain distributions distraction force was 29.1 N, below the experimental mean of 37.9±5.4 N. Preliminary comparison of regional strain results to histological tissue damage distribution also demonstrated poor biofidelity of the distraction simulations, with none of the experimental features evident and showing tissue damage decreasing with increasing strain, therefore further analysis of the distraction results was omitted. 3.3 Internal strain distributions Mid-sagittal images of the deformed spinal cords with internal strain distributions highlight differences between the contusion and dislocation mechanisms (see Figure 3.13). Some of the most striking features were the high strain (>0.16) ‘tails’ seen dorsocaudal and rostroventral to the dislocation injury epicentre (see Figure 3.13b), showing the local regions subjected to tension due to locally isolated dynamic rotation of the cord at the epicentre caused by opposing C4/5 vertebral motion. These rostrocaudal asymmetries about the epicentre were evident in the affected dorsal, ventromedial, and ventrolateral white matter region strain plots (Figure 3.15f, 3.15h, 3.15i). Tissue damage in both ventral white matter regions also showed increases rostral compared to caudal to the epicentre (Figure 3.15c, 3.15d), though the caudal increase in the dorsal white matter was not evident in the experimental results (Figure 3.15a). (a) Contusion (b) Dislocation Figure 3.13: Distribution of maximum principal strain during contusion (a) and dislocation (b) injury simulations. Sagittal slices shown at midline and ±1 mm. Contusion plots show a relatively small very high strain zone (>0.4 maximum principal strain) localized to the dorsal surface and extending less than 1 mm lateral to the midline, while the very high strain zone for dislocation zone encompasses the full depth of the cord and extends beyond 1 mm laterally, demonstrating differences in mechanical injury severity and extent between the two mechanisms. Other notable features are the high strain (>0.16) ‘tails’ observed dorsocaudal and rostroventral to the dislocation injury epicentre, showing the local regions subjected to tension due to local dynamic rotation of the cord at the epicentre. 3.4 Correlation with histology Strong overall correlations between maximum principal strain and tissue damage indicated the model’s biofidelity and corresponding utility. For the pooled white matter regions, maximum principal strain showed significant (p < 0.0001) correlations with tissue damage for both contusion 41 3.4. Correlation with histology 70 70 60 60 Permeable Axons/mm Permeable Axons/mm (R2 = 0.86, Figure 3.14a) and dislocation (R2 = 0.54, Figure 3.14b) mechanisms, with more damage with increasing strain (see Table 3.1). For the ventral gray matter, maximum principal strain correlated strongly with tissue damage for both contusion and dislocation mechanisms (both R2 = 0.93, Figures 3.14c and 3.14d). Maximum principal strain distributions (see Figure 3.15) were similar in nature to the distributions of membrane permeability for contusion and dislocation mechanisms, with central peaks flanking the injury epicentre and decreasing tails toward the caudal and rostral extremes. All non-pooled regions yielded high correlation coefficients for contusion (R2 between 0.90 and 0.96), as did four of the five regions for dislocation (R2 between 0.61 and 0.96), with a notable exception being the dorsal white matter (R2 = 0.38). Interestingly, in the dorsal white matter the minimum, or third, principal strain showed much higher correlation with tissue damage (R2 = 0.83) for dislocation than did maximum principal strain (see Figure E.2). 50 40 30 20 10 0 0 50 40 30 20 10 0.1 0.2 ε1 0.3 0 0 0.4 70 70 60 60 50 50 40 30 10 10 0.2 ε1 0.4 0.3 0.4 30 20 0.1 0.3 40 20 0 0 0.2 ε1 (b) Dislocation WM % of cells % of cells (a) Contusion WM 0.1 0.3 0.4 0 0 (c) Contusion GM 0.1 0.2 ε1 (d) Dislocation GM Figure 3.14: Correlations of maximum principal strain with tissue damage for contusion and dislocation mechanisms in the white and gray matter. The complete results are presented in Appendices D and E. 42 3.5. Regional distribution of strain and tissue damage Table 3.1: Correlation coefficients for maximum principal strain and tissue damage within cord regions Contusion Cord region Ventral gray Dorsal white Lateral white Ventromedial white Ventrolateral white Pooled white Dislocation Correlation coefficient, R2 p-value Correlation coefficient, R2 p-value 0.93 0.94 0.96 0.90 0.91 0.86 1.5E-06 5.9E-07 1.1E-07 8.0E-06 5.6E-06 2.1E-19 0.93 0.38 0.96a 0.82 0.61 0.54b 6.3E-06 0.056 9.5E-7 3.1E-04 0.0075 6.8E-08 a This value is improved substantially from that presented previously (R2 = 0.38 and p-value = 0.056) by Russell et al. [95], due to correction of the data point at +1 mm for dislocation in the lateral white matter after discovery and correction of a data exporting error affecting that point. b This value is improved slightly from that presented previously (R2 = 0.52 and p-value = 1.8E-07) by Russell et al. [95], due to correction of the data point at +1 mm for dislocation in the lateral white matter. 3.5 Regional distribution of strain and tissue damage Key differences between the contusion and dislocation mechanisms lie in the size and shape of the central zones of very high maximum principal strain (>0.4), with the contusion zone located only near the dorsal surface and extending less than 1 mm lateral to the midline, while the dislocation zone encompassed the full depth of the cord and extended beyond 1 mm laterally. This corresponded to much higher average peak strains at epicentre for dislocation than contusion in the lateral white matter (Figure 3.15g13 ) and ventral gray matter (Figure 3.15j), and was reflected in the experimental dislocation results (Figure 3.15b, 3.15e) by especially low epicentre counts due to primary axotomy and widespread necrosis in these areas respectively. 13 The plot in Figure 3.15g differs slightly from that presented previously by Russell et al. [95], due to correction of the data point at +1 mm for dislocation (red) in the lateral white matter after discovery and correction of a data exporting error affecting that point. 43 Dorsal white matter 70 Lateral white matter Contusion Dislocation 60 Ventromedial white matter Ventrolateral white matter 70 70 70 60 60 60 50 50 50 Ventral gray matter 100 30 40 30 % of cells 30 40 Axons/mm Axons/mm Axons/mm 40 Axons/mm 80 50 40 30 20 20 20 20 10 10 10 10 60 40 20 4 2 0 −2 −4 Rostral distance from epicenter (mm) 0 6 −6 (a) 0 6 −6 (b) Dorsal white matter 4 2 0 −2 −4 Rostral distance from epicenter (mm) 0 6 −6 (c) Lateral white matter 4 2 0 −2 −4 Rostral distance from epicenter (mm) 0 6 −6 (d) Ventromedial white matter 0.4 0.4 0.4 0.4 0.4 0.3 ε1 0.5 ε1 0.5 ε1 0.5 0.3 0.3 0.3 0.2 0.2 0.2 0.2 0.2 0.1 0.1 0.1 0.1 0.1 0 6 4 2 0 −2 −4 Rostral distance from epicenter (mm) (f ) −6 0 6 4 2 0 −2 −4 Rostral distance from epicenter (mm) (g) −6 0 6 4 2 0 −2 −4 Rostral distance from epicenter (mm) (h) −6 0 6 4 2 0 −2 −4 Rostral distance from epicenter (mm) (i) −6 Ventral gray matter 0.5 0.3 4 2 0 −2 −4 Rostral distance from epicenter (mm) (e) Ventrolateral white matter 0.5 ε1 ε1 4 2 0 −2 −4 Rostral distance from epicenter (mm) −6 0 6 4 2 0 −2 −4 Rostral distance from epicenter (mm) −6 (j) Figure 3.15: Rostrocaudal distributions of experimentally measured tissue damage and computed maximum principal strain are plotted at 1 mm intervals for both contusion (blue, solid) and dislocation (red, broken). (a-e) Tissue damage as measured by counts of cells permeable to dextran is generally higher for dislocation than contusion, with peaks near injury epicentre trailing down to lower levels rostral and caudal. Local dips directly at injury epicentre for dislocation injuries are the result of immediate hemorrhagic necrosis (these points were excluded from statistical analyses and correlations). (f-j) The average peak maximum principal strains for each cord region show highest strains near injury epicentre, as expected, and higher peak strains for dislocation than contusion. Much higher average peak strains at epicentre for dislocation than contusion in the lateral white matter (g) and ventral gray matter (j) regions reflect the larger lateral extent of the very high strain zone for dislocation as seen in Figure 3.13, and are supported by especially low epicentre cell counts due to rampant necrosis in these areas (b,e). In addition, rostrocaudal asymmetries about the dislocation epicentre corresponding to the high strain tails from Figure 5 are seen in the affected dorsal (f), ventromedial (h), and ventrolateral (i) white matter strain plots, matching similar patterns of tissue damage in the ventral white matter (c,d), though the caudal increase in the dorsal white matter (f) is not evident in the experimental results (a). Error bars denote standard error for experimental cell permeability counts (a-e), and standard deviation of peak strains for all elements in the simulation regions (f-j). 3.5. Regional distribution of strain and tissue damage 0 6 44 Chapter 4 Discussion In the following sections I highlight notable results and relate them to relevant literature, and then identify and address the inherent limitations of my work. 4.1 Injury simulation I developed a finite element model of contusion and dislocation injuries in the rat cervical spine to address the lack of knowledge regarding the effect of injury mechanism on primary mechanical damage patterns in the spinal cord. It is the first multi-mechanism computational model to be based on experimental injury mechanisms in an animal model, and thus has the advantages of more direct validation and comparison with histological tissue damage. In addition, this work demonstrated the feasibility of using the Smoothed Particle Hydrodynamics method to model the cerebrospinal fluid during impact, which may be useful in large animal or human SCI FE models (that involve larger subdural spaces and CSF volumes compared to the rat) in which the CSF has been shown to play an important role to cushion impact to the cord [55, 86]. Overall, the model demonstrated its versatility to simulate both contusion and dislocation injury mechanisms with good biofidelity. The hyperviscoelastic material properties of the spinal cord yielded a realistic contact force during contusion, and the material properties of the intervertebral disc and spinal ligaments (including failure strain limits) resulted in dislocation force profiles that were similar to those measured in animal models. Note that the applied dislocation forces (experimental peak range ∼17-38 N) are much higher than direct cord contusion forces (experimental peak range ∼0.75-2.1 N) as they are instead applied to the entire vertebral column. The simulated force history for dislocation, however, deviated somewhat from the experimental corridors, including a prolonged toe region of increasing stiffness at the start of displacement that indicates the involved disc and ligaments were not behaving stiffly enough in this initial phase – possibly due to inaccurate pre-tension in ligaments, omission of muscle attachments, or inability of the linear elastic material properties of the disc and ligaments to model behavior accurately in this regime. Because of this, the current ligamentous cervical spine model cannot be considered fully validated and should not be used without further refinement to model more general, external perturbations to the spine such as rear or head first impacts. This limitation on the biofidelity of gross spinal column forces during dislocation, though, did not affect the time course or amount of spinal cord deformation in our study as this was determined by contact with displacementcontrolled vertebrae. Indeed, maximum principal strain was shown to correlate well to tissue damage for both contusion and dislocation cervical injury mechanisms, extending previous findings for thoracic weight-drop contusion by Maikos et al. [68]. 4.1.1 Interpretation of cord strain patterns and correlations The correlations between maximum principal strain and tissue damage presented here did not suggest any obvious damage thresholds – such as a minimum strain level below which no damage was observed – that could be used as an injury criterion. However, all regions subjected to at least 0.1 maximum principal strain corresponded to elevated average levels of tissue damage, while averages for regions less than 0.1 strain varied between baseline and moderate levels of damage. The variation below 0.1 strain was especially high within the dislocation results, possibly due to 45 4.1. Injury simulation less repeatability for this mechanism. This loose “threshold” of 10% strain is slightly below lower bounds on tissue damage thresholds of 13-19% found previously [4, 103]. Furthermore, steeper slopes for the ventral gray matter correlations compared to the white matter are suggestive of a lower injury tolerance to tensile strain in the gray matter, though it should be noted that the correlations for gray matter are less conclusive due to the small number of samples in that region in the current study; future studies with a focus on better quantifying such differences in injury tolerance are certainly warranted. Correlations of maximum principal strain with tissue damage were high in all regions studied for contusion and in most regions for dislocation, with the exception of poor correlation in the dorsal white matter. In fact, within the dorsal white matter the minimum principal strain showed much higher correlation with tissue damage that maximum principal strain; this may indicate a compressive tissue failure mechanism may play an important role in the dorsal white matter during dislocation. Interestingly, although the computational results for contusion injuries show peak strains beneath the tip, histology following contusion shows damage focused in the gray matter. This discrepancy is due to the lower injury threshold of the highly vascularized gray matter in comparison to white matter, as found previously by Maikos and Shreiber [66] in weight drop contusion injuries. At high enough strains, one would expect to see primary damage in the white matter. Indeed this is the case. The models predict greater strain in the lateral white matter during dislocation injury (>0.4 in Figure 3.15g) compared to the dorsal white matter during contusion (approximately 0.3 in Figure 3.15f). Accordingly, the dorsal white matter is often spared during contusion whereas primary damage of the lateral white matter is common in our animal models of fracture-dislocation (see loss of axons at epicentre in Figure 3.15b). While I did find a correlation of maximum principal strain with the tissue damage results of Choo et al. [15], the pattern of strain from my cervical simulations is quite different from the pattern of tissue damage found for thoracolumbar dislocation results by Clarke et al. [18] (Figure 4.1). Unlike the cervical dislocation experiments, the anterior thoracolumbar dislocations were performed at a slower velocity (peak of 0.22 m/s compared to 1 m/s) and with an additional vertebra in between the fixed and displaced vertebrae. Elevated tails of maximum principal strain (green) are observed in the white matter opposite to vertebral contact, while the thoracolumbar dislocation experiments demonstrated tails of tissue damage in the white matter closest to the contacting vertebrae. This may indicate important differences in local cord deformation between the two dislocation injury mechanisms, likely due to the absence of the localized vertebral pinching (as present between the fixed C4 and displaced C5 in the dislocations by Choo et al. [14] which I simulated) in the more distributed and slower dislocations performed by Clarke et al. [18]. In my simulations, this vertebral pinching was observed to tightly constrain the cord at injury epicentre, resulting in rotation of the cord at epicentre as the nearby vertebral surfaces slid past each other. This dynamic cord rotation was resisted by inertia along the cord length, resulting in tension along the surfaces of the white matter and the tails of elevated maximum principal strain. The distributed thoracolumbar dislocation, which includes an additional vertebra between the fixed and displaced vertebrae, did not produce this localized pinching mechanism. This difference highlights the wide range of cord deformation and corresponding tissue damage possible for different injury mechanisms, even between two examples of anterior dislocation. In addition to correlating strain and tissue damage and investigating injury tolerances, the model is also a useful tool to compare injury severity between mechanisms. For the mechanisms I studied, which were developed previously by Choo et al. [14], dislocation appears much more severe in both peak maximum principal strain intensity and extent. This bears some similarity to the results of Greaves et al. [40], whose quasi-static human model showed deeper and wider extent of von Mises strain for dislocation than contusion. Such comparisons and predictions of injury mechanism severity will be useful for further development of consistent and well characterized injury 46 4.1. Injury simulation (a) Simulated maximum principal strain pattern (b) Experimental dislocation tissue damage pattern Figure 4.1: Pattern of maximum principal strain in cervical dislocation is different from experimental tissue damage in distributed thoracolumbar dislocation. A sagittal slice showing the pattern of maximum principal strain (a) is compared to a diagrammatic depiction of the tissue damage pattern (dark shading) observed by Clarke et al. [18] (b). Elevated tails of maximum principal strain (green) are observed in the white matter opposite to vertebral contact, while the thoracolumbar dislocation experiments demonstrated tails of tissue damage in the white matter closest to the contacting vertebrae. This may indicate important differences in local cord deformation between the two dislocation injury mechanisms. In my simulations, vertebral pinching was observed to tightly constrain the cord at injury epicentre, resulting in rotation of the cord at epicentre as the nearby vertebral surfaces slid past each other. This dynamic cord rotation was resisted by inertia along the cord length, resulting in tension along the surfaces of the white matter and the tails of elevated maximum principal strain. The distributed dislocation, which includes an additional vertebra between the fixed (T12) and displaced (L1) vertebrae, did not include this localized pinching mechanism. [Diagram in (b) adapted from Clarke et al. [18].] protocols and, alongside behavioral survival studies, can increase our understanding of differences in functional deficits and treatment goals between mechanisms. 4.1.2 Force history delay and subarachnoid space An interesting observation from contusion simulation results was a 1.8 ms delay in simulation force peaks relative to experimental results. This may be due to the different subarachnoid space between the dura and cord (filled with cerebrospinal fluid) at cervical and thoracic levels. Maikos et al. [68]’s model used a thickness of ∼50-80 µm, while Franconi et al. [33] show the thickness at ∼300 µm, which is what my model dura was based on. Note that after specifying the offset for the dura, its size had to be reduced in some places to conform to the spinal column geometry and avoid penetrating vertebrae or vertebral discs. Another phenomenon related to the timing of peak contusion force, noticed during close examination of the experimental results of Choo et al. [15] and recent experiments by Bhatnagar (personal correspondence, 2010), was that the peak force occurred noticeably before the peak cord compression. This effect may be due to very short term viscoelastic relaxation during the dynamic cord compression, or perhaps due to some tissue damage or failure occurring prior to peak compression. A preliminary investigation to attempt to more accurately simulate this effect by reducing the viscoelastic relaxation time constant yielded modest improvements, but could not explain the majority of this viscoelastic effect. If relaxation is the cause, it is likely that a fully nonlinear viscoelastic model may be necessary to capture this nuance of spinal cord behaviour during dynamic deformation [114]. 47 4.2. Material model investigation 4.2 Material model investigation By plotting experimental spinal cord stress-strain results alongside results corresponding to several hyperviscoelastic models, using the method of Goh et al. [38], I was able to elucidate large differences in material behaviour. In particular, these plots demonstrated that the cord material parameters used in both the model by Maikos et al. [68] and in my own model exhibit a much more linear response in tension than the highly hyperelastic experimental response recorded during spinal cord tensile testing by Fiford and Bilston [30]. I found, however, that trial simulations using the material parameters from the experimental data fit resulted in a much too stiff cord with contact forces far exceeding those measured experimentally – this is likely the reason that the material parameters calibrated by Maikos et al. [68] to match their experimental weight-drop results deviated substantially from those fit to the Fiford and Bilston [30] data. The strain plots also revealed the different behaviour corresponding to the differing hyperelastic nonlinearity parameter values of α = 4.7 (as used by Maikos et al. [68]) and -4.7 (as proposed by Miller and Chinzei [74]), with behaviour for the latter falling somewhere between that for α = 4.7 and the fit to the Fiford and Bilston [30] data. Finally, the full set of hyperviscoelastic material parameters proposed by Miller and Chinzei [74], as derived from their testing of porcine brain tissue samples at relatively slow strain rates (0.0064 and 0.64 s−1 ), yielded stress-strain curves far less stiff than each of the other plotted models. A key message from these results is that further spinal cord dynamic tissue testing and complementary simulation work are needed to bridge the current gap between the methodologies and allow for a fully biofidelic spinal cord material model for all strain rates and magnitudes. I also used the hyperviscoelastic fitting algorithm to attempt fitting all of the data from Fiford and Bilston [30] for several peak strains and strain rates. Despite allowing for viscoelastic relaxation during the loading ramp which is typically neglected in quasi-linear viscoelastic (QLV) modelling, this investigation yielded the same conclusion as found by Fiford and Bilston [30], namely that a linear viscoelastic model (ie. one with strain-independent viscoelasticity, in this case via a Prony series) cannot explain the different steady-state stress levels reached after relaxation from fixed strain levels reached at varying strain rates. Unfortunately, while the algorithm I employed is ideal for determining material constants to correspond to current finite element software hyperviscoelastic material models that employ Prony series, it could not be used for more general nonlinear viscoelastic models that may be employed in the future, such as the recent fully nonlinear viscoelastic modelling method proposed by Troyer and Puttlitz [114]. Extending my material model investigation from theoretical modelling to explicit finite element simulation, via a simplified tissue coupon model, allowed for further exploration of the precise behaviour of the hyperviscoelastic material models available in PAM-CRASH. Tensile tissue coupon simulations highlighted the nature of the different viscoelastic implementations used in each of Material Type 37 (Viscoelastic Ogden Rubber for Solid Elements, G-Based Viscous Response) and Material Type 38 (Viscoelastic Ogden Rubber for Solid Elements, Ogden-Based Viscous Response), allowing me to confidently choose Material 38 for injury simulations14 to be consistent with the theoretical hyperviscoelastic models described in the literature [38, 67, 68, 73, 104]. The importance of using double precision arithmetic during FE simulation was also evidenced by the tissue coupon viscoelastic relaxation results, and is important to note for further simulations. 14 12.5 mm weight-drop contusion simulations further demonstrated the magnitude of differences in injury simulation results that could be expected for the two different material implementations, with Material 38 yielding a peak impactor displacement of 1.7 mm compared to 1.4 mm for Material 37. 48 4.3. Limitations 4.3 Limitations Several limitations of the current work suggest possible improvements for the future and are discussed below. 4.3.1 Spinal cord material properties The hyperviscoelastic cord properties proposed by Maikos et al. [68] – and further validated by our study – model spinal cord behavior quite well but were based initially on material testing data and then adjusted to better match experimental behavior. Further material testing and modelling of the rat cord is necessary, including investigation of white and gray matter inhomogeneity and anisotropy of the tissue to determine the importance of such factors in modelling. White and gray matter inhomogeneity Furthermore, there is no current consensus on possible differences in white and gray matter material properties [3, 20, 50, 68]. I omitted such differences due to the ongoing uncertainty in their exact numerical values. Relative differences in the material properties of white and gray matter would cause the two components to deform by different magnitudes, with greater deformation in the softer material than the stiffer material; this difference in deformation would result in an additional shearing strain at the interface between the gray and white matter [106]. Viscoelastic properties In particular, the viscoelastic properties – extracted by Maikos et al. [68] from dynamic brain tissue tests by Mendis et al. [70] – seem to yield good results in simulating cord behaviour in high velocity impacts (with peak velocities of 0.489-0.690 m/s in Maikos et al.’s weight drop contusions and 1 m/s in the current study), but a more detailed characterization of rat cord viscoelasticity may be required to model its behavior well over a wider range of impact velocities. In fact, the Ogden and Prony hyperviscoelastic model we employed to simulate the cord properties may ultimately prove inadequate to accurately model cord behavior in all desired scenarios, and development of more complicated material models may be necessary that include fully nonlinear viscoelasticity and a mechanism for dynamic tissue failure during simulation. Because of the difficulty of accurately testing soft tissues like the cord to determine material properties, some combination of material testing and FE simulation to match test behavior may be required to achieve further characterization, as proposed by Morriss et al. [77]. Further testing considerations should also include careful attention to the influence of preconditioning on mechanical test results of the cord at various strain rates and magnitudes, an issue recently highlighted by Cheng et al. [13]. 4.3.2 Strain direction My analysis of the maximum principal strain during spinal cord injury simulation focused on the strain distribution patterns and quantification of average peak strains within regions of the cord. This analysis did not investigate principal strain directions associated with specific points within the cord, which could be ascertained from plotting the direction vectors from FE results for select elements in regions of interest. Such an approach has not been taken in other FE models of SCI, and, while outside the scope of my thesis, may yield novel information on the complex deformations of the cord. Detailed analysis of strain direction could be especially useful for investigating localized damage patterns for specific anatomic structures within the cord, considering anisotropy of those structures when interpreting dominant strain directions in the area. 49 4.3. Limitations 4.3.3 Cerebrospinal fluid validation The Smoothed Particle Hydrodynamics approach I used to model the cerebrospinal fluid is novel, and qualitatively seems to allow a much more realistic simulation of fluid behaviour in the spine than for methods used previously Kleiven [60], Maikos et al. [68], Scifert et al. [98]. However, specific validation of the behaviour of the CSF SPH model has not yet been performed. This is partly due to a current limitation of the SPH method in PAM-CRASH, namely that the solver does not yield smooth pressure distribution results15 , thus precluding accurate comparisons against experimental measurements such as those performed recently by Jones et al. [56] using a porcine SCI model. Further, since the rat spine contains a relatively thin layer of CSF 16 , and since the presence of the CSF seemed to have a modest effect on the simulation results as demonstrated in Figure 3.9, specific validation of the CSF behaviour itself was not deemed central or necessary for the current work. 15 ESI have indicated that they have a development version of the PAM-CRASH solver that has addressed the pressure distribution issue, and that they should be able to provide us with such a version in the future. 16 The thickness of the CSF layer in my cervical rat model was an average of 0.4 mm compared with the corresponding human value of 3.35 mm, and bovine of 1.5 mm [86, 87]. The thickness of the CSF layer in the porcine spine is similar to that in the human [56]. 50 Chapter 5 Conclusion In this final chapter I state the conclusions of my project with respect to my objectives and summarize the contributions my thesis project has made to the field. I then provide recommendations for future related work and finish with a reflection on the relevance of my thesis to the overall goals of better prevention and treatment of spinal cord injuries. 5.1 Conclusions The main results for each of my research objectives were: 1. Objective: Create a biofidelic dynamic, nonlinear finite element model of the rat cervical spine. My first objective was accomplished with the creation of a novel FE model of the rat cervical spine based on geometry obtained from magnetic resonance imaging. I based material properties for model components on values from previous experimental or modelling literature to achieve model biofidelity wherever possible. In particular, my model incorporated the spinal cord white and gray matter, dura mater, cerebrospinal fluid, spinal ligaments, intervertebral discs, a rigid indenter and vertebrae, and failure criteria for ligaments and vertebral endplates. I used a hyperviscoelastic Ogden-Prony material model for the spinal cord based on that used in a similar model of the thoracic rat spine by Maikos et al. [68] – the only rat cord material property calibration at similar injury velocities. Model creation was a necessary stepping stone for my other objectives, with specific biofidelity of the model for various applications assessed during the validation stage. 2. Objective: Simulate spinal cord contusion experiments at impact velocities of different orders of magnitude and compare to experimental results. I simulated the 3 mm/s and 300 mm/s contusion injuries performed experimentally by Sparrey et al. [107], using both hyperviscoelastic spinal cord properties and linear elastic cord properties for reference. I used a medium finite element mesh size (edge length ∼0.3 mm) for these simulations which had long durations of up to one third of a second for the 300 mm/s velocity. Results showed poor correspondence with experimental behaviour (see Validation below), so further analysis was not performed. 3. Objective: Simulate dynamic spinal cord injury experiments for contusion and dislocation injury mechanisms and compare FE strains to tissue damage. I simulated high speed (1 m/s) contusion and dislocation injuries between vertebral levels C3 and C6 to match previous experiments by Choo et al. [14]. For these short duration (on the order of several milliseconds) simulations I used a fine mesh size (edge length ∼0.15 mm) to obtain more detailed strain distribution results. I calculated average peak maximum principal strains for several cord regions at the injury epicentre and at 1 mm intervals from +5 mm rostral to -5 mm caudal to the lesion, and compared peak strains to tissue damage measured previously via axonal permeability to 10 kD fluorescein-dextran. Linear regression of tissue damage against peak maximum principal strain for pooled data within white matter regions showed significant (p < 0.0001) correlations that are similar for both contusion (R2 = 0.86) and dislocation (R2 = 0.54). 51 5.2. Contributions 4. Objective: Validate the FE model by comparing computed injury forces to experimentally measured values. I performed a first step of model validation by recreating a weight-drop simulation from Maikos et al. [68] and demonstrated similar strain patterns and magnitudes to their results. In this step I also assessed model performance with different finite element mesh sizes, rejecting the coarse mesh size (edge length ∼0.6 mm) results that deviated from the similar results for both medium (∼0.3 mm) and fine (∼0.15 mm) mesh sizes. I attempted validation of the contusion injury velocity simulations by comparing resulting force-displacement curves with the experimental curves. While the 300 mm/s injury velocity simulation results showed a slightly stiffer force-displacement response than for 3 mm/s, the difference was much smaller than that observed experimentally and the model cannot be considered validated. This indicates that the current material properties do not accurately model cord deformation for injury velocities in this 3-300 mm/s range. Additional contusion velocity simulations using linear elastic cord properties deviated even further from experimental results, with no difference between 300 mm/s and 3 mm/s simulations, highlighting the need for hyperviscoelastic properties to model cord behaviour over a range of injury velocities. I validated contusion and dislocation injury mechanisms by comparing force time histories with experimental corridors. The simulated peak contusion force of 1.4 N was very close to the experimental mean of 1.5±0.4 N (± SD), while the simulated dislocation peak of 17.6 N was within the experimental range but below the mean of 24.7±5.7 N. Contusion force followed history closely resembled the experimental traces, indicating appropriate contact forces from the cord material model during deformation. While both experimental and simulated dislocation force traces demonstrated multiple local peaks that indicate sequential failure of local soft tissue components (ie. the intervertebral disc and spinal ligaments), the simulation demonstrated lower stiffness during the initial displacement and a more abrupt failure than the experiments; this indicates that the material properties of the ligaments and intervertebral disc should be refined before my model can be used to simulate external perturbations (that do not prescribe vertebral displacement directly) such as impact to the spinal column during rear or head-first impact. Additional simulation of a distraction injury mechanism could not be validated, indicating that substantial refinements to the model should be made to accurately simulate distraction. 5.2 Contributions Various aspects of my thesis project yield novel contributions. Foremost of these is the fact that this represents the first application of a dynamic model of spinal cord injury to multiple injury mechanisms. I show that maximum principal strain results for contusion and dislocation injury mechanisms correlate well with axonal damage. I further find that distinct strain distribution patterns help distinguish damaging cord deformations for the two mechanisms. My model also enabled the first finite element attempt at investigating the effect of injury velocity on contusion SCI, by simulating contusion velocities over several orders of magnitude (3 mm/s to 300 mm/s). While the results were inconclusive due to lack of validation against experimental results, my investigation did highlight the limitations of current hyperviscoelastic material models. In particular, this work demonstrates that the viscoelastic material properties currently used do not model cord behaviour well over such a wide range of strain rates. The use of Smoothed Particle Hydrodynamics (SPH) in my model to simulate cerebrospinal fluid flow during SCI is also a novel approach. SPH is a recent addition to dynamic finite element software and allows for efficient and straightforward fluid-structure interaction within the model. 52 5.3. Recommendations for future work A further contribution arises from my development of code to plot stress-strain behaviour of hyperviscoelastic material models and to optimize parameters to fit experimental data. I use this code to compare and contrast the wide range of behaviour predicted by different material parameters found in the literature. Comparison of model behaviour is difficult without the use of such analytical code, and there is no similar comparison in the literature to date. 5.3 Recommendations for future work The following are some recommendations I make for future work related to my thesis project: ❼ Specific refinement of the FE model to enable biofidelic simulation of distraction. Notably, nerve roots (which serve to locally anchor spinal cord segments and were omitted from the model to date) should be added to better simulate the distraction mechanism. One possibility for modelling the nerve roots could be as a series of bar elements, similar to the spinal ligaments currently in the model, with one end of each bar coincident with a spinal cord node and the other end fixed relative to the vertebra at that segment. Also, because distraction injury mechanism is more distributed along the length of the spine compared to the local extent of contusion or dislocation, a longer section of the cervical spine than only levels C3-C6 would likely be more suitable for distraction simulation. ❼ Simulation of experiments measuring cerebrospinal fluid pressures during SCI. In particular the current model, or a simplified derivative version (perhaps based on the porcine spine with simplified geometry), could be used to simulate experiments performed by Jones [54]. Validation of the Smoothed Particle Hydrodynamics CSF model would then enable further simulation variations to explore the interaction between the CSF and cord during SCI. ❼ Use of the model to investigate the effect of different indenter shapes. Different laboratories have used various indenter shapes for contusion SCI experimental models, and simulation of these variations would help to understand any differences in injury pattern. Simulations within the injury velocity range of 0.6-1 m/s, as validated for the current model, would be straightforward to adapt to different indenter shapes, though injury mechanisms outside of this range would require further material property refinement and validation. ❼ A combined experimental and simulation approach to refining spinal cord material properties, especially cord viscoelasticity. There has yet to be a detailed material testing study of spinal cord tissue hyperviscoelasticity over a wide range of strain and strain rates, including up to clinically relevant high strains and rates. Spinal cord tissue testing to date has typically been tensile only, but a biaxial tension methodology may also be useful in better characterizing the tissue. Furthermore, if the dynamic viscoelastic testing of the spinal cord does demonstrate that linear viscoelasticity is substantially inadequate, a custom material model incorporating nonlinear viscoelasticity could be developed. It should be noted that the work to implement, test, and validate such a custom material model would likely constitute a major portion of a thesis project. 53 5.4. Concluding statement 5.4 Concluding statement An improved understanding of the differences in injury patterns between mechanisms – afforded by combined experimental and computational approaches – has the potential to influence future treatment and prevention approaches. For example, therapies might specifically target the central cavitation injury pattern and spared white matter rim associated with contusion injuries, or a full cord width dislocation injury with some rostral and caudal white matter damage, rather than attempting to treat more generic damage patterns. Better stratification of the patient population by factors such as injury mechanism may better identify the strengths and weaknesses of different SCI therapies and lead to improved clinical trial outcomes overall. Furthermore, with continued improvements to the computational modelling of the spinal cord during traumatic injury, the cord may one day be added to full-body models used to aid design and safety testing of products ranging from automobiles to helmets, raising the profile of SCI in safety standards and promoting preventative strategies. 54 Bibliography [1] Spinal Cord Injury Facts and Figures at a Glance. Technical report, National Spinal Cord Injury Statistical Center, 2009. [2] T. E. Anderson. A controlled pneumatic technique for experimental spinal cord contusion. J Neurosci Methods, 6(4):327–333, Nov 1982. [3] K. B. Arbogast and S. S. Margulies. Material characterization of the brainstem from oscillatory shear tests. J Biomech, 31(9):801–807, Sep 1998. [4] A. C. Bain and D. F. Meaney. Tissue-level thresholds for axonal damage in an experimental model of central nervous system white matter injury. J Biomech Eng, 122(6):615–622, Dec 2000. [5] D. M. Basso, M. S. Beattie, and J. C. Bresnahan. A sensitive and reliable locomotor rating scale for open field testing in rats. J Neurotrauma, 12(1):1–21, Feb 1995. [6] D. M. Basso, M. S. Beattie, and J. C. Bresnahan. Graded histological and locomotor outcomes after spinal cord contusion using the nyu weight-drop device versus transection. Exp Neurol, 139(2):244–256, Jun 1996. doi: 10.1006/exnr.1996.0098. URL http://dx.doi.org/10.1006/ exnr.1996.0098. [7] L. E. Bilston and L. E. Thibault. The mechanical properties of the human cervical spinal cord in vitro. Ann Biomed Eng, 24(1):67–74, 1996. [8] A. R. Blight and V. Decrescito. Morphometric analysis of experimental spinal cord injury in the cat: the relation of injury intensity to survival of myelinated axons. Neuroscience, 19(1): 321–341, Sep 1986. [9] J. C. Bresnahan, M. S. Beattie, FD Todd, 3rd, and D. H. Noyes. A behavioral and anatomical analysis of spinal cord injury produced by a feedback-controlled impaction device. Exp Neurol, 95(3):548–570, Mar 1987. [10] G. L. Chang, T. K. Hung, and W. W. Feng. An in-vivo measurement and analysis of viscoelastic properties of the spinal cord of cats. J Biomech Eng, 110(2):115–122, May 1988. [11] H. Cheng, Y. Cao, and L. Olson. Spinal cord repair in adult paraplegic rats: partial restoration of hind limb function. Science, 273(5274):510–513, Jul 1996. [12] Shaokoon Cheng, Elizabeth C. Clarke, and Lynne E. Bilston. Rheological properties of the tissues of the central nervous system: a review. Med Eng Phys, 30(10):1318–1337, Dec 2008. doi: 10.1016/j.medengphy.2008.06.003. URL http://dx.doi.org/10.1016/j.medengphy. 2008.06.003. [13] Shaokoon Cheng, Elizabeth C. Clarke, and Lynne E. Bilston. The effects of preconditioning strain on measured tissue properties. J Biomech, 42(9):1360–1362, Jun 2009. doi: 10.1016/j. jbiomech.2009.03.023. URL http://dx.doi.org/10.1016/j.jbiomech.2009.03.023. 55 Bibliography [14] Anthony M. Choo, Jie Liu, Clarrie K. Lam, Marcel Dvorak, Wolfram Tetzlaff, and Thomas R. Oxland. Contusion, dislocation, and distraction: primary hemorrhage and membrane permeability in distinct mechanisms of spinal cord injury. J Neurosurg Spine, 6(3):255–266, Mar 2007. doi: 10.3171/spi.2007.6.3.255. URL http://dx.doi.org/10.3171/spi.2007.6. 3.255. [15] Anthony M. Choo, Jie Liu, Marcel Dvorak, Wolfram Tetzlaff, and Thomas R. Oxland. Secondary pathology following contusion, dislocation, and distraction spinal cord injuries. Exp Neurol, 212(2):490–506, Aug 2008. doi: 10.1016/j.expneurol.2008.04.038. URL http: //dx.doi.org/10.1016/j.expneurol.2008.04.038. [16] Anthony Min-Te Choo, Jie Liu, Zhuowei Liu, Marcel Dvorak, Wolfram Tetzlaff, and Thomas R. Oxland. Modeling spinal cord contusion, dislocation, and distraction: characterization of vertebral clamps, injury severities, and node of ranvier deformations. J Neurosci Methods, 181(1):6–17, Jun 2009. doi: 10.1016/j.jneumeth.2009.04.007. URL http://dx.doi.org/10.1016/j.jneumeth.2009.04.007. [17] Elizabeth C. Clarke. Spinal cord mechanical properties. In Lynne E. Bilston, editor, Neural Tissue Biomechanics, volume 3 of Studies in Mechanobiology, Tissue Engineering and Biomaterials, pages 25–40. Springer Berlin Heidelberg, 2011. ISBN 978-3-642-13890-4. URL http://dx.doi.org/10.1007/8415/_2010/_15. [18] Elizabeth C. Clarke, Anthony M. Choo, Jie Liu, Clarrie K. Lam, Lynne E. Bilston, Wolfram Tetzlaff, and Thomas R. Oxland. Anterior fracture-dislocation is more severe than lateral: a biomechanical and neuropathological comparison in rat thoracolumbar spine. J Neurotrauma, 25(4):371–383, Apr 2008. doi: 10.1089/neu.2007.0421. URL http://dx.doi.org/10.1089/ neu.2007.0421. [19] Elizabeth C. Clarke, Shaokoon Cheng, and Lynne E. Bilston. The mechanical properties of neonatal rat spinal cord in vitro, and comparisons with adult. J Biomech, 42(10):1397– 1402, Jul 2009. doi: 10.1016/j.jbiomech.2009.04.008. URL http://dx.doi.org/10.1016/j. jbiomech.2009.04.008. [20] Brittany Coats and Susan S. Margulies. Material properties of porcine parietal cortex. J Biomech, 39(13):2521–2525, 2006. doi: 10.1016/j.jbiomech.2005.07.020. URL http://dx. doi.org/10.1016/j.jbiomech.2005.07.020. [21] J. V. Coumans, T. T. Lin, H. N. Dai, L. MacArthur, M. McAtee, C. Nash, and B. S. Bregman. Axonal regeneration and functional recovery after complete spinal cord transection in rats by delayed treatment with transplants and neurotrophins. J Neurosci, 21(23):9334–9344, Dec 2001. [22] D Kacy Cullen and Michelle C. LaPlaca. Neuronal response to high rate shear deformation depends on heterogeneity of the local strain field. J Neurotrauma, 23(9):1304–1319, Sep 2006. doi: 10.1089/neu.2006.23.1304. URL http://dx.doi.org/10.1089/neu.2006.23.1304. [23] K. K. Darvish and J. R. Crandall. Nonlinear viscoelastic effects in oscillatory shear deformation of brain tissue. Med Eng Phys, 23(9):633–645, Nov 2001. [24] Virtual Performance Solutions - Explicit Solver Notes Manual. ESI Group, 2010. [25] Angela Farry and David Baxter. The Incidence and Prevalence of Spinal Cord Injury in Canada. Technical report, Rick Hansen Institute and Urban Futures, 2010. 56 Bibliography [26] M. G. Fehlings and C. H. Tator. The effect of direct current field polarity on recovery after acute experimental spinal cord injury. Brain Res, 579(1):32–42, May 1992. [27] M. G. Fehlings and C. H. Tator. The relationships among the severity of spinal cord injury, residual neurological function, axon counts, and counts of retrogradely labeled neurons after experimental spinal cord injury. Exp Neurol, 132(2):220–228, Apr 1995. [28] M. G. Fehlings, C. H. Tator, and R. D. Linden. The relationships among the severity of spinal cord injury, motor and somatosensory evoked potentials and spinal cord blood flow. Electroencephalogr Clin Neurophysiol, 74(4):241–259, 1989. [29] R. J. Fiford, L. E. Bilston, P. Waite, and J. Lu. A vertebral dislocation model of spinal cord injury in rats. J Neurotrauma, 21(4):451–458, Apr 2004. doi: 10.1089/089771504323004593. URL http://dx.doi.org/10.1089/089771504323004593. [30] R.J. Fiford and L.E. Bilston. The mechanical properties of rat spinal cord in vitro. Journal of Biomechanics, 38(7):1509–1515, 2005. [31] Rodney J. Fiford and Lynne E. Bilston. The mechanical properties of rat spinal cord in vitro. J Biomech, 38(7):1509–1515, Jul 2005. doi: 10.1016/j.jbiomech.2004.07.009. URL http://dx.doi.org/10.1016/j.jbiomech.2004.07.009. [32] J.R. Flynn and P.S. Bolton. Measurement of the vertebral canal dimensions of the neck of the rat with a comparison to the human. The Anatomical Record, 2007. [33] F. Franconi, L. Lemaire, L. Marescaux, P. Jallet, and J.J. Le Jeune. In vivo quantitative microimaging of rat spinal cord at 7 T. Magnetic Resonance in Medicine, 44(6):893–898, 2000. [34] Y. Fung. Biomechanics: mechanical properties of living tissues, volume 12. Springer, 1993. [35] J. A. Galbraith, L. E. Thibault, and D. R. Matteson. Mechanical and electrical responses of the squid giant axon to simple elongation. J Biomech Eng, 115(1):13–22, Feb 1993. [36] Beth Galle, Hui Ouyang, Riyi Shi, and Eric Nauman. Correlations between tissue-level stresses and strains and cellular damage within the guinea pig spinal cord white matter. J Biomech, 40(13):3029–3033, 2007. doi: 10.1016/j.jbiomech.2007.03.014. URL http://dx. doi.org/10.1016/j.jbiomech.2007.03.014. [37] Donna M. Geddes-Klein, Kimberly B. Schiffman, and David F. Meaney. Mechanisms and consequences of neuronal stretch injury in vitro differ with the model of trauma. J Neurotrauma, 23(2):193–204, Feb 2006. doi: 10.1089/neu.2006.23.193. URL http://dx.doi.org/ 10.1089/neu.2006.23.193. [38] S.M. Goh, M.N. Charalambides, and J.G. Williams. Determination of the constitutive constants of non-linear viscoelastic materials. Mechanics of Time-Dependent Materials, 8(3): 255–268, 2004. [39] Henry Gray. Anatomy of the Human Body. Philadelphia: Lea & Febiger; Bartleby.com, 2000. www.bartleby.com/107/. [2007/08/16], 1918. [40] Carolyn Y. Greaves, Mohamed S. Gadala, and Thomas R. Oxland. A three-dimensional finite element model of the cervical spine with spinal cord: an investigation of three injury mechanisms. Ann Biomed Eng, 36(3):396–405, Mar 2008. doi: 10.1007/s10439-008-9440-0. URL http://dx.doi.org/10.1007/s10439-008-9440-0. 57 Bibliography [41] C.Y. Greaves. Spinal Cord Injury Mechanisms: A Finite Element Study. Master’s thesis, University of British Columbia, 2004. [42] Eunice Greene. Anatomy of the Rat, pages 19,149. The American Philosophical Society, Philadelphia, 1935. [43] J. A. Gruner. A monitored contusion model of spinal cord injury in the rat. J Neurotrauma, 9(2):123–6; discussion 126–8, 1992. [44] A. Guha, C. H. Tator, L. Endrenyi, and I. Piper. Decompression of the spinal cord improves recovery after acute experimental spinal cord compression injury. Paraplegia, 25(4):324–339, Aug 1987. doi: 10.1038/sc.1987.61. URL http://dx.doi.org/10.1038/sc.1987.61. [45] Johnson Ho and Svein Kleiven. Dynamic response of the brain with vasculature: a threedimensional computational study. J Biomech, 40(13):3006–3012, 2007. doi: 10.1016/j. jbiomech.2007.02.011. URL http://dx.doi.org/10.1016/j.jbiomech.2007.02.011. [46] Johnson Ho and Svein Kleiven. Can sulci protect the brain from traumatic injury? J Biomech, 42(13):2074–2080, Sep 2009. doi: 10.1016/j.jbiomech.2009.06.051. URL http: //dx.doi.org/10.1016/j.jbiomech.2009.06.051. [47] A. Howell. Anatomy of the Wood Rat, pages 124–127. The Williams & Wilkins Company, Baltimore, 1926. [48] T. K. Hung, G. L. Chang, J. L. Chang, and M. S. Albin. Stress-strain relationship and neurological sequelae of uniaxial elongation of the spinal cord of cats. Surg Neurol, 15(6): 471–476, Jun 1981. [49] Jung Keun Hyun and Hae-Won Kim. Clinical and experimental advances in regeneration of spinal cord injury. J Tissue Eng, 2010:650857, 2010. doi: 10.4061/2010/650857. URL http://dx.doi.org/10.4061/2010/650857. [50] K. Ichihara, T. Taguchi, Y. Shimada, I. Sakuramoto, S. Kawano, and S. Kawai. Gray matter of the bovine cervical spinal cord is mechanically more rigid and fragile than the white matter. J Neurotrauma, 18(3):361–367, Mar 2001. doi: 10.1089/08977150151071053. URL http://dx.doi.org/10.1089/08977150151071053. [51] Kazuhiko Ichihara, Toshihiko Taguchi, Itsuo Sakuramoto, Shunichi Kawano, and Shinya Kawai. Mechanism of the spinal cord injury and the cervical spondylotic myelopathy: new approach based on the mechanical features of the spinal cord white and gray matter. J Neurosurg, 99(3 Suppl):278–285, Oct 2003. [52] L. B. Jakeman, Z. Guan, P. Wei, R. Ponnappan, R. Dzwonczyk, P. G. Popovich, and B. T. Stokes. Traumatic spinal cord injury produced by controlled contusion in mouse. J Neurotrauma, 17(4):299–319, Apr 2000. ¨ Oguz. Shape differences in the cervical and upper [53] DR Johnson, TJ McAndrew, and O. thoracic vertebrae in rats (Rattus norvegicus) and bats (Pteropus poiocephalus): can we see shape patterns derived from position in column and species membership? Journal of Anatomy, 194(02):249–253, 1999. [54] C.F. Jones. Cerebrospinal fluid mechanics during and after experimental spinal cord injury. PhD thesis, University of British Columbia, 2011. 58 Bibliography [55] Claire F. Jones, Shannon G. Kroeker, Peter A. Cripton, and Richard M. Hall. The effect of cerebrospinal fluid on the biomechanics of spinal cord: an ex vivo bovine model using bovine and physical surrogate spinal cord. Spine (Phila Pa 1976), 33(17):E580–E588, Aug 2008. doi: 10.1097/BRS.0b013e31817ecc57. URL http://dx.doi.org/10.1097/BRS. 0b013e31817ecc57. [56] Claire F. Jones, Jae H T. Lee, Brian K. Kwon, and Peter A. Cripton. Development of a largeanimal model to measure dynamic cerebrospinal fluid pressure during spinal cord injury. J Neurosurg Spine, Apr 2012. doi: 10.3171/2012.3.SPINE11970. URL http://dx.doi.org/ 10.3171/2012.3.SPINE11970. [57] P. A. Kearney, S. A. Ridella, D. C. Viano, and T. E. Anderson. Interaction of contact velocity and cord compression in determining the severity of spinal cord injury. J Neurotrauma, 5(3): 187–208, 1988. [58] J.H. Kim, T.W. Tu, P.V. Bayly, and S.K. Song. Impact speed does not determine severity of spinal cord injury in mice with fixed impact displacement. Journal of Neurotrauma, 26(8): 1395–1404, 2009. [59] Hideyuki Kimpara, Yuko Nakahira, Masami Iwamoto, Kazuo Miki, Kazuhiko Ichihara, Shunichi Kawano, and Toshihiko Taguchi. Investigation of anteroposterior head-neck responses during severe frontal impacts using a brain-spinal cord complex fe model. Stapp Car Crash J, 50:509–544, Nov 2006. [60] Svein Kleiven. Influence of impact direction on the human head in prediction of subdural hematoma. J Neurotrauma, 20(4):365–379, Apr 2003. doi: 10.1089/089771503765172327. URL http://dx.doi.org/10.1089/089771503765172327. [61] N. R. Krenz and L. C. Weaver. Sprouting of primary afferent fibers after spinal cord transection in the rat. Neuroscience, 85(2):443–458, Jul 1998. [62] Brian K. Kwon, Tom R. Oxland, and Wolfram Tetzlaff. Animal models used in spinal cord regeneration research. Spine (Phila Pa 1976), 27(14):1504–1510, Jul 2002. [63] Xin-Feng Li and Li-Yang Dai. Three-dimensional finite element model of the cervical spinal cord: preliminary results of injury mechanism analysis. Spine (Phila Pa 1976), 34(11):1140– 1147, May 2009. doi: 10.1097/BRS.0b013e31819e2af1. URL http://dx.doi.org/10.1097/ BRS.0b013e31819e2af1. [64] Paul Lu, Armin Blesch, Lori Graham, Yaozhi Wang, Ramsey Samara, Karla Banos, Verena Haringer, Leif Havton, Nina Weishaupt, David Bennett, Karim Fouad, and Mark H. Tuszynski. Motor axonal regeneration after partial and complete spinal cord transection. J Neurosci, 32(24):8208–8218, Jun 2012. doi: 10.1523/JNEUROSCI.0308-12.2012. URL http://dx.doi.org/10.1523/JNEUROSCI.0308-12.2012. [65] E. Lucas. Measuring in vivo internal spinal cord deformations during experimental spinal cord injury using a rat model, radiography, and fiducial markers. Master’s thesis, University of British Columbia, 2010. [66] Jason T. Maikos and David I. Shreiber. Immediate damage to the blood-spinal cord barrier due to mechanical trauma. J Neurotrauma, 24(3):492–507, Mar 2007. doi: 10.1089/neu.2006. 0149. URL http://dx.doi.org/10.1089/neu.2006.0149. [67] Jason T. Maikos, Ragi A I. Elias, and David I. Shreiber. Mechanical properties of dura mater from the rat brain and spinal cord. J Neurotrauma, 25(1):38–51, Jan 2008. doi: 10.1089/neu.2007.0348. URL http://dx.doi.org/10.1089/neu.2007.0348. 59 Bibliography [68] Jason T. Maikos, Zhen Qian, Dimitris Metaxas, and David I. Shreiber. Finite element analysis of spinal cord injury in the rat. J Neurotrauma, 25(7):795–816, Jul 2008. doi: 10.1089/neu. 2007.0423. URL http://dx.doi.org/10.1089/neu.2007.0423. [69] Thomas W. McAllister, James C. Ford, Songbai Ji, Jonathan G. Beckwith, Laura A. Flashman, Keith Paulsen, and Richard M. Greenwald. Maximum principal strain and strain rate associated with concussion diagnosis correlates with changes in corpus callosum white matter indices. Ann Biomed Eng, 40(1):127–140, Jan 2012. doi: 10.1007/s10439-011-0402-6. URL http://dx.doi.org/10.1007/s10439-011-0402-6. [70] K. K. Mendis, R. L. Stalnaker, and S. H. Advani. A constitutive relationship for large deformation finite element modeling of brain tissue. J Biomech Eng, 117(3):279–285, Aug 1995. [71] GA Metz, A. Curt, H. van de Meent, I. Klusman, ME Schwab, and V. Dietz. Validation of the weight-drop contusion model in rats: a comparative study of human spinal cord injury. Journal of Neurotrauma, 17(1):1–17, 2000. [72] K. Miller and K. Chinzei. Constitutive modelling of brain tissue: experiment and theory. J Biomech, 30(11-12):1115–1121, 1997. [73] Karol Miller. Method of testing very soft biological tissues in compression. J Biomech, 38 (1):153–158, Jan 2005. doi: 10.1016/j.jbiomech.2004.03.004. URL http://dx.doi.org/10. 1016/j.jbiomech.2004.03.004. [74] Karol Miller and Kiyoyuki Chinzei. Mechanical properties of brain tissue in tension. J Biomech, 35(4):483–490, Apr 2002. [75] J.J. Monaghan. An introduction to sph. Computer Physics Communications, 48(1):89–96, 1988. [76] J.J. Monaghan. Simulating free surface flows with SPH. Journal of Computational Physics, 110:399–399, 1994. [77] Leith Morriss, Adam Wittek, and Karol Miller. Compression testing of very soft biological tissues using semi-confined configuration–a word of caution. J Biomech, 41(1):235– 238, 2008. doi: 10.1016/j.jbiomech.2007.06.025. URL http://dx.doi.org/10.1016/j. jbiomech.2007.06.025. [78] S. Muskopf. The biology corner. Online Website. [2007/07/28], 2007. biologycorner.com/worksheets/rat_external.html. URL www. [79] D. H. Noyes. Correlation between parameters of spinal cord impact and resultant injury. Exp Neurol, 95(3):535–547, Mar 1987. [80] D. H. Noyes. Electromechanical impactor for producing experimental spinal cord injury in animals. Med Biol Eng Comput, 25(3):335–340, May 1987. [81] R. J. Oakland, R. M. Hall, R. K. Wilcox, and D. C. Barton. The biomechanical response of spinal cord tissue to uniaxial loading. Proc Inst Mech Eng H, 220(4):489–492, May 2006. [82] R.W. Ogden. Large deformation isotropic elasticity-on the correlation of theory and experiment for incompressible rubberlike solids. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 326(1567):565–584, 1972. 60 Bibliography [83] R.W. Ogden, G. Saccomandi, and I. Sgura. Fitting hyperelastic models to experimental data. Computational Mechanics, 34(6):484–502, 2004. [84] S. Osher and J. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988. [85] Dwight Parkinson. Immediate neurocognitive effects of concussion: graded contusion model of the mouse spinal cord using a pneumatic impact device. Neurosurgery, 52(6):1505, Jun 2003. [86] Cecilia Persson, Stewart W D. McLure, Jon Summers, and Richard M. Hall. The effect of bone fragment size and cerebrospinal fluid on spinal cord deformation during trauma: an ex vivo study. J Neurosurg Spine, 10(4):315–323, Apr 2009. doi: 10.3171/2009.1.SPINE08286. URL http://dx.doi.org/10.3171/2009.1.SPINE08286. [87] Cecilia Persson, Jon Summers, and Richard M. Hall. The effect of cerebrospinal fluid thickness on traumatic spinal cord deformation. J Appl Biomech, 27(4):330–335, Nov 2011. [88] L.A. Piegl and W. Tiller. The NURBS book. Springer Verlag, 1997. [89] Michael T. Prange and Susan S. Margulies. Regional, directional, and age-dependent properties of the brain undergoing large deformation. J Biomech Eng, 124(2):244–252, Apr 2002. [90] K.P. Quinn and B.A. Winkelstein. Cervical facet capsular ligament yield defines the threshold for injury and persistent joint-mediated neck pain. Journal of Biomechanics, 40(10):2299– 2306, 2007. [91] A. Ramn-Cueto, G. W. Plant, J. Avila, and M. B. Bunge. Long-distance axonal regeneration in the transected adult rat spinal cord is promoted by olfactory ensheathing glia transplants. J Neurosci, 18(10):3803–3815, May 1998. [92] A. S. Rivlin and C. H. Tator. Effect of duration of acute spinal cord compression in a new acute cord injury model in the rat. Surg Neurol, 10(1):38–43, Jul 1978. [93] HGQ Rowett. The Rat as a Small Mammal, pages 11–13. John Murray Ltd., London, third edition, 1974. [94] C. Russell. Modeling the rat cervical spine. Co-op work term report, University of British Columbia, 2007. [95] Colin M. Russell, Anthony M. Choo, Wolfram Tetzlaff, Tae-Eun Chung, and Thomas R. Oxland. Maximum principal strain correlates with spinal cord tissue damage in contusion and dislocation injuries in the rat cervical spine. J Neurotrauma, 29(8):1574–1585, May 2012. doi: 10.1089/neu.2011.2225. URL http://dx.doi.org/10.1089/neu.2011.2225. [96] Stephen W. Scheff, Alexander G. Rabchevsky, Isabella Fugaccia, John A. Main, and James E Lumpp, Jr. Experimental modeling of spinal cord injury: characterization of a force-defined injury device. J Neurotrauma, 20(2):179–193, Feb 2003. doi: 10.1089/08977150360547099. URL http://dx.doi.org/10.1089/08977150360547099. [97] G. Schwartz and M. G. Fehlings. Evaluation of the neuroprotective effects of sodium channel blockers after spinal cord injury: improved behavioral and neuroanatomical recovery with riluzole. J Neurosurg, 94(2 Suppl):245–256, Apr 2001. 61 Bibliography [98] Jeffrey Scifert, Koji Totoribe, Vijay Goel, and Jan Huntzinger. Spinal cord mechanics during flexion and extension of the cervical spine: a finite element study. Pain Physician, 5(4): 394–400, Oct 2002. [99] Toshitaka Seki, Kazutoshi Hida, Mitsuhiro Tada, Izumi Koyanagi, and Yoshinobu Iwasaki. Graded contusion model of the mouse spinal cord using a pneumatic impact device. Neurosurgery, 50(5):1075–81; discussion 1081–2, May 2002. [100] J. Sethian. Level Set Methods and Fast Marching Methods. Cambridge University Press, 1999. [101] Mehdi Shafieian, Kurosh K. Darvish, and James R. Stone. Changes to the viscoelastic properties of brain tissue after traumatic axonal injury. J Biomech, 42(13):2136–2142, Sep 2009. doi: 10.1016/j.jbiomech.2009.05.041. URL http://dx.doi.org/10.1016/j.jbiomech.2009.05. 041. [102] R. Shi and A. R. Blight. Compression injury of mammalian spinal cord in vitro and the dynamics of action potential conduction failure. J Neurophysiol, 76(3):1572–1580, Sep 1996. [103] D.I. Shreiber, A.C. Bain, and D.F. Meaney. In vivo thresholds for mechanical injury to the blood-brain barrier. In 41st Stapp Car Crash Conference Proceedings, volume 315, page 277. SAE International, 1997. [104] J. G. Snedeker, P. Niederer, F. R. Schmidlin, M. Farshad, C. K. Demetropoulos, J. B. Lee, and K. H. Yang. Strain-rate dependent material properties of the porcine and human kidney capsule. J Biomech, 38(5):1011–1021, May 2005. doi: 10.1016/j.jbiomech.2004.05.036. URL http://dx.doi.org/10.1016/j.jbiomech.2004.05.036. [105] Carolyn J. Sparrey and Tony M. Keaveny. Compression behavior of porcine spinal cord white matter. J Biomech, 44(6):1078–1082, Apr 2011. doi: 10.1016/j.jbiomech.2011.01.035. URL http://dx.doi.org/10.1016/j.jbiomech.2011.01.035. [106] Carolyn J. Sparrey, Geoffrey T. Manley, and Tony M. Keaveny. Effects of white, grey, and pia mater properties on tissue level stresses and strains in the compressed spinal cord. J Neurotrauma, 26(4):585–595, Apr 2009. doi: 10.1089/neu.2008.0654. URL http://dx.doi. org/10.1089/neu.2008.0654. [107] C.J. Sparrey, A.M. Choo, J. Liu, W. Tetzlaff, and T.R. Oxland. The distribution of tissue damage in the spinal cord is influenced by the contusion velocity. Spine, 33(22):E812, 2008. [108] David P. Stirling, Kourosh Khodarahmi, Jie Liu, Lowell T. McPhail, Christopher B. McBride, John D. Steeves, Matt S. Ramer, and Wolfram Tetzlaff. Minocycline treatment reduces delayed oligodendrocyte death, attenuates axonal dieback, and improves functional outcome after spinal cord injury. J Neurosci, 24(9):2182–2190, Mar 2004. doi: 10.1523/JNEUROSCI. 5275-03.2004. URL http://dx.doi.org/10.1523/JNEUROSCI.5275-03.2004. [109] B. T. Stokes. Experimental spinal cord injury: a dynamic and verifiable injury device. J Neurotrauma, 9(2):129–31; discussion 131–4, 1992. [110] B.T. Stokes, D.H. Noyes, and D.L. Behrmann. An electromechanical spinal injury technique with dynamic sensitivity. Journal of Neurotrauma, 9(3):187–195, 1992. [111] Peter K. Stys. White matter injury mechanisms. Curr Mol Med, 4(2):113–130, Mar 2004. [112] P. G. Sullivan, A. G. Rabchevsky, P. C. Waldmeier, and J. E. Springer. Mitochondrial permeability transition in cns trauma: cause or effect of neuronal cell death? J Neurosci Res, 79(1-2):231–239, 2005. doi: 10.1002/jnr.20292. URL http://dx.doi.org/10.1002/ jnr.20292. 62 Bibliography [113] Charles H. Tator. Review of treatment trials in human spinal cord injury: issues, difficulties, and recommendations. Neurosurgery, 59(5):957–82; discussion 982–7, Nov 2006. doi: 10.1227/01.NEU.0000245591.16087.89. URL http://dx.doi.org/10.1227/01.NEU. 0000245591.16087.89. [114] Kevin L. Troyer and Christian M. Puttlitz. Nonlinear viscoelasticty plays an essential role in the functional behavior of spinal ligaments. J Biomech, 45(4):684–691, Feb 2012. doi: 10.1016/j.jbiomech.2011.12.009. URL http://dx.doi.org/10.1016/j.jbiomech.2011.12. 009. [115] David C. Viano, Ira R. Casson, Elliot J. Pellman, Liying Zhang, Albert I. King, and King H. Yang. Concussion in professional football: brain responses by finite element analysis: part 9. Neurosurgery, 57(5):891–916; discussion 891–916, Nov 2005. [116] DC Viano and P. Lovsund. Biomechanics of brain and spinal-cord injury: analysis of neuropathologic and neurophysiologic experiments. Crash Prevention and Injury Control, 1(1), 1999. [117] C. Watson, G. Paxinos, and G. Kayalioglu. The spinal cord: a Christopher and Dana Reeve Foundation text and atlas. Academic Press, 2008. [118] TAG Wells. The Rat: A Practical Guide, pages 10–14. Heinemann Educational Books Ltd., London, 1964. [119] R.T. Whitaker. A level-set approach to 3D reconstruction from range data. International Journal of Computer Vision, 29(3):203–231, 1998. [120] Bruce D. Wingerd. Rat Dissection Manual, pages vii–viii,9–10. Johns Hopkins, Baltimore, 1988. [121] Sang Jun Yeo, Sung Nam Hwang, Seung Won Park, Young Baeg Kim, Byung Kook Min, Jeong Taik Kwon, and Jong Sik Suk. Development of a rat model of graded contusive spinal cord injury using a pneumatic impact device. J Korean Med Sci, 19(4):574–580, Aug 2004. [122] N. Yoganandan, S. Kumaresan, and F. A. Pintar. Geometric and mechanical properties of human cervical spine ligaments. J Biomech Eng, 122(6):623–629, Dec 2000. [123] N. Yoganandan, S. Kumaresan, and F.A. Pintar. Geometric and mechanical properties of human cervical spine ligaments. Journal of Biomechanical Engineering, 122:623, 2000. [124] PA Yushkevich, J. Piven, HC Hazlett, RG Smith, S. Ho, JC Gee, and G. Gerig. User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability. Neuroimage, 31(3):1116–28, 2006. [125] Paul A. Yushkevich, Joseph Piven, Heather Cody Hazlett, Rachel Gimpel Smith, Sean Ho, James C. Gee, and Guido Gerig. User-guided 3d active contour segmentation of anatomical structures: significantly improved efficiency and reliability. Neuroimage, 31(3):1116–1128, Jul 2006. doi: 10.1016/j.neuroimage.2006.01.015. URL http://dx.doi.org/10.1016/j. neuroimage.2006.01.015. [126] Yi Ping Zhang, Darlene A. Burke, Lisa B E. Shields, Sergey Y. Chekmenev, Toros Dincman, Yongjie Zhang, Yiyan Zheng, Rebecca R. Smith, Richard L. Benton, William H. DeVries, Xiaoling Hu, David S K. Magnuson, Scott R. Whittemore, and Christopher B. Shields. Spinal cord contusion based on precise vertebral stabilization and tissue displacement measured by combined assessment to discriminate small functional differences. J Neurotrauma, 25(10): 63 1227–1240, Oct 2008. doi: 10.1089/neu.2007.0388. URL http://dx.doi.org/10.1089/neu. 2007.0388. [127] Qiliang Zhu, Michael Prange, and Susan Margulies. Predicting unconsciousness from a pediatric brain injury threshold. Dev Neurosci, 28(4-5):388–395, 2006. doi: 10.1159/000094165. URL http://dx.doi.org/10.1159/000094165. [128] SC Zhu and A. Yuille. Region competition: unifying snakes, region growing, and Bayes/MDL for multiband image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 18(9):884–900, 1996. 64 Appendix A Images of Rat Vertebrae from Literature (a) Nervous System (b) Left Lateral Aspect of Atlas and Epistropheus (c) Anterior Aspect of Fifth Cervical Vertebra Figure A.1: Cervical vertebrae and nervous system [Illustrations and caption text reproduced from Greene [42].] 65 Appendix A. Images of Rat Vertebrae from Literature (a) The Atlas. Anterior view. (b) The Atlas. Dorsal view. (c) The Atlas. Posterior view. (d) The Axis. Anterior view. (e) The Axis. Lateral view. (f ) The Axis. Posterior view. (g) Fourth Cervical Anterior view. Vertebra.(h) Fourth Cervical Dorsal view. (j) Fourth Thoracic Anterior view. Vertebra.(i) Sixth Cervical Anterior view. Vertebra.(k) Fourth Thoracic Lateral view. Vertebra. Vertebra. Figure A.2: Cervical and thoracic vertebrae [Illustrations reproduced from Wells [118].] 66 Appendix A. Images of Rat Vertebrae from Literature (a) The Atlas. Anterior view. (d) The Axis. Anterior view. (b) The Atlas. Dorsal view. (e) The Axis. Dorsal view. (c) The Atlas. Lateral view. (f ) The Axis. Lateral view. (g) Third Cervical Vertebra. Anterior and Lateral views. Figure A.3: Cervical vertebrae [Illustrations reproduced from Rowett [93].] 67 Appendix A. Images of Rat Vertebrae from Literature (a) C1 (b) C2 (c) C3 (d) C4 (e) C5 (f ) C6 (g) C7 (h) T1 (i) T2 Figure A.4: Cervical and thoracic vertebrae [Illustrations reproduced with permission from Johnson et al. [53].] 68 Appendix B Procedure for extracting simulation results for solid elements in each predefined spinal cord zone slice 1. In Visual-Viewer open desired simulation .THP time history results file. 2. Choose SOLID Entity Type and TIME X Component at the top of Import and Plot window. 3. Choose CC Group by Component. 4. Under Entities selection box, click Expand. 5. In Entity List window type zone slice name into filter, such as “dorsal epicentre” or “dorsal +1mm”, and click enter17 . 6. Click purple filled box to Select All. 7. In Import and Plot window select all three principal stresses and all three Solid Auxiliary variables (principal strains) in the Components box (hold CTRL to select multiple components). 8. Click Plot. 9. Select All Pages tab in Explorer. 10. Right-click on one of the selected plot groups and choose Export Curves. 11. Choose CSV (Comma delimited) and Use Curve Titles. 12. Type an appropriate file name such as “dorsal +1mm” and click Save to write the .csv file. 13. Click Delete to delete all the selected plot groups. Repeat steps 5-13 for each zone slice. 17 Note that “lateral” and “ventrolateral” are not distinguishable at this step, and the correct set must be manually selected instead of using Select All after performing filter. 69 Appendix C MATLAB code C.1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 Hyperviscoelastic curve fitting - hvstress.m function [ stress ] = hvstress( params, vars ) % hvstress Computes hyperviscoelastic stress from strain history % S = hvstress( params, vars ) returns the hyperviscoelastic % stress given the Ogden and Prony parameters and an arbitrary strain % history over time. This function is designed so that it may be used % with the nonlinear least squares curve fitting function lsqcurvefit to % optimize the hyperviscoelastic Ogden and Prony parameters against a known strain % history and corresponding stresses (this is why n,m are passed as % global variables rather than part of params, and also why timestep dt % is calculated externally rather than during each function call). % % [params] is a row vector containing all Ogden terms in the order mu, % alpha followed by all Prony terms in the order g, tau. % % [vars] is a matrix containing three row vectors of equal length: % t = vars(1,:) discrete timepoints % lambda = vars(2,:) streth ratio at each timepoint % dt = vars(3,:) timestep between timepoints % % n must be a global variable equal to the desired number of Ogden term % pairs. % % m must be a global variable equal to the desired number of Prony term % pairs. % % Written June 29, 2010 by Colin Russell 27 28 29 30 31 32 % n = number of Ogden term pairs % m = number of Prony term pairs global n m; 33 34 35 36 37 38 39 40 41 for i=1:n mu(i) = params(2*i−1); alpha(i) = params(2*i); end for i=1:m g(i,1) = params(2*i−1+2*n); tau(i,1) = params(2*i+2*n); end 42 43 44 45 t = vars(1,:); lambda = vars(2,:); dt = vars(3,:); 46 47 %% Compute stress 48 70 C.2. Hyperviscoelastic optimized parameter display - hvlabels.m 49 50 51 52 53 % Hyperelastic stress s0 = zeros(1,length(t)); for i=1:n s0 = s0 + (2*mu(i))*(lambda.ˆ(alpha(i)−1)−lambda.ˆ(−alpha(i)/2 −1)); end 54 55 56 57 58 59 % Return only hyperelastic stress if no Prony series terms used if m==0 stress = s0; return; end 60 61 ginf = 1−sum(g); 62 63 64 65 66 67 68 69 70 % Compute discretized convolution integral to incorporate viscoelasticity stress = zeros(1,length(t)); h = zeros(m,length(t)); for i=1:length(t)−1 h(:,i+1) = exp(−dt(i)*(tau.ˆ−1)).*h(:,i) +... g.*tau.*( (1−exp(−dt(i)*(tau.ˆ−1))) )*( s0(i+1)−s0(i) )/dt(i); stress(i+1) = ginf*s0(i+1) + sum(h(:,i)); end 71 72 end C.2 1 2 3 Hyperviscoelastic optimized parameter display - hvlabels.m function [ returntable ] = hvlabels( N, M, params ) % hvlabels Displays Ogden and Prony parameters in tabulated form, with labels. % Written June 9, 2010 by Colin Russell 4 5 table = cell(4,2*M+2*N); 6 7 8 9 10 11 12 13 14 for i=1:N mu(i) = params(2*i−1); alpha(i) = params(2*i); end for i=1:M g(i,1) = params(2*i−1+2*N); tau(i,1) = params(2*i+2*N); end 15 16 17 18 G0 = sum(mu.*alpha); G = g*G0; Ginf = G0 − sum(G); 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 table{3,1} = ['G0 (MPa)']; table{4,1} = G0; for i=1:N table{1,2*i−1} = ['mu' num2str(i) ' (MPa)']; table{1,2*i} = ['alpha' num2str(i)]; table{2,2*i−1} = mu(i); table{2,2*i} = alpha(i); end for i=1:M table{1,2*i−1+2*N} = ['g' num2str(i)]; table{1,2*i+2*N} = ['tau' num2str(i) ' (ms)']; table{2,2*i−1+2*N} = g(i); table{2,2*i+2*N} = tau(i); table{3,2*i−1+2*N} = ['G' num2str(i) ' (MPa)']; 71 C.3. Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m table{4,2*i−1+2*N} = G(i); 34 35 end 36 37 disp(table) 38 39 40 41 42 %% Return optional output if nargout>0 returntable = table; end 43 44 end C.3 Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m 1 2 clc; clear; close all; 3 4 global m n t relax; 5 6 7 % Units: Stress (MPa) % Time (ms) 8 9 10 11 12 % % n m n m = = = number of Ogden term pairs = number of Prony term pairs 1; 4; 13 14 15 % Relaxation time (ms) t relax = 1800e3; 16 17 18 % Approx. ratio of steady state to initial stress at 5%, 0.2/s relax ratio = 0.045/0.076; 19 20 21 % Set G0 equal to initial stress over 5% strain at 0.2/s rate G0 = 0.076/0.05; 22 23 24 25 % Load Fiford and Bilston data % t in ms, s and sSD in MPa load fbdata; 26 27 28 29 30 % 5% max strain, 0.2/s e max = 0.05; e rate = 0.2e−3; % (per ms) t1 = e max/e rate; 31 32 33 34 35 36 37 38 39 40 41 42 43 44 % % % % % % % % % % % % % % Piecewise cubic Hermite interpolation, with logarithmic spacing during % relaxation [ans,t1L] = max(s 5 02); t1 = t 5 02(t1L); xx = [0:0.25:t1−0.25, logspace(log10(t1),log10(t relax),1000)]; yy = interp1(t 5 02 ,s 5 02 ,xx,'cubic'); figure hold on plot(xx,yy,'.r') plot(t 5 02 ,s 5 02) hold off t 5 02 = xx; s 5 02 = yy; 45 46 72 C.3. Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m 47 48 49 50 51 52 53 54 55 t = t 5 02; lambda 5 02 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); % true strain = log(lambda 5 02); % s 5 02 = s 5 02./(1+true strain); % Convert true stress to nominal stress dt 5 02 = zeros(1,length(t)); for i=1:length(t)−1 dt 5 02(i) = t(i+1)−t(i); end dt 5 02(end) = dt 5 02(end−1); 56 57 58 59 60 61 62 63 64 65 66 67 68 % 5% max strain, 0.02/s e max = 0.05; e rate = 0.02e−3; % (per ms) t1 = e max/e rate; t = t 5 002; lambda 5 002 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 5 002 = zeros(1,length(t)); for i=1:length(t)−1 dt 5 002(i) = t(i+1)−t(i); end dt 5 002(end) = dt 5 002(end−1); 69 70 71 72 73 74 75 76 77 78 79 80 % 5% max strain, 0.002/s e max = 0.05; e rate = 0.002e−3; % (per ms) t1 = e max/e rate; t = t 5 0002; lambda 5 0002 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 5 0002 = zeros(1,length(t)); for i=1:length(t)−1 dt 5 0002(i) = t(i+1)−t(i); end dt 5 0002(end) = dt 5 0002(end−1); 81 82 83 84 85 86 87 88 89 90 91 92 % 3.5% max strain, 0.2/s e max = 0.035; e rate = 0.2e−3; % (per ms) t1 = e max/e rate; t = t 35 02; lambda 35 02 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 35 02 = zeros(1,length(t)); for i=1:length(t)−1 dt 35 02(i) = t(i+1)−t(i); end dt 35 02(end) = dt 35 02(end−1); 93 94 95 96 97 98 99 100 101 102 103 104 % 3.5% max strain, 0.002/s e max = 0.035; e rate = 0.002e−3; % (per ms) t1 = e max/e rate; t = t 35 0002; lambda 35 0002 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 35 0002 = zeros(1,length(t)); for i=1:length(t)−1 dt 35 0002(i) = t(i+1)−t(i); end dt 35 0002(end) = dt 35 0002(end−1); 105 106 107 108 % 2% max strain, 0.2/s e max = 0.02; e rate = 0.2e−3; % (per ms) 73 C.3. Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m 109 110 111 112 113 114 115 116 t1 = e max/e rate; t = t 2 02; lambda 2 02 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 2 02 = zeros(1,length(t)); for i=1:length(t)−1 dt 2 02(i) = t(i+1)−t(i); end dt 2 02(end) = dt 2 02(end−1); 117 118 119 120 121 122 123 124 125 126 127 128 % 2% max strain, 0.002/s e max = 0.02; e rate = 0.002e−3; % (per ms) t1 = e max/e rate; t = t 2 0002; lambda 2 0002 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 2 0002 = zeros(1,length(t)); for i=1:length(t)−1 dt 2 0002(i) = t(i+1)−t(i); end dt 2 0002(end) = dt 2 0002(end−1); 129 130 131 132 133 134 135 136 137 138 139 140 141 142 %% Test model at higher strain and strain rate % 20% max strain, 80/s t 20 80 = [0:0.05:2.5]; e max = 0.2; e rate = 80e−3; t1 = e max/e rate; t = t 20 80; lambda 20 80 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 20 80 = zeros(1,length(t)); for i=1:length(t)−1 dt 20 80(i) = t(i+1)−t(i); end dt 20 80(end) = dt 20 80(end−1); 143 144 145 146 147 148 149 150 151 152 153 154 155 % 20% max strain, 0.8/s t 20 08 = [0:1:250]; e max = 0.2; e rate = 0.8e−3; t1 = e max/e rate; t = t 20 08; lambda 20 08 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 20 08 = zeros(1,length(t)); for i=1:length(t)−1 dt 20 08(i) = t(i+1)−t(i); end dt 20 08(end) = dt 20 08(end−1); 156 157 158 159 160 161 162 163 164 165 166 167 168 % 5% max strain, 80/s t 5 80 = [0:0.05:2.5]; e max = 0.05; e rate = 80e−3; t1 = e max/e rate; t = t 5 80; lambda 5 80 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 5 80 = zeros(1,length(t)); for i=1:length(t)−1 dt 5 80(i) = t(i+1)−t(i); end dt 5 80(end) = dt 5 80(end−1); 169 170 % 5% max strain, 0.8/s 74 C.3. Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m 171 172 173 174 175 176 177 178 179 180 181 t 5 08 = [0:1:250]; e max = 0.05; e rate = 0.8e−3; t1 = e max/e rate; t = t 5 08; lambda 5 08 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 5 08 = zeros(1,length(t)); for i=1:length(t)−1 dt 5 08(i) = t(i+1)−t(i); end dt 5 08(end) = dt 5 08(end−1); 182 183 184 185 186 187 188 189 %% Concatenated variables % s = [s 5 02 , s 5 002 , s 5 0002, s 35 02 , s 35 % t = [t 5 02 , t 5 002 , t 5 0002, t 35 02 , t 35 % lambda = [lambda 5 02, lambda 5 002, lambda 5 % lambda 35 0002, lambda 2 02, lambda 2 0002]; % dt = [dt 5 02 , dt 5 002, dt 5 0002, dt 35 02, 0002, s 2 02 , s 2 0002]; 0002, t 2 02 , t 2 0002]; 0002, lambda 35 02,... dt 35 0002, dt 2 02 , dt 2 0002]; 190 191 192 193 194 % % % % s = [s 5 t = [t 5 lambda = dt = [dt 02 , s 35 02 , s 2 02]; 02 , t 35 02 , t 2 02]; [lambda 5 02, lambda 35 02, lambda 2 02]; 5 02 , dt 35 02, dt 2 02]; % % % % s = [s 5 t = [t 5 lambda = dt = [dt 02 , s 5 002 , s 5 0002]; 02 , t 5 002 , t 5 0002]; [lambda 5 02, lambda 5 002, lambda 5 0002]; 5 02 , dt 5 002, dt 5 0002]; % % % % s = [s 5 t = [t 5 lambda = dt = [dt 002]; 002]; [lambda 5 002]; 5 002]; 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 s = [s 5 t = [t 5 lambda = dt = [dt 02]; 02]; [lambda 5 02]; 5 02]; 210 211 212 213 214 215 % % % % s t e e = [s 5 0002]; = [t 5 0002]; max = [e max 5 0002]; rate = [e rate 5 0002]; 216 217 218 219 %% Concatenate hv variables vars = [t; lambda; dt]; 220 221 222 %% Curve fitting 223 224 vfun = @hvstress; 225 226 227 228 229 230 231 232 % Initialize model parameters with starting guesses that are spaced out for i=1:n % Set initial mu i params(2*i−1) = 200e−3; % Set initial alpha i params(2*i) = 2*i; end 75 C.3. Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m 233 234 235 236 237 238 for i=1:m % Set initial g i params(2*i−1+2*n) = 0.9; % Set initial tau i params(2*i+2*n) = 10ˆ(i−1); end 239 240 241 242 % default upper and lower bounds for model parameters lb = zeros(1,length(params)); ub = inf*ones(1,length(params)); 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 % Set smarter upper and lower bounds for i=1:n lb(2*i−1) = 0; ub(2*i−1) = inf; lb(2*i) = 0; ub(2*i) = inf; end % for i=1:n % Try negative mu, alpha % lb(2*i−1) = −inf; % ub(2*i−1) = 0; % lb(2*i) = −inf; % ub(2*i) = 0; % end for i=1:m lb(2*i−1+2*n) = 0; ub(2*i−1+2*n) = inf; lb(2*i+2*n) = 0.1; ub(2*i+2*n) = inf; end % % Try Miller value of alpha % for i=1:n % lb(2*i−1) = −inf; % ub(2*i−1) = 0; % lb(2*i) = −4.71; % ub(2*i) = −4.7; % end %% Assign Mendis/Maikos Prony series parameters for 8ms %% time constant factors since Fiford data is not dense enough in initial %% strain ramp to fit this short term relaxation constant for i=1 lb(2*i−1+2*n) = 0.5282; ub(2*i−1+2*n) = 0.5283; lb(2*i+2*n) = 8; ub(2*i+2*n) = 8.1; end % for i=2 % lb(2*i−1+2*n) = 0.3018; % ub(2*i−1+2*n) = 0.3019; % lb(2*i+2*n) = 150; % ub(2*i+2*n) = 150.1; % end % for i=3:m % lb(2*i−1+2*n) = 0; % ub(2*i−1+2*n) = 1−.5282−.308; % lb(2*i+2*n) = 10ˆ(i−1); % ub(2*i+2*n) = inf; % end 291 292 % randparams = lb+(ub−lb).*rand(1,length(params)); 293 294 % Setting lsqcurvefit algorithm options 76 C.3. Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m 295 296 297 298 299 300 301 302 303 304 305 306 307 % options = optimset('MaxFunEvals',1e9,'MaxIter',30000,'Display',... % 'notify','TolFun',1e−14,'TolX',1e−14,'PlotFcns',{[]} ); options = optimset('MaxFunEvals',1e9,'MaxIter',1e6,'Display','notify',... 'TypicalX',params,'TolFun',1e−16,'TolX',1e−16,... 'PlotFcns',{[@optimplotfval]} ); % PlotFcns % Plots various measures of progress while the algorithm executes, select % from predefined plots or write your own. Specifying @optimplotx plots the % current point; @optimplotfunccount plots the function count; % @optimplotfval plots the function value; @optimplotconstrviolation plots % the maximum constraint violation; @optimplotresnorm plots the norm of the % residuals; @optimplotstepsize plots the step size; % @optimplotfirstorderopt plots the first−order of optimality. 308 309 310 311 312 313 314 315 316 %lsqcurvefit function tic [solvedconstants,resnorm,residual,output] = lsqcurvefit(vfun,params,... vars,s,[lb],[ub],options); toc 317 318 319 disp 'Fiford and Bilston fit constants' hvlabels(n,m,solvedconstants); 320 321 322 323 324 325 326 327 328 329 330 331 %% figure title(['Fit to Fiford and Bilston (2005) data, 5% max strain'... ' and 0.2/s strain rate']) ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(solvedconstants,[t 5 02; lambda 5 02; dt 5 02]),'−r') legend('5% 0.2/s','5% 0.2/s fit') hold off 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 figure title('Fit to Fiford and Bilston (2005) data, 0.2/s strain rate') ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(solvedconstants,... [t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 35 02/1000,s 35 02 ,'.g') plot(t 35 02/1000,vfun(solvedconstants,... [t 35 02; lambda 35 02; dt 35 02]),'−g') plot(t 2 02/1000,s 2 02 ,'.k') plot(t 2 02/1000,vfun(solvedconstants,... [t 2 02; lambda 2 02; dt 2 02]),'−k') legend('5% 0.2/s','5% 0.2/s fit','3.5% 0.2/s','3.5% 0.2/s fit',... '2% 0.2/s','2% 0.2/s fit') hold off 350 351 352 353 354 355 356 figure title('Fit to Fiford and Bilston (2005) data, 5% max strain') ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') 77 C.3. Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m 357 358 359 360 361 362 363 364 365 366 367 plot(t 5 02/1000,vfun(solvedconstants,... [t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 5 002/1000,s 5 002 ,'.g') plot(t 5 002/1000,vfun(solvedconstants,... [t 5 002; lambda 5 002; dt 5 002]),'−g') plot(t 5 0002/1000,s 5 0002,'.k') plot(t 5 0002/1000,vfun(solvedconstants,... [t 5 0002; lambda 5 0002; dt 5 0002]),'−k') legend('5% 0.2/s','5% 0.2/s fit','5% 0.02/s','5% 0.02/s fit',... '5% 0.002/s','5% 0.002/s fit') hold off 368 369 370 371 372 373 374 375 376 377 378 379 380 381 %% Test model at high strain rates figure title('Test model at 5% strain, at 80/s and 0.8/s strain rates') ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(solvedconstants,[t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 5 80/1000,vfun(solvedconstants,[t 5 80; lambda 5 80; dt 5 80]),'−g') plot(t 5 08/1000,vfun(solvedconstants,[t 5 08; lambda 5 08; dt 5 08]),'−k') xlim([−0.025 0.25]) legend('5% 0.2/s','5% 0.2/s fit','5% 80/s prediction','5% 0.8/s prediction') hold off 382 383 384 385 386 %% Predictions with Maikos constants n=1; m=2; 387 388 389 390 391 392 393 394 maikosconstants = % maikosconstants % maikosconstants % maikosconstants % maikosconstants % maikosconstants % maikosconstants [−40.04e−3, −4.7, 0.5282, 8, 0.3018, 150]; = [32e−3, 4.7, 0.5282, 8, 0.3018, 150]; = [16e−3, 4.7, 0.5282, 8, 0.3018, 150]; = [40e−3, 4.7, 0.5282, 8, 0.3018, 150];%using mu=G0/alpha =````````````````` [32e−3, 4.7, 0.6609, 8, 0.3777, 150]; = [200e−3, 4.7, 0.5282, 8, 0.3018, 150]; = [32e−3, 4.7]; 395 396 397 disp 'Maikos constants' hvlabels(n,m,maikosconstants); 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 figure title(['Maikos constants (2008) predictions for Fiford and Bilston'... ' 0.2/s strain rate data (2005)']) ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(maikosconstants,[t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 35 02/1000,s 35 02 ,'.g') plot(t 35 02/1000,vfun(maikosconstants,[t 35 02; lambda 35 02; dt 35 02]),'−g') plot(t 2 02/1000,s 2 02 ,'.k') plot(t 2 02/1000,vfun(maikosconstants,[t 2 02; lambda 2 02; dt 2 02]),'−k') legend('5% 0.2/s','5% 0.2/s fit','3.5% 0.2/s','3.5% 0.2/s fit','2% 0.2/s','2% 0.2/s fit') hold off 413 414 415 416 417 418 figure title(['Maikos constants (2008) predictions for Fiford and Bilston'... ' 5% max strain data (2005)']) ylabel('Stress (MPa)') xlabel('Time (s)') 78 C.3. Hyperviscoelastic fitting of Fiford data - hyperviscofitfinal.m 419 420 421 422 423 424 425 426 427 428 hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(maikosconstants,[t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 5 002/1000,s 5 002 ,'.g') plot(t 5 002/1000,vfun(maikosconstants,[t 5 002; lambda 5 002; dt 5 002]),'−g') plot(t 5 0002/1000,s 5 0002,'.k') plot(t 5 0002/1000,vfun(maikosconstants,[t 5 0002; lambda 5 0002; dt 5 0002]),'−k') legend('5% 0.2/s','5% 0.2/s fit','5% 0.02/s','5% 0.02/s fit',... '5% 0.002/s','5% 0.002/s fit') hold off 429 430 431 432 433 434 435 436 437 438 439 440 441 figure title('Test Maikos model at 5% strain, at 80/s and 0.8/s strain rates') ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(maikosconstants,[t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 5 80/1000,vfun(maikosconstants,[t 5 80; lambda 5 80; dt 5 80]),'−g') plot(t 5 08/1000,vfun(maikosconstants,[t 5 08; lambda 5 08; dt 5 08]),'−k') xlim([−0.025 0.25]) legend('5% 0.2/s','5% 0.2/s fit','5% 80/s prediction','5% 0.8/s prediction') hold off 442 443 444 445 %% Predictions with Miller constants millerconstants = [842e−6/−4.7, −4.7, 0.45, 0.5e3, 0.365, 50e3]; 446 447 448 disp 'Miller constants' hvlabels(n,m,millerconstants); 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 figure title(['Miller constants for brain (2002) predictions for Bilston'... ' 0.2/s strain rate data']) ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(millerconstants,[t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 35 02/1000,s 35 02 ,'.g') plot(t 35 02/1000,vfun(millerconstants,[t 35 02; lambda 35 02; dt 35 02]),'−g') plot(t 2 02/1000,s 2 02 ,'.k') plot(t 2 02/1000,vfun(millerconstants,[t 2 02; lambda 2 02; dt 2 02]),'−k') legend('5% 0.2/s','5% 0.2/s fit','3.5% 0.2/s','3.5% 0.2/s fit',... '2% 0.2/s','2% 0.2/s fit') hold off 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 figure title(['Miller constants for brain (2008) predictions for Fiford and Bilston'... ' 5% max strain data (2005)']) ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(millerconstants,[t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 5 002/1000,s 5 002 ,'.g') plot(t 5 002/1000,vfun(millerconstants,[t 5 002; lambda 5 002; dt 5 002]),'−g') plot(t 5 0002/1000,s 5 0002,'.k') plot(t 5 0002/1000,vfun(millerconstants,[t 5 0002; lambda 5 0002; dt 5 0002]),'−k') legend('5% 0.2/s','5% 0.2/s fit','5% 0.02/s','5% 0.02/s fit',... '5% 0.002/s','5% 0.002/s fit') hold off 79 C.4. Hyperviscoelastic algorithm validation scripts 481 482 483 484 485 486 487 488 489 490 491 492 493 figure title('Test Miller model at 5% strain, at 80/s and 0.8/s strain rates') ylabel('Stress (MPa)') xlabel('Time (s)') hold on plot(t 5 02/1000,s 5 02 ,'.r') plot(t 5 02/1000,vfun(millerconstants,[t 5 02; lambda 5 02; dt 5 02]),'−r') plot(t 5 80/1000,vfun(millerconstants,[t 5 80; lambda 5 80; dt 5 80]),'−g') plot(t 5 08/1000,vfun(millerconstants,[t 5 08; lambda 5 08; dt 5 08]),'−k') xlim([−0.025 0.25]) legend('5% 0.2/s','5% 0.2/s fit','5% 80/s prediction','5% 0.8/s prediction') hold off C.4 C.4.1 1 Hyperviscoelastic algorithm validation scripts Recreation of Greaves hyperelastic plot - greavesogdentest.m clc; clear; close all; 2 3 global n m; 4 5 6 n=1; m=0; 7 8 9 mu = 0.09/29.52; alpha = 29.52; 10 11 params = [mu; alpha]; 12 13 14 15 16 17 18 19 20 21 22 23 timestep = 1; e rate fast = 0.24e−3; % strain rate per ms e max = 0.11; t end = e max/e rate fast; t fast = 0:timestep:t end; lambda fast = e rate fast * t fast+1; dt fast = zeros(1,length(t fast)); for i=1:length(t fast)−1 dt fast(i) = t fast(i+1)−t fast(i); end dt fast(end) = dt fast(end−1); 24 25 26 vfun = @hvstress; 27 28 29 30 31 32 33 34 figure title('Plot of Bilston (1996) data fit from Greaves thesis') ylabel('Stress (MPa)') xlabel('Strain') hold on plot(lambda fast−1,vfun(params,[t fast; lambda fast; dt fast]),'−') hold off C.4.2 1 Recreation of Miller hyperviscoelastic plots - millertesthv.m clc; clear; close all; 2 3 global n m; 4 80 C.4. Hyperviscoelastic algorithm validation scripts 5 6 n=1; m=2; 7 8 timestep = 0.1; 9 10 K = 1.583; 11 12 13 14 15 16 17 alpha = −4.7; mu = 842e−6/alpha; g1 = 0.45; tau1 = 0.5e3; g2 = 0.365; tau2 = 50e3; 18 19 20 params = [mu; alpha; g1; tau1; g2; tau2]; posparams = [−mu; −alpha; g1; tau1; g2; tau2]; 21 22 % n=0; params = [mu; alpha]; % Test hyperelastic without visco 23 24 25 26 27 28 29 30 31 32 33 34 % Fast tension e rate fast = 6.4e−4; % strain rate per ms e max = (1.6−1)/K; t end = e max/e rate fast; t fastT = 0:timestep:t end; lambda fastT = K*(e rate fast *t fastT)+1; dt fastT = zeros(1,length(t fastT)); for i=1:length(t fastT)−1 dt fastT(i) = t fastT(i+1)−t fastT(i); end dt fastT(end) = dt fastT(end−1); 35 36 37 38 39 40 41 42 43 44 45 46 % Fast compression e rate fast = 6.4e−4; % strain rate per ms e min = 0.6; t end = (1−e min)/e rate fast; t fastC = 0:timestep:t end; lambda fastC = 1−e rate fast * t fastC; dt fastC = zeros(1,length(t fastC)); for i=1:length(t fastC)−1 dt fastC(i) = t fastC(i+1)−t fastC(i); end dt fastC(end) = dt fastC(end−1); 47 48 49 50 51 52 53 54 55 56 57 58 59 % Slow tension e rate slow = 6.4e−6; % strain rate per ms e max = (1.3−1)/K; t end = e max/e rate slow; t slowT = 0:100*timestep:t end; lambda slowT = K*(e rate slow *t slowT)+1; dt slowT = zeros(1,length(t slowT)); for i=1:length(t slowT)−1 dt slowT(i) = t slowT(i+1)−t slowT(i); end dt slowT(end) = dt slowT(end−1); 60 61 62 63 64 65 66 % Slow compression e rate slow = 6.4e−6; % strain rate per ms e min = 0.7; t end = (1−e min)/e rate slow; t slowC = 0:100*timestep:t end; lambda slowC = 1−e rate slow * t slowC; 81 C.4. Hyperviscoelastic algorithm validation scripts 67 68 69 70 71 dt slowC = zeros(1,length(t slowC)); for i=1:length(t slowC)−1 dt slowC(i) = t slowC(i+1)−t slowC(i); end dt slowC(end) = dt slowC(end−1); 72 73 74 75 76 77 78 79 80 81 82 83 % Very Slow compression e rate veryslow = 6.4e−9; % strain rate per ms e min = 0.7; t end = (1−e min)/e rate veryslow; t veryslowC = 0:1e5*timestep:t end; lambda veryslowC = 1−e rate veryslow *t veryslowC; dt veryslowC = zeros(1,length(t veryslowC)); for i=1:length(t veryslowC)−1 dt veryslowC(i) = t veryslowC(i+1)−t veryslowC(i); end dt veryslowC(end) = dt veryslowC(end−1); 84 85 86 87 % vfun = @hvtensilestress; vfun = @hvstress; 88 89 hvlabels(n,m,params) 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 figure title('Plot of Miller and Chinzei (2002) data fit') ylabel('Stress (Pa)') xlabel('Stretch ratio, lambda') hold on grid plot(lambda fastT,(1e6)*vfun(params,[t fastT; lambda fastT; dt fastT]),... lambda fastC,(1e6)*vfun(params,[t fastC; lambda fastC; dt fastC]),'g') plot(lambda slowT,(1e6)*vfun(params,[t slowT; lambda slowT; dt slowT]),'−−',... lambda slowC,(1e6)*vfun(params,[t slowC; lambda slowC; dt slowC]),'−−g') plot(lambda veryslowC,(1e6)*vfun(params,[t veryslowC; lambda veryslowC;... dt veryslowC]),'−.g') plot(lambda fastT,(1e6)*vfun(posparams,[t fastT; lambda fastT; dt fastT]),'−k',... lambda fastC,(1e6)*vfun(posparams,[t fastC; lambda fastC; dt fastC]),'−r') plot(lambda slowT,(1e6)*vfun(posparams,[t slowT; lambda slowT; dt slowT]),'−−k',... lambda slowC,(1e6)*vfun(posparams,[t slowC; lambda slowC; dt slowC]),'−−r') plot(lambda veryslowC,(1e6)*vfun(posparams,[t veryslowC; lambda veryslowC;... dt veryslowC]),'−.r') legend('Fast Tension {\alpha = −4.7}','Fast Compression {\alpha = −4.7}',... 'Slow Tension {\alpha = −4.7}','Slow Compression {\alpha = −4.7}',... 'Very Slow Compression {\alpha = −4.7}','Fast Tension {\alpha = +4.7}',... 'Fast Compression {\alpha = +4.7}','Slow Tension {\alpha = +4.7}',... 'Slow Compression {\alpha = +4.7}','Very Slow Compression {\alpha = +4.7}',... 'Location','Best') hold off C.4.3 1 Recreation of Snedeker hyperviscoelastic plots - snedeckertest.m clc; clear; close all; 2 3 global n m; 4 5 6 n=2; m=4; 7 8 9 % Mu divided by two because Snedecker's use of Ogden strain energy did not % include leading 2 multiplier 82 C.5. Material model comparison 10 11 12 13 14 15 16 17 18 19 20 21 mu1 = 0.2/2; alpha1 = 15; mu2 = 4.2/2; alpha2 = 7.5; g1 = 0.12; tau1 = 5; g2 = 0.4; tau2 = 1; g3 = 0.08; tau3 = 0.1e3; g4 = 0.13; tau4 = 5e3; 22 23 params = [mu1; alpha1; mu2; alpha2; g1; tau1; g2; tau2; g3; tau3; g4; tau4]; 24 25 26 27 28 29 30 31 32 33 34 35 timestep = 1e−3; e rate fast = 250e−3; % strain rate per ms e max = 0.325; t end = e max/e rate fast; t fast = 0:timestep:t end; lambda fast = e rate fast * t fast+1; dt fast = zeros(1,length(t fast)); for i=1:length(t fast)−1 dt fast(i) = t fast(i+1)−t fast(i); end dt fast(end) = dt fast(end−1); 36 37 38 39 40 41 42 43 44 45 46 47 timestep = 100; e rate slow = 0.005e−3; e max = 0.325; t end = e max/e rate slow; t slow = 0:timestep:t end; lambda slow = e rate slow * t slow+1; dt slow = zeros(1,length(t slow)); for i=1:length(t slow)−1 dt slow(i) = t slow(i+1)−t slow(i); end dt slow(end) = dt slow(end−1); 48 49 50 vfun = @hvstress; 51 52 53 54 55 56 57 58 59 60 figure title('Plot of Snedecker (2005) data fit') ylabel('Stress (MPa)') xlabel('Strain') hold on plot(lambda fast−1,vfun(params,[t fast; lambda fast; dt fast]),'−') plot(lambda slow−1,vfun(params,[t slow; lambda slow; dt slow]),'−−') legend('Fast','Slow') hold off C.5 Material model comparison 1 2 clc; clear; close all; 3 4 global m n t relax; 5 6 % Units: Stress (MPa) 83 C.5. Material model comparison 7 % Time (ms) 8 9 10 n = 1; m = 4; 11 12 13 % Relaxation time (ms) t relax = 1800e3; 14 15 16 % Approx. ratio of steady state to initial stress at 5%, 0.2/s relax ratio = 0.045/0.076; 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 % Set G0 equal to initial stress over 5% strain at 0.2/s rate G0 = 0.076/0.05; for i=1:n % Set initial mu i params(2*i−1) = 200e−3; % Set initial alpha i params(2*i) = 2*i; end for i=1:m % Set initial g i params(2*i−1+2*m) = 0.9; % Set initial tau i params(2*i+2*m) = 10ˆ(i−1); end 32 33 34 35 36 37 % Load Fiford and Bilston data % t in ms, s and sSD in MPa load fbdata; 38 39 40 41 42 % 5% max strain, 0.2/s e max = 0.05; e rate = 0.2e−3; % (per ms) t1 = e max/e rate; 43 44 45 46 47 48 49 50 51 52 t = t 5 02; lambda 5 02 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); % true strain = log(lambda 5 02); % s 5 02 = s 5 02./(1+true strain); % Convert true stress to nominal stress dt 5 02 = zeros(1,length(t)); for i=1:length(t)−1 dt 5 02(i) = t(i+1)−t(i); end dt 5 02(end) = dt 5 02(end−1); 53 54 55 56 57 58 59 60 61 62 63 64 65 % 5% max strain, 0.02/s e max = 0.05; e rate = 0.02e−3; % (per ms) t1 = e max/e rate; t = t 5 002; lambda 5 002 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 5 002 = zeros(1,length(t)); for i=1:length(t)−1 dt 5 002(i) = t(i+1)−t(i); end dt 5 002(end) = dt 5 002(end−1); 66 67 68 % 5% max strain, 0.002/s e max = 0.05; 84 C.5. Material model comparison 69 70 71 72 73 74 75 76 77 e rate = 0.002e−3; % (per ms) t1 = e max/e rate; t = t 5 0002; lambda 5 0002 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 5 0002 = zeros(1,length(t)); for i=1:length(t)−1 dt 5 0002(i) = t(i+1)−t(i); end dt 5 0002(end) = dt 5 0002(end−1); 78 79 80 81 82 83 84 85 86 87 88 89 % 3.5% max strain, 0.2/s e max = 0.035; e rate = 0.2e−3; % (per ms) t1 = e max/e rate; t = t 35 02; lambda 35 02 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 35 02 = zeros(1,length(t)); for i=1:length(t)−1 dt 35 02(i) = t(i+1)−t(i); end dt 35 02(end) = dt 35 02(end−1); 90 91 92 93 94 95 96 97 98 99 100 101 % 3.5% max strain, 0.002/s e max = 0.035; e rate = 0.002e−3; % (per ms) t1 = e max/e rate; t = t 35 0002; lambda 35 0002 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 35 0002 = zeros(1,length(t)); for i=1:length(t)−1 dt 35 0002(i) = t(i+1)−t(i); end dt 35 0002(end) = dt 35 0002(end−1); 102 103 104 105 106 107 108 109 110 111 112 113 % 2% max strain, 0.2/s e max = 0.02; e rate = 0.2e−3; % (per ms) t1 = e max/e rate; t = t 2 02; lambda 2 02 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 2 02 = zeros(1,length(t)); for i=1:length(t)−1 dt 2 02(i) = t(i+1)−t(i); end dt 2 02(end) = dt 2 02(end−1); 114 115 116 117 118 119 120 121 122 123 124 125 % 2% max strain, 0.002/s e max = 0.02; e rate = 0.002e−3; % (per ms) t1 = e max/e rate; t = t 2 0002; lambda 2 0002 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 2 0002 = zeros(1,length(t)); for i=1:length(t)−1 dt 2 0002(i) = t(i+1)−t(i); end dt 2 0002(end) = dt 2 0002(end−1); 126 127 128 %% Describe strain history 129 130 timestep = 0.001; 85 C.5. Material model comparison 131 132 133 134 135 136 137 138 139 140 141 142 % 60% max tensile strain, 100/s e max = 0.6; e rate = 100e−3; t1 = e max/e rate; t 60T 100 = [0:timestep:t1]; t = t 60T 100; lambda 60T 100 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 60T 100 = zeros(1,length(t)); for i=1:length(t)−1 dt 60T 100(i) = t(i+1)−t(i); end dt 60T 100(end) = dt 60T 100(end−1); 143 144 145 146 147 148 149 150 151 152 153 154 155 % 40% max compressive strain, 100/s e max = 0.4; e rate = 100e−3; t1 = e max/e rate; t 40C 100 = [0:timestep:t1]; t = t 40C 100; lambda 40C 100 = (t<=t1).*(1−e rate *t)+(t>t1).*(e rate *t1+1); dt 40C 100 = zeros(1,length(t)); for i=1:length(t)−1 dt 40C 100(i) = t(i+1)−t(i); end dt 40C 100(end) = dt 40C 100(end−1); 156 157 158 159 160 161 162 163 164 165 166 167 168 169 timestep = timestep*100; % 60% max tensile strain, 1/s e max = 0.6; e rate = 1e−3; t1 = e max/e rate; t 60T 1 = [0:timestep:t1]; t = t 60T 1; lambda 60T 1 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 60T 1 = zeros(1,length(t)); for i=1:length(t)−1 dt 60T 1(i) = t(i+1)−t(i); end dt 60T 1(end) = dt 60T 1(end−1); 170 171 172 173 174 175 176 177 178 179 180 181 182 % 40% max tensile strain, 1/s e max = 0.4; e rate = 1e−3; t1 = e max/e rate; t 40C 1 = [0:timestep:t1]; t = t 40C 1; lambda 40C 1 = (t<=t1).*(1−e rate *t)+(t>t1).*(e rate *t1+1); dt 40C 1 = zeros(1,length(t)); for i=1:length(t)−1 dt 40C 1(i) = t(i+1)−t(i); end dt 40C 1(end) = dt 40C 1(end−1); 183 184 185 186 187 188 189 190 191 192 timestep = timestep*5; % 60% max tensile strain, 0.2/s e max = 0.6; e rate = 0.2e−3; t1 = e max/e rate; t 60T 02 = [0:timestep:t1]; t = t 60T 02; lambda 60T 02 = (t<=t1).*(e rate *t+1)+(t>t1).*(e rate *t1+1); dt 60T 02 = zeros(1,length(t)); 86 C.5. Material model comparison 193 194 195 196 for i=1:length(t)−1 dt 60T 02(i) = t(i+1)−t(i); end dt 60T 02(end) = dt 60T 02(end−1); 197 198 199 200 201 202 203 204 205 206 207 208 209 % 40% max tensile strain, 0.2/s e max = 0.4; e rate = 0.2e−3; t1 = e max/e rate; t 40C 02 = [0:timestep:t1]; t = t 40C 02; lambda 40C 02 = (t<=t1).*(1−e rate *t)+(t>t1).*(e rate *t1+1); dt 40C 02 = zeros(1,length(t)); for i=1:length(t)−1 dt 40C 02(i) = t(i+1)−t(i); end dt 40C 02(end) = dt 40C 02(end−1); 210 211 vfun = @hvstress; 212 213 214 215 216 217 218 219 220 %% Master plot of experimental data and model predictions maikosconstants = [40.04e−3, 4.7, 0.5282, 8, 0.3018, 150]; negmaikosconstants = [−40.04e−3, −4.7, 0.5282, 8, 0.3018, 150]; % negmaikosconstants = [−40.04e−3, −4.7, 0.5282, 1, 0.3018, 150]; % test with faster relaxation 1ms time constant millerconstants = [842e−6/−4.7, −4.7, 0.45, 0.5e3, 0.365, 50e3]; fifordconstants = [7.7e−3, 49.9342, 0.5282, 8, 0.094, 256.0977, 0.051,... 38183, 0.0488, 456240]; 221 222 223 224 225 226 227 228 229 disp maikosconstants hvlabels(1,2,maikosconstants); disp negmaikosconstants hvlabels(1,2,negmaikosconstants); disp millerconstants hvlabels(1,2,millerconstants); disp fifordconstants hvlabels(1,4,fifordconstants); 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 figure title('Master stress−strain plot of experimental data and model predictions') ylabel('Stress, {\sigma} (MPa)') xlabel('Stretch Ratio, {\lambda}') ylim([−1.5 0.5]) xlim([0.6 1.4]) % grid hold on currentPlot = []; % Fast strain rate 100/s n=1; m=2; currentPlot = [currentPlot;plot(lambda 60T 100,vfun(maikosconstants,... [t 60T 100; lambda 60T 100; dt 60T 100]),'r')]; plot(lambda 40C 100,vfun(maikosconstants,[t 40C 100; lambda 40C 100; dt 40C 100]),'r') currentPlot = [currentPlot;plot(lambda 60T 100,vfun(negmaikosconstants,... [t 60T 100; lambda 60T 100; dt 60T 100]),'g')]; plot(lambda 40C 100,vfun(negmaikosconstants,[t 40C 100; lambda 40C 100; dt 40C 100]),'g') currentPlot = [currentPlot;plot(lambda 60T 100,vfun(millerconstants,... [t 60T 100; lambda 60T 100; dt 60T 100]),'b')]; plot(lambda 40C 100,vfun(millerconstants,[t 40C 100; lambda 40C 100; dt 40C 100]),'b') n=1; m=4; currentPlot = [currentPlot;plot(lambda 60T 100,vfun(fifordconstants,... [t 60T 100; lambda 60T 100; dt 60T 100]),'k')]; 87 C.5. Material model comparison 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 plot(lambda 40C 100,vfun(fifordconstants,[t 40C 100; lambda 40C 100; dt 40C 100]),'k') % Slow strain rate 1/s n=1; m=2; currentPlot = [currentPlot;plot(lambda 60T 1,vfun(maikosconstants,... [t 60T 1; lambda 60T 1; dt 60T 1]),'−−r')]; plot(lambda 40C 1,vfun(maikosconstants,[t 40C 1; lambda 40C 1; dt 40C 1]),'−−r') currentPlot = [currentPlot;plot(lambda 60T 1,vfun(negmaikosconstants,... [t 60T 1; lambda 60T 1; dt 60T 1]),'−−g')]; plot(lambda 40C 1,vfun(negmaikosconstants,[t 40C 1; lambda 40C 1; dt 40C 1]),'−−g') currentPlot = [currentPlot;plot(lambda 60T 1,vfun(millerconstants,... [t 60T 1; lambda 60T 1; dt 60T 1]),'−−b')]; plot(lambda 40C 1,vfun(millerconstants,[t 40C 1; lambda 40C 1; dt 40C 1]),'−−b') n=1; m=4; currentPlot = [currentPlot;plot(lambda 60T 1,vfun(fifordconstants,... [t 60T 1; lambda 60T 1; dt 60T 1]),'−−k')]; plot(lambda 40C 1,vfun(fifordconstants,[t 40C 1; lambda 40C 1; dt 40C 1]),'−−k') hold off legend(currentPlot,'100/s − Maikos parameters','100/s − Negative Maikos parameters',... '100/s − Miller parameters (brain)','100/s − Fit to Fiford and Bilston data',... '1/s − Maikos parameters','1/s − Negative Maikos parameters',... '1/s − Miller parameters (brain)','1/s − Fit to Fiford and Bilston data',... 'Location','SouthEast') saveas(gcf,['masterplot.eps'],'psc2') 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 figure title('Master stress−strain plot of experimental data and model predictions (zoom)') ylabel('Stress, {\sigma} (MPa)') xlabel('Stretch Ratio, {\lambda}') ylim([−0.1 0.1]) xlim([0.9 1.1]) % grid hold on currentPlot = []; % Fast strain rate 100/s n=1; m=2; currentPlot = [currentPlot;plot(lambda 60T 100,vfun(maikosconstants,... [t 60T 100; lambda 60T 100; dt 60T 100]),'r')]; plot(lambda 40C 100,vfun(maikosconstants,[t 40C 100; lambda 40C 100; dt 40C 100]),'r') currentPlot = [currentPlot;plot(lambda 60T 100,vfun(negmaikosconstants,... [t 60T 100; lambda 60T 100; dt 60T 100]),'g')]; plot(lambda 40C 100,vfun(negmaikosconstants,[t 40C 100; lambda 40C 100; dt 40C 100]),'g') currentPlot = [currentPlot;plot(lambda 60T 100,vfun(millerconstants,... [t 60T 100; lambda 60T 100; dt 60T 100]),'b')]; plot(lambda 40C 100,vfun(millerconstants,[t 40C 100; lambda 40C 100; dt 40C 100]),'b') n=1; m=4; currentPlot = [currentPlot;plot(lambda 60T 100,vfun(fifordconstants,... [t 60T 100; lambda 60T 100; dt 60T 100]),'k')]; plot(lambda 40C 100,vfun(fifordconstants,[t 40C 100; lambda 40C 100; dt 40C 100]),'k') % Slow strain rate 1/s n=1; m=2; currentPlot = [currentPlot;plot(lambda 60T 1,vfun(maikosconstants,... [t 60T 1; lambda 60T 1; dt 60T 1]),'−−r')]; plot(lambda 40C 1,vfun(maikosconstants,[t 40C 1; lambda 40C 1; dt 40C 1]),'−−r') currentPlot = [currentPlot;plot(lambda 60T 1,vfun(negmaikosconstants,... [t 60T 1; lambda 60T 1; dt 60T 1]),'−−g')]; plot(lambda 40C 1,vfun(negmaikosconstants,[t 40C 1; lambda 40C 1; dt 40C 1]),'−−g') currentPlot = [currentPlot;plot(lambda 60T 1,vfun(millerconstants,... [t 60T 1; lambda 60T 1; dt 60T 1]),'−−b')]; plot(lambda 40C 1,vfun(millerconstants,[t 40C 1; lambda 40C 1; dt 40C 1]),'−−b') n=1; m=4; currentPlot = [currentPlot;plot(lambda 60T 1,vfun(fifordconstants,... 88 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 [t 60T 1; lambda 60T 1; dt 60T 1]),'−−k')]; plot(lambda 40C 1,vfun(fifordconstants,[t 40C 1; lambda 40C 1; dt 40C 1]),'−−k') % Fiford and Bilston data, 5% strain 0.2/s currentPlot = [currentPlot;plot(lambda 5 02,s 5 02 ,'.−','MarkerSize',5)]; currentPlot = [currentPlot;plot(lambda 60T 02,vfun(fifordconstants,... [t 60T 02; lambda 60T 02; dt 60T 02]),':k')]; plot(lambda 40C 02,vfun(fifordconstants,[t 40C 02; lambda 40C 02; dt 40C 02]),':k') hold off legend(currentPlot,'100/s − Maikos parameters','100/s − Negative Maikos parameters',... '100/s − Miller parameters (brain)','100/s − Fit to Fiford and Bilston data',... '1/s − Maikos parameters','1/s − Negative Maikos parameters',... '1/s − Miller parameters (brain)','1/s − Fit to Fiford and Bilston data',... '0.2/s − Fiford and Bilston data','0.2/s − Fit to Fiford and Bilston data',... 'Location','SouthEast') saveas(gcf,['masterplot zoom.eps'],'psc2') C.5.1 Model parameters - output from hvlabels.m maikosconstants 'mu1 (MPa)' [ 0.0400] 'G0 (MPa)' [ 0.1882] 'alpha1' 'g1' 'tau1 (ms)' 'g2' 'tau2 (ms)' [4.7000] [ 0.5282] [ 8] [ 0.3018] [ 150] [] 'G1 (MPa)' [] 'G2 (MPa)' [] [] [ 0.0994] [] [ 0.0568] [] negmaikosconstant s 'mu1 (MPa)' 'alpha1' 'g1' 'tau1 (ms)' 'g2' 'tau2 (ms)' [ -0.0400] [-4.7000] [ 0.5282] [ 8] [ 0.3018] [ 150] 'G0 (MPa)' [] 'G1 (MPa)' [] 'G2 (MPa)' [] [ 0.1882] [] [ 0.0994] [] [ 0.0568] [] millerconstants 'mu1 (MPa)' 'alpha1' 'g1' [-1.7915e-004] [-4.7000 ] [ 0.4 'G0 (MPa)' [ 'G1 (MPa) [ 8.4200e-004 ] [ [3.7890e- 'tau1 ( ms)' 'g2' 'tau 500] [ 500] [ 0.3650] [ ' [] 'G2 (M Pa)' 004] [] [3.073 3e-004] 2 (ms)' 50000] [] [] ﬁfordconstants 'mu1 (MPa)' 'alpha1' 'g1' 'tau1 (ms)' 'g2' 'tau2 (ms)' 'g3' 'tau3 (ms)' 'g4' 'tau4 (ms)' [ 0.0077] [49.9342] [ 0.5282] [ 8] [ 0.0940][ 256.0977] [ 0.0510] [ 38183] [ 0.0488] [ 456240] 'G0 (MPa)' [] 'G1 (MPa)' [] 'G2 (MPa)' [] 'G3 (MPa)' [] 'G4 (MPa)' [] [ 0.3845] [] [ 0.2031] [] [ 0.0361] [] [ 0.0196] [] [ 0.0188] [] C.6 C.6.1 1 FE data extraction and linear regression of tissue damage versus maximum principal strain Reading of exported FE data from .csv files and saving to .mat zonesliceread.m %% Read stress and strain simulation results from .csv and save to .mat 2 3 clc; clear; 4 5 6 7 8 %% Specify region to process ie. 'negMaikos WM DC CT', and folder where %% data reside region = 'Maikos WM LC DL'; folder = '\Dislocation\Maikos\'; 9 89 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 for i=−5:5 if i>0 file name = [region ' +' num2str(i) 'mm.csv']; data name = [region ' data pos' num2str(abs(i))]; labels name = [region ' labels pos' num2str(abs(i))]; else if i<0 file name = [region ' ' num2str(i) 'mm.csv']; data name = [region ' data neg' num2str(abs(i))]; labels name = [region ' labels neg' num2str(abs(i))]; else file name = [region ' ' num2str(i) 'mm.csv']; data name = [region ' data ' num2str(i)]; labels name = [region ' labels ' num2str(i)]; end end eval(['[' data name ',' labels name ']' '= xlsread([cd folder file name]);']); disp(['reading data: ' region ' ' num2str(i) 'mm']) end clear file name data name labels name i; 29 30 save([region '.mat']) C.6.2 1 2 3 Calculation of mean peak values and SD and saving to .mat mechanismFEdatasave.m %% Load stress and strain simulation results from all zone slices, %% calculate average peak values, and save to mechanismFEdata.mat %% and mechanismFEdataSD.mat 4 5 clc; clear; close all; 6 7 %% Calculate mean of peaks 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 % Contusion with Maikos properties [Maikos WM DC CT strain1,Maikos WM DC CT strain2,Maikos WM DC CT strain3,... Maikos WM DC CT stress1,Maikos WM DC CT stress2,... Maikos WM DC CT stress3] = avgpeak('Maikos WM DC CT'); [Maikos WM LC CT strain1,Maikos WM LC CT strain2,Maikos WM LC CT strain3,... Maikos WM LC CT stress1,Maikos WM LC CT stress2,... Maikos WM LC CT stress3] = avgpeak('Maikos WM LC CT'); [Maikos WM VM CT strain1,Maikos WM VM CT strain2,Maikos WM VM CT strain3,... Maikos WM VM CT stress1,Maikos WM VM CT stress2,... Maikos WM VM CT stress3] = avgpeak('Maikos WM VM CT'); [Maikos WM VL CT strain1,Maikos WM VL CT strain2,Maikos WM VL CT strain3,... Maikos WM VL CT stress1,Maikos WM VL CT stress2,... Maikos WM VL CT stress3] = avgpeak('Maikos WM VL CT'); [Maikos GM CT strain1,Maikos GM CT strain2,Maikos GM CT strain3,... Maikos GM CT stress1,Maikos GM CT stress2,... Maikos GM CT stress3] = avgpeak('Maikos GM CT'); 25 26 27 28 29 30 31 32 33 34 35 36 % Contusion with negative Maikos properties [negMaikos WM DC CT strain1,negMaikos WM DC CT strain2,negMaikos negMaikos WM DC CT stress1,negMaikos WM DC CT stress2,... negMaikos WM DC CT stress3] = avgpeak('negMaikos WM DC CT'); [negMaikos WM LC CT strain1,negMaikos WM LC CT strain2,negMaikos negMaikos WM LC CT stress1,negMaikos WM LC CT stress2,... negMaikos WM LC CT stress3] = avgpeak('negMaikos WM LC CT'); [negMaikos WM VM CT strain1,negMaikos WM VM CT strain2,negMaikos negMaikos WM VM CT stress1,negMaikos WM VM CT stress2,... negMaikos WM VM CT stress3] = avgpeak('negMaikos WM VM CT'); [negMaikos WM VL CT strain1,negMaikos WM VL CT strain2,negMaikos WM DC CT strain3,... WM LC CT strain3,... WM VM CT strain3,... WM VL CT strain3,... 90 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 37 38 39 40 41 negMaikos WM VL CT stress1,negMaikos WM VL CT stress2,... negMaikos WM VL CT stress3] = avgpeak('negMaikos WM VL CT'); [negMaikos GM CT strain1,negMaikos GM CT strain2,negMaikos GM CT strain3,... negMaikos GM CT stress1,negMaikos GM CT stress2,... negMaikos GM CT stress3] = avgpeak('negMaikos GM CT'); 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 % Dislocation with Maikos properties [Maikos WM DC DL strain1,Maikos WM DC DL strain2,Maikos WM DC DL strain3,... Maikos WM DC DL stress1,Maikos WM DC DL stress2,... Maikos WM DC DL stress3] = avgpeak('Maikos WM DC DL'); [Maikos WM LC DL strain1,Maikos WM LC DL strain2,Maikos WM LC DL strain3,... Maikos WM LC DL stress1,Maikos WM LC DL stress2,... Maikos WM LC DL stress3] = avgpeak('Maikos WM LC DL'); [Maikos WM VM DL strain1,Maikos WM VM DL strain2,Maikos WM VM DL strain3,... Maikos WM VM DL stress1,Maikos WM VM DL stress2,... Maikos WM VM DL stress3] = avgpeak('Maikos WM VM DL'); [Maikos WM VL DL strain1,Maikos WM VL DL strain2,Maikos WM VL DL strain3,... Maikos WM VL DL stress1,Maikos WM VL DL stress2,... Maikos WM VL DL stress3] = avgpeak('Maikos WM VL DL'); [Maikos GM DL strain1,Maikos GM DL strain2,Maikos GM DL strain3,... Maikos GM DL stress1,Maikos GM DL stress2,... Maikos GM DL stress3] = avgpeak('Maikos GM DL'); 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 % Dislocation with negative Maikos properties [negMaikos WM DC DL strain1,negMaikos WM DC DL strain2,negMaikos WM DC DL strain3,... negMaikos WM DC DL stress1,negMaikos WM DC DL stress2,... negMaikos WM DC DL stress3] = avgpeak('negMaikos WM DC DL'); [negMaikos WM LC DL strain1,negMaikos WM LC DL strain2,negMaikos WM LC DL strain3,... negMaikos WM LC DL stress1,negMaikos WM LC DL stress2,... negMaikos WM LC DL stress3] = avgpeak('negMaikos WM LC DL'); [negMaikos WM VM DL strain1,negMaikos WM VM DL strain2,negMaikos WM VM DL strain3,... negMaikos WM VM DL stress1,negMaikos WM VM DL stress2,... negMaikos WM VM DL stress3] = avgpeak('negMaikos WM VM DL'); [negMaikos WM VL DL strain1,negMaikos WM VL DL strain2,negMaikos WM VL DL strain3,... negMaikos WM VL DL stress1,negMaikos WM VL DL stress2,... negMaikos WM VL DL stress3] = avgpeak('negMaikos WM VL DL'); [negMaikos GM DL strain1,negMaikos GM DL strain2,negMaikos GM DL strain3,... negMaikos GM DL stress1,negMaikos GM DL stress2,... negMaikos GM DL stress3] = avgpeak('negMaikos GM DL'); 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 % Distraction with Maikos properties [Maikos WM DC DT strain1,Maikos WM DC DT strain2,Maikos WM DC DT strain3,... Maikos WM DC DT stress1,Maikos WM DC DT stress2,... Maikos WM DC DT stress3] = avgpeak('Maikos WM DC DT'); [Maikos WM LC DT strain1,Maikos WM LC DT strain2,Maikos WM LC DT strain3,... Maikos WM LC DT stress1,Maikos WM LC DT stress2,... Maikos WM LC DT stress3] = avgpeak('Maikos WM LC DT'); [Maikos WM VM DT strain1,Maikos WM VM DT strain2,Maikos WM VM DT strain3,... Maikos WM VM DT stress1,Maikos WM VM DT stress2,... Maikos WM VM DT stress3] = avgpeak('Maikos WM VM DT'); [Maikos WM VL DT strain1,Maikos WM VL DT strain2,Maikos WM VL DT strain3,... Maikos WM VL DT stress1,Maikos WM VL DT stress2,... Maikos WM VL DT stress3] = avgpeak('Maikos WM VL DT'); [Maikos GM DT strain1,Maikos GM DT strain2,Maikos GM DT strain3,... Maikos GM DT stress1,Maikos GM DT stress2,... Maikos GM DT stress3] = avgpeak('Maikos GM DT'); 93 94 95 96 97 98 % Distraction with negative Maikos properties [negMaikos WM DC DT strain1,negMaikos WM DC DT strain2,negMaikos WM DC DT strain3,... negMaikos WM DC DT stress1,negMaikos WM DC DT stress2,... negMaikos WM DC DT stress3] = avgpeak('negMaikos WM DC DT'); [negMaikos WM LC DT strain1,negMaikos WM LC DT strain2,negMaikos WM LC DT strain3,... 91 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 99 100 101 102 103 104 105 106 107 108 109 negMaikos WM LC DT stress1,negMaikos WM LC DT stress2,... negMaikos WM LC DT stress3] = avgpeak('negMaikos WM LC DT'); [negMaikos WM VM DT strain1,negMaikos WM VM DT strain2,negMaikos WM VM DT strain3,... negMaikos WM VM DT stress1,negMaikos WM VM DT stress2,... negMaikos WM VM DT stress3] = avgpeak('negMaikos WM VM DT'); [negMaikos WM VL DT strain1,negMaikos WM VL DT strain2,negMaikos WM VL DT strain3,... negMaikos WM VL DT stress1,negMaikos WM VL DT stress2,... negMaikos WM VL DT stress3] = avgpeak('negMaikos WM VL DT'); [negMaikos GM DT strain1,negMaikos GM DT strain2,negMaikos GM DT strain3,... negMaikos GM DT stress1,negMaikos GM DT stress2,... negMaikos GM DT stress3] = avgpeak('negMaikos GM DT'); 110 111 112 113 % Save as .mat file save('mechanismFEdata.mat') clear; 114 115 %% Calculate S.D. of peaks 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 % Contusion with Maikos properties [Maikos WM DC CT strain1SD,Maikos WM DC CT strain2SD,Maikos WM DC CT strain3SD,... Maikos WM DC CT stress1SD,Maikos WM DC CT stress2SD,... Maikos WM DC CT stress3SD] = sdpeak('Maikos WM DC CT'); [Maikos WM LC CT strain1SD,Maikos WM LC CT strain2SD,Maikos WM LC CT strain3SD,... Maikos WM LC CT stress1SD,Maikos WM LC CT stress2SD,... Maikos WM LC CT stress3SD] = sdpeak('Maikos WM LC CT'); [Maikos WM VM CT strain1SD,Maikos WM VM CT strain2SD,Maikos WM VM CT strain3SD,... Maikos WM VM CT stress1SD,Maikos WM VM CT stress2SD,... Maikos WM VM CT stress3SD] = sdpeak('Maikos WM VM CT'); [Maikos WM VL CT strain1SD,Maikos WM VL CT strain2SD,Maikos WM VL CT strain3SD,... Maikos WM VL CT stress1SD,Maikos WM VL CT stress2SD,... Maikos WM VL CT stress3SD] = sdpeak('Maikos WM VL CT'); [Maikos GM CT strain1SD,Maikos GM CT strain2SD,Maikos GM CT strain3SD,... Maikos GM CT stress1SD,Maikos GM CT stress2SD,... Maikos GM CT stress3SD] = sdpeak('Maikos GM CT'); 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 % Contusion with negative Maikos properties [negMaikos WM DC CT strain1SD,negMaikos WM DC CT strain2SD,negMaikos WM DC CT strain3SD,... negMaikos WM DC CT stress1SD,negMaikos WM DC CT stress2SD,... negMaikos WM DC CT stress3SD] = sdpeak('negMaikos WM DC CT'); [negMaikos WM LC CT strain1SD,negMaikos WM LC CT strain2SD,negMaikos WM LC CT strain3SD,... negMaikos WM LC CT stress1SD,negMaikos WM LC CT stress2SD,... negMaikos WM LC CT stress3SD] = sdpeak('negMaikos WM LC CT'); [negMaikos WM VM CT strain1SD,negMaikos WM VM CT strain2SD,negMaikos WM VM CT strain3SD,... negMaikos WM VM CT stress1SD,negMaikos WM VM CT stress2SD,... negMaikos WM VM CT stress3SD] = sdpeak('negMaikos WM VM CT'); [negMaikos WM VL CT strain1SD,negMaikos WM VL CT strain2SD,negMaikos WM VL CT strain3SD,... negMaikos WM VL CT stress1SD,negMaikos WM VL CT stress2SD,... negMaikos WM VL CT stress3SD] = sdpeak('negMaikos WM VL CT'); [negMaikos GM CT strain1SD,negMaikos GM CT strain2SD,negMaikos GM CT strain3SD,... negMaikos GM CT stress1SD,negMaikos GM CT stress2SD,... negMaikos GM CT stress3SD] = sdpeak('negMaikos GM CT'); 150 151 152 153 154 155 156 157 158 159 160 % Dislocation with Maikos properties [Maikos WM DC DL strain1SD,Maikos WM DC DL strain2SD,Maikos WM DC DL strain3SD,... Maikos WM DC DL stress1SD,Maikos WM DC DL stress2SD,... Maikos WM DC DL stress3SD] = sdpeak('Maikos WM DC DL'); [Maikos WM LC DL strain1SD,Maikos WM LC DL strain2SD,Maikos WM LC DL strain3SD,... Maikos WM LC DL stress1SD,Maikos WM LC DL stress2SD,... Maikos WM LC DL stress3SD] = sdpeak('Maikos WM LC DL'); [Maikos WM VM DL strain1SD,Maikos WM VM DL strain2SD,Maikos WM VM DL strain3SD,... Maikos WM VM DL stress1SD,Maikos WM VM DL stress2SD,... Maikos WM VM DL stress3SD] = sdpeak('Maikos WM VM DL'); 92 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 161 162 163 164 165 166 [Maikos WM VL DL strain1SD,Maikos WM VL DL strain2SD,Maikos WM VL DL strain3SD,... Maikos WM VL DL stress1SD,Maikos WM VL DL stress2SD,... Maikos WM VL DL stress3SD] = sdpeak('Maikos WM VL DL'); [Maikos GM DL strain1SD,Maikos GM DL strain2SD,Maikos GM DL strain3SD,... Maikos GM DL stress1SD,Maikos GM DL stress2SD,... Maikos GM DL stress3SD] = sdpeak('Maikos GM DL'); 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 % Dislocation with negative Maikos properties [negMaikos WM DC DL strain1SD,negMaikos WM DC DL strain2SD,negMaikos WM DC DL strain3SD,... negMaikos WM DC DL stress1SD,negMaikos WM DC DL stress2SD,... negMaikos WM DC DL stress3SD] = sdpeak('negMaikos WM DC DL'); [negMaikos WM LC DL strain1SD,negMaikos WM LC DL strain2SD,negMaikos WM LC DL strain3SD,... negMaikos WM LC DL stress1SD,negMaikos WM LC DL stress2SD,... negMaikos WM LC DL stress3SD] = sdpeak('negMaikos WM LC DL'); [negMaikos WM VM DL strain1SD,negMaikos WM VM DL strain2SD,negMaikos WM VM DL strain3SD,... negMaikos WM VM DL stress1SD,negMaikos WM VM DL stress2SD,... negMaikos WM VM DL stress3SD] = sdpeak('negMaikos WM VM DL'); [negMaikos WM VL DL strain1SD,negMaikos WM VL DL strain2SD,negMaikos WM VL DL strain3SD,... negMaikos WM VL DL stress1SD,negMaikos WM VL DL stress2SD,... negMaikos WM VL DL stress3SD] = sdpeak('negMaikos WM VL DL'); [negMaikos GM DL strain1SD,negMaikos GM DL strain2SD,negMaikos GM DL strain3SD,... negMaikos GM DL stress1SD,negMaikos GM DL stress2SD,... negMaikos GM DL stress3SD] = sdpeak('negMaikos GM DL'); 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 % Distraction with Maikos properties [Maikos WM DC DT strain1SD,Maikos WM DC DT strain2SD,Maikos WM DC DT strain3SD,... Maikos WM DC DT stress1SD,Maikos WM DC DT stress2SD,... Maikos WM DC DT stress3SD] = sdpeak('Maikos WM DC DT'); [Maikos WM LC DT strain1SD,Maikos WM LC DT strain2SD,Maikos WM LC DT strain3SD,... Maikos WM LC DT stress1SD,Maikos WM LC DT stress2SD,... Maikos WM LC DT stress3SD] = sdpeak('Maikos WM LC DT'); [Maikos WM VM DT strain1SD,Maikos WM VM DT strain2SD,Maikos WM VM DT strain3SD,... Maikos WM VM DT stress1SD,Maikos WM VM DT stress2SD,... Maikos WM VM DT stress3SD] = sdpeak('Maikos WM VM DT'); [Maikos WM VL DT strain1SD,Maikos WM VL DT strain2SD,Maikos WM VL DT strain3SD,... Maikos WM VL DT stress1SD,Maikos WM VL DT stress2SD,... Maikos WM VL DT stress3SD] = sdpeak('Maikos WM VL DT'); [Maikos GM DT strain1SD,Maikos GM DT strain2SD,Maikos GM DT strain3SD,... Maikos GM DT stress1SD,Maikos GM DT stress2SD,... Maikos GM DT stress3SD] = sdpeak('Maikos GM DT'); 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 % Distraction with negative Maikos properties [negMaikos WM DC DT strain1SD,negMaikos WM DC DT strain2SD,negMaikos WM DC DT strain3SD,... negMaikos WM DC DT stress1SD,negMaikos WM DC DT stress2SD,... negMaikos WM DC DT stress3SD] = sdpeak('negMaikos WM DC DT'); [negMaikos WM LC DT strain1SD,negMaikos WM LC DT strain2SD,negMaikos WM LC DT strain3SD,... negMaikos WM LC DT stress1SD,negMaikos WM LC DT stress2SD,... negMaikos WM LC DT stress3SD] = sdpeak('negMaikos WM LC DT'); [negMaikos WM VM DT strain1SD,negMaikos WM VM DT strain2SD,negMaikos WM VM DT strain3SD,... negMaikos WM VM DT stress1SD,negMaikos WM VM DT stress2SD,... negMaikos WM VM DT stress3SD] = sdpeak('negMaikos WM VM DT'); [negMaikos WM VL DT strain1SD,negMaikos WM VL DT strain2SD,negMaikos WM VL DT strain3SD,... negMaikos WM VL DT stress1SD,negMaikos WM VL DT stress2SD,... negMaikos WM VL DT stress3SD] = sdpeak('negMaikos WM VL DT'); [negMaikos GM DT strain1SD,negMaikos GM DT strain2SD,negMaikos GM DT strain3SD,... negMaikos GM DT stress1SD,negMaikos GM DT stress2SD,... negMaikos GM DT stress3SD] = sdpeak('negMaikos GM DT'); 218 219 220 221 % Save as .mat file save('mechanismFEdataSD.mat') clear; 93 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain Helper function to calculate mean peak values - avgpeak.m 1 2 3 4 5 6 7 8 function [avgpeak strain1,avgpeak strain2,avgpeak strain3,avgpeak stress1,... avgpeak stress2,avgpeak stress3] = avgpeak(region) %% avgpeak −− Calculates and returns average peak strains and stresses from %% zoneslice region % [avgpeak strain1,avgpeak strain2,avgpeak strain3,avgpeak stress1,... % avgpeak stress2,avgpeak stress3] = avgpeak(region) % Can call with only one output, in that case only avgpeak strain1 will be % returned. 9 10 load([region '.mat']); 11 12 13 %% Extract max princ. strain, stress from data and calculate average peak %% stress/strain 14 15 16 17 % Figure out # of elements in region by reading in the column % labels and counting how many start with 'Stress First' numels = length(cell2mat(strfind(eval([region ' labels 0']),'Stress First'))); 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 for i=−5:5 if i>0 data name = [region ' data pos' num2str(abs(i))]; else if i<0 data name = [region ' data neg' num2str(abs(i))]; else data name = [region ' data ' num2str(i)]; end end peak = eval(['sign(max(' data name ')).*max(abs(' data name '));']); avgpeak strain1(i+6) = mean(peak(2+3*numels:1+4*numels)); if nargout <= 1 continue end avgpeak strain2(i+6) = mean(peak(2+4*numels:1+5*numels)); avgpeak strain3(i+6) = mean(peak(2+5*numels:1+6*numels)); 34 avgpeak stress1(i+6) = mean(peak(2:1+numels)); avgpeak stress2(i+6) = mean(peak(2+1*numels:1+2*numels)); avgpeak stress3(i+6) = mean(peak(2+2*numels:1+3*numels)); 35 36 37 38 39 end 40 41 end Helper function to calculate SD of peak values - sdpeak.m 1 2 3 4 5 6 7 8 function [sdpeak strain1,sdpeak strain2,sdpeak strain3,sdpeak stress1,... sdpeak stress2,sdpeak stress3] = sdpeak(region) %% avgpeak −− Calculates and returns average and S.D. of peak strains and stresses from %% zoneslice region % [avgpeak strain1,avgpeak strain2,avgpeak strain3,avgpeak stress1,... % avgpeak stress2,avgpeak stress3] = avgpeak(region,folder) % Can call with only one output, in that case only avgpeak strain1 will be % returned. 9 10 load([region '.mat']); 11 12 13 %% Extract max princ. strain, stress from data and calculate S.D. of peak %% stresses/strains 94 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 14 15 16 17 % Figure out # of elements in region by reading in the column % labels and counting how many start with 'Stress First' numels = length(cell2mat(strfind(eval([region ' labels 0']),'Stress First'))); 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 for i=−5:5 if i>0 data name = [region ' data pos' num2str(abs(i))]; else if i<0 data name = [region ' data neg' num2str(abs(i))]; else data name = [region ' data ' num2str(i)]; end end peak = eval(['sign(max(' data name ')).*max(abs(' data name '));']); sdpeak strain1(i+6) = std(peak(2+3*numels:1+4*numels)); if nargout <= 1 continue end sdpeak strain2(i+6) = std(peak(2+4*numels:1+5*numels)); sdpeak strain3(i+6) = std(peak(2+5*numels:1+6*numels)); 34 sdpeak stress1(i+6) = std(peak(2:1+numels)); sdpeak stress2(i+6) = std(peak(2+1*numels:1+2*numels)); sdpeak stress3(i+6) = std(peak(2+2*numels:1+3*numels)); 35 36 37 38 end 39 40 end C.6.3 1 2 Plotting of regional FE stress and strain as function of slice position zonesliceplot.m %% Load stress and strain simulation results from .mat and plot as function %% of slice position 3 4 clc; clear; close all; 5 6 7 8 % Load FE data load mechanismFEdata.mat; load mechanismFEdataSD.mat; 9 10 11 12 13 14 %% Specify default plot settings set(0,'DefaultAxesLineStyle','−') set(0,'DefaultLineLineWidth',2.5) set(0,'DefaultAxesFontSize',14) set(0,'DefaultAxesColorOrder',[0 0 1;1 0 0;0 0.7 0; 0.5 0.5 0]) 15 16 17 savefigs = 1; % 1 = save figures closefigs = 1; % 1 = close figures after saving 18 19 20 21 %% Plot max princ. stresses and strains as function of slice position distance = −5:5; param = {'strain1','strain2','strain3','stress1','stress2','stress3'}; 22 23 24 25 26 27 28 29 30 % Maikos properties for i=1:6 sliceplot('Maikos sliceplot('Maikos sliceplot('Maikos sliceplot('Maikos sliceplot('Maikos end WM DC',param{i},savefigs,closefigs) WM LC',param{i},savefigs,closefigs) WM VM',param{i},savefigs,closefigs) WM VL',param{i},savefigs,closefigs) GM',param{i},savefigs,closefigs) 95 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 31 32 33 34 35 36 37 38 39 % % % % % % % % Negative Maikos properties for i=1:1 sliceplot('negMaikos WM DC',param{i},savefigs,closefigs) sliceplot('negMaikos WM LC',param{i},savefigs,closefigs) sliceplot('negMaikos WM VM',param{i},savefigs,closefigs) sliceplot('negMaikos WM VL',param{i},savefigs,closefigs) % sliceplot('negMaikos GM',param{i},savefigs,closefigs) end Helper function to format slice distance plots - sliceplot.m 1 2 3 function [] = sliceplot(region,param,savefigs,closefigs) %% sliceplot −− Plots FE data as function of slice position % [] = sliceplot(region,param,savefigs,closefigs) 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 %% Format param as TeX symbol switch param case {'strain1'} paramsym = '{\epsilon} case {'strain2'} paramsym = '{\epsilon} case {'strain3'} paramsym = '{\epsilon} case {'stress1'} paramsym = '{\sigma} 1 case {'stress2'} paramsym = '{\sigma} 2 case {'stress3'} paramsym = '{\sigma} 3 end 1'; 2'; 3'; (MPa)'; (MPa)'; (MPa)'; 20 21 22 23 24 25 26 27 %% Assign data corresponding to region and param to variable 'paramdata' paramdataCT = evalin('base',[region ' CT ' param]); paramdataDL = evalin('base',[region ' DL ' param]); paramdataDT = evalin('base',[region ' DT ' param]); paramSDCT = evalin('base',[region ' CT ' param 'SD']); paramSDDL = evalin('base',[region ' DL ' param 'SD']); paramSDDT = evalin('base',[region ' DT ' param 'SD']); 28 29 30 31 32 distance = −5:5; paramdata = [paramdataCT; paramdataDL; paramdataDT]'; distance = [distance; distance; distance]'; paramSD = [paramSDCT; paramSDDL; paramSDDT]'; 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 %% Plot figure figure hold on title([strrep(region,' ','\ ') ': ' paramsym]) ylabel(paramsym) xlabel('<<Caudal Distance from epicenter (mm) switch param case {'strain1'} ylim([0 0.6]) case {'strain2'} ylim([0 0.3]) case {'strain3'} ylim([−0.4 0]) case {'stress1'} ylim([0 0.11]) case {'stress2'} Rostral>>') 96 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 50 51 52 53 54 55 ylim([0 0.15]) case {'stress3'} ylim([0 0.9]) end errorbar(distance,paramdata,paramSD) hold off 56 57 58 59 60 61 62 if savefigs saveas(gcf,[cd '\Figures\Distance\' region ' ' param ' dist.eps'],'psc2') end if closefigs close(gcf) end C.6.4 1 Pooled and regional correlation plots of FE stress and strain with tissue damage (dislocation epicentre points omitted) zoneslicecorrplotOMIT.m %% Load stress and strain simulation results from .mat and display in correlation plots 2 3 clc; clear; close all; 4 5 6 7 8 % Load Choo data load choodextrancounts.mat; % Load FE data load mechanismFEdata.mat; 9 10 11 12 13 14 %% Specify default plot settings set(0,'DefaultAxesLineStyle','−.o') set(0,'DefaultLineLineWidth',2) set(0,'DefaultAxesFontSize',20) set(0,'DefaultAxesColorOrder',[0 0 1;1 0 0;0 0.7 0; 0.5 0.5 0]) 15 16 17 18 plotlinreg = 1; % 1 = add linear regression lines to correlation plots savefigs = 1; % 1 = save figures closefigs = 1; % 1 = close figures after saving 19 20 21 22 23 24 25 %% Concatenated data param = {'strain1','strain2','strain3','stress1','stress2','stress3'}; % Maikos properties for i=1:6 corrplotconcatOMIT('Maikos WM',param{i},plotlinreg,savefigs,closefigs) end 26 27 28 29 30 % % Negative Maikos properties % for i = 1:6 % corrplotconcatOMIT('negMaikos WM',param{i},plotlinreg,savefigs,closefigs) % end 31 32 33 34 35 36 37 38 39 40 41 42 %% Regional plots % Maikos properties for i = 1:6 corrplotOMIT('Maikos WM DC',param{i},plotlinreg,savefigs,closefigs) corrplotOMIT('Maikos WM LC',param{i},plotlinreg,savefigs,closefigs) corrplotOMIT('Maikos WM VM',param{i},plotlinreg,savefigs,closefigs) corrplotOMIT('Maikos WM VL',param{i},plotlinreg,savefigs,closefigs) corrplotOMIT('Maikos GM',param{i},plotlinreg,savefigs,closefigs) end % % % Negative Maikos properties 97 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 43 44 45 46 47 48 49 % for i = 1:6 % corrplotOMIT('negMaikos WM DC',param{i},plotlinreg,savefigs,closefigs) % corrplotOMIT('negMaikos WM LC',param{i},plotlinreg,savefigs,closefigs) % corrplotOMIT('negMaikos WM VM',param{i},plotlinreg,savefigs,closefigs) % corrplotOMIT('negMaikos WM VL',param{i},plotlinreg,savefigs,closefigs) % corrplotOMIT('negMaikos GM',param{i},plotlinreg,savefigs,closefigs) % end Helper function to format pooled WM plots - corrplotconcatOMIT.m 1 2 3 4 5 6 7 8 9 10 11 12 function [] = corrplotconcatOMIT(region,param,plotlinreg,savefigs,closefigs) %% corrplotconcatOMIT −− Calculates and plots lumped data correlation %% scatter plots of FE result parameter vs. cellular permeability counts. %% Omits data points at injury epicenter (for dislocation) because %% they are underestimates of cellular damage due to necrosis. % [] = corrplotconcat(region,param,plotlinreg,savefigs,closefigs) % region = 'Maikos WM' e.g. % param = 'strain1','stress1', etc. (1−3) % plotlinreg = 0,1 (0 will prevent linear correlation results being shown % on plot) % savefigs = 1 −> save figures % closefigs = 1 −> close figures after saving 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 %% Format param as TeX symbol switch param case {'strain1'} paramsym = '{\epsilon} 1'; case {'strain2'} paramsym = '{\epsilon} 2'; case {'strain3'} paramsym = '{\epsilon} 3'; case {'stress1'} paramsym = '{\sigma} 1'; case {'stress2'} paramsym = '{\sigma} 2'; case {'stress3'} paramsym = '{\sigma} 3'; end 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 %% Plot figure figure hold all title(['Scatter plot lumped data ' strrep(region,' ','\ ') ': '... paramsym ' vs. cell count']) ylabel('Axons/mm') if strfind(param,'stress') xlabel([paramsym ' (MPa)']) ymin = 0; ymax = 70; xmin = −0.1; xmax = 0.2; else xlabel(paramsym) ymin = 0; ymax = 70; xmin = −0.4; xmax = 0.6; end ylim([ymin ymax]) xlim([xmin xmax]) 47 48 49 50 51 % Contusion mech = ' CT'; %% Calculate and concatenate mean of Choo data for all WM zones mean choo = [evalin('base','mean(choo WM DC CT)')... 98 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 52 53 evalin('base','mean(choo WM LC CT)')... evalin('base','mean(choo WM VM CT)') evalin('base','mean(choo WM VL CT)')]; 54 55 56 57 58 59 60 %% Assign concatenated data corresponding to all WM zones and one param %% to variable 'paramdata' paramdata = [evalin('base',[region(1:end) ' DC' mech ' ' param])... evalin('base',[region(1:end) ' LC' mech ' ' param])... evalin('base',[region(1:end) ' VM' mech ' ' param])... evalin('base',[region(1:end) ' VL' mech ' ' param])]; 61 62 63 64 65 %% Sort concatenated data pairs from least strain to highest strain %% (prevents multiple lines during regression plotting) [paramdata,index] = sort(paramdata); mean choo = mean choo(index); 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 %% Optionally perform linear regression if(plotlinreg) [R,p] = corrcoef(paramdata,mean choo); [P,ErrorEst] = polyfit(paramdata,mean choo,1); [bestfit,delta] = polyval(P,paramdata,ErrorEst); end scatter(paramdata,mean choo,'MarkerFaceColor','b') if(plotlinreg) text(xmin+0.05*(xmax−xmin),ymax−0.05*(ymax−ymin),... ['Rˆ2=' num2str(R(2)ˆ2) ' p=' num2str(p(2))],'Color','b') plot(paramdata,bestfit,'−b') plot(paramdata,bestfit+2*delta,'−−b',... paramdata,bestfit−2*delta,'−−b') end 81 82 83 84 85 86 87 % Dislocation mech = ' DL'; %% Calculate and concatenate mean of Choo data for all WM zones mean choo = [evalin('base','mean(choo WM DC DL)')... evalin('base','mean(choo WM LC DL)')... evalin('base','mean(choo WM VM DL)') evalin('base','mean(choo WM VL DL)')]; 88 89 90 91 92 93 94 %% Assign concatenated data corresponding to all WM zones and one param %% to variable 'paramdata' paramdata = [evalin('base',[region(1:end) ' DC' mech ' ' param])... evalin('base',[region(1:end) ' LC' mech ' ' param])... evalin('base',[region(1:end) ' VM' mech ' ' param])... evalin('base',[region(1:end) ' VL' mech ' ' param])]; 95 96 97 98 99 100 101 102 103 104 105 %% Omit injury epicenter datapoints corresponding to necrotic zones %% (dislocation) mean choo for i=1:4 necrotic stress strain(i) = paramdata(1−i+6+(i−1)*(11)); mean choo(1−i+6+(i−1)*(11)) = []; paramdata(1−i+6+(i−1)*(11)) = []; end % mean choo % necrotic stress strain 106 107 108 109 110 %% Sort concatenated data pairs from least strain to highest strain %% (prevents multiple lines during regression plotting) [paramdata,index] = sort(paramdata); mean choo = mean choo(index); 111 112 113 %% Optionally perform linear regression if(plotlinreg) 99 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain [R,p] = corrcoef(paramdata,mean choo); [P,ErrorEst] = polyfit(paramdata,mean choo,1); [bestfit,delta] = polyval(P,paramdata,ErrorEst); 114 115 116 117 118 119 120 121 122 123 124 125 end scatter(paramdata,mean choo,'MarkerFaceColor','r') if(plotlinreg) text(xmin+0.05*(xmax−xmin),ymax−0.15*(ymax−ymin),... ['Rˆ2=' num2str(R(2)ˆ2) ' p=' num2str(p(2))],'Color','r') plot(paramdata,bestfit,'−r') plot(paramdata,bestfit+2*delta,'−−r',... paramdata,bestfit−2*delta,'−−r') end 126 127 128 129 130 131 132 % Distraction mech = ' DT'; %% Calculate and concatenate mean of Choo data for all WM zones mean choo = [evalin('base','mean(choo WM DC DT)')... evalin('base','mean(choo WM LC DT)')... evalin('base','mean(choo WM VM DT)') evalin('base','mean(choo WM VL DT)')]; 133 134 135 136 137 138 139 %% Assign concatenated data corresponding to all WM zones and one param %% to variable 'paramdata' paramdata = [evalin('base',[region(1:end) ' DC' mech ' ' param])... evalin('base',[region(1:end) ' LC' mech ' ' param])... evalin('base',[region(1:end) ' VM' mech ' ' param])... evalin('base',[region(1:end) ' VL' mech ' ' param])]; 140 141 142 143 144 %% Sort concatenated data pairs from least strain to highest strain %% (prevents multiple lines during regression plotting) [paramdata,index] = sort(paramdata); mean choo = mean choo(index); 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 %% Optionally perform linear regression if(plotlinreg) [R,p] = corrcoef(paramdata,mean choo); [P,ErrorEst] = polyfit(paramdata,mean choo,1); [bestfit,delta] = polyval(P,paramdata,ErrorEst); end scatter(paramdata,mean choo,'MarkerFaceColor',[0 0.7 0]) if(plotlinreg) text(xmin+0.05*(xmax−xmin),ymax−0.25*(ymax−ymin),... ['Rˆ2=' num2str(R(2)ˆ2) ' p=' num2str(p(2))],'Color',[0 0.7 0]) plot(paramdata,bestfit,'−','Color',[0 0.7 0]) plot(paramdata,bestfit+2*delta,'−−',... paramdata,bestfit−2*delta,'−−','Color',[0 0.7 0]) end 160 161 hold off 162 163 164 165 166 167 168 169 if savefigs saveas(gcf,[cd '\Figures\CorrelationsOMIT\' region ' ' param... ' corr.eps'],'psc2') end if closefigs close(gcf) end 170 171 end Helper function to format regional plots - corrplotOMIT.m 100 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 1 2 3 4 5 6 7 8 9 10 function [] = corrplotOMIT(region,param,plotlinreg,savefigs,closefigs) %% corrplotOMIT −− Calculates and plots correlation scatter plots of FE result %% parameter vs. cellular permeability counts. Omits data points at injury epicenter %% (for dislocation) because they are underestimates of cellular damage due to necrosis. % [] = corrplot(region,param,plotlinreg,savefigs,closefigs) % param = 'strain1','stress1', etc. (1−3) % plotlinreg = 0,1 (0 will prevent linear correlation results being shown % on plot) % savefigs = 1 −> save figures % closefigs = 1 −> close figures after saving 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 %% Format param as TeX symbol switch param case {'strain1'} paramsym = '{\epsilon} 1'; case {'strain2'} paramsym = '{\epsilon} 2'; case {'strain3'} paramsym = '{\epsilon} 3'; case {'stress1'} paramsym = '{\sigma} 1'; case {'stress2'} paramsym = '{\sigma} 2'; case {'stress3'} paramsym = '{\sigma} 3'; end 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 %% Plot figure figure hold all title(['Scatter plot ' strrep(region,' ','\ ') ': ' paramsym ' vs. cell count']) if strfind(region,' GM') % special case for changing yaxis label when plotting GM data ylabel('% of cells') else ylabel('Axons/mm') end if strfind(param,'stress') xlabel([paramsym ' (MPa)']) ymin = 0; ymax = 70; xmin = −0.1; xmax = 0.2; else xlabel(paramsym) ymin = 0; ymax = 70; xmin = −0.4; xmax = 0.6; end ylim([ymin ymax]) xlim([xmin xmax]) 48 49 50 51 52 53 54 % Contusion mech = ' CT'; %% Calculate mean of Choo data for current region % Find character index where 'Maikos' appears in region name MaikosI = strfind(region,'Maikos'); mean choo = evalin('base',['mean(choo ' region(MaikosI+7:end) mech ')']); 55 56 57 %% Assign data corresponding to region and param to variable 'paramdata' paramdata = evalin('base',[region mech ' ' param]); 58 59 60 61 62 %% Sort data pairs from least strain to highest strain %% (prevents multiple lines during regression plotting) [paramdata,index] = sort(paramdata); mean choo = mean choo(index); 101 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 %% Optionally perform linear regression if(plotlinreg) [R,p] = corrcoef(paramdata,mean choo); [P,ErrorEst] = polyfit(paramdata,mean choo,1); [bestfit,delta] = polyval(P,paramdata,ErrorEst); end %% Make scatter plot scatter(paramdata,mean choo,'MarkerFaceColor','b') if(plotlinreg) text(xmin+0.05*(xmax−xmin),ymax−0.05*(ymax−ymin),... ['Rˆ2=' num2str(R(2)ˆ2) ' p=' num2str(p(2))],'Color','b') plot(paramdata,bestfit,'−b') plot(paramdata,bestfit+2*delta,'−−b',... paramdata,bestfit−2*delta,'−−b') end 79 80 81 82 83 84 85 % Dislocation mech = ' DL'; %% Calculate mean of Choo data for current region % Find character index where 'Maikos' appears in region name MaikosI = strfind(region,'Maikos'); mean choo = evalin('base',['mean(choo ' region(MaikosI+7:end) mech ')']); 86 87 88 %% Assign data corresponding to region and param to variable 'paramdata' paramdata = evalin('base',[region mech ' ' param]); 89 90 91 92 93 94 95 96 97 98 99 %% Omit injury epicenter datapoints corresponding to necrotic zones %% (dislocation) mean choo for i=1:1 necrotic stress strain(i) = paramdata(1−i+6+(i−1)*(11)); mean choo(1−i+6+(i−1)*(11)) = []; paramdata(1−i+6+(i−1)*(11)) = []; end % mean choo % necrotic stress strain 100 101 102 103 104 %% Sort data pairs from least strain to highest strain %% (prevents multiple lines during regression plotting) [paramdata,index] = sort(paramdata); mean choo = mean choo(index); 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 %% Optionally perform linear regression if(plotlinreg) [R,p] = corrcoef(paramdata,mean choo); [P,ErrorEst] = polyfit(paramdata,mean choo,1); [bestfit,delta] = polyval(P,paramdata,ErrorEst); end %% Make scatter plot scatter(paramdata,mean choo,'MarkerFaceColor','r') if(plotlinreg) text(xmin+0.05*(xmax−xmin),ymax−0.15*(ymax−ymin),... ['Rˆ2=' num2str(R(2)ˆ2) ' p=' num2str(p(2))],'Color','r') plot(paramdata,bestfit,'−r') plot(paramdata,bestfit+2*delta,'−−r',... paramdata,bestfit−2*delta,'−−r') end 121 122 123 124 % Distraction mech = ' DT'; %% Calculate mean of Choo data for current region 102 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 125 126 127 % Find character index where 'Maikos' appears in region name MaikosI = strfind(region,'Maikos'); mean choo = evalin('base',['mean(choo ' region(MaikosI+7:end) mech ')']); 128 129 130 %% Assign data corresponding to region and param to variable 'paramdata' paramdata = evalin('base',[region mech ' ' param]); 131 132 133 134 135 %% Sort data pairs from least strain to highest strain %% (prevents multiple lines during regression plotting) [paramdata,index] = sort(paramdata); mean choo = mean choo(index); 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 %% Optionally perform linear regression if(plotlinreg) [R,p] = corrcoef(paramdata,mean choo); [P,ErrorEst] = polyfit(paramdata,mean choo,1); [bestfit,delta] = polyval(P,paramdata,ErrorEst); end %% Make scatter plot scatter(paramdata,mean choo,'MarkerFaceColor',[0 0.7 0]) if(plotlinreg) text(xmin+0.05*(xmax−xmin),ymax−0.25*(ymax−ymin),... ['Rˆ2=' num2str(R(2)ˆ2) ' p=' num2str(p(2))],'Color',[0 0.7 0]) plot(paramdata,bestfit,'−','Color',[0 0.7 0]) plot(paramdata,bestfit+2*delta,'−−',... paramdata,bestfit−2*delta,'−−','Color',[0 0.7 0]) end 152 153 hold off 154 155 156 157 158 159 160 161 if savefigs saveas(gcf,[cd '\Figures\CorrelationsOMIT\' region ' ' param... ' corr.eps'],'psc2') end if closefigs close(gcf) end 162 163 end C.6.5 1 2 Regional time history plots of FE stress and strain - zonecurveplot.m %% Load stress and strain simulation curves from .mat and plot as function %% of time 3 4 clc; clear; close all; 5 6 7 8 9 %% Specify default plot settings set(0,'DefaultAxesLineStyle','−') set(0,'DefaultLineLineWidth',0.5) set(0,'DefaultAxesFontSize',20) 10 11 12 savefigs = 1; % 1 = save figures closefigs = 1; % 1 = close figures after saving 13 14 15 16 17 18 %% Plot max princ. stresses and strains as function of time mech = {'CT','DL','DT'}; prop = {'Maikos','negMaikos'}; zone = {'WM DC','WM LC','WM VM','WM VL','GM'}; param = {'strain1','strain2','strain3','stress1','stress2','stress3'}; 19 103 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain 20 21 22 23 24 25 26 27 28 29 30 for i=3:3 % Cycle through mechanisms for j=1:1 % Cycle through Maikos,negMaikos for k=1:5 % Cycle through zones for l=1:6 % Cycle through params region = [prop{j} ' ' zone{k} ' ' mech{i}]; curveplot(region,param{l},savefigs,closefigs) end end end end Helper function to format time history plots - curveplot.m 1 2 3 4 function [] = curveplot(region,param,savefigs,closefigs) %% curveplot −− Plots zone slice time histories for chosen parameter % [] = curveplot(region,param,savefigs,closefigs) % Plots raw parameter curves as function of time (ms). 5 6 load([region '.mat']); 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 %% Format param as TeX symbol switch param case {'strain1'} paramsym = '{\epsilon} 1'; multip = 3; % numels multiplier for indexing into data case {'strain2'} paramsym = '{\epsilon} 2'; multip = 4; case {'strain3'} paramsym = '{\epsilon} 3'; multip = 5; case {'stress1'} paramsym = '{\sigma} 1 (MPa)'; multip = 1; case {'stress2'} paramsym = '{\sigma} 2 (MPa)'; multip = 2; case {'stress3'} paramsym = '{\sigma} 3 (MPa)'; multip = 3; end 29 30 31 32 33 %% Plot zone slice curves % Figure out # of elements in region by reading in the column % labels and counting how many start with 'Stress First' numels = length(cell2mat(strfind(eval([region ' labels 0']),'Stress First'))); 34 35 36 37 38 39 40 41 42 43 for i=−5:5 if i>0 data name = [region ' data pos' num2str(abs(i))]; else if i<0 data name = [region ' data neg' num2str(abs(i))]; else data name = [region ' data ' num2str(i)]; end end 44 45 46 47 figure hold on title([strrep(region,' ','\ ') ' ' num2str(i) 'mm']) 104 C.6. FE data extraction and linear regression of tissue damage versus maximum principal strain ylabel(paramsym) xlabel('Time (ms)') switch param case {'strain1'} ylim([0 0.6]) case {'strain2'} ylim([0 0.3]) case {'strain3'} ylim([−0.4 0]) case {'stress1'} ylim([0 0.11]) case {'stress2'} ylim([0 0.15]) case {'stress3'} ylim([0 0.9]) end eval(['plot(' data name '(:,1),' data name... '(:,2+multip*numels:1+(1+multip)*numels))']); hold off 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 if savefigs saveas(gcf,[cd '\Figures\Time\' region ' ' param... ' time ' num2str(i) 'mm.eps'],'psc2') end if closefigs close(gcf) end 68 69 70 71 72 73 74 75 end 76 77 end 105 Appendix D Histological and FE parameters as function of distance from injury epicentre Section D.1 shows the distributions of cell membrane permeability as a function of distance from injury epicentre for each of the contusion (blue), dislocation (red), and distraction (green) injury mechanisms. In section D.2 the same colour scheme and regional layout are used to plot each of the three principal strains and three principal stresses as a function of distance from epicentre. Peak values for each element were calculated and then averaged within regions to yield the plotted mean peak values for each region. Each parameter is presented on a separate page. 106 D.1. Histological data D.1 Histological data Choo Dextran counts − WM DC Choo Dextran counts − WM LC 70 70 Contusion Dislocation Distraction 60 60 50 Axons/mm Axons/mm 50 40 30 40 30 20 20 10 10 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 0 −6 6 −4 <<Caudal (a) Dorsal white 6 4 Rostral>> 6 Choo Dextran counts − WM VL 70 60 60 50 50 Axons/mm Axons/mm Choo Dextran counts − WM VM 40 30 40 30 20 20 10 10 −4 <<Caudal 4 Rostral>> (b) Lateral white 70 0 −6 −2 0 2 Distance from epicenter (mm) −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 0 −6 −4 <<Caudal (c) Ventro-medial white −2 0 2 Distance from epicenter (mm) (d) Ventro-lateral white Choo Dextran counts − GM 100 90 80 % of cells 70 60 50 40 30 20 10 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 (e) Gray matter Figure D.1: Distance plots – Choo cellular permeability. Error bars denote the standard error of the mean. 107 D.2. FE results D.2 FE results Maikos_WM_LC: ε1 0.5 0.5 0.4 0.4 ε1 ε1 Maikos_WM_DC: ε1 0.3 0.3 0.2 0.2 0.1 0.1 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 0 −6 6 −4 <<Caudal (a) Dorsal white −2 0 2 Distance from epicenter (mm) 0.5 0.4 0.4 ε1 ε1 4 Rostral>> 6 Maikos_WM_VL: ε1 0.5 0.3 0.3 0.2 0.2 0.1 0.1 −4 <<Caudal 6 (b) Lateral white Maikos_WM_VM: ε1 0 −6 4 Rostral>> −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 0 −6 −4 <<Caudal (c) Ventro-medial white −2 0 2 Distance from epicenter (mm) (d) Ventro-lateral white Maikos_GM: ε1 0.5 ε 1 0.4 0.3 0.2 0.1 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 (e) Gray matter Figure D.2: Distance plots – first principal strain (Maikos properties). Mean peak values plotted for each region, with error bars to denote the standard deviation within region. 108 D.2. FE results Maikos_WM_LC: ε2 0.25 0.25 0.2 0.2 ε2 ε2 Maikos_WM_DC: ε2 0.15 0.15 0.1 0.1 0.05 0.05 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 0 −6 6 −4 <<Caudal (a) Dorsal white −2 0 2 Distance from epicenter (mm) 0.25 0.2 0.2 ε2 ε2 4 Rostral>> 6 Maikos_WM_VL: ε2 0.25 0.15 0.15 0.1 0.1 0.05 0.05 −4 <<Caudal 6 (b) Lateral white Maikos_WM_VM: ε2 0 −6 4 Rostral>> −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 0 −6 −4 <<Caudal (c) Ventro-medial white −2 0 2 Distance from epicenter (mm) (d) Ventro-lateral white Maikos_GM: ε2 0.25 ε2 0.2 0.15 0.1 0.05 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 (e) Gray matter Figure D.3: Distance plots – second principal strain (Maikos properties). Mean peak values plotted for each region, with error bars to denote the standard deviation within region. 109 D.2. FE results Maikos_WM_LC: ε3 0 −0.05 −0.05 −0.1 −0.1 −0.15 −0.15 ε3 ε3 Maikos_WM_DC: ε3 0 −0.2 −0.2 −0.25 −0.25 −0.3 −0.3 −0.35 −0.35 −0.4 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> −0.4 −6 6 (a) Dorsal white −4 <<Caudal −0.05 −0.05 −0.1 −0.1 −0.15 −0.15 ε3 ε3 0 −0.2 4 Rostral>> 6 −0.2 −0.25 −0.25 −0.3 −0.3 −0.35 −0.35 −2 0 2 Distance from epicenter (mm) 6 Maikos_WM_VL: ε3 0 −4 <<Caudal 4 Rostral>> (b) Lateral white Maikos_WM_VM: ε3 −0.4 −6 −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 −0.4 −6 (c) Ventro-medial white −4 <<Caudal −2 0 2 Distance from epicenter (mm) (d) Ventro-lateral white Maikos_GM: ε3 0 −0.05 −0.1 ε3 −0.15 −0.2 −0.25 −0.3 −0.35 −0.4 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 (e) Gray matter Figure D.4: Distance plots – third principal strain (Maikos properties). Mean peak values plotted for each region, with error bars to denote the standard deviation within region. 110 D.2. FE results Maikos_WM_LC: σ1 (MPa) 0.1 0.1 0.09 0.09 0.08 0.08 0.07 0.07 σ1 (MPa) σ1 (MPa) Maikos_WM_DC: σ1 (MPa) 0.06 0.05 0.06 0.05 0.04 0.04 0.03 0.03 0.02 0.02 0.01 0.01 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 0 −6 6 −4 <<Caudal (a) Dorsal white −2 0 2 Distance from epicenter (mm) 0.1 0.1 0.09 0.08 0.08 0.07 0.07 σ1 (MPa) σ1 (MPa) 4 Rostral>> 6 Maikos_WM_VL: σ1 (MPa) 0.09 0.06 0.05 0.06 0.05 0.04 0.04 0.03 0.03 0.02 0.02 0.01 0.01 −4 <<Caudal 6 (b) Lateral white Maikos_WM_VM: σ1 (MPa) 0 −6 4 Rostral>> −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 0 −6 −4 <<Caudal (c) Ventro-medial white −2 0 2 Distance from epicenter (mm) (d) Ventro-lateral white Maikos_GM: σ1 (MPa) 0.1 0.09 0.08 σ1 (MPa) 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 (e) Gray matter Figure D.5: Distance plots – first principal stress (Maikos properties). Mean peak values plotted for each region, with error bars to denote the standard deviation within region. 111 D.2. FE results Maikos_WM_DC: σ2 (MPa) Maikos_WM_LC: σ2 (MPa) σ2 (MPa) 0.1 σ2 (MPa) 0.1 0.05 0 −6 0.05 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 0 −6 6 −4 <<Caudal (a) Dorsal white −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 4 Rostral>> 6 (b) Lateral white Maikos_WM_VM: σ2 (MPa) Maikos_WM_VL: σ2 (MPa) σ2 (MPa) 0.1 σ2 (MPa) 0.1 0.05 0 −6 0.05 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 0 −6 −4 <<Caudal (c) Ventro-medial white −2 0 2 Distance from epicenter (mm) (d) Ventro-lateral white Maikos_GM: σ2 (MPa) σ2 (MPa) 0.1 0.05 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 (e) Gray matter Figure D.6: Distance plots – second principal stress (Maikos properties). Mean peak values plotted for each region, with error bars to denote the standard deviation within region. 112 D.2. FE results Maikos_WM_LC: σ3 (MPa) 0.9 0.8 0.8 0.7 0.7 0.6 0.6 σ3 (MPa) σ3 (MPa) Maikos_WM_DC: σ3 (MPa) 0.9 0.5 0.4 0.5 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 0 −6 6 −4 <<Caudal (a) Dorsal white −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 Maikos_WM_VL: σ3 (MPa) 0.9 0.8 0.8 0.7 0.7 0.6 0.6 σ3 (MPa) σ3 (MPa) Maikos_WM_VM: σ3 (MPa) 0.5 0.4 0.5 0.4 0.3 0.3 0.2 0.2 0.1 0.1 −4 <<Caudal 6 (b) Lateral white 0.9 0 −6 4 Rostral>> −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 0 −6 −4 <<Caudal (c) Ventro-medial white −2 0 2 Distance from epicenter (mm) (d) Ventro-lateral white Maikos_GM: σ3 (MPa) 0.9 0.8 0.7 3 σ (MPa) 0.6 0.5 0.4 0.3 0.2 0.1 0 −6 −4 <<Caudal −2 0 2 Distance from epicenter (mm) 4 Rostral>> 6 (e) Gray matter Figure D.7: Distance plots – third principal stress (Maikos properties). Mean peak values plotted for each region, with error bars to denote the standard deviation within region. 113 Appendix E Correlation plots Experimental measurements made by Choo et al. [14] of spinal cord cell membrane permeability are plotted against corresponding stress and strain parameters from FE simulation results for each injury mechanism. Mean data points are matched according to cross-sectional region and axial slice position. Section E.1 shows summary plots wherein data from all four white matter regions have been pooled together for each injury mechanism, while section E.2 shows correlations within each region. As in Appendix D, blue represents contusion, red dislocation, and green distraction. Linear regression lines of best fit are plotted with solid lines in the colour corresponding to each mechanism, along with 95% confidence intervals in dashed lines. R2 values are also listed for each mechanism colour, with corresponding p values. 114 E.1. Correlations for pooled white matter regions E.1 Correlations for pooled white matter regions Scatter plot lumped data Maikos_WM: ε1 vs. cell count 70 Scatter plot lumped data Maikos_WM: σ1 vs. cell count 70 2 R =0.85789 p=2.11e−019 60 60 R2=0.53963 p=6.8125e−008 2 R =0.30593 p=9.8693e−005 50 Axons/mm Axons/mm 50 40 30 10 10 0 ε 0.2 0.4 R2=0.14961 p=0.0094964 30 20 −0.2 R2=0.65913 p=2.0578e−010 40 20 0 −0.4 2 R =0.71589 p=4.7845e−013 0 −0.1 0.6 −0.05 0 1 (a) First principal strain 60 R2=0.56135 p=4.901e−009 60 2 R =0.31503 p=0.00016467 50 Axons/mm Axons/mm R =0.35994 p=1.6756e−005 40 30 10 10 0 ε 0.2 0.4 R2=0.48679 p=5.6267e−007 R2=0.39027 p=5.8289e−006 30 20 −0.2 R2=0.55001 p=8.4542e−009 40 20 0 −0.4 0 −0.1 0.6 −0.05 0 2 (c) Second principal strain 70 2 R =0.81747 p=4.1447e−017 60 2 R =0.58513 p=9.088e−009 R2=0.45487 p=5.1776e−007 50 Axons/mm Axons/mm 50 40 30 10 10 0 ε3 0.2 (e) Third principal strain 0.4 0.6 R2=0.3131 p=7.8545e−005 2 R =0.35462 p=5.0445e−005 R2=0.39374 p=5.1515e−006 30 20 −0.2 0.2 40 20 0 −0.4 0.15 Scatter plot lumped data Maikos_WM: σ3 vs. cell count 3 60 0.05 0.1 σ2 (MPa) (d) Second principal stress Scatter plot lumped data Maikos_WM: ε vs. cell count 70 0.2 Scatter plot lumped data Maikos_WM: σ2 vs. cell count 70 2 50 0.15 (b) First principal stress Scatter plot lumped data Maikos_WM: ε2 vs. cell count 70 0.05 0.1 σ1 (MPa) 0 −0.1 −0.05 0 0.05 0.1 σ3 (MPa) 0.15 0.2 (f ) Third principal stress Figure E.1: Parameter correlations – pooled data for all white matter regions (Maikos properties) 115 E.2. Regional correlations E.2 Regional correlations Scatter plot Maikos_WM_DC: ε1 vs. cell count 70 60 Scatter plot Maikos_WM_DC: σ1 vs. cell count 70 R2=0.94455 p=5.8951e−007 60 R2=0.38336 p=0.056283 2 R =0.12102 p=0.29449 50 Axons/mm Axons/mm 50 40 30 10 10 0 ε1 0.2 0.4 0 −0.1 0.6 (a) First principal strain Axons/mm 50 R2=0.8983 p=9.216e−006 60 R2=0.56372 p=0.01233 R2=0.18208 p=0.19059 50 40 30 10 10 0 ε2 0.2 0.4 0 −0.1 0.6 (c) Second principal strain 70 2 60 2 R =0.83255 p=0.00023111 R2=0.00022818 p=0.96484 50 Axons/mm Axons/mm 50 40 30 10 10 0 ε3 0.2 (e) Third principal strain R2=0.00037223 p=0.9551 −0.05 0 0.05 0.1 σ2 (MPa) 0.15 0.2 0.4 0.6 R2=0.4978 p=0.015279 R2=0.78873 p=0.000598 R2=0.0034989 p=0.86284 30 20 −0.2 R2=0.69755 p=0.0026321 40 20 0 −0.4 R2=0.91692 p=3.6808e−006 Scatter plot Maikos_WM_DC: σ3 vs. cell count R =0.98598 p=1.1905e−009 60 0.2 (d) Second principal stress Scatter plot Maikos_WM_DC: ε3 vs. cell count 70 0.15 30 20 −0.2 0.05 0.1 σ1 (MPa) 40 20 0 −0.4 0 Scatter plot Maikos_WM_DC: σ2 vs. cell count 70 Axons/mm 60 −0.05 (b) First principal stress Scatter plot Maikos_WM_DC: ε2 vs. cell count 70 R2=0.33611 p=0.06156 30 20 −0.2 R2=0.48453 p=0.025362 40 20 0 −0.4 R2=0.64316 p=0.002984 0 −0.1 −0.05 0 0.05 0.1 σ3 (MPa) 0.15 0.2 (f ) Third principal stress Figure E.2: Parameter correlations – white matter dorsal column region (Maikos properties) 116 E.2. Regional correlations Scatter plot Maikos_WM_LC: ε1 vs. cell count 70 60 Scatter plot Maikos_WM_LC: σ1 vs. cell count 70 R2=0.96188 p=1.0842e−007 60 R2=0.95705 p=9.4722e−007 2 R =0.32694 p=0.066091 50 Axons/mm Axons/mm 50 40 30 10 10 0 ε1 0.2 0.4 R2=0.14589 p=0.24637 30 20 −0.2 R2=0.94521 p=2.5191e−006 40 20 0 −0.4 R2=0.96375 p=8.6364e−008 0 −0.1 0.6 −0.05 (a) First principal strain Axons/mm 50 R2=0.94141 p=7.5611e−007 60 R2=0.80055 p=0.00047221 R2=0.040151 p=0.55469 50 40 30 10 10 0 ε2 0.2 0.4 R2=0.81935 p=0.0001269 R2=0.83671 p=0.00020854 R2=0.4973 p=0.015354 30 20 −0.2 0 −0.1 0.6 −0.05 (c) Second principal strain Axons/mm 50 70 R2=0.94548 p=5.4617e−007 60 2 R =0.89674 p=3.2461e−005 50 30 0.2 R2=0.33307 p=0.063033 2 R =0.51866 p=0.018824 30 20 10 10 0 ε3 0.2 (e) Third principal strain 0.4 0.6 R =0.3638 p=0.049476 40 20 −0.2 0.15 2 R2=0.4764 p=0.018731 40 0 −0.4 0.05 0.1 σ2 (MPa) Scatter plot Maikos_WM_LC: σ3 vs. cell count Axons/mm 60 0 (d) Second principal stress Scatter plot Maikos_WM_LC: ε3 vs. cell count 70 0.2 40 20 0 −0.4 0.15 Scatter plot Maikos_WM_LC: σ2 vs. cell count 70 Axons/mm 60 0.05 0.1 σ1 (MPa) (b) First principal stress Scatter plot Maikos_WM_LC: ε2 vs. cell count 70 0 0 −0.1 −0.05 0 0.05 0.1 σ3 (MPa) 0.15 0.2 (f ) Third principal stress Figure E.3: Parameter correlations – white matter lateral column region (Maikos properties) 117 E.2. Regional correlations Scatter plot Maikos_WM_VM: ε1 vs. cell count 70 60 Scatter plot Maikos_WM_VM: σ1 vs. cell count 70 R2=0.90146 p=7.9868e−006 60 R2=0.82001 p=0.00031034 2 R =0.27144 p=0.10031 50 Axons/mm Axons/mm 50 40 30 10 10 0 ε1 0.2 0.4 R2=0.069772 p=0.43252 30 20 −0.2 R2=0.87655 p=6.6912e−005 40 20 0 −0.4 R2=0.88405 p=1.6739e−005 0 −0.1 0.6 −0.05 (a) First principal strain Axons/mm 50 R2=0.94084 p=7.8991e−007 60 R2=0.10827 p=0.35322 R2=0.81654 p=0.00013622 50 40 30 10 10 0 ε2 0.2 0.4 R2=0.73323 p=0.00076605 R2=0.41691 p=0.043742 R2=0.39155 p=0.039468 30 20 −0.2 0 −0.1 0.6 −0.05 (c) Second principal strain Axons/mm 50 70 R2=0.9496 p=3.8271e−007 60 2 R =0.48932 p=0.024348 50 30 0.2 R2=0.89772 p=9.4577e−006 2 R =0.26013 p=0.13204 30 20 10 10 0 ε3 0.2 (e) Third principal strain 0.4 0.6 R =0.64803 p=0.0027968 40 20 −0.2 0.15 2 R2=0.56162 p=0.0079282 40 0 −0.4 0.05 0.1 σ2 (MPa) Scatter plot Maikos_WM_VM: σ3 vs. cell count Axons/mm 60 0 (d) Second principal stress Scatter plot Maikos_WM_VM: ε3 vs. cell count 70 0.2 40 20 0 −0.4 0.15 Scatter plot Maikos_WM_VM: σ2 vs. cell count 70 Axons/mm 60 0.05 0.1 σ1 (MPa) (b) First principal stress Scatter plot Maikos_WM_VM: ε2 vs. cell count 70 0 0 −0.1 −0.05 0 0.05 0.1 σ3 (MPa) 0.15 0.2 (f ) Third principal stress Figure E.4: Parameter correlations – white matter ventro-medial region (Maikos properties) 118 E.2. Regional correlations Scatter plot Maikos_WM_VL: ε1 vs. cell count 70 60 Scatter plot Maikos_WM_VL: σ1 vs. cell count 70 R2=0.90878 p=5.6242e−006 60 R2=0.61175 p=0.0075047 2 R =0.18879 p=0.18174 50 Axons/mm Axons/mm 50 40 30 10 10 0 ε1 0.2 0.4 R2=0.13639 p=0.26367 30 20 −0.2 R2=0.8011 p=0.00046694 40 20 0 −0.4 R2=0.87732 p=2.1645e−005 0 −0.1 0.6 −0.05 (a) First principal strain Axons/mm 50 R2=0.86485 p=3.3649e−005 60 R2=0.11132 p=0.3461 R2=0.60015 p=0.0051118 50 40 30 10 10 0 ε2 0.2 0.4 R2=0.81721 p=0.00013395 R2=0.43048 p=0.039378 R2=0.3974 p=0.037597 30 20 −0.2 0 −0.1 0.6 −0.05 (c) Second principal strain Axons/mm 50 70 R2=0.9139 p=4.3264e−006 60 2 R =0.35992 p=0.066725 50 30 0.2 R2=0.81353 p=0.00014677 2 R =0.19075 p=0.20694 30 20 10 10 0 ε3 0.2 (e) Third principal strain 0.4 0.6 R =0.53035 p=0.011043 40 20 −0.2 0.15 2 R2=0.65107 p=0.002685 40 0 −0.4 0.05 0.1 σ2 (MPa) Scatter plot Maikos_WM_VL: σ3 vs. cell count Axons/mm 60 0 (d) Second principal stress Scatter plot Maikos_WM_VL: ε3 vs. cell count 70 0.2 40 20 0 −0.4 0.15 Scatter plot Maikos_WM_VL: σ2 vs. cell count 70 Axons/mm 60 0.05 0.1 σ1 (MPa) (b) First principal stress Scatter plot Maikos_WM_VL: ε2 vs. cell count 70 0 0 −0.1 −0.05 0 0.05 0.1 σ3 (MPa) 0.15 0.2 (f ) Third principal stress Figure E.5: Parameter correlations – white matter ventro-lateral region (Maikos properties) 119 E.2. Regional correlations Scatter plot Maikos_GM: ε1 vs. cell count 70 60 Scatter plot Maikos_GM: σ1 vs. cell count 70 R2=0.9314 p=1.5446e−006 60 R2=0.93132 p=6.2576e−006 2 R =0.47121 p=0.019659 50 % of cells % of cells 50 40 30 10 10 0 ε1 0.2 0.4 R2=0.23626 p=0.12955 30 20 −0.2 R2=0.5005 p=0.022108 40 20 0 −0.4 R2=0.95834 p=1.6184e−007 0 −0.1 0.6 −0.05 (a) First principal strain % of cells 50 R2=0.78539 p=0.00028016 60 R2=0.85825 p=0.00011727 R2=0.71094 p=0.0011127 50 40 30 10 10 0 ε2 0.2 0.4 R2=0.50224 p=0.014634 R2=0.71675 p=0.002004 R2=0.68202 p=0.0017365 30 20 −0.2 0 −0.1 0.6 −0.05 (c) Second principal strain % of cells 50 70 R2=0.93926 p=8.901e−007 60 2 R =0.9236 p=9.6174e−006 50 30 0.2 R2=0.41352 p=0.032818 2 R =0.92466 p=9.0869e−006 30 20 10 10 0 ε3 0.2 (e) Third principal strain 0.4 0.6 R =0.72978 p=0.00081316 40 20 −0.2 0.15 2 R2=0.70028 p=0.0013171 40 0 −0.4 0.05 0.1 σ2 (MPa) Scatter plot Maikos_GM: σ3 vs. cell count % of cells 60 0 (d) Second principal stress Scatter plot Maikos_GM: ε3 vs. cell count 70 0.2 40 20 0 −0.4 0.15 Scatter plot Maikos_GM: σ2 vs. cell count 70 % of cells 60 0.05 0.1 σ1 (MPa) (b) First principal stress Scatter plot Maikos_GM: ε2 vs. cell count 70 0 0 −0.1 −0.05 0 0.05 0.1 σ3 (MPa) 0.15 0.2 (f ) Third principal stress Figure E.6: Parameter correlations – gray matter region (Maikos properties) 120 Appendix F Embedded 3D Model of Segmented Geometry of the Rat Cervical Spine Click on the image below to enter the 3D model interface. Elements of the model may be highlighted by clicking on them, and can be hidden or made transparent by using the model tree in the navigation panel at the left. The model may be rotated by clicking and holding the left mouse button, and then moving the mouse. Click and hold the right mouse button to zoom. For other options, see the 3D toolbar above the model. This model requires Adobe Acrobat Reader version 7 or higher, available at www.adobe.com/products/acrobat/readstep2.html. Click here to cycle through a set of predefined views. 121
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A nonlinear finite element model of the rat cervical...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A nonlinear finite element model of the rat cervical spine : validation and correlation with histological… Russell, Colin Macdonald 2012
pdf
Page Metadata
Item Metadata
Title | A nonlinear finite element model of the rat cervical spine : validation and correlation with histological measures of spinal cord injury |
Creator |
Russell, Colin Macdonald |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | Researchers and clinicians do not currently use the heterogeneity of the primary mechanism of spinal cord injury (SCI) to tailor treatment strategies because the effects of these distinct patterns of acute mechanical damage on long-term neuropathology have not been fully investigated. Computational modelling of SCI enables the analysis of mechanical forces and deformations within the spinal cord tissue that are not visible experimentally. I created a dynamic, hyperviscoelastic three-dimensional finite element (FE) model of the rat cervical spine and simulated contusion and dislocation SCI mechanisms. I investigated the relationship between maximum principal strain and previously published tissue damage patterns, and compared primary injury patterns between mechanisms. My model incorporates the spinal cord white and gray matter, dura mater, cerebrospinal fluid, spinal ligaments, intervertebral discs, a rigid indenter and vertebrae, and failure criteria for ligaments and vertebral endplates. High-speed (1 m/s) contusion and dislocation injuries were simulated between vertebral levels C3 and C6 to match previous animal experiments, and average peak maximum principal strains were calculated for several regions at the injury epicentre and at 1 mm intervals from +5 mm rostral to -5 mm caudal to the lesion. I compared average peak principal strains to tissue damage measured previously via axonal permeability to 10 kD fluorescein-dextran (Choo, 2007). Linear regression of tissue damage against peak maximum principal strain for pooled data within white matter regions yields significant (p < 0.0001) correlations that are similar for both contusion (R² = 0.86) and dislocation (R² = 0.54). With additional simulations of cord contusion injuries at lower injury velocities of 3 and 300 mm/s, I found that current material properties used to model the cord are not biofidelic within this velocity range. By fitting existing experimental cord material testing data and plotting alongside the material properties used in several related models, I further demonstrated the remaining divide between experimental data and computational models. My model enhances our understanding of the differences in injury patterns between SCI mechanisms, and provides further evidence for the link between principal strain and tissue damage. Furthermore, my results speak to a continued need to test cord material properties at a range of strains and strain rates to better refine cord hyperviscoelastic properties. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-08-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0073131 |
URI | http://hdl.handle.net/2429/43119 |
Degree |
Master of Applied Science - MASc |
Program |
Biomedical Engineering |
Affiliation |
Applied Science, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_fall_russell_colin.pdf [ 23.66MB ]
- Metadata
- JSON: 24-1.0073131.json
- JSON-LD: 24-1.0073131-ld.json
- RDF/XML (Pretty): 24-1.0073131-rdf.xml
- RDF/JSON: 24-1.0073131-rdf.json
- Turtle: 24-1.0073131-turtle.txt
- N-Triples: 24-1.0073131-rdf-ntriples.txt
- Original Record: 24-1.0073131-source.json
- Full Text
- 24-1.0073131-fulltext.txt
- Citation
- 24-1.0073131.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0073131/manifest