During swimming in dogfish sharks, Squalus acanthias, both the intervertebral joints and the vertebral centra undergo significant strain. To investigate this system, unique among vertebrates, we cyclically bent isolated segments of 10 vertebrae and nine joints. For the first time in the biomechanics of fish vertebral columns, we simultaneously characterized non-linear elasticity and viscosity throughout the bending oscillation, extending recently proposed techniques for large-amplitude oscillatory shear (LAOS) characterization to large-amplitude oscillatory bending (LAOB). The vertebral column segments behave as non-linear viscoelastic springs. Elastic properties dominate for all frequencies and curvatures tested, increasing as either variable increases. Non-linearities within a bending cycle are most in evidence at the highest frequency, 2.0 Hz, and curvature, 5 m−1. Viscous bending properties are greatest at low frequencies and high curvatures, with non-linear effects occurring at all frequencies and curvatures. The range of mechanical behaviors includes that of springs and brakes, with smooth transitions between them that allow for continuously variable power transmission by the vertebral column to assist in the mechanics of undulatory propulsion.

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, κ0 (m−1), 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 κ0 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 κ0 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 κ0. From the tests on sharks, we would expect at high amplitudes of κ0 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.

List of symbols and abbreviations
     
  • A

    cross-sectional area (m2)

  •  
  • e

    elastic Chebyshev coefficient (MPa)

  •  
  • E

    Young's modulus (MPa)

  •  
  • E

    linear viscoelastic storage modulus (MPa)

  •  
  • E

    linear viscoelastic loss modulus (MPa)

  •  
  • (EI)

    composite flexural stiffness (N m2)

  •  
  • (EI)′

    elastic storage flexural stiffness (N m2)

  •  
  • first-harmonic flexural storage stiffness (N m2)

  •  
  • local flexural stiffness at the largest strain (N m2)

  •  
  • local flexural stiffness at minimum strain (N m2)

  •  
  • (EI)″

    flexural loss (N m2)

  •  
  • first-harmonic flexural loss modulus (N m2)

  •  
  • f

    bending frequency

  •  
  • F

    force (N)

  •  
  • g

    gauge length (m)

  •  
  • G

    linear viscoelastic modulus (MPa)

  •  
  • G

    linear viscoelastic storage in shear (MPa)

  •  
  • G

    linear viscoelastic loss in shear (MPa)

  •  
  • I

    second moment of area (m4)

  •  
  • L

    specimen/segment length

  •  
  • LAOB

    large amplitude oscillatory bending

  •  
  • LAOS

    large amplitude oscillatory shear

  •  
  • LV

    linear viscoelastic

  •  
  • M

    moment (N m)

  •  
  • MTS

    materials testing system

  •  
  • N

    axial tension (N)

  •  
  • NLV

    non-linear viscoelastic

  •  
  • P

    pivot point

  •  
  • Q

    grip point

  •  
  • r

    moment arm length (m)

  •  
  • R

    resilience (%)

  •  
  • Rc

    radius of curvature (m)

  •  
  • S

    pivot point

  •  
  • T

    torque around the pivot point (N m)

  •  
  • v

    viscous Chebyshev coefficient (MPa s)

  •  
  • V

    vertical force (N)

  •  
  • Wb

    work to bend (µJ)

  •  
  • Wr

    work to recoil (µJ)

  •  
  • x0

    displacement amplitude

  •  
  • α

    moment arm angle during bending (rad)

  •  
  • γ

    shear strain (%)

  •  
  • δ

    phase shift due to viscous forces

  •  
  • η

    viscosity (N m2 s)

  •  
  • first harmonic flexural viscosity (N m2 s)

  •  
  • local flexural viscosity at the largest strain rate (N m2 s)

  •  
  • local flexural viscosity at minimum strain rate (N m2 s)

  •  
  • κ

    curvature (m−1)

  •  
  • κ0

    curvature amplitude (m−1)

  •  
  • rate of curvature (m−1 s−1)

  •  
  • σ

    stress (N m−2)

  •  
  • ω

    angular frequency (rad s−1)

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.

Study organisms

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).

Fig. 1.

Study samples, custom bending rig and free-body diagram. (A) Vertebral column sections used for testing consisted of 10 vertebrae and nine intervertebral joints taken from the precaudal region of the body just below the first dorsal fin. Neural arches were removed for testing. (B) Custom-made bending rig showing the calibration sample (white bar). The linear actuator's motion, delivered by an MTS Tytron 250, was converted into a bending moment by virtue of the moment arm (green). Vertebral column samples were kept hydrated with elasmobranch Ringer’s solution at 22°C. (C) Free-body diagrams of the bending apparatus full system and (D) a single moment arm and end of specimen. The displacement of pivot point S is controlled, and the force F=Px is measured. Positive force and displacement are defined by the arrow directions. P and S, pivot points; Q and R, specimen attachment points; T, torque around the pivot point; M, moment; N, axial tension; V, vertical force; r moment arm length; α, moment arm angle during bending.

Fig. 1.

