Automatic control: the vertebral column of dogfish sharks behaves as a continuously variable transmission with smoothly shifting functions

Among vertebrates, the vertebral column of sharks functions in a novel way: during swimming, both intervertebral joints and vertebral centra undergo strain of similar magnitude (Porter et al., 2014). As the body bends, centra do too, storing and releasing mechanical power of the order of 10% of the shark's total metabolic power output. If this is true for cartilaginous fishes in general, then the loss of bone – which permits the centra to undergo large strains – may be an adaptation for enhanced spring-like behavior of the vertebral column. It is unlikely that the vertebral column behaves as a simple spring, as it is a composite of different materials – mineralized and unmineralized, tessellated and areolar cartilages – built into different structures (Dean and Summers, 2006). Given the functional novelty and adaptive importance of the shark vertebral column, our goal was to investigate its mechanical behavior during bending using an approach that is novel to comparative biomechanics: non-linear viscoelasticity.

Linear viscoelastic (LV) properties have been used to characterize the mechanical behavior of backbones – notochords and vertebral columns – as they undergo cyclic bending at curvatures and frequencies that approximate those observed as the fish swim (Long, 1992, 1995, 2002, 2011a). Unlike simpler static or quasi-static tests, cyclic tests vary time-dependent inputs in order to characterize dynamic time-dependent viscoelastic behavior (Summers and Long, 2006). For a given combination of cyclic strain and strain rate, the structure's mechanical behavior can be decomposed into independent viscous and elastic properties. The loss modulus, E″ (MPa), characterizes viscosity, which is proportional to the energy dissipated as heat. The storage modulus, E′ (MPa), characterizes elasticity, which is proportional to the mechanical work returned in recoil. Strictly, E″ and E′ are material properties measured from homogeneous structures; as biological structures are inherently heterogeneous, we caution that they represent ‘apparent’ material properties of the composite structure rather than material properties per se (Long et al., 2011a).

For sharks, LV properties of the bending vertebral column show functionally important variation across species and inputs. Foremost, as the amplitude of bending curvature, κ₀ (m⁻¹), increases, the apparent E′ increases in blacktip shark, Carcharinus limbatus, bonnethead shark, Sphryna tiburo, and dogfish shark, Squalus acanthias (Long et al., 2011a). While the apparent E″ does not vary with increasing κ₀ in C. limbatus and S. tiburo, it does in S. acanthias. As bending frequency, f (Hz), increases in tests done only on S. acanthias, E′ but not E″ increases. Mechanical properties and biochemical constituents (mineral, water, collagen, proteoglycan) of vertebral centra vary significantly among species (Porter et al., 2006, 2007). Any one of these differences in elastic mechanical properties may result in dramatic interspecific differences in the viscoelastic response as well. We can infer that the vertebral columns get stiffer (higher E′) without more energy loss (constant E″) as swimming speed, when controlled by tailbeat frequency alone, increases. If κ₀ increases with anatomical position, acceleration or turns, then the backbone stiffens (higher E′) and loses more energy (higher E″).

While it is illuminating to know that the mechanical behavior of the backbone changes depending on the locomotion of the shark, the LV approach makes an assumption that is demonstrably false: E′ and E″ are constant within a bending cycle across all values of κ₀. From the tests on sharks, we would expect at high amplitudes of κ₀ that the within-cycle E′ would increase, distorting the work loop of the bending moment, M (N m), and curvature, κ, from a symmetrical ellipse to one down-turned and up-turned at either end. In the extreme, within-cycle changes might create a ‘neutral zone’ where E′ is zero for low values of κ and then at the extremes provide a mechanical brake by a sudden increase in E′ (Nowroozi and Brainerd, 2014); quasi-static tests have detected neutral zones in the vertebral columns of longnose gar, Lepisosteus osseus (Long et al., 1996) and striped bass, Morone saxatalis (Nowroozi and Brainerd, 2013). If E″ increases with curvature, it could further alter the shape of the work loop, creating more energy loss at the extremes of the κ range. These predictions, with clear functional implications, can be tested using a non-linear viscoelastic (NLV) approach that explicitly accounts for changes in mechanical properties within a cycle.

