## ABSTRACT

Flapping insect flight is a complex and beautiful phenomenon that relies on fast, active control mechanisms to counter aerodynamic instability. To directly investigate how freely flying *Drosophila* *melanogaster* control their body pitch angle against such instability, we perturbed them using impulsive mechanical torques and filmed their corrective maneuvers with high-speed video. Combining experimental observations and numerical simulation, we found that flies correct for pitch deflections of up to 40 deg in 29±8 ms by bilaterally modulating their wings' front-most stroke angle in a manner well described by a linear proportional–integral (PI) controller. Flies initiate this corrective process only 10±2 ms after the perturbation onset, indicating that pitch stabilization involves a fast reflex response. Remarkably, flies can also correct for very large-amplitude pitch perturbations – greater than 150 deg – providing a regime in which to probe the limits of the linear-response framework. Together with previous studies regarding yaw and roll control, our results on pitch show that flies' stabilization of each of these body angles is consistent with PI control.

## INTRODUCTION

From walking humans to flying insects, many fascinating forms of bio-locomotion are contingent upon robust stabilization control. Implementing this control is particularly difficult in the case of small, flapping-wing insects, as flapping flight is inherently subject to rapidly divergent aerodynamic instabilities (Faruque and Humbert, 2010; Gao et al., 2011; Liu et al., 2010; Pérez-Arancibia et al., 2011; Ristroph et al., 2013; Sun, 2014; Sun et al., 2007; Sun and Xiong, 2005; Taylor and Thomas, 2003; Taylor and Żbikowski, 2005; Xu and Sun, 2013; Zhang and Sun, 2011, 2010). As such, flying insects have evolved stabilization techniques relying on reflexes that are among the fastest in the animal kingdom (Beatus et al., 2015) and robust to the complex environment that insects must navigate (Combes and Dudley, 2009; Dickerson et al., 2012; Ortega-Jimenez et al., 2013; Ravi et al., 2013; Vance et al., 2013).

In particular, pitching instability is a prominent obstacle for flight control in flapping insects. Analytical and numerical modeling suggest that, for many two-winged insects (e.g. flies), periodic flapping couples with longitudinal body motion to produce rapidly growing oscillations of the body pitch angle (Chang and Wang, 2014; Ristroph et al., 2013; Sun, 2014; Sun et al., 2007). This oscillatory instability can be understood as the result of differential drag on the wings due to longitudinal body motion (Ristroph et al., 2013; Sun et al., 2007). For example, if a fly pitches down while hovering, its re-directed lift propels its body forward, causing an increased drag on the wings during the forward stroke relative to the backward stroke. Because the centers of pressure of the fly's wings are always above the body center of mass during normal flapping, this drag asymmetry generates a torque that pitches the fly up. Rather than acting as a restoring torque, the drag – together with the body inertia – pitches the fly up, beyond its initial pitch orientation. The fly then begins to move backwards, and oscillation ensues in the opposite direction. This mechanism results in an undulating instability of the body pitch angle, which doubles over a time scale of ∼9 wingbeats (Sun, 2014). Mitigating the effects of this instability requires flies to actively adjust their wing motion on time scales faster than the growth of these oscillations.

Our work builds upon an already rich corpus of literature on insect flight control, a sizable portion of which addresses the pitch degree of freedom. Experimental studies subjecting tethered insects to both mechanical pitching perturbations and visual pitching stimuli (Dickinson, 1999; Nalbach, 1994; Sherman and Dickinson, 2004, 2003; Taylor and Thomas, 2003; Zanker, 1990) have elucidated stereotyped kinematic responses for pitch correction, including manipulation of wingstroke angle, stroke plane orientation, wingbeat frequency and body configuration. However, tethered insects do not constitute a closed-loop feedback system, as changes to their wing kinematics do not affect their body orientation (Roth et al., 2014). Moreover, in the case of tethered flies, it has been shown that the wing kinematics are qualitatively different from those in free flight (Bender and Dickinson, 2006; Fry et al., 2005). Thus, free-flight studies are necessary for a comprehensive understanding of pitch control. Significant analysis has been performed on freely flying insects executing voluntary maneuvers (Bergou et al., 2010; Ennos, 1989; Fry et al., 2003; Ristroph et al., 2011, 2009) or responding to visual stimuli (Cheng et al., 2011; Muijres et al., 2014; Tammero and Dickinson, 2002; Windsor et al., 2014), but the general challenge of systematically inducing mechanical perturbations on untethered insects has traditionally been a barrier to the study of stabilization reflexes. Some notable exceptions to this include methods of mechanical perturbation using air-flow vortices (Combes and Dudley, 2009; Ortega-Jimenez et al., 2013; Ravi et al., 2013) or gusts of wind (Vance et al., 2013). However, such fluid-impulse methods are difficult to tune, and are thus not ideal for inducing the fast, precise mechanical perturbations that are required for a quantitative understanding of pitch control.

To achieve the necessary speed and precision for measurements of body pitch control, we used a perturbation scheme that has previously been applied to analyzing control of the yaw (Ristroph et al., 2010) and roll (Beatus et al., 2015) degrees of freedom. We mechanically perturbed free-flying *Drosophila* *melanogaster* by gluing small magnetic pins to their dorsal thoracic surfaces and applying short bursts (5–8 ms) of a vertical magnetic field that pitched their bodies up or down. As the flies corrected their orientation, we measured their body and wing kinematics using high-speed video (Fig. 1A).

#### List of symbols and abbreviations

