While many studies involving surgical procedures have been conducted to study the pectoralis (Biewener et al., 1992; Dial and Biewener, 1993; Biewener et al., 1998; Hedrick et al., 2003; Tobalske et al., 2003; Ellerby and Askew, 2007; Tobalske and Biewener, 2008; Jackson and Dial, 2011a; Robertson and Biewener, 2012), few quantified the change in flight performance due to the surgery (Tobalske et al., 2005). While we primarily analyzed the second wingbeat that the doves made after takeoff after undergoing surgery, we recognize that this may constrain our understanding of dove flight in general. Hence, we first analyzed the effects of the surgical procedure and muscle recording cable on the overall flight dynamics for one dove for which we made recordings before as well as after the surgical procedures.

For this one dove, surgery consistently changed the way it took off and flew (Figure 1—figure supplement 2). After the surgery, the dove pushed off from the takeoff perch with similar vertical force (vertical impulse before: 0.39±0.03; after: 0.34±0.03), but substantially less forward force (horizontal impulse before: 0.81±0.03; after: 0.56±0.05; gravity was subtracted from vertical impulses; all impulses were scaled by bodyweight and integration time). During flight, the dove produced a smaller braking force after the surgery (horizontal impulse before: –0.09±0.02; after: –0.03±0.01). Additionally, the dove used more wingbeats to fly the same distance, resulting in a higher flapping frequency (before: 8.28±0.19 Hz; after: 9.20±0.15 Hz), and a different first stroke. The first downstroke was smaller in magnitude (maximum force scaled by bodyweight before: 2.18±0.13; after: 1.09±0.05) and began 14.4±4.6 milliseconds before release from the perch rather than 54.9±12.7 milliseconds after release from the perch.

The effects of the suspended recording cable were less pronounced than the effects of surgery, although still noticeable. The mass of the entire cable was 26.1 grams, but the extra mass moved by the dove was only an average of 10.9±3.8 grams (6.4% ± 2.2% of bodyweight) based on bodyweight measurements recorded while the dove was sitting on the takeoff perch. When the dove was carrying the extra mass of the cable, it pushed off from the takeoff perch in a more forward direction (vertical impulse scaled by bodyweight before: 0.34±0.03; after scaled by bodyweight plus cable: 0.29±0.05; horizontal impulse scaled by bodyweight before: 0.56±0.05; after scaled by bodyweight plus cable: 0.64±0.03). Once airborne, the dove produced greater vertical impulse (before: –0.28±0.04; after: –0.23±0.06), and greater forward impulse (before: –0.03±0.01; after: 0.02±0.03), while flapping more frequently (before: 9.20±0.15 Hz; after: 9.62±0.22 Hz). However, all of these flight differences are small compared to the differences caused by the surgery.

When we compare the aerodynamic forces produced before and after surgery by dove three (Figure 1—figure supplement 2 and Supplementary file 1b), we find an increase in flapping frequency after surgery which mirrors the differences in flight performance reported for a small passerine (Tobalske et al., 2005). Additionally, after surgery, we consistently observed a modified takeoff pattern: the dove pushed off with less forward force, and its first downstroke was earlier and produced less aerodynamic force. Overall, the flight performance of the dove was substantially altered by the surgery, and to a lesser degree, by the weight increase due to the attached recording cable. However, these changes are analogous to changing the task the dove needs to accomplish, similar to changing flight distance, speed, or angle. Our results therefore quantify this particular, post-surgery flight condition, rather than flight in general.

The four doves used somewhat different flight strategies from one another (Figure 2B, Figure 1—figure supplement 3), yet parameters during the second downstroke were similar across 20 flights (Figures 1A and 2, Supplementary file 1c). During takeoff, dove 1 pushed off with the most force, while dove 4 pushed off with the least force in both the vertical (vertical impulse for dove 1: 0.34±0.04; dove 4: 0.23±0.06) and horizontal directions (horizontal impulse for dove 1: 0.77±0.06; dove 4: 0.44±0.04; gravity was subtracted from vertical impulses; impulses herein are dimensionless as they were scaled by bodyweight (N) and integration time (s)). Because dove 1 pushed off with greater forward force, it could start braking earlier in the flight than dove 4, which still needed to produce forward thrust during the second stroke after takeoff (horizontal impulse for dove 1: –0.09±0.01; dove 4: 0.17±0.03). This may explain the difference in stroke plane angle during the second downstroke for the two doves (Figure 2B; dove 1: 31.6°±2.5°; dove 4: 45.07°±2.8°), even though there is not substantial variation in the other kinematic angles among the doves (Figure 2C): whereas the aerodynamic force magnitude needed was similar for both doves due to weight support, the force direction was different. The vertical force during the second stroke was additionally influenced by the timing and number of strokes the doves took (number of strokes between takeoff and landing for dove 1: 4; dove 2: 5; dove 3: 4; dove 4: 5), but, overall, we found that the vertical aerodynamic force was relatively consistent during the second stroke (vertical impulse for all doves: 0.86±0.09; gravity not subtracted here for clarity).

During the second stroke, the pectoralis muscle performance was consistent across the 20 flights (Figure 1A, Figure 1—figure supplement 3). The electrical activation of the pectoralis began during the middle portion of the upstroke, 7.07±5.10 milliseconds before the pectoralis began shortening, which occurred 10.54±2.86 milliseconds before the beginning of the downstroke (defined by kinematic stroke reversal: see stroke angle in Figure 2C). The electrical activation then lasted 37.63±7.66 milliseconds and ended 68.57±3.38 milliseconds before the pectoralis completed contracting, which occurred 11.06±1.99 milliseconds before the end of the downstroke. The pectoralis shortened for 59.62±3.84 milliseconds during the stroke which lasted 102.45±6.59 milliseconds (downstroke: 61.80±4.73 milliseconds). Furthermore, we see that the timing of the pectoralis contraction corresponds with the majority of aerodynamic force production.