In the NLV approach (Ewoldt et al., 2008a,b), linearity can be checked by viewing the raw work loops, by quantifying higher harmonics in the response, and by quantifying local viscoelastic properties throughout the beat cycle. As this NLV approach has been used to deform materials in shear, it has been dubbed large-amplitude oscillatory shear (LAOS). Here, by adjusting the underlying mathematics to account for the differences in strain, we adapted LOAS to large-amplitude oscillatory bending (LAOB). Using LAOB, mechanical properties of biological materials in bending can be described both in the classic LV cycle-averaged sense and in the new NLV within-cycle sense. We predict that this NLV approach will be a powerful tool for improving the measurement accuracy of the mechanical properties of shark backbones, herein, and any biological material more generally in bending. In particular, the assumption of linearity will be shown to be false for the material and curvature amplitudes considered here.

Fresh and frozen Squalus acanthias Linneaus (Squalidea: Squaliformes) vertebral columns were obtained from Mt Desert Island Biological Laboratory, Salisbury Cove, ME, USA. Animals were 77.2 cm (male), 74.0 cm (male) and 86.7 cm (female) total length. This species is usually demersal, sometimes found schooling up in the water column, and can make long seasonal migrations cued by water temperature (Compagno et al., 2005). We removed segments of 10 centra and nine intervertebral capsules from the precaudal region of the vertebral column in register with the anterior dorsal fin (Fig. 1A). Segment lengths (mm) and radii (mm) from the three individuals were as follows: 62.72 and 3.55; 52.64 and 2.95; and 53.12 and 3.3. We chose vertebrae from this region of the vertebral column because our previous work showed vertebral morphology from this region significantly predicted body curvature during routine turning (Porter et al., 2009). The neural arch cartilages were removed as previous mechanical testing showed that the neural arches do not bear significant load during compression (Porter and Long, 2010).

We used information about the swimming kinematics of S. acanthias to design the inputs for our mechanical tests (see below). Using S. acanthias individuals of similar size to those used in this study, Maia and Wilga (2016) trained the sharks to swim in a flow tank at speeds of 0.50 and 0.75 body lengths per second. At those two relatively slow speeds, individuals had an average tailbeat (bending) frequency, f, of 0.88 and 1.20 Hz. As the variance was not reported, we estimated standard deviation from videos we had of S. acanthias swimming (Porter et al., 2014). For four animals of similar size, swimming in a still-water tank, f ranged from 0.51 to 0.88 Hz, with an average f of 0.70±0.168 Hz (±1 s.d.). From these two sets of data we estimated a range of f during steady swimming of 0.5 Hz (lowest value measured) to 1.5 Hz (mean from Maia and Wilga, 2016; +2 s.d. from Porter et al., 2014). During unsteady escape maneuvers in smaller S. acanthias (mean total length of 58 cm), Domenici et al. (2004) measured half-cycle durations that, when doubled, yielded a range of transient f from 2.4 to 6.8 Hz. Using 2.0 Hz as a possible upper limit of f for steady swimming, we tested the vertebral columns over a range of steady cyclic f of 0.25–2.0 Hz. For maximal curvature amplitude, κ₀, we used a value, 5.0 m⁻¹, slightly above the mean value of 4.5 m⁻¹ reported for steadily swimming S. acanthias in the anatomical region from which we took our samples (region 3 in Porter et al., 2014).

Our goal was to measure the LV and NLV mechanics of the vertebral columns in bending. LV and NLV properties of a structure determine how it behaves mechanically. The term ‘behavior’ here indicates how a structure responds when external forces act upon it. In developing a reaction force by Newton's third law, structures change shape, a reconfiguration process that can be measured as strain, the ratio of the change in shape to the original shape. As reconfiguration involves movement of the structure, one may describe the system's dynamics, the interplay of force and motion, using Newtonian equations of motion. The total force can be partitioned into separate forces proportional to displacement (or strain), velocity (or strain rate) and acceleration: elastic, viscous and inertial forces, respectively (for review, see Den Hartog, 1934). In a linear equation of motion, each type of force has an associated coefficient that describes the magnitude of that force: stiffness, damping and mass. If the coefficients for stiffness and damping are constant with respect to changes in strain and strain rate, then the system is linearly viscoelastic (Lakes, 1998). If those coefficients change, then the system is non-linearly viscoelastic.