Study samples, custom bending rig and free-body diagram. (A) Vertebral column sections used for testing consisted of 10 vertebrae and nine intervertebral joints taken from the precaudal region of the body just below the first dorsal fin. Neural arches were removed for testing. (B) Custom-made bending rig showing the calibration sample (white bar). The linear actuator's motion, delivered by an MTS Tytron 250, was converted into a bending moment by virtue of the moment arm (green). Vertebral column samples were kept hydrated with elasmobranch Ringer’s solution at 22°C. (C) Free-body diagrams of the bending apparatus full system and (D) a single moment arm and end of specimen. The displacement of pivot point S is controlled, and the force F=Px is measured. Positive force and displacement are defined by the arrow directions. P and S, pivot points; Q and R, specimen attachment points; T, torque around the pivot point; M, moment; N, axial tension; V, vertical force; r moment arm length; α, moment arm angle during bending.

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, κ0, we used a value, 5.0 m−1, slightly above the mean value of 4.5 m−1 reported for steadily swimming S. acanthias in the anatomical region from which we took our samples (region 3 in Porter et al., 2014).

Mechanical testing

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−1), 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, x0 (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 x0 (order: 5, 10 and 15 mm) and each f (order: 0.25–2.00 Hz).

Calculation of bending moment and curvature amplitude

The mechanical testing instrument imposes a known linear displacement x and measures a force F=Px (where P is the pivot point; Fig. 1C). This transfers the linear F to bending-dominated deformation by placing the sample at a distance r away from the line of force, connected with pivoted arms. The equations to calculate M and curvature amplitude, κ0, from force and displacement are derived here, with careful attention to axial compression and friction in the pivots. We show that a sufficiently large r is required for bending to dominate, and the apparatus is designed so that axial deformation is negligible and bending deformation is dominant in the resulting kinematics. The dominant role of bending exists under certain geometric conditions of the sample cross-section and moment arm of the instrument. The calculations in this section will be used to show that bending dominates when:
(1)
where r is the moment arm, A is the cross-sectional area and I is the second moment of area of the specimen.
The quasistatic loading balance for the entire setup gives the following set of equations:
(2)
where g is the gauge length (neutral bending axis length) of the specimen. The torques TP and TS originate from the possibility of friction at the pivot points P and S. If the friction is equal with TP=TS, or if we assume an idealized frictionless pivot with TP=TS=0, then the vertical forces will become zero, Py=Sy=0. For our analysis, we will assume that Py=Sy=0.
To find the load applied to the specimen, consider the free-body diagram (Fig. 1D). For small displacements, a moment balance around the specimen attachment point Q shows that M=FrTP, where M is the bending moment at point Q, r is the moment arm, F is the force measured by the instrument and TP is the torque at pivot point P. The torque TP is generally unknown, and in order to make calculations, we will assume an idealized frictionless pivot with TP=0, in which case we estimate the bending moment as:
(3)
In the presence of pivot friction, the actual applied moment will be less than the value calculated by Eqn 3. The bending moment will be constant across the length of the specimen only if the vertical force is zero, Py=0, which we will assume. The shear force in the specimen is also zero if the vertical loads are zero, V=Py=0. The axial tension in the specimen is given by N=−F.
Curvature of the specimen is induced by moving the pivot point S a distance x, which results in each moment arm turning through an angle α as shown in Fig. 1. The angle is related to the displacement and moment arm by:
(4)
The curvature of the specimen is determined by the relative angle of the two moment arms. Two imaginary lines along the lengths of the moment arms would intersect at an angle 2α. This angle is defined by the arc length g (at the neutral axis) and the radius of curvature Rc:
(5)
We combine Eqns 4 and 5 to get a single equation for the curvature κ=1/Rc:
(6)
For small displacements, x<<2r, the arcsine function simplifies and the curvature is linearly proportional to the instrument displacement:
(7)
With the expressions for bending moment, Eqn 3, and curvature amplitude, Eqn 6, we can define the apparent flexural stiffness EI with the expression:
(8)
For a homogeneous material, E is the Young's modulus and I is the second moment of area about the neutral bending axis. For our composite structure, we will simply use the composite flexural stiffness (EI) to describe the response, and extend this characterization measure to include non-linear viscoelastic effects with oscillatory input κ(t)=κ0sin(ωt) and measured output M(t0,ω).
To validate the assumption of negligible axial compression, the results given above can be used to compare the bending and axial contributions to the total displacement x=xbend+xaxial. For small displacements, the bending contribution curvature is related to the bending moment M=(EI)κ by classic beam theory. Using this with Eqn 7:
(9)
The moment M=Fr is used from Eqn 3 to give the result:
(10)
The axial deformation contribution to the displacement x is given by uniaxial compression:
(11)
where A is the load-bearing cross-sectional area. For bending to dominate, we desire xbend/xaxial>>1. The results of Eqns 10 and 11 give the criteria that:
(12)
for bending to dominate, as noted in Eqn 1. This depends only on the specimen cross-section geometry and moment arm of the instrument. As r→0, there is no bending, but with a sufficiently long moment arm, all x displacement goes to bending. For our tests, the range of r is great enough so that bending dominates; the average value for the experiments here is Ar2/I≈1500>>1.

Mechanical analysis

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)″.