*C*_{D}- coefficient of drag
*C*_{L}- coefficient of lift
- CI
- confidence interval
- CoM
- center of mass
- CoP
- center of pressure
*F*_{D}- drag force
*F*_{L}- lift force
*K*_{i}- PI controller integral gain
*K*_{p}- PI controller proportional gain
- PI
- proportional–integral
- second moment of wing area
- RMS
- root mean square
*S*- wing area
*U*_{t}- wing tip velocity (lab frame)
- α
- wing angle of attack
- Δ
*T* - controller latency time
- η
_{w} - wing pitch angle
- θ
_{b} - body pitch angle
- θ
_{w} - wing elevation angle
- ρ
_{0} - air density
- ρ
_{b} - body roll angle
- τ
_{p} - pitch torque
- τ
_{w} - wing torque
- front wingstroke angle
- back wingstroke angle
- body yaw angle
- ϕ
_{w} - wingstroke angle

We recorded perturbation events with amplitudes typically ranging from 5 to 40 deg for both pitching up and pitching down. For these perturbations, flies recover 90% of their pitch orientation within ∼30 ms. Moreover, we found that the corrective process is initiated ∼2 wingbeats (∼10 ms) after the onset of the impulsive torque; such a short latency time indicates that this corrective process is a reflexive behavior largely governed by input from the halteres, the flies' rate-gyro-like mechanical sensing organs (Dickinson, 1999; Nalbach, 1994; Pringle, 1948).

To generate corrective pitching torques, flies bilaterally modulate their wings' front-most stroke angle, i.e. they flap their wings more or less in the front to pitch up or down. This corrective mechanism is in general agreement with previous findings on active body pitch stabilization in fruit flies (Chang and Wang, 2014; Dickinson, 1999; Taylor, 2001; Zanker, 1990). We show that flies' modulation of front stroke angle is well modeled by a linear, continuous, proportional–integral (PI) controller with a time delay, Δ*T=*6±1.7 ms (mean±s.d.). Such simple feedback controller models have been successfully used to describe a wide variety of sensorimotor behavior in many species, including flies (Beatus et al., 2015; Elzinga et al., 2012; Ristroph et al., 2010), hawkmoths (Cheng et al., 2011; Dyhr et al., 2013), electric fish (Cowan and Fortune, 2007; Madhav et al., 2013; Roth et al., 2011), pigeons (Lin et al., 2014) and cockroaches (Cowan et al., 2006; Lee et al., 2008; for comprehensive reviews, see Cowan et al., 2014; Roth et al., 2014; Tytell et al., 2011). Our results indicate that pitch control in fruit flies is an extremely fast and robust process, which can be accurately modeled by a simple controller for a wide range of perturbation amplitudes. Moreover, we found that flies are capable of correcting for pitch deflections of 150 deg or more, a perturbation regime in which the linear controller theory begins to break down. Together with previous results on how flies control yaw (Ristroph et al., 2010) and roll (Beatus et al., 2015), the analysis of pitch control presented here addresses a missing piece in our understanding of how flies control each of their body angles individually.

## MATERIALS AND METHODS

### Fly preparation

We performed each experiment using common fruit flies (*D. melanogaster* Meigen, females) from an out-bred laboratory stock. Individual flies were anesthetized at 0–4°C, at which point we carefully glued 1.5–2 mm long, 0.15 mm diameter ferromagnetic pins to their notum (dorsal thoracic surface), oriented so that the pin lay in the fly's sagittal plane. The pin is shown in Fig. 1A (false-colored blue) and Fig. 2 (images). Control experiments with untreated flies showed that the addition of the pin did not qualitatively alter flies' flight kinematics. When attached, the pins added ∼20% (∼0.2 mg) to the mass and pitch moment of inertia of the fly, which falls within the range of their natural body mass variations. Moreover, the pin contributed negligibly to off-diagonal components of the flies' inertia tensors, and therefore did not introduce any coupling between the rotational degrees of freedom of the body. The primary effect of the added pin mass was a dorsal shift of ∼0.2 mm in the flies' center of mass, which is accounted for in the calculation of aerodynamic torque. Following Card and Dickinson (2008), the center of mass for untreated flies was assumed to be located at the centroid of the fly body, halfway between the head and tail along the body axis.

### Videography and mechanical perturbation

Once 15–30 flies had been prepared as above, we released them into a transparent cubic filming chamber of side length 13 cm. On the top and bottom of the chamber were attached two horizontally oriented Helmholtz coils, which produced a vertical magnetic field. The central region of the chamber was filmed by three orthogonal high-speed cameras (Phantom V7.1) at 8000 frames s^{−1}. The cameras were calibrated using a direct linear transformation scheme detailed in Lourakis and Argyros (2009) and Theriault et al. (2014). When flies entered the filming volume, an optical trigger simultaneously signaled the cameras to record and supplied a 5–8 ms current pulse to the Helmholtz coils. We varied both the strength and duration of the magnetic pulse produced by the coils across experiments. Maximal field strengths reached ∼10^{−2} T, and most experiments were performed with a pulse that lasted 5.8 ms. The magnetic pulse exerted a torque on the ferromagnetic pin, pitching the fly either up or down.