The key to measuring dynamic mechanical behavior is to apply a known force to a structure and measure the structure's change in shape. Sinusoidal motion is one type of common input, a cyclic loading that can occur in translation (tension and/or compression), torque or bending. We focused on cyclic bending as that is the primary motion of reconfiguration in vertebral columns of swimming fish. If a structure is bent within a two-dimensional plane, like a lateral bend in a swimming fish, then curvature, κ, describes the change in shape at any point along the column. Curvature, κ (m⁻¹), is the inverse of the radius of curvature, which is taken from the circle sized to fit perfectly along the bent midline at any point. Force, F (N), that causes the κ develops a moment, M (N m), through a lever or moment arm, r.

To cyclically bend vertebral columns, we designed a rig that would allow us to translate the uni-axial loading of our materials testing system (MTS, Tytron 250) into bending loads imposed by end-loading (Fig. 1B; Long et al., 2011b). We were able to change the amount of curvature in the vertebral column segment in two ways: by increasing the length of the moment arm or by increasing the actuator amplitude of the machine (Fig. 1C). Please see Long et al. (2011b) for a full description and validation of the bending rig.

Mechanical tests were conducted using a 50 N single-axis load cell at room temperature. Segments were bathed in 22°C elasmobranch Ringer’s solution (Forster et al., 1972) after each frequency sweep to maintain tissue hydration. A total of 45 mechanical tests were performed on each segment. We tested segments at five frequencies, f (0.25, 0.50, 1.00, 1.50 and 2.00 Hz), and three actuator displacement amplitudes, x₀ (5, 10 and 15 mm). Additionally, we adjusted the moment arm, r, of our rig to three lengths (57.15, 62.35 and 67.75 mm). The order of segment testing was identical. Each segment was tested at each r (order: 67.75, 62.35, 57.15 mm), each x₀ (order: 5, 10 and 15 mm) and each f (order: 0.25–2.00 Hz).

Oscillatory viscoelastic characterization is based on the idea that stress σ is a function of both strain and strain rate. For example, in shear σ is a function of both shear strain γ and shear strain rate ⁠. In oscillatory deformation, strain and strain rate are 90 deg out of phase, and the corresponding stress signal can be decomposed into elastic and viscous components. In this framework, intrinsic material properties can be calculated. For shear, linear viscoelastic moduli G are G′, the storage modulus, and G″, the loss modulus; for bending the analogous terms are E′, the storage modulus, and E″, the loss modulus; the tangent of the phase shift, δ, or loss angle, is the ratio of loss and storage moduli. Additional measures are required for non-linear viscoelasticity at larger strain amplitudes.

For composite structures with complex geometry, using intrinsic material properties may be inappropriate, as calculations of σ and strain require strong assumptions about spatially dependent material properties. In our case, we chose to characterize the composite structure response in terms of flexural stiffness, EI, which relates the measured bending moment M to the imposed bending curvature κ, defined for a purely elastic response as in Eqn 8. For a viscoelastic structure, the bending moment is a function of both the curvature κ and the rate of curvature ⁠. Based on this, we will characterize the composite in terms of both an elastic storage flexural stiffness (EI)′ and a viscous loss flexural stiffness (EI)″. EI will also serve as the framework for non-linear measures. In all cases, there is a direct conceptual mapping from G′ and G″ to (EI)′ and (EI)″.

The framework of (EI)′ and (ηI)′=(EI)″/ω can be extended to non-linear responses. NLV characterization is widely applied in oscillatory shear, where the term large-amplitude oscillatory shear (LAOS) is often used in the literature.

Here, we extend the ideas of LAOS to LAOB, i.e. large-amplitude oscillatory bending. By analogy, we can extend measures based on G′ and G″ to (EI)′ and (EI)″. In particular, all of the non-linear viscoelastic measures proposed by Ewoldt et al. (2008b) directly map from viscoelastic moduli to viscoelastic flexural stiffnesses. To complete the analogy, consider that shear stress σ is replaced by bending moment M, shear strain γ is replaced by curvature κ, and therefore G′ is replaced by (EI)′, and η′=G″/ω is replaced by (ηI)′=(EI)″/ω.

The raw oscillatory waveforms are conveniently viewed as parametric curves of bending as force (N) versus displacement (mm) (Fig. 2), in analogy with curves of stress versus strain for LAOS. LV responses appear as tilted ellipses in these coordinates. A non-linear response is non-elliptical.

The raw oscillatory data were quantitatively analyzed in Matlab using a custom-written program called MITlaos (Ewoldt et al., 2007, 2008a,b, 2009). The program calculates several quantitative descriptions in terms of viscoelastic bending stiffnesses, by analogy to the framework of LAOS (Ewoldt et al., 2008b). This includes (i) cycle-averaged first-harmonic measures, (ii) higher-harmonic components, and (iii) local measures of both elasticity and viscous dissipation at particular locations in the oscillatory cycle, e.g. at the extremes of maximum and minimum curvature. Our results will outline how these measures depend on both ω and κ_0.

Higher harmonics describe non-elliptical distortion of waveforms. Third-harmonics can be used to physically interpret the elastic stiffening/softening and viscous thinning/thickening throughout a particular oscillatory cycle (Ewoldt et al., 2008b). We will report the third-harmonics related to elasticity and viscosity, e₃ and ν₃, respectively.

For these measures, the subscript M refers to minimal strain (or strain rate) and L to large strain (or strain rate). The ratios of these measures indicate the degree of non-linearity within a cycle. The degree of elastic non-linearity can be described by the ratio of the elastic modulus at largest strain and the elastic modulus at minimum strain, ⁠, which we call the elastic index of non-linearity. The viscous index of non-linearity is similarly calculated as ⁠.

The results from each modulus and ratio were analyzed as the dependent responses to changes in the continuous main effects of f and κ₀, and the f×κ₀ interaction term in a mixed model ANCOVA with length of the segment, L, as covariate (JMP version 5.0.1a, SAS Institute; Zar, 1999).

NLV characterization generates a high-dimensional set of data, which we organize as follows. First, representative oscillatory waveforms are shown in detail individually and across the two-dimensional domain {f,κ₀} (Fig. 2). Waveforms are then quantified in several ways. The most basic measures are cycle-averaged (first-harmonic) elasticity and viscosity, used in Figs 3–5. Advanced measures then capture non-linearity within a single beat cycle, using local or regional measures rather than whole-cycle measurements (Figs 3,4,6), and higher harmonics (Fig. 7). An additional perspective in terms of bending work and recovered work is given in Fig. 8. Finally, Fig. 9 summarizes the most important results as response surfaces dependent on {f,κ₀}.

A single test at fixed frequency and amplitude is shown as a parametric plot of force and displacement (Fig. 2). The enclosed area indicates dissipated energy. Force and displacement directly map to moment and curvature such that slopes have dimensions proportional to [M/κ₀] and at certain locations can be interpreted as stiffness (N m²). Note that the first-harmonic (average) quantity will always report a smaller bending stiffness than the maximum stiffness experienced during the cycle. For the data in Fig. 2, the maximum stiffness occurs at large instantaneous curvature and (Fig. 2D).

Asymmetry is apparent at large instantaneous curvature in Fig. 2. That is, differences are seen in the magnitudes of M(t) between positive and negative curvature. This oscillatory asymmetry is typically not considered with shear characterization, even in LAOS. Most materials are shear-symmetric (Hyun et al., 2011), although methods exist for quantifying shear asymmetry and this involves even-harmonics (Yu et al., 2009). In our study, bending asymmetry may arise from the inherent biological geometry, or material distribution, or imperfect vertebral column alignment with respect to the neutral bending axis of the test instrument. We consider the observed asymmetry to be a sub-dominant feature, and therefore the quantitative measures presented in all figures will neglect asymmetry. That is, we will focus on the already complex symmetric aspects of the responses.

Representative force and displacement data for multiple tests are shown (Fig. 2), in the two-dimensional space of {f,κ₀} following Pipkin (1972). Each parametric plot is force (N) and displacement (mm), as in Fig. 10. The results are from a 10 centra long segment of vertebral column from S. acanthias. One can qualitatively interpret the response curves even before reporting quantitative descriptions of viscoelastic behavior (Fig. 10). Each curve shows one cycle of the test. Curves are more elliptically shaped at lower frequencies and amplitudes, although they are not precisely elliptical (and therefore not precisely LV). The elliptical deviation is pronounced in tests at higher frequencies and amplitudes, indicating that the mechanical behavior of the vertebral column is non-linear under these circumstances. In general, the curves show a predominantly elastic response, as there is little dissipative hysteresis. With limited dissipation, the local slope at any location of force versus displacement approximately corresponds to local bending elasticity. For the non-linear (non-elliptical curves), the slope tends to be large at the maximum instantaneous bending curvature. Thus, the maximum instantaneous (tangent) stiffness is experienced at maximum curvature within the cycle.

It was verified that compressive loads are negligible (Eqn 12) and do not contribute to the interpretation of LAOB results. Tests were conducted with different moment arm length r to verify that measured storage, ⁠, and loss stiffness, ⁠, were independent of the test geometry (Fig. S1). Non-linearity was also assessed by considering the dependence of these properties on curvature amplitude ⁠. A mixed model ANCOVA was applied, including r and MTS amplitude x₀, as well as the r×x₀ interaction. Both and changed significantly (P<0.0001) in the ANCOVA models (F_3,14=44.61 and F_3,14=26.87, respectively; Fig. S1). Amplitude x₀ was a significant main effect (P<0.0001), indicating non-linearity, but there was no significant effect of r or an x₀×r interaction, indicating that true bending properties are being measured and validating our assumption that compressive loads are negligible in these tests.

We found significant differences in bending elasticity ⁠, and for multiple frequencies f and amplitudes κ₀ (Fig. 3; Table S1). Mixed model ANCOVA for the linear flexural stiffness, ⁠, flexural stiffness at the largest strain, ⁠, and flexural stiffness at minimum strain, ⁠, are all statistically significant (P<0.0001; Fig. 3). Frequency, f, and curvature amplitude, κ₀, are significant main effects increasing (Fig. 3A). Additionally, for the linear flexural stiffness, ⁠, there was a significant interaction f×κ₀. Flexural stiffness at the largest strain, ⁠, increases significantly with increasing f, κ₀ and the interaction f×κ₀ (Fig. 3B). Flexural stiffness at minimum strain, ⁠, increases with increasing f and κ₀; however, the f×κ₀ interaction is not significant (Fig. 3C). Flexural stiffness values are largest at high f and κ₀. Vertebral column segment length was a significant covariate for these measures of bending elasticity.

Viscous effects are quantified in terms of ⁠, and (Fig. 4; Table S1). The first harmonic viscous modulus, ⁠, the viscous modulus at the largest strain, ⁠, and the viscous modulus at minimum strain, ⁠, were all statistically significant in regression models (P<0.0001; Fig. 4). Both f and κ₀ were significant model effects as well as the f×κ₀ interaction influencing and (Fig. 4). Both and are largest at the lowest f, while is an order of magnitude lower than those moduli. Vertebral column segment length was a significant covariate for these measures of viscous moduli.

Both elastic energy storage (stiffness) and viscous energy dissipation (viscosity) increase with increasing curvature amplitude (Figs 3 and 4; Table S1). But which property is increasing faster? The answer lies in the ratio and its dependence on κ₀. While dissipation increases with κ₀ (Fig. 5A), the ratio decreases with κ₀ (Fig. 5B). Thus, although dissipation is increasing, the stiffness increases comparatively more at large curvature amplitude, κ₀. The average (first-harmonic) elastic loss modulus and loss angle both have significant ANCOVA models (P<0.0001; Table S1; Fig. 5). increases in response to increasing κ₀ only (Fig. 5A). However, both f and κ₀ significantly impact tan δ₁ (Fig. 5B). Vertebral column segment length was a significant covariate for but not for tanδ₁.

Non-linearity within a beat cycle can be shown by comparing local measures as and (Fig. 6; Table S1). A ratio equal to 1 is the limit of a LV response. Flexural stiffness tends to increase at large instantaneous curvature, ⁠, indicating stiffening within a cycle. The elastic index of non-linearity has a significant ANCOVA model. Both f and κ₀ are significant effects influencing (Fig. 6A). The flexural dissipation is always non-linear within a cycle, for the available tests, tending to be more prominent at small instantaneous deformation rates, ⁠. The {f,κ₀} dependence of the viscous index of non-linearity was not significant (P<0.0001; Fig. 6). Vertebral column segment length was a significant covariate for but not for ⁠.

Higher harmonics are also used to quantify non-linearity within a cycle; Fig. 7 shows third-harmonic coefficients of the non-linear response (Ewoldt et al., 2008b). A value of 0 is the limit of a linear viscoelastic response. The initial deviation from linearity shows elastic stiffening e₃>0, and viscous thinning ν₃<0, consistent with local measures (Fig. 7). The elastic, e₃, and viscous, ν₃, intra-cycle non-linearity measures change significantly in the ANCOVA models (Table S1; Fig. 7). For e₃ the only significant effect was κ₀ (Fig. 7A). Both f and κ₀ were significant main effects for ν₃, but there was not a significant f×κ₀ interaction effect (Fig. 7B). Vertebral column segment length was a significant covariate for e₃ but not for ν₃.

Vertebral column work to bend, W_b, increases significantly with both f and κ₀ (Table S1; Fig. 8A). There is also a significant f×κ₀ interaction effect. Mechanical work recovered by recoil, W_r, also increases significantly with both f and κ₀ and there is a significant f×κ₀ interaction effect (Fig. 8B). Vertebral column segment length was a significant covariate.

The hypothesis that the vertebral column of sharks functions as a spring is too simple (Porter et al., 2014). Based on its mechanical behavior, the vertebral column may serve as both a spring and a brake. The transition between those two behaviors in an isolated vertebral column is continuous in the motor control space of tailbeat frequency, f, and bending curvature amplitude, κ₀ (Figs 9 and 10). We can reasonably speculate that to cruise faster with an elastic spring, the spiny dogfish, S. acanthias, would increase f and κ. These changes in bending motion will increase the vertebral column's flexural stiffness, ⁠, and reduce its flexural viscosity ⁠. More bending moment, M, can then be transmitted by the vertebral column as it bends and still more can be transmitted in recoil as it straightens. To slow down with an energy-dissipating brake, the shark would decrease f while maintaining a larger magnitude κ₀. These changes in bending motion will decrease the vertebral column's flexural stiffness, ⁠, and maintain its flexural viscosity ⁠. Kinetic energy is thus dissipated and bending motion is attenuated.

However, we need to proceed cautiously. While mechanical properties can address the capacities of the vertebral column in isolation, they do not necessarily determine its function during swimming. Any structure must be studied as part of the integrated, behaving animal to which it belongs. Ideally, vertebral column function should be studied in vivo. We have done initial work in vivo to determine loading in S. acanthias, implanting sonomicrometry crystals in the vertebral centra (Porter et al., 2014). This in vivo work is missing the direct measurement of the mechanical interactions of the vertebral column with the rest of the body. We would want to measure the strains within and loading of connective tissues attached to the vertebral column, the force activity patterns of muscles, metabolic rate of the whole shark, and the dynamic pressure field in the surrounding fluid. No one to our knowledge has attempted an experiment of this complexity in a living fish, probably because it would require, simultaneously and in a flow tank, surgical implantation of vast arrays of EMG electrodes, pressure transducers, strain gauges and sonomicrometry crystals; digital particle image velocimetry; and oxygen consumption measurements. New closed-loop physiological preparations in fish may offer an important intermediate approach between ideally complete and isolated-structure experiments (Tytell et al., 2014).

Without an appropriate model that assesses the mechanical behavior of the vertebral column in the integrated force-coupled context, we are left to speculate on its functional importance in swimming S. acanthias. Thinking about the considerations for isolated testing described above, we designed our experiments to mimic the f and κ of steady swimming (see Materials and methods). In addition, the precaudal section of the vertebral column that we tested undergoes considerable bending during steady swimming as determined from external kinematics (Porter et al., 2014). That same study showed that the vertebral column experiences mechanical loads in vivo sufficient to cause strains in the intervertebral joints and centra. But we do not know the nature of the three-dimensional stress fields that cause those strains, nor do we know how those stress fields cause strains in other tissues of the body. We can make some rough estimates regarding vertebral column energetics relative to the rest of the body.

Previously (Porter et al., 2014), we estimated from oxygen consumption measurements on S. acanthias (Brett and Blackburn, 1978) that the vertebral column returns about 12% of the metabolic power cost of swimming as elastic recoil. Another way to assess energetics is by using a hydrodynamic approach. While hydrodynamic power has not been measured or modeled in S. acanthias, Tytell and Lauder (2004) did both in eel, Anguilla rostrata. Both dogfish and eel generate undulations along the whole length of their bodies, allowing for elongated-body theory estimates to be used for comparison. In elongated-body theory, total hydrodynamic power is proportional to virtual mass, a property of the momentum transfer done by the body that is, in turn, proportional to the square of the dorso-ventral span of the caudal fin (Lauder and Tytell, 2006). As neither eels nor dogfish have simple trailing edges, we made rough estimates of 2 cm for the 20 cm long eels (Tytell and Lauder, 2004) and 19 cm for the 71 cm long dogfish of this study using a span to length ratio of 0.26 (Flammang, 2010). The scaling factor of virtual mass is thus 89, by which we multiply the eel's total hydrodynamic power, 191 µW, at a speed of 20 cm s⁻¹, to get a rough estimate of 17 µW in dogfish. The mechanical work of elastic recoil (Fig. 8, using a value of 130 µJ at f=1.0 Hz and κ₀=4 m⁻¹) is then multiplied by a factor of 11.4 to scale the small section of vertebral column to its full length, assuming that the vertebral column is 90% of the total length of the shark and invariant in its shape or properties. At 1.0 Hz, this yields a total recoil power of the vertebral column of 1.5 µW, about 9% of the total mechanical power. This first approximation suggests that the vertebral column of S. acanthias might function as a spring. Keeping in mind that we lack sufficient understanding of the whole-body mechanics to know the details of how this might work, we discuss the possible functional capacities of the vertebral column made possible by its mechanical properties.

To our knowledge, this is the first time that within-cycle NLV properties have been measured during cyclic loading in the vertebral columns of any species of fish. Thus, with LOAB, we are in a position to begin to complicate models of how vertebral columns might function (Fig. 10A) and to probe actual and possible spaces of mechanical behavior (Fig. 10B). The NLV approach shows that, strikingly, flexural stiffness increases with increasing κ within a cycle, as determined by the ratio of (EI)′ at the largest and minimum strains being greater than 1 (Fig. 6A). The flexural viscosity changes within a cycle but in the opposite direction with respect to ⁠, as determined by the ratio of (ηI)′ at the largest and minimum strains being less than 1 (Fig. 6B). These within-cycle changes are assigned by LOAB to properties of higher odd harmonics beyond the first harmonic properties of LV, with the third-harmonic coefficients, e₃ and ν₃ (Fig. 7).

Instead of thinking of vertebral columns as springs or brakes, they are both, and may be thought of as a type of continuously variable transmission. By avoiding the shifting of gears, mechanical continuously variable transmissions provide continuous and variable power flow over a range of speeds by adjusting the physical engagement of pulleys or hydraulic displacements (Singh and Nair, 1992). Physical engagement is adjusted in the vertebral column of S. acanthias by varying the relative strain occurring in intervertebral joints and vertebral centra (Porter et al., 2014). Remarkably, both joints and centra undergo significant strain during swimming and in ex vivo tests. This complicated mechanical behavior is matched by the vertebral column's complex structure, with the joints composed of fibrous elastic capsules containing fluid and the centra composed of calcified cartilage in biconic structures.

We expect that previous studies examining only the first-harmonic LV dynamic mechanical properties [ and ] may be underestimating the maximum stiffness of biological materials during cyclic loading when using large inputs. We found that the storage modulus and the elastic modulus at the largest strain were statistically similar (P<0.001; Fig. 3). By only taking the storage modulus into account and ignoring the modulus at minimum strain ⁠, we are ignoring potential variability in stiffness throughout the loading cycle. Previous research that examines exclusively first-harmonic mechanical properties may have underestimated, if not completely overlooked, some of the complexities and variability of vertebral column mechanics and the role of the vertebral column during swimming. Using the LAOB framework to determine and (Fig. 3) allows us to examine the mechanical variability of vertebral column stiffness. For example, increases as both testing f and κ₀ increase; however, only increasing κ₀ significantly increases (Table S1). Thus, is not frequency dependent in this testing regime, and the stiffness of the vertebral column at minimum strain will remain the same regardless of tailbeat frequency, but will change with body curvature.

In these experiments, we varied the mechanical loading on each specimen by varying f and κ₀. The length of the moment arm, r, the length of the specimen, L, and the amplitude of the MTS, x₀, are all variables contributing to κ₀. As the bending tested here is likely the result of pure bending and a superimposed and additive compressive load, we expected that increasing moment arm length would decrease the compressive loading on the specimen during testing. However, we found that changing the length of the moment arm does not significantly impact either or ⁠, but rather changing the MTS x₀ significantly increases those mechanical properties (Fig. S1).

We also acknowledge that there may be some additional frictional forces measured here that are not included in these stress calculations. We find that our force deformation curves have a parallelogram shape at small deformations, for example at f of 0.25–0.5 Hz and 5 mm x₀ (Fig. 2). We hypothesized that the pivot in the bending rig was causing potential frictional forces during these experiments. The material itself may actually cause these frictional forces. Hearle et al. (1969) suggested force displacement curves of this shape are the result of internal friction between fibers.

Passive mechanical non-linearity has been observed in many systems. Strain stiffening is the most common mechanical non-linearity, as observed in many tissues (Fung, 2010) such as artery walls (Shadwick, 1999) and the cornea (Hjortdal, 1995), as well as biological gels (Shah and Janmey, 1997; Storm et al., 2005; Ewoldt et al., 2009, 2011; Yao et al., 2013). Biomolecules themselves are inherently non-linear, including DNA (Perkins et al., 1995) and proteins (Rief et al., 1997; Kellermayer et al., 1997). Strain stiffening has been associated with both inherent molecular non-linearities and non-linear geometric effects. Inherent molecular non-linearities arise from non-linear force–extension curves, e.g. of cytoskeletal and extracellular proteins and filaments (Gardel et al., 2004; Storm et al., 2005; Lin et al., 2010). Non-linear geometric effects include buckling and non-affine deformation (Head et al., 2003; Onck et al., 2005; Chandran and Barocas, 2006; Heussinger and Frey, 2006).

In an attempt to circumscribe the possible mechanical functions of the vertebral column of sharks during swimming, we performed dynamic oscillatory bending tests on isolated sections of the vertebral column of spiny dogfish sharks, S. acanthias. Conducted over a range of physiologically realistic values of tailbeat frequency and midline curvature, these tests were analyzed using a novel technique for vertebral columns: non-linear viscoelasticity. The non-linear mechanical behavior detected includes changes in dissipative and elastic properties as frequency changes between cycles and as curvature changes within a cycle. These mechanical behaviors transition across frequency and curvature in a way that suggests that the vertebral column may function primarily as an elastic spring at high frequencies and as a brake at low frequencies.

Stephen Kajiura graciously provided S. acanthias vertebral columns for this experiment. Josh de Leeuw and Nicole Doorly wrote the C script to batch convert the force displacement data from the MTS into stress and strain. Two anonymous reviewers provided important feedback on this manuscript.