For our displacement controlled tests, κ oscillates as a sinusoid. We note that although displacement x(t) is explicitly controlled, which causes κ(t), we have shown that κ directly follows the imposed displacement x(t) for our conditions Ar2/I≈1500>>1. For steady oscillations, a reference time is arbitrary except when calculating signs of higher harmonics in the response (Ewoldt, 2013). We choose to represent the oscillating curvature amplitude as:
(13)
where κ­0 is the amplitude of curvature and ω is the angular frequency. The measured bending moment in the linear viscoelastic regime is then:
(14)
which is phase shifted by an angle δ due to viscoelasticity. In terms of viscoelastic flexural stiffness, this can be re-written as:
(15)
where
(16)
and
(17)
where η is the viscosity.

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.

Fig. 2.

Mechanical behavior of the dynamic bending system. (A) Work loop of five cycles for ultra high-weight polyethylene (NIST reference material 84567), bent at a frequency of 1.0 Hz and with a machine displacement amplitude of ±2.0 mm. Gage length was 114.5 mm, width 10.3 mm, thickness 3.7 mm. (B) The same five cycles over time, with force as the black line and displacement of the actuator as the green line. (C) Close up of the fifth cycle, showing the crossing of the x-axis. The difference in crossing times is the phase lag, δ, in seconds. (D) Work loop of five cycles for a shark vertebral column motion segment, bent a frequency of 1.0 Hz and with a machine displacement amplitude of ±5.0 mm. Non-linearities are clearly present in the deviations from a simple ellipsoidal shape. Yellow line indicates the slope at minimum strain, which is proportional to (EI)′M. The blue line represents the slope at largest strain, which is proportional to (EI)′L. (E) The same five cycles over time.

Fig. 2.

Mechanical behavior of the dynamic bending system. (A) Work loop of five cycles for ultra high-weight polyethylene (NIST reference material 84567), bent at a frequency of 1.0 Hz and with a machine displacement amplitude of ±2.0 mm. Gage length was 114.5 mm, width 10.3 mm, thickness 3.7 mm. (B) The same five cycles over time, with force as the black line and displacement of the actuator as the green line. (C) Close up of the fifth cycle, showing the crossing of the x-axis. The difference in crossing times is the phase lag, δ, in seconds. (D) Work loop of five cycles for a shark vertebral column motion segment, bent a frequency of 1.0 Hz and with a machine displacement amplitude of ±5.0 mm. Non-linearities are clearly present in the deviations from a simple ellipsoidal shape. Yellow line indicates the slope at minimum strain, which is proportional to (EI)′M. The blue line represents the slope at largest strain, which is proportional to (EI)′L. (E) The same five cycles over time.

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.

First-harmonic measures are the most basic measures used to describe non-linear oscillatory responses and are typically used by commercial instrument software. First-harmonics represent average measures (Ewoldt et al., 2008b; Ewoldt, 2013), as they are integrated over an entire period of raw waveform data. For LAOB, we will use the cycle-averaged first-harmonic elasticity and dissipation . The first-harmonic loss is proportional to the total energy lost in a single cycle, whereas can be related to work stored in a cycle. We calculated mechanical work to bend (Long et al., 2011a), Wb (µJ), per cycle done on the system, where L (m) is the specimen length, κ0 (m−1), as:
(18)
The mechanical work recovered as elastic recoil, Wr (µJ), takes into account the work to bend, Wb, and resilience, R (%), where:
(19)
The first-harmonic phase represents a relative measure of dissipated and stored energy, .

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, e3 and ν3, respectively.

Local non-linear measures, a regional measure at minimum or largest strains rather than a whole-cycle measure, will give the clearest description of non-linear viscoelasticity in this study: local elastic bending stiffness at minimal curvature amplitude:
(20)
and largest curvature amplitude:
(21)
The complementary viscous measures are best interpreted with parametric curves of moment M(t) versus rate of curvature . Local dissipative measures are defined at locations of minimal rate of curvature:
(22)
and largest rate-of-curvature:
(23)

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 .

Statistical analysis

The results from each modulus and ratio were analyzed as the dependent responses to changes in the continuous main effects of f and κ0, and the f×κ0 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 {f0} (Fig. 2). Waveforms are then quantified in several ways. The most basic measures are cycle-averaged (first-harmonic) elasticity and viscosity, used in Figs 35. 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 {f0}.

Fig. 3.

Flexural stiffness of the bending vertebral column shows significant amplitude-dependent non-linearity over a range of frequencies, f. This amplitude-dependent non-linearity (P<0.0001; see Table S1) is shown by the cycle-averaged measure of flexural stiffness, , as well as local measures of flexural stiffness at the largest and minimum instantaneous strain, and . (A) First-harmonic increased with increasing curvature amplitude, κ0, and f. At higher f, the κ0 effect was amplified, as indicated by a significant κ0×f interaction. (B) also increased with increasing κ0 and f. At higher f, the κ0 effect was amplified, as indicated by a significant κ0×f interaction. (C) also increased with increasing κ0 and f. However, did not vary with the κ0×f interaction. For a given trial, the inequality of and indicates elastic non-linearity within a beat cycle, i.e. is stiffening at large curvature (see Fig. 6). These data are from 145 different mechanical tests and each point represents a mechanical test at each f and each displacement amplitude, x0, for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 3.

