## ABSTRACT

Both hummingbirds and insects flap their wings to hover. Some insects, like fruit flies, improve efficiency by lifting their body weight equally over the upstroke and downstroke, while utilizing elastic recoil during stroke reversal. It is unclear whether hummingbirds converged on a similar elastic storage solution, because of asymmetries in their lift generation and specialized flight muscle apparatus. The muscles are activated a quarter of a stroke earlier than in larger birds, and contract superfast, which cannot be explained by previous stroke-averaged analyses. We measured the aerodynamic force and kinematics of Anna's hummingbirds to resolve wing torque and power within the wingbeat. Comparing these wingbeat-resolved aerodynamic weight support measurements with those of fruit flies, hawk moths and a generalist bird, the parrotlet, we found that hummingbirds have about the same low induced power losses as the two insects, lower than that of the generalist bird in slow hovering flight. Previous analyses emphasized how bird flight muscles have to overcome wing drag midstroke. We found that high wing inertia revises this for hummingbirds – the pectoralis has to coordinate upstroke to downstroke reversal while the supracoracoideus coordinates downstroke to upstroke reversal. Our mechanistic analysis aligns with all previous muscle recordings and shows how early activation helps furnish elastic recoil through stroke reversal to stay within the physiological limits of muscles. Our findings thus support Weis-Fogh's hypothesis that flies and hummingbirds have converged on a mechanically efficient wingbeat to meet the high energetic demands of hovering flight. These insights can help improve the efficiency of flapping robots.

## INTRODUCTION

Hovering is one of the most energy-consuming flight modes in rotorcraft and animal flight (Lasiewski, 1963; Weis-Fogh, 1972; Suarez et al., 1991; Chai and Dudley, 1995; Leishman, 2006). Hummingbirds (*Trochilidae*) are the only vertebrates that can continuously hover in still air. Unlike other birds, hummingbirds invert their wings on the upstroke like insects – a remarkable example of behavioral convergence despite profound differences in body plans (Warrick et al., 2012). In 1972, Weis-Fogh analyzed the aeromechanics of hovering *Drosophila* and hummingbirds (Weis-Fogh, 1972). His first principle model predicts (i) that hummingbirds and flies must have converged on harnessing elastic recoil to beat their wings back and forth in resonance and (ii) that they need to lift their body weight symmetrically over the downstroke and upstroke to save energy. These predictions have been confirmed for *Drosophila* (Pringle, 1957; Alexander and Bennet-Clark, 1977; Dickinson and Lighton, 1995; Dickinson et al., 1999; Muijres et al., 2014); however, testing them for hummingbirds has resulted in a multi-decade quandary based on stroke-averaged muscle power models (Weis-Fogh, 1972; Ellington, 1984). Many studies report two different stroke-averaged power requirements for hummingbirds: one under the assumption of perfect elastic storage and one assuming zero elastic storage (Wells, 1993a,b; Chai and Dudley, 1995, 1996; Chai and Millard, 1997; Welch et al., 2007; Ortega-Jimenez and Dudley, 2012; Kim et al., 2014).

Previous *in vivo* airflow measurements and computational fluid dynamics (CFD) simulations have shown that asymmetry in the wing chamber (Warrick et al., 2005) restricts the ability of hummingbirds to support their body weight equally during the downstroke and upstroke (Warrick et al., 2005, 2009; Wolf et al., 2013; Song et al., 2014, 2015). The aerodynamic force was, however, calculated based on shed vorticity or CFD models; it has never been measured directly throughout the wingbeat. In contrast, hummingbird muscle function has been established within a wingbeat through *in vivo* electromyography (Hagiwara et al., 1968; Altshuler et al., 2010b; Tobalske et al., 2010; Donovan et al., 2012; Mahalingam and Welch, 2013), sonomicrometry (Tobalske et al., 2010) and isometric force measurements (Hagiwara et al., 1968), but the findings are unique relative to those for other larger birds. Hummingbirds activate their flight muscles at midstroke, roughly a quarter period in advance of stroke reversal (Altshuler et al., 2010b; Tobalske et al., 2010; Donovan et al., 2012; Mahalingam and Welch, 2013), whereas larger birds activate their flight muscles around stroke reversal (Biewener et al., 1992; Tobalske and Biewener, 2008), which cannot be fully reconciled via scaling (Tobalske, 2016). Further, once activated, the hummingbird's pectoralis muscle reaches peak isometric muscle force (Hagiwara et al., 1968) in under 8 ms – much faster than the pectoralis of any other bird (e.g. ∼32 ms in starlings: Biewener et al., 1992; ∼39 ms in pigeons: Biewener et al., 1998; and ∼71 ms in mallards: Williamson et al., 2001).

To reconcile these unusual attributes of muscle contractile dynamics, and to determine how hummingbirds could rely on aerodynamic symmetry and elastic recoil to hover efficiently as hypothesized for *Drosophila* (Weis-Fogh, 1972; Dickinson and Lighton, 1995), we analyzed the instantaneous aerodynamic force they generate within a wingbeat. We used a new sensitive aerodynamic force platform (AFP; Fig. 1A) to mechanically integrate the Navier–Stokes control-surface integral for the aerodynamic force (Lentink et al., 2015; Hightower et al., 2017; Lentink, 2018) generated by a hummingbird. We then compared our vertical aerodynamic force measurements for hummingbirds with data for fruit flies (following Weis-Fogh), hawk moths (an insect with similar flight behavior; Willmott and Ellington, 1997; Wakeling and Ellington, 1977) and parrotlets (a generalist bird not capable of sustained hovering; Chin and Lentink, 2017) to better understand the specific aerodynamic adaptations hummingbirds rely on for energy-efficient hovering. Finally, combining the novel, direct aerodynamic force measurements with detailed wing kinematics analyses and aeromechanical models, we determined the torque and power the hummingbird flight muscles must generate to coordinate the wingbeat during hovering.

## MATERIALS AND METHODS

### Hummingbirds