Using the 3D angular momentum balance model (Figure 6), we can compute the time-resolved moment vector (torques) generated by the pectoralis muscle during the downstroke. Due to the complexity of 3D moments, we consider each of the three components of muscle torque on the wing separately.

First, we consider the moment in the roll axis (x_b points in the cranial direction and is aligned with the thoracic vertebrae of the dove; Figure 6B, Figure 6—figure supplement 1B). During the majority of the downstroke, the flight muscles pulled ventrally on the wing (average ventral moment during downstroke: 0.140±0.013 N-m) to oppose lift and drag, whose world z-axis components each point vertically up (Figure 3—figure supplement 1C). Just before the end of the downstroke, and 1.97±1.91 milliseconds before the pectoralis began lengthening, inertial effects from decelerating the wing were associated with a small dorsal moment (average dorsal moment at the end of the downstroke: 0.057±0.026 N-m). As is the case for power, inertia dominated during the upstroke.

Second, we consider the moment in the pitch axis (y_b points laterally towards the distal wing, perpendicular to the thoracic vertebrae; Figure 6—figure supplement 1C), in which the effects of aerodynamics dominated over inertia. When we examine the lift and drag forces (Figure 3, Figure 3—figure supplement 1), we again see that the vertical direction was the most important during the downstroke because forces in the lateral direction created no y_b moment, and the horizontal forces were small (15.5% ± 4.1% of the stroke-averaged net aerodynamic force during the downstroke). At the beginning of the downstroke, the wing was at the posterior of the bird, creating a positive aerodynamic torque in the y_b axis, which was balanced by a negative muscle moment in the y_b axis (flight muscle acting to supinate the wing). At the end of the downstroke, the wing was at the anterior of the bird, with a negative aerodynamic torque. Hence, the total muscle moment in the pitching axis acted to first supinate (average supination moment during beginning of downstroke: 0.047±0.010 N-m), then pronate (average pronation moment during end of downstroke: 0.068±0.026 N-m), the wing during the downstroke. Since the pectoralis is thought mainly to pronate the wing, it appears that other flight muscles are driving the small supination moment at the beginning of the downstroke.

Finally, we consider the moment in the yaw axis (z_b points in the dorsal direction; Figure 6—figure supplement 1D), where aerodynamics again dominate the downstroke. To gain a more holistic view of the relationship between the moment in the yaw axis and the horizontal, lateral, and vertical world components of lift and drag, it is helpful to examine Figure 6—figure supplement 1E–G, where we plot the 3D moments in the world coordinate system. The aerodynamic moment in the vertical world axis was small throughout the downstroke (average magnitude: 0.029±0.009 N-m), because the vertical lift and drag create no moment in the vertical direction (Figure 6—figure supplement 1G). Therefore, the yaw aerodynamic moment originated primarily from the moment in the horizontal world axis, rotated about the lateral world axis by the amount that the dove pitches its body forward. The net effect was to generate an average muscle moment of 0.063±0.017 N-m in the caudal direction during the downstroke.

As is the case for the power analysis, the effects of the aerodynamics dominated during the downstroke, whereas inertia dominates during the upstroke and near stroke reversal. Furthermore, for all three components of the aerodynamic moment, we find that the components of lift and drag in the vertical world direction were dominant.

To measure the electrical activation and strain of the left and right pectoralis muscles, we surgically implanted EMG electrodes and sonomicrometry crystals using standard methods for the pectoralis of birds (Biewener et al., 1998; Tobalske et al., 2005). The doves were anesthetized using isoflurane by inhalation. Then, each pair of sonomicrometry crystals (2.0 mm; Sonometrics, Inc, London, ON, Canada) and each custom-made, fine-wire bipolar EMG hook electrode (0.5 mm bared tips with 2 mm spacing; California Fine Wire, Inc, Grover Beach, CA, USA) were implanted parallel to the fascicle axis of the mid-anterior region of sternobrachial portion of pectoralis at a depth of approximately 4 mm beneath the superficial fascia of the muscle. Before suturing (4–0 coated Vicryl) closed all of the incisions, the sonomicrometry signals were tested to ensure good signal quality. Then, all of the electrodes were sutured to superficial connective tissue a few millimeters away from the exit point of the pectoralis, to prevent movement independent of the fascicles. We left a small loop in the EMG electrode at the surface of the muscle to limit low-frequency noise. The electrode and transducer wires were tunneled subscutaneously along a narrow path from the ventral incision, along the lateral side of the dove, to the dorsal surface where the wires terminated into a back plug. The back plug was secured to the dove’s back by suturing the base of the plug to intervertebral ligaments, and it was custom-made (prior to surgery) using minature connectors and epoxy. The skin was sutured closed around the back plug. After the experiments, the doves were euthanized using an overdose of isoflurane. All training, surgical procedures, and experimental procedures were approved by Stanford’s Administrative Panel on Laboratory Animal Care (APLAC-30905 and APLAC-31426) and followed established methods and protocols (Biewener et al., 1998; Heers et al., 2016; Tobalske et al., 2003; Ellerby and Askew, 2007; Tobalske and Biewener, 2008; Jackson and Dial, 2011a; Robertson and Biewener, 2012).