Flexural stiffness of the bending vertebral column shows significant amplitude-dependent non-linearity over a range of frequencies, f. This amplitude-dependent non-linearity (P<0.0001; see Table S1) is shown by the cycle-averaged measure of flexural stiffness, , as well as local measures of flexural stiffness at the largest and minimum instantaneous strain, and . (A) First-harmonic increased with increasing curvature amplitude, κ0, and f. At higher f, the κ0 effect was amplified, as indicated by a significant κ0×f interaction. (B) also increased with increasing κ0 and f. At higher f, the κ0 effect was amplified, as indicated by a significant κ0×f interaction. (C) also increased with increasing κ0 and f. However, did not vary with the κ0×f interaction. For a given trial, the inequality of and indicates elastic non-linearity within a beat cycle, i.e. is stiffening at large curvature (see Fig. 6). These data are from 145 different mechanical tests and each point represents a mechanical test at each f and each displacement amplitude, x0, for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 4.

Flexural viscosity of the bending vertebral column shows significant amplitude-dependent non-linearity over a range of frequencies, f. This amplitude-dependent non-linearity (P<0.0001; see Table S1) is shown by the cycle-averaged measure of flexural viscosity, , as well as local measures of flexural viscosity at large and minimum instantaneous strain rate, and . (A) First-harmonic , increased with increasing κ0 and decreased with increasing f. At higher f, the κ0 effect was attenuated, as indicated by a decreasing slope and a significant κ0×f interaction. (B) also increased with increasing κ0 and decreased with increasing f. At higher f, the κ0 effect was attenuated, as indicated by a decreasing slope and a significant κ0×f interaction. (C) also increased with increasing κ0 and decreased with increasing f. The inequality of and indicates viscous non-linearity within a beat cycle, i.e. is viscous thickening at large rates of curvature κ (see Fig. 7). These data are from 145 different mechanical tests and each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 4.

Flexural viscosity of the bending vertebral column shows significant amplitude-dependent non-linearity over a range of frequencies, f. This amplitude-dependent non-linearity (P<0.0001; see Table S1) is shown by the cycle-averaged measure of flexural viscosity, , as well as local measures of flexural viscosity at large and minimum instantaneous strain rate, and . (A) First-harmonic , increased with increasing κ0 and decreased with increasing f. At higher f, the κ0 effect was attenuated, as indicated by a decreasing slope and a significant κ0×f interaction. (B) also increased with increasing κ0 and decreased with increasing f. At higher f, the κ0 effect was attenuated, as indicated by a decreasing slope and a significant κ0×f interaction. (C) also increased with increasing κ0 and decreased with increasing f. The inequality of and indicates viscous non-linearity within a beat cycle, i.e. is viscous thickening at large rates of curvature κ (see Fig. 7). These data are from 145 different mechanical tests and each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 5.

Flexural loss and loss angle of the bending vertebral column show significant non-linearity. (A) Flexural loss, , increases significantly only with curvature, κ0 (P<0.0001; Table S1). (B) The relative loss ratio, , decreases significantly with both κ0 and frequency, f (P<0.0001; Table S1). Note that (where ω is angular frequency); thus, and (Fig. 4A) are not independent responses. Each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 5.

Flexural loss and loss angle of the bending vertebral column show significant non-linearity. (A) Flexural loss, , increases significantly only with curvature, κ0 (P<0.0001; Table S1). (B) The relative loss ratio, , decreases significantly with both κ0 and frequency, f (P<0.0001; Table S1). Note that (where ω is angular frequency); thus, and (Fig. 4A) are not independent responses. Each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 6.

Non-linearity within a beat cycle as quantified by ratios of local flexural stiffness and local flexural viscosity, for a range of amplitudes, κ0, and frequencies, f. A ratio of 1 is the limit of a linear viscoelastic response. (A) Flexural stiffness tends to increase at large instantaneous curvature, , indicating stiffening within a cycle. The ratio increased significantly with increasing κ0 (P=0.01) and f (P<0.0001; Table S1). (B) The flexural dissipation tends to be more prominent at small instantaneous deformation rates, . The viscous index of non-linearity did not change significantly with increasing κ0 and increasing f (Table S1). Each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 6.

Non-linearity within a beat cycle as quantified by ratios of local flexural stiffness and local flexural viscosity, for a range of amplitudes, κ0, and frequencies, f. A ratio of 1 is the limit of a linear viscoelastic response. (A) Flexural stiffness tends to increase at large instantaneous curvature, , indicating stiffening within a cycle. The ratio increased significantly with increasing κ0 (P=0.01) and f (P<0.0001; Table S1). (B) The flexural dissipation tends to be more prominent at small instantaneous deformation rates, . The viscous index of non-linearity did not change significantly with increasing κ0 and increasing f (Table S1). Each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 7.

Non-linearity within a beat cycle as quantified by third-harmonic coefficients of the oscillatory bending moment, over a range of amplitudes, κ0, and frequencies, f. A value of 0 is the limit of a linear viscoelastic response. The initial deviation from linearity shows (A) elastic stiffening e3>0 and (B) viscous thinning v3<0, consistent with local measures in Fig. 6. See Table S1 for quantitative significance regarding dependence on κ0, f, and the κ0×f interaction. Each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 7.