Anna's hummingbirds, *Calypte anna* (Lesson 1829), were studied at Stanford University, CA, USA, in March 2017. The procedures were approved by the Stanford Administrative Panel on Laboratory Animal Care and carried out under appropriate Federal and California scientific collecting permits. Over a period of 4 days, we captured and trained six male Anna's hummingbirds to hover in front of a stationary feeder (Fig. 1B). Each bird was then placed in a custom-developed AFP with a flight volume of 50×50×50 cm^{3}, a perch and an artificial flower, with *ad libitum* access to a 25% sugar solution (Fig. 1A). The AFP was used to measure body weight support in flight and an instrumented perch was used to measure resting, take-off and landing ground reaction forces (Fig. 1C). In addition, we used five synchronized high-speed cameras to measure wing kinematics (Fig. 2A). The experiment was performed at a field station to minimize acoustic noise from external sources. Five flights per bird were recorded (∼4 s each), which included feedings ranging from 72 to 173 wingbeats. The aerodynamic force measurements from one flight of birds 3 and 4 were not analyzed, because of a lack of synchrony between the high-speed cameras and AFP. The kinematics from all five flights of all six birds were analyzed. After five recorded flights, each bird was marked with a temporary non-toxic dot of paint and released at the location of capture on the same day. We calculated bird mass based on the average mass recorded by the instrumented perch before and after each flight (4.9±0.4 g, mean±s.d.; *N*=6 birds), and found that all birds gained 4–13% of their initial mass during the experiment.

### Time-resolved *in vivo* force measurements

We measured the aerodynamic force *in vivo* by developing a new AFP (Fig. 2A) that is sensitive enough to mechanically integrate the high-frequency low-amplitude pressure field generated by a hummingbird within a wingbeat (Fig. 2B). The working principle and precision of AFPs have been demonstrated using analytical solutions of the Navier–Stokes control surface formulation (Lentink, 2018), design guidelines for high-frequency measurements (Hightower et al., 2017), a quad copter to demonstrate that independent ground-truth unsteady force measurements indeed line up within sensor resolution (Lentink et al., 2015), and *in vivo* recordings of parrotlets to demonstrate the effectiveness of the platform for working with freely behaving animals (Chin and Lentink, 2017). The platform we present here follows the same working principles, but measures aerodynamic pressure force at a faster sampling rate (6 versus 1 kHz) and the resonance frequency of the structure is at least 2.5 times faster (260 versus 102 Hz). The aerodynamic forces (Fig. 2C) generated by the hummingbird hovering at the center of the volume are transmitted to the walls of the platform by pressure waves that travel at the speed of sound (Hightower et al., 2017), which takes about 3% of the wingbeat period (∼0.7 ms). The vertical aerodynamic force is measured with a bottom and a top plate, which are custom-built carbon fiber sandwich panels (50×50 cm^{2}; 498 and 496 g; KVE Composites Group, Den Haag, The Netherlands). The plates integrated the vertical pressure force over the ceiling (Fig. 2D) and floor (Fig. 2E) of the flight volume, which is enclosed by vertical acrylic side walls. The acrylic sidewalls are not instrumented, and therefore we ignored the small shear forces due to boundary layer flow development along the walls, which constitute only a few percent of the total vertical force (Fig. 2F) in this setup, at most (Lentink et al., 2015). Thin strips of Saran wrap cover the small gaps between the plates and side walls to prevent air leakage while minimizing mechanical coupling (Hightower et al., 2017). Each plate is supported with three Vee Blocks (VB-375-SM, Bal-tec, Los Angeles, CA, USA) and corresponding spherical contact surfaces fixed to the plate, which results in six statically defined contact points per plate (kinematically constrained). Tensioned elastic bands secure the connection of the inverted top plate to the Vee Blocks by counteracting gravity. Each Vee Block is mounted on a Nano 43 six-axis force/torque sensor (6 kHz sampling rate, silicon strain gage based, with SI-9-0.125 calibration, 2 mN resolution, ATI Industrial Automation, Pittsburgh, PA, USA). Two additional Nano 43 sensors measured forces on the artificial flower as well as the perch (Fig. 2B,C) to determine the bird's weight (averaged over 0.2 s before and after each flight while the bird was resting). The average rms noise level in the force measured before takeoff was 4.5 mN (9.4% of body weight) for the AFP and 4.2 mN (8.7% of body weight) for the instrumented perch. The noise is due to ambient background noise and is close to sensor resolution (2 mN). The average absolute difference between the unloaded AFP before and after the flight was 1.2 mN (2.5% of body weight) while the average difference between the perch forces was 0.8 mN (1.7% of body weight). These are a factor of two lower than sensor resolution (2 mN) and are therefore inconsequential (differences may be caused by hummingbird urination and sensor drift). Vertical forces were filtered offline (8th order digital low-pass Butterworth filter with cut-off frequency of 180 Hz; ∼4.4 times the wingbeat frequency; unfiltered traces shown as light colors in Fig. 2C–F). The 180 Hz cut-off frequency was chosen to filter out the natural frequencies of the AFP, which were obtained by measuring the impulse response of the system. The structural natural frequencies were measured by tapping the top and bottom plate with a stiff carbon fiber rod (Fig. 3A–D). A comparison of the frequency spectrum of the force trace generated by a hummingbird and the structural natural frequency spectrum of the AFP is shown in Fig. 3E. The large peak near twice the wingbeat frequency of ∼79 Hz demonstrates that the downstroke and upstroke each generate a vertical force pulse (which doubles the number of force extremes within a beat and thus the frequency).

### High-speed filming

Five high-speed cameras, four grayscale Phantom Miro M310 and one color Phantom Miro LC310 (Vision Research, Wayne, NJ, USA), were synchronized and positioned around the flight chamber as shown in Figs 1A and 2A. After the bird hovered in front of the feeder, the cameras (2000 Hz; 490 μs exposure; 1280×800 pixel resolution) were post-triggered, resulting in flight recordings (Fig. 4) of 8310 frames (∼4 s). Camera-ready signals were also sent to the DAQ (USB-6210, National Instruments, Austin, TX, USA) to synchronize forces and kinematics. In total, five high-speed recordings of each bird were taken while it was hovering at the feeder.

### 3D camera calibration

The high-speed cameras were 3D calibrated by sweeping a wand with two black spheres at a fixed distance (2.99 mm diameter at 43.70 mm distance apart) through the volume horizontally and vertically. A single black sphere was also dropped through the volume to define the gravitational *z*-axis direction. A direct linear transformation (DLT) calibration was performed using easyWand5 (Hedrick, 2008), which resulted in wand quality scores of 0.40, 0.58, 0.59, 0.37, 0.63, 0.57 and 0.53 (a score of less than 1 indicates a good calibration) (Hedrick, 2008; Theriault et al., 2014).

### Automated wingbeat segmentation for calculating average force within a beat