Of the movies collected from several independent fly batches using the above method, we selected 18 to analyze in full; 16 additional movies were partially analyzed to collect more data on pitch correction time and to observe correction for very large-amplitude perturbations (∼150 deg). We chose movies to fully analyze based on the criteria that (i) the fly was in view of all three cameras for sufficiently long to observe pre- and post-perturbation kinematics, (ii) the perturbation primarily affected the fly's pitch orientation, (iii) the fly was not performing any maneuver other than correction, and (iv) we sampled a wide range of perturbation magnitudes for both pitching up and pitching down across our data set. To glean kinematic data from the raw footage, we used a custom-developed image analysis algorithm detailed elsewhere (Beatus et al., 2015; Ristroph et al., 2009). This 3D hull reconstruction algorithm provided a kinematic description of 12 degrees of freedom for the fly (body orientation and center of mass position, as well as three Euler angles for each wing).

## RESULTS

### Body and wing kinematics during pitch correction

Representative kinematics for pitch perturbation events are shown in Figs 1 and 2. Fig. 1B shows definitions for the body Euler angle coordinates – yaw (), roll (ρ_{b}) and pitch (θ_{b}); Fig. 1D–F shows time series of these Euler angles before and after the application of a 5.8 ms magnetic pulse (yellow strip) for 18 perturbation events with at least six different flies. Before the perturbation, flies typically maintained a pitch angle of roughly 50 deg. Perturbations deflected the pitch angle by as much as 40 deg either up or down. While the flies' yaw and roll angles were sometimes altered by the perturbation, pitch was the most consistently and significantly affected degree of freedom immediately following the application of the pulse (at *t*≈6 ms). Highlighted curves in Fig. 1D–F show an event in which the fly was pitched up by 25 deg, attaining its maximal angular deflection at 15 ms after the onset of the perturbation (Movie 1). By 29 ms it had corrected for 90% of the pitch deflection. The maximum pitch velocity from the perturbation was 2400 deg s^{−1}.

The wing kinematics for two representative perturbations, one in which the fly was pitched up by 25 deg (highlighted in Fig. 1D–F) and another in which the fly was pitched down by 23 deg, are shown in Fig. 2. Wing Euler angles – stroke (φ_{w}), pitch (η_{w}) and elevation (θ_{w}) – are defined in Fig. 1C. In general, the wing kinematics we observed following the perturbation were left/right symmetric. For the pitch-up event, ∼10 ms, or 2 wingbeats, after the onset of the perturbation, the minima of the wingstroke angles shifted upward for both the left and right wing (Fig. 2A, orange arrows). During a given wingbeat cycle, the minimum of the stroke angle for each wing corresponds to its front-most position. We refer to the average of the front-most positions for the left and right wings as the front wingstroke angle, . By *t*=15 ms, the fly in Fig. 2A–D had increased its from its pre-perturbation value by ∼25 deg. Physically, this means that the fly was significantly reducing the amplitude of its ventral stroke, i.e. flapping less forward; the duration of this increase in was 3 wingbeats. We did not observe any shifts in the fly's back-most stroke angle, , during the correction maneuvers.

Conversely, for the pitch-down event in Fig. 2E–H, ∼10 ms after the perturbation onset, the pitched-down fly began to decrease its . This corresponds to the fly increasing the amplitude of its ventral stroke, i.e. flapping further forward. Again, there appeared to be little to no change in the back stroke angle.

### Aerodynamic forces and torques

Intuitively, the relationship between front stroke angle and pitching torque can be understood as follows. To within a good approximation, the net aerodynamic force generated by a flapping wing is directed perpendicular to the wing's surface (Dickinson et al., 1999; Sane and Dickinson, 2001), so that portions of the wingstroke during which the wing is in the front half of the stroke plane (φ_{w}≤90 deg) generate pitch-up torques, while portions in the back half (φ_{w}≥90 deg) generate pitch-down torques (Fig. 2D,H, Fig. 3). During non-maneuvering flight, these torques cancel over a wingstroke. By biasing a wingstroke so that it spends a smaller fraction of the stroke period in the forward position, a smaller pitch-up torque is generated during that cycle, such that the net pitch torque will be directed downward. Conversely, by increasing the front stroke angle, and thus increasing the portion of the stroke spent in the front position, flies can generate a net pitch-up torque over the course of a full stroke. This can be observed in Fig. 2 (orange arrows) and Fig. 4, in which active adjustments to front stroke angle resulted in net corrective pitching torques over the course of individual wingbeats.

To quantify the effect of changing the front stroke angle, we calculated the pitching torque generated by the wings during the maneuvers shown in Fig. 2 using the full 3D fly kinematics. To calculate the aerodynamic force generated by the wings, we used a quasi-steady aerodynamic force model that was previously calibrated on a mechanical, scaled-up fly model (Dickinson et al., 1999; Sane and Dickinson, 2001). This model gives the lift (*F*_{L}) and drag (*F*_{D}) forces generated by the wings as:
(1)
(2)

where *C*_{L} and *C*_{D} are the wing's lift and drag coefficients, respectively, and are given as functions of wing angle of attack (α) by Sane and Dickinson (2001); *S* is the wing area; is the non-dimensional radius of the second moment of wing area (given as 0.313 by Cheng et al., 2009); ρ_{0} is the density of air; and *U*_{t} is the wingtip velocity as measured in the lab frame. Drag is directed anti-parallel to the wing tip velocity, and lift is perpendicular to both drag and the wing span vector. The total aerodynamic force is the vector sum of the lift and drag forces. While this is a simple method for calculating forces on flapping wings, we found that it quantitatively captures the relevant force production for both pre- and post-perturbation wing kinematics. We tested the effect of adding rotational forces to the aerodynamic model (Sane and Dickinson, 2001), which should give the next largest contribution to force and account for the force peaks at wingstroke reversal, but found negligible changes to our calculation results.

