Finite element and deformation analyses predict pattern of bone failure in loaded zebrafish spines

The spine is the central skeletal support structure in vertebrates consisting of repeated units of bone, the vertebrae, separated by intervertebral discs (IVDs) that enable the movement of the spine. Spinal pathologies such as idiopathic back pain, vertebral compression fractures and IVD failure affect millions of people worldwide. Animal models can help us to understand the disease process, and zebrafish are increasingly used as they are highly genetically tractable, their spines are axially loaded like humans, and they show similar pathologies to humans during ageing. However, biomechanical models for the zebrafish are largely lacking. Here, we describe the results of loading intact zebrafish spinal motion segments on a material testing stage within a micro-computed tomography machine. We show that vertebrae and their arches show predictable patterns of deformation prior to their ultimate failure, in a pattern dependent on their position within the segment. We further show using geometric morphometrics which regions of the vertebra deform the most during loading, and that finite-element models of the trunk subjected reflect the real patterns of deformation and strain seen during loading and can therefore be used as a predictive model for biomechanical performance.


Introduction
The spine consists of a repeated pattern of motion segments (MSs) of bony vertebrae separated by intervertebral discs (IVDs) that enable movement. Back pain and IVD degeneration affect millions of people worldwide [1,2], and vertebral compression fractures are a frequent feature of osteoporosis [3]. Biomechanical pathologies of the spine are underpinned by genetic, physiological and environmental pathways that together damage IVD, muscle and the bone, changing the mechanics of the system. Animal models, typically rodents, are frequently used to study mechanisms of spinal pathology [4]. However, quadrupeds are disadvantageous for studying the human spine as gravitational load acts perpendicular to their axial skeleton. Zebrafish are increasingly used as a model for human disease, due to their genetic tractability. Unlike quadrupeds, but similar to humans under gravity (figure 1a), their spine is antero-posteriorly loaded as a result of swimming through viscous water [5]. Zebrafish are well established as models for skeletogenesis, pathology and ageing [6], and develop spinal pathologies in response to altered genetics [7] and ageing [8]. However, the biomechanics of the zebrafish spine are comparatively poorly characterized.
Finite-element analysis (FEA) has proven a pivotal tool in the study of biomechanical subjects [9], and offers a method for biomechanically characterizing the zebrafish spine, including intact MSs. This technique digitally models an object of known material properties using a series of linked nodes of known    royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 16: 20190430 number and geometry, that can be subjected to a wide variety of forces outputting the predicted geometry, strain and deformation. Results can be validated by comparison with the results of loading experiments in which a sample is loaded ex vivo [10,11]. FEA has been used in zebrafish to test contributions of shape and material properties in joint morphogenesis [12,13] and to study strain patterns in a single vertebra [14].
Here, we describe a novel integrated experimental platform that brings together imaging, modelling and real-world validation to explore the biomechanics of intact zebrafish spinal MSs. We generated an FEA model of the spine, which we validated with a loading experiment using a high-precision material testing stage (MTS) under set loading regimes using micro-computed tomography (µCT). Three-dimensional geometric morphometrics (3D-GM) was used to explore patterns of deformation seen in each vertebra during loading. Comparison of results demonstrated that our FEA model accurately predicted the relative patterns of deformation and strain experienced by real samples loaded ex vivo.

Zebrafish samples
One-year-old, wild-type (WT) zebrafish were fixed in 4% paraformaldehyde and dehydrated to 70% EtOH. MSs were acquired by making two cuts in the trunk, between the morphologically homogeneous vertebrae 18 and 24 of a total of 33 vertebrae [5] (figure 1a-c).

In vitro vertebral loading experiment
Loading experiments were conducted using a custom-built material testing stage (MTS2) in a Bruker SKYSCAN 1272 µCT system. Radiographic visualization of each MS (n = 3) was performed and if required, vertebrae were trimmed to retain three complete vertebrae and associated IVDs (figure 1b-d). Samples were stabilized (anterior-up) in the MTS2 using cyanoacrylate glue. The MTS2 was programmed to perform a sequential series of seven scans at a series of increasing loads (table 1), using 60 keV X-ray energy, 50 W current, 5 µm isotropic voxel size and a 0.25 mm aluminium filter. A total of 1501 projections were collected during a 180°rotation, with 400 ms exposure time. Reconstructions were performed using NRecon (v. 1.7.1.0). Surfaces of vertebrae, muscle and IVDs in each dataset were generated using Avizo (Avizo v. 8; Vizualisation Sciences Group) (figure 1c-e and table 1) and linear measurements of IVDs and MS lengths made using the '3D Measurement' tool. Vertebrae surfaces were further processed in Meshlab (table 2).