We determined the average vertical aerodynamic force trace during a wingbeat by analyzing stationary hovering flight at the feeder (Movie 1). Steady hovering flight was defined as the bird hovering at the feeder with its bill in the artificial flower, which we automatically detected with a custom-written MATLAB (R2015b, MathWorks, Sunnyvale, CA, USA) script. We excluded the first three and last three wingbeats at the feeder to ensure stationary hovering. Wingbeat transitions were also determined automatically by tracking the horizontal position of the centroid of the birds' silhouettes (side view, Movie 2). These video data were used to segment the force measurements and automatically identify the individual force trace for each wingbeat. Next, the average vertical force trace was determined for a wingbeat, which was subsequently normalized with the average weight of the bird (as measured by the perch). Upstroke weight support was calculated by integrating the vertical aerodynamic force trace over the upstroke and dividing it by the integral of vertical aerodynamic force over the entire wingbeat (upstroke impulse/total wingbeat impulse). Further, we calculated the flapping frequency, the relative duration of the downstroke and the total weight support for each wingbeat (Table S1). When we applied the maximum and minimum stroke angle definition for wingbeat transition, it did not change the transitions more than one frame (∼2% of wingbeat period) or upstroke weight support percentage more than 1% (tracking the kinematics of five wingbeats per bird). Finally, we determined the normalized weight support trace within a complete wingbeat for all six hummingbirds; for this, we first interpolated every force sample to 1000 samples per wingbeat, then calculated the average force trace (Fig. 5A) over the wingbeats (*n*=568–789) of all birds (*N*=6), totaling 4094 wingbeats. Average values for each bird are given in Table S1.

### Wing kinematics measurement

We determined the kinematics of the right wing by selecting and analyzing five wingbeats from each bird (*N*=6 birds; *n*=5 wingbeats per bird; 30 wingbeats total). From each of the five independent flight recordings per bird, one wingbeat was randomly selected from a subset of wingbeats that matched the average wingbeat frequency. Next, we projected a 3D linearly twisted planform of the right wing of the bird in the five camera views (Fig. 4) using the DLT calibration (Hedrick, 2008) (custom MATLAB GUI). The 3D representation of the wing was defined by seven parameters and best fitted across all five views by hand. Three parameters defined the 3D position (*x*, *y*, *z*) of the wing, and three Euler angles defined the rotation about the *x*-, *y*- and *z*-axes. The seventh parameter defined linear twist about the long axis of the wing. These seven parameters were manually adjusted to best fit the twisted and projected hummingbird wing outline in all five camera views. Finally, we filtered the manually tracked wing kinematics parameters with a 4th order digital low-pass Butterworth filter with a cut-off frequency of 400 Hz (∼10 times the wingbeat frequency).

### Wing kinematics parameters

The measured wing kinematics were converted to the standard decomposition of wing stroke, deviation and angle of attack (Sane and Dickinson, 2001) (Fig. 6A,C). We improved the description of wing kinematics by including linear wing twist (Fig. 6B,C), considering the kinematics of 100 blade element segments (Dickson et al., 2008; Leishman, 2006) (see Appendix, Fig. A1) and by determining the projected vertical ‘actuator disk’ area of the wing (Fig. S1). The wing planform (Fig. 4) fitted well throughout the wingbeat cycle, because the upstroke to downstroke span ratio is approximately constant for hovering hummingbirds (Tobalske et al., 2007). The wing tip location was used to calculate the stroke and deviation angle with respect to the horizontal stroke plane. We determined the orthogonal stroke and deviation angular velocities using the ‘perfect smoother’ developed by Eilers (2003). This smoother optimized the time derivatives we needed to determine the wing angular velocity vector (see Appendix, Fig. A1A). Eilers' smoother (Eilers, 2003) was also used to determine the wing angular acceleration vector (see Appendix, Fig. A1B). Subsequently, we divided the wing planform into 100 wing segments, equally spaced down the span of the wing along the midline. The wing's angular velocity and acceleration vectors were crossed with the position vector of each wing segment to obtain the local velocity and acceleration over the wingbeat (see Appendix, section A2). The local aerodynamic angle of attack of the wing was defined as the angle between the chord vector and the wing velocity vector. The chord vector at each spanwise wing element pointed from the midspan line to the leading edge of the wing. As the wing is twisted, the angle of attack varied along the span (see Appendix, Fig. A1C). Because the wing is inverted on the upstroke, the angle of attack with respect to the wing was negative (Kruyt et al., 2014). This inversion allowed us to use an angle of attack-specific drag coefficient for both the downstroke and upstroke (see Fig. S2 and Appendix, section A3). The angle of attack was ill defined at stroke reversal, because wing velocity approached zero. We therefore calculated the instantaneous local rotational velocity (see Appendix, Fig. A1D) using the rotational velocity of the root and tip of the wing (see Appendix, Eqn A4e), which accounted for angular velocity caused by wing twist and eliminated error caused by the velocity singularity. Finally, we calculated the area swept by the wing by creating a point cloud based on the position of the wing outline over the whole stroke. These points were projected into the horizontal stroke plane, which we used to calculate the vertical ‘actuator disk’ area of the right wing shown in Fig. S1 (using an ‘alpha hull’ with a probe radius of 10 mm in a custom-written MATLAB script). To also account for the left wing, we multiplied the actuator disk area by a factor of two.

### Force and power model for hovering flight