From these lift and drag forces, we calculated aerodynamic torques exerted by the wings on the body (Fig. 2D,H). The moment arm for the torque is given by the vector from the fly's center of mass to the wing's center of pressure, assumed to be in the chord center, 70% along the length of the wing's span (Fig. 3) (Cheng et al., 2009; Fry et al., 2005). During the active correction of the pitch-up perturbation, the fly's wings generated less upward torque, resulting in a net downward bias for the torque on the body (orange arrows in Fig. 2D). Similarly, the pitched-down fly in Fig. 2H generated net upward pitching torque during active correction (orange arrows in Fig. 2H). The corrective torques for these two events are well correlated with the measured modulations of front stroke angle (Fig. 4A,C), a trend that we observed across all perturbation events (see Results, ‘The corrective effect of stroke angle over a range of perturbations’).

Interestingly, after the flies generated a corrective torque for 2–3 wingbeats, we also observed a few wingbeats in which they generated net torque in the opposite direction (Fig. 2D,H, Fig. 4A,C). As with the initial corrective torque, this subsequent counter-torque arose from modulations of the front stroke angle, evident in Fig. 4A,C. The counter-torque acted to brake the corrective pitching motion, mitigating the overshoot in body pitch angle caused by the initial correction response. This allows for faster correction times, as the initial corrective maneuver can generate larger torques, and thus more quickly return the fly to pitch angles near its original orientation. In movies that allowed us to track the fly for long periods after the perturbation, we observed that the front stroke angle and the net aerodynamic torque often oscillated with decaying amplitude and a period of ∼3–4 wingbeats.

Importantly, passive damping of pitch motion contributed negligibly to the correction maneuvers we observed, in contrast to the case of yaw stability (Cheng et al., 2009; Hedrick et al., 2009; Hesselberg and Lehmann, 2007; Sun, 2014; Warrick et al., 2012). The characteristic time scale at which passive pitch damping becomes significant has been estimated as ∼80 ms (Ristroph et al., 2013), much longer than the entire correction maneuver. Taken together, our results indicate that pitch correction for flies is an active process involving modulation of .

### The importance of stroke angle relative to other degrees of freedom

To assess the effect of front stroke angle modulation on body pitch correction, we calculated the changes to body pitch angle, Δθ_{b}, generated by changes in wing kinematics. We isolated the corrective effect of each of the three wing kinematic variables by first identifying wingstrokes that correspond to both non-maneuvering (no net torque) and corrective flight. We calculated the Δθ_{b} resulting from all eight combinations of corrective (red) and non-corrective (blue) kinematics (color code in Fig. 5). For example, a triplet of squares with color combination red–blue–blue corresponds to a wingbeat with the stroke angle taken from a maneuvering wingbeat, and both wing pitch and elevation angles taken from a non-maneuvering wingbeat; its coordinate along the horizontal axis gives the calculated change in body pitch angle resulting from such a wingbeat. Our calculation assumed rigid wings attached to a stationary body with the geometry as in Fig. 3, and we determined the net change in body pitch angle using numerical integration and assuming an initial condition with zero pitch angular velocity.

We performed this analysis for data from two different perturbation events: the pitch-up and pitch-down events in Figs 2 and 4. The grouping of the points in Fig. 5 indicates that body pitch correction is most closely associated with changes to wingstroke angle. The red–blue–blue point, corresponding to corrective stroke angle but non-maneuvering wing pitch and elevation, achieves at least 60% of the body pitch correction, consistent with Muijres et al. (2014). Moreover, the only combinations that account for more than 30% of the total correction include the stroke angle from a maneuvering wingbeat (points of the form red–*x*–*x*). These results motivate a minimal model for body pitch stabilization that considers only variations in front stroke angle to drive pitch correction.

### The corrective effect of stroke angle over a range of perturbations

To further flesh out the relationship between front stroke angle and corrective torque, we plotted the maximum measured corrective pitch acceleration generated by the fly in each of the 18 maneuvers as a function of the corresponding change in front stroke angle measured at that time, (Fig. 6A). The maximum pitch acceleration was measured at the extremum of θ_{b}, using a quadratic polynomial fit. The plot demonstrates a strong correlation between changes in front stroke angle and corrective acceleration (linear *R*^{2}=0.87). Across our data set of 18 perturbation events, flies increased (flap less forward) to pitch themselves down, and decreased (flap further forward) to pitch themselves up.

The correlation between and corrective pitch acceleration in Fig. 6A is also predicted by a calculation based on the quasi-steady aerodynamic force model in Eqns 1 and 2 (Fig. 6A, gray line). To calculate the aerodynamic pitch torques, we used a simplified wing kinematic model similar to that in Chang and Wang (2014) in which only the front stroke angle is varied (see Appendix). We averaged the computed torques over a wingbeat, and divided by the moment of inertia to obtain pitch acceleration. The calculation relied only on the wing kinematics and fly morphology (Cheng et al., 2009), and had no fitted parameters. The calculation quantitatively reproduced the measured pitch acceleration, further corroborating a model for pitch control that includes only modulation of .

To rule out an alternative corrective mechanism, based on modulation of back stroke angle, we plotted corrective pitch acceleration as a function of (Fig. 6B) and found no discernible correlation between these two variables.

### Correction time scales

We also analyzed the pitch correction time scales across our data set in terms of both the response latency and the overall duration of the correction maneuver. Fig. 6C shows a histogram of latency times for each corrective maneuver, defined as the time between the onset of the perturbation and the first measurable change in the front stroke angle ( deg) for all 18 perturbation events. The mean latency time was 9.9±2.1 ms (mean±s.d., *N=*18), corresponding to ∼2 wingbeats. Fig. 6D plots the total correction time for each maneuver as a function of the maximum body pitch deflection in each perturbation event. We define the correction time as the time between the onset of the perturbation and the fly correcting 90% of the pitch deflection. The mean correction time was 29±8 ms (mean±s.d., *N*=32). Finally, we found that the correction time is weakly correlated with the perturbation amplitude (linear *R*^{2}=0.093), which is consistent with a linear control model.

### Control-theory model

We used a control-theoretic framework to describe the flies' strategy for pitch stabilization. In particular, we modeled actuated changes to the front stroke angle as the output of a PI controller with time delay Δ*T*, for which the input is body pitch velocity (block diagram in Fig. 7A). The response is given by:
(3)

Eqn 3 states that adjustment of the front stroke angle () at a given time *t* is given by a linear combination of the body pitch angle deviation from the pre-perturbation orientation (Δθ_{b}) and body pitch velocity (), both measured at an earlier time *t*−Δ*T*. The parameters *K*_{p} and *K*_{i} are the proportionality constants that determine the relative weights of body pitch angle and pitch velocity. Note that the same controller could be termed a proportional–derivative (PD) controller if the input to the controller was chosen to be the body pitch angle. Because the fly halteres are known to measure body angular velocities (Dickinson, 1999; Nalbach, 1994; Pringle, 1948), we choose the PI nomenclature. We did not consider controller models that depend on angular acceleration (proportional–integral–derivative, PID) based on previous studies that have shown flies' corrective pitch response to be insensitive to angular acceleration (Dickinson, 1999).

Using measured values for, Δθ_{b} and , we fitted for the parameters *K*_{p}, *K*_{i} and Δ*T*. The three parameters were fitted for each movie individually, using one data point per wingstroke. The fit was performed by evaluating Eqn 3 on a dense 3D grid in parameter space and finding the global minimum for the sum of squared residuals between the control model and real data. The results of a controller fit for the pitch-up event in Figs 1, 2 and 4 are shown in Fig. 7B (orange dots, blue line). We found excellent quantitative agreement between the controller fit and our measured data with a root mean square (RMS) error of 1.9 deg, which is of the order of the measurement uncertainty for (Ristroph et al., 2009). The controller model captures not only the sharp rise in in response to the perturbation but also the subsequent decrease in corresponding to the braking counter-torque that slows the fly's downward pitching motion (see Results, ‘Aerodynamic forces and torques’). The fast rise time of the response can be attributed to the proportional term.

We applied the same fitting process to nine movies. We found the values of fitted control parameters (Table 1) to be *K*_{i}=0.3±0.15, *K*_{p}=7±2.1 ms and Δ*T=*6±1.7 ms (means±s.d.), with an average RMS fitting error of 3.0 deg. The mean value of Δ*T* corresponds to ∼1 wingbeat. Fig. 6C shows the region corresponding to mean Δ*T*±1 s.d. (highlighted in orange) compared with measured latency times. The confidence interval (CI) size is large relative to the fitted control parameters (*>*50% in some cases). The large confidence intervals, combined with the accuracy of the fit, indicate that the controller output is insensitive to the choice of model parameters, i.e. that the controller gains and time delay do not require fine tuning.

### Numerical simulation

To corroborate our experimental evidence for the PI controller, we performed a dynamical simulation of a mechanically perturbed fruit fly. The simulation solved the equations of motion for the pitch, longitudinal and vertical degrees of freedom, assuming the fly's geometry, simplified wing kinematics and the quasi-steady aerodynamic force model detailed above (for details, see Appendix). The body pitch angle over time for simulated flies implementing different control strategies is shown in Fig. 8. The four simulated control schemes shown are: (i) proportional–integral (PI, blue), (ii) proportional (P, green), (iii) integral (I, orange) and (iv) no control (red). To determine parameters for the simulated controllers, we fitted the parameters of each model to experimental data. To mimic experimental conditions, we imposed a 5 ms external mechanical torque on the simulated flies (yellow strip), with magnitude comparable to our real system.

Among the four tested models, PI control is the only one that is consistent with the fast, robust pitch control that we observed experimentally (Fig. 8A). We found that flies with no control or I control exhibited large, rapid oscillations of body pitch angle (Fig. 8B), while flies with P control exhibited slightly smaller, long time scale oscillations (Fig. 8A). With I, P and no control, the simulated fly failed to remain aloft and rapidly lost altitude. In contrast, simulated flies implementing PI control corrected their orientation over time scales similar to those in our experimental data, maintained pitch stability over long times, and remained aloft. The general features of each control scheme showed little sensitivity to the values of the control parameters, in agreement with the large confidence intervals of the fitted control parameters obtained above (see ‘Control-theory model’).

### Extreme perturbations

In addition to the 18 perturbation events analyzed in full (Figs 1, 6), we examined two large-amplitude perturbation events of two additional flies. Snapshots and time courses of body pitch angle are shown in Fig. 9 for a pitch-up and a pitch-down event, both with maximum pitch deflection greater than 130 deg (Movie 2). Remarkably, both flies performed successful correction maneuvers, although they were not in-frame long enough to observe them returning to their original orientation. The correction time for both large-amplitude events (*>*50 ms) was longer than the correction times shown in Fig. 6, which can be attributed to the fact that the controlled quantity is biologically constrained: front stroke angle is limited, for instance, by the angle at which the body or the other wing obstructs a wing's forward motion. When we input the body pitch kinematics for the events in Fig. 9 into our PI controller model, the model predicted changes to in excess of 100 deg – a value that is physiologically impossible in the forward direction and not observed in the backward direction. Assuming the flies' corrective response is bounded by deg, a physiologically reasonable estimate, our numerical simulation predicted a response time of ∼70 ms for a perturbation Δθ_{b}=−150 deg, in excellent agreement with the experimental data.