Finite-element analysis
An MS surface mesh was created based on a 1-year-old WT specimen µCT scanned using a Nikon XTH 225ST μCT system as described under two conditions: (a) native state and (b) contrastenhanced following 14 day incubation in 2.5% phosphomolybdemic acid [16]. Scan (a) was used to segment vertebrae (V18-V24), and scan (b) to segment IVDs. The resulting binary labels from scans (a) and (b) were saved as 8-bit tiff stacks, manually registered in 3D space in Avizo ('Trackball' tool) and algorithmically combined ('Algebra' tool), creating a single volume of separate materials representing three vertebrae and four IVDs (figure 1d,e and table 2). A 500 µm thick cylinder was created contacting the anterior-most IVD perpendicular to the model axis, to mimic the stainless-steel compressive plate and distribution of forces applied during loading (figure 1f ).
The complete vertebral surface mesh was imported into Simpleware ScanIP (v. 2018.12, Synopsys Inc.) to create an FE model. The model consisted of 1 054 187 linear tetrahedral elements joined at 257 392 nodes comprising four material types: vertebral bone, annulus-fibrosus, nucleus-pulposus and stainless steel (figure 1d-f, table 2). The model was analysed in Abaqus (2018 version). A custom datum coordinate system was created centred on the antero-posterior axis of the model, and a concentrated force applied to the central node of the anterior face of the compressive plate. This loading case was repeated in each of seven steps of a multi-step analysis, with load values matching the increments applied in the MTS (table 1). The model was constrained in two locations using boundary conditions, at the base of the posterior-most IVD (constrained in three axes) and at the top of the compressive plate (constrained in two axes), allowing movement along the model's antero-posterior axis (figure 1f ). Deformed meshes from each step were exported as surface files and analysed using 3D-GM for quantitative comparison between relative and absolute patterns of deformation predicted by FEA and observed in MTS data.

Three-dimensional geometric morphometrics
Three-dimensional geometric morphometrics analysis of vertebral deformation was performed using the 'Geomorph' package for the R statistics software [17]. For each loading experiment, we used the first scan (1 N load) to create a template of 3D coordinates for 22 fixed three-dimensional landmarks (figure 2a-c) linked by 300 surface sliding semi-landmarks (using the 'buildtemplate' function). By assigning the same landmarks in each scan (using the 'digitsurface' function), we compared the first scan with subsequent scans of the same vertebra using generalized Procrustes analysis (allowing semilandmarks to 'slide' in order to remove arbitrary spacing). Resulting shape variables were subjected to principal component analysis (PCA) to identify the principal patterns of variation between scans of the same vertebra, and isolate trends in deformation with increasing compressive load.

Vertebral motion segments fail under loading of 12-16 N at positions of maximum von Mises strain
To test the range of compressive loads that the MS could resist until failure, we subjected an MS to exponentially increasing compressive forces from 1 to 100 N. This specimen failed at 16 N whereupon the central vertebra fractured midcentrum. A primary loading regime between 1 and 16 N was thus established (table 1) for the three primary specimens; occupying the elastic, plastic and failure regions of the compressive loading profile of a typical MS. Failure was considered when at least one vertebral centrum fractured across the axis (e.g. figure 1j,l ). All samples failed between 12 and 16 N upon shallow angle fracture in the central vertebra, with the smallest specimen (specimen 3) failing at the lowest force (figure 1g,h). This is higher than maximum aquatic forces experienced during swim training by Fiaz et al. [5], which reached approximately 9.5 N. Minor differences in mounting orientation created differences in linear deformation between right and left sides, but specimens follow similar patterns. Prior to failure, linear measurements show an increase in IVD antero-posterior thickness (        We found characteristic patterns of deformation and strain in response to compressive loading of zebrafish vertebrae. Three-dimensional results from MTS data follow distinct trends for each vertebrae between the three specimens (figure 2d,i,n), showing consistent dorsoventral compression, and lateral compression that is reversed at higher loads potentially due to elastic rebound of the IVD and fracturing along the zygopophyses that occurs at these loads ( figure 2). This relative pattern is shared between each specimen, although specimen 3 experiences this at lower loads than specimens 1-2, before failing at 12 N. Fractures are observed where the arches and zygopophyses contact the centrum, at loads that precede the failure of the segment (figure 2f,h,k,m,p,r). Comparison with FEA data (blue points in figure 2d,i,n) suggests that the FE model accurately predicts these patterns (figure 2d,i,n), and that patterns of deformation could explain the first signs of damage prior to failure. In both datasets, the anterior vertebra undergoes most deformation, particularly posterior deformation of the arches (figure 2e-h). The central vertebrae and arches show strong torsion (figure 2j-m), increasing through the loading regime leading to the failure of the segment ( figure 1l,o). The posterior vertebra shows the least deformation and is most isotropic in pattern (figure 2or), potentially due to protection offered by the anterior IVDs. Comparison with ex vivo loading of vertebral MSs validates the accuracy of our FEA model for predicting patterns of deformation and strain across these structures. This offers a step towards a digital 'sandbox' approach to modelling the effects of genetic, physiological and morphological properties on the reaction and resistance of vertebral MSs to loading. Inputting specific properties of vertebral samples into a validated FE model will allow their effects on the biomechanics of the spine to be quantitatively tested in silico, allowing the relative contributions of shape and material properties to be explored and empirically tested. This will aid comparison of mechanical performance between different model systems. As an advantage of the zebrafish system is the wealth of mutants modelling human disease genetics [18], comparisons of mechanical performance between genotype and phenotype will be possible. In the longer term, this approach may give insight into biomechanical aspects of spinal pathology, allowing identification of 'at risk sites' in the spine. This could provide a basis for more specific or earlier interventions than those commonly employed.