A quasi-steady wing element model (Sane and Dickinson, 2002; Dickson et al., 2008; Song et al., 2015) was used to estimate the instantaneous aerodynamic force (Fig. 6D,E) and power (Fig. 6F,G). The aerodynamic force calculation included translational lift and drag components, as well as rotational and added mass forces (Fig. 6D; see Appendix, section A4). Drag contributed to vertical weight support, because the stroke plane is U-shaped, which points the drag force out of the horizontal plane. Translational lift and profile drag calculations were based on wing spinner force measurements with prepared *C. anna* wings (Kruyt et al., 2014) (Fig. S2A). Rotational force was calculated with our measured wing rotational velocity in combination with a rotational force coefficient derived based on the axis of rotation of a hummingbird wing (Song et al., 2015) (see Appendix, Eqn A4c). Finally, the added mass force was calculated based on flat plate theory (Dickson et al., 2008; Sedov, 1980) in combination with our measured wing acceleration (see Appendix, Eqn A4f). We integrated the vertical force contributions over all 100 spanwise wing elements. Aerodynamic torque (Fig. 7) was calculated by crossing the aerodynamic force at each blade element with its radial position from the shoulder (see Appendix, section A8). A similar quasi-steady wing element model was developed to calculate the aerodynamic power needed to hover (Fig. 6F,G). This model included profile, induced, rotational and added mass power (see Appendix, section A5). Compared with earlier experimental studies, we improved the profile and induced power calculations as follows: we calculated profile power by separating measured translational drag (Kruyt et al., 2014) into its induced and profile drag components, so that we could use profile drag coefficient as a function of angle of attack to calculate profile drag and power (Fig. S2B; see Appendix, Eqn A3b), instead of assuming a constant profile drag coefficient (Ellington, 1984; Chai and Dudley, 1995, 1996; Altshuler et al., 2010a). As our unique experimental setup measures vertical force *in vivo*, we were also able to calculate the induced power (Leishman, 2006) by using the actual vertical force generated by the hummingbird (Fig. 5A). Aerodynamic power due to rotation and added mass were calculated by dot multiplying the forces with the local wing velocity vector and integrating them for all wing elements over the whole stroke, which was similar to the force calculation above. Inertial power is calculated based on wing mass and radius of gyration for *C. anna* wings (Fig. 6G; see Appendix, section A7). We determined the body mass-specific power (Fig. 6F,G) by dividing the power by the mass of each bird, and the muscle mass-specific power (Fig. 8) by dividing the power by the mass of the associated flight muscle. Each muscle was defined to be involved from the moment the power dropped below zero through its storage phase and release until the power once again dropped below zero as shown in Fig. 8B. This definition agrees with the EMG firing and muscle tension traces (Fig. 7D) as well as the sign of the torque (Fig. 7B). We adopted EMG traces from the pectoralis and supracoracoideus (Fig. 7D) measured for Anna's hummingbirds (Altshuler et al., 2010b), because these measurements were also made for hovering flight. Whereas these particular EMG signals are representative for hovering Anna's hummingbirds, flight muscle EMG traces have multiple peaks during maximal load-lifting trials in Anna's hummingbirds (Altshuler et al., 2010b) and may vary somewhat between species (Tobalske et al., 2010). Combined, the two flight muscles represent approximately 25% of the hummingbird's body mass (Chai and Dudley, 1995, 1996), with the pectoralis and supracoracoideus representing approximately 17% (Clark, 2009) and 8%, respectively. Models and calculations are presented in full detail in the Appendix.

## RESULTS

### An active upstroke allows for efficient hovering

Our instantaneous force measurements show that Anna's support 28±1% of their bodyweight during the upstroke (Fig. 5A; *N*=6 birds; *n*=568–789 wingbeats per bird). Upstroke weight support is thus appreciable but asymmetric in hummingbirds as previously estimated based on flow measurements (Warrick et al., 2005, 2009; Wolf et al., 2013) and CFDs models (Song et al., 2014, 2015). Comparing the *in vivo* hummingbird (Fig. 5A) and *in vivo* Pacific parrotlet (Chin and Lentink, 2017) (Fig. 5B) force measurements, we see that hummingbirds generate appreciably more weight support during the upstroke than generalist birds in slow or hovering flight. Whereas hummingbirds also generate more upstroke support than hawk moths (Fig. 5D), they do not generate as much upstroke support as fruit flies (Fig. 5C). On average, we measure 99±3% of the vertical aerodynamic force that Anna's hummingbirds need to support their bodyweight in flight. We derived simple aerodynamic theory (see Appendix, section A1) to show how generating more equally distributed weight support over the wingbeat reduces the induced power requirements. Because of the non-linearity of the induced power equation, minor deviations from generating constant vertical force have small effects, while larger deviations have an increasingly greater induced power penalty (in the form of a temporal cost factor; see Appendix, Eqn A1e; Fig. 5F). This explains why we find that hummingbirds incur a somewhat lower penalty (1.25±0.03) than *Drosophila* (1.29), because weight support fluctuates more in Drosophila, whereas both incur a lower cost than the Pacific parrotlet (1.35). Overall, the hummingbird induced power penalty is similarly low in hummingbirds, *Drosophila* and hawk moths, compared with that in Pacific parrotlets (Fig. 5F). This mechanistically underpins why hummingbirds converged closer towards generating more upstroke weight support, relative to their generalist bird ancestors (McGuire et al., 2014), to reduce the cost of sustained hovering flight.

### Quasi-steady model predicts high inertial requirements

The quasi-steady aerodynamic force model (Fig. 6; *N*=6 birds; *n*=5 wingbeats per bird) for Anna's hummingbirds accounts for 64±4% of bodyweight support and predicts that the upstroke contributes 39±1% to the total quasi-steady vertical force (Fig. 6E). Stroke-averaged vertical force is composed of ∼84% translational lift, ∼0% translational drag, ∼5% rotational force and ∼11% added mass force (Fig. 6D). Our detailed 3D kinematic wing model showed that hummingbird wings twist up to −38 deg during the downstroke and up to 62 deg during the upstroke (Fig. 6C). Dynamic changes in wing twist affects both the magnitude and radial distribution of rotational and added mass forces acting on the wing during pronation and supination. The model's capacity to predict most of the weight support in hummingbirds is similar to predictions for *Drosophila* (Dickinson et al., 1999; Dickson et al., 2008) but below 100% because of additional effects not captured by a quasi-steady model. Therefore, aerodynamic power must also be underestimated, because induced power is proportional to vertical aerodynamic force to the power of 1.5 (see Appendix, Eqn A1b). To overcome these underestimates (Weis-Fogh, 1972; Chai and Dudley, 1995; Kruyt et al., 2014), as well as the over-prediction of upstroke weight support, we calculated the induced power (Fig. 6F) directly based on the measured *in vivo* vertical force (Fig. 5A). The total aerodynamic power consists of ∼63% induced power, while profile power represents ∼23%. The power due to rotational force and to added mass force is relatively small, ∼7% and ∼7%, respectively (Fig. 6F). In total, Anna's hummingbirds need a stroke-averaged aerodynamic power of 35±1 W kg^{−1} of body mass, with peak values of ∼100 W kg^{−1} (Fig. 6G). Whereas these sustained power values are amongst the highest reported for skeletal muscle, even more power is needed to overcome wing inertia at stroke reversal (Wells, 1993b). With the inertia added, the hummingbird flight muscles, the pectoralis and supracoracoideus, need to deliver a peak power of almost 200 and 150 W kg^{−1} of body mass, respectively, to beat the wings back and forth. This power value cannot be explained with existing theories unless elastic storage in the flight apparatus is assumed (Weis-Fogh, 1972; Chai and Dudley, 1995, 1996; Ellington, 1984; Wells, 1993b) (Fig. 6G).