Recordings were made by connecting a shielded cable (flexible wires sheathed approximately every 30 cm; diameter 7 mm) to the back plug on the dorsal side of the dove. We loosely suspended the cable above the dove (Figure 1—figure supplement 1), and the length of the cable that was suspended was approximately 76 mm, weighing 26.1 g. The cable was connected to an Ultrasound Dimension Gauge (UDG; Sonometrics, Inc, London, ON, Canada) to measure the sonomicrometry signals, and to a differential AC amplifier (Model 1700; A-M Systems, Sequim, WA, USA) to measure the EMG signals. All signals were recorded at 10,000 Hz using a Digidata 1550 A A/D converter (Axon Instruments, Union City, CA, USA). Sonomicrometry signals were converted into fiber lengths, ${\displaystyle L_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t},\mathrm{f}}}$, by calibrating at 0 mm, 5 mm, and 10 mm, and assuming the speed of sound transmission is 1590 m∙s^–1 when the muscle is at 37 degrees Celsius (Marsh, 2016). Then, these fiber lengths were converted into fiber strain, ${\displaystyle {\gamma }_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t},\mathrm{f}}}$, using Equation 4. We corrected for a time delay of 2ms (2.0% of the stroke) based on the time delay in the UDG, whereas the EMG signals did not have any time lag. To compare the EMG signals across all flights, we filtered the rectified signal with a second order butterworth filter with a cutoff frequency of 50 Hz. While we recorded EMG signals for both the left and right pectoralis muscles, the right signal was unreliable, so we only include data from the left EMG.

The 2D AFP setup (Lentink, 2018; Lentink et al., 2015) was validated by tethering a quadcopter to an instrumented beam and then comparing vertical and horizontal forces measured by the beam with forces measured by the AFP, as reported in (Lentink et al., 2015). Over a period of 10 seconds and a sampling rate of 1000 Hz, the impulse ratio and the mean force ratio of the 2D AFP to the beam were both 1.00±0.02 (n=10 trials) in the vertical direction and 1.00±0.01 in the horizontal direction (n=20 trials). Additionally, for a flight that starts and ends at rest, we expect that the total vertical impulse imparted by the legs and wings should equal full bodyweight support (Chin and Lentink, 2017), and that the net horizontal impulse should equal zero. Integrating the forces from takeoff to landing for the dove’s flight, we measured a vertical impulse ratio (impulse from legs and wings divided by impulse due to bodyweight, bw) of 1.03±0.02, and a horizontal impulse (impulse from legs and wings divided by bodyweight) of 0.02±0.02 bw∙s for 15 flights where the doves did not have an external cable attached to the back. For the 20 flights where the doves had a cable attached, we measured a vertical impulse ratio of 1.01±0.02, and a horizontal impulse of 0.08±0.02 bw∙s. The horizontal impulse is greater than zero when the cable is attached, because the dove did positive work to pull the cable forward.

To compare all flights and understand the underlying dynamics, we reduced the highly resolved 3D surface and kinematic measurements, in combination with the measured aerodynamic forces, into simplified aerodynamic and kinematic parameters which summarize the time-resolved state of entire wing.

To facilitate this simplification, we divided the wing into blade elements spaced 1 mm apart in the spanwise direction. For each blade element (J total blade elements, index j), we computed: (1) the surface area, ${\displaystyle S_{\mathrm{j}}}$, (2) the chordline vector from the trailing edge to the leading edge, (3) the velocity vector at the quarter-chord location (25% chord length behind the leading edge; measured relative to the Newtonian reference frame, N), ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{v}}_{}^{\mathrm{j}}}$, which includes the velocity components of both the flapping and morphing motion of the wing, and (4) the angle of attack induced by wing motion, ${\displaystyle {\alpha }_{\mathrm{j}}}$, which is the angle between the chordline and velocity vectors (Figure 2A and C).

We computed many of the simplified wing parameters using the spanwise blade element measurements and the 3D kinematics measurements. From the spanwise blade element measurements, we computed representative wing measurements which are the weighted averages of the blade element values (derived in Deetjen et al., 2020):

The weight for each blade element is the surface area times the velocity magnitude squared. To precisely track the position and orientation of the wing, we defined the stroke plane angle, ${\displaystyle {\varphi }_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{k}\mathrm{e}}}$, stroke angle, ${\displaystyle {\theta }_{\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{k}\mathrm{e}}}$, deviation angle, ${\displaystyle {\theta }_{\mathrm{d}\mathrm{e}\mathrm{v}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}}$, and twist angle, ${\displaystyle {\theta }_{\mathrm{t}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{t}}}$, of the wing (Figure 2A). The stroke plane was linearly fit using the path swept by the tip of the ninth primary feather across the full wingbeat, where a positive stroke plane angle corresponds to a stroke plane which is pitched down. The stroke angle and deviation angle are defined relative to the stroke plane according to Figure 2A, and the twist angle is the angle between the chordline vector at the root and the chordline vector along the span of the wing. Finally, we also computed time-resolved parameters to describe wing morphing. The semi-wing surface area, ${\displaystyle S_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$, is the area of a single wing:

which we used to calculate the aspect ratio, ${\displaystyle A{R}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$, of the wing:

where the wing radius, ${\displaystyle b_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$, is the length of the V-shaped path from the shoulder joint to the wrist joint to the ninth primary feather. The folding ratio, ${\displaystyle F{R}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$, of the wing:

depends on the span vector, ${\displaystyle {\overrightarrow{\mathbf{r}}}_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}}}$, which is the vector from the shoulder joint to the ninth primary feather. The max function finds the maximum value across the stroke.

To compute the lift and drag forces and the lateral (right to left) aerodynamic forces produced by the wing (Figure 1C), we utilized the measured vertical and horizontal aerodynamic forces, in addition to the 3D surface and kinematic measurements. From the derivation in Hedrick et al., 2003 (Deetjen et al., 2020), the directions of lift, ${\displaystyle {\hat{\mathbf{L}}}_{\mathrm{L}}}$, and drag, ${\displaystyle {\hat{\mathbf{D}}}_{\mathrm{L}}}$, for the entire left wing are:

To compute the magnitude of lift, $L$, and drag, $D$, on both wings, we included the net measured horizontal, ${\displaystyle F_{\mathrm{A}\mathrm{F}\mathrm{P},\mathrm{x}}}$, and vertical, ${\displaystyle F_{\mathrm{A}\mathrm{F}\mathrm{P},\mathrm{z}}}$, aerodynamic forces:

where the x and z subscripts indicate the components of the lift and drag unit vectors, and $\theta$ is the yaw angle of the dove. A positive yaw angle corresponds with the dove yawing to the left. Finally, we summed the lift and drag vectors to compute the total 3D aerodynamic force vector produced by the wing, ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$:

Near stroke reversal, these equations encounter a singularity, so we smoothed the angle of attack, and lift and drag forces in these regions using the `perfect smoother' (Eilers, 2003) with smoothing weights, $W$, based on the same approach as in Deetjen et al., 2020:

with $c_0=0.15$ and $c_1=0.35$.

To compute the time-resolved power that the dove’s muscles must have produced, we balanced the muscle power with the effects of the aerodynamic forces and inertia.

To do this, we first needed to decide what the most effective scope to analyze power is: power balance of a single wing of the dove, or power balance of the entire dove. To aid in this decision, we examined the definition of power, ${\displaystyle {}_{}{}^{\mathrm{N}}P_{}^{\mathrm{B}}}$, for a system B (a collection of point masses and rigid bodies: scope of the power analysis) in the Newtonian reference frame N:

where ${\displaystyle {}_{}{}^{\mathrm{N}}K_{}^{\mathrm{B}}}$ is the kinetic energy of system B in reference frame N. Specifying what system B encompasses was needed to clarify how terms relating to the muscles, aerodynamics, and inertia would be incorporated into this high-level equation. One possible way to define system B was as a system that contains only point masses and rigid bodies which model a single wing of the dove. This definition is advantageous because it eliminates the term related to the kinetic energy of the body, which involves the second derivative of the position of the body. However, the disadvantage of this definition is that an extra power term appeared which we could not measure. That is, we had no method to measure the reaction force on the shoulder joint, which contributed to power because it was an external force acting on the wing. The only alternative way to eliminate this term was to assume that the shoulder joint itself moved in a Newtonian reference frame (constant velocity and no rotation), which is a large over-simplification. Instead, we defined system B to contain the entire dove, and accepted the difficulty in computing the second derivative of the position of the body, which could add noise to the results. With the scope of the power balance decided, we expanded each side of Equation S13 to incorporate the effects of the muscles, aerodynamics, and inertia.

The kinetic energy on the right side of Equation S13 is made up of the distributed masses of the wings and body. This term expands to:

where ${\displaystyle m_{\mathrm{i}}}$ is the mass of the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ of ${\displaystyle {\mathbf{I}}_{\mathrm{W},\mathrm{B}}}$ point masses in the body and wings, ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{v}}_{}^{\mathrm{i}}}$ is the velocity vector of the center of mass of the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ point mass in the Newtonian reference frame, ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathit{\omega }}_{}^{{\mathrm{B}}_{\mathrm{i}}}}$ is the angular velocity of the reference frame of the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ body in system B relative to the Newtonian reference frame, and ${\displaystyle {\overrightarrow{\overrightarrow{\mathbf{I}}}}^{{\mathrm{B}}_{\mathrm{i}}}}$ is the moment of inertia of the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ body in system B relative to the center of mass of that body. Since we modeled the dove using point masses, Equation S14 simplifies to:

The power term on the left side of Equation S13 is made up of the external power exerted on the dove and the internal power produced by the dove. Assuming that all the internal power produced by the dove was accounted for by the two pectorales and supracoracoideuses, the dove’s power expands to:

where ${\displaystyle -{}_{}{}^{\mathrm{N}}P_{}^{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}}$ is the power done by external aerodynamic forces, ${\displaystyle -{}_{}{}^{\mathrm{N}}P_{}^{\mathrm{g}}}$ is the power done by external gravitational forces, ${\displaystyle {}_{}{}^{\mathrm{N}}P_{}^{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$ is the internal power produced by each pectoralis muscle, and ${\displaystyle {}_{}{}^{\mathrm{N}}P_{}^{\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{r}\mathrm{a}}}$ is the internal power produced by each supracoracoideus muscle. We derived the aerodynamic power in Deetjen et al., 2020:

where ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}}$ is the overall aerodynamic force on the wing, and ${\displaystyle {\overrightarrow{\mathbf{v}}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$ is the representative wing velocity vector computed using , Equation S1 . The gravitational power is the summation of the gravity forces on each point mass dotted with their velocity vectors:

where ${\displaystyle \overrightarrow{\mathbf{g}}}$ is the gravity vector. This equation can be simplified, showing that only the velocity in the z direction is relevant for gravitational power:

where $g$ is the magnitude of gravity and ${\displaystyle {}_{}{}^{\mathrm{N}}v_{\mathrm{z}}^{\mathrm{i}}}$ is the vertical component of the velocity vector.

To compare the effects on required muscle power from terms related to inertia versus aerodynamics, we grouped all of the terms related to inertia together, and named the resulting term inertial power:

Hence, the total muscle power simply equals the inertial power plus the aerodynamic power:

which could be computed using the aerodynamic forces and kinematics that we measured. During the mid-downstroke, the pectoralis primarily produces the power, so we eliminated the supracoracoideus term.

In order to convert the muscle power into force and stress, we incorporated the strain measured in vivo and other parameters measured during dissection. The power produced by the pectoralis equals the total force generated by the pectoralis measured at the attachment point, ${\displaystyle F_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$, times the rate of change of the length of the entire pectoralis, ${\displaystyle {\dot{L}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$:

Since the pectoralis is not uniform along its length, we computed the rate of change of its length using measurements taken from the in vivo experiment and during the dissection:

where ${\displaystyle {\overline{\alpha }}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t},\mathrm{f}}}$ is the average measured muscle fiber angle (dissection), ${\displaystyle {\overline{L}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t},\mathrm{f}}}$ is the average measured muscle fiber length (dissection), and ${\displaystyle {\dot{\gamma }}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t},\mathrm{f}}}$ is the measured strain rate of the muscle fibers assuming uniform strain throughout the muscle (in vivo). To convert the muscle force to the average muscle stress, ${\displaystyle {\sigma }_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$, we calculated the Physiological Cross-Sectional Area (PCSA) (Alexander, 1983), which is the area of the cross section of the pectoralis perpendicular to its fibers:

For muscle density, we use ${\displaystyle {\rho }_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}=1060}$ kg∙m^–3 (Tobalske and Biewener, 2008), and ${\displaystyle m_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$ is the average mass of the left and right pectorales of each dove measured during dissection.

To compute the time-resolved 3D moment vector (torque) that the dove’s muscles must have produced, we balanced the muscle moment with the effects of the aerodynamic forces and inertia. We used this result, together with the pectoralis force magnitude, computed by balancing 1D power, to compute the pull angle of the pectoralis on the humerus during downstroke.

Like the power balance derivation, we needed to decide what the most effective scope to analyze angular momentum is: a single wing of the dove, or the entire dove. We decided to analyze the angular momentum of a single wing of the dove because fortunately, whereas the reaction force at the shoulder joint is problematic for the power balance, it cancels out for the angular momentum balance. The base form of the 3D angular momentum balance equation is:

where ${\displaystyle {\overrightarrow{\mathbf{M}}}^{\mathrm{W}/\mathrm{O}}}$ is the sum of the moments on the wing (W: a system of particles and rigid bodies) about the shoulder joint (point O: fixed on the wing), ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{H}}_{}^{\mathrm{W}/\mathrm{O}}}$ is the angular momentum of the wing about the shoulder joint in the Newtonian reference frame, ${\displaystyle \frac{{}_{}{}^{\mathrm{N}}d_{}^{}}{dt}}$ is the time derivative in the Newtonian reference frame, ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{v}}_{}^{\mathrm{O}}}$ is the velocity of the shoulder joint in the Newtonian reference frame, and ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{L}}_{}^{\mathrm{W}}}$ is the linear momentum of the wing in the Newtonian reference frame. Reaction forces at the shoulder joint do not create any moment about that shoulder joint, so we canceled that term. We eliminated any reaction torques by modelling the shoulder joint as a freely-rotating ball and socket joint with negligible friction. We also assumed negligible internal wing forces, such as muscles which bend the wrist or elbow. With the scope of the angular momentum balance decided, we expanded each side of Equation S27 to incorporate the effects of the muscles, aerodynamics, and inertia.

The terms on the right side of Equation S27 all relate to the distributed mass of the wing. The linear momentum of the wing in the Newtonian reference frame is:

where ${\displaystyle m_{\mathrm{i}}}$ is the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ of ${\displaystyle {\mathrm{I}}_{\mathrm{W}}}$ bodies (or point masses) on the wing, and ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{v}}_{}^{\mathrm{i}}}$ is the velocity vector of the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ body in the Newtonian reference frame. The angular momentum of the wing about the shoulder joint in the Newtonian reference frame is:

where ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathit{\omega }}_{}^{{\mathrm{W}}_{\mathrm{i}}}}$ is the angular velocity of the reference frame of the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ body in system W (the wing) relative to the Newtonian reference frame, ${\displaystyle {\overrightarrow{\overrightarrow{\mathbf{I}}}}^{{\mathrm{W}}_{\mathrm{i}}}}$ is the moment of inertia of the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ body in system W relative to the center of mass of that body, and ${\displaystyle {\overrightarrow{\mathbf{r}}}^{\mathrm{i}/\mathrm{O}}}$ is the position vector from the shoulder joint to the ${\displaystyle {\mathrm{i}}^{\text{th}}}$ body in system W. Because we modeled the dove’s wings using point masses, we simplified to:

Next, we expanded the left-hand side of Equation S27 which is the sum of the moments on the wing. This is made up of aerodynamic, gravitational, and muscle torques:

where ${\displaystyle {\overrightarrow{\mathbf{M}}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}^{\mathrm{W}/\mathrm{O}}}$ is the aerodynamic moment, ${\displaystyle {\overrightarrow{\mathbf{M}}}_{\mathrm{g}}^{\mathrm{W}/\mathrm{O}}}$ is the gravitational moment, ${\displaystyle {\overrightarrow{\mathbf{M}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}^{\mathrm{W}/\mathrm{O}}}$ is the pectoralis moment, and ${\displaystyle {\overrightarrow{\mathbf{M}}}_{\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{r}\mathrm{a}}^{\mathrm{W}/\mathrm{O}}}$ is the supracoracoideus moment. Each of these terms equals the cross product of the vector from the shoulder joint to the location where the force acts crossed with the force vector. Hence, the gravitational moment term can be written as:

This concludes the expansion of all terms related to mass in the angular momentum balance, and to summarize them, we defined a parameter called the inertial moment:

By grouping terms in this manner, the total muscle moment equals the inertial moment plus the aerodynamic moment:

similar to Equation S22 for the power balance.

The aerodynamic moment was challenging to compute because aerodynamic forces acted continuously along the wing. In order to solve for the aerodynamic moment simply in terms of the aerodynamic force on the wing and known parameters, we modeled the wing using blade elements:

where ${\displaystyle {\overrightarrow{\mathbf{r}}}^{\mathrm{j}/\mathrm{O}}}$ is the position vector from the shoulder joint to the quarter-chord location (a quarter of the way from the leading edge to the trailing edge) on the ${\displaystyle {\mathrm{j}}^{\text{th}}}$ of J blade elements, and ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o},\mathrm{j}}}$ is the aerodynamic force on the ${\displaystyle {\mathrm{j}}^{\text{th}}}$ blade element. We expanded ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o},\mathrm{j}}}$ to:

where the ${\displaystyle {\mathrm{j}}^{\text{th}}}$ blade element has surface area, ${\displaystyle S_{\mathrm{j}}}$, lift coefficient, ${\displaystyle C_{\mathrm{L},\mathrm{j}}}$, lift direction, ${\displaystyle {\hat{\mathbf{L}}}_{\mathrm{j}}}$, drag coefficient, ${\displaystyle C_{\mathrm{D},\mathrm{j}}}$, drag direction, ${\displaystyle {\hat{\mathbf{D}}}_{\mathrm{j}}}$, and velocity ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{v}}_{}^{\mathrm{j}}}$. Assuming that the lift and drag coefficients and directions are spanwise invariant (${\displaystyle C_{\mathrm{L},\mathrm{j}}=C_{\mathrm{L}}}$, ${\displaystyle {\hat{\mathbf{L}}}_{\mathrm{j}}=\hat{\mathbf{L}}}$, ${\displaystyle C_{\mathrm{D},\mathrm{j}}=C_{\mathrm{D}}}$, and ${\displaystyle {\hat{\mathbf{D}}}_{\mathrm{j}}=\hat{\mathbf{D}}}$ for all j), this reduces to:

and further simplifies to:

In order to incorporate the total measured aerodynamic force on the wing, ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}}$ into the equation for the aerodynamic moment, we simplified the aerodynamic moment to:

where ${\displaystyle {\overrightarrow{\mathbf{R}}}^{\mathrm{A}\mathrm{M}/\mathrm{O}}}$ is the position vector from the shoulder joint to the aerodynamic moment center, and ${\displaystyle {\overrightarrow{\mathbf{T}}}_{\mathrm{O}}}$ accounts for extra torque due to paired, non-aligned forces. We assumed ${\displaystyle \left|{\overrightarrow{\mathbf{T}}}_{\mathrm{O}}\right|=0}$ based on the assumption that all of the lift and drag forces point in the same direction, and primarily act at locations along the wing that form a straight line starting at the shoulder joint. Using blade elements, we expanded ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}}$, giving:

and then simplified by assuming spanwise invariance in lift and drag coefficients and direction:

By comparing Equations S39 and S42, we solved for ${\displaystyle {\overrightarrow{\mathbf{R}}}^{\mathrm{A}\mathrm{M}/\mathrm{O}}}$:

The full expansion of the aerodynamic moment is:

which closely resembles Equation S17 for aerodynamic power. By incorporating this equation for the aerodynamic moment, Equation S35 can be used to solve for the moment that the flight muscles together must generate.

In order to compute the angle that the pectoralis pulls on the humerus, we used the computed pectoralis force magnitude from the 1D power balance and the muscle moment from the 3D angular momentum balance. We started by expanding the pectoralis moment:

where ${\displaystyle {\overrightarrow{\mathbf{r}}}^{\mathrm{P}/\mathrm{O}}}$ is the vector from the shoulder joint to the deltopectoral crest where the pectoralis pulls on the humerus, and where ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$ is the total pectoralis force. While the directions of ${\displaystyle {\overrightarrow{\mathbf{r}}}^{\mathrm{P}/\mathrm{O}}}$ and ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$ were unknown and changed continuously through the stroke, we knew the magnitudes of both: ${\displaystyle \left|{\overrightarrow{\mathbf{r}}}^{\mathrm{P}/\mathrm{O}}\right|}$ from dissection and ${\displaystyle \left|{\overrightarrow{\mathbf{F}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}\right|}$ from the 1D power balance. Hence, we used the following equation, derived from Equation S46, to solve for the effective pull ratio, ${\displaystyle \sin {\theta }_{\mathrm{p}}}$, and the pectoralis pull angle on the humerus, ${\displaystyle {\theta }_{\mathrm{p}}}$:

The next section details how we used this result to compute the full 3D pectoralis force vector.

To determine the direction in which the pectoralis pulls on the humerus, we overlaid the skeleton of a dove on the measured 3D surface (Figure 6, Figure 6—figure supplement 1). Using Computed Tomography (CT), we scanned a dove skeleton of similar mass (170.2 g) to the four doves in this study, and then scaled the dimensions for each of the four doves to match the size of the dove we CT scanned. The scaling factor was the one third power of the ratio of each of the four dove’s masses to the CT-scanned dove’s mass. We assumed that the bones in the body of the dove (ribs, sternum, furcula, coracoid, and scapula) do not move relative to each other during flight. In reality, they do move, but not substantially compared to the movement of the wings (Baier et al., 2013; Heers et al., 2016). We also assumed that the bones in the wing consist of three rigid bone segments connected to each other: (1) the brachium (humerus), (2) the antebrachium (radius and ulna), and (3) the manus (radiale +ulnare + carpometacarpus +phalanges). The positions and orientations of the reference frames (joint axes) associated with each bone segment were determined based on the inertial axes of the bones and the joint anatomy (Baier et al., 2013; Heers et al., 2016). We assumed that all wing joints were free to rotate with 3 degrees of freedom and no translational degrees of freedom (measurements show that translation of the joints relative to the reference frame is minimal for modeling purposes) (Baier et al., 2013; Heers et al., 2016). The positioning of each joint was determined based on the centers of rotation (Stowers, 2017) when previously-measured kinematics in chukars (Alectoris chukar; kinematics measured using X-Ray Reconstruction of Moving Morphology, i.e., XROMM) (Heers et al., 2016) were applied to the dove model. Muscle paths were modeled based on dissection of the CT-scanned dove. Pins were pushed into the body through the pectoralis and supracoracoideus in multiple locations along the central tendons, and pictures were taken in situ for multiple wing positions. Then the pectoralis and supracoracoideus were removed from the body, and the positions of the pins were re-photographed to visualize the paths of the tendons with respect to the skeleton. These paths were used to construct a simplified musculoskeletal model in SIMM (Software for Interactive Musculoskeletal Modeling), using wrapping surfaces to prevent the muscles from passing through bone.

We overlaid the musculoskeletal model, made up of four skeletal sections (one for the body and three for the left wing) onto the 3D surface of the dove during each flight in four stages.

1. First, we determined the position and orientation of the skeleton relative to the 3D surface of the body. We matched the orientation of the vertebral column with the pitch and yaw angles of 3D surface fit of the body (assuming zero roll), and we matched the plane of symmetry of the skeletal and 3D surface bodies. We positioned the skeletal shoulder joint as close as possible to the average position (relative to the body of the dove) of the manually tracked shoulder joint, subject to the previous constraints.

2. Second, we rotated the humerus about the shoulder joint (rigidly attached to the dove’s body) and the radius/ulna about the elbow joint (rigidly attached to the end of the humerus) until the end of the radius/ulna (the wrist joint) best matched the manually tracked wrist of the flying doves. Since rotation about the long axis of the humerus and the long axis of the radius/ulna do not impact the final position of the wrist from a purely geometric perspective (in reality, long axis rotation of the humerus partially determines how much abduction or adduction occurs at the elbow, which in turn determines the position of the wrist), this left two rotational axes of consequence for each bone segment: hence four rotational parameters to fit. Constraining the wrist joint to match between the skeletal and 3D surface models generates three constraints. For the fourth constraint, we rotated the humerus and radius/ulna so the deltopectoral crest and dorsal surfaces of the radius and ulna were approximately in the plane of the 3D surface of the wing. We then solved for these four rotational angles using a nonlinear least-squares regression, and smoothed the result (Eilers, 2003), to find the orientations of the humerus and radius/ulna.

3. Third, we oriented the manus such that the most distal phalanx pointed toward the ninth primary wingtip (which was tracked).

   4. Fourth, we allowed the humerus to rotate in the long axis so that the moment generated by the pectoralis best matched the required muscle moment, ${\displaystyle {\overrightarrow{\mathbf{M}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}^{\mathrm{W}/\mathrm{O}}}$ when the pectoralis was shortening.

Using the overlaid skeleton, the pectoralis force magnitude, the pectoralis moment, and the pectoralis pull angle, we computed the time-resolved 3D force vector of the pectoralis. Examining Equation S46, the overlaid skeleton provides the missing information necessary to determine the vector from the shoulder joint to the deltopectoral crest where the pectoralis pulls on the humerus, ${\displaystyle {\overrightarrow{\mathbf{r}}}^{\mathrm{P}/\mathrm{O}}}$. With this information, we solved for the direction of the pectoralis force, ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}/\left|{\overrightarrow{\mathbf{F}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}\right|}$:

where ${\displaystyle R\left(\mathbf{a},\theta \right)}$ is a function which outputs a rotation matrix corresponding to a rotation about axis $a$ by an angle of $\theta$. In other words, the direction of ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$ is found by rotating the negative direction of ${\displaystyle {\overrightarrow{\mathbf{r}}}^{\mathrm{P}/\mathrm{O}}}$ about the direction of ${\displaystyle {\overrightarrow{\mathbf{M}}}_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}^{\mathrm{W}/\mathrm{O}}}$ by an angle of ${\displaystyle -{\theta }_{\mathrm{p}}}$. Note that in order to compute the direction of ${\displaystyle {\overrightarrow{\mathbf{r}}}^{\mathrm{P}/\mathrm{O}}}$, only the orientation of the humerus is needed. The other skeletal segments (radius/ulna and manus) are computed for visualization purposes in the figures.

To address the modeling challenge of disentangling pectoralis and supracoracoideus muscle power during stroke-reversal, we considered the case where there is no elastic energy storage, and the muscles only generate or absorb power. Based on the sonomicrometry measurements of the pectoralis and total computed muscle power, each time-step during the stroke fell into one of four categories: (1) Positive total muscle power while the pectoralis was shortening. This was the case for much of the downstroke and indicated that the pectoralis was generating power. However, near stroke reversal, low pectoralis muscle velocity prevented high levels of power generation (see Equation S23). Hence, we restricted the maximum force production of pectoralis in stroke-reversal regions based on the maximum force produced by the pectoralis during the mid-downstroke. Any leftover total power is categorized as supracoracoideus power generation. (2) Negative total muscle power while the pectoralis was shortening. This was the case at the end of the upstroke and indicated that the supracoracoideus was acting as a brake by absorbing power. (3) Positive total muscle power while the pectoralis was lengthening. This was the case for much of the upstroke and indicated that the supracoracoideus was generating power. (4) Negative total muscle power while the pectoralis was lengthening. This was the case during the latter half of the upstroke when the wing was decelerating and indicated that the pectoralis was acting as a brake and absorbing power. However, similar to case 1, if pectoralis muscle velocity was too low to absorb the entirety of the total muscle power, we attributed the remainder to supracoracoideus power absorption.

To model energy storage in the supracoracoideus tendon, we considered how different amounts and timing of elastic energy storage in the tendon would impact the power distribution between the pectoralis and supracoracoideus. The two tuning parameters we varied are: (1) the proportion of supracoracoideus power which was elastically stored and released, and (2) the time period over which energy was stored in the supracoracoideus tendon. When we increased the first tuning parameter, we proportionally shifted the distribution of supracoracoideus power from generated to released (upstroke), and from absorbed to stored (downstroke). The power generated by the pectoralis was correspondingly increased to ensure sufficient energy for the supracoracoideus tendon to store during the downstroke. To account for hysteresis effects (Tobalske and Biewener, 2008) in the supracoracoideus tendon, the integrated extra power produced by the pectoralis must exceed the energy released by the supracoracoideus tendon by 7%. In order to estimate the temporal distribution of energy storage in the supracoracoideus, we assumed that the supracoracoideus lengthens according to measurements of ascending flight in pigeons (Figure 2 of Tobalske and Biewener, 2008). Specifically, we matched the timing of the pectoralis strain between our measurements and their measurements, in order to match the supracoracoideus strain and find when the supracoracoideus finishes lengthening. This timing, combined with our timing tuning parameter, defined the region of the stroke over which the supracoracoideus stores energy. Next, we assumed that the shape of the force-length curve of the supracoracoideus tendon matched the curve given in Millard et al., 2013. With this timing, strain, and force shape information, we calculated the shape of the power curve for the supracoracoideus tendon as a function of stroke percentage. We then scaled the magnitude of the power curve according to the amount of energy the supracoracoideus needed to store. By varying these two tuning parameters, we examined the effects of elastic energy storage in the supracoracoideus tendon on muscle performance.

We scaled our aerodynamic, inertial, and muscle measurements for doves across multiple extant birds, to estimate comparative patterns of power output during slow flight.

For each extant bird analyzed, we scaled our results for doves by modifying all relevant variables (see Supplementary file 2 and Supplementary file 3 for details), and plugging them directly into the relevant equations. We scaled the following variables based on extant data of 27 birds measured in literature (Berg and Rayner, 1995): body mass, ${\displaystyle m_{\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}}}$, wing mass (single wing), ${\displaystyle m_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$, wing moment of inertia (single wing), ${\displaystyle I_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$, wingspan (distance from shoulder joint to wingtip; single wing), ${\displaystyle \left|{\overrightarrow{\mathbf{r}}}_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}}\right|}$, distance from the shoulder joint to the center of gravity of the wing (single wing), ${\displaystyle \left|{\overrightarrow{\mathbf{r}}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g},\mathrm{c}\mathrm{g}}\right|}$, wing area (single wing), ${\displaystyle S_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}}$, and wingbeat frequency, $f$. For four additional birds where data was missing, we used data from other sources (Jackson and Dial, 2011a; Greenewalt, 1962), scaled isometrically to match the body mass. For pectoralis mass (single pectoralis), ${\displaystyle m_{\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}}}$, we isometrically scaled data in literature to match the body mass (Tobalske and Biewener, 2008; Tobalske et al., 2005; Greenewalt, 1962). Finally, we assumed the following relationships for the following variables: time step for integration, $\Delta t$:

aerodynamic force, ${\displaystyle {\overrightarrow{\mathbf{F}}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}}$:

wing velocity (distributed), ${\displaystyle {\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{v}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}}^{\mathrm{i}}}}$:

and body velocity, ${\displaystyle {}_{}{}^{\mathrm{N}}\overrightarrow{\mathbf{v}}_{\mathrm{b}\mathrm{o}\mathrm{d}\mathrm{y}}^{}}$ equivalent to the doves. We calculated the point mass distribution in the wing, ${\displaystyle m_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g},\mathrm{i}}}$, using a least-squares fit to match wing mass, inertia, and center of gravity data. Additionally, to compare pectoralis electrical activation timing to our scaling parameters in Figure 7D, we used the same procedure for 17 birds (Supplementary file 2b and Supplementary file 3b-c) (Tobalske and Biewener, 2008; Biewener et al., 1992; Jackson and Dial, 2011a; Dial and Biewener, 1993; Hedrick et al., 2003; Tobalske et al., 2005; Berg and Rayner, 1995; Tobalske et al., 2017; Altshuler et al., 2010; Tobalske et al., 2010; Hedrick et al., 2004; Ingersoll and Lentink, 2018; Jackson et al., 2011b; Donovan et al., 2013; Greenewalt, 1962; Dial et al., 1997).

To interpret the scaling relationships, we traced the effects of scaling parameters in the aerodynamic power and moment equations. For aerodynamic power, we started with Equations S1 and S17 combined for some bird:

Next, for each variable, we replaced it with an equivalent expression by utilizing the scaling relationships:

where variables with a tilde represent a second bird. When we reduced this expression, we recovered the equation for aerodynamic power for the second bird:

which led us to the scaling relationship for aerodynamic power:

We then rearranged this equation to form interpretable parameters: