Finite element modelling predicts changes in joint shape and cell behaviour due to loss of muscle strain in jaw development

Abnormal joint morphogenesis is linked to clinical conditions such as Developmental Dysplasia of the Hip (DDH) and to osteoarthritis (OA). Muscle activity is known to be important during the developmental process of joint morphogenesis. However, less is known about how this mechanical stimulus affects the behaviour of joint cells to generate altered morphology. Using zebrafish, in which we can image all joint musculoskeletal tissues at high resolution, we show that removal of muscle activity through anaesthetisation or genetic manipulation causes a change to the shape of the joint between the Meckel's cartilage and Palatoquadrate (the jaw joint), such that the joint develops asymmetrically leading to an overlap of the cartilage elements on the medial side which inhibits normal joint function. We identify the time during which muscle activity is critical to produce a normal joint. Using Finite Element Analysis (FEA), to model the strains exerted by muscle on the skeletal elements, we identify that minimum principal strains are located at the medial region of the joint and interzone during mouth opening. Then, by studying the cells immediately proximal to the joint, we demonstrate that biomechanical strain regulates cell orientation within the developing joint, such that when muscle-induced strain is removed, cells on the medial side of the joint notably change their orientation. Together, these data show that biomechanical forces are required to establish symmetry in the joint during development.


Introduction
The development of reciprocal interlocking joints at cartilage elements is central to ensuring normal skeletal function (Nowlan et al., 2010). Processes that disrupt joint shape formation can cause abnormal loading and joint function (Michaeli et al., 1997), e.g. hip shape correlates strongly with risk of osteoarthritis (Jacobsen and Sonne-Holm, 2005). The initial formation of the cartilage template from mesenchymal cell condensations, mostly replaced by bone, is well understood (DeLise et al., 2000;Thorogood and Hinchliffe, 1975). However, how the early joint structures undergo morphogenesis to form their mature shape remains less clear (Pacifici et al., 2005).
Zebrafish, with relevant fluorescent transgenic lines (Apschner et al., 2014;Hammond and Moro, 2012), allow dynamic imaging of the musculoskeletal system at cellular resolution. Zebrafish are, therefore, a useful model to examine how mechanical loading from muscle impacts on cartilage behaviour. By 48 h post-fertilisation (hpf), mesenchymal cells have condensed to form the mandibular arches (Eames et al., 2013;Kimmel et al., 1995). At 53-55 hpf, the cartilaginous elements of the Meckel's cartilage (MC), palatoquadrate (PQ) and ceratohyal appear, as does the adductor mandibulae jaw musculature, with the intermandibularis anterior and protractor hyoideus, identifiable by 62 hpf (Schilling and Kimmel, 1997). Larval zebrafish use the protractor hyoideus to constrict the buccal chamber of the jaw and adductor mandibulae to close the mouth (Diogo et al., 2008;Hernandez et al., 2002). The joint between the MC and PQ (Figs. 1A and 2A) is described as the jaw joint in (Talbot et al., 2010) and referred to as such hereafter. In the joint structure, the retroarticular process (RAP) of the MC protrudes ventrally to interlock with the PQ (Miller et al., 2003), typically cavitation of this joint occurs at around 7 dpf (http://zfatlas.psu.edu/).
Many studies have linked absence of muscle activity with abnormal joint shaping and fusions of articular surfaces in long bones. Early studies used Decamethonium Bromide and botulinum toxin to generate paralysis in chick models; leading to a flattening of articular surfaces and a failure of joint cavitation (Drachman and Sokoloff, 1966;Murray and Drachman, 1969). More recently, (Roddy et al., 2011b) found that rigid paralysis of chicks during early development caused the knee joint to appear flattened.
Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com Additionally, genetic models in mice, such as Splotch mutants that lack limb muscle, exhibit abnormal limb joint shaping (Kahn et al., 2009;Nowlan et al., 2010). While muscle force has been shown to be necessary for normal joint morphogenesis and chondrocyte intercalation (Shwartz et al., 2012), it remains largely unclear how cells within the joint interpret such forces to bring about changes in behaviour. Little is known concerning the effects of paralysis on jaw joint morphology, though the temporomandibular joint region shows signs of adaptation when the mechanical environment is altered (Enomoto et al., 2014;McNamara and Carlson, 1979). Finite element (FE) models simulating the biophysical response to muscle-induced loading have been used to investigate joint development, endochondral ossification, and joint development, including developmental dysplasia of the hip (DDH) (Carter and Wong, 1988;Heegaard et al., 1999;Nowlan et al., 2008;Roddy et al., 2011a;Shefelbine and Carter, 2004). Thus far, developmental FE-models have focused primarily on the femur, particularly in humans and chick. Whilst there are a handful of studies using FE to deduce the mechanical performance of extant shark jaws (Ferrara et al., 2011) and other jawed vertebrates (Rayfield, 2007), FE-modelling has not been used to explore the mechanobiology of the developing zebrafish jaw pre-cavitation.
Here, we document the process of joint morphogenesis in wild type zebrafish jaws from the time of first muscle activity to generation of the refined interlocking joint shape. Then, by removing muscle activity pharmacologically and genetically, we quantify the timing and extent of the response to lack of muscle activity on joint morphogenesis. Using FE analysis (FEA) we identify the locations of muscle-induced strain acting on the zebrafish jaw cartilage throughout normal joint morphogenesis. To understand the mechanobiological changes that underpin joint shape we quantify differences in cellular orientation between wild type, anaesthetised and mutant models.

