The opercular mouth-opening mechanism of largemouth bass functions as a 3D four-bar linkage with three degrees of freedom

The four-bar linkage, a mechanical device consisting of four rigid links interconnected to form a closed loop, is an important model of the transmission of motion and force in musculoskeletal systems. Many biomechanical systems have been modeled as four-bar linkages, including portions of the jaw apparatus in some lizards (Frazzetta, 1962; Metzger, 2002) and birds (Bock, 1964; Van Gennip and Berkhoudt, 1992), the striking appendages of mantis shrimps (Patek et al., 2007), the ligaments of the human knee (Greene, 1983), the thorax and wings of flying insects (Wootton, 2009), and multiple skeletal elements in the skulls of fishes (Ballintijn, 1969; Westneat, 1990). In addition, four-bar models have been used to evaluate functional hypotheses by comparing simulated and in vivo kinematics (Westneat, 1991; Van Wassenbergh et al., 2005; Roos et al., 2009), to measure how force and motion are transmitted through musculoskeletal systems (Aerts and Verraes, 1984; Adriaens et al., 2001; Van Wassenbergh et al., 2013), and to examine the distribution and evolution of functional diversity (Westneat, 1995; Alfaro et al., 2004, 2005; Wainwright et al., 2004; Hulsey and García de León, 2005).

In biomechanical studies, four-bar models have been applied most extensively to the skulls of fishes. Relative to other vertebrates, the heads of ray-finned fishes are extremely mobile, containing over 15 movable skeletal elements (Gregory, 1933; Schaeffer and Rosen, 1961; Alexander, 1967; Lauder, 1983; Ferry-Graham and Lauder, 2001), joined together by a structurally and functionally diverse set of articulations (Liem, 1980; Lauder, 1982; Westneat, 2006). In most fishes, these mobile elements are commonly subdivided into a series of interconnected mechanisms that can be modeled as four-bar linkages: the hyoid-pectoral mechanism (Muller, 1987), the anterior jaw mechanism (Westneat, 1990), and the opercular mechanism (Fig. 1; Ballintijn, 1969). These mechanisms enable the skull to convert force inputs from relatively few muscles into rapid, three-dimensional (3D) expansion of the mouth cavity to generate suction and engulf prey items during feeding (Muller, 1989; Carroll, 2004; Carroll and Wainwright, 2006; Camp and Brainerd, 2015; Camp et al., 2015).

Previous applications of the four-bar linkage to biological systems have primarily used a 2D, or planar, four-bar (e.g. Aerts and Verraes, 1984; Westneat, 1990; Alfaro et al., 2004; Van Wassenbergh et al., 2005; Patek et al., 2007; Roos et al., 2009). In the 2D four-bar, all the joints are hinges with their rotational axes oriented in parallel, confining all rotations to a single plane (Fig. 2B). The 2D four-bar also has a single degree of freedom (DoF), the total number of inputs needed to fully characterize the position and orientation of each link. As a consequence, the rotation of one link is sufficient to resolve the conformation of the entire linkage. Despite its simplicity, the 2D four-bar model is generally accurate in reproducing the motion of several biological systems, including the striking appendage of mantis shrimp (Patek et al., 2007; McHenry et al., 2012), the hyoid depression mechanism in catfish (Van Wassenbergh et al., 2005; with some deviation observed in one individual), and the anterior jaw mechanisms in wrasses (Westneat, 1990, 1991, 1994).

However, one system in which a 2D four-bar model has proven consistently inaccurate is the opercular mechanism in fishes, a coupling of the operculum, suspensorium, lower jaw and interoperculum that depresses the lower jaw (Fig. 1; Ballintijn, 1969). Although originally proposed as a 2D four-bar based on anatomical observations (Anker, 1974; Elshoud-Oldenhave and Osse, 1976; Aerts and Verraes, 1984), subsequent tests have shown substantial error between in vivo rotation of the lower jaw and the rotation predicted by a 2D four-bar: ranging from 20 to 50% in wrasses (Westneat, 1990, 1994), over 45% in catfish after initial mouth opening (Van Wassenbergh et al., 2005), and up to 50% in largemouth bass (Camp and Brainerd, 2015). At the same time, lower jaw depression is temporally correlated with opercular rotation (Westneat, 1990, 1994; Van Wassenbergh et al., 2005; Camp and Brainerd, 2015) and surgical severance of the linkage results in disruption of lower jaw depression (Liem, 1970; Durie and Turingan, 2004), suggesting some form of mechanical coupling.

The inaccuracy of the 2D four-bar in predicting opercular mechanism motion may be due to any of several simplifying assumptions, including planar and single-DoF motion. In largemouth bass the operculum moves with three rotational DoF and the lower jaw flares out laterally while depressing ventrally (Camp and Brainerd, 2015), violating the 2D four-bar assumptions of planar motion and single-axis rotations. In addition, the interoperculum (coupler) link is formed not by a single bone, but rather by a serial arrangement of the interoperculomandibular (IOM) ligament, interoperculum bone and suboperculum bone (Fig. 1). During jaw depression, tension could stretch the ligamentous component of this link, violating the more fundamental four-bar assumption of rigid links. If stretching of the interoperculum link contributes relatively little to 2D four-bar model error, then a higher DoF, 3D four-bar would provide an accurate model for the opercular mechanism. While a 3D four-bar model has been proposed previously for biological systems (Olsen and Westneat, 2016), such a model has never been tested against in vivo kinematics.

In this paper we present new 3D four-bar linkage models with varying joint configurations and DoF. We fit each of these models to the in vivo 3D motion of the opercular mechanism in largemouth bass (Micropterus salmoides Lacepède 1802), collected using X-ray reconstruction of moving morphology (XROMM), to test whether the opercular mechanism functions as a four-bar linkage and, if so, to quantify how many degrees of freedom are needed to characterize its motion. We also compare the fit of these models at two different levels of variation, for each feeding strike and across multiple feeding strikes from the same individual. We evaluate the fit of each model based on a comparison with an error benchmark derived from inter-individual variation and ANCOVA (analysis of covariance) tests of the relationship between output and input rotation. If the bass opercular mechanism functions as a four-bar linkage, then the fit of at least one of the four-bar models should either be significantly less than our established error benchmark or produce a ratio of output to input rotation that does not differ significantly from in vivo. We discuss the evaluation and comparison of multiple models, the use of four-bar models to simulate opercular mechanism kinematics in fishes, and the prediction of motion from morphology in the skulls of fishes more broadly.

We measured the in vivo 3D kinematics of the opercular mechanism in largemouth bass during suction feeding using XROMM data collected as part of a previous study (Camp and Brainerd, 2015) (Supplementary Information 1). XROMM generates accurate 3D skeletal animations from biplanar X-ray videos and 3D digital bone models (Brainerd et al., 2010). Radio-opaque beads were surgically implanted in the neurocranium and left-side operculum, suspensorium, interoperculum (Bass2 and Bass3 only) and lower jaw (Fig. S3) of three largemouth bass (standard lengths of 307, 287 and 316 mm). Suction feeding strikes on live goldfish prey were then recorded using biplanar high-speed X-ray video for a total of 28 feeding strikes across all individuals (nine strikes were recorded from Bass1 and Bass2, and 10 strikes from Bass3). The 3D motion of each radio-opaque bead was reconstructed from these videos, with a mean marker-tracking precision of 0.11 mm (Camp and Brainerd, 2015). The bead motions were combined with polygonal-mesh models of the skeletal elements generated from computed tomography (CT) scans using the XROMM Maya Tools package (available at xromm.org) to create 3D skeletal animations (the methods are described in more detail in Camp and Brainerd, 2015). All experimental procedures were approved by the Brown University Institutional Animal Care and Use Committee (IACUC).

We used the XROMM animations to measure the motion of three non-collinear points from each of the three main skeletal elements of the opercular mechanism (Fig. S1): the lower jaw, suspensorium and operculum. These points (referred to here as ‘fit points’) were subsequently used in fitting the four-bar models. Of these nine fit points, five were located at the joints of the hypothesized opercular mechanism (Fig. 1A): one suspensorium-associated point at the suspensorium–lower jaw joint, one lower jaw-associated point on the angular process where the IOM ligament attaches to the lower jaw, one operculum-associated point where the interoperculum attaches to the operculum, and two points at the suspensorium–operculum joint (one associated with the suspensorium and one associated with the operculum). The locations of the remaining four fit points were selected to minimize collinearity and maximize spread among the three points associated with each element so that the fitted four-bar models could adequately capture the full 3D rotation and translation of each element. Point positions were then measured during the strike as each point moved with its associated polygonal-mesh bone model.

Motion data were recorded from one individual at 300 frames s⁻¹ and from two individuals at 500 frames s⁻¹. To analyse all data in a consistent time frame, the xyz-coordinates of the fit points at 500 frames s⁻¹ were converted to 300 frames s⁻¹ by linear interpolation using a custom function (Supplementary Information 1). To remove the motion of the entire fish's body within the tank, all motion was measured relative to a fixed operculum. The operculum was chosen as the reference because it has previously been shown to move the least relative to the body, among the opercular mechanism elements (Camp and Brainerd, 2015). Early analyses confirmed that results did not differ substantially whether the operculum or suspensorium, the traditional fixed link of the opercular four-bar (e.g. Liem, 1980), was chosen as the fixed link.

To characterize the rotations and translations of the opercular mechanism in an anatomically relevant way, we measured the 6-DoF motion of the lower jaw (output link) and suspensorium (input link) using joint coordinate systems (JCS). A JCS consists of two anatomical coordinate systems, one that remains fixed to a bone and another that moves with a second bone, allowing motion to be described as a set of translations and Euler rotations along and about the three JCS axes (Grood and Suntay, 1983; Brainerd et al., 2010). The lower jaw JCS was placed at the quadratomandibular joint and used to measure motion of the lower jaw relative to the suspensorium (Fig. 3A). The suspensorium JCS was placed at the suspensorium–operculum joint and used to measure motion of the suspensorium relative to the operculum (Fig. 3D). The JCS axes follow the right-hand rule and correspond to the major body axes: positive rostrocaudal (RC) translation corresponds to rostral translation; positive dorsoventral (DV) axis rotation and translation correspond to lateral rotation and dorsal translation, respectively; and positive transverse axis rotation and translation correspond to elevation and rightward translation, respectively. For each JCS, translations and rotations were calculated relative to their initial value prior to the onset of the suction feeding strike (67 ms prior to peak gape). These measurements were used only to describe the 3D motion of the lower jaw and suspensorium with anatomically relevant axes. For fitting all four-bar models, we estimated best-fit centers and axes of rotation, as described in the ‘Model fitting’ section, below.

Two sources of variability in our kinematic dataset that are not captured directly in the JCS analysis or the four-bar models are changes in length of the coupler link and bone deformation. We measured each of these additional sources of variability to quantify their potential contribution to model fit error. The coupler link of the opercular mechanism is not formed by a single rigid link but rather by the interoperculum bone and the IOM ligament in series (Fig. 1). As a consequence, stretching or slackening of the ligament, folding between the bone and ligament, or deformation of the interoperculum bone could cause the coupler link to change length. We quantified each of these sources of variation by measuring length changes within the interoperculum bone (as the distance between the interoperculum–operculum joint and the IOM ligament attachment on the interoperculum), length changes in the entire link (as the distance between the interoperculum–operculum joint and the attachment of the IOM ligament on the lower jaw), and, for the two individuals in which the interoperculum bone was animated, length changes in the IOM ligament (as the distance between landmarks at its attachments on the dentary and interoperculum). To measure bone deformation we selected the two radio-opaque beads in each skeletal element positioned at the greatest distance from one another, thereby maximizing our ability to capture bending or deformation of the bone. We then measured the distance between the beads throughout each strike. If the skeletal element is a completely rigid body, the change in the average distance between the beads over time should be minimal (similar to the marker-tracking precision of 0.11 mm) and occur randomly throughout the strike. However, a larger and consistent change in the distance over time suggests that the skeletal element is not a rigid body and that there is deformation within the bone (operculum, interoperculum) or motion between bones making up the skeletal element (suspensorium, lower jaw).

To test whether the opercular mechanism moves in a manner consistent with a four-bar linkage we fitted five different four-bar models (Fig. 2), varying in joint configurations and degrees of freedom, to the fit points measured from the animated bone models (see ‘X-ray reconstruction of moving morphology’ section, above). The simplest model is the traditional 1-DoF 2D four-bar, consisting of four parallel hinge joints (also known as ‘revolute joints’ or ‘R-joints’), denoted by the name ‘RRRR’ (Fig. 2B). The next more complex model is a 1-DoF 3D four-bar with two non-parallel hinge joints and two ball-and-socket joints (also known as ‘spherical’ joints), denoted by the name ‘RRSS’ (Fig. 2C). The next two more complex four-bar models (RSSS and SRSS) each gain two additional degrees of freedom by exchanging one of the hinge joints in the RRSS model (at the suspensorium–operculum joint and the suspensorium–lower jaw joint, respectively) with a ball-and-socket joint (Fig. 2D,E). Lastly, the highest DoF model (SSSS) is a four-bar consisting entirely of ball-and-socket joints (Fig. 2F).

