STUDIES ON THE DEVELOPMENT AND FRACTURE MECHANICS OF CORTICAL BONE By Trevor S. DeLand A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Engineering Mechanics – Master of Science 2013 ABSTRACT STUDIES ON THE DEVELOPMENT AND FRACTURE MECHANICS OF CORTICAL BONE By Trevor S. DeLand Bone is similar to other engineering materials in that it is a composite tissue consisting of an organic collagen matrix, a hydroxyapatite mineral phase, and water. This connective tissue is unique in that it is constantly remodeling its structure in response to applied loads. However, high mechanical stresses resulting from trauma can exceed the bone’s material strength and cause fracture. Due to bone’s unique structure, fractures often form distinct patterns and forensic scientists can analyze these patterns to determine the cause of an injury. This thesis explores two major topics in bone biomechanics. First, the effect of exercise on skeletal mechanical properties is studied, focusing on the contributions of material and geometrical properties to whole-bone strength and stiffness. Using these data, a model is developed to provide a predictive metric for whole-bone mechanical properties from a non-invasive computed tomography scan. Second, the fracture patterns of different bones under various traumatic scenarios are studied. One study advances the development of an infant porcine model for fracture in the developing human skull. In particular, the effects of fall height and impact interface on fracture severity are examined. A second study explores the fracture patterns generated from bending failure of long bones and introduces a novel mechanism for the formation of these fractures. Copyright by TREVOR S. DELAND 2013 ACKNOWLEDGMENTS I would like to begin by thanking my advisor, Dr. Roger Haut, for providing me the opportunity to perform this work and the guidance to become a better researcher. I would also like to thank everyone at the OBL whose help and friendship made the past two years much more enjoyable, especially Cliff, Keith, Mark, and Jerrod. As ever, my thanks to my family to whom I owe all that I have accomplished. Finally, my deepest gratitude to Abby: whose unwavering love and support have gotten me through the bad times and helped me to enjoy the good. Thank you all. iv TABLE OF CONTENTS LIST OF TABLES ................................................................................................................... vii LIST OF FIGURES ................................................................................................................ viii CHAPTER ONE AN INTRODUCTION TO BONE ............................................................................................. 1 REFERECES..................................................................................................................... 9 CHAPTER TWO MECHANICAL PROPERTIES OF PULLET BONES FROM DIFFERENT HOUSING SYSTEMS ............................................................................................................................ 11 INTRODUCTION ........................................................................................................... 11 METHODS .................................................................................................................... 12 RESULTS....................................................................................................................... 17 DISCUSSION................................................................................................................. 20 ACKNOWLEDGEMENTS ............................................................................................... 24 APPENDIX .................................................................................................................... 25 REFERENCES ................................................................................................................ 27 CHAPTER THREE DEVELOPMENT OF A MODEL FOR IN VIVO ASSESSMENT OF WHOLE-BONE STRENGTH AND STIFFNESS ............................................................................................... 31 INTRODUCTION ........................................................................................................... 31 METHODS .................................................................................................................... 33 Finite Element Simulation for Inverse Analyses ................................................... 34 Quantitative CT and Linear Regression ................................................................. 38 RESULTS....................................................................................................................... 40 DISCUSSION................................................................................................................. 46 APPENDIX .................................................................................................................... 51 REFERENCES ................................................................................................................ 54 CHAPTER FOUR THE EFFECTS OF FALL HEIGHT AND INTERFACE COMPLIANCE ON FRACTURE PATTERNS IN THE DEVELOPING PORCINE SKULL .............................................................. 64 INTRODUCTION ........................................................................................................... 64 METHODS .................................................................................................................... 67 RESULTS....................................................................................................................... 74 Energy Effects ....................................................................................................... 74 Interface Effects .................................................................................................... 77 v Fracture Distribution and GIS ............................................................................... 80 Finite Element Modeling ....................................................................................... 84 DISCUSSION................................................................................................................. 85 ACKNOWLEDGEMENTS ............................................................................................... 91 APPENDIX .................................................................................................................... 92 REFERENCES ................................................................................................................ 97 CHAPTER FIVE PATTERNS OF FRACTURE DUE TO THREE-POINT BENDING IN LONG BONES ................. 101 INTRODUCTION ......................................................................................................... 101 METHODS .................................................................................................................. 103 RESULTS..................................................................................................................... 106 DISCUSSION............................................................................................................... 109 APPENDIX .................................................................................................................. 114 REFERENCES .............................................................................................................. 121 CHAPTER SIX CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK .................................... 125 vi LIST OF TABLES Table 2.1. Mechanical, material, and geometrical properties of the tibia for each housing condition. AV = aviary system, CC = conventional cage.......................... 26 Table 2.2. Mechanical, material, and geometrical properties of the humerus for each housing condition. AV = aviary system, CC = conventional cage. ................ 26 Table 3.1. Raw data from Chapter 3. ................................................................................ 59 Table 4.1. Revisited raw data from low-energy drops to a rigid surface (Powell et al. 2013). ............................................................................................................... 93 Table 4.2. Raw data from high energy drop impacts to a rigid surface............................ 94 Table 4.3. Raw data from high-energy drop impacts to a compliant surface (Carpet 1). ............................................................................................................. 95 Table 4.4. Raw data from high-energy drop impacts to a compliant surface (Carpet 2). ............................................................................................................. 96 Table 5.1. Comparison of fracture characteristics for paired bones. Left femora were loaded on the anterior surface and right femora were loaded on the posterior surface. When identifying the gross fracture type: T=transverse, O=oblique ............................................................................................................ 108 vii LIST OF FIGURES Figure 1.1. Cross-sectional view of a bone showing the sponge-like lattice of trabecular bone surrounded by dense cortical bone. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this thesis. .................................................................. 2 Figure 1.2. Cross-section diagram of a long bone showing the three regions and their composition .................................................................................................... 3 Figure 1.3. Anatomy of the infant (top) and adult (bottom) human skull ......................... 4 Figure 1.4. Fracture patterns of long bones resulting from tensile, compressive, torsion, and bending loads. Redrawn from Carter and Spengler (1982). .............. 6 Figure 2.1. Fixture used to align and secure test specimens for potting ......................... 13 Figure 2.2. Four-point bending fixture with pivoting cups holding the potted bone ends. Vertical displacement of the hydraulic actuator acted through the pivoting crossbar (rear) on the cups, causing rotation of the bone ends. ........... 14 Figure 2.3. Four-point bending mechanical test setup and geometry ............................. 16 Figure 2.4. Beam theory scenario and deflections of the test section ends from the four-point bending setup in Figure 2.3........................................................... 17 Figure 2.5. Differences in geometrical properties between housing systems. Tibiae and humeri from AV pullets had a larger cross-sectional area and second moment of inertia than CC pullets. A * denotes a significant difference between means (P<0.05). ................................................................... 18 Figure 2.6. Representative moment-rotation curves showing the mechanical behavior of the bones in bending up to failure of AV (dotted) and CC (solid) bones. Stiffness was determined from the slope of the initial, linear portion of the curves. ................................................................................. 19 viii Figure 2.7. Mechanical testing differences between housing systems. AV tibiae and humeri withstood higher moments and were stiffer than CC bones. A * denotes a significant difference between means (P<0.05). .............................. 19 Figure 2.8. Comparison of bone material properties calculated from the geometry and mechanical data. A * denotes a significant difference between means (P<0.05). ..................................................................................... 20 Figure 2.9. Bone cross-sections reconstructed using average measurements for AV (dashed line) and CC (solid line) humeri and tibiae ........................................ 22 Figure 3.1. Beam theory scenario used to model the current study showing the application of moments to the bone test section ends and the resultant deformations, measured by the rotation at the ends. ......................................... 36 Figure 3.2. Sample rotation-moment curve showing the selection of the moment used for FE analysis loading conditions (O) based on half the bending moment at failure (X). For all specimens, this point fell within the linear region of the moment-rotation curve. ................................................................. 41 Figure 3.3. a) Topographical contour plot of the natural log of the objective function, Z, Equation (3.2), (out of the page) versus A and B. Lower “elevations” (dark blue) indicate that FE models using that correlation more closely matched experimental observations. It appeared that there were a number of correlations lying along a line that minimized the function (traced by a white line). b) Plot of the objective function, Z, versus A along the white line showing the concavity of the solution. The A and B combination that minimized the objective function (resulted in a model that most closely matched the experimental result) is identified by an X........................................................................................................................ 42 Figure 3.4. Correlation between the stiffness determined from mechanical tests and the stiffness of the FE bone models generated from CT using Equation (3.8). The same relationship was used for all bones in the study (both humerus and tibia). ..................................................................................... 44 Figure 3.5. a) Linear fit of bending stiffness versus the metric generated from QCT and beam theory following Equation (3.5). b) Linear fit of failure ix moment versus the section modulus measured from QCT following Equation (3.7)........................................................................................................ 45 Figure 4.1. Experimental setup of the drop tower and porcine head specimen ............. 68 Figure 4.2. Photo detail of the floating rod and specimen mounting apparatus attached to the drop trolley ................................................................................. 69 Figure 4.3. a) The specimen and drop trolley were suspended above the impact interface and held in place by solenoid clamps. b) When released, the trolley and specimen mounting rod fell together until the trolley struck the padded stop. The floating rod and specimen continued downward, striking the instrumented impact interface.......................................................... 71 Figure 4.4. Right-side and occipital-oriented views used for marking fracture patterns. The impact site is indicated by an X ...................................................... 72 Figure 4.5. Domed shell used to model impact loading on the right parietal. A pressure was applied to the center of the outer surface with constraints modeling the rim resting on a frictionless table. .................................................. 74 Figure 4.6. Peak impact force versus age for the current (high-energy) and previous (low-energy) studies .............................................................................. 75 Figure 4.7. Time to peak load versus age for high and low energy impacts .................... 76 Figure 4.8. Total bone and diastatic fracture length versus age for high and low energy impacts to a rigid surface .......................................................................... 77 Figure 4.9. Peak impact force versus age for high energy impacts to three different interfaces. Carpet 1 = thin loop pile commercial carpet. Carpet 2 = cut pile carpet with pad underlay ...................................................................... 78 Figure 4.10. Time to peak load versus age for high-energy impacts to three different interfaces. Carpet 1 = thin loop pile commercial carpet. Carpet 2 = cut pile carpet with pad underlay ...................................................................... 79 x Figure 4.11. Total bone and diastatic fracture length versus age for high energy impacts to three different interfaces. Carpet 1 = thin loop pile commercial carpet. Carpet 2 = cut pile carpet with pad underlay ........................................... 80 Figure 4.12. Geographic Information System maps of the a) low energy rigid b) high energy rigid c) high energy carpet 1 d) high energy carpet 2 impacts ......... 82 Figure 4.13. Maximum principal stresses generated in the shell as a result of the boundary and loading conditions of Figure 4.5. Arrows indicate the relative magnitude and orientation of the maximum principal stress with larger stresses signified by large, red arrows and smaller stresses by short, blue arrows. The largest tensile stresses were formed around the periphery of the shell. In addition, tensile stress increased directly under the central load application area. ......................................................................... 84 Figure 5.1. Diagram of a typical butterfly fracture of a long bone as a result of bending loads. ..................................................................................................... 102 Figure 5.2. Technique used to pot the bone ends with the femur passing through the bottom cup to leave the knee intact. ........................................................... 104 Figure 5.3. Three-point bending setup on the servohydraulic testing machine. ........... 105 Figure 5.4. Oblique (a) and Transverse (b) gross fracture patterns. The anvil impacted the top surface of the bone at the site marked by the arrow, initiating tensile failure on the bottom surface. Notice that the hairline fractures branch away from the gross fracture on the tensile side (bottom half) of the bone. ................................................................................................ 107 Figure 5.5. Sample load-displacement plot of a three-point bending test of a femur up to failure .............................................................................................. 108 Figure 5.6. Fracture profile of compact tensile test specimen with an initial crack made transverse to the long axis of the bone. Redrawn from Behiri and Bonfield, 1989 ..................................................................................................... 110 Figure 5.7. Composition of cortical bone showing the orientation of osteons parallel to the long axis. ...................................................................................... 112 xi Figure 5.8. Specimen 12-1441L-AP. Top: lateral view. Bottom: medial view. For all photographs, the bone was oriented with the anvil striking the top surface................................................................................................................. 115 Figure 5.9. Specimen 12-1441R-PA. Top: lateral view. Bottom: medial view. ............... 115 Figure 5.10. Specimen 12-1447L-AP. Top: lateral view. Bottom: medial view. ............. 116 Figure 5.11. Specimen 12-1447R-PA. Top: lateral view. Bottom: medial view. ............. 116 Figure 5.12. Specimen 12-1466L-AP. Top: lateral view. Bottom: medial view. ............. 117 Figure 5.13. Specimen 12-1466R-PA. Top: lateral view. Bottom: medial view. ............. 117 Figure 5.14. Specimen 13-1026L-AP. Top: lateral view. Bottom: medial view. ............. 118 Figure 5.15. Specimen 13-1026R-PA. Top: lateral view. Bottom: medial view. ............. 118 Figure 5.16. Specimen 13-1065L-AP. Top: lateral view. Bottom: medial view. ............. 119 Figure 5.17. Specimen 13-1065R-PA. Top: lateral view. Bottom: medial view. ............. 119 Figure 5.18. Specimen 13-1180L-AP. Top: lateral view. Bottom: medial view. ............. 120 Figure 5.19. Specimen 13-1180R-PA. Top: lateral view. Bottom: medial view. ............. 120 xii CHAPTER ONE AN INTRODUCTION TO BONE The primary function of the skeletal system is to provide a rigid structure to support the body. The skeleton keeps the body upright against the pull of gravity and provides a frame of articulating linkages that work with muscles to produce a wide range of motions. Owing to its inherently mechanical nature and close parallels to other engineering structures, bone is a popular subject of biomechanics research. Similar to other composite materials, bone is composed of a matrix of fibers, which give flexibility and tensile strength, and a mineral phase that stiffens the bone and provides compressive strength. The matrix is organic, composed largely of Type I collagen; and the mineral phase consists almost entirely of hydroxyapatite. Additionally, bone contains a large amount of water, some of which is bonded to collagen fibers and some of which is free to flow in the matrix. As a result of its composite structure, bone is stronger in compression than it is in tension (Reilly and Burstein 1975). In addition, the presence of water changes the response of bone tissue at different strain rates, a phenomenon known as viscoelasticity (Carter and Hayes 1976). On the macroscopic scale, bone is classified into two general categories based on porosity. Highly porous bone has an architecture that resembles a sponge and is commonly referred to as trabecular, (also spongy or cancellous) bone. The second, more 1 dense type of bone has a linear, organized structure and is referred to as cortical (compact) bone (Figure 1.1). Figure 1.1. Cross-sectional view of a bone showing the sponge-like lattice of trabecular bone surrounded by dense cortical bone. For interpretation of the references to color in this and all other figures, the reader is referred to the electronic version of this thesis. Long bones of the appendicular skeleton contain both cortical and trabecular bone. Around the outer edge of these bones is a layer of cortical bone. Inside this shell at either end of the bone is a region of trabecular bone. In the middle section of the bone, the cortex thickens and, in the absence of trabecular bone, forms a hollow tube. Thus, long bones can be divided into three distinct regions based on their general composition. The sections containing trabecular bone at either end are called epiphyses and the middle, tubular portion of solid cortical bone is called the diaphysis (Figure 1.2). 2 Figure 1.2. Cross-section diagram of a long bone showing the three regions and their composition The skull is unique in that it – in addition to providing structural support – serves as a protective shell for the brain. In early developmental stages, the skull begins as a series of bony plates that are connected by soft tissue at joints called sutures forming gaps known as fontanels. As the skull grows and matures, these sutures narrow and harden and the bones become fused and interlocked, forming a solid structure (Figure 1.3). Because of these developmental changes, the mechanical properties of the skull vary significantly between infants and adults (Margulies and Thibault 2000). 3 Figure 1.3. Anatomy of the infant (top) and adult (bottom) human skull A unique feature of bone is that it is a dynamic structure; it is constantly being “turned over” by metabolic cells that resorb existing bone (osteoclasts) and deposit new bone (osteoblasts). This reconstruction process is ongoing throughout the life of the bone and is thought to occur as an automated remodeling/repair mechanism for small cracks that, if not repaired, might lead to whole-bone failure (Cowin 1989). In addition to this passive remodeling process, external mechanical stimuli have also been found to trigger bone remodeling. Woo et al. (1981) documented that, after a year of regimented exercise, the femoral diaphysis of pigs became significantly larger with the bone growing inward, shrinking the medullary cavity. In another study (Woo et al. 1976), mechanical 4 loading was effectively removed from intact dog femora by means of rigid internal fixation. It was documented that unloading the bone in this manner resulted in significant atrophy and cortical thinning. Through these studies, it has become clear that the degree of mechanical loading on a bone can induce different modeling pathways: depositing new bone to help support increased loads and resorbing existing bone where structure is no longer needed. While intermediate mechanical loads can activate modeling pathways to further develop bone, high forces generated in traumatic scenarios can stress the bone to levels that exceed its strength resulting in bone fracture. Due to its composite viscoelastic structure, the way in which bone fractures varies with both orientation and rate of the applied load. As a result, the pattern of fracture in long bones can be used to determine which type of loading likely caused the fracture (Carter and Spengler 1982) (Figure 1.4). Also, for certain fracture patterns, it is even possible to identify the direction in which the load was applied (Symes et al. 2012). In addition, the degree of comminution or fracture fragmentation can give some idea of the applied energy or rate of loading causing the fracture. 5 Figure 1.4. Fracture patterns of long bones resulting from tensile, compressive, torsion, and bending loads. Redrawn from Carter and Spengler (1982). Similar techniques can be used to assess the cause of fracture in the skull based on the fracture pattern. As is the case in other bones, an increase in impact energy to the skull has been found to significantly increase the degree of fracture comminution (Powell et al. 2012). Such information is of great importance in the forensics community where the patterns of skull fracture in pediatric cases can be examined to suggest a mechanism of injury and determine whether the trauma was accidental or inflicted (Hobbs 1984). The research presented in this thesis explores the biomechanics of bone; studying how bones grow and change in response to exercise and mechanical loading as well as the effect of impact loading conditions on the fracture process. In Chapter 2, mechanical testing is used to assess the stiffness and strength of the humerus and tibia from young hens reared in different environments. It was hypothesized that the housing system that allowed more room for movement and exercise would result in stiffer and stronger 6 bones. Chapter 3 builds on the work of Chapter 2 and describes the development of models using finite element and linear regression analyses in conjunction with computed tomography (CT) scan data to predict the mechanical characteristics of bone. The study proposes a method for in vivo assessments of bone strength and stiffness through the use of high-resolution micro-CT data and predictions from FE and linear regression analyses based on simple beam theory. Chapter 4 involves the study of bone fracture through a continuation of previous work in the Orthopaedic Biomechanics Laboratories examining fracture patterns in the developing skull using an infant porcine surrogate model. Specifically, the effects of fall height and impact interface compliance on the patterns of cranial fracture are evaluated as a function of specimen age. It was hypothesized that increasing the impact surface compliance would have an effect similar to reducing the fall height. Chapter 5 examines another mechanism of bone fracture, performing controlled 3-point bending tests on human femora. The hypothesis of this study was that a specific pattern of fracture would develop that could be used to determine the direction of the applied lateral load to the long bone. This thesis provides new insight across multiple disciplines within the field of bone biomechanics. The work presented in Chapters 2 and 3 will help researchers better understand the factors that are responsible for whole bone characteristics such as stiffness and strength. In addition, a model is introduced for predicting these characteristics that may be of value in conducting longitudinal studies of bone development and remodeling under various in vivo loading conditions. The study of fracture patterns presented in Chapters 4 and 5 may be of significant benefit in the 7 forensics community by furthering the understanding and some of the science behind bone fracture mechanisms and how the information can be used to determine injury causation. 8 REFERENCES 9 REFERENCES Carter, D. R., and W. C. Hayes. 1976. "Bone compressive strength - influence of density and strain rate". Science 194 (4270):1174-1176. Carter, DR, and DM Spengler. 1982. "Biomechanics of fracture". In Bone in Clinical Orthopaedics, edited by G. Sumner-Smith. Philadelphia: Saunders. Cowin, Stephen Corteen. 1989. Bone Mechanics. Boca Raton, FL: CRC Press. Hobbs, C. J. 1984. "Skull fracture and the diagnosis of abuse". Archives of Disease in Childhood 59 (3). Margulies, S. S., and K. H. Thibault. 2000. "Infant skull and suture properties: Measurements and implications for mechanisms of pediatric brain injury". Journal of Biomechanical Engineering 122 (4):364-371. Powell, Brian J., Nicholas V. Passalacqua, Timothy G. Baumer, Todd W. Fenton, and Roger C. Haut. 2012. "Fracture patterns on the infant porcine skull following severe blunt impact". Journal of Forensic Sciences 57 (2). Reilly, Donald T, and Albert H Burstein. 1975. "The elastic and ultimate properties of compact bone tissue". Journal of Biomechanics 8 (6):393-405. Symes, Steven A, Ericka N L'Abbé, Erin N Chapman, Ivana Wolff, and Dennis C Dirkmaat. 2012. "Interpreting Traumatic Injury to Bone in Medicolegal Investigations". In A Companion to Forensic Anthropology, edited by D. C. Dirkmaat: Blackwell Publishing Ltd. Woo, SL, WH Akeson, RD Coutts, LADD Rutherford, DAVID Doty, GF Jemmott, and DAVID Amiel. 1976. "A comparison of cortical bone atrophy secondary to fixation with plates with large differences in bending stiffness". The Journal of bone and joint surgery. American volume 58 (2):190. Woo, SL, Steve C Kuei, David Amiel, MA Gomez, WC Hayes, FC White, and WH Akeson. 1981. "The effect of prolonged physical training on the properties of long bone: a study of Wolff's Law". The Journal of Bone and Joint Surgery American Volume 63 (5):780. 10 CHAPTER TWO MECHANICAL PROPERTIES OF PULLET BONES FROM DIFFERENT HOUSING SYSTEMS INTRODUCTION Commercial chicken egg laying operations typically house birds in blocks of cages designed to optimize population density and maximize egg production. Recently, the welfare of the animals housed in these systems has been questioned and new, alternative systems have been introduced to alleviate these concerns (Jendral et al. 2008). One such alternative is referred to as a “cage-free” aviary system. Both aviary and conventional systems provide the housed birds with ad libitum access to feed and fresh water. Conventional cages typically house birds on a single level at a high population density and provide no opportunity for enrichment (Tauson 1998). Aviary systems provide an open, communal setting with more opportunity for movement and exercise. Hens can forage and dust bathe on the ground, have access to perches, and lay eggs in isolated nest boxes (Lay et al. 2011). Previous studies have examined the effect of such alternative housing systems on laying hen behavior and welfare (Jendral et al. 2008; Leyendecker et al. 2002; Leyendecker et al. 2005; Newman and Leeson 1998). A common point of interest for these studies has been the determination of whole-bone strength to assess predisposition to fracture. This testing is typically done by three-point bending failure of long bones, usually humeri and tibiae. Such studies have found that birds from aviary systems develop 11 higher whole-bone strength (i.e., tolerate higher bending loads before failing) than bones from birds raised in conventional cages (Newman and Leeson 1998; Leyendecker et al. 2002). Laying hen bone quality is most commonly assessed by correlating measured ash content or CT density with strength and Young’s modulus (Currey 1969; Currey 1988; Carter and Hayes 1976; Leyendecker et al. 2002; Leyendecker et al. 2005). Some studies on laying hens have used whole-bone mechanical test results in combination with beam theory equations to approximate the Young’s modulus and ultimate strength of the bone. These studies found that there is little to no difference in material properties between housing conditions. Instead it was found that mechanical property differences are attributable changes in cortical geometries (Newman and Leeson 1998; Rath et al. 2000). Expanding on this work, the current study examined the mechanical properties of humeri and tibiae from 16-week old pullets reared in conventional cages and aviary systems using whole-bone four-point bending tests. The hypothesis of the current study was that bones from birds raised in an aviary system would be stiffer and fail at higher loads than bones from conventional cage birds. Changes in mechanical properties were expected to result only from changes in bone geometry with material properties remaining constant between housing systems. METHODS At 16 weeks of age, 120 white leghorn pullets from each housing system – conventional cage (CC) and aviary (AV) – were killed by cervical dislocation and stored at -20°C. Two 12 days prior to mechanical testing, left legs and wings were thawed at room temperature. The tibiae and humeri were harvested and cleaned of all soft tissue. The bones were wrapped in saline soaked gauze and kept moist throughout all preparations and testing procedures. A uniform mid-diaphysis section – 20 mm long for the humerus and 30 mm for the tibia – was used in the tests and the remaining ends of the bones were potted in cups with polyester resin (Martin Senour Fibre Strand Plus 6371, Sherwin-Williams; Cleveland, OH). A custom rig secured the bones in alignment with the cups while the resin cured (Figure 2.1). After potting, the outer anterior-posterior (AP) and mediallateral (ML) dimensions of the bone were measured with digital calipers (Absolute IP 66, Mitutoyo Corp.; Kawasaki, Kanagawa, Japan) at the ends and center of the test section. Potting Mold Cups Bone Test Section Positioning Clamp Polyester Resin Figure 2.1. Fixture used to align and secure test specimens for potting The potted specimens were installed in a four-point bending fixture mounted in a servohydraulic testing machine (model 1331, Instron; Norwood, MA). Potted bone ends were installed into pivoting cups. A balanced force was applied to each cup from the actuator through a pivoting crossbar (Figure 2.2). The setup applied equal bending 13 moments to each end of the specimen, uniformly loading the test section in pure bending. An actuator preload of 2 N was applied to eliminate residual system compliance before bone failure was induced with a 1 Hz, 10 mm haversine displacement. Load and displacement output of the actuator were recorded at 5000 Hz with a 100-lb load transducer (model 1500ASK-100, Interface; Scottsdale, AZ) and a 6inch linear variable differential transformer mounted on the actuator (model HR 3000, Measurement Specialties; Hampton, VA). The tibiae were oriented with the lateral surface loaded in tension and the humeri with the posterior surface loaded in tension to replicate the loading orientation most likely to cause fracture in vivo. Actuator Crossbar Pivoting Cups Figure 2.2. Four-point bending fixture with pivoting cups holding the potted bone ends. Vertical displacement of the hydraulic actuator acted through the pivoting crossbar (rear) on the cups, causing rotation of the bone ends. After fracture, the bone fragments were reassembled in order to measure anterior, posterior, medial, and lateral cortical thicknesses at the fracture site. Outer dimensions 14 and diaphysial thicknesses were used to approximate the cortical cross-section, assumed to be a hollow ellipse. The material properties of the whole bones were determined based on classical beam theory with the exposed test section modeled as a uniform beam with moments applied to each end. The computations required the bone’s geometrical resistance to bending, called the second moment of area, to be computed using the expression 𝜋𝑎0 𝑏𝑏0 3 𝜋𝑎1 𝑏𝑏1 3 𝐼𝐼 = − 4 4 (2.1) where 𝑎 and 𝑏𝑏 were the radii parallel and perpendicular to the neutral axis of the bone respectively with subscripts 0 and 1 denoting the outer and inner dimensions of the bone respectively (Beer, Johnston, and Dewolf 2002, inside back cover). The applied moment, 𝑀𝑀, and resultant bone-end rotation angle, 𝜃𝜃 (Figure 2.4), were computed from the force-displacement data and fixture geometry using the equations 𝑀𝑀 = 𝜃𝜃 = 𝐹𝐹𝑟𝑟 2 𝑑 𝑟𝑟 (2.2) (2.3) where 𝐹𝐹 was the force measured by the load cell, 𝑑 was the actuator displacement, and 𝑟𝑟 was the distance between the load application and pivot points on the rigid potting cups (Figure 2.3). Whole-bone bending stiffness 𝐾 was determined by the slope of a line fit to the initial, linear portion of the moment-rotation curve. Mechanical stiffness and geometry information were substituted into the beam theory equations to estimate the 15 material properties. Young’s modulus for the whole bone was determined with the equation 𝐾𝐿𝐿 2𝐼𝐼 𝐸𝐸 = (2.4) where L was the length of the test section (Beer, Johnston, and Dewolf 2002, 218) (Figure 2.3). The material strength of each bone was determined based on a computation of the maximum bending stress (force per unit area) in the bone using the expression 𝑆= 𝑀𝑀f 𝑏𝑏0 𝐼𝐼 (2.5) where 𝑀𝑀f was the maximum moment applied at each end of the specimen (Beer, Johnston, and Dewolf 2002, 217) (Figure 2.3). Statistical comparisons of geometrical, mechanical, and material properties of the pullet bones between housing conditions were performed using Student’s t-tests. Results were considered significant for P < 0.05. 𝑟𝑟 𝐹𝐹/2 2𝑏𝑏0 𝐿𝐿 𝐹𝐹/2 𝑟𝑟 Figure 2.3. Four-point bending mechanical test setup and geometry 16 𝑀𝑀 𝜃𝜃 𝐿𝐿 𝑀𝑀 𝜃𝜃 Figure 2.4. Beam theory scenario and deflections of the test section ends from the fourpoint bending setup in Figure 2.3 RESULTS Aviary ML (7.05±0.32 mm) and AP (5.92±0.28 mm) outer dimensions of the humerus were 5.8% and 4.0% larger than those of the CC humeri (6.66±0.28 mm and 5.69±0.22 mm) respectively (P<0.001). Cortical thicknesses of the AV humeri (0.71±0.09 mm) were 33.6% greater than those of the CC humeri (0.53±0.07 mm) (P<0.001). In contrast, there were no significant differences in the ML (P=0.067) and AP (P=0.224) outer bone dimensions of the tibia between housing conditions. Cortical thickness of the AV tibia (0.84±0.07 mm) was, however, on average 14% larger than that of the CC birds (0.74±0.08 mm) (P<0.001). Cross-sectional areas of AV tibiae and humeri were 13.2% and 36.4% greater than those of CC cross-sections respectively (P<0.001) (Figure 2.5). Second moments of area, 𝐼𝐼, of AV tibiae were also 12% greater than those of CC tibiae (P<0.001) (Figure 2.5). Second moments of area of AV humeri were 39.8% greater than those of the CC humeri (P<0.001). 17 10 5 * * 70 AV CC * 60 (mm4) 15 80 CC Second moment of area, 𝐼𝐼 Cross-sectional area (mm2) 20 50 AV * 40 30 20 10 0 0 Tibia Humerus Tibia Humerus Figure 2.5. Differences in geometrical properties between housing systems. Tibiae and humeri from AV pullets had a larger cross-sectional area and second moment of inertia than CC pullets. A * denotes a significant difference between means (P<0.05). An overlay of representative moment-rotation data of a bone from each housing condition illustrates the effect of maximum moment and stiffness on bone mechanics up to failure (Figure 2.6). The AV tibia maximum moment, 𝑀𝑀f , was 12.5% greater than that of the CC tibia (P<0.001) (Figure 2.7). The maximum moment withstood by the AV humerus was 43.9% higher than that of the CC humerus (P<0.001). AV tibia stiffness, 𝐾, was 11.6% greater than that of the CC tibia (P<0.001). The stiffness of the AV humeri were 33.6% greater than that of the CC humeri (P<0.001) (Figure 2.7). 18 4000 Tibia 5000 Humerus Moment (N-mm) Moment (N-mm) 6000 3000 4000 2000 3000 2000 CC Tibia 1000 AV Tibia 1000 AV Humerus CC Humerus 0 0 0 0 0.03 0.06 0.09 0.12 0.15 Rotation (rad) 0.03 0.06 0.09 Rotation (rad) 0.12 Figure 2.6. Representative moment-rotation curves showing the mechanical behavior of the bones in bending up to failure of AV (dotted) and CC (solid) bones. Stiffness was determined from the slope of the initial, linear portion of the curves. 5000 * 4000 40 CC * AV 3000 2000 1000 Stiffness, K (N-mm/rad) Maximum Moment, Mf (N-mm) 6000 0 35 * 30 CC * AV 25 20 15 10 5 0 Tibia Humerus Tibia Humerus Figure 2.7. Mechanical testing differences between housing systems. AV tibiae and humeri withstood higher moments and were stiffer than CC bones. A * denotes a significant difference between means (P<0.05). Aviary tibia bone strength, 𝑆, was 3.7% greater than that of the CC tibia (P=0.012). The strength of the AV humerus was 6.3% greater than that of the CC humerus (P=0.002) (Figure 2.8). There was no significant difference in Young’s modulus, 𝐸𝐸, of the tibia between housing conditions (P=0.38). In contrast, Young’s modulus of CC humeri was 4.6% greater than that of AV humeri (P=0.005) (Figure 2.8). 19 350 CC * 300 250 200 150 100 50 * 20 Young's modulus, 𝐸𝐸 (GPa) Bone strength, S (MPa) 400 AV 15 CC AV * 10 5 0 0 Tibia Humerus Tibia Humerus Figure 2.8. Comparison of bone material properties calculated from the geometry and mechanical data. A * denotes a significant difference between means (P<0.05). DISCUSSION Measurements of bone thickness taken after fracture showed that bones from birds reared in aviary systems developed a thicker cortex than birds from cage systems. Such hypertrophy is a well-established response of bone that has been exposed to increased mechanical loads as a result of exercise (Jones et al. 1977; Woo et al. 1981; Saville and Whyte 1969; Carter 1982). A suggested explanation for this hypertrophy is that bone constantly remodels itself, depositing new bone in response to high strains or induced microdamage (Carter 1984). The findings of the current study suggest that the hypertrophy of aviary-reared pullet bones was a result of the increased opportunity for exercise and therefore additional mechanical loading provided to these birds. Conversely, the sedentary lifestyle of pullets reared in conventional cages may have also contributed to some degree of bone resorption and a reduced cortical cross-section (Whitehead and Fleming 2000). 20 These changes in cross-sectional area were not the same between the tibia and humerus. In the tibia there was no significant difference in the outer (periosteal) dimensions between housing conditions. The larger cross-sectional area of aviary tibiae was due to a reduction in the medullary canal (Figure 2.9). While such inward growth increases the second moment of area, 𝐼𝐼, the addition of bone more proximal to the neutral axis means that second moment of area differences were not as pronounced as cross-sectional area differences. The primary function of endosteal growth may be to increase a bone’s axial rigidity, as has been observed in response to in vivo dynamic longitudinal compressive loading (Robling et al. 2002). Thus, the differences in AV tibial cross-sections in the current study may suggest that these bones experienced more axial loading than CC tibial cross-sections. In contrast, in the AV humerus there was an increased cross-sectional area that was largely due to increased periosteal diameters, with the endosteal dimensions remaining largely unchanged (Figure 2.9). In addition to increasing axial rigidity, this outward expansion of the cortex greatly increased the second moment of area, 𝐼𝐼. As a result, the second moment of area of the humeri increased more dramatically than the cross-sectional area under AV versus CC housing conditions. An increased second moment of area is a characteristic response in bone that has been subjected to additional bending or torsional loads (Bass et al. 2002; Ducher, Daly, and Bass 2009). These findings suggest that changes in loading patterns of the AV system differ between the humerus and tibia, resulting in different growth patterns in each bone. 21 Tibia Humerus Figure 2.9. Bone cross-sections reconstructed using average measurements for AV (dashed line) and CC (solid line) humeri and tibiae Geometrical and mechanical property changes were more dramatic in the humerus than in the tibia. The percentage increases in each parameter of the AV humeri over CC humeri were 2-3 times higher than the changes observed in the tibia. While this is likely due, in part, to different modes of loading, there may be other factors that resulted in unequal growth between bones. Previous studies have shown that cages equipped with short perches yield stronger tibiae but cause little change in the bones of the wing (Whitehead 2004). However, when birds are housed in aviary systems, the opportunity to fly has dramatically increased wing bone breaking force (Knowles and Broom 1990; Fleming et al. 1994). Such was the case in the current study: CC birds had limited opportunity to load their legs but not being able to fly was a stark contrast to birds raised in the open aviary system. The act of flight would dramatically increase the mechanical forces on the humerus more than in the tibia and may have provided the rationale for an increased effect of the AV system on the mechanical and geometrical properties of the AV humeri over those of the CC birds. 22 Bone strength and Young’s modulus were calculated from the geometry and loaddisplacement data to assess the material properties of the bones. These results matched well with literature established values of strength ranging from 200-300 MPa and Young’s modulus averaging 12 GPa (Reilly and Burstein 1975). Although the differences between groups for some of these quantities were statistically significant, the relatively small percentage changes suggest that these differences may not have significant clinical/practical implications. Previous studies have assessed the material properties of exercised bones and found that there is no significant change in their strength or modulus (Woo et al. 1981; Haapasalo et al. 2000). The material property differences identified as statistically significant in the current study were likely due to the large sample sizes that boosted the sensitivity of these statistical tests. Future studies might be conducted with smaller sample sizes and still yield data of practical significance. Multiple simplifying assumptions were made in calculating the material properties of the bones from mechanical test data. Firstly, the beam theory equations used in these calculations approximated the bone geometry as an elliptical cross-section that did not vary along the test section. Although this provides a good approximation for comparison between housing conditions, any potential geometrical variations between conditions might be explored in the future with the use of a computed tomography-generated finite element model simulation of each bone (Niebur et al. 2000). Another limitation of the current study was that the material properties were considered to be constant and isotropic throughout the entire bone. Again, this assumption may provide a good basis for comparison between housing conditions but does not reflect any potential variations 23 that may exist between housing conditions. Finite elements could account for such variations in material properties within each bone, defining the modulus of each element based on the specific CT density of the bone at that point (Wirtz et al. 2000). Additionally, finite elements could be used to model the orthotropic elastic properties of bone, assigning different elastic moduli to different material orientations. However, using such detailed finite element techniques would have been prohibitive for the large sample sizes used in the current study. In conclusion, the results of the current study suggest that the improved load-bearing capability and stiffness of bones from pullets reared in aviary systems were largely due to increased bone quantity, identified by geometry changes that varied between bones, and not significant changes in bone quality. This cortical hypertrophy is a universal response to increased levels of exercise and is likely due to bone-forming mechanisms working to minimize local strains induced by increased activity and mechanical stresses. In addition, the distribution patterns of supplemental bone growth shown between the tibia and the humerus of AV birds may be useful in identifying the predominant mode of loading that was altered by the current AV system of housing. ACKNOWLEDGEMENTS This study is part of the Coalition for a Sustainable Egg Supply (CSES) project. The author would like to thank Mr. Ross Dudgeon for his work developing the method and assisting with preparation and testing and Mr. Cliff Beckett for his technical support and assistance. 24 APPENDIX 25 Table 2.1. Mechanical, material, and geometrical properties of the tibia for each housing condition. AV = aviary system, CC = conventional cage. Mean, (SD) AV CC Failure Moment Rotation (N-m) (rad) Stiffness (N-m/rad) Stress (MPa) Modulus (GPa) Inertia (mm4) Area (mm2) Thickness (mm) Diameter ML AP (mm) (mm) 5144.36 0.16 33.91 316.30 13.69 59.91 14.16 0.84 6.66 5.77 (487.42) (0.02) (3.68) (32.28) (1.21) (8.53) (1.22) (0.07) (0.29) (0.25) 4573.93 0.15 30.39 304.88 13.74 53.51 12.50 0.74 6.58 5.73 (466.65) (0.02) (3.99) (33.53) (1.51) (7.94) (1.25) (0.08) (0.29) (0.30) Table 2.2. Mechanical, material, and geometrical properties of the humerus for each housing condition. AV = aviary system, CC = conventional cage. Mean, (SD) AV CC Failure Moment Rotation (N-m) (rad) Stiffness (N-m/rad) Stress (MPa) Modulus (GPa) Inertia (mm4) Area (mm2) Thickness (mm) Diameter ML AP (mm) (mm) 3636.22 0.12 29.13 244.18 10.25 46.20 12.86 0.71 7.05 5.92 (443.82) (0.02) (4.02) (29.18) (1.21) (7.19) (1.47) (0.09) (0.32) (0.28) 3636.22 0.12 29.13 244.18 10.25 46.20 12.86 0.71 7.05 5.92 (443.82) (0.02) (4.02) (29.18) (1.21) (7.19) (1.47) (0.09) (0.32) (0.28) 26 REFERENCES 27 REFERENCES Bass, SL, Leanne Saxon, RM Daly, CH Turner, AG Robling, E Seeman, and Stephen Stuckey. 2002. "The effect of mechanical loading on the size and shape of bone in pre‐, peri‐, and postpubertal girls: a study in tennis players". Journal of Bone and Mineral Research 17 (12):2274-2280. Beer, Ferdinand P, E Russell Johnston, and John T Dewolf. 2002. Mechanics of Materials. Third ed: McGraw-Hill. Carter, D. R., and W. C. Hayes. 1976. "Bone compressive strength - influence of density and strain rate". Science 194 (4270):1174-1176. Carter, Dennis R. 1984. "Mechanical loading histories and cortical bone remodeling". Calcified Tissue International 36 (1):S19-S24. Carter, DR. 1982. "The relationship between in vivo strains and cortical bone remodeling". Critical Reviews in Biomedical Engineering 8 (1):1. Currey, J. D. 1988. "The effect of porosity and mineral content on the Young's modulus of elasticity of compact bone". Journal of Biomechanics 21 (2):131-139. Currey, JD. 1969. "The mechanical consequences of variation in the mineral content of bone". Journal of Biomechanics 2 (1):1-11. Ducher, Gaele, Robin M Daly, and Shona L Bass. 2009. "Effects of repetitive loading on bone mass and geometry in young male tennis players: a quantitative study using MRI". Journal of Bone and Mineral Research 24 (10):1686-1692. Fleming, RH, CC Whitehead, D Alvey, NG Gregory, and LJ Wilkins. 1994. "Bone structure and breaking strength in laying hens housed in different husbandry systems". British Poultry Science 35 (5):651-662. Haapasalo, H, S Kontulainen, H Sievanen, P Kannus, M Jarvinen, and I Vuori. 2000. "Exercise-induced bone gain is due to enlargement in bone size without a change in volumetric bone density: a peripheral quantitative computed tomography study of the upper arms of male tennis players". Bone 27 (3):351-358. Jendral, M. J., D. R. Korver, J. S. Church, and J. J. R. Feddes. 2008. "Bone mineral density and breaking strength of white leghorns housed in conventional, modified, and commercially available colony battery cages". Poultry Science 87 (5):828-837. 28 Jones, Henry H, James D Priest, Wilson C Hayes, Carol Chin Tichenor, and Donald A Nagel. 1977. "Humeral hypertrophy in response to exercise". Journal of Bone and Joint Surgery American Volume 59 (2):204-208. Knowles, TG, and DM Broom. 1990. "Limb bone strength and movement in laying hens from different housing systems". Veterinary Record 126 (15):354-356. Lay, DC, RM Fulton, PY Hester, DM Karcher, JB Kjaer, JA Mench, BA Mullens, RC Newberry, CJ Nicol, and NP O’Sullivan. 2011. "Hen welfare in different housing systems". Poultry Science 90 (1):278-294. Leyendecker, M., H. Hamann, J. Hartung, G. Glunder, U. Neumann, C. Surie, J. Kamphues, and O. Distl. 2002. "Bone breaking strength, eggshell stability and production traits of laying hens kept in battery cages, furnished cages and an aviary housing system". Paper read at Proceedings of the 7th World Congress on Genetics Applied to Livestock Production, Montpellier, France, August, 2002. Session 4. Leyendecker, M., H. Hamann, J. Hartung, J. Kamphues, U. Neumann, C. Surie, and O. Distl. 2005. "Keeping laying hens in furnished cages and an aviary housing system enhances their bone stability". British Poultry Science 46 (5):536-544. Newman, S., and S. Leeson. 1998. "Effect of housing birds in cages or an aviary system on bone characteristics". Poultry Science 77 (10):1492-1496. Niebur, Glen L, Michael J Feldstein, Jonathan C Yuen, Tony J Chen, and Tony M Keaveny. 2000. "High-resolution finite element models with tissue strength asymmetry accurately predict failure of trabecular bone". Journal of Biomechanics 33 (12):15751583. Rath, NC, GR Huff, WE Huff, and JM Balog. 2000. "Factors regulating bone maturity and strength in poultry". Poultry Science 79 (7):1024-1032. Reilly, Donald T, and Albert H Burstein. 1975. "The elastic and ultimate properties of compact bone tissue". Journal of Biomechanics 8 (6):393-405. Robling, Alexander G, Felicia M Hinant, David B Burr, and Charles H Turner. 2002. "Improved bone structure and strength after long‐term mechanical loading is greatest if loading is separated into short bouts". Journal of Bone and Mineral Research 17 (8):1545-1554. Saville, Paul D, and Michael P Whyte. 1969. "9 Muscle and Bone Hypertrophy: Positive Effect of Running Exercise in the Rat". Clinical Orthopaedics and Related Research 65:81-88. Tauson, R. 1998. "Health and production in improved cage designs". Poultry Science 77 (12):1820-1827. 29 Whitehead, CC. 2004. "Overview of bone biology in the egg-laying hen". Poultry Science 83 (2):193-199. Whitehead, CC, and RH Fleming. 2000. "Osteoporosis in cage layers". Poultry Science 79 (7):1033-1041. Wirtz, Dieter Christian, Norbert Schiffers, Thomas Pandorf, Klaus Radermacher, Dieter Weichert, and Raimund Forst. 2000. "Critical evaluation of known bone material properties to realize anisotropic FE-simulation of the proximal femur". Journal of Biomechanics 33 (10):1325-1330. Woo, SL, Steve C Kuei, David Amiel, MA Gomez, WC Hayes, FC White, and WH Akeson. 1981. "The effect of prolonged physical training on the properties of long bone: a study of Wolff's Law". The Journal of Bone and Joint Surgery American Volume 63 (5):780. 30 CHAPTER THREE DEVELOPMENT OF A MODEL FOR IN VIVO ASSESSMENT OF WHOLE BONE STRENGTH AND STIFFNESS INTRODUCTION The previous chapter examined the whole-bone mechanical and material properties of pullets from different housing systems. However, whole-bone mechanical testing requires the bird to be killed. To provide the opportunity for longitudinal studies and reduce the number of animals sacrificed, a process for assessing the mechanical and material properties of a bone in vivo was sought in the current study. Previous investigations have used computed tomography (CT) to generate 3-D finite element (FE) models of human femora to assess weaknesses in the bone that could lead to increased fracture risk (Keyak et al. 1997; Lotz, Cheal, and Hayes 1991). Simulated loads can be applied to such FE models that return the estimated stress distributions in the bone, predicting where and how the bone might fail (Keyak et al. 1997; Lotz, Cheal, and Hayes 1991). In order to generate an accurate FE model, appropriate material properties are essential. Bone is not a homogeneous material and such variations must be taken into account to accurately model its material properties (Currey 1988). Previous studies have used computer algorithms to assign element modulus values based on CT density using established relationships for bone (Van Rietbergen et al. 1995; Niebur et al. 2000). These relationships are developed using inverse analyses wherein the parameters of a 31 numerical model are adjusted to make it match an experimental result. Such relationships are commonly developed to account for the porosity of trabecular bone and the partial volume effects that come with medical imaging (Lotz, Gerhart, and Hayes 1990; Ciarelli et al. 1991). However, cortical bone is often assumed to be of uniform density and is assigned a single set of material properties (Rho, Hobatho, and Ashman 1995; Snyder and Schneider 1991). A goal of the current study was to develop a CT density-modulus relationship specific to cortical bone that would account for spatial distributions of density instead of assuming a single, homogeneous material property. A second method of predicting whole-bone strength has used quantitative CT (QCT) or dual energy X-ray absorptiometry (DXA) techniques. Clinical assessment using QCT or DXA makes key measures of bone geometry or density and correlates these measures with experimentally determined strength or stiffness (Cody et al. 1999; Ferretti et al. 1993; Martin, Severns, and Kabo 2004). QCT and DXA methods are less complex than FE modeling in that predictions are done statistically and do not require specialized solver software. However, FE models have been found to be more accurate than QCT or DXA in predicting failure loads (Cody et al. 1999; Crawford, Cann, and Keaveny 2003). The objective of the current study was to assess the efficacy of using FE and QCT models to predict the stiffness and failure moment of laying-hen bones. Using inverse analyses, a relationship between CT density and elastic modulus was to be determined for cortical bone by matching numerical data to the experimental mechanical tests performed in Chapter 2. It was hypothesized that there would be a single relationship between CT 32 density and elastic modulus that could be used in finite element models of both the humerus and the tibia of birds from various housing conditions. A second technique would implement QCT principles to develop models for bone stiffness and failure load based on measurements made directly from the CT scan. It was hypothesized that, since material properties have been shown to remain constant between housing environments, a simple measurement of bone geometry and average density from a CT scan would be sufficient to predict the stiffness and failure load of a bone from any housing environment. METHODS Tibiae and humeri of five white leghorn hens from each of three housing systems [conventional cage (CC), enriched cage (EN), and aviary (AV)] were collected at 78 weeks of age. The three environments provided varying opportunity for movement and exercise and, as Chapter 2 documents, may significantly alter cortical geometry. A uniform section of each bone was identified for testing and its boundaries were marked with drops of epoxy used previously in Chapter 2. These markers were used for consistent orientation of the bone during scanning and mechanical testing. Bones were then wrapped in saline-soaked gauze, sealed in a plastic bags, and scanned with a microCT scanner (eXplore Locus, GE Healthcare; Cleveland, OH). The long axis of each bone was aligned with the scanner bore to acquire axial slices and images were obtained with an in-plane resolution of 46 μm at a slice thickness of 46 μm. After scanning, the bones were potted and tested as described in Chapter 2. Raw CT scanner data were 33 reconstructed and converted into DICOM images using MicroView (Parallax Innovations; Ilderton, Ontario) to crop and align the tubular bone test section with the image axes. Finite Element Simulation for Inverse Analyses To generate the FE model, image slices were imported into a medical image processor (Mimics 15.01, Materialise; Leuven, Belgium). A thresholding mask was created to reproduce the cortical geometry by selecting all voxels with a CT density higher than 2000 Hounsfield units (Hu). This threshold level was selected so that the dimensions of the CT model matched the physical outer diameters taken with calipers at the ends and center of the test section. An FE mesh was automatically generated from this mask with each element composed of a 3-pixel square spanning 10 slices (a total of 90 voxels per element). This grouping produced elements measuring 138x138x460 μm with the long dimension lying along the longitudinal axis of the bone. The density of each element was calculated from the average of the 90 voxels contained within the element. Mesh geometry and element densities were loaded into Abaqus CAE (v6.11, Dassault Systemes; Vélizy-Villacoublay, France) for analysis. Boundary conditions were applied to duplicate the simply supported experimental setup (Figure 3.1). Nodes on each face were constrained such that both ends were free to pivot and one end was free to translate along the bone axis. For all specimens, the range of loading up to one-half the failure moment fell within the linear portion of the moment-rotation curve in each test (i.e. deformation was within the elastic range), so a linear, elastic FE model was appropriate for the study (Figure 3.2). 34 Based on previous studies of cortical and trabecular bone (Carter and Hayes 1976; Rho, Hobatho, and Ashman 1995), a relationship between CT density, 𝐶𝑇, and Young’s modulus, 𝐸, of each element was assumed to take the form 𝐸 = 𝐴 × 𝐶𝑇 B (3.1) where A and B were correlation coefficients that needed to be determined. Equation (3.1) was applied to each element, computing the element’s Young’s modulus based on its CT density. An inverse analysis technique was used to find the values of A and B that made the FE model bone-end rotation match that of each mechanical test. For this analysis procedure, an objective function, Z, of the form 30 2 𝑍 = � �𝜃𝜃model i − 𝜃𝜃exp � i i=1 (3.2) was defined where 𝜃𝜃model was the averaged rotation of the FE model’s ends and 𝜃𝜃exp was the experimental rotation at half the failure moment (Figure 3.2). Squares of the differences were summed over the 30 bones in this study. 35 𝑀𝑀 𝑀𝑀 𝜃𝜃 𝜃𝜃 Figure 3.1. Beam theory scenario used to model the current study showing the application of moments to the bone test section ends and the resultant deformations, measured by the rotation at the ends. The inverse analysis process began by evaluating the FE model and Equation (3.2) for an array of moduli computed using a range of A and B combinations. This process was automated using MATLAB (R2012a, Mathworks; Natick, MA) to iterate over a grid of A and B combinations (see Appendix A). For each set of correlation coefficients (A and B), a Python script (v3.3.2, Python Software Foundation; Beaverton, OR) was called to open the Abaqus model, redefine the modulus of each element from its density, solve for the static displacement, and read the rotation of each bone section end. In this way, a contour plot of the objective function against A and B was generated to assess the nature of the solution domain (i.e., determine if the solution was concave and if there was a global minimum). Since the FE model was linear, an expression was developed that calculated the response of the model for a given combination of A and B based on the known results 36 for another combination (see Appendix B). In this way, the plotting and optimization routines did not require that the FE simulation be run for each new modulus correlation. This technique produced the same results as evaluating the model at each point but significantly reduced computation time. Since the objective function appeared to be concave over the domain, an optimization routine was developed using the “fmincon” function in MATLAB employing an interiorpoint algorithm to find the correlation constants of Equation (3.1) that minimized the objective function in Equation (3.2). The “fmincon” routine is a gradient-based iterative solver and requires an initial guess for the parameters to be optimized. The solver operated by approximating the gradient of the objective function at the initial guess point using a finite difference method. Using this gradient, another point was picked “downhill” from the first, the function’s gradient was re-evaluated, and another iteration point was then selected. This process continued until the change in objective function (function tolerance) or input parameter values (step tolerance) from one point to the next fell below a specified tolerance, indicating a minimum of the objective function. To test the sensitivity, accuracy, and consistency of the solution, the optimization routine was run for an array of initial guesses. In addition, the effect of the objective function and iteration step tolerances were varied such that the algorithm consistently converged on a single solution in as few of iterations as possible. Once the optimal combination of A and B was found, the bending stiffness of each model was calculated 37 and compared to the bending stiffness computed from mechanical tests. Linear regression analysis was performed and considered significant for a p-value less than 0.05. Quantitative CT and Linear Regression Cropped and aligned micro-CT image slices were measured for QCT metrics using ImageJ (Wayne Rasband, NIH; Bethesda, MD) supplemented with the BoneJ plugin (Doube et al. 2010). For each slice, the cortex geometry was identified by a threshold that only included pixels with densities higher than 2000 Hu. Measurements of average CT density, cross-sectional area, major and minor principal second moments of area, and section moduli (second moment of area divided by the distance to the extreme fiber from the neutral axis) were taken for each slice. Geometrical measurements from each slice were then averaged over the length of the test section. To compute the volumetric CT density, ����, the average CT density from each slice, 𝐶𝑇i , was weighted with the area 𝐶𝑇 of that slice, 𝐴i , using the equation ���� = � 𝐶𝑇 i 𝐶𝑇i 𝐴i 𝑛𝐴̅ (3.3) where 𝐴̅ was the average cross-sectional area of all 𝑛 slices. Predictive values for initial stiffness and failure moment were adapted from the beam theory equations used in Chapter 2. Equation (2.4) was solved for the whole-bone stiffness, 𝐾, given by the equation 38 𝐾= 2𝐸𝐼 𝐿 (3.4) where E was the Young’s modulus, I was the second moment of area, and L was the test length (20mm for humeri, 30 mm for tibiae). Approximating the Young’s modulus of the bone by substituting ���� into Equation (3.1) and merging the constants into a single 𝐶𝑇 proportionality allowed Equation (3.4) to be expressed as ���� B 𝐼 𝐶𝑇 𝐾∝ 𝐿 (3.5) where the principal second moment of area, 𝐼, depended on the experimental configuration. The humerus was loaded with bending about the minor second moment of area and bending of the tibia was about the major second moment of area. To determine the failure moment, 𝑀𝑀f , Equation (2.5) was rewritten, yielding 𝐼𝜎 𝑀𝑀f = f 𝑏0 (3.6) where 𝜎f was the failure strength and 𝑏0 was the outer radius of the bone perpendicular to the neutral axis. Failure strength was assumed to be constant across all specimens. Replacing the strength with a proportionality and substituting the section modulus, 𝑆, for 𝐼 ⁄ 𝑏0 , Equation (3.6) was rewritten as 𝑀𝑀f ∝ 𝑆 (3.7) Again, the section modulus value (major or minor) depended on the orientation of the bone in the mechanical test. For the humerus, the minor axis value was used and, for 39 the tibia, the major axis value was used. The right hand sides of Equations (3.5) and (3.7), as measured from CT data, were correlated to the experimental whole-bone stiffness and failure moment, respectively, using linear regression analysis. Correlations were considered significant for a P-value less than 0.05. RESULTS Five bones (2 AV tibiae, one AV humerus, one CC humerus, and 1 EN tibia) were fractured during the harvesting and cleaning process, making them unsuitable for mechanical testing. The bones from 78 week-old specimens used in this study were of the same size and produced mechanical data similar to the 16 week-old pullets described previously. For the 25 bones that were tested in the current study, the point considered for the FE analysis fell in the linear portion of the moment-rotation curve as demonstrated in Figure 3.2. 40 8000 Bending Moment, M (N-mm) 7000 6000 5000 4000 3000 2000 1000 0 0 0.05 0.1 End Rotation, θ (rad) 0.15 0.2 Figure 3.2. Sample rotation-moment curve showing the selection of the moment used for FE analysis loading conditions (O) based on half the bending moment at failure (X). For all specimens, this point fell within the linear region of the moment-rotation curve. The objective function, Equation (3.2), was evaluated over a matrix of CT-modulus correlations from Equation (3.1) whose rows and columns corresponded to predefined values of A and B (Figure 3.3). This preliminary plot showed that the solution was indeed concave with a single combination of A and B that would minimize Equation (3.2). 41 a b -4 10 Figure 3.3. a) Topographical contour plot of the natural log of the objective function, Z, Equation (3.2), (out of the page) versus A and B. Lower “elevations” (dark blue) indicate that FE models using that correlation more closely matched experimental observations. b) Plot of the objective function, Z, versus A along the white line showing the concavity of the solution. The A and B combination that minimized the objective function (resulted in a model that most closely matched the experimental result) is identified by an X. 42 An optimization routine was written to automatically search for the minimum value of the objective function by iterating on A and B. In order to provide good convergence, a -09 function tolerance of 1.0x10 -10 and iteration step tolerance of 1.0x10 were required. Using these tolerances, the optimization solver converged on the same point and Equation (3.1) took the form 𝐸 = 37.97𝐶𝑇 0.7596 (3.8) regardless of the location of the initial guess. On average, the algorithm required 80 iterations to converge. Using this optimal combination of A and B, the stiffness of the FE model was computed and plotted against the experimentally determined stiffness for each tested bone (Figure 3.4). There was a strong correlation between the model and experimental 2 stiffness values (P<0.001, R =0.8723). The linear fit of the data showed a slope close to unity and a near-zero intercept (Figure 3.4). 43 90 FE Model Stiffness (N-m/rad) 80 y = 1.0643x + 0.237 R² = 0.8723 70 60 50 40 30 20 10 0 0 20 40 60 Experimental Stiffness (N-m/rad) 80 100 Figure 3.4. Correlation between the stiffness determined from mechanical tests and the stiffness of the FE bone models generated from CT using Equation (3.8). The same relationship was used for all bones in the study (both humerus and tibia). Following Equation (3.5), linear regression analysis of QCT-measured geometry and 2 mechanical properties showed a significant correlation (P<0.001, R =0.8814) (Figure 3.5a). Likewise, following Equation (3.7), there was a strong correlation between the 2 section modulus and failure moment (P<0.001, R =0.9539) (Figure 3.5b). 44 a 90 80 y = 0.0617x + 4.078 Stiffness (N-m/rad) 70 R² = 0.8814 60 50 40 30 20 10 0 0 b 400 600 CTB * I / L 800 1000 1200 8 7 Failure Moment (N-m) 200 y = 0.3549x - 0.2665 R² = 0.9539 6 5 4 3 2 1 0 0 5 10 15 Section Modulus, S, (mm3) 20 25 Figure 3.5. a) Linear fit of bending stiffness versus the metric generated from QCT and beam theory following Equation (3.5). b) Linear fit of failure moment versus the section modulus measured from QCT following Equation (3.7). 45 DISCUSSION The current study examined two techniques for estimating the mechanical properties of whole laying-hen bones using micro-CT data. Through inverse analysis, a single set of constants was found to relate CT density to Young’s modulus for CT-generated finite element models as reported by Equation (3.8). By design, this optimized relationship resulted in a strong correlation between the FE model and experimental stiffnesses. Part of the success of the FE model was likely due to the high resolution micro-CT geometrical data that was required to accurately measure the small features and thin cortices of the bones (Prevrhal, Engelke, and Kalender 1999). Such high-resolution images, in turn, allowed development of a very fine mesh that was able to accurately represent the irregular geometry of the cortex. A second method of estimating mechanical properties from the micro-CT data involved calculating key geometrical measurements of the bones as based on beam theory. The study indicated that whole-bone bending stiffness strongly correlated with the second moment of area, and the failure moment was strongly correlated with the section modulus of the cortex. These findings confirm the conclusion of the previous chapter that the strength and stiffness of laying-hen bones were determined largely by bone cross-sectional geometry, versus their material properties. On average, QCT measurements of the second moment of area predicted the experimental stiffness with the same degree of accuracy as the FE model. These findings contrast with previous studies that reported improved predictions using FE models over 46 QCT (Crawford, Cann, and Keaveny 2003; Cody et al. 1999). A possible explanation for this inconsistency is that the QCT methods used the full resolution of the CT scan, whereas the FE model grouped multiple voxels together to form each element. This grouping was required to keep the number of elements and size of the FE models to reasonable levels but effectively reduced the resolution of the model. It is possible that developing a model with an element for every voxel would improve the results but would have exceeded the capabilities of the machines used to run the current analyses. Another reason there was no improvement in the prediction of stiffness using FE analysis could be a function of the type of bone considered in the study. The previous work compared specimens from a vertebral body (Crawford, Cann, and Keaveny 2003) and the proximal femur (Cody et al. 1999), structures that derive much of their strength from trabecular bone. The spatial distribution of density in trabecular bone is an essential component of its structure and is best accounted for through the use of CT; a single measure of geometry or average density cannot provide such information. However, the long bone diaphyses studied here were comprised entirely of cortical bone. Thus, the similar performance of QCT and FE techniques in the current study suggests that the mechanical properties of cortical bone can be characterized solely by average geometry and CT density without consideration to the distribution of density. Applying the density-modulus relation developed using the FE model to QCT predictions for stiffness allowed both bone types of different moduli to be analyzed together and improved the strength of the correlation. Through further study, it may be possible to develop a relationship correlating bone strength to CT density similar to the modulus47 density correlation used in the current study. In this way, the average density of each bone could be accounted for when developing a predictive metric for failure moment from QCT data and improve upon the correlation in the current study, which was made on geometry alone. A major limitation of the current study is the simplicity of the finite element model. Although the model accurately reproduced the geometry and density distribution, the analysis was limited to a single static displacement. Cortical bone is not a simple linear elastic solid; its cellular structure makes it a viscoelastic material whose behavior is dependent on strain rate (Carter and Hayes 1977). Since all tests in the current study were performed at the same speed, however, it has been assumed that such viscoelastic effects between bones and environments were negligible. However, it is likely that for loading rates significantly different than that used in the current study, a potential change in Young’s modulus would alter the stiffness of the bones and require a new set of correlation coefficients to account for such changes. There are obvious advantages and limitations of both techniques presented here. While more computationally intensive, FE models can be developed to predict the behavior of bone structures under a wide range of loading and boundary conditions to simulate physiological scenarios. Perhaps the most beneficial aspect of FE models is that, once their behavior has been tested and validated, subsequent models can be created and analyzed without additional validations. On the other hand, QCT models correlating a specific mechanical characteristic with geometrical and densiometric aspects require a 48 full dataset to calibrate the response. However, such QCT models are an attractive substitute for the time and labor-intensive testing methods used in Chapter 2, where the test conditions between environments and bones were the same. Although the current study used micro-CT data on ex vivo samples, recent advances in the medical imaging field have introduced in vivo micro-CT with resolutions matching the current study. In future studies, it may be possible to “virtually” assess a live hen’s bone stiffness and breaking strength using micro-CT generated metrics. Also, since the animal wouldn’t need to be sacrificed to perform a mechanical test, the same animal could be used in longitudinal studies to track bone growth and development over time. In conclusion, the current study examined two techniques for assessing whole bone mechanical properties using micro-CT. It was found that FE and QCT techniques predicted the bone bending stiffness to a similar degree. In addition, QCT data were able to predict the failure moment from geometrical properties alone. These findings have significant implications in poultry science where ex vivo mechanical testing is the standard practice for determining whole bone properties. Following the methods of the current study, a simple measurement of geometry from an in vivo micro-CT could provide researchers with an accurate assessment of bone strength and stiffness. A suggestion for future work would be to develop a non-linear FE model capable of predicting bone behavior up to failure and add a CT density metric to the QCT correlation for failure moment. An additional area of study should assess the feasibility of using in vivo micro-CT to predict bone mechanics by performing mechanical tests on 49 bones that were imaged before sacrificing the animal. In this way, future studies could make accurate predictions of laying-hen bone strength and stiffness from an in vivo micro-CT scan. 50 APPENDIX B 51 LINEARIZATION METHOD Equation (3.1), which was originally used to find the modulus of a single element based on its CT density, can be applied to the model to express the stiffness of a whole bone. In this case, the bulk modulus, 𝐸0 , of the bone can be written as 𝐸0 = 𝐴𝐶𝑇0 B (3.9) where, again, 𝐴 and 𝐵 are constants and 𝐶𝑇0 is defined as the bulk CT density, a unique property of each bone that must be determined using the bone’s finite element model. In order to calculate the bulk CT density, a simplifying assumption is made that restricts the model to small deformations. In this case, the deformation is the rotation of each end of the test section,𝜃𝜃, and can be related to the bulk modulus using the proportionality 𝐸0 ∝ 1 𝜃𝜃 (3.10) Substituting Equation (B.9) for the modulus gives 𝐴𝐶𝑇0 B ∝ 1 𝜃𝜃 (3.11) This relation was checked by solving for the end rotation using a single 𝐵 and two values of 𝐴. It was confirmed that changing the bulk modulus had an inverse effect on the rotation (eg. doubling the value of 𝐴 reduced 𝜃𝜃 by half). 52 Thus, Equation (B.11) could be used to solve for 𝐶𝑇0 of a particular model by using two different values of 𝐵, namely (3.12) 1 𝐴𝐶𝑇0 𝐵2 ∝ 𝜃𝜃2 (3.13) 1 𝜃𝜃2 𝐵1−𝐵2 𝐶𝑇0 = � � 𝜃𝜃1 and 1 𝐴𝐶𝑇0 𝐵1 ∝ 𝜃𝜃1 (3.14) Dividing Equation (B.12) by (B.13) and simplifying gives the expression In this way, 𝐶𝑇0 can bet determined from two executions of the FE model. Once 𝐶𝑇0 has been calculated for a bone model, solving for the end rotation of a new combination of 𝐴 and 𝐵 no longer required the full FE model. Instead, the end rotation, 𝜃𝜃i , for the new combination, 𝐴i and 𝐵i , could be found from the FE model result, 𝜃𝜃1 , for a known combination, 𝐴1 and 𝐵1 , using 𝐴 𝜃𝜃i = 𝜃𝜃1 1 𝐶𝑇 �B1−Bi � 𝐴i 53 (3.15) CODE %% setup.m % Matlab script to call a python script to run abaqus models with % predetermined parameters as to calculate the bulk CT density of each % specimen. %% Define names of specimen to be modeled, read data from .xls file load xldata % cell containing specimen numbers and moments/rotations at % half-failure (generated using xlsread function) %% Compute the bulk CT density for each specimen A0 = 2; % set single value for A B0 = [0.5,1.5]; % set initial values for B for i = 1:size(xldata,1) % iterate over all specimens in xldata specnum = xldata{i,1}; % define the specimen ID number for j = 1:2 fileID = fopen('inputs.dat','w'); fprintf(fileID,'%g\r\n%g\r\n%s',A0,B0(j),specnum); fclose('all'); ! abaqus cae nogui=runmodel.py theta0(i,j) = str2double(fileread('rotation.dat')); delete('abaqus.rpy*', '*.dat','Job-1.*','*.jnl') end CT(i) = (theta0(i,1)/theta0(i,2))^(1/(B0(2)-B0(1))); % compute the % bulk CT for each specimen end save allModelsSetup A0 B0 CT theta0 xldata 54 # runmodel.py #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Python script to call a single abaqus .cae file in the same folder # and change elastic modulus values according to density and # parameters specified in inputs.dat. After running the job, find the # maximum and minimum rotations (both will occur at the ends) and then # average their magnitudes. Finally, write this averaged value to # rotation.dat to be accessed by MATLAB. #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Author: # Trevor DeLand # Revisions: # 2013-Feb-27...Created script # 2013-Mar-08...Merged with readdata.py for streamlining # 2013-Apr-02...Match model file to specimen listed in inputs.dat; # error exits; removed specific naming requirements # for model, material, and job subroutines # 2013-Jun-13...Made changes for reading full specimen name instead # of just specimen number from input file #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Load modules for added functionality with Abaqus and Command Prompt from abaqus import * from abaqusConstants import * import material from odbAccess import * import os import job #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Read a and b values and specimen number from inputs.dat (written by objfun.m) f = open('inputs.dat','r') inputs = [] for line in f: # generate inputs vector from lines of inputs.dat inputs.append(line) f.close() # close inputs.dat a = float(inputs[0]) b = float(inputs[1]) specnum = str(inputs[2]) # formats specimen number as string #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Find the current directory and the name of the matching abaqus model within folder = os.getcwd() filename = [] for files in os.listdir('.'): # list of files in current directory if files.endswith('.cae') and files.find(specnum) >= 0: # if find doesn't # match specnum in files, returns -1 filename = files if filename == []: # exit if no matching file is found errmsg = 'Error: No ABAQUS model found for specimen %s' % (specnum) os.system('echo. && echo ' + errmsg) os._exit(0) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Open the .cae model database matching excel file in the current directory filepath = folder + '\\' + filename openMdb(filepath) 55 modellist = mdb.models.keys() modelname = [] for model in modellist: if model.find(specnum) >= 0: modelname = model if modelname == []: errmsg = 'Error: Cannot locate the matching model in database %s' \ % (filename) os.system('echo. && echo ' + errmsg) os._exit(0) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Redefine the Young's modulus using a, b, and the density model = mdb.models[modelname] matlist = model.materials # repository: can only be indexed by material name matnamelist = matlist.keys() # generate list of material names for indexing for matname in matnamelist: # set to iterate over all materials material = matlist[matname] density = material.density.table[0][0] # read material's average CT density modulus = a*density**b # calculate new modulus from density using a & b properties = (modulus, 0.3) # define new modulus and poisson ratio table material.Elastic(table=(properties,)) # set new material properties #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Find and run the job jobname = mdb.jobs.keys()[0] # retreive the first job name for repository index job = mdb.jobs[jobname] job.submit(consistencyChecking=OFF) job.waitForCompletion() #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Read the results from the generated output database parameter = 'UR' # declare which field value to analyze. eg UR = rotation j = 0 # data entry to examine. eg 0 = UR1, 1 = UR2, 2 = UR3 # Access the output database from the job odbPath = jobname + '.odb' odb = openOdb(path=odbPath) lastStep = odb.steps.keys()[-1] # find name of last step for repository index results = odb.steps[lastStep].frames[-1].fieldOutputs[parameter].values maxU = 0 # initialize min and max values minU = 0 # Search for the maximum and minimum values (max ur > 0 and min ur < 0) if node.data[j] > maxU: maxU = node.data[j] if node.data[j] < minU: minU = node.data[j] urmodel = (maxU - minU)/2 # calculate the average rotation #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Write the averaged rotation data to rotation.dat for MATLAB access outputFile = open('rotation.dat','w+') outputFile.write(str(urmodel)) # write argument needs to be a string outputFile.close() 56 %% ExpandPlot.m % After running setup.m, extrapolate the results of the Abaqus model over a % range of values and plot the results %% Fill out the results over the grid load allModelsSetup difference = 0; for i = 1:size(xldata,1) %% Set values for A and B abRes = 2000; a = linspace(25,50,abRes); b = linspace(0.725,0.8123,abRes); thetamodel = zeros(abRes); %% Run the objective function over the specified A and B values for j = 1:abRes for k = 1:abRes thetamodel(j,k) = theta0(i,1)*A0/a(j)*CT(i)^(B0(1)-b(k)); end end %% Calculate the summed squares of the difference difference = difference + (thetamodel - xldata{i,3}).^2; end %% Locate the absolute minimum [colmins, ibmins] = min(difference'); [absmin, iamin] = min(colmins); amin = a(iamin); bmin = b(ibmins(iamin)); fprintf('Global Minimum:\nA = %g \nB = %g \n\n',amin,bmin); figure plot(a,colmins); %% Plot the results %surf(a,b,log(difference'),'LineStyle','none'); axis tight %xlabel('A'); ylabel('B'); figure contour(a,b,log(difference'),15); axis tight hold on xlabel('A'); ylabel('B'); plot(amin,bmin,'ro') 57 %% Optimization_fmincon.m % Routine to read relevant info from experimental data (.xls file) and % perform an inverse analysis optimization using fmincon %% Initialize clear clc load allModelsSetup t1 = datestr(now); data = {A0 B0 theta0 CT xldata}; % cell of constants % Create anonymous fxn for passing data into objfun without iterating f = @(x)objfun(x,data); % Set the desired options for the optimization options = optimset('fmincon'); options.TolFun = 1e-06; options.TolX = 1e-06; options.FinDiffType = 'forward'; options.PlotFcns = @optimplotfval; options.Algorithm = 'interior-point'; options.Display = 'iter-detailed'; %% Run the optimization x0 = [5 1]; % Optimization starting point Y = [0.0476 0.67]; % Linear constraints: x(dot)Y <= z z = 1; [xopt,mindiff,exitflag,output] = fmincon(f,x0,[],[],[],[],[0,0],... [50,1.5],[],options); t2 = datestr(now); fprintf(' A = %g \n B = %g \n Difference = %g \n',xopt(1),xopt(2),... mindiff); % display the optimization results fprintf('\n Started: %s \n Finished: %s \n',t1,t2); % display timer 58 Table 3.1. Raw data from Chapter 3. Specimen CT Density (HU) AV1130HL AV1130TL AV1138HL AV1138TL AV1140HL AV1172HL AV1172TL CC1014HL CC1014TL CC1035HL CC1035TL CC1071TL CC1076HL CC1076TL CC1087HL CC1087TL EN1251HL EN1251TL EN1278HL EN1278TL EN1302HL EN1302TL 2812 3192 2816 3189 2918 2940 3306 3001 3203 2888 3331 3303 2995 3312 2906 3171 2933 3206 2994 3398 2910 3207 QCT Measurements Second Moment of Area Section Modulus 4 3 (mm ) (mm ) 45.1 14.1 73.5 20.0 38.0 12.4 56.8 16.4 30.2 10.2 41.9 13.4 71.6 19.6 28.9 10.2 45.9 13.3 21.0 8.00 44.1 14.0 47.0 14.2 30.5 10.4 68.3 18.1 32.2 11.1 58.2 17.0 28.4 9.92 50.6 14.6 22.4 8.18 43.1 13.7 34.1 11.7 58.6 16.5 59 CT Model Stiffness B (CT *I/L) 941 1126 795 870 647 905 1127 633 706 447 698 739 668 1077 690 887 612 777 490 693 730 901 Experimental Results Stiffness (N-m/rad) Stiffness (N-m/rad) Failure Moment (N-m) 68.6 84.7 58.8 62.7 45.1 62.9 77.0 45.0 47.9 27.8 47.2 53.7 46.3 74.8 50.2 63.4 49.1 55.8 33.0 48.9 55.4 64.6 63.7 76.3 51.3 61.8 46.0 48.5 69.4 37.1 48.9 30.3 50.2 51.6 37.6 68.4 45.9 61.1 34.9 59.1 36.5 50.0 49.8 61.8 4.54 7.37 4.25 5.29 3.52 4.07 6.25 3.37 4.07 2.83 4.91 5.05 3.57 6.27 3.68 5.39 3.22 4.88 2.96 4.60 3.69 5.86 Table 3.1 Cont’d EN1321HL EN1340HL EN1340TL 2976 2857 3270 21.4 31.5 64.0 7.90 10.7 17.6 467 665 998 60 33.6 49.0 75.6 35.5 47.6 69.1 2.37 3.19 6.39 REFERENCES 61 REFERENCES Carter, D. R., and W. C. Hayes. 1976. "Bone compressive strength - influence of density and strain rate". Science 194 (4270):1174-1176. Carter, Dennis R, and Wilson C Hayes. 1977. "The compressive behavior of bone as a two-phase porous structure". The Journal of Bone and Joint Surgery. American Volume 59 (7):954-962. Ciarelli, MJ, SA Goldstein, JL Kuhn, DD Cody, and MB Brown. 1991. "Evaluation of orthogonal mechanical properties and density of human trabecular bone from the major metaphyseal regions with materials testing and computed tomography". Journal of Orthopaedic Research 9 (5):674-682. Cody, Dianna D, Gary J Gross, Fu J Hou, Horace J Spencer, Steven A Goldstein, and David P Fyhrie. 1999. "Femoral strength is better predicted by finite element models than QCT and DXA". Journal of Biomechanics 32 (10):1013-1020. Crawford, R Paul, Christopher E Cann, and Tony M Keaveny. 2003. "Finite element models predict in vitro vertebral body compressive strength better than quantitative computed tomography". Bone 33 (4):744-750. Currey, J. D. 1988. "The effect of porosity and mineral content on the Young's modulus of elasticity of compact bone". Journal of Biomechanics 21 (2):131-139. Doube, Michael, Michał M Kłosowski, Ignacio Arganda-Carreras, Fabrice P Cordelières, Robert P Dougherty, Jonathan S Jackson, Benjamin Schmid, John R Hutchinson, and Sandra J Shefelbine. 2010. "BoneJ: Free and extensible bone image analysis in ImageJ". Bone 47 (6):1076-1079. Ferretti, José Luis, Ricardo Francisco Capozza, Nélida Mondelo, and José Rubén Zanchetta. 1993. "Interrelationships between densitometric, geometric, and mechanical properties of rat femora: inferences concerning mechanical regulation of bone modeling". Journal of Bone and Mineral Research 8 (11):1389-1396. Keyak, Joyce H, Stephen A Rossi, Kimberly A Jones, and Harry B Skinner. 1997. "Prediction of femoral fracture load using automated finite element modeling". Journal of Biomechanics 31 (2):125-133. Lotz, JC, EJ Cheal, and WC Hayes. 1991. "Fracture prediction for the proximal femur using finite element models: Part I--Linear analysis". Journal of Biomechanical Engineering 113 (4):353. 62 Lotz, Jeffrey C, Tobin N Gerhart, and Wilson C Hayes. 1990. "Mechanical properties of trabecular bone from the proximal femur: a quantitative CT study". Journal of Computer Assisted Tomography 14 (1):107-114. Martin, Daniel E, Anne E Severns, and J Michael Kabo. 2004. "Determination of mechanical stiffness of bone by pQCT measurements: correlation with non-destructive mechanical fourpoint bending test data". Journal of biomechanics 37 (8):1289-1293. Niebur, Glen L, Michael J Feldstein, Jonathan C Yuen, Tony J Chen, and Tony M Keaveny. 2000. "High-resolution finite element models with tissue strength asymmetry accurately predict failure of trabecular bone". Journal of Biomechanics 33 (12):1575-1583. Prevrhal, Sven, Klaus Engelke, and Willi A Kalender. 1999. "Accuracy limits for the determination of cortical width and density: the influence of object size and CT imaging parameters". Physics in medicine and biology 44 (3):751. Rho, JY, MC Hobatho, and RB Ashman. 1995. "Relations of mechanical properties to density and CT numbers in human bone". Medical engineering & physics 17 (5):347-355. Snyder, Susan M, and Erich Schneider. 1991. "Estimation of mechanical properties of cortical bone by computed tomography". Journal of Orthopaedic Research 9 (3):422-431. Van Rietbergen, B, Harrie Weinans, R Huiskes, and Anders Odgaard. 1995. "A new method to determine trabecular bone elastic properties and loading using micromechanical finiteelement models". Journal of Biomechanics 28 (1):69-81. 63 CHAPTER FOUR THE EFFECTS OF FALL HEIGHT AND INTERFACE COMPLIANCE ON FRACTURE PATTERNS IN THE DEVELOPING PORCINE SKULL INTRODUCTION Trauma is the leading cause of hospitalization for children worldwide (Lam and Mackersie 1999). Three quarters of these cases and 70% of injury-related deaths come as a result of head trauma (Lam and Mackersie 1999). Falls are the leading cause of childhood head trauma and the third-leading cause of death in children aged 1-4 years (Reece and Sege 2000; Greenes and Schutzman 1999; Hobbs 1984; Hall et al. 1989). Head injuries are also very prevalent in cases of abuse, accounting for 80% of mortalities in battered children (Case et al. 2001). An estimated 30% of pediatric head trauma cases result in skull fracture with the majority of fractures occurring in the parietal bone (Schutzman et al. 2001; Hobbs 1984; Billmire and Myers 1985). Forensic investigations often examine characteristics of these fractures in an attempt to distinguish accidental from inflicted injury. However, reaching a verdict based on skull fracture alone is often questionable as both cases may produce similar fractures (Billmire and Myers 1985). The purpose of the current study was to gain a deeper understanding of the fracture mechanics of the skull through controlled simulation of traumatic head injuries. 64 Early studies have shown that simulated falls onto the occiput using infant cadaver surrogates generate multiple fractures that often cross sutures (Weber 1984, 1985). However, due to availability and ethical concerns, such experiments are rarely feasible today. Studies using anthropomorphic dummies have documented accelerations and loading rates of fracture-inducing impact scenarios but provide no insight into the type of fractures generated under such conditions (Coats and Margulies 2008; Prange et al. 2003). Recent studies by our laboratory have demonstrated the utility of an infant porcine model in simulating impacts to the developing skull and assessing the resultant fracture patterns. An initial study showed that, in terms of mechanical development of the skull, a pig’s age in days closely correlates with months in a human (Baumer et al. 2009; Baumer et al. 2010). Furthermore, it has been documented that increasing the energy of impact to an entrapped head increases the number and length of cranial fractures (Baumer et al. 2010; Powell et al. 2012). At levels equivalent to high-energy entrapped conditions, controlled free-fall head drop impacts result in fractures resembling low-energy entrapped head impacts (Powell et al. 2013). As a result of these porcine model tests, it has been concluded that, in addition to impact energy, the method used to constrain the head during impact has a profound effect on the degree of fracture (Powell et al. 2013; Baumer et al. 2010; Powell et al. 2012). Common to all the infant porcine studies, fractures in the impacted bone appear to originate at the suture and propagate radially inward toward the centrally-located site of impact (Baumer et al. 2010; Powell et al. 2012; Powell et al. 2013). This phenomenon is thought to occur due to “out-bending” of the bone; a theory based on fracture 65 patterns observed in the adult human skull (Gurdjian, Lissner, and Webster 1947; Gurdjian, Webster, and Lissner 1950; Gurdjian, Webstef, and Lissner 1953). This outbending effect was explored in a study by Powell that used finite element analysis to approximate the skull as a simple sphere and compared the stresses developed during impacts under different constraint conditions (Powell 2010). Powell’s finite element model, however, was limited in that he treated the skull as a solid, homogeneous body. In an infant, the sutures that connect bones of the skull are different than bone and may be more cartilaginous depending on the degree of maturation (Baumer et al. 2009; Coats and Margulies 2003; Coats and Margulies 2006). An additional goal of the current study was to develop a model that might account for discontinuity at the sutures by simulating a cranial impact to a single bone of the skull and applying boundary conditions at its periphery that might model constraint offered to the bone by compliant sutures. Fracture-inducing falls in children are most often from heights of 6 feet or less and the resultant fractures are linear and uncomplicated (Reece and Sege 2000). Higher falls can cause complex fractures that could easily be misdiagnosed as abuse (Reece and Sege 2000; Hobbs 1984; Greenes and Schutzman 1999). In addition, a more compliant interface has been shown to dramatically reduce the incidence of skull fracture for a given impact energy (Weber 1984, 1985).The impetus for the current study was to examine the individual effects of fall height and interface compliance on skull fractures using the infant porcine model. 66 The hypotheses of this study were threefold. First, raising the drop height of controlled head drops would increase the frequency and length of fractures on the infant porcine skull. Second, high-energy drops would induce fractures in bones away from that which was impacted. Third, impacts to a compliant interface would have the same effect as lowering the drop height, that is to reduce the frequency and length of cranial fractures for a given level of impact energy. In addition, a simple finite element model was developed to illustrate a mechanism for the formation of fractures at the edge of the impacted bone adjacent to the sutures. METHODS A total of 87 pigs aged 2-19 days old were collected by a local hog farm and frozen within 12 hours of death at -20°C. All animals died of natural causes. While frozen, heads were separated via cervical dislocation using an upright bandsaw and stored at 20°C. Specimens were thawed at room temperature for 24 hours prior to testing. The heads were screened for pre-existing head trauma. To generate a single impact, a custom drop tower was built. A trolley was free to travel along the length of a 2.5-meter vertical rail on linear bearings (Figure 4.1). An electromagnetic solenoid clamp held the trolley in place and allowed easy adjustment of the drop height. If the required drop height exceeded the height of the tower, a spring at the top of the tower could be used to increase the velocity of the dropped head. 67 Figure 4.1. Experimental setup of the drop tower and porcine head specimen A second set of bearings connected a small floating rod to the trolley which could be locked together via a second solenoid clamp. Specimens were strapped to a padded aluminum plate that attached to the end of the floating rod via an adjustable ball-andsocket clamp (Figure 4.2). Specimens were oriented so as to impact the center of the parietal bone. In its entirety, this apparatus allowed specimens to be dropped from a wide range of heights and impacted at a consistent orientation while minimizing the mass attached directly to the head. 68 Solenoid Clamp Drop Trolley Floating Rod Adjustable Positioning Clamp Padded Mounting Plate Velcro Straps Figure 4.2. Photo detail of the floating rod and specimen mounting apparatus attached to the drop trolley Impact energy for these experiments was set to be double that of the previous study of Powell et al. (Powell et al. 2013). Specimens 2-9 days old were impacted at 13 J and specimens 10-18 days old at 23 J by adjusting the drop height determined by using the following equation: ℎ= 𝑈 𝑚𝑔 (4.1) where 𝑈 was the specified drop energy, 𝑚 was the mass of the head, and 𝑔 was the gravitational acceleration. If the value from Equation (1) was less than the drop tower height, the trolley was positioned and held such that the specimen was at a height h 69 above the impact surface. If the height of the tower was not sufficient, the spring at the top supplemented the gravitational potential energy. With the free end of the spring pushing against the trolley, the position of the fixed end could be adjusted to result in the required spring compression determined by 2𝑀𝑔�ℎ − ℎ0 � 𝑥= � 𝐾 (4.2) where 𝑀 was the combined mass of the trolley and specimen, ℎ was the required drop height determined using Equation (4.1), ℎ0 was the maximum height of the drop tower, and 𝐾 was the spring stiffness. An electronic controller initiated the drop by simultaneously releasing both clamps (Figure 4.3a). Upon impact, the trolley struck a padded surface and quickly came to rest. With the trolley-mounted clamp open, the floating rod continued downward and the head struck the impact surface (Figure 4.3b). The rigid surface was a 1.5 cm thick aluminum disk 20 cm in diameter (N=32). A more compliant interface was provided by covering the disk with thin, loop pile commercial-grade carpet, referred to as carpet 1 (N=23). A second level of interface compliance was provided by covering the disk with carpet cushion (6-lb, 7/16-inch virgin foam rebond, Carpenter) and a cut pile polyester carpet (1/2-inch pile height, 29.1 oz/sq yard, Frontier Park) referred to as carpet 2 (N=32). The compliance of each carpet was tested using a Clegg Impact Tester (Model 95049A, Lafayette Instrument Co.; Lafayette, IN) which measured the peak deceleration of a 2.25 kg missile dropped from a height of 45 cm on to the surface. A 2500-lb load transducer (model 1010AF-2.5K, Interface; Scottsdale, AZ) mounted immediately behind 70 the aluminum disk recorded the reaction force of the impact at 10,000 Hz (Figure 4.1). To guarantee a single strike, a voltage comparator circuit monitored the transducer output and reactivated the trolley-mounted clamp, catching the floating rod when the impact reaction force returned to zero to prevent a second impact. Figure 4.3. a) The specimen and drop trolley were suspended above the impact interface and held in place by solenoid clamps. b) When released, the trolley and specimen mounting rod fell together until the trolley struck the padded stop. The floating rod and specimen continued downward, striking the instrumented impact interface. After impact, the scalp and other soft tissues were carefully removed from the skull. Any diastatic fractures were photographed and marked on a standard diagram (Figure 4.4). All remaining soft tissue was removed by boiling the skulls in water and detergent for four hours. After air-drying, the bones were re-assembled and glued together. Bone 71 fractures were added to the diagram and length of each fracture was recorded to the nearest millimeter using a flexible measuring tape. Figure 4.4. Right-side and occipital-oriented views used for marking fracture patterns. The impact site is indicated by an X Fracture diagrams were digitized and entered into a Geographic Information System (GIS), which overlaid multiple patterns from a group of specimens for simple means of assessing the overall patterns. The current study included a total of 3 groups based on the impact interface (rigid, carpet 1, carpet 2). A previous study by Powell et al. (Powell et al. 2013) used the same methods but with half the drop energy onto a rigid interface. These data were revisited here for a fourth group. Two factor linear regression analyses were used to examine differences in the responses of peak impact force, impact loading duration, and total bone fracture length with respect to specimen age and impact energy using the lower-energy data of Powell et al. (2013) and the high energy data of the current study. A second set of two-factor linear regression analyses examined the responses of peak impact force, impact loading duration, and total bone fracture length to specimen age and impact interface for high- 72 energy impacts. Statistical results were considered significant for a p-value less than 0.05. To study impact stress distributions in the parietal bone theoretically, the right parietal bone was modeled as a circular domed shell in Abaqus CAE (v6.11, Dassault Systemes; Vélizy-Villacoublay, France). Based on measurements from dissected skulls, the model dome measured 60 mm across and 10 mm high with a radius of curvature of 120 mm. Isotropic elastic material properties were assigned with a shell thickness of 1.5 mm. An impact was simulated by applying a static pressure to the center of the outer surface of 2 the shell with the area of contact set at 200 mm based on unpublished pressure film data from a previous study (Powell et al. 2013). The circular edge of the shell was constrained against the impact force but allowed circumferential expansion to model the potential effect of the lack of mechanical constraint afforded by sutures on the periphery of the parietal bone (Figure 4.5). 73 Figure 4.5. Domed shell used to model impact loading on the right parietal. A pressure was applied to the center of the outer surface with constraints modeling the rim resting on a frictionless table. RESULTS Energy Effects Peak forces generated by the current set of high-energy impacts showed a significant correlation with age (P<0.001) (Figure 4.6). Additionally, these forces were significantly larger than those generated by Powell et al. (2013) referenced here as low-energy (P<0.001) (Figure 4.6). 74 Low-Energy 2000 High-Energy Peak Impact Force (N) 1800 1600 1400 1200 1000 800 600 400 200 0 0 5 10 Age (Days) 15 20 Figure 4.6. Peak impact force versus age for the current (high-energy) and previous (low-energy) studies Time to peak load was measured up to the point at which the impact force began to return to zero. There was no correlation between time to peak load and specimen age (P=0.698). Time to peak load for the current high-energy impacts was significantly longer than the low-energy impacts of Powell et al. (2013) (P<0.001) (Figure 4.7). 75 Low Energy 8.0 High-Energy Time to Peak Load (ms) 7.0 6.0 5.0 4.0 3.0 2.0 1.0 0.0 0 5 10 15 20 Age (Days) Figure 4.7. Time to peak load versus age for high and low energy impacts There was no correlation between age and total fracture length (P=0.996) (Figure 4.8). Total fracture length was significantly greater in the current high-energy impacts than in the previous low-energy work of Powell et al. (2013) (P<0.001). Diastatic fractures were present in 42% of low-energy and 16% of high-energy impacts (Appendix C). 76 Low-Energy Total Fracture Length (mm) 300 High-Energy 250 200 150 100 50 0 0 5 10 15 20 Age (Days) Figure 4.8. Total bone and diastatic fracture length versus age for high and low energy impacts to a rigid surface Interface Effects Results of the Clegg impact test showed that carpet 1 generated a peak deceleration of 413 G and carpet 2 generated a peak deceleration of 152 G. As was the case when impacting a rigid surface, peak impact forces generated by the two carpeted interfaces showed a significant increase with specimen age (P<0.001). The most compliant carpet and pad surface (carpet 2) generated significantly higher impact forces for every age than the commercial carpet (carpet 1) and rigid interfaces (P<0.001). There were no significant differences between the forces generated by carpet 1 and the rigid interface. (Figure 4.9). 77 Rigid 3000 Carpet 1 Carpet 2 Peak Impact Force (N) 2500 2000 1500 1000 500 0 0 5 10 Age (Days) 15 20 Figure 4.9. Peak impact force versus age for high energy impacts to three different interfaces. Carpet 1 = thin loop pile commercial carpet. Carpet 2 = cut pile carpet with pad underlay There was no significant correlation between time to peak load and age (P=0.692) or impact interface (P=0.631) (Figure 4.10). 78 9 Rigid Carpet 1 Carpet 2 Time to Peak Load (ms) 8 7 6 5 4 3 2 1 0 0 5 10 Age (Days) 15 20 Figure 4.10. Time to peak load versus age for high-energy impacts to three different interfaces. Carpet 1 = thin loop pile commercial carpet. Carpet 2 = cut pile carpet with pad underlay Across the 3 interfaces there was no significant correlation between total fracture length and specimen age (P=0.231). Impacts to the rigid surface produced the most fracture with an average total fracture length of 118±73 mm. Commercial carpet (carpet 1) produced the second most fracture (54±43 mm) and had one unfractured specimen. Carpet with padding (carpet 2) produced the least total fracture length (20.3±25.4 mm) and had 7 non-fractured specimens. Total fracture length for impacts to the rigid interface was significantly higher than both carpet 1 (P=0.002) and carpet 2 (P<0.001). Carpet 1 produced a significantly higher total fracture length than carpet 2 (P=0.009) (Figure 4.11). Diastatic fracture frequencies were 16% for both rigid and carpet 1 impacts and 13% for carpet 2 impacts (Appendix C). 79 Total Fracture Length (mm) 350 Rigid Carpet 1 Carpet 2 300 250 200 150 100 50 0 0 5 10 Age (Days) 15 20 Figure 4.11. Total bone and diastatic fracture length versus age for high energy impacts to three different interfaces. Carpet 1 = thin loop pile commercial carpet. Carpet 2 = cut pile carpet with pad underlay Fracture Distribution and GIS The GIS maps of the revisited low-energy data (Powell et al. 2013) showed that fractures in the right parietal commonly initiated at its edge along the coronal suture. Other, less frequent initiation sites were concentrated on the ventral aspect, along the squamosal suture, and the caudal edge along the lamboidal suture bordering the occipital bone. Generally, these fractures were linear and grew inwards toward the center of the bone. There was substantial damage to the right frontal bone with many fractures aligning with those formed in the right parietal. The occipital bone was undamaged, except in one case (Figure 4.12a). As the GIS maps show, doubling the impact energy onto a rigid surface dramatically increased the extent of fracturing. Again, fracture initiation sites were concentrated 80 around the edges of the right parietal, forming along the coronal, squamosal, and lamboidal sutures (Figure 4.12b). Damage to the frontal bone was more pronounced than for the low energy impacts with many fractures crossing the coronal suture from the right parietal. However, as was the case in other groups, very few fractures formed along the sagittal suture. The most distinguishing characteristic of high energy impacts versus low energy impacts was the extensive damage to the occipital bone with many fractures passing through the center of the bone. In addition, a number of high-energy rigid specimens showed fractures in the left parietal bone; which initiated primarily at the lamboidal and coronal sutures and rarely at the saggital suture. There was only one case where a fracture crossed the sagittal suture between the left and right parietals (Figure 4.12b). Although the total amount of bone fracture was reduced, head impacts onto the thin carpet 1 interface developed fracture patterns similar to the high energy rigid group. Fracture initiation sites were concentrated in three areas: one at the center of the right coronal suture, another at the coronal-squamosal suture intersection, and a third on the lamboidal suture (Figure 4.12c). High energy carpet 2 impact damage patterns were similar to the low energy rigid impacts of Powell et al, after accounting for the reduced sample size. Compared to low energy rigid impacts, carpet 2 impacts showed similar initiation sites around the right parietal but with more occipital and less frontal bone damage as well as a few specimens with bilateral fractures (Figure 4.12d). 81 a b Figure 4.12. Geographic Information System maps of the a) low energy rigid b) high energy rigid c) high energy carpet 1 d) high energy carpet 2 impacts 82 Figure 4.12 Cont’d c d 83 Finite Element Modeling Analysis of the principal stress distributions in the domed shell under a centralized pressure loading showed that the largest stresses were tensile in nature and located around the periphery of the shell oriented in the circumferential direction. However, there were also large tensile stresses formed directly under the applied load (Figure 4.13). Note that only the relative magnitudes of these stresses are indicated. These data are intended only to illustrate the relative stress distribution under these conditions and are not representative of the loading levels required to initiate fracture. Figure 4.13. Maximum principal stresses generated in the shell as a result of the boundary and loading conditions of Figure 4.5. Arrows indicate the relative magnitude and orientation of the maximum principal stress with larger stresses signified by large, red arrows and smaller stresses by short, blue arrows. The largest tensile stresses were formed around the periphery of the shell. In addition, tensile stress increased directly under the central load application area. 84 DISCUSSION To better understand the effects of fall height on fracture patterns in the developing skull, the current study performed fall-induced impacts to the center of the right parietal bone of infant porcine heads. It was hypothesized that raising the impact energy versus a previous study by Powell (2013) would increase the frequency and length of fractures in the impacted right parietal bone. Doubling the drop height was found to increase the number of fracture initiation sites around the edge of the right parietal (Figure 4.12a, b) as well as the total length of fractures (Figure 4.8). When a crack forms, stress is relieved in the surrounding material, preventing the nucleation of subsequent cracks (Martin, Burr, and Sharkey 1998; Freund 1998). However, when the impact energy is increased and loading occurs more quickly, multiple cracks can form at remote sites before this stress shielding takes place (Martin, Burr, and Sharkey 1998). This effect was observed in a previous study by Powell et al. in which increasing the impact energy to an entrapped head increased the number of fracture initiation sites in the impacted right parietal (Powell et al. 2012). The patterns observed in the current study suggest the same is true in fall-induced impacts. In the impacted bone, increased impact energy from a higher fall height is dissipated by means of initiating new cracks and forming branches, as opposed to elongating simple linear cracks. Consistent with the observations of previous studies (Baumer et al. 2010; Powell et al. 2012; Powell et al. 2013), fractures in the impacted right parietal bone appeared to initiate at the suture boundaries and propagate inward toward the center of the bone (Figure 4.12). This effect has been previously attributed to “out-bending” of the skull 85 away from the impact point that generate large tensile stresses oriented circumferentially and produce radial oriented fractures in the skull (Gurdjian, Lissner, and Webster 1947; Gurdjian, Webster, and Lissner 1950). A study by Powell examined the stresses generated in a spherical model following a simulated drop onto a rigid surface (Powell 2010). He documented large tensile stresses formed at the point of impact, but the model did not actually examine the state of stress in the parietal bone itself. Finite element analysis performed in the current study modeled this bone as an isolated body supported around the edges while allowing for radial and circumferential expansions, as might be allowed by the surrounding compliant sutures (Figure 4.5). All of the maximum principal stresses from the analysis were tensile in nature and oriented primarily in a circumferential direction (Figure 4.13). Similar to the study of Powell (2010), large tensile stresses were formed directly under the impact site due to inbending of the shell under the load. However, in the current study, the largest tensile stresses were oriented circumferentially around the shell’s edge. This may explain the remote nature of the fractures observed in the current and previous studies (Baumer et al. 2010; Powell et al. 2012; Powell et al. 2013) (Figure 4.13). It is possible that impacting the center of the parietal bone generates circumferential tensile stresses sufficiently large to pull the bone apart around the edges, initiating radial cracks at the bone-suture boundary. In addition, there could be underlying geometrical or material features of the skull that help localize these stresses at the periphery of the parietal bone and cause fracture to initiate in consistent locations as observed in the distribution patterns above (Figure 4.12). However, such detail was not accounted for in the current modeling 86 effort. The current model is a gross simplification of the porcine skull anatomy and only meant to generate a potential explanation for fracture initiation around the periphery of this bone under a centralized impact to the bone. In reality, the mechanical properties and fracture processes are likely much more complex and beyond the scope of this study. The current investigation also hypothesized that increasing the drop energy would increase the level of damage to bones away from the impact site. This hypothesis was confirmed as the frequency of fracture in the occipital and right frontal bones was documented in the current study to increase with impact energy level (Figure 4.12a, b). As was the case with the impacted right parietal, fractures in bones away from the impact site always intersected the edge of the bone at a suture. In the occipital and right frontal bones such fractures were concentrated in regions adjacent to the right parietal. It stands to reason then that the insult applied to the right parietal was somehow transferred to the adjacent bones. Since sutures have been found to stiffen at high strain rates (Coats and Margulies 2006), it is possible that the circumferential stresses around the periphery of the right parietal bone illustrated above were transferred by the sutures to the surrounding bones. A second possible mechanism was that there was sufficient deformation of the parietal bone upon impact that the right frontal and occipital bones impacted the surface directly. In either case, the stresses in these neighboring bones would have then been highest in areas immediately adjacent to (i.e., sharing a suture with) the impacted bone. When combined with the discontinuity provided by the bone edge, fractures in the bones surrounding the impacted bone 87 would tend to form along the suture closest to the impact site as observed in the current study. In addition to producing fractures in the occipital and right parietal bones, the current study documented a number of specimens with fractures in the left parietal, opposite the impact site. Such bilateral fractures have previously been used as a requisite for diagnosing multiple blows to the head, an injury suggestive of abusive trauma (Meservy et al. 1987; Billmire and Myers 1985). However, their prevalence in the current study suggests that bilateral fractures can form as a result of a single impact from a high fall and may be a poor characteristic for distinguishing abuse from falls. To analyze the effect of impact interface on the pattern of cranial fractures, the current study also performed controlled high energy head drops onto two different carpeted interfaces. Drops to the more compliant carpet resulted in significantly higher peak impact forces compared to the commercial-grade carpet and rigid aluminum surface (Figure 4.9). This trend was likely due to the softer carpet and padding effectively increasing the area of contact over the skull. The area of contact may have been large enough that the impact force was distributed to other bones of the skull before stresses grew large enough in the right parietal to cause fracture, if fracture occurred. The commercial carpet interface was likely not sufficiently complaint to produce this same distribution, resulting in peak forces that were similar to rigid surfaces. The final hypothesis of the study was that increasing the impact interface compliance would have the same effect on the fracture patterns as decreasing the drop energy. It 88 was documented that impacts onto a carpeted surface, even if the carpet provided minimal compliance, resulted in significantly less total skull bone fracture across all ages when compared with impacts to a rigid interface (Figure 4.11, Figure 4.12b, c, d). For the more compliant carpet, the degree of fracture matched levels from the low energy rigid drops (Figure 4.12a, d). Since bone fracture is a means of dissipating impact energy, it follows that energy absorbed through the deformation of a compliant interface would have reduced the impact energy imparted to the skull and distributed the impact force, both lessening the bone stresses to reduce the degree of cranial fracture. Interestingly, high-energy impacts to all surfaces had a consistent frequency of diastatic fracture which was lower than that of low-energy impacts to a rigid surface. This trend may be due to the fact that sutures tend to stiffen at high strain rates (Coats and Margulies 2006). In the current study, the high-energy impacts may have increased the loading rates and stiffened the sutures regardless of the interface compliance, reducing the degree of diastatic fracture compared to the slower loading rates of the low-energy impacts. Thus, the presence of diastatic fractures may be a good indicator of the height of the fall regardless of interface type. In summary, this study confirmed that increasing the impact energy through a raised fall height increases the incidence of cranial bone fracture in controlled drop tests of infant porcine heads. Upon further examination, this extra fracturing was shown to be generated by the formation and branching of new fractures as a means of dissipating the additional impact energy (Figure 4.12a, b). Furthermore, the type of impact 89 interface was shown to have a significant effect on the degree of cranial bone fracture with compliant carpeted surfaces attenuating some of the impact energy and distributing the impact force over the skull to help limit the extent of cranial fracturing. For a highly compliant interface, this attenuation was found to have an effect on the total bone fracture length that was similar to reducing the fall height. Finally, a new theoretical model suggested that radial fractures in the impacted parietal bone may have been generated by large circumferential tensile stresses developed at the edge of the bone. A major limitation of this study was that the impact energy levels were based on a previous study that increased drop height with age to generate consistent fractures. A suggestion for future work would be to keep the energy constant across all ages in order to isolate the effect of age and skull development on fracture formation. Using the results from the current, previous, and future studies, the authors hope to build a generalized database for correlating an impact scenario with a fracture pattern for the developing porcine skull. Such knowledge could then be applied to patterning human case studies to develop a tool for aiding medical examiners and forensic scientists in distinguishing accident from abuse, potentially by utilizing techniques commonly used in machine learning. Development of such a human database and machine learning algorithm may be able to associate the resultant pattern of pediatric cranial fracture patterns with the impact load scenario to help determine injury causation in a quantitative manner. 90 ACKNOWLEDGEMENTS This study was funded by a grant from the National Institute of Justice. The author would like to thank Mr. Ed Reed and Ms. Star Lewis (of Reed-McKenzie Farms, Decatur, MI) for supplying the porcine specimens and Mr. Cliff Beckett for his technical support. 91 APPENDIX 92 Table 4.1. Revisited raw data from low-energy drops to a rigid surface (Powell et al. 2013). Specimen Age (days) Mass (g) HD001 HD002 HD003 HD004 HD005 HD006 HD007 HD008 HD009 HD010 HD011 HD012 HD013 HD014 HD015 HD016 HD017 HD018 HD019 HD020 HD021 HD022 HD023 HD024 HD025 HD026 HD027 HD028 HD029 HD030 HD031 5 5 2 9 5 10 10 9 13 13 8 8 15 15 16 16 6 6 14 14 3 3 17 4 4 7 7 11 11 12 12 500 500 350 555 500 500 500 555 990 990 430 430 800 800 870 870 395 395 720 720 380 380 685 375 375 350 350 605 605 635 635 Impact Max Duration Load (N) (ms) 320 751 587 639 662 906 808 728 832 949 718 790 1226 1251 1341 1193 800 711 1439 889 777 689 1208 1228 790 1121 990 1227 987 1247 1190 3.1 1.7 3.3 2.1 1.3 1.3 2.3 2.2 2.0 0.2 2.0 0.8 2.4 1.9 1.1 1.6 4.5 2.2 1.2 2.7 3.7 4.7 1.6 1.9 2.1 2.5 1.6 1.4 1.3 4.0 1.7 93 Bone Fracture (mm) Diastatic Fracture (mm) 7 0 3 0 18 0 61 14 29 33 13 22 51 11 30 16 12 17 57 19 66 15 12 15 19 43 13 26 24 38 16 18 0 29 0 0 0 10 0 0 0 0 0 20 0 0 0 25 10 0 11 10 0 0 0 0 49 15 0 23 22 17 Total Fracture Length (mm) 25 0 32 0 18 0 71 14 29 33 13 22 71 11 30 16 37 27 57 30 76 15 12 15 19 92 28 26 47 60 33 Table 4.2. Raw data from high energy drop impacts to a rigid surface. Specimen Age (days) Mass (g) Max Load (N) Impact Duration (ms) Bone Fracture (mm) Diastatic Fracture (mm) Total Fracture (mm) HD101 HD102 HD103 HD104 HD105 HD106 HD107 HD108 HD109 HD110 HD111 HD112 HD113 HD114 HD115 HD116 HD117 HD118 HD119 HD120 HD121 HD122 HD123 HD124 HD125 HD126 HD127 HD128 HD129 HD130 HD131 HD132 6 5 8 6 5 11 7 3 8 10 13 13 15 18 7 17 8 14 6 16 3 9 13 13 15 4 8 9 7 9 14 12 380 370 450 354 385 613 423 445 490 550 636 720 615 605 520 785 467 691 317 607 475 587 576 592 641 456 496 530 423 397 596 589 1073 923.1 1271 1014 1040 1410 1203 1136 1253 1302 1654 1668 1505 1366 1248 1493 1380 1606 1078 1696 528.6 1064 1445 1652 1501 1353 1139 1400 1377 1237 1553 1662 3.8 6.3 5.7 6.2 7.6 4.7 5.9 5.4 5.7 6.2 5 4.4 6 6.9 4.3 5.7 4.3 5.2 4.9 4.9 6.9 4.6 4.4 4.2 5.6 4.3 4.4 4.4 4.2 4.2 3.4 3.7 190 291 64 54 178 157 89 119 82 142 279 33 263 179 110 53 244 120 145 36 14 19 118 93 65 45 63 152 132 158 85 65 0 0 0 0 0 8 0 37 0 0 0 0 0 16 0 0 0 23 0 0 0 0 0 0 0 0 0 0 0 9 0 0 190 291 64 54 178 165 89 156 82 142 279 33 263 195 110 53 244 143 145 36 14 19 118 93 65 45 63 152 132 167 85 65 94 Table 4.3. Raw data from high-energy drop impacts to a compliant surface (Carpet 1). Specimen Age Mass (g) HD161 HD162 HD163 HD164 HD165 HD166 HD167 HD168 HD169 HD170 HD171 HD172 HD173 HD174 HD175 HD176 HD177 HD178 HD179 HD180 HD181 HD182 HD183 HD184 HD185 HD186 HD187 HD191 HD192 HD193 HD194 HD195 9 18 5 15 12 18 16 7 9 8 4 8 7 16 3 16 7 9 5 6 9 16 15 17 5 18 6 16 15 14 14 11 590 789 399 652 575 603 789 474 595 647 420 419 395 603 320 683 450 493 515 516 496 746 715 630 382 694 535 700 611 594 578 550 Max Load (N) Impact Duration (ms) Bone Fracture (mm) Diastatic Fracture (mm) 1415 1661 1172 1750 1702 1905 1986 1133 1210 1234 1131 1115 1186 1827 950.5 N/A 1045 1227 1504 1320 1367 1738 1615 1655 1074 1800 1471 1739 1186 1600 1159 1459 4.9 6.6 5.5 5.1 5.1 5 5.2 3.9 3.8 5.7 4.3 4.7 3.9 4.5 6.1 N/A 4.6 6.2 4.2 3.8 3.9 2.5 4.8 5.3 4 6.3 3.4 5.7 7.2 5.1 4.2 4.3 9 18 58 146 68 50 26 0 32 47 103 57 133 151 140 50 35 48 4 54 73 39 18 11 55 3 22 36 9 2 66 42 0 0 0 0 26 0 0 0 0 0 0 0 10 0 22 13 0 0 0 0 0 15 0 0 0 0 0 0 0 0 0 0 95 Total Fracture Length (mm) 9 18 58 146 94 50 26 0 32 47 103 57 143 151 162 63 35 48 4 54 73 54 18 11 55 3 22 36 9 2 66 42 Table 4.4. Raw data from high-energy drop impacts to a compliant surface (Carpet 2). Specimen Age (days) Mass (g) Max Load (N) Impact Duration (ms) Bone Fracture (mm) Diastatic Fracture (mm) HD141 HD142 HD143 HD144 HD145 HD146 HD147 HD148 HD149 HD150 HD151 HD152 HD153 HD154 HD155 HD156 HD157 HD158 HD159 HD160 HD188 HD189 HD190 12 15 19 17 3 2 5 9 8 18 6 13 2 14 3 15 4 13 17 9 7 9 4 575 694 860 862 314 363 388 494 445 758 393 613 355 650 395 701 370 707 792 522 467 423 420 1854 1856 1683 1965 1460 1466 1390 1458 1785 2281 1160 2113 1863 2147 1724 2427 1269 1956 1902 1444 1763 1145 1098 5.3 4.3 4.5 4.5 4.8 5.9 6.6 6.9 3.6 7.2 5.6 4.8 4.7 3.6 3.1 3.5 6.1 5.2 4.5 6.1 2.8 6.5 6.2 68 0 0 0 17 60 8 27 0 0 26 20 101 0 12 5 19 26 7 21 0 33 17 0 0 0 0 0 18 0 0 0 0 0 17 0 0 0 0 0 0 0 0 0 21 0 96 Total Fracture Length (mm) 68 0 0 0 17 78 8 27 0 0 26 37 101 0 12 5 19 26 7 21 0 54 17 REFERENCES 97 REFERENCES Baumer, Timothy G., Nicholas V. Passalacqua, Brian J. Powell, William N. Newberry, Todd W. Fenton, and Roger C. Haut. 2010. "Age-dependent fracture characteristics of rigid and compliant surface impacts on the infant skull - a porcine model". Journal of Forensic Sciences 55 (4). Baumer, Timothy G., Brian J. Powell, Todd W. Fenton, and Roger C. Haut. 2009. "Age dependent mechanical properties of the infant porcine parietal bone and a correlation to the human". Journal of Biomechanical Engineering-Transactions of the Asme 131 (11). Billmire, M. E., and P. A. Myers. 1985. "Serious head injuries in infants - Accident or abuse". Pediatrics 75 (2). Case, M. E., M. A. Graham, T. C. Handy, J. M. Jentzen, J. A. Monteleone, and Co Natl Assoc Med Examiners Ad Hoc. 2001. "Position paper on fatal abusive head injuries in infants and young children". American Journal of Forensic Medicine and Pathology 22 (2):112-122. Coats, B, and Susan S Margulies. 2003. "Characterization of Pediatric Porcine Skull Properties During Impact". IRCOBI Conf.:57-66. Coats, Brittany, and Susan S Margulies. 2008. "Potential for head injuries in infants from low-height falls: Laboratory investigation". Journal of Neurosurgery: Pediatrics 2 (5):321-330. Coats, Brittany, and Susan S. Margulies. 2006. "Material properties of human infant skull and suture at high rates". Journal of Neurotrauma 23 (8). Freund, L Benjamin. 1998. Dynamic fracture mechanics: Cambridge university press. Greenes, D. S., and S. A. Schutzman. 1999. "Clinical indicators of intracranial injury in head-injured infants". Pediatrics 104 (4). Gurdjian, E. S., J. E. Webster, and H. R. Lissner. 1950. "The mechanism of skull fracture". Radiology 54 (3). 98 Gurdjian, ES, HR Lissner, and JE Webster. 1947. "The mechanism of production of linear skull fracture; further studies on deformation of the skull by the stresscoat technique". Surgery, Gynecology & Obstetrics 85 (2):195. Gurdjian, ES, JE Webstef, and HR Lissner. 1953. "Observations on prediction of fracture site in head injury". Radiology 60 (2):226-235. Hall, John R, Hernan M Reyes, Maria Horvat, Janet L Meller, and Robert Stein. 1989. "The mortality of childhood falls". The Journal of Trauma and Acute Care Surgery 29 (9):1273-1275. Hobbs, C. J. 1984. "Skull fracture and the diagnosis of abuse". Archives of Disease in Childhood 59 (3). Lam, W. H., and A. Mackersie. 1999. "Paediatric head injury: incidence, aetiology and management". Paediatric Anaesthesia 9 (5):377-385. Martin, R Bruce, David B Burr, and Neil A Sharkey. 1998. Skeletal tissue mechanics: Springer Verlag. Meservy, C. J., R. Towbin, R. L. McLaurin, P. A. Myers, and W. Ball. 1987. "Radiographic characteristics of skull fractures resulting from child abuse". American Journal of Roentgenology 149 (1):173-175. Powell, Brian J. 2010. Fracture Patterns of the developing skull attributable to different impact scenarios [MS Thesis]. Michigan State University, East Lansing, MI. Powell, Brian J, Nicholas V Passalacqua, Todd W Fenton, and Roger C Haut. 2013. "Fracture characteristics of entrapped head impacts versus controlled head drops in infant porcine specimens". Journal of Forensic Sciences. Powell, Brian J., Nicholas V. Passalacqua, Timothy G. Baumer, Todd W. Fenton, and Roger C. Haut. 2012. "Fracture patterns on the infant porcine skull following severe blunt impact". Journal of Forensic Sciences 57 (2). Prange, M. T., B. Coats, A. C. Duhaime, and S. S. Margulies. 2003. "Anthropomorphic simulations of falls, shakes, and inflicted impacts in infants". Journal of Neurosurgery 99 (1):143-150. Reece, R. M., and R. Sege. 2000. "Childhood head injuries - Accidental or inflicted?". Archives of Pediatrics & Adolescent Medicine 154 (1). 99 Schutzman, S. A., P. Barnes, A. C. Duhaime, D. Greenes, C. Homer, D. Jaffe, R. J. Lewis, T. G. Luerssen, and J. Schunk. 2001. "Evaluation and management of children younger than two years old with apparently minor head trauma: Proposed guidelines". Pediatrics 107 (5). Weber, W. 1984. "Experimental study of skull fracture in infants". Zeitschrift Fur Rechtsmedizin-Journal of Legal Medicine 92 (2). ------. 1985. "Biomechanical fragility of skull fractures in infants". Zeitschrift Fur Rechtsmedizin-Journal of Legal Medicine 94 (2). 100 CHATPER FIVE PATTERNS OF FRACTURE DUE TO THREE-POINT BENDING IN LONG BONES INTRODUCTION Clinical observations show that long bones failed under bending forces often fracture in a characteristic “butterfly” pattern resembling a Y-shaped profile with the removal of a wedge-shaped piece of bone referred to as a butterfly fragment (Figure 5.1) (Carter and Spengler 2002; Martin, Burr, and Sharkey 1998; Galloway 1999). The reason for this unique pattern is explained in the current literature using basic solid mechanics and bone strength asymmetry; with cortical bone being strongest in compression, weakest in shear, and intermediate in tension (Martin, Burr, and Sharkey 1998). The fracture is thought to initiate as a tensile failure opposite the impacted side in the form of a transverse crack. Once the crack reaches the compressed side of the neutral or centroidal axis, with the stresses not large enough to cause compressive failure, it has been suggested that shear stresses along the principal shear planes (oriented at 45° angles) exceed the bone’s shear strength, causing the crack to bifurcate and propagate in those directions (Figure 5.1) (Sharir, Barak, and Shahar 2008; Carter and Spengler 2002). Following this theory, the presence of a butterfly fragment is often used to diagnose the direction of impact loading for a bending fracture (Carter and Spengler 2002; Martin, Burr, and Sharkey 1998; Passalacqua and Fenton 2012; Symes et al. 2012; Rockhold and Hermann 1999). Such evidence of impact direction can be a useful tool for medical examiners or forensic scientists attempting to reconstruct a traumatic event. 101 Figure 5.1. Diagram of a typical butterfly fracture of a long bone as a result of bending loads. An extensive series of transverse impacts to human leg bones by Kress et al. found that butterfly fractures, although prevalent, were not the exclusive form of gross bone fracture due to bending (Kress et al. 1995). The study reported a high incidence of oblique and some transverse oriented fractures. Such patterns are traditionally associated with failure due to axial compression or tension, respectively (Kress et al. 1995; Carter and Spengler 2002). Additionally, photographs of fractures from that study show butterfly fragments that were much larger than traditionally depicted and having a branch point that appeared to be located on the tensile side of the bone, forming an inverted V-, as opposed to a Y-, shaped pattern (Kress et al. 1995). Kress et al. and other studies have reported experimental cases of the butterfly fragment forming in the reversed orientation with the fragment opening toward the tensile side (Thomas and Simmons 2011; Martens et al. 1986; Kress et al. 1995). This inconsistent fragment orientation was so prevalent in some of those studies that it was concluded that butterfly fragment orientation was not a good indicator of impact direction (Thomas and Simmons 2011; Martens et al. 1986). A study of bending fracture 102 in dry human bones was performed in response to this inconsistency of butterfly fragment orientation and found that incomplete hairline fractures gave a better indication of loading direction than gross fracture pattern alone (Fenton et al. 2012). When these findings were applied to the work of Thomas and Simmons (2011), the revisited fractures were matched with the correct loading orientation in 100% of the specimens (Reber 2013). The hypothesis of the current study was that, through controlled three-point bending tests, simulated impact to an upright leg with an applied normal body weight would produce a butterfly fracture pattern that could consistently be used to determine impact loading direction. METHODS Six pairs of human legs were dislocated at the head of the femur and stored at -20°C. All specimens were male aged 52 to 65 years. The legs were thawed at room temperature two days prior to testing. The femora were then cleaned of all soft tissue to a point just above the knee. The proximal end of each femur was potted in a polyester resin (Martin Senour Fibre Strand Plus 6371, Sherwin-Williams; Cleveland, OH). A similar technique was used to pot around the distal portion of the bone just above the knee, as the knee joint itself was to be used in another set of experiments (Figure 5.2). The exposed section of bone was wrapped with gauze saturated with phosphate-buffered saline solution and kept moist throughout the cleaning and potting procedures to avoid potential changes in the fracture pattern due to drying (Burstein et al. 1972). 103 Bone Wrapped in Saline-Soaked Gauze Potting Cup Molds Polyester Resin Potting Material Figure 5.2. Technique used to pot the bone ends with the femur passing through the bottom cup to leave the knee intact. The resin pots served as attachment points to install the bone into a 3-point bending fixture mounted on a servohydraulic material testing machine (MTS; Eden Prairie, MN) (Figure 5.3). The fixture allowed both ends to pivot freely and an X-Y table allowed translations of the distal end of the bone. Two large springs applied a static axial preload of 450 N to simulate an upright standing posture with the leg supporting one-half normal body weight. A 30 mm-diameter solid steel rod served as the impact anvil and was oriented perpendicular to the bone axis and a preload of 50 N was applied to the bone mid-diaphysis to eliminate any potential residual system compliance. Failure was induced by a position-controlled 2 Hz, 20 mm haversine displacement of the anvil into the bone. Force and displacement data were recorded at 10,000 Hz by an actuatormounted load transducer (3210AF-5K, 5000lb capacity, Interface; Scottsdale, AZ) and 104 linear variable differential transformer (LMT-711P35, ±3.5” stroke, G.L. Collins Corp.; Long Beach, CA). Energy absorbed by the bones to failure was calculated as the area under the load-displacement curve up to the peak load. Right legs were oriented so that the anvil loaded the posterior surface of the femur and left legs loaded the anterior surface of contralateral pairs of femora from each cadaver. Actuator Impact Anvil Load cell Axial Preload Spring X-Y Translation Table Pivoting Potting Cups Figure 5.3. Three-point bending setup on the servohydraulic testing machine. After fracture, potting was removed and the femora were boiled to remove remaining oils and soft tissue. The bones were then reassembled to photographically document fracture patterns. Each fracture was categorized using the following convention, based on the primary characteristic of the complete fracture orientation: transverse, oblique, 105 or butterfly. In addition, the presence and orientation of any incomplete “hairline” fractures were also documented in the current study. Any potential differences in maximum load and displacement at failure were documented for orientation dependence using a paired t-test between the left and right femora for each specimen. Results were considered statistically significant for a p-value less than 0.05. RESULTS Gross fracture assessments on these isolated human femora yielded 7 transverse fractures and 5 oblique fractures. In all femora at least one hairline fracture occurred. In grossly oblique fractures the hairline crack branched off away from the initial transverse fracture on the tensile side of the bone, resulting in what was classified as an “incomplete butterfly” pattern (Figure 5.4a). The branch point of these hairline fractures were always located on the tensile side of the neutral axis and the cracks extended from the tensile, transverse oriented crack and propagated toward the impacted surface of the bone (Figure 5.4a). In cases of grossly transverse patterns of fracture, hairline fractures would curve predictably away from the transverse crack and propagate along the long axis of the bone (Figure 5.4b). 106 Gross Oblique Fracture Forming Incomplete Butterfly a Incomplete Hairline Crack Crack curves to parallel long b Gross Transverse Fracture Figure 5.4. Oblique (a) and Transverse (b) gross fracture patterns. The anvil impacted the top surface of the bone at the site marked by the arrow, initiating tensile failure on the bottom surface. Notice that the hairline fractures branch away from the gross fracture on the tensile side (bottom half) of the bone. The load-displacement plots to failure of the bone were rather linear with a slightly curved orientation downward, but failure was abrupt in all cases (Figure 5). Time to bone failure was less than 150 milliseconds for all bones. Failure loads ranged from 3.99 to 9.20 kN with a mean of 6.10 kN and standard deviation of 1.70 kN. Failure displacement ranged from 5.83 to 13.24 mm with a mean of 9.62 and a standard deviation of 2.25 mm. Failure energy ranged from 14.96 to 66.30 J with a mean of 34.21 107 J and a standard deviation of 15.71 J (Table 5.1). There was no significant difference in the failure load (P=0.257), displacement (P=0.223), or energy to failure (P=0.242) between loading orientations. Bending Load (N) 6000 5000 4000 3000 2000 1000 0 0 2 4 6 Displacement (mm) 8 10 Figure 5.5. Sample load-displacement plot of a three-point bending test of a femur up to failure Table 5.1. Comparison of fracture characteristics for paired bones. Left femora were loaded on the anterior surface and right femora were loaded on the posterior surface. When identifying the gross fracture type: T=transverse, O=oblique Specimen 1 2 3 4 5 6 Left Right Left Right Left Right Left Right Left Right Left Right Failure Load (N) Failure Displacement (mm) Time to Failure (ms) Failure Energy (J) Gross Fracture Type 4714 6007 7325 9195 4501 5113 5642 5322 5325 3992 7050 9001 5.83 10.14 6.41 10.90 7.13 8.62 12.04 9.81 10.83 9.72 10.75 13.24 86.5 122 92.3 192.2 97.5 109.3 137.4 120.3 128.5 118.1 128.1 146.5 14.96 36.95 24.63 57.82 17.20 24.05 40.87 29.01 34.54 22.50 41.71 66.30 T O T T O T T T O T O O 108 DISCUSSION The current study conducted controlled three-point bending tests to fracture human femora and examined the resultant fracture patterns. While none of the specimens resulted in a fracture pattern as shown in Figure 5.1 with a complete wedge of bone being removed on the impacted side of the bone, 5/12 cases yielded oblique fractures with hairline cracks that suggested an incomplete butterfly wedge fracture pattern. In these cases, the orientation of the incomplete wedge was consistent with the loading orientation, making the direction of impact clear. In the remaining specimens, gross transverse bone fractures were generated; but in each case, characteristic hairline cracks were formed that branched and curved off the transverse-oriented fracture that indicated the direction of loading. Thus, the hypothesis of the current study was validated by the generation of these characteristic hairline crack patterns in all specimens. The consistent location of the incomplete butterfly fragment branch point on the tensile side of the neutral (centroidal) axis of the bone would contradict the solid mechanicsbased mechanism of bone fracture often quoted in the engineering literature, which attributes the branching and angling of a butterfly fracture to shear failure of bone on the compression side. Other studies on bone fracture mechanics have performed crack propagation tests on compact and single edge notched specimens machined from bovine femoral and tibial cortical bone (Melvin and Evans 1973; Behiri and Bonfield 1989; Melvin 1993). These studies have documented that when an initial notch is made in a transverse orientation (perpendicular to the bone’s long axis), tensile loading at 109 rates that result in a catastrophic failure of the bone sample produces crack bifurcation with two cracks rapidly propagating diagonally away from the initial notch (Melvin and Evans 1973; Melvin 1993; Behiri and Bonfield 1989) (Figure 5.6). These findings matched the whole-bone patterns of the current study that occurred on the tensile side of the bones. However, no mechanism was given in the above studies to explain the reason for this branching. Osteon Orientation Tensile Load Application Points Diagonal Crack Branch Figure 5.6. Fracture profile of compact tensile test specimen with an initial crack made transverse to the long axis of the bone. Redrawn from Behiri and Bonfield, 1989 Dynamic fracture propagation studies on brittle materials have shown branching patterns resembling river deltas for cracks in plates under pure tension (mode 1) loading (Ha and Bobaru 2010; Ramulu and Kobayashi 1985; Ravi-Chandar and Knauss 1984). A proposed mechanism for branching involves the formation of microcracks ahead of the main crack tip, directing propagation and, occasionally, causing instabilities that produce bifurcation (Ravi-Chandar and Knauss 1984, 1984). Since bone behaves as a brittle solid at high rates of loading (Bonfield 1987), this dynamic fracture theory lends itself nicely to the current study. In addition, the cellular microstructure of bone contains various 110 lacunae and other defects that could have aided in the development of microcracks ahead of the propagating fracture (Vashishth, Tanner, and Bonfield 2000). A second characteristic observed in the current study was that fractures often curved to parallel the long axis of the bone. This curving mirrors the findings of bone fracture mechanics research on machined specimens which found that, regardless of the initial crack angle and loading direction, fractures propagated parallel to the bone’s long axis, along the orientation of osteons (Behiri and Bonfield 1989). These same fracture mechanics tests have found that the fracture toughness of compact bone is significantly lower when the crack is propagated parallel to, as opposed to across, the bone’s long axis (Melvin and Evans 1973; Melvin 1993; Behiri and Bonfield 1989). Such transversely isotropic material properties are a characteristic of bone resulting from its cellular nature and the alignment of osteons along the long axis of the bone (Reilly and Burstein 1975) (Figure 5.7). Thus, the curving observed in the current study may be a result of cracks following the path of least resistance, running between osteons instead of across them as they propagated through the bone. 111 Figure 5.7. Composition of cortical bone showing the orientation of osteons parallel to the long axis. In addition to material anisotropy, loading rate and energy level have been shown to influence fracture patterns by changing the degree of comminution in the fractured bone (Beardsley et al. 2002). Such effects were observed in the above fracture mechanics studies; transverse specimens loaded slowly (0.01 mm/min) steadily propagated a single crack through the material. However, increasing the loading rates to 0.02 mm/min was enough to destabilize the crack propagation, causing catastrophic failure and branching of the crack (Behiri and Bonfield 1989). In the current study, loading rates approached 5000 mm/min. This excess impact energy was likely transferred to the bone through the initiation of multiple fractures. However, once these cracks had formed, the strain energy stored in the deformed bone had been released. After this unloading, it is possible that the speed of the impact anvil provided insufficient energy to propagate all the cracks simultaneously and only one continued on to cause gross fracture. 112 In conclusion, the current study showed that, although not identical to the fracture pattern described in biomechanics textbooks, 3-point bending of human femora produced characteristic patterns of fracture that are helpful in diagnosing impact direction. Grossly oblique fractures were formed with an incomplete hairline fracture that branched at an angle and the two cracks diverged as they propagated toward the impact site. Transverse gross fractures were also always accompanied with incomplete cracks that branched off at an angle. Another common characteristic was that gross and incomplete fractures initiated at a diagonal but curved to parallel the long axis as they moved away from the tensile initiation site. By observing these patterns of branching, diverging, and curving the loading orientation was made clear in all specimens. However, the explanation for the formation of these fractures given in biomechanical texts was inconsistent with the patterns observed in the current study. A static, strength-of-materials explanation for fracture branching on the compressive side of the bone under 3-point bending is likely an over-simplification of the complex mechanisms involved in fracture initiation and propagation. A new explanation, based on previous studies of dynamic fracture propagation and branching of a crack in brittle, transversely isotropic bone was introduced here and better explained the mechanics behind the formation of butterfly fractures initiating on the tensile side of long bones under 3-point bending. 113 APPENDIX 114 Figure 5.8. Specimen 12-1441L-AP. Top: lateral view. Bottom: medial view. For all photographs, the bone was oriented with the anvil striking the top surface Figure 5.9. Specimen 12-1441R-PA. Top: lateral view. Bottom: medial view. 115 Figure 5.10. Specimen 12-1447L-AP. Top: lateral view. Bottom: medial view. Figure 5.11. Specimen 12-1447R-PA. Top: lateral view. Bottom: medial view. 116 Figure 5.12. Specimen 12-1466L-AP. Top: lateral view. Bottom: medial view. Figure 5.13. Specimen 12-1466R-PA. Top: lateral view. Bottom: medial view. 117 Figure 5.14. Specimen 13-1026L-AP. Top: lateral view. Bottom: medial view. Figure 5.15. Specimen 13-1026R-PA. Top: lateral view. Bottom: medial view. 118 Figure 5.16. Specimen 13-1065L-AP. Top: lateral view. Bottom: medial view. Figure 5.17. Specimen 13-1065R-PA. Top: lateral view. Bottom: medial view. 119 Figure 5.18. Specimen 13-1180L-AP. Top: lateral view. Bottom: medial view. Figure 5.19. Specimen 13-1180R-PA. Top: lateral view. Bottom: medial view. 120 REFERENCES 121 REFERENCES Beardsley, Christina L, Christine R Bertsch, J Lawrence Marsh, and Thomas D Brown. 2002. "Interfragmentary surface area as an index of comminution energy: proof of concept in a bone fracture surrogate". Journal of Biomechanics 35 (3):331-338. Behiri, JC, and W Bonfield. 1989. "Orientation dependence of the fracture mechanics of cortical bone". Journal of Biomechanics 22 (8):863-872. Bonfield, W. 1987. "Advances in the fracture mechanics of cortical bone". Journal of Biomechanics 20 (11):1071-1081. Burstein, Albert H, John D Currey, Victor H Frankel, and Donald T Reilly. 1972. "The ultimate properties of bone tissue: the effects of yielding". Journal of biomechanics 5 (1):35-44. Carter, DR, and DM Spengler. 2002. "Biomechanics of fracture". Bone in Clinical Orthopaedics: a Study in Comparative Osteology:261-285. Fenton, Todd W, Ashley E Kendell, Trevor S DeLand, and Roger C Haut. 2012. "Determination of impact direction based on fracture patterns in human long bones". Paper read at American Academy of Forensic Sciences, at Atlanta, GA. Galloway, Alison. 1999. Broken bones: anthropological analysis of blunt force trauma: Charles C. Thomas Springfield. Ha, Youn Doh, and Florin Bobaru. 2010. "Studies of dynamic crack propagation and crack branching with peridynamics". International Journal of Fracture 162 (1-2):229-244. Kress, TA, DJ Porta, JN Snider, PM Fuller, JP Psihogios, WL Heck, SJ Frick, and JF Wasserman. 1995. "Fracture patterns of human cadaver long bones". Paper read at Proceedings of the 1995 International IRCOBI Conference on the Biomechanics of Impact, September 13-15, 1995, Brunnen, Switzerland. Martens, M, Remy Van Audekercke, P De Meester, and JC Mulier. 1986. "Mechanical behaviour of femoral bones in bending loading". Journal of Biomechanics 19 (6):443454. 122 Martin, R Bruce, David B Burr, and Neil A Sharkey. 1998. Skeletal Tissue Mechanics: Springer Verlag. Melvin, John W. 1993. "Fracture mechanics of bone". Journal of Biomechanical Engineering 115 (4B):549. Melvin, JW, and FG Evans. 1973. "Crack propagation in bone". Paper read at Biomechanics Symposium ASME. Passalacqua, Nick V, and Todd W Fenton. 2012. "Developments in Skeletal Trauma: Blunt-Force Trauma". In A Companion to Forensic Anthropology, edited by D. C. Dirkmaat: Blackwell Publishing Ltd. Ramulu, M, and AS Kobayashi. 1985. "Mechanics of crack curving and branching—a dynamic fracture analysis". International Journal of Fracture 27 (3-4):187-201. Ravi-Chandar, K, and WG Knauss. 1984. "An experimental investigation into dynamic fracture: II. Microstructural aspects". International Journal of Fracture 26 (1):65-80. ------. 1984. "An experimental investigation into dynamic fracture: III. On steady-state crack propagation and crack branching". International Journal of Fracture 26 (2):141154. Reber, Samantha L. 2013. "Interpreting the injury mechanisms of blunt force trauma from butterfly fracture patterns". Paper read at American Academy of Forensic Sciences, at Washington, DC. Reilly, Donald T, and Albert H Burstein. 1975. "The elastic and ultimate properties of compact bone tissue". Journal of Biomechanics 8 (6):393-405. Rockhold, Lauren A, and Nicholas P Hermann. 1999. "A Case Study of a Vehicular Hitand-Run Fatality: Direction of Force". In Broken Bones: Anthropological Analysis of Blunt Force Trauma, edited by C. C. Thomas. Springfield, IL. Sharir, Amnon, Meir Max Barak, and Ron Shahar. 2008. "Whole bone mechanics and mechanical testing". The Veterinary Journal 177 (1):8-17. Symes, Steven A, Ericka N L'Abbé, Erin N Chapman, Ivana Wolff, and Dennis C Dirkmaat. 2012. "Interpreting Traumatic Injury to Bone in Medicolegal Investigations". In A Companion to Forensic Anthropology, edited by D. C. Dirkmaat: Blackwell Publishing Ltd. 123 Thomas, Tammy S, and Tal Simmons. 2011. "The relationship between directionality of force and the formation of butterfly fractures". Paper read at American Academy of Forensic Sciences: Annual Scientific Meeting, at Chicago, IL. Vashishth, D., K. E. Tanner, and W. Bonfield. 2000. "Contribution, development and morphology of microcracking in cortical bone during crack propagation". Journal of Biomechanics 33 (9):1169-1174. 124 CHAPTER SIX CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK This thesis has documented a body of work on the development and fracture mechanics of bone. A study on the effect of exercise on whole-bone mechanics was introduced based on a current study of chicken bones from various commercial egg-laying housing systems. A subsequent study then developed a method for predicting these mechanical properties based on computed tomography (CT) data along with FE and linear regression models based on simple beam theory. Finally, the patterns of bone fracture resulting from various traumas were explored and the interpretation of these patterns and how they pertain to forensic investigations was discussed. In Chapter 2, whole-bone mechanical tests documented that pullets reared in aviary environments that had more room for movement and exercise developed bones that were, on the whole, stiffer and stronger than those from birds housed in cages with restricted movement. The differences in mechanical properties were attributed to changes in the cortical geometry as there was little difference in the bone material properties between environments. In addition, the differences between environments were found to be more pronounced in the humerus than the tibia. This result was attributed to the type of activity each environment might have allowed. While conventional cages permitted walking and mechanical loading of the legs, there was likely insufficient room for exercising the wings. However, the open aviary environment 125 allowed the birds to fly, an action that was likely responsible for the more dramatic changes observed in the humerus compared to the tibia. Chapter 3 extended the work of Chapter 2 in developing a method to assess the mechanical strength and stiffness of laying hen bones in non-destructive in vivo mechanical tests. Using high-resolution micro-CT, accurate measurements of bone geometry and density distribution were obtained. Finite element (FE) analysis was used to generate linear elastic models of each bone and used an inverse analysis to find the optimum relationship between micro-CT density and Young’s modulus. Using this relationship, the stiffness of the FE models matched experimental results with reasonable accuracy. A second technique using the micro-CT measurements for density and geometry along with simple beam theory was shown to predict whole bone stiffness with the same degree of accuracy as the FE model. Similarly, these data were used to predict whole-bone failure strength, which produced excellent correlations with experimental data. Due to the similar performances of the FE-based and beam theory models, it was concluded that it was unnecessary to account for density and material variation within a bone as a simple measure of average bone density was sufficient. A suggestion for future work would be to refine the beam theory-based model for bone failure strength and improve its predictive accuracy. The use of these techniques in an in vivo scenario seems appropriate and needs to be validated in future studies once an in vivo micro-CT system becomes available for use in these projects. 126 Chapter 4 detailed the use of an infant porcine model for pediatric skull fracture to study the effects of fall height and interface compliance on fracture patterns. It was shown that an increased fall height resulted in significantly higher impact forces, longer impact durations, and more cranial damage as indicated by total fracture length. Using diagrammatic representations of the resultant fracture patterns, the study additionally showed that high falls dramatically increased the incidence of fracture in cranial bones adjacent to that which was impacted and resulted in substantial damage to the side of the skull opposite the impact. Using different carpets for impact interfaces, this body of work documented that increasing interface compliance significantly reduced the degree of cranial bone fracture. Following impact onto a highly compliant carpet interface, the cranial fracture pattern matched that due to an impact from a shorter fall height onto a rigid surface. The study showed that in order to accurately determine the cause of a specific cranial fracture pattern, there are a number of variables that must be considered including the alleged height of the fall and characteristics of the potential impacted surface. Using these and other data from the infant porcine model, future work should focus on developing machine learning techniques to help analyze fracture patterns and automatically classify them according to the impact scenario. This tool could then be adapted to human case studies for the development of a tool to help medical examiners identify the cause of a trauma using quantitative versus current qualitative techniques. Chapter 5 focused on fracture patterns developed on long bones of the appendage as a result of applied bending loads. This thesis documented that there were characteristic 127 features of the resultant fractures that consistently identified the direction of the applied lateral load. In addition, a more phenomenological mechanism of fracture development was presented as adapted from the field of dynamic fracture mechanics rather than quasi-static beam theory to account for the location and characteristic pattern of long bone fracture. A suggestion for future work may be to study more laboratory-controlled impact orientations and loading scenarios on long bones, such as pure bending and varying axial load components, to help verify characteristic fracture patterns that might prove helpful to the forensics community for the determination of injury causation. 128