Zebrafish husbandry/zebrafish lines
Zebrafish were maintained as described (Westerfield, 2000). All zebrafish experiments were approved by the local ethics committee and the Home Office (Project license number 30/2863). The Tg(Col2a1aBAC:mcherry) zebrafish line has been previously described (Hammond and Moro, 2012;Hammond and Schulte-Merker, 2009) and labels all chondrocytes ( Fig. 2A). The line Tg(smhyc:EGFP)i104 labels all slow twitch muscle fibres (Elworthy et al., 2008). myod fh261 mutants have   been previously described and lack craniofacial muscle, with the exception of the sternohyoideus (Hinits et al., 2011).

Counting mouth movements and measuring jaw displacement
Zebrafish were anaesthetised with MS222 (Tricaine methanesulfonate), (Sigma) and mounted laterally onto coverslips in 3% agarose. The agarose surrounding the head was removed and Danieau buffer (Westerfield, 2000) flushed over the coverslip until jaw movements resumed. The number of mouth movements in 1 min was recorded from at least six fish per timepoint. Movies of jaw movements were made using fish labelled with Tg(Col2a1aBAC:mcherry) on a Leica SP8 confocal microscope. Jaw movement was quantified as the positional change at the anterior end of the MC measured (mm) from individual movie frames.

Jaw imaging
Confocal images of approximately 100 slices at 1.3 mm intervals were taken of jaws at 3, 4 and 5 dpf labelled with Tg(Col2a1aBAC:mcherry) using an SP8 Leica confocal.

Outline traces
The curve drawing tool in powerpoint was used to produce the outlines of the jaw joint labelled with Tg(Col2a1aBAC:mcherry) from stacks of confocal images; each specimen was drawn in a separate colour.

Interzone measurements
The interval between the MC and PQ cartilage elements of the jaw joint on their medial and lateral sides (typically the smallest and largest gaps between cartilage across the joint, respectively) were measured (from both the left and right joints) from confocal images of the jaw using Leica LAS AF Lite software. Negative values were recorded for overlapped cartilage elements.

Geometry and Finite Element Analysis
A 3D 512 Â 512 pixel resolution representation of the stained larval zebrafish jaw at 3, 4 and 5 dpf was produced using confocal microscopy. These datasets were imported into Avizo (Version 7.0.0 FEI Visualization Science Group) for material segmentation and digital reconstruction to produce 3D models of the cartilage elements. The models were imported into Hypermesh (Version 10, Altair Engineering) to generate a solid mesh of approximately 1.5 million tetrahedral linear elements (Supp. Fig. 1, Supp. Table 1). Due to the absence of known zebrafish cartilage material properties, a range of Young's moduli for the cartilage and interzone material were tested on the 5 dpf model to determine which properties generated jaw displacements, (measured from the tip of the MC, Fig. 1C), that best matched those observed in Table 1. When 50% of maximum calculated muscle force was applied to the models, all models that had a cartilage Young's modulus of 1.1 MPa fell within a physiological range of jaw displacement (Table 1). Young's modulus of 1.1 MPa for cartilage (as previously reported for unmineralised embryonic murine cartilage (Tanck et al., 2004)) and 0.75 MPa for the interzone were chosen for the models (Fig.1). These material properties were assigned hereafter, with a Poisson's ratio of 0.25, previously described for embryonic unmineralised murine cartilage (Tanck et al., 2000).

Muscle force measurement and calculation
To model strains acting on the jaw, we calculated muscle forces during jaw opening (protractor hyoideus and intermandibularis) and closing movements (adductor mandibulae). The Tg(smhyc:EGFP)i104 (slow twitch skeletal muscle) transgenic line was used to determine muscle attachment points. However, as jaw muscles are composed of a mixture of fast and slow twitch fibres (Hernandez et al., 2005), we performed immunohistochemistry for the antibody A4.1025, which detects both slow and fast twitch fibres (Dangoor et al., 1990), allowing us to count the actual number of muscle fibres in each muscle (Table 2).
Cross-sectional area of each fibre was used to determine the force that each muscle group exerts at 40 nN/mm 2 (Table 2), the force for skeletal larval zebrafish muscle previously described (Iorga et al., 2011). Muscle attachments were added to the FE-models for mouth closure (adductor mandibulae) and opening (protractor hyoideus and intermandibularis) using confocal datasets for reference for attachment sites (Supp. Figs. 1, 1A and B). Models were constrained from movement at the junction between the hyosymplectic and palatoquadrate where the jaw attaches the skull and the ceratohyal (Supp. Table 1, Supp. Fig. 1). The left and right side of the models were constrained in y and z, whilst the ceratohyal was constrained in x, y and z. The geometrically linear models were imported from Hypermesh into Abaqus FE-software (v6.10.2 Simulia, Dassault Systèmes) for analysis. Maximum principal strain, (tension), and minimum principal strain (compression) were recorded.

Cell orientation measurements
Cell orientation analysis on the Meckel's cartilage joint was performed using confocal images of joints labelled with II-II6B3 antibody (expanded methodology in Supp. Fig. 2). Cells from three joints per condition were identified by thresholding in ImageJ (Supp. Table 2, Supp. Fig. 2B). The angle of the joint from the horizontal axis was recorded (Supp. Fig. 2B). An automated readout of cell orientation was produced (Supp. Fig. 2C). These data were exported into Excel, and angles corrected in relation to the joint axis, (Supp. Fig. 2C). Joint cell orientation was plotted in circular histograms using PAST software (Hammer et al., 2001).

Statistics
Statistics were performed using SPSS software. Kruskal-Wallis tests were performed to compare the size of the cartilage interval at the medial and lateral edges of the joint following each anaesthesia treatment (Figs. 2F and 3C). The Kruskal-Wallis test was used to make multi-comparisons between non-normal data.

Jaw mobility
The first stage at which jaw opening was reliably present was 72 hpf (Fig. 2B). Subsequently, there was a significant increase in jaw motility between 72 and 96 hpf, whereas there was no significant difference between 96 hpf and 120 hpf. Mouth opening had therefore reached its maximal frequency by 96 hpf. Table 1 The maximum jaw displacement at 5 dpf in 1 min. The average jaw displacement is 37.2 mm.

Jaw morphology
Confocal imaging and the resulting 3D digital reconstructions show that the interval between cartilage elements (the joint interzone) changes in morphology across the time period in the control fish (Fig. 2). At 3 dpf there is a smaller interval between the Meckel's cartilage (MC) and the palatoquadrate (PQ) on the medial side than on the opposing lateral side (Fig. 2C, C'), i.e. the joint interzone is asymmetrically shaped at 3 dpf. By 5dpf the interzone is similar in width across its extent (Fig. 2 C-E'). To quantify this change in morphology we measured the interval between the MC and the PQ at the medial and lateral side at each time point (Fig. 2F). Kruskal-Wallis tests showed a significant difference in the relative size of the interzone at its greatest and smallest extent at 3 dpf, which decreases over time, such that by 5 dpf the difference in the size of the interval across the joint is no longer significant and the interval is roughly uniform across its extent. myod mutants, (which lacks all jaw movement), have a similar morphology to wild type at 3 dpf (Supp. Fig. 4) but a significant difference between the cartilage element interval on the medial and lateral side at 5 dpf compared to 5 dpf wild type (Fig. 2F, G and G'). We also tested the size of the cartilage interval of the myod compared to younger wild type fish; myod fish at 5 dpf retain a joint interzone shape comparable with that of a 3 day old wild type ( Fig. 2F and G, Supp. Table 3).
Fish anaesthetised from 3-5 dpf have a significant difference in the extent of the cartilage interval between the medial and lateral side compared to wild type fish at 5 dpf (Fig. 2F, H, and H'), but are not significantly different to the myod mutants ( Fig. 2F-H, Supp. Table 3). Therefore, it is the activity rather than the presence of muscle which leads to joint morphogenesis.

Variable anaesthesia and joint morphology
Kruskal-Wallis tests show that anaesthesia from 3 to 4 dpf does not lead to any significant difference in joint spacing between the medial and lateral sides at 5 dpf relative to control fish ( Fig. 3A1-2, C). Anaesthesia from 4 to 5 dpf, 4 to 5þ dpf, 3 to 5 þ dpf or 5 dpf alone significantly alters the size of the interzone between the medial and lateral sides (Fig. 3A1, 3-6, C).
Outline traces indicate that the joint surface of the MC is more plastic than the surface of the palatoquadrate (Fig. 3Bi, ii). Anaesthesia at 3-4 dpf alone did not lead to a dramatic shape change of the MC (Fig. 3Bi-iii), as reflected in the interzone interval analysis (Fig. 3C), however, later anaesthesia leads to a change in joint morphology such that the MC approaches or overlaps the surface of the PQ on the medial side (Fig. 3Biii).

Functional strains within the developing cartilage jaw
FE-models demonstrate maximum principal strains present upon jaw closure (generated by adductor mandibulae muscle contraction) (Fig. 4A-I) and minimum principal strains at 3-5 dpf (Fig. 4J-R). We also generated models of jaw opening via contraction of the protractor hyoideus and the intermandibularis (Diogo et al., 2008) (Fig. 5). From these models it is apparent jaw closure causes maximum principal strains to be located medially and laterally at the interzone and minimum principal strains are located laterally (Fig. 4 D''-F'', M''-O''). During jaw opening, compressive strains (minimum principal strains) are recorded on the medial side of the MC joint and at the interzone (Fig. 5G-J). Relevant interzone strains were always between þ3500 mstrain (maximum) and À 5000 mstrain (minimum), in line with other models (Nowlan et al., 2008), (Supp. Fig. 5).

Joint cell orientation
Joint morphology changes from an initially asymmetric shape to an increasingly symmetrical shape from 3−5 dpf (Fig. 2) and in the absence of muscle contraction joint shape remains asymmetric with the medial side of the MC joint protruding over the PQ (Figs. 2 and 3). Minimum strains are located on the medial side during mouth opening. We, therefore, tested whether cells change their orientation over time to reflect the pattern of strain on the joint and interzone (Figs. 4 and 5). As the joint surface of the MC is more plastic than the surface of the PQ (Fig. 3), we only considered cells from this element. During normal joint morphogenesis, cells on the joint surface of the MC did not show a significant change in cell orientation, apart from those located at the medial side (Fig. 6A, Supp. Table 4). However, for zebrafish that were anaesthetised for 3-5 dpf (Fig. 6B) or myod mutants (Fig. 6C), a significant change in orientation in cells on the medial side of the joint was observed at 5 dpf compared to the 5 dpf control medial cell orientations (Supp. Table 4).

Discussion
Here we present data describing the morphology of the developing zebrafish jaw joint over time, and the effect of removing muscle loads and subsequent biomechanical strain on joint morphology. Normal, wild type joint morphology changes from an initial flattened shape to an interlocking morphology between 3 and 5 dpf. Building on previous work from mouse and chick models showing that joint formation and integrity require muscle activity during development (Drachman and Sokoloff, 1966;Kahn et al., 2009;Nowlan et al., 2010;Roddy et al., 2011b), we show muscle activity is critical for refining the morphology of the zebrafish jaw joint. Critically, our myod mutant and anaesthetisation experiments demonstrate that muscle activity is required for joint integrity, rather than the presence of muscle per se. By performing temporal experiments in which we anaesthetise fish for differing time periods we show early muscle activity (from 3 to 4 dpf) is dispensable for normal joint morphogenesis, but movement during later time periods is required for normal joint shape. Therefore the joint retains a degree of plasticity such that loss of early muscle activity does not lead to significantly altered morphology, similar to the situation in chick (Nowlan et al., 2014). In part this may be because muscle activity does not peak until day 4, therefore fewer mechanical cycles will   have been experienced during this period. Alternatively, recovery of muscle activity could be sufficient to recover altered joint shape. The MC surface of the jaw joint is more affected by muscle paralysis than the PQ, suggesting the PQ is less mechanosensitive than the MC. Differential mechanosensitivity of bones has previously been described in limbs (Nowlan et al., 2008) and the mandible (Enomoto et al., 2014).
Our models represent the first FE models for the zebrafish craniofacial skeletal system. Zebrafish are increasingly used as a model for biomechanics, but so far these studies have focused mainly on the ontogeny of swimming behaviour (Fiaz et al., 2012;Fontaine et al., 2008;Green et al., 2011;Li et al., 2012), though studies exist that describe the onset of different jaw movements for respiration and feeding (Hernandez et al., 2002(Hernandez et al., ., 2005 and the description of early cell behaviour in forming the jaw elements is well described (Eames et al., 2013;Talbot et al., 2010). Our models show that the strains occurring in the interzone, which drive the changes to zebrafish joint shape are similar in magnitude (up to þ3500 m strain) to those seen in comparable models such as embryonic chick limb (Nowlan et al., 2008), and those used in vitro to elicit cell behaviour changes that drive joint cavitation (Dowthwaite et al., 2003(Dowthwaite et al., ., 1999. Higher strains were generated at constraints and muscle attachment points, but are unlikely to have a significant impact on the model, as the location rather than magnitude of strain is primarily studied. We show that cell orientation within the joint is significantly affected in fish subjected to a period of immobility. This demonstrates that the cells are altering their behaviour in a strain dependant fashion. Cells on the medial side of the joint change their orientation over time partially accounting for the natural change in joint shape; whereas, in joints where biomechanical strains are reduced, medial cells fail to adopt the correct orientation, leading to an overlapped, non-functional joint. The role of strain in mediating cell orientation in vitro is well characterised for many cell types including chondrocytes (Ghezzi et al., 2014) and ex vivo (Clark et al., 2003).
The FE models show that strains during mouth closure are higher than those during opening. During jaw closure minimum Fig. 6. Cell orientations in wild type, anaesthetised and myod mutant zebrafish jaw joints. Orientation angle of chondrocytes in (A) 3-5 dpf control; (B) 3-5 þ dpf anaesthetised; and (C) myod mutant zebrafish in the Meckel's cartilage element of the joint, plotted on circular histograms (rose plots), where 0°lies on the medial side of the joint and 180°at the lateral side of the joint. n¼ 3 joints per experimental condition (1 or 2 refers to number of joints per blue wedge). The number of cells analysed per condition are listed in Supp. Table 2. Histogram bins equal 20°. The red line marks mean orientation and the green line marks the 95% confidence interval.
(compressive) principal strains are higher than maximum strains and predominantly seen laterally, loss of strain during immobilisation could lead to a change in cellular processes on the lateral side, with a secondary effect on the medial side. Alternatively, changes to strain location may arise in the morphologically altered immobilised jaws.
During mouth opening peak minimum strains are located at the medial side of the joint interzone. The loss of compression at the medial side of the joint and interzone, could explain the change in chondrocyte behaviour that leads to a switch in cell orientation at the medial side in anaesthetised and myod mutant zebrafish. The change in orientation, while significantly contributing to the change in joint shape does not completely explain it, suggesting that there is an additional component, which could be a change to rates of differentiation, proliferation, or migration of chondroprogenitors of the joint. Indeed the failure of chondrocytes to mature fully and intercalate has been previously reported (Shwartz et al., 2012).
We believe that adding 3-dimensional modelling of the joints using FEA to this emerging field will further increase the utility of the zebrafish model for biomechanical studies. Better description of the ontology of the skeletal system that relates zebrafish to higher vertebrate models (Dahdul et al., 2012) will also facilitate future comparisons between model organisms for skeletal biology. Excitingly, use of the zebrafish model with its genetic amenability, opens up the possibility of directly manipulating mechanosensitive genes using transgenic tools to switch genes on or off in the skeletal system, unpicking the relationship between genes and biomechanical influences in shaping a joint.

Conflict of interest statement
All authors of this manuscript state that they have no conflicts of interest to report. The funders of the study (all of whom are listed in the acknowledgements section of the manuscript) were in no way involved in the experimental design or in the writing of the manuscript.