## DISCUSSION

### Wing inertia determines muscle activation timing

To better understand hummingbird muscle function, we calculated (see Appendix, section A8) the net torque around the shoulder joint, wingbeat resolved based on wing drag and inertia (Fig. 7A). The net torque switches direction at midstroke for both flight muscles, because wing inertia is high and the aerodynamic drag is insufficient to decelerate the wing (Fig. 7B). The antagonistic flight muscle pair has to generate the instantaneous net torque needed to overcome wing inertia and drag, which forms an aero-mechanical torque loop throughout the wingbeat (Fig. 7C). The effect of wing inertia dominating aerodynamic drag is that it tilts the torque loop, as predicted based on first principles by Weis-Fogh (1972). This causes wing torque to peak at stroke reversal, instead of midstroke, and dictates the effective work that the pectoral and supracoracoideus muscle–tendon units need to deliver together. This ‘internal’ net muscle-generated torque×angular displacement has to be equal to (or larger than) the ‘external’ net wing torque×angular displacement (depending on friction and other internal losses; Weis-Fogh, 1972). We can determine which muscle–tendon unit has to perform which effective work, because the pectoralis delivers the net torque (muscle force×lever arm) required for upstroke to downstroke reversal, while the supracoracoideus coordinates downstroke to upstroke reversal. The antagonistic flight muscle pair must accommodate this peak inertial torque requirement to reverse the stroke direction. The pectoralis (main downstroke muscle; Altshuler et al., 2010b; Tobalske et al., 2010; Donovan et al., 2012; Mahalingam and Welch, 2013) is the primary muscle that can pull to generate the positive braking torque required to coordinate the upstroke to downstroke reversal. Similarly, the supracoracoideus (main upstroke muscle; Altshuler et al., 2010b; Tobalske et al., 2010; Donovan et al., 2012; Mahalingam and Welch, 2013) has to be activated well before the start of the upstroke to generate the negative braking torque required for downstroke to upstroke reversal (Fig. 7B). Indeed, it has been shown (Altshuler et al., 2010b; Donovan et al., 2012) that Anna's hummingbirds activate the pectoralis and supracoracoideus a quarter of a wingbeat in advance of the upstroke–downstroke and downstroke–upstroke transition, respectively (Fig. 7D) – approximately a quarter-stroke earlier than found in other birds (Altshuler et al., 2010b; Tobalske et al., 2010; Mahalingam and Welch, 2013). The activation appears to be well timed, because Hagiwara et al. (1968) found that the pectoralis of Anna's hummingbirds reaches peak isometric force in under 8 ms after activation *in situ*, regardless of whether it is activated by a single pulse (twitch) or a pulse train at 23 or 53 Hz [a range that encompasses Anna's hummingbird wingbeat frequency of 41 Hz, although Hagiwara et al. (1968) even showed noticeable force peaks with tetanus from stimulations of up to 300 Hz; Fig. S3]. Comparing the time to peak isometric force of Anna's hummingbirds (8 ms; 4.9 g body mass) with values for the Etruscan shrew extensor digitorum longus (11 ms; 1.8 g body mass), the fastest mammal locomotory muscle recorded to date (Peters et al., 1999; Jürgens, 2002), it becomes clear that scaling effects cannot fully explain the speed of Anna's hummingbird pectoralis muscle. It is reasonable to assume the pectoralis and supracoracoideus are similarly fast in Anna's hummingbirds, because the fiber type of both flight muscles is uniformly fast oxidative glycolytic (Donovan et al., 2012; Welch and Altshuler, 2009). These fibers enable Anna's hummingbird pectoralis to reach peak force faster than the so-called ‘superfast’ muscles that control the dove's trill (Elemans et al., 2004). Superfast muscles in other animals are not generally known to be primary motors for locomotion (Fuxjager et al., 2016), because they generate low specific power at high frequencies (Elemans et al., 2004, 2008; Fuxjager et al., 2016; Rome, 2006). Similar early activation in frog muscle–tendon units corresponds to a more spring-like work loop that stores elastic energy, albeit imperfectly (Sawicki et al., 2015). Assuming isometric muscle measurements are a reasonable approximation to estimate the time between activation and peak muscle force *in vivo*, we find that the resulting peak muscle force measured by Hagiwara et al. (1968) (Fig. 7D) aligns with the peak torque (muscle force×moment arm) demand at stroke reversal that we calculated (Fig. 7B). It would require a further 6 ms delay (¼ wingbeat) in muscle force development due to contractile kinetics for the muscle peak to instead coincide with maximum wing velocity, much longer than the 2 ms uncertainty we found in the definition of stroke reversal. Whereas torque peaks at stroke reversal, power demand (Fig. 6G) is approximately zero, because wing angular velocity reaches zero. This is consistent with the near-zero instantaneous strain rate in the pectoralis muscle (low power) at stroke reversal in Rufous hummingbirds (Tobalske et al., 2010).

### Recoil in the flight apparatus explains peak instantaneous power

In the absence of elastic storage, a hummingbird would have to generate 195±11 and 347±29 W kg^{−1} of muscle mass in the pectoralis and supracoracoideus, respectively (averaged over a full wingbeat, assuming negative power is free; Fig. 8A). By storing all negative power elastically and releasing it to accelerate the wing again, the hummingbird would only have to generate 105±11 and 220±16 W kg^{−1} of power in the pectoralis and supracoracoideus, respectively (Fig. 8A). Weis-Fogh and Alexander (1977) estimated the theoretical maximum is about 250 W kg^{−1} muscle mass based on hummingbirds' physiological limits, while Pennycuick and Rezende (1984) estimated it to be as high as 430 W kg^{−1} muscle mass. Josephson (1993) suspected these are overestimates due to unrealistic muscle contractile parameter assumptions and a failure to account for activation/inactivation time. Although these power outputs are already very high for vertebrates, the flight muscles are capable of delivering over three times more profile power (Altshuler et al., 2004) during maximum load-lifting trails (Chai and Millard, 1997; Altshuler et al., 2004; Skandalis et al., 2017). Using additional load-lifting data from Anna's hummingbirds (Altshuler et al., 2010b; Segre et al., 2015), we confirm this estimate and calculate that the induced and inertial power also triple during burst muscle performance (see Appendix, section A10). Further, the nominal instantaneous total power peaks in hovering hummingbirds are between 1000 and 2000 W kg^{−1} muscle mass for each muscle (Fig. 8B), exceeding peak power during burst takeoff in quail (1121 W kg^{−1}; Askew and Marsh, 2001) and matching peak power in projecting chameleon tongues, which use elastic storage to achieve it (1892 W kg^{−1}; Anderson and Deban, 2010). It is therefore more parsimonious to assume peak power is achieved by releasing stored elastic energy, especially because this extraordinary power output has to be sustained over many thousands of contraction cycles during foraging flight. Whereas physiological limits do not enable muscles to deliver such high power over so many cycles, springs can generate extreme instantaneous power with similar double frequency traces using near-zero work per cycle. Another way to visualize the opportunity for elastic storage in hummingbirds is to plot the aero-mechanical torque loop in the reference frame of each muscle. During forward flight, a generalist starling (at 13.7 m s^{−1}; Biewener et al., 1992), magpie (at 6 m s^{−1}; Dial et al., 1997), budgerigar (at 4 m s^{−1}; Ellerby and Askew, 2007) and pigeon during slow flight (at 3.9 m s^{−1}; Tobalske and Biewener, 2008) have a relatively upright pectoralis work loop with minimal negative work during the end of the upstroke. The supracoracoideus, however, generates more negative work, which may be stored elastically in the supracoracoid tendon (Tobalske and Biewener, 2008). In contrast, the torque loop of hovering hummingbirds is tilted as a result of high wing inertia, allowing a larger area of negative work while decelerating the wing (Fig. 8C,D). They can efficiently store the high instantaneous negative work and release it again as positive work via elastic recoil (Roberts et al., 1997; Dickinson et al., 2000; Biewener, 2003). Elastic energy can be stored in tendons (Mahalingam and Welch, 2013; Alexander, 2002) and the muscle itself (Alexander and Bennet-Clark, 1977; Dickinson et al., 2000; Sawicki et al., 2015). Finally, because wing inertia dominates aerodynamics and our predictions for inertial effects are grounded in established dynamics theory (see Appendix, Fig. A2 and section A9), our understanding of how hummingbirds use elastic recoil to tune their wingbeat are robust to aerodynamic modeling errors.

### Conclusion

Our analysis shows that hummingbirds and *Drosophila* are about equally aerodynamically efficient based on the induced power required to sustain hovering. Further, three lines of wingbeat-resolved evidence show that the most parsimonious explanation for the extreme instantaneous performance of the hummingbird flight apparatus is elastic recoil. Only recoil can reconcile (i) the instantaneous peak torque at stroke reversal (Fig. 7B), which requires exceptional muscle work, with (ii) the early activation (Altshuler et al., 2010b; Donovan et al., 2012) and fast contractile dynamics (Hagiwara et al., 1968) of Anna's hummingbird flight muscle (Fig. 7D), and (iii) the peak instantaneous power of well over 1000 W kg^{−1} for each flight muscle–tendon unit (Fig. 8B), in concert. Overall these findings support Weis-Fogh's hypothesis that *Drosophila* and hummingbirds must have converged on similarly low induced power and inertial cost to sustain hovering. The finding that recoil is critical to overcome the inertia of flapping wings at low energetic cost in hummingbirds shows how current bird-sized flapping aerial robots (Lentink et al., 2009; Keennon et al., 2012) could be made dramatically more energy efficient by adding a spring in their flap. Future studies can build upon these results by exploring the diversity in the flight apparatus of birds (Tobalske, 2016) and bats (Riskin et al., 2012; Konow et al., 2017), both in laboratory settings and in more realistic, ecological settings with animals entering the platform voluntarily to feed (see ruggedized platform in Movie 3).

### APPENDIX

Our hummingbird aerodynamics models and equations developed here integrate and expand the existing quasi-steady aerodynamics model for insect flight (Dickson et al., 2008) and the blade element model for rotorcraft (Leishman, 2006).

#### A1: induced power penalty

*P*

_{ind}, delivered by the flight muscles to the wings and is proportional to vertical force raised to the power of 1.5. This cost function intuitively shows that it is more efficient to generate constant thrust throughout the stroke because any decrease in vertical force during one part of the wingbeat will need to be compensated for by an equivalent increase in force during another part of the wingbeat, which will cost more power than is saved. We expanded the corresponding ‘actuator disk’ theory (Leishman, 2006; Chai and Dudley, 1995; Weis-Fogh, 1972; Norberg, 2012) to determine the analytical fitness function that relates the actual average induced power, , to the theoretical minimum induced power, , as a function of instantaneous weight support over the wingbeat. The average induced power, , was calculated (Leishman, 2006; Ellington, 1984) by multiplying the robot or animal's weight,

*W*, by the velocity of the induced flow,

*v*

_{ind}: where ρ is the density of air and

*A*is the area swept by the wings. This allowed us to describe as: where κ is the induced power factor that accounts for tip losses, non-uniform inflow and other non-ideal effects (Leishman, 2006). A value of κ=1.2 is typically assumed for bird flapping flight (Norberg, 2012; Pennycuick, 2008; Tobalske et al., 2003). Ellington (1984) broke this induced power factor up into a spatial correction factor, σ, and a temporal correction factor, τ: Given we measured the temporal distribution of the vertical force directly, we can quantify the additional penalty that an animal incurs by generating a time-varying wake. To decouple, we define κ

_{σ}to be the spatial cost factor (with κ

_{σ}=1 corresponding to a uniform wake) and κ

_{τ}to be the temporal cost factor (with κ

_{τ}=1 corresponding to a constant wake). Using Ellington's estimates (Ellington, 1984), we assume a value of σ=0.1, which leads to κ

_{σ}=1.1. We can then describe the ideal average induced power as: and the induced power including losses due to temporal force fluctuation as: where

*T*is the wingbeat period and

*F*(

*t*) is the aerodynamic vertical force, which must be equal to the animal’s weight on average to enable sustained hovering. The induced power cost due to time-varying vertical force can be quantified by calculating the ratio with the ideal induced power: This cost factor is equal to 1 if

*F*(

*t*)=

*W*, but will be greater than 1 if the vertical force varies during the wingbeat. As we measured the term

*F*(

*t*)/

*W in vivo*, we can directly calculate the induced power factor that the animal incurs. The calculated induced power temporal cost factor for Anna's hummingbirds, parrotlets,

*Drosophila*and hawk moths is shown in Fig. 5F.

#### A2: wing element velocity and acceleration

*r*

_{e}is the distance to each wing element. This allowed the velocity, , at each element to be calculated by: and the acceleration, , to be calculated by:

#### A3: profile drag coefficient

*C*

_{L}, and drag,

*C*

_{D}, coefficients of spinning

*C. anna*wings from Kruyt et al. (2014) (Fig. S2A) in which

*C*

_{D}included both profile and induced drag as a function of angle of attack. Based on these coefficients, we separated the total drag power into its induced (see Appendix, Eqn A1b) and profile components as follows: where the vertical force is ,

*R*is the wing radius, and are the second and third moments of area (Kruyt et al., 2014), and

*S*is the surface area of the wing. Replacing

*F*and rearranging this equation to get an expression for

*C*

_{D,prof}, we obtained:

We next substituted the specific wing geometry parameters from Kruyt et al. (2014) (, , *S*=680 mm^{2}, *A*=7260 mm^{2}), and assumed (Norberg, 2012; Pennycuick, 2008; Tobalske et al., 2003) κ=1.2, to calculate the profile drag coefficient as a function of angle of attack (Fig. S2B).

#### A4: quasi-steady vertical force

*c*

_{e}is the chord length of the wing element, Δ

*r*is the width of the element and α

_{e}is the angle of attack of the element. Similarly, the translational profile drag force of each wing element was calculated by: The induced drag is assumed to act perpendicular to the induced flow, so does not contribute to the vertical force. The rotational force at each wing element was calculated by (Dickson et al., 2008): where

*C*

_{r}is the rotational force coefficient and ω

_{e}is the rotational velocity of each element. The theoretical value for the rotational force coefficient was calculated by (Dickson et al., 2008; Fung, 2008; Sane and Dickinson, 2002): where is the non-dimensional axis of rotation. Song et al. (2015) calculated an average value of for five representative chord locations for hovering ruby-throated hummingbird kinematics. We rounded this value to (which represents the estimated accuracy of this number for

*C. anna*) to calculate the rotational force coefficient for Anna's hummingbirds in our study,

*C*

_{r}=0.8. To determine the rotational force, we calculated the rotational velocity, ω

_{e}, by first finding the rotational velocity of the wing base, ω

_{b}, and then the rotational velocity of the wing tip with respect to the base due to wing twist, ω

_{t}. The rotational velocity at each wing element (Fig. A1D) was then: where

*R*is the radius of the wing. The added mass force at each wing element was thus calculated as: This term accounted for the additional air the wing had to accelerate as it flapped. The rotational and added mass forces acted normal to the back surface of the wing. As the wing was twisted, this normal vector changed with the radial position of the wing element. We accounted for wing twist in our vertical force calculation, which sums up the vertical component of the translational, rotational and added mass forces in a custom-written MATLAB script.

#### A5: quasi-steady power model

*P*

_{P,e}, was calculated by: where

*C*

_{D,prof}is the profile component of the drag coefficient calculated above (see Appendix, Eqn A3b) and plotted in Fig. S2B. The instantaneous induced power was calculated more precisely by substituting the measured

*in vivo*vertical force in established actuator disk theory (Leishman, 2006; Ellington, 1984; Chai and Dudley, 1995, 1996; Altshuler et al., 2010a; Norberg, 2012). The induced power for both wings combined,

*P*

_{I}(

*t*), was thus calculated as: The rotational and added mass power required at each wing element were calculated by dot multiplying the force vectors by the incoming air velocity :

*P*

_{aero}(

*t*), we summed the power contribution from each wing element, multiplied by 2 to account for the two wings, and added the induced power based on the

*in vivo*measured vertical force: The calculation of induced power based on

*in vivo*force measurement allows us to quantitatively estimate the total instantaneous aerodynamic power.

#### A6: simplified wingbeat-averaged aerodynamic power estimate

*v*(

*t*) is the wing velocity, is the third moment of area and is the stroke angular velocity. By assuming a harmonic stroke, we can write: where Φ (in radians) is the stroke amplitude (half of the total swept stroke angle) and

*f*is the wingbeat frequency (1/

*T*). The integral then becomes: and the average profile power over the wingbeat for one wing is: We found that a value for of 0.16 best matched our stroke-averaged high-fidelity profile power calculation. This equation is similar to the one derived by Ellington (1984) with the exception that the stroke is assumed to be harmonic, which approximates hummingbirds well (Fig. 6C). The average induced power, , for the bird can be calculated based on the weight,

*W*, and the actuator disk area swept by the wing,

*A*=2 Φ

*R*

^{2}. If vertical force is assumed to be equal to weight throughout the stroke, the power will be underestimated. To correct for this, we use the spatial and temporal cost factors introduced above (see Appendix, section A1): We use κ

_{σ}=1.1 and κ

_{τ}=1.25 to best match our stroke-averaged high-fidelity induced power calculation. As profile and induced power alone make up 86% of the modeled total aerodynamic power, , we approximated the aerodynamic power without detailed kinematics as:

#### A7: inertial power

*P*

_{inertial}, is the dot product of the wing torque and angular velocity, (Fig. A1A). Torque is the moment of inertia about the wing base,

*I*, times the angular acceleration, (Fig. A1B). The inertial power of one hummingbird wing is thus: The moment of inertia of the wing can be calculated by dividing the wing into small chord strips and measuring their masses,

*m*, and distances,

_{i}*r*, from the wing base: The non-dimensional radius of gyration, , is defined such that: where

_{i}*m*

_{w}is the mass of one wing and

*R*is the wing length. Using the following equation, we calculated an value (0.265) using the radial distance and mass of 10 wing chord strips from a male Anna's hummingbird (Ortega-Jimenez and Dudley, 2012) (

*N*=1): As a second measure, we calculated the non-dimensional value for another species of hummingbird (

*Amazilia fimbriata fluvaitilis*) based on Weis-Fogh's inertial wing measurements. By combining Weis-Fogh's (1972) published wing moment of inertia, wing mass and wing length, we calculated the to be 0.276, which is only 4% larger than our Anna's hummingbird value: In addition, we compared the measurements from Wells (1993b) for broad-tailed (

*Selasphorus platycercus*: 0.263±0.017,

*N*=4) and rufous (

*Selasphorus rufus*: 0.248±0.007,

*N*=4) hummingbirds and found that our Anna's hummingbird calculations were within the same range (1% and 6% smaller, respectively). We then used the left and right wing mass (0.1135 and 0.0916 g) and the body mass (4.16 g;

*N*=1 provided by D. L. Altshuler at The University of British Columbia, Canada) to calculate the average (single) wing mass as a fraction of body mass, which is 2.47%. This percentage was used in combination with our measured bird mass and wing length to estimate the wing mass and calculate each moment of inertia.

The wing rotation origin was calculated by fitting a line through the midspan of the wing at every instant of the wingbeat. The best-fit intersection of all of these lines was defined to be the wing rotation origin that the wing radius was measured from, as well as the rotation angles. This is a more accurate location than the physical shoulder, as it is mathematically the best-fit origin of rotation for the wing planform. Over each wingbeat, we assumed a constant wing length and thus moment of inertia, which is fair for hovering hummingbirds (Tobalske et al., 2007).

#### A8: quasi-steady torque model

*z*component of these torques as we summed the contributions from each element: where is a unit vector of [0 0 1]. Lastly, we calculated the torque due to the induced drag,

*Q*

_{I}, which acted parallel to the induced flow (vertical

*z*direction). As the induced power is torque times the

*z*component of angular velocity, we calculated the induced torque as follows: where the factor of 1/2 accounted for the torque of a single wing. This equation introduced singularities at stroke transitions where angular velocity is zero. Because of measurement error, the induced power and angular velocity do not simultaneously approach zero at wingbeat transitions. Small phase delay errors between power and angular velocity cause the singularity, which we negated numerically through a linear fit between 40% and 60% of the wingbeat for the downstroke to upstroke transition and between 90% and 10% of the wingbeat for the upstroke to downstroke transition. Because induced torque is insignificant relative to inertial torque (Fig. 7B), this simplification does not change our conclusions. The aerodynamic torque was calculated by summing the lift, drag, rotational, added mass and induced components: Finally, for our model, we calculated the inertial torque in the vertical direction by multiplying the

*z*component of the angular acceleration by the moment of inertia: To further simplify, by assuming a harmonic stroke in the horizontal plane, we can insert the stroke angular acceleration, , explained above (see Appendix, Eqn A7h), to get:

#### A9: simplified wingbeat-resolved power and torque model

*W*is the bird's weight. The induced power can then be estimated by: where the absolute value is needed to ensure the induced power is positive as the induced flow is always downwards through the actuator disk. The induced torque for a single wing can then be estimated by: where the negative half of the sine function can be incorporated into the piecewise function so the absolute value can be removed. The other significant aerodynamic contribution is from profile drag for a sinusoidal stroke:

*C*

_{D,down}, and upstroke,

*C*

_{D,up}. The absolute value sign can be incorporated into the piecewise factor to ensure the profile power is always positive (drag always opposes velocity): The profile torque is estimated by crossing the profile drag force with the distance out that the force acts on the wing, . Because the velocity of the wing changes direction during the upstroke, the negative factor can be once again placed in the piecewise function: Combining this with the inertial power and torque estimates, one can compare the inertial and aerodynamic torque and power requirements with simple parameters. Similar to the piecewise breakdown of the aerodynamic force, the inertial power and torque can be separated into two halves with differing moments of inertia,

*I*, representing a fully extended downstroke and upstroke. This can be extended for generalist birds, which fold their wing on the upstroke, using an inertia reduction factor of

*n*:

Finally, we can compare these piecewise harmonic models with the detailed blade element models, which shows the aerodynamic and inertial variables are remarkably well approximated (Fig. A2). We thus conclude that the hummingbird wingbeat can be modeled analytically within reasonable approximation. We used the following values averaged over six Anna's hummingbirds: bodyweight *W*=0.048 N, body mass *m*=4.89 g, upstroke force fraction χ=0.28, wingbeat frequency *f*=41 Hz, spatial cost factor κ_{σ}=1.1, density of air ρ=1.2 kg m^{−3}, swept area *A*=0.0072 m^{2}, stroke amplitude Φ=1.29 radians, wing surface area *S*=0.000565 m^{2}, third moment of area , wing radius *R*=0.0514 m, downstroke drag coefficient *C*_{D,down}=0.162, upstroke drag coefficient *C*_{D,up}=0.162, wing inertia *I*=2.23×10^{−8} kg m^{2}, upstroke reduction factor *n*=1.

#### A10: maximum power estimates under load-lifting trials

^{3}and

*f*

^{3}, so the profile power margin,

*M*

_{profile}, is: Induced power scales with

*W*

^{1.5}and Φ

^{−0.5}, so the induced power margin,

*M*

_{induced}, is: Inertial power scales with Φ

^{2}and

*f*

^{3}, so the inertial power margin,

*M*

_{inertial}, is: Using these equations and the Anna's hummingbird values from the literature, we estimate that the profile, induced and inertial power increase by 3.2, 3.2 and 2.7 times, respectively, during maximum load lifting (which we round to a factor of 3 increase in muscle power for burst performance).

## Acknowledgements

We thank, B. Hightower, D. D. Chin, B. Gomez, N. R. Chiariello and T. Fukami for assistance during the field work; R. Dudley and C. D. Mendenhall for hummingbird handling instructions; F. T. Muijres for Robofly data; D. L. Altshuler for muscle and wing mass data and critically reading the manuscript; and D. D. Chin, A. H. Heers, D. B. Quinn and E. R. Tucci for detailed manuscript feedback.

## Footnotes

**Author contributions**

Conceptualization: R.I., D.L.; Methodology: R.I., D.L.; Software: R.I.; Validation: R.I., D.L.; Formal analysis: R.I., D.L.; Investigation: R.I., D.L.; Resources: R.I., D.L.; Data curation: R.I.; Writing - original draft: R.I., D.L.; Writing - review & editing: R.I., D.L.; Visualization: R.I., D.L.; Supervision: D.L.; Project administration: D.L.; Funding acquisition: D.L.

**Funding**

This research is supported by National Science Foundation CAREER Award 1552419 and the King Abdulaziz City for Science and Technology Center of Excellence for Aeronautics and Astronautics at Stanford.

**Data availability**

Data and Matlab code are available from the figshare digital repository: https://doi.org/10.6084/m9.figshare.6083393

## References

**Competing interests**

The authors declare no competing or financial interests.