Building Finite Element Models to Investigate Zebrafish Jaw Biomechanics.

Skeletal morphogenesis occurs through tightly regulated cell behaviors during development; many cell types alter their behavior in response to mechanical strain. Skeletal joints are subjected to dynamic mechanical loading. Finite element analysis (FEA) is a computational method, frequently used in engineering that can predict how a material or structure will respond to mechanical input. By dividing a whole system (in this case the zebrafish jaw skeleton) into a mesh of smaller 'finite elements', FEA can be used to calculate the mechanical response of the structure to external loads. The results can be visualized in many ways including as a 'heat map' showing the position of maximum and minimum principal strains (a positive principal strain indicates tension while a negative indicates compression. The maximum and minimum refer the largest and smallest strain). These can be used to identify which regions of the jaw and therefore which cells are likely to be under particularly high tensional or compressional loads during jaw movement and can therefore be used to identify relationships between mechanical strain and cell behavior. This protocol describes the steps to generate Finite Element models from confocal image data on the musculoskeletal system, using the zebrafish lower jaw as a practical example. The protocol leads the reader through a series of steps: 1) staining of the musculoskeletal components, 2) imaging the musculoskeletal components, 3) building a 3 dimensional (3D) surface, 4) generating a mesh of Finite Elements, 5) solving the FEA and finally 6) validating the results by comparison to real displacements seen in movements of the fish jaw.


Introduction
Finite Element (FE) modelling is an engineering technique that can computationally calculate and map the magnitude and location of strains acting on a structure 1 . The model consists of the 3D structure, represented by a mesh of "Finite Elements", and the end result of the analysis is governed by a number of factors including the structure and number of elements in the mesh, the magnitude and location of the mechanical loads and the material properties. Material properties describe certain aspects of a material's behavior under a given type of load; Young's modulus (E) describes the elasticity of the material while Poisson's ratio describes the proportional decrease in the width of a material to its length when a sample is stretched. FE modelling can be used to calculate a variety of variables including displacement, stress, pressure and strain acting on the model by taking into account the unique input data about the structure's shape, location and magnitude of loads and the specific material properties.
FE modeling is widely used in engineering 2 and increasingly for orthopedic 3 and paleontological applications 4 . In development biomechanical forces are known to act as a stimulus in many cells to activate cell responses [5][6][7][8] and it is useful to predict both the relative positions and magnitudes of mechanical stimuli within developing organ systems, however, currently FE modeling has been little used for zebrafish development.
Both cartilage and bone have been shown to be mechanosensitive materials. For example, in vitro compression has been found to activate chondrogenic pathways, whilst tension has been shown to be necessary for bone formation 9 . FE analysis (FEA) has been exploited to model strains acting on biological specimens, including those acting on skeletal elements during bone formation 10 . Other development applications include its use to predict the shape of a joint after it has been exposed to theoretical biomechanical forces 11,12 and to show the pattern of strains present during chick knee joint morphogenesis 8 .
This protocol is aimed at sharing the experience of generating 3-dimensional surfaces, meshes and Finite Element models from confocal images with a view to understanding the mechanics of developing tissues. We also show ways of validating the FE models though capturing real joint displacement information in vivo. While we use the zebrafish jaw as an exemplar the same techniques could be used on any small biological system for which 3D information on the structure of the musculoskeletal system can be obtained by confocal or multiphoton imaging. 3. Permeabilize larva in 0.25% trypsin in PBT on ice for 5-6 min. Wash in 4x PBT for 5 min each. 4. Block for 2-3 hr in 5% serum in PBT. 5. Incubate larva in the recommended dilution of rabbit anti-type 2 collagen and mouse anti-myosin antibodies in 5% serum in PBT for 1 hr at room temp or overnight at 4 °C. NOTE: The recommended dilution range is normally on the antibody data sheet. Choose antibodies that are raised against different species to one another and that are also different to the tissue. 6. Wash larva 6x for 15 mins in PBT. 7. Block for 1-2 hr in 5% serum in PBT. 8. Incubate in secondary antibodies in the dark. Use fluorescently labelled anti-mouse (550) and anti-rabbit (488) secondary antibodies at an appropriate dilution in 5% serum and PBT for the specific antibody. 9. Wash 6X for 10 minutes each in PBT and image on a 10X confocal microscope as soon as possible.

Imaging Musculoskeletal Geometry
1. Mount the larvae ventrally on a coverslip in tepid 0.3-0.5% low melting point (LMP) agarose in Danieau's solution 16 . NOTE: Transgenic fish will need to be sedated in 0.02% MS222 (Tricaine methanesulfonate, pH 7) before mounting and during imaging. 2. Take a confocal image stack of the region of interest using the 10X objective lens and the approximately 2.5X digital zoom. Prepare images of the green and red channel using the 488 nm and 561 nm lasers respectively. Image at a 512 x 512 pixel resolution with an interval between z planes of 1.3 µm and 3 line averages. The resulting stack will comprise of approximately 100 z slices. 3. Export the data as a tiff series. Maximum projections of the muscle and cartilage elements from a 5dpf zebrafish larva are shown in i.e., cartilage and joint. Select the cartilage region of the image ( Figure 2C, white signal, purple outline) using the magic wand tool. Use the brush tool to remove noise from the outlines. NOTE: If using magic wand tool, click 'All slices'. 4. Select the joint region with the brush tool and assign to a joint component ( Figure 2C, blue outline) 5. Smooth multiple slices at once by selecting segmentation in the top menu and smooth labels. Right click on the image and select generate surface to produce a 3D surface render of the component ( Figure 2D). 6. Click on the surface and save data as a hmascii file for import into meshing software.

Calculating the Muscle Forces to Be Used in the FE Model
1. Count the number of muscle fibers from confocal images of smyhc:GFP transgenic zebrafish ( Figure 1A, arrowhead, 1C) and measure the diameter of the fibers to calculate their cross sectional area (πr 2 ). 2. Identify suitable force per muscle unit area from the literature. The maximal muscle force generated per unit area for larval zebrafish skeletal muscles (40 nN/µm 2 ) were used 17 . 3. Calculate the forces for each anatomical muscle group by multiplying the number of fibers and their area by the force per unit area. See Table 1.

Generating a Mesh
1. Import the 3D model generated in section 2 (above) into a software package capable of generating a finite element mesh. 2. Generate a2D mesh of the cartilage and joint surfaces by using the shrink wrap tool under the 2D menu. Choose an appropriate element size. Note: Use an element size between 1.5-2.5. If necessary, generate a range of differently sized 2D surface meshes to carry out 3D mesh optimization (Section 4.4). 3. Carry out mesh quality checks found under the '2D>Tools>Check elements' panel to check for duplicated elements, insertions and penetrations in the mesh. Fix dihedral angles by using the utility tab in model tree. 4. Generate a 3D mesh from the 2D surface meshes of differing element size using the 3D>Tetramesh subpanel.
NOTE: Compare the results of different mesh sizes and select the FE model with the lowest mesh size that converges after further simulations and does not compromise feature definition. The example in Figure 3 contains 1.5 million tetrahedral elements for the lower jaw cartilages and had a 2D element size of 2.0. 5. Transform the mesh so jaw model is to scale as per confocal stack using the Geometry>Distance subpanel.
NOTE: Ensure the cartilage and joint components are connected in the model by exporting a merged model or by using ties.