## DISCUSSION

### Front stroke angle as the controlled quantity

We showed that front stroke angle modulation is the primary mechanism for body pitch control in fruit flies, consistent with previous experiments (Dickinson, 1999; Ristroph et al., 2013; Taylor, 2001; Zanker, 1990), and in the same spirit as other proposed mechanisms that include modulation of the mid-stroke angle (Chang and Wang, 2014). Calculated aerodynamic forces predicted that changes to the back stroke angle could produce corrective pitching torques in the same way that changes to front stroke angle do; the fact that we did not observe this in the data hints that morphological constraints favor modulation of the front stroke angle. Our computational results (Figs 6, 8) showed that a minimal model, which only incorporates changes to front stroke angle, and uses control parameters extracted from fits to our measurements, is the simplest linear, continuous model capable of stabilizing the body pitch angle on time scales similar to those observed in the experiments.

Kinematic variables other than may also contribute to pitch correction. Previous studies have associated changes to stroke plane elevation (Zanker, 1988), wingbeat frequency (Dickinson, 1999) and body posture (Taylor, 2001) with pitch correction. In particular, we observed transient alterations in both wing pitch and elevation angle during corrective maneuvers. Fig. 5 suggests that, when combined with modulation of front stroke angle, changes to wing pitch and elevation angle can account for up to 40% of body pitch correction, consistent with Muijres et al. (2014). The detailed role of these kinematic variables in pitch control and whether they are actively or passively actuated remains unknown (Bergou et al., 2010).

### Long-term stabilization

Importantly, the PI controller model presented here, which assumes a ‘dead reckoning’ method of angle estimation based on angular velocity input, cannot account for pitch control on long time scales because of integration and sensor errors affecting the estimation of the absolute pitch angle. Long-term pitch control in this framework requires direct measurement of the pitch angle, as could be achieved by the visual system at longer time scales (see Dyhr et al., 2013 for a discussion of time delays in the visual response system). An interplay between the haltere and visual systems (as in Huston and Krapp, 2009; Sherman and Dickinson, 2004) would thus be necessary for comprehensive pitch stability. An intriguing alternative to the dead reckoning assumption is that the fly could implement a model-based estimator to make measurements of its absolute pitch angle using only information about its pitch velocity (N. Cowan, personal communication). Such a model-based estimator is indeed possible: the linearized, flying insect pitch dynamics in Ristroph et al. (2013) give rise to an observability matrix that is full rank, indicating that body pitch angle should be observable with only angular velocity feedback. Whether this method for control is actually used by flies would be the subject of further research. The PI model presented here, however, can completely and accurately account for the fly's fast reflex response, which stabilizes it against rapid pitch perturbations.

### Discrete versus continuous control models

The periodic motion of wing flapping introduces inherent discreteness to insect flight. For processes occurring on time scales comparable to a wingstroke period – like the perturbations and maneuvers we reported here – we expect discrete effects to be more pronounced. In particular, modulations of front stroke angle can, by definition, occur only once per wingbeat. Because perturbations can be induced at any time during the wingbeat, but the actuated kinematics are discretely constrained, latency times for correction depend on the phase of the perturbation relative to the wingstroke. Latency times will be bounded from below by the flies' neural response time, but could potentially be as much as 1 wingbeat longer as a result of the phase of the perturbation within the wingbeat.

Measured latency times can also depend on discrete sensing. The temporal sampling resolution with which flies can measure mechanical perturbations is likely determined by motion of their halteres, the rate-gyro sensory organs used in the fast perturbation response (Dickinson, 1999). Dipteran halteres beat at the wing frequency and use Coriolis forces to measure body angular velocities (Nalbach, 1994; Pringle, 1949). The largest sensitivity to mechanical perturbations is likely to occur at times during the fly's mid-stroke, when Coriolis forces on the halteres are the largest (Nalbach, 1994). Sensing at discrete times introduces a second relevant phase for correction latency time: the phase of the perturbation relative to sensing. Similar to discrete actuation, discrete sensing would lead to latency times longer than the neural response time. Moreover, even during the fly's mid-stroke, its halteres only have finite sensitivity. It is likely that there exists some threshold for angular velocities that are large enough to elicit a control response (Fox and Daniel, 2008).

The continuous PI controller model does not account for the aforementioned effects of discrete actuation, discrete sampling or sensing threshold. Hence, the measured latency time should constitute an upper bound for the delay time that we obtain from the controller model. Indeed, the time delay from our controller model (Δ*T=*6±1.7 ms) is roughly 1 wingbeat shorter than the measured latency time (9.9±2.1 ms).

Despite the inherent discreteness of the fly control systems, the continuous PI controller model quantitatively captures the behavior of flies in response to pitch perturbations (Fig. 7). This quantitative agreement leads to an interesting open question: under what conditions does it become necessary to use a discrete controller model to describe flight stabilization? Previous studies on the legged locomotion of cockroaches have shown quantitative consistency between discrete and continuous control models in the context of wall following (Lee et al., 2008). To address this question in the context of flight stabilization would require precise perturbation timing, in order to probe the short time scales at which discretization becomes relevant. Such an analysis could provide significant insight into the timing and thresholding of fruit fly reflexes.

### Physiological basis for pitch control

Both the mechanism and timing of the pitch correction indicate a likely candidate muscle for control actuation: the first basalare muscle (b1), as suggested by previous studies (Chang and Wang, 2014; Fayyazuddin and Dickinson, 1999). Among dipteran flight control muscles, b1 is unique in that it is active during every wingstroke (Heide, 1983; Heide and Götz, 1996; Miyan and Ewing, 1985), which would allow for the wingbeat time scale pitch control that we observed. Moreover, b1 activity is strongly correlated with modulations of ventral stroke amplitude, i.e. changes in (Dickinson and Tu, 1997; Walker et al., 2014). In blowflies, b1 activity is also correlated with changes in elevation angle during the ventral stroke (Balint and Dickinson, 2004), which could explain the slight shifts in elevation angle that we observed during correction (Fig. 2B,F). Our results indirectly support the hypothesis that the b1 muscle is responsible for pitch control through the regulation of . Testing flies with disabled or altered b1 muscles could provide an avenue for confirming the role of b1 in the body pitch control process.

### Linear control of body orientation

In addition to the results on fruit flies reported here, PI control has also been identified in pitch control for hawkmoths (Cheng et al., 2011; Windsor et al., 2014). The anatomical similarities found across species suggest that pitch instability is an obstacle faced by many flapping insects (Ristroph et al., 2013; Sun, 2014; Sun et al., 2007); a natural question raised by these collective findings is whether PI control is a generic feature of pitch stabilization in flying insects. Beyond flying insects, what we refer to as PI control has also been observed in fast obstacle avoidance in pigeons (Lin et al., 2014) and antenna-based wall following in cockroaches (Cowan et al., 2006; Lee et al., 2008). Future research on the ubiquity of PI control could have fascinating implications for the evolution of flight stabilization and sensorimotor control mechanisms.

In the larger context of comprehensive flight stabilization, our results, together with previous studies on yaw (Ristroph et al., 2010) and roll (Beatus et al., 2015) control in fruit flies, show that the strategies flies use to control each of their body Euler angles can be modeled as PI controllers. However, the overarching structure in which these three individual controllers are embedded is still unknown. Given the non-commutativity of rotations in 3D, the relationship between controllers that measure different angular coordinates is likely to be non-trivial. For example, a previous study has shown two cases where the fly's response to certain perturbation could not be explained by a superposition of linear controllers (Beatus et al., 2015). First, in response to perturbations that simultaneously affected both the roll and yaw angles, preferential correction for roll over yaw was observed, hinting that the control of these two degrees of freedom is coupled. Second, it was shown that the response of flies to extreme perturbations consisting of multiple rotations along roll cannot be explained by a linear control model. Taken together, these results hint at a complex and intriguing control architecture used by flies to stabilize their orientation that may depend on the amplitude and timing of the perturbations along each axis.

Finally, an understanding of the relationship between control of different Euler angles could have profound implications for how the fly encodes information about its body orientation. In the case of vision, organism-specific demands have spurred the development of novel, specialized neural structures in both mammals (Hafting et al., 2005; Yartsev et al., 2011) and insects (Ofstad et al., 2011; Seelig and Jayaraman, 2013). Pioneering work on information processing from halteres has suggested similar morphology/function relationships for the gyroscopic rate sensing in insects (Fox et al., 2010). Connecting such analyses with the resultant control structure observed in free-flight behavioral experiments could provide a window into the most basic ways in which flies sense and interpret the world.

## Acknowledgements

We thank Grace (Li) Chi and Andy Clark for providing flies; Ty Hedrick for advice on camera calibration; and Andy Ruina and the Cohen group members for useful conversations. We especially thank Professor Noah Cowan for his insights and the important role he has played in this and previous manuscripts helping us develop and apply control theory models to our systems.

## Appendix

### Simplified wing kinematics

For both the calculation of aerodynamic torques in Fig. 6A and the dynamical simulation in Fig. 8, we use an analytic form for simplified wing kinematics taken from Chang and Wang (2014). These kinematics closely resemble the motion of real fly wings, but are simple enough to write down concisely: (A 1) (A 2) (A 3)

The wing Euler angles are defined in Fig. 1C. The terms in Eqns A1–3 are defined as follows: φ_{0}, η_{0} and θ_{0} are angle offsets; φ_{m}, η_{m} and θ_{m} are amplitudes; *K* and *C* are waveform parameters – *K* tunes the stroke angle from pure sine wave to triangle wave, while *C* tunes the wing pitch angle from sinusoid to square wave; ω is the wingbeat frequency; and δ_{η} and δ_{θ} are phase offsets.

For the parameter values used in the calculation in Fig. 6A, see Table 2. Note that in the main text we refer to front and back stroke angles ( and ), which are related to φ_{m} and φ_{0} by the linear relationships:
(A 4)
(A 5)

### Numerical simulation

Using the simplified wing kinematics above and the quasi-steady aerodynamic model detailed in the main text (Eqns 1 and 2), our numerical simulation solves the Newton–Euler equations for vertical (*z*) and forward (*x*) motion of the body center of mass, as well as pitch rotational (θ_{b}) motion:
(A 6)
(A 7)
(A 8)

where *M* is the fly body mass, *I*_{pitch} is its pitch moment of inertia, ** g** is the gravitational acceleration,

**F**

_{w}is the wing aerodynamic force and

**r**

_{CoP}is the vector from the fly's center of mass (CoM) to the wing center of pressure (CoP).