Non-linearity within a beat cycle as quantified by third-harmonic coefficients of the oscillatory bending moment, over a range of amplitudes, κ0, and frequencies, f. A value of 0 is the limit of a linear viscoelastic response. The initial deviation from linearity shows (A) elastic stiffening e3>0 and (B) viscous thinning v3<0, consistent with local measures in Fig. 6. See Table S1 for quantitative significance regarding dependence on κ0, f, and the κ0×f interaction. Each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 8.

Mechanical work, W, of cartilaginous vertebral columns changes with curvature, κ0, and frequency, f. See Table S1. (A) Work to bend, Wb, shows significant changes (P<0.0001) with changing inputs. Wb increases with increasing curvature amplitude, κ0, and frequency, f, with the affects of κ0 amplified at higher f, as indicated by a significant κ0×f interaction. (B) Mechanical work recovered as elastic recoil, Wr, also changes significantly (P<0.0001) over the range of curvature amplitudes, κ0, and frequencies, f, examined here. Each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 8.

Mechanical work, W, of cartilaginous vertebral columns changes with curvature, κ0, and frequency, f. See Table S1. (A) Work to bend, Wb, shows significant changes (P<0.0001) with changing inputs. Wb increases with increasing curvature amplitude, κ0, and frequency, f, with the affects of κ0 amplified at higher f, as indicated by a significant κ0×f interaction. (B) Mechanical work recovered as elastic recoil, Wr, also changes significantly (P<0.0001) over the range of curvature amplitudes, κ0, and frequencies, f, examined here. Each point represents a mechanical test at each f and each x0 for each specimen. Linear trend lines are fitted to these data at each frequency.

Fig. 9.

A continuously variable transmission as a function of frequency and curvature amplitude, {f0}, relevant to swimming. The first-harmonic average measures of elastic energy storage (A, bending elasticity) and viscous energy loss (B, bending viscosity) are normalized by their maximum value within the observed domain and color-coded according to the z-coordinate height. Data from Fig. 3A and Fig. 4A.

Fig. 9.