Finite Element Model Construction
1. Using commercial finite element (FE) software, create a FE model. Using the 3D muscle and cartilage labeled confocal stacks generated in section 1 as a guide, assign nodes that correspond to muscle attachment points. Create a vector between two nodes representing the origin and insertion of each muscle (Figure 3). 2. Create a load collector of type 'history' to apply a 'C load' for each muscle. Specify the magnitude in Newtons (calculated in step 3.3) and assign the associated vector. Figure 3 shows the attachment points for the adductor mandibulae (AM), protractor hyoideus (PH) and intermandibularis (IM). NOTE: For these jaw muscles, maximum contractile force is distributed between the origin and insertion so that only 50% of each load is applied at each site. 3. Assign appropriate elastic isotropic material properties as determined by the literature. Young's modulus for cartilage and the interzone in this model were 1.1 MPa and 0.25 MPa respectively and Poisson's ratio was 0.25 for both 18,19 . 4. Create a load collector of type 'boundary' to apply initial constraints on the model. Go to the tab Analysis>Constraints and in the create subpanel, pick the nodes on the model you wish to constrain. Select the degrees of freedom (DOF) that limit movement of the model to the best approximation of its natural range of motion. NOTE: The model in Figure 3 was constrained in all axes of motion (DOF: 1, 2, 3 represent x, y and z, respectively) at the ceratohyal to anchor it in space at a mid-point in the model and in the y and z axis at the point where the palatoquadrate connects to the rest of the zebrafish skull (Figure 3, Table 1). The model must be constrained in all three DOF in a least one node. 5. Create a 'Load step', for each type of movement you wish to simulate (i.e., opening, closing), under the analysis menu and select all the relevant loads (made in Section 5.2) and constraints (made in Section 5.4) to simulate this movement. Select 'Static' from the drop down menu when it appears. 6. Export model including the mesh, loads, constraints and material properties in an appropriate file format, in this case ".inp" format. 7. Load model into FE analysis software. Create and execute a job for the model using the Job module. 8. Analyze output for stress, strain, displacement, etc. found in the results tab and visualization menu (Figure 4 and 5). 6. Validation of Jaw Deformation/Displacement Distances 1. Select 3-6 Tg(Col2a1aBAC:mcherry) transgenic zebrafish. 2. Lightly anaesthetize the larvae with 0.02% MS222 until they cease to respond to touch but their hearts are still beating. 3. Mount the larvae laterally (while anaesthetized) on coverslips in tepid 1% LMP agarose (made up in Danieau's solution). 4. Remove the agarose from around the head and jaw with forceps. 5. Flush fresh Danieau's solution (with no MS222) over the head of the larva to remove anesthesia using a Pasteur pipette until normal mouth movements resume. 6. Use movie capture software to take bright-field high-speed videos of mouth movements. Take movies of around a minute duration at the highest frame rate, or sufficient to record multiple cycles of jaw opening. 7. Choose frames that show the jaw open to its maximum displacement. Measure the distance between the anterior tip of the Meckel's cartilage and the upper jaw (tip of ethmoid plate) in µm. 8. Calculate the average displacement from multiple fish larvae. 9. Extract displacement data from the model. Use the average displacement calculated in 6.8 to verify model displacement behavior (Figure 4).
If the model deviates from expected, run sensitivity analysis by sequentially altering material properties and muscle loads.
The FE models once run can be used to display a range of data, such as stress (Figure 5A), minimum and maximum principal strain ( Figure  5B-K). These results are three-dimensional so the model can be magnified to see fine patterns of detail (Figure 5E, 5I) rotated to obtain relevant views (Figure 5F, 5G, 5J, 5K) and digitally sectioned (Figure 5E', 5E'', 5I', 5I'') to show how the patterns of stress, strain or pressure change throughout the model. It is also possible to extract quantitative data from the model (not shown). By verifying the model and using the most accurate material properties, loads and mesh shape the FE model can be used to explore the best estimate of the mechanical environment being experienced by cells during that window of development. The results of the model can be directly compared to changes in cellular behavior and gene expression 20 .

Discussion
Finite Element models have been used to relate the areas of skeletal elements that are under strain with those undergoing bone formation 10 , as well as to map the areas under strain during endochondral ossification and joint morphogenesis 8,12,21 . Other studies have also been able to apply theoretical growth models to replicate changes during joint development 11,12 . Here we show the protocol for building FE models for a relatively simple system, the zebrafish jaw 20 . Unlike alternative methods of collecting raw images for the FE models, such as CT scanning 22 , confocal imaging of transgenic lines or immunostained zebrafish allows for multiple tissues to be studied. It can, therefore, provide direct information on muscle attachment points in relation to cartilage. Among vertebrate models zebrafish are particularly amenable to genetic and pharmacological manipulation. The generation of FE models for zebrafish craniofacial cartilage now opens up the possibility of further study of the interplay between biomechanics and genetics in joint morphogenesis.
There are a number of critical steps to the process of creating an FE model; the first is generating an accurate three-dimensional representation of the system. This requires imaging at high enough resolution to clearly define boundaries. Note that even with high-resolution imaging to make a good surface one may have to smooth out some regions. Another critical step is defining the correct placement of the load and correct constraints. An insufficiently constrained model will fail to solve and incorrect placement of the loads will cause abnormal movement.
Some processing of the raw data (Figure 2) is necessary as a surface generated from the raw data would be difficult to mesh (Figure 2B). We filtered the data using a Gaussian filter ( Figure 2C) and we carried out some manual smoothing of the curves to produce a set of clean outlines that can be converted into a 3D surface. Too much smoothing can produce a "melted" surface that has lost many of its features. Choosing the correct element size is an iterative process as choosing too small an element size creates too large a mesh which is computationally intensive. However, choosing too large an element size will produce a mesh which fails to recapitulate the correct shape of the structure. The correct mesh had the smallest element size that captured the correct shape of the jaw and converged on a correct solution, verified using the jaw displacement. It may also be necessary to modify the material properties or load calculations to better emulate the correct displacement as different ages and species will have substantially different properties.
It is important to remember that there are always limitations to a hypothetical model and assumptions made to run FE models. When only modeling one or a small number of samples it is critical to ensure that a representative sample is chosen as there are likely to be small variations between individuals. As only some of the jaw elements and muscles were included, the model is a simplified version of the zebrafish craniofacial musculoskeletal system. Therefore, constraints had to be positioned to account for where the modelled jaw elements would connect with the rest of the skull and the model was artificially constrained in the center to fix it in 'space'. This artificial constraint did not impact on the interpretation drawn from the models as the ceratohyal itself was not analyzed. The inclusion of more of the craniofacial structure, especially other jaw opening muscles such as the sternohyals and its attached cartilage 23 , could have added to the model, but limitations include the ability of larger models to run in the Finite Element software.
Another limitation is that we have not modeled ligament insertion, though this could be achieved by the insertion of springs 8 . One other assumption made in this case was that the model would behave linearly. The magnitudes of strains on the models were comparable to those in published models and applied to in vitro cells 10,24 , with strains being below +3,500 and above -5,000 µɛ apart from constraint and muscle attachment points. Therefore, the strains at the relevant regions of the model were deemed within a range acceptable for a linear model. Cartilage does not behave entirely as a linear material and has previously been modelled as a poroelastic material, which enabled analysis of the fluid behavior in the model 25 . Spreading the muscle attachment points amongst a cluster of local nodes would distribute the peak forces and more accurately represent the muscle insertion for certain muscles.
Use of FE allows an assessment of the strains and stresses acting on a structure. As a technique it is frequently used in many bioscience disciplines including orthopedics, paleontology and more recently developmental biology. Here we describe how to build FEs for the zebrafish lower jaw. In the future these models could be extended to look at the whole jaw, including the palate. Similar techniques could be used to model spinal biomechanics in fish, which to date have mostly been studied by kinematic means.

Disclosures
The authors have nothing to disclose.