**r**_{CoP} incorporates both the motion of the wing relative to the hinge and a fixed distance between the body CoM and the wing hinge, and can thus be written in the fly body frame as (apostrophe denotes body frame). The vector from the body CoM to the wing hinge, **r**_{hinge}, is a fixed length estimated from our 3D reconstructions, and is given in the lab frame in Table 2. *R* is the wing span, and the unit vector corresponding to the direction of the wing span relative to the hinge can be written as:
(A 9)

Note that these are expressions for **r**_{CoP}, the moment arm vector in a frame with axes fixed to the body; for our calculations, we used a fixed-orientation frame, and thus calculated **r**_{CoP}=**Q****r**_{CoP}, where **Q** is a rotation matrix. We determine **F**_{w} using the quasi-steady model in Eqns 1 and 2 of the main text, with the wingtip velocity written in the lab frame as:
(A 10)

We calculated the aerodynamic force **F**_{w}, including the counter-flapping torque (CFT) terms described in Hedrick et al. (2009). Note, however, that our equations of motion do not include drag on the body, which we assumed to be negligible. We chose to solve the non-linear equations of motion, rather than the linearized dynamics in Ristroph et al. (2013).

Control is implemented by adjusting the front stroke angle of the prescribed wing kinematics according to Eqn 3. The *K*_{i} and *K*_{p} parameters were determined by the fit to the experimental data. The inputs for the controller in each wingbeat – the body pitch angle and pitch velocity – were taken as their mean values during the previous wingbeat. This scheme represents a time delay of 1 wingbeat while avoiding the effects of the inherent small-scale pitch oscillations. Before we apply the perturbation, we let the simulated fly stabilize to a steady-state body pitch angle of 45 deg. The perturbation is then applied, with magnitude roughly corresponding to the accelerations observed in the experiments. In simulation runs that tested controller models other than PI, we let the system stabilize at θ_{b}=45 deg using a PI controller and only then applied the perturbation and simultaneously changed the controller type.

### Analysis of data from previous studies

Following Beatus et al. (2015), we used the PI control model to predict the pitch response of tethered fruit flies previously published by Dickinson (1999). In Dickinson (1999), flies were tethered to a gimbal apparatus that oscillated about different rotational axes. The left wingstroke amplitude of the flies, ɸ_{left}, was measured using photodetectors that recorded the wings' shadows. As noted by Dickinson (1999), flies do not adjust their back stroke angle during pitch correction, so stroke amplitude is a good proxy for , the quantity we measured in free-flight experiments.

In one of the measurements reported in Dickinson (1999), pitch perturbations were imposed so that: θ_{b}(*t*)=*A*sin(ω*t*) and , with *A=*25 deg, period *T=*0.63 s and maximum pitch velocity 250 deg s^{−1}. The left wingstroke amplitude was plotted against both pitch angle and pitch velocity (fig. 3A,C in Dickinson, 1999). Using standard image-processing techniques, we extract the data from these plots.

The pitch oscillations in the tethered experiments have period of 630 ms, which is much longer than the observed pitch correction latency times from our experiments (≈10 ms), so we consider the controller time delay negligible. We then write the form for our controller, now in terms of left wingstroke amplitude, as: (A 11)

where the left wingstroke amplitude, ɸ_{left}(*t*), and mean stroke amplitude, ɸ_{mean}, are related to by the linear relationship: . Considering only the left wing does not reduce the generality of this analysis, as pitch correction is left/right symmetric. We manually fitted for the control parameters from the data. The fitted parameters obtained are *K*_{i}=0.3 and *K*_{p}=8 ms, comparable to the parameters from the main text.

The predictions of the PI controller fit are shown in Figs S1–3. The output of the PI controller is plotted as a function of both pitch angle and pitch velocity in Fig. S1, yielding an ellipse in both cases (*R*^{2}=0.761 and 0.842, respectively). The linear model for ɸ_{left} as a function of gives *R*^{2}=0.556. We also plotted the PI controller prediction in the 3D space whose axes are (θ_{b}, , ɸ_{left}) (Fig. S2). The PI controller predicts an inclined ellipse in this space, the projections of which onto the horizontal axes yield the plots in Fig. S1. The inclination of the ellipse shows that the corrective response depends on both pitch angle and pitch velocity, i.e. the measured data in Dickinson (1999) are consistent with a PI controller.

Additionally, we show the predicted output of the PI controller plotted as a function of pitch acceleration in Fig. S3. Consistent with Dickinson (1999), Fig. S3 shows that the fly's corrective response can be quantitatively captured without including information about the pitch acceleration. Fig. S3 suggests that pitch acceleration is unimportant in determining the flies' corrective wing kinematics, and thus excludes a PID controller model.

## FOOTNOTES

**Competing interests**The authors declare no competing or financial interests.

**Author contributions**S.C.W. and L.C. conducted experiments and collected data. T.B. and S.C.W. developed code for data analysis and simulation. S.C.W. analyzed the data. S.C.W., T.B. and I.C. wrote the paper.

**Funding**This work was supported in part by a National Science Foundation (NSF) DMR award (no. 1056662) and in part by an Army Research Office (ARO) award (no. 61651-EG). S.C.W. was supported by a National Defense Science and Engineering Graduate (NDSEG) Fellowship. T.B. was supported by the Cross Disciplinary Post-Doctoral Fellowship of the Human Frontier Science Program.

**Supplementary information**Supplementary information available online at http://jeb.biologists.org/lookup/suppl/doi:10.1242/jeb.122622/-/DC1

- © 2015. Published by The Company of Biologists Ltd