A continuously variable transmission as a function of frequency and curvature amplitude, {f0}, relevant to swimming. The first-harmonic average measures of elastic energy storage (A, bending elasticity) and viscous energy loss (B, bending viscosity) are normalized by their maximum value within the observed domain and color-coded according to the z-coordinate height. Data from Fig. 3A and Fig. 4A.

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 [M0] and at certain locations can be interpreted as stiffness (N m2). 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 {f0} 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 x0, as well as the r×x0 interaction. Both and changed significantly (P<0.0001) in the ANCOVA models (F3,14=44.61 and F3,14=26.87, respectively; Fig. S1). Amplitude x0 was a significant main effect (P<0.0001), indicating non-linearity, but there was no significant effect of r or an x0×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 κ0 (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, κ0, are significant main effects increasing (Fig. 3A). Additionally, for the linear flexural stiffness, , there was a significant interaction f×κ0. Flexural stiffness at the largest strain, , increases significantly with increasing f, κ0 and the interaction f×κ0 (Fig. 3B). Flexural stiffness at minimum strain, , increases with increasing f and κ0; however, the f×κ0 interaction is not significant (Fig. 3C). Flexural stiffness values are largest at high f and κ0. 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 κ0 were significant model effects as well as the f×κ0 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 κ0. While dissipation increases with κ0 (Fig. 5A), the ratio decreases with κ0 (Fig. 5B). Thus, although dissipation is increasing, the stiffness increases comparatively more at large curvature amplitude, κ0. 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 κ0 only (Fig. 5A). However, both f and κ0 significantly impact tan δ1 (Fig. 5B). Vertebral column segment length was a significant covariate for but not for tanδ1.

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 κ0 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 {f0} 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 e3>0, and viscous thinning ν3<0, consistent with local measures (Fig. 7). The elastic, e3, and viscous, ν3, intra-cycle non-linearity measures change significantly in the ANCOVA models (Table S1; Fig. 7). For e3 the only significant effect was κ0 (Fig. 7A). Both f and κ0 were significant main effects for ν3, but there was not a significant f×κ0 interaction effect (Fig. 7B). Vertebral column segment length was a significant covariate for e3 but not for ν3.

Vertebral column work to bend, Wb, increases significantly with both f and κ0 (Table S1; Fig. 8A). There is also a significant f×κ0 interaction effect. Mechanical work recovered by recoil, Wr, also increases significantly with both f and κ0 and there is a significant f×κ0 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, κ0 (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 κ0. 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.

Fig. 10.

Mechanical behavior of the vertebral column. (A) Non-linear viscoelastic responses (NLV; blue, raw data) compared with first-harmonic linear viscoelastic (LV) responses (gray, smooth ellipses) in Squalus acanthias. Each hysteresis curve is moment (M) versus curvature (κ) normalized by the maximum absolute value of the measured moment (vertical axis) and imposed curvature (horizontal axis). At low frequencies, the vertebral column functions as a brake or damper, dissipating energy (energy lost is proportional to the area circumscribed by the curve). At high frequencies, the vertebral column functions more as a spring, storing and releasing energy. For a given frequency, increasing curvature increases the NLV behavior. (B) Hypothetical mechanical behaviors, not seen in sharks, but possible by NLV. When the third-order coefficient of elasticity is positive, elastic stiffening may create a neutral zone. Thus, NLV properties offer an explanation for neutral zones of bending reported in some bony fish species. When the third-order coefficient of elasticity is negative, the opposite occurs, with elastic softening like that reported in blue marlin.

Fig. 10.

Mechanical behavior of the vertebral column. (A) Non-linear viscoelastic responses (NLV; blue, raw data) compared with first-harmonic linear viscoelastic (LV) responses (gray, smooth ellipses) in Squalus acanthias. Each hysteresis curve is moment (M) versus curvature (κ) normalized by the maximum absolute value of the measured moment (vertical axis) and imposed curvature (horizontal axis). At low frequencies, the vertebral column functions as a brake or damper, dissipating energy (energy lost is proportional to the area circumscribed by the curve). At high frequencies, the vertebral column functions more as a spring, storing and releasing energy. For a given frequency, increasing curvature increases the NLV behavior. (B) Hypothetical mechanical behaviors, not seen in sharks, but possible by NLV. When the third-order coefficient of elasticity is positive, elastic stiffening may create a neutral zone. Thus, NLV properties offer an explanation for neutral zones of bending reported in some bony fish species. When the third-order coefficient of elasticity is negative, the opposite occurs, with elastic softening like that reported in blue marlin.

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−1, 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 κ0=4 m−1) 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, e3 and ν3 (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 κ0 increase; however, only increasing κ0 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 κ0. The length of the moment arm, r, the length of the specimen, L, and the amplitude of the MTS, x0, are all variables contributing to κ0. 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 x0 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 x0 (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).

Conclusions

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.

Author contributions

M.E.P. and J.H.L. designed the experiments and collected data. All authors analyzed these data using an approach developed by R.H.E. All authors contributed to manuscript preparation and editing prior to submission.

Funding

This work was supported by the US National Science Foundation [DBI-0442269, IOS-0922605, and grant no. 1344227, INSPIRE, Special Projects].

Brett
,
J. R.
and
Blackburn
,
J. M.
(
1978
).
Metabolic rate and energy expenditure of the spiny dogfish, Squalus acanthias
.
J. Fish. Res. Board Can.
35
,
816
-
821
.
Chandran
,
P. L.
and
Barocas
,
V. H.
(
2006
).
Affine versus non-affine fibril kinematics in collagen networks: theoretical studies of network behavior
.
J. Biomech. Eng.
128
,
259
-
270
.
Compagno
,
L.
,
Dando
,
M.
and
Fowler
,
S.
(
2005
).
Sharks of the World
.
New Jersey
:
Princeton University Press
.
Dean
,
M. N.
and
Summers
,
A. P.
(
2006
).
Mineralized cartilage in the skeleton of chondrichthyan fishes
.
Zoology
109
,
164
-
168
.
Den Hartog
,
J. P.
(
1934
).
Mechanical Vibrations
.
New York
:
Dover Publications
.
Domenici
,
P.
,
Standen
,
E. M.
and
Levine
,
R. P.
(
2004
).
Escape manoeuvres in the spiny dogfish (Squalus acanthias)
.
J. Exp. Biol.
207
,
2339
-
2349
.
Ewoldt
,
R. H.
(
2013
).
Defining nonlinear rheological material functions for oscillatory shear
.
J. Rheol.
57
,
177
.
Ewoldt
,
R. H.
,
Winter
,
P.
and
McKinley
,
G. H.
(
2007
).
MITlaos version 2.1 Beta for MATLAB
.
Cambridge, MA
:
self-published
.
Ewoldt
,
R. H.
,
Hosoi
,
A. E.
and
McKinley
,
G. H.
(
2008a
).
An ontology for large amplitude oscillatory shear flow. The XV International Congress on Rheology: The Society of Rheology 80th Annual Meeting, Monterey (California)
.
AIP Conf. Proc.
1027
,
1135
.
Ewoldt
,
R. H.
,
Hosoi
,
A. E.
and
McKinley
,
G. H.
(
2008b
).
New measures for characterizing nonlinear viscoelasticity in large amplitude oscillatory shear
.
J. Rheol.
52
,
1427
.
Ewoldt
,
R. H.
,
Hosoi
,
A. E.
and
McKinley
,
G. H.
(
2009
).
Nonlinear viscoelastic biomaterials: meaningful characterization and engineering inspiration
.
Integr. Comp. Biol.
49
,
40
-
50
.
Ewoldt
,
R. H.
,
Winegard
,
T. M.
and
Fudge
,
D. S.
(
2011
).
Non-linear viscoelasticity of hagfish slime
.
Int. J. Non Linear Mech.
46
,
627
-
636
.
Flammang
,
B. E.
(
2010
).
Functional morphology of the radialis muscle in shark tails
.
J. Morphol.
271
,
340
-
352
.
Forster
,
R. P.
,
Goldstein
,
L.
and
Rosen
,
J. K.
(
1972
).
Intrarenal control of urea reabsorption by renal tubules of the marine elasmobranch, Squalus acanthias
.
Comp. Biochem. Physiol. A Physiol.
42
,
3
-
12
.
Fung
,
Y.-C.
(
2010
).
Biomechanics: Mechanical Properties of Living Tissues
.
New York
:
Springer
.
Gardel
,
M. L.
,
Shin
,
J. H.
,
MacKintosh
,
F. C.
,
Mahadevan
,
L.
,
Matsudaira
,
P.
and
Weitz
,
D. A.
(
2004
).
Elastic behavior of cross-linked and bundled actin networks
.
Science
304
,
1301
-
1305
.
Head
,
D.
,
Levine
,
A.
and
MacKintosh
,
F.
(
2003
).
Deformation of cross-linked semiflexible polymer networks
.
Phys. Rev. Lett.
91
,
108102
.
Hearle
,
J. W. S.
,
Grosberg
,
P.
and
Backer
,
S.
(
1969
).
Structural Mechanics of Fibers, Yarns, and Fabrics
.
New York
:
Wiley-Interscience
.
Heussinger
,
C.
and
Frey
,
E.
(
2006
).
Stiff polymers, foams, and fiber networks
.
Phys. Rev. Lett.
96
,
017802
.
Hjortdal
,
J. Ø.
(
1995
).
Extensibility of the normo-hydrated human cornea
.
Acta Ophthalmol. Scand.
73
,
12
-
17
.
Hyun
,
K.
,
Wilhelm
,
M.
,
Klein
,
C. O.
,
Cho
,
K. S.
,
Nam
,
J. G.
,
Ahn
,
K. H.
,
Lee
,
S. L.
,
Ewoldt
,
R. H.
and
McKinley
,
G. H.
(
2011
).
A review of nonlinear oscillatory shear tests: analysis and application of Large Amplitude Oscillatory Shear (LAOS)
.
Prog. Polym. Sci.
36
,
1697
-
1753
.
Kellermayer
,
M. S. Z.
,
Smith
,
S. B.
,
Granzier
,
H. L.
and
Bustamante
,
C.
(
1997
).
Folding-unfolding transitions in single titin molecules characterized with laser tweezers
.
Science
276
,
1112
-
1116
.
Lakes
,
R. S.
(
1998
).
Viscoelastic Solids
.
New York
:
CRC Press
.
Lauder
,
G. V.
and
Tytell
,
E. D.
(
2006
).
Hydrodynamics of undulatory propulsion
.
Fish Physiol.
23
,
425
.
Lin
,
Y.-C.
,
Yao
,
N. Y.
,
Broedersz
,
C. P.
,
Herrmann
,
H.
,
MacKintosh
,
F. C.
and
Weitz
,
D. A.
(
2010
).
Origins of elasticity in intermediate filament networks
.
Phys. Rev. Lett.
104
,
058101
.
Long
,
J. H.
, Jr
(
1992
).
Stiffness and damping forces in the intervertebral joints of blue marlin, (Makaira nigricans)
.
J. Exp. Biol.
162
,
131
-
155
.
Long
,
J. H.
, Jr
(
1995
).
Morphology, mechanics, and locomotion: the relation between the notochord and swimming motions in sturgeon
.
Environ. Biol. Fish.
44
,
199
-
211
.
Long
,
J. H.
, Jr,
,
Hale
,
M. E.
,
McHenry
,
M. J.
and
Westneat
,
M. W.
(
1996
).
Functions of fish skin: flexural stiffness and steady swimming of longnose gar Lepisosteus osseus
.
J. Exp. Biol.
199
,
2139
-
2151
.
Long
,
J. H.
, Jr
,
Koob-Emunds
,
M.
,
Sinwell
,
B.
and
Koob
,
T. J.
(
2002
).
The notochord of hagfish, Myxine glutinosa: viscoelastic properties and mechanical functions during steady swimming
.
J. Exp. Biol.
205
,
3819
-
3831
.
Long
,
J. H.
, Jr
,
Koob
,
T. J.
,
Irving
,
K.
,
Combie
,
K.
,
Engel
,
V.
,
Livingston
,
N.
,
Lammert
,
A.
and
J.
and
Schumacher
, (
2006
).
Biomimetic evolutionary analysis: testing the adaptive value of vertebrate tail stiffness in autonomous swimming robots
.
J. Exp. Biol.
209
,
4732
-
4746
.
Long
,
J. H.
, Jr
,
Koob
,
T. J.
,
Schaefer
,
J. T.
,
Summers
,
A. P.
,
Bantilan
,
K.
,
Grotmol
,
S.
and
Porter
,
M. E.
(
2011a
).
Inspired by sharks: a biomimetic skeleton for the flapping, propulsive tail of an aquatic robot
.
Mar. Technol. Soc. J.
45
,
119
-
129
.
Long
,
J. H.
, Jr
,
Krenitsky
,
N.
,
Roberts
,
S. F.
,
Hirokawa
,
J.
,
de Leeuw
,
J.
and
Porter
,
M. E.
(
2011b
).
Testing biomimetic structures in bioinspired robots: how vertebrae control the stiffness of the body and the behavior of fish-like swimmers
.
Integr. Comp. Biol.
51
,
158
-
175
.
Maia
,
A.
and
Wilga
,
C. A.
(
2016
).
Dorsal fin function in spiny dogfish during steady swimming
.
J. Zool.
298
,
139
-
149
.
Nowroozi
,
B. N.
and
Brainerd
,
E. L.
(
2013
).
X-ray motion analysis of the vertebral column during the startle response in striped bass, Morone saxatilis
.
J. Exp. Biol.
216
,
2833
-
2842
.
Nowroozi
,
B. N.
and
Brainerd
,
E. L.
(
2014
).
Importance of mechanics and kinematics in determining the stiffness contribution of the vertebral column during body-caudal-fin swimming in fishes
.
Zoology
117
,
28
-
35
.
Onck
,
P.
,
Koeman
,
T.
,
van Dillen
,
T.
and
van der Giessen
,
E.
(
2005
).
Alternative explanation of stiffening in cross-linked semiflexible networks
.
Phys. Rev. Lett.
95
,
178102
.
Perkins
,
T. T.
,
Smith
,
D. E.
,
Larson
,
R. G.
and
Chu
,
S.
(
1995
).
Stretching of a single tethered polymer in a uniform flow
.
Science
268
,
83
-
87
.
Pipkin
,
A. C.
(
1972
).
Lectures on Viscoelasticity Theory
, Vol.
7
.
New York
:
Springer
.
Porter
,
M. E.
and
Long
,
J. H.
, Jr
(
2010
).
Vertebrae in compression: mechanical behavior of arches and centra in the gray smooth-hound (Mustelus californicus)
.
J. Morphol.
271
,
366
-
375
.
Porter
,
M. E.
,
Beltrán
,
J. L.
,
Koob
,
T. J.
and
Summers
,
A. P.
(
2006
).
Material properties and biochemical composition of mineralized vertebral cartilage in seven elasmobranch species (Chondrichthyes)
.
J. Exp. Biol.
209
,
2920
-
2928
.
Porter
,
M. E.
,
Koob
,
T. J.
and
Summers
,
A. P.
(
2007
).
The contribution of mineral to the material properties of vertebral cartilage from the smooth-hound shark Mustelus californicus
.
J. Exp. Biol.
210
,
3319
-
3327
.
Porter
,
M. E.
,
Roque
,
C. M.
and
Long
,
J. H.
, Jr
(
2009
).
Turning maneuvers in sharks: predicting body curvature from vertebral morphology
.
J. Morphol.
270
,
954
-
965
.
Porter
,
M. E.
,
Diaz
,
C.
, Jr
,
Sturm
,
J. J.
,
Grotmol
,
S.
,
Summers
,
A. P.
and
Long
,
J. H.
, Jr
(
2014
).
Built for speed: strain in the cartilaginous vertebral columns of sharks
.
Zoology
117
,
19
-
27
.
Rief
,
M.
,
Gautel
,
M.
,
Oesterhelt
,
F.
,
Fernandez
,
J. M.
and
Gaub
,
H. E.
(
1997
).
Reversible unfolding of individual titin immunoglobulin domains by AFM
.
Science
276
,
1109
-
1112
.
Shadwick
,
R. E.
(
1999
).
Mechanical design in arteries
.
J. Exp. Biol.
202
,
3305
-
3313
.
Shah
,
J. V.
and
Janmey
,
P. A.
(
1997
).
Strain hardening of fibrin gels and plasma clots
.
Rheol. Acta
36
,
262
-
268
.
Singh
,
T.
and
Nair
,
S. S.
(
1992
).
A mathematical review and comparison of continuously variable transmissions. SAE Technical Paper 922107
.
Storm
,
C.
,
Pastore
,
J.
,
MacKintosh
,
F.
,
Lubensky
,
T. C.
and
Janmey
,
P. A.
(
2005
).
Nonlinear elasticity in biological gels
.
Nature
435
,
191
-
194
.
Summers
,
A. P.
and
Long
,
J. H.
, Jr.
(
2006
).
Skin and bones, sinew and gristle: the mechanical behavior of fish skeletal tissues
. In
Fish Biomechanics
, Vol.
23
(ed.
R. E.
Shadwick
and
G. V.
Lauder
), pp.
141
-
177
.
New York
:
Academic Press
.
Tytell
,
E. D.
and
Lauder
,
G. V.
(
2004
).
The hydrodynamics of eel swimming
.
J. Exp. Biol.
207
,
1825
-
1841
.
Tytell
,
E. D.
,
Hsu
,
C.-Y.
and
Fauci
,
L. J.
(
2014
).
The role of mechanical resonance in the neural control of swimming in fishes
.
Zoology
117
,
48
-
56
.
Yao
,
N. Y.
,
Broedersz
,
C. P.
,
Depken
,
M.
,
Becker
,
D. J.
,
Pollak
,
M. R.
,
MacKintosh
,
F. C.
and
Weitz
,
D. A.
(
2013
).
Stress-enhanced gelation: a dynamic nonlinearity of elasticity
.
Phys. Rev. Lett.
110
,
018103
.
Yu
,
W.
,
Wang
,
P.
and
Zhou
,
C.
(
2009
).
General stress decomposition in nonlinear oscillatory shear flow
.
J. Rheol.
53
,
215
.
Zar
,
J. H.
(
1999
).
Biostatistical Analysis
.
New Jersey
:
Prentice Hall
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information