Rather than use the traditional approach of testing how well a predefined model (e.g. based on anatomical measurements) performs against in vivo kinematics (e.g. Westneat, 1990, 1994; Van Wassenbergh et al., 2005; Roos et al., 2009; Camp and Brainerd, 2015) we used a best-fit approach. A best-fit approach tests how well a given set of model assumptions (e.g. a 1-DoF, 2D four-bar) performs against in vivo kinematics by searching for a set of parameters that result in the lowest model error. A best-fit approach is also less sensitive to the uncertainty or choice of a particular set of model parameters. We fitted each model to the in vivo kinematics during a 67 ms time window (21 frames) prior to peak gape using the ‘fitMechanism’ function, newly added to the R package ‘linkR’ (http://cran.r-project.org/package=linkR; available for download from the CRAN repository; see ‘Data availability’ section). The more basic kinematic simulation (i.e. moving a linkage through possible conformations, also known as ‘position analysis’) is performed by the linkR function ‘animateMechanism’, which uses an incremental geometric assembly approach (Olsen and Westneat, 2016), similar in nature to, though independently developed from, Kramer's degrees of freedom analysis (Kramer, 1992). The fitMechanism function acts as a wrapper to animateMechanism, managing the parameter optimization based on the kinematic simulation results returned by animateMechanism.

Starting with an initial estimate of the four-bar joint coordinates, the orientation of any hinge axes and the reference pose (i.e. position and orientation) of the points associated with each skeletal element (here the ‘fit points’), the fitMechanism function optimizes the four-bar input parameters (e.g. the magnitude of input link rotation) at each time point to best align the model fit points and in vivo fit points (the number of input parameters at each time point is equal to the DoF of the four-bar). Following this, the fitMechanism function uses these input parameters to optimize the position of the lower jaw–IOM ligament joint (three parameters in 3D, two parameters in 2D), effectively finding an optimal length for the coupler link, which was found to stretch during the strike. During the optimization the joint was not allowed to translate more than 2 mm from its initial position, corresponding to the mean lengthening observed in the coupler link. Next, the fitMechanism function optimizes any hinge joint axes across all time points (three parameters per axis, corresponding to a 3D vector). Then the optimized input parameters and joint axes are used to optimize the reference pose of each element's fit points across all time points relative to the four-bar joint coordinates (six parameters per element, corresponding to a 3D translation and rotation). Each optimization step is based on the root mean square (RMS) error, which as a result of the square root imposes a disproportionately greater penalty on larger error values relative to a simple mean error.

All optimizations are performed using the R ‘stats’ package function ‘nlminb’, a non-linear optimization function that uses gradients to search a multiple-parameter space for a combination of parameter values that minimize a particular error function (http://www.R-project.org/). Because each optimization step affects the error of the preceding step this sequence of optimizations is repeated until consecutive error changes are less than 0.1%, in a manner analogous to Procrustes superimposition in shape analysis (Rohlf, 1999). Each four-bar model was fitted at two different levels of variation: by strike and by individual. For the former, hinge joint axes and the fit point reference poses were optimized for each strike separately. For the latter, hinge joint axes and the fit point reference poses were optimized for all strikes from the same individual. Fitting by individual was done by concatenating all strikes from the same individual into one sequence. As the model fitting is only positional (time independent), jumps between the start and end of consecutive strikes, and in fact the order of the time points themselves, has no effect on the model fit. To reduce the time required for optimization, hinge joint axes and fit point reference poses were optimized using a subset of 50 frames (24–26% of the total frames), evenly spaced across the concatenated frames; input parameter optimization and model fit error calculation were still performed using all frames and the use of more frames for hinge axis and reference pose optimization did not substantially change model fit errors.

For the 2D four-bar model, no optimization of the hinge joint axes is needed as all the axes are by definition orthogonal to the linkage plane. In addition, as we are simplifying the interopercular link as a coupler link in tension, we ignore long-axis rotations of the interopercular link and have not included any interoperculum-associated fit points. Thus, although the anatomical joints between the lower jaw and IOM ligament and between the interoperculum and operculum are not in fact ball-and-socket joints, treating them simply as joints with three rotational degrees of freedom enables the interopercular link to function as a simple constant-length coupler link. Given that all motion is measured relative to a fixed operculum, the only fit points included in the RMS error are the three fit points associated with the lower jaw and three fit points associated with the suspensorium.

To evaluate and compare the fit of each model we used two approaches: comparison against an error benchmark and ANCOVA tests of the relationship between linkage output and input rotation. For the first approach, we compared model fit errors to a set of Procrustes consensus distances. One intended use of these four-bar models is to compare musculoskeletal function across species. Thus, we would like to obtain model fit errors that do not exceed the variation among individuals so that comparisons of model results reflect true interspecific differences and not simply model errors. One measure of inter-individual variation is the distance of homologous points from the corresponding mean or consensus points among individuals.

To calculate this, Procrustes superimposition was first used to find a mean coordinate set for the four opercular mechanism joints across all three individuals; this process finds an optimal alignment of multiple coordinate sets by minimizing differences in centroid size, translation and rotation (Rohlf, 1999). Next, distances were measured between the joint coordinates of each individual and the corresponding mean positions at five evenly spaced times through the expansive phase of a strike, creating a sample of 60 consensus distances. Model error was calculated by finding the mean distance between the model and in vivo fit point positions at peak gape (positional error). The ‘lm’ function in the ‘stats’ R package was then used to test whether the mean model fit errors at peak gape were significantly greater than the mean consensus distance using the following model formula: d∼src+indiv. Consensus distance and error values were then modeled as the predicted variable (d), predicted by the categorical factors ‘src’ (whether the values corresponded to consensus distance or error) and ‘indiv’ (which individual values were measured from). This test takes into account any potential individual-specific effects, such as if results from only a single individual were driving the significant difference. A more complex model that included interactions between ‘src’ and ‘indiv’ produced the same conclusions and did not provide a consistently better fit [determined by difference in Akaike Information Criterion (AIC) score]; thus, the simple non-interaction model was used.

For the second approach, we used ANCOVA tests to compare the in vivo relationship between output rotation (total lower jaw depression) and input rotation (total suspensorium elevation) to that predicted by each model, as done in previous evaluations of four-bar model error (Westneat, 1990, 1991, 1994). Although the relationship between output and input rotations in linkages is non-linear, it is often sufficiently linear for the limited range of motion under consideration to use a linear regression (e.g. Westneat, 1990, 1991, 1994). We performed the ANCOVA tests using the ‘lm’ function in the ‘stats’ R package, using the following formula: ljaw∼susp×src+indiv. Lower jaw or output rotation (‘ljaw’) was predicted by suspensorium or input rotation (‘susp’), with the categorical factors ‘src’ (whether values were taken from in vivo or the model) and ‘indiv’ (which individual values were measured from). We then used the ‘anova’ function (‘stats’ package) to test for significant differences between the in vivo kinematics and the four-bar model predictions, including a difference in slope, or a difference in either the intercept or mean of the output or input rotation. The inclusion of an individual factor accounts for any potential individual-specific effects. A more complex model that included interaction terms for the factor ‘indiv’ produced the same conclusions, and for all but one model had a higher AIC score (worse fit); thus, the simpler model without interaction terms for ‘indiv’ was used.

Although an ANCOVA can detect model rotational errors, it cannot necessarily detect model positional errors (i.e. predicted lower jaw rotation may be accurate but its predicted position is inaccurate). Thus, the positional error benchmark and ANCOVA tests provide complementary evaluations of model error. In addition to the benchmark and ANCOVA approaches we also calculated the angular error in predicted total suspensorium and lower jaw rotation separately. Lower jaw rotation was measured as the rotation of the vector defined from the IOM ligament–lower jaw attachment to the lower jaw–suspensorium joint. Suspensorium rotation was measured as the rotation of the vector from the lower jaw–suspensorium joint to the suspensorium–operculum joint. Vector rotations were measured as the change between the start and end of each strike. As we had no clear benchmark for angular errors, this analysis lacked significance statistics and served instead to provide a measure of model performance more intuitive in relation to opercular mechanism function. All P-value adjustments for multiple comparisons were performed using the ‘p.adjust’ function in the R ‘stats’ package.

Lower jaw and suspensorium motion showed substantial deviation from the sagittal plane (i.e. deviation from rotation about the transverse axis) during mouth opening (Figs 3 and 4). While lower jaw and suspensorium rotations were greatest about the transverse axis of each JCS, with mean (±s.e.m.) rotations of −36.5±1.2 and 12.1±0.35 deg at peak gape, respectively, there were also substantial rotations about the RC and DV axes (Fig. 4C,F). The rotational components of the lower jaw about the RC and DV axes were −12.5±1.2 and −8.56±0.8 deg, respectively, while those of the suspensorium about the RC and DV axes were 4.19±0.7 and −4.5±1.2 deg, respectively (Fig. 4C,F). Mean translations of the suspensorium at the suspensorium–operculum joint (relative to the operculum) were minimal (generally <1 mm) along all three axes (Fig. 4E and Fig. S2). Lower jaw translations were minimal along the RC and DV axes, but the quadroarticular joint of the lower jaw translated laterally a mean distance of 2.9±0.3 mm relative to the suspensorium (Fig. 4B and Fig. S2).

The interoperculum (coupler) link showed substantial length changes during mouth opening and closing (Fig. 5) in the two individuals where the interoperculum bone was included in the XROMM animation (N=19 strikes). Defined by joints positioned at the ventral portion of the suboperculum and the articular process of the lower jaw (Fig. 1), the interoperculum link lengthened by a mean±s.e.m. of 1.8±0.2 mm or 5.3±0.5% strain at peak gape, relative to its length prior to the start of the strike (Fig. 5). We also measured the length of the IOM ligament as the distance between its attachment sites on the interoperculum and the articular process. The IOM ligament lengthened by a mean of 0.6±0.1 mm or 5.1±0.8% strain at peak gape, relative to its length prior to the strike (Fig. 5). As expected, the interoperculum bone itself underwent relatively small length changes corresponding to less than 2% strain.

For two of three individuals (Bass2 and Bass3), the lower jaw, suspensorium and operculum mean inter-bead distances changed consistently over time, with the greatest distance generally coinciding with peak gape (Figs S3, S4). The lack of change in the distance for Bass1 is probably due to the beads being implanted closer together (e.g. all lower jaw beads were placed within the dentary). Although this pattern of inter-bead distance changes suggests that these skeletal elements are not completely rigid bodies and may be deforming, the magnitude of these changes was relatively small. Mean inter-bead distances changed by less than 0.3 mm (0.8% initial length) for the lower jaw, 0.4 mm (1.7% initial length) for the operculum, and 0.5 mm (1% initial length) for the suspensorium across all three individuals (N=28 strikes).

As expected, the positional error, how well the model fit points overlapped with the in vivo fit points, decreased as a function of the model's DoF (Fig. 6A, Table 1; Movie 1). However, differing joint configurations for the same DoF also influenced fit error. For example, the RSSS and SRSS models, although both having three DoF, had significantly different errors by strike and by individual (non-overlapping notches between boxes in Fig. 6 generally correspond to significantly different median values). For the RRRR model fit by strike, the model fitting failed to converge for two strikes, probably due to the poor fit of the RRRR model more generally; thus, summary statistics for the RRRR model by strike include 26 rather than 28 strikes. Also as expected, model fit errors were higher when fitted to all strikes across an individual than when fitted to each strike separately. However, we also found that whether fitting by strike or by individual, at least three DoF are needed to obtain an error at or below the inter-individual benchmark (Table 2). By strike and by individual the 3-DoF model SRSS is at or below the benchmark. Thus, for both levels of variation the SRSS model represents a best-fit model, that is, the simplest model (having the fewest DoF) at or below the error threshold.

As with positional error, the mean errors for suspensorium and lower jaw rotations were generally greater when fitted by individual than when fitted by strike (Fig. 6B,C, Table 1). There was also a general decrease in error with increasing model DoF, with the 3-DoF SRSS model and the 5-DoF SSSS models having the lowest mean errors (<1.5 deg). Angular errors measured relative to the total rotation of each body were generally greater for the suspensorium than for the lower jaw. All the models tended to underestimate suspensorium rotation while there was no consistent tendency to under- or overestimate lower jaw rotation.

The linear models used for each ANCOVA test fit moderately well, with R² values ranging from 0.64 to 0.92 (Fig. 7). The ANCOVA tests identified one model as differing significantly from in vivo by strike (RRRR) and three models by individual (RRRR, RRSS and RSSS), consistent with a lower model fit error when fitting each strike separately. Thus, the ANCOVA tests indicate that by strike one DoF is sufficient to accurately represent the in vivo suspensorium–lower jaw relationship, while by individual at least three DoF are necessary. The results of the ANCOVA tests were generally consistent with those of positional error with the exception of the RRSS and RSSS models by strike. When the benchmark and ANCOVA results are considered together, both the 3-DoF SRSS model and 5-DoF SSSS model satisfy the desired level of both positional and rotational accuracy.

For nearly three decades the precise constraints guiding the motion of the opercular mechanism in fishes have remained a frustrating puzzle. Anatomical observations and surgical disruption suggested that this mechanism was analogous to a four-bar linkage (Anker, 1974; Elshoud-Oldenhave and Osse, 1976; Aerts and Verraes, 1984). Yet previous planar four-bar models of the kinematics have been unsuccessful in accurately matching the evident temporal correlation between opercular rotation and lower jaw depression in several species of ray-finned fishes (Westneat, 1990, 1994; Van Wassenbergh et al., 2005; Camp and Brainerd, 2015). The results of this study provide resolution to this question by demonstrating that, for at least one species of ray-finned fish, the opercular mechanism functions as a 3D four-bar linkage with three degrees of freedom. A comparison of multiple four-bar models identifies a four-bar model with a ball-and-socket joint at the suspensorium–operculum joint, a hinge joint at the suspensorium–lower jaw joint, and a link of constant length connecting the lower jaw to the operculum (referred to here as the SRSS model) as the simplest model to accurately match the position and rotation of the suspensorium and lower jaw across multiple strikes (Figs 6, 7, Tables 1, 2; Movie 1).

A 3D, multiple-DoF four-bar has not been used previously to model a biomechanical system. Prior tests of four-bar models against biomechanical systems have primarily evaluated 2D, 1-DoF models (e.g. Aerts and Verraes, 1984; Westneat, 1990; Alfaro et al., 2004; Van Wassenbergh et al., 2005; Patek et al., 2007; Roos et al., 2009). Although other studies have evaluated the fit of 3D multi-joint models against in vivo or passive kinematics, including models with more than one DoF, none of the models has taken the form of a 3D four-bar linkage (Reinbolt et al., 2005; Di Gregorio et al., 2007; Chang and Pollard, 2008; Franci et al., 2009; Cerveri et al., 2010; Fohanno et al., 2014; Laitenberger et al., 2015). Thus, this study expands the versatility of four-bar models in simulating biomechanical systems and demonstrates that their utility is not limited to planar or single-DoF systems.

The error in predicted lower jaw depression by the 2D four-bar (RRRR) model in this study, 8% by individual (Table 1), appears less than that of a previous study (Camp and Brainerd, 2015), which obtained 2D four-bar model errors of 50% using the same kinematic dataset. However, the lower jaw positional error is high: 43% relative to the total mean displacement of the lower jaw and suspensorium fit points. In addition, the optimization approach used in this study is fundamentally different from the previous approach, which compared observed kinematics with a 2D four-bar model derived from anatomical measurements. Optimizing each model for each individual will result in lower model fit errors but can also shift error to different mechanism elements: for example, angular errors for the 2D four-bar shifted from the lower jaw to the suspensorium, which deviates 25% by individual from in vivo rotation (Table 1). The results of this study confirm that even with optimization, a 2D four-bar remains an inadequate model for the bass opercular mechanism.

The inability of previous 2D four-bar models to match opercular mechanism kinematics is probably due to the fact that these models were both two-dimensional and permitted only a single DoF. A 3D, 1-DoF four-bar model is inaccurate in predicting the position and rotation of the opercular mechanism across multiple strikes (the RRSS model in Figs 6,7). Thus, simply using a 3D, rather than 2D, four-bar is insufficient to produce an accurate four-bar model for this system. Rather, an increase in the rotational degrees of freedom at the suspensorium–operculum joint is needed, particularly to match the mechanism's range of motion across multiple strikes. These extra degrees of freedom imply a role for active motor control in determining the mechanism's particular trajectory. A previous study found that 95% of the power required for mouth opening in largemouth bass comes from the axial muscles that run down the length of the body, while the cranial muscles (the levator arcus palatini, dilator operculi and levator operculi muscles) contribute relatively little power (Camp et al., 2015). These cranial muscles may play a role as fine-tune controllers, directing the relatively imprecise force generated by the axial muscles through this high-mobility system to more precisely achieve a particular mouth-opening trajectory. These extra degrees of freedom also imply that suspensorium rotation does not necessarily result in lower jaw rotation, as would occur in a 1-DoF linkage.

A 3D, 3-DoF four-bar model performs well in spite of several experimental and anatomical sources of variation. Firstly, if skeletal joints function as pure hinge or ball-and-socket joints, there should be no translation of the two bones relative to each other. Thus, non-zero translations of the suspensorium relative to the operculum (suspensorium–operculum joint) or the lower jaw relative to the suspensorium (lower jaw joint) measured from the XROMM animations reflect either animation errors, that the joint coordinate system was not placed at the exact center of joint rotation, or both. These sources of error probably explain the small (generally <1 mm) and highly variable (within and among individuals) translations measured along all axes of the suspensorium–operculum joint and along the RC and DV axes of the lower jaw joint. However, we measured relatively large (mean peak of 3 mm) and consistent lateral translation of the lower jaw joint relative to the suspensorium. We propose that this translation is the result of deformations of the lower jaw and/or suspensorium – both of which are composed of multiple bones – as well as sliding at this hinge joint, possibly as a result of hydrodynamic forces or those applied by the upper jaw. While the source of these joint translations is intriguing, they do not significantly affect the model fit (Fig. 7, Table 1).

Secondly, the coupler link of the opercular mechanism, formed by the interoperculum and IOM ligament, stretches by 5% (2 mm) of its initial length at peak gape (Fig. 5), violating the four-bar model assumption of rigid link lengths. While the coupler link length was optimized as a part of the model fitting, its length remained constant throughout the simulation. Thus, dynamic coupler length changes are not essential to achieve accurate model performance. Thirdly, changes in the distance between radio-opaque beads embedded in the same skeletal element reached as high as 0.5 mm in some bones (Fig. S3), suggesting the lower jaw, suspensorium and operculum may not be fully rigid skeletal elements. As the fit points for each element were exported from a rigid body model these errors are absorbed into displacement of each element, which may not necessarily be accommodated by the four-bar joint constraints. Relative to these sources of variation, the mean consensus distance of 1.1 mm is a relatively conservative error benchmark: it is 10 times greater than the bead tracking precision, two times greater than error attributable to bone deformation, half the length change of the coupler link, and three times less than the translation at the lower jaw–suspensorium joint.

Studies evaluating kinematic models against in vivo motion typically consider a single model. Such a model may be based solely on the anatomy of the individual or species (e.g. Hoese and Westneat, 1996; Bergert and Wainwright, 1997; Patek et al., 2007; Roos et al., 2009; Camp and Brainerd, 2015) or based initially on anatomy and subsequently optimized to fit motion data, as is more common in studies of human biomechanics (e.g. Sommer and Miller, 1980; Reinbolt et al., 2005; Chang and Pollard, 2008; Franci et al., 2009; Cerveri et al., 2010; Fohanno et al., 2014). However, consideration of a single kinematic model does not allow selection of a best-fit model from among several hypothesized models or the evaluation of how model assumptions influence its accuracy. While some previous studies have compared multiple kinematic models (e.g. Van Wassenbergh et al., 2005; Di Gregorio et al., 2007; Cerveri et al., 2008; Laitenberger et al., 2015), these studies compared models having the same or similar degrees of freedom.

This study advances biomechanical model fitting approaches by evaluating multiple models of varying degrees of freedom in relation to in vivo kinematics. There is an a priori expectation that models with greater DoF will have lower error: extra degrees of freedom can accommodate noise (inherent in in vivo data) rather than ‘true’ motion. For this reason, one cannot simply designate as the ‘best-fit model’ that which has the lowest error. A common approach to comparing models of varying complexity is to use a statistic that takes into account both the error and the number of parameters, such as AIC (Myung, 2000; Johnson and Omland, 2004). However, the use of such a statistic in comparing mechanism models is not entirely straightforward. AIC and similar measures depend on a likelihood function to evaluate not simply model error but the likelihood of a particular set of parameter values given the joint probability distribution of the parameters. Such a function is not easily calculable for a mechanism model and we know of no such calculation in the literature. In addition, one must assign relative weights to balance the number of parameters against model error, the designation of which, at this point, would be rather arbitrary.

Rather, we have identified the best-fit model as the simplest model that satisfies two different error criteria relevant to how mechanism models are used in the field of biomechanics: positional error at or below an inter-individual variation benchmark and rotational error measured as the ratio of output to input link rotation. Thus, although both the 3-DoF SRSS and the 5-DoF SSSS models satisfy both error criteria, the SRSS model is identified as the best-fit model as it is the simplest of these two models. These error criteria provide complementary measures of model performance. Positional error (Fig. 6) evaluates the full displacement error (translation and rotation) of each element rather than simply rotation. However, positional error is also sensitive to the choice and distribution of the fit points and significance depends on the choice of error benchmark. For this study, the mean distance to the Procrustes consensus shape was chosen as the benchmark because it is measured in a manner comparable to positional error (i.e. distance from model to in vivo positions) and represents the mean shape deviation among individuals. As a common objective in the field of comparative biomechanics is to compare model predictions among species, model error should not exceed the shape differences among individuals of the same species.

However, a comparison of output to input link rotation using an ANCOVA test (Fig. 7) relies on intuitive line-fitting to detect significant differences between model predictions and in vivo and is thus not dependent on the choice of a particular error benchmark (Westneat, 1990, 1991, 1994). In addition, an evaluation of output–input rotation error relates directly to the mechanical advantage of a linkage, a metric commonly used in biomechanics, evolution and ecology to understand organismal function and functional diversity (Aerts and Verraes, 1984; Westneat, 1995; McGee et al., 2013; Hu and Albertson, 2014; Jamniczky et al., 2014). However, a lack of a significant difference in output–input rotation does not necessarily indicate the model has correctly predicted an element's position: a model may correctly predict the rotation of an element while incorrectly predicting its position. For example, no significant difference is detected between model and in vivo output–input rotation for the RRSS model by strike (Fig. 7B), yet the positional error significantly exceeds the error benchmark (Fig. 6A). In general, seemingly minor model errors based on ANCOVA tests (e.g. RRRR model results in Fig. 7A,F) have among the highest positional errors, underscoring the importance of considering multiple error assessments.

By comparing the fit of models varying in their DoF for the same behavior at two different levels of variation (i.e. by strike and by individual) we are able to ask, for the first time, how DoF are partitioned within a biomechanical system. For example, are multiple DoF manifest within the trajectory of each feeding strike or do they arise solely as a result of variation in the trajectories taken across multiple feeding strikes? Our finding that a 1-DoF model is inadequate even at the level of a particular strike (Table 2) supports the former, with individual strikes exhibiting deviations from a 1-DoF mechanism. However, the substantially better model fits by strike (versus by individual; Figs 6,7) show that strike-to-strike variation also contributes to the opercular mechanism's multiple degrees of freedom. Ultimately, it is the fit by individual that is meaningful in determining the best-fit model: all motion constraints acting on the opercular mechanism may not be evident during a particular strike; however, they should become increasingly apparent as the number of observed strikes increases. By that same principle, our results represent minimum DoF for this mechanism as we have only made observations of its motion during a particular behavior and for a particular number of strikes. However, the fact that a 3-DoF model is the best-fit model whether considering one or nine to 10 strikes is reassuring that these results are not overly sensitive to the number of strikes examined.

Planar four-bar models of the opercular mechanism have been used extensively to understand evolution, development and ecomorphology across many groups of fishes (Aerts and Verraes, 1984; Westneat, 1995; McGee et al., 2013; Hu and Albertson, 2014; Jamniczky et al., 2014). More broadly, planar four-bar models have been used to predict motion from morphology in fish and other animals (Hoese and Westneat, 1996; Bergert and Wainwright, 1997; Patek et al., 2007; Roos et al., 2009). The use of these models has been essential to understanding the evolution of functional diversity because it enables rapid comparisons across a broad sampling of species. It is not possible to make any general conclusions as to how previous applications of the 2D four-bar opercular model should be reinterpreted. If model predictions are consistently inaccurate across species then interspecific differences or trends may still hold, despite being inaccurate in the case of each individual species. In addition, while the poor fit of 2D four-bar models to the opercular mechanism in species as varied as wrasses (Westneat, 1990, 1994), catfish (Van Wassenbergh et al., 2005) and bass (Camp and Brainerd, 2015; this paper) suggests that these results are broadly applicable, the 3D kinematics of the opercular mechanism have yet to be examined outside of largemouth bass. There may be substantial diversity in the function of the opercular mechanism across fish species, including mechanisms that operate with less – or more – DoF than in largemouth bass.

While previous studies in several fish species have shown a 2D four-bar to be a poor predictor of opercular mechanism motion, these studies had raised the possibility that a 3D, 1-DoF four-bar may provide an adequate model. Although the results of this study demonstrate that 1-DoF four-bar models, whether 2D or 3D, are inadequate for the opercular mechanism, they also recover the utility of morphology in predicting the motion of this mechanism. A 3-DoF four-bar model may not be able to predict the particular trajectory of the opercular mechanism, but it can predict its range of motion, and species may have sufficiently different ranges of motion to still draw important conclusions about functional diversity from morphology. Range of motion itself may prove to be a trait relevant to organismal performance.

The four-bar model, despite its wide range of possible conformations and joint configurations, is still one of the simplest linkage models in terms of number of links and connectivity. Thus, the four-bar represents a small subset of a diverse universe of potential mechanisms that can be tested against biomechanical systems. The incorporation of additional kinetic skeletal elements within fish heads will probably uncover new joint types and mechanism configurations. In addition, given the implied role of active motor control in multiple-DoF systems to determine particular motion trajectories through all the potential conformations, the integration of mechanism models with neuromuscular data may yield insights into the co-evolution between neuromotor control and musculoskeletal mobility. This study provides the basic framework essential for such endeavors and a broadly applicable approach for uncovering the wondrous and varied mechanisms within motor systems.

The authors thank Patricia Hernández, Sam Van Wassenbergh and Mark Westneat for helpful discussions.