## ABSTRACT

Flies fly at a broad range of speeds and produce sophisticated aerial maneuvers with precisely controlled wing movements. Remarkably, only subtle changes in wing motion are used by flies to produce aerial maneuvers, resulting in little directional tilt of the aerodynamic force vector relative to the body. Therefore, it is often considered that flies fly according to a helicopter model and control speed mainly via force vectoring by body pitch change. Here, we examined the speed control of bluebottle flies using a magnetically levitated (MAGLEV) flight mill, as they fly at different body pitch angles and with different augmented aerodynamic damping. We identified wing kinematic contributors to the changes of estimated aerodynamic force through testing and comparing two force-vectoring models – a constant force-vectoring model and a variable force-vectoring model – while using Akaike's information criterion for selection of the best-approximating model. The results show that the best-approximating variable force-vectoring model, which includes the effects of wing kinematic changes, yields a considerably more accurate prediction of flight speed, particularly in higher velocity range, as compared with those of the constant force-vectoring model. Examination of the variable force-vectoring model reveals that, in the flight mill tethered flight, flies use a collection of wing kinematic variables to control primarily the force magnitude, while the force direction is also modulated, albeit to a smaller extent compared with those due to the changes in body pitch. The roles of these wing kinematic variables are analogous to those of throttle, and collective and cyclic pitch of helicopters.

## INTRODUCTION

Flies are eminent miniature flyers that exercise stable and agile flight over a large flight envelop (Beatus et al., 2015; Fry et al., 2003; Muijres et al., 2014). This aerial success hinges partly on flies' ability to precisely control subtle wing movement through regulating the firing rate and timing of steering muscles (Dickinson and Tu, 1997; Lindsay et al., 2017), despite their wings being a difficult locomotor apparatus to control neuromuscularly (Balint, 2004; Deora et al., 2015) or to emulate in engineering designs (Keennon et al., 2012; Ma et al., 2013; Roll et al., 2013). Previous research has shown that although flies are able to produce large modulation of wing motion through the clutch and gearing mechanisms at the wing hinge (Deora et al., 2015), even subtle changes of wing motion are sufficient to produce a large maneuvering moment for a fly to execute rapid maneuvers, for example during saccades (Fry et al., 2003), evasive maneuvers (Muijres et al., 2014) and recovery from aerial stumbles (Ristroph et al., 2010). However, such subtle changes only result in small directional changes of the aerodynamic force vector relative to the body; therefore, flies maneuver mostly according to the helicopter model (Medici and Fry, 2012; Muijres et al., 2014), i.e. they change direction mainly by vectoring their body and the averaged aerodynamic force vector together towards the intended flight direction, while maintaining their relative angle (i.e. the angle between the body and the aerodynamic force vector) approximately constant.

The helicopter model may also apply to forward flight, as the flight speed of flies and other insects is well known to tightly correlate with body pitch angle (David, 1978; Dudley and Ellington, 1990; Meng and Sun, 2016; Willmott and Ellington, 1997), suggesting that the tilt of the aerodynamic force vector might also be small during forward flight. However, key questions remain: what wing kinematic variables do flies use to control flight speed?; how do these variables vary with body pitch and thrust force?; and to what degree do flies change the magnitude and direction of aerodynamic forces while flying at different speeds? While the answers to these questions remain elusive, they are of critical importance for insect flight research and also for inspiring novel engineered flight, especially considering that flies and other insects fly at a broad range of speeds and produce large linear acceleration while foraging, chasing mates and escaping from predators (Collett and Land, 1975; Dudley, 2002). For example, locusts (*Nomadacris septemfasciata*) reach a speed of 13 m s^{−1} in a wind tunnel (Waloff, 1972), drone flies (*Eristalis tenax*) reach 8.5 m s^{−1} (Meng and Sun, 2016), hawkmoths (*Manduca sexta*) reach 5 m s^{−1} (Willmott and Ellington, 1997) and bumblebees (*Bombus terrestris*) reach 4.5 m s^{−1} (Dudley and Ellington, 1990); dragonflies (*Plathemis lydia*) accelerate near 2 ** g** in prey interception flights (Mischiati et al., 2014) and bluebottle flies (

*Calliphora vicina*) demonstrate 3

**acceleration in a free-flight chamber (Bomphrey et al., 2009).**

*g*Forward flight of insects is commonly studied in laboratory settings using wind tunnels and flight mills. Wind tunnel studies range from the observation of body and wing kinematics in free flight (Azuma and Watanabe, 1988; David, 1978; Dudley and Ellington, 1990; Meng and Sun, 2016; Willmott and Ellington, 1997) and tethered flight settings (Vogel, 1966), identifying visual control principles (Baird et al., 2005; Fry et al., 2009; Medici and Fry, 2012; Srinivasan et al., 1996), and kinematic data-driven modeling of flight control and stabilization (Fuller et al., 2014). Flight mills – devices that approximate continuous forward flight in a confined space through restricting an insect to a circular flight path around a pivot joint – are commonly used to determine the traveling distance of insects and their dispersal potentials (Attisano et al., 2015; Ranius, 2006; Ribak et al., 2017). Although flight mills are rarely used to study other aspects of insect forward flight, they have the potential to provide more naturalistic visual and proprioceptive sensory feedback than wind tunnel experiments. This is important because vision plays a key role in regulating forward flight speed. Previous studies have found that many insects (e.g. honeybees and flies) can robustly extract their ground speed (or retinal slip velocity) from visual patterns of varying spatial and temporal frequencies (David, 1982; Fry et al., 2009). Flies have also been shown in the wind tunnel experiments to sometimes maintain a preferred ground speed invariant to substantial changes in airspeed (David, 1982) (putatively detected by air flow sensors, e.g. antenna: Fuller et al., 2014). Therefore, flies, and possibly other insects also, are able to fly at their preferred speed independent of the aerodynamic power requirement if it is within their locomotor capacity. However, it is unknown whether such behavior can be reproduced in the flight mill experiments, how insects control their flight speed in the flight mill, or what happens when the limit of their locomotor capacity is reached.

*A*- area
- AIC
- Akaike’s information criterion
- AICc
- small sample size-corrected AIC
- AoA
- angle of attack
- BIC
- Bayesian information criterion
- average damping coefficient
*C*_{D}- drag coefficient
- D0–D2
- no damper, medium damper and large damper
*F*- resultant force
*F*_{lift}- lift force
*F*_{thrust}- thrust force
*I*- moment of inertia
*k*- constant of integration
*K*- number of parameters in a candidate model
*K*,_{jF}*K*_{j}_{χ}- standardized regression coefficients of wing kinematic variables contributing to force magnitude and force direction
*l*- radius of flight mill shaft
*m*- total number of trials
*M*- total number of objects
*n*- wingbeat frequency
*r*- radial distance
*R*- wing length
*R*^{2}- coefficient of determination
- RMSE
- root mean squared error
- RSS
- residual sum of squares
*t*- time
- dimensionless time
*T*_{d}/*T*_{u}- ratio of downstroke to upstroke duration
*v*- body forward velocity
*w*- Akaike weight
*w*_{+}- summation of Akaike weights
*X*_{b}- body roll axis
*X*,_{jF}*X*_{j}_{χ}- wing kinematic variables contributing to force magnitude and force direction
*X*_{w}- wing normal axis
- (
*X*,*Y*,*Z*) - global coordinate frame
- (
*X*_{b},*Y*_{b},*Z*_{b}) - body coordinate frame
- (
*X*_{s},*Y*_{s},*Z*_{s}) - stroke plane coordinate frame
- (
*X*_{w},*Y*_{w},*Z*_{w}) - wing coordinate frame
*Y*_{b}- body pitch axis
*Y*_{w}- wing spanwise axis
*Z*_{b}- body yaw axis
*Z*_{w}- wing chordwise axis
- β
- stroke plane angle relative to horizontal
- β+χ
- stroke plane angle relative to body long axis
- γ
- angle-pin angle
- Δ
_{i} - AIC difference of
*i*th model - Θ
- deviation amplitude
- mean deviation angle
- θ
_{0}, θ_{si}, θ_{ci} - constant, sine, and cosine Fourier coefficients of stroke deviation
- ρ
- air density
- τ
_{D} - torque due to aerodynamic drag
- τ
_{T} - torque due to thrust
- Φ
- stroke amplitude
- mean stroke angle
- φ
_{0}, φ_{si}, φ_{ci} - constant, sine and cosine Fourier coefficients of stroke position
- χ
- body pitch angle
- χ
_{F} - aerodynamic force angle from body longitudinal axis
- χ
_{0} - angle between total force vector and body pitch axis
- Ψ
- rotation amplitude
- mean rotation angle
- ψ
_{0}, ψ_{si}, ψ_{ci} - constant, sine and cosine Fourier coefficients of wing rotation
- ω
- flight mill angular velocity

In this study, we examined the speed control of bluebottle flies (*Calliphora vomitoria*) in forward flight using a novel magnetically levitated (MAGLEV) flight mill. The MAGLEV flight mill eliminates the mechanical friction of pivots and permits systematic manipulation of a fly's body pitch angle and aerodynamic damping; thereby, it enabled us to study the details of speed control and force vectoring and the corresponding wing kinematic control in forward flight. In particular, we tested and compared a constant force-vectoring model (without considering the wing kinematics) and a variable force-vectoring model (considering the wing kinematics) to determine the wing kinematic contributors to the changes in magnitude and direction of aerodynamic forces.

## MATERIALS AND METHODS

### MAGLEV flight mill apparatus

The MAGLEV flight mill apparatus was composed of four main components (Fig. 1A): (1) three magnetically levitated permanent magnets as the pivot joint; (2) a horizontally rotating shaft with attached fly and damper; (3) inner and outer enclosing walls with grating patterns; and (4) three high-speed video cameras (Fastcam Mini UX100, Photron, Tokyo, Japan). Magnetic levitation was achieved through two electromagnet actuators that stabilized the vertical position of the permanent magnets (i.e. the pivot joint) and rotating shaft using positional feedback provided by two linear Hall-effect sensors (A1321, Allegro microsystem, LLC, Worcester, MA, USA). The first Hall-effect sensor was placed slightly above the permanent magnets to measure the total magnetic field of the permanent magnets and electromagnets combined. The second Hall-effect sensor was attached to the rim of the top electromagnet to separate the noise (magnetic field of electromagnets) from the first Hall-effect sensor. The strength of the magnetic field was then transformed to distance as a proximity signal. A proportional integral derivative (PID) controller computed the current compensations for the two electromagnets to keep the pivot pin/rotating shaft vertically stable. All sensor readings and computations were processed with a microcontroller (Uno, Arduino, Italy).

The shaft was made of a 2×254 mm (diameter×length) carbon fiber rod, which was sandwiched between the permanent magnets. On one end of the shaft, a magnetic metal angle-pin connected a bluebottle fly to the shaft through two micro-permanent magnets (1.58×3.18 mm, diameter×length; Fig. 1A,B) glued to the fly's dorsal side. Caution was taken in the gluing process to minimize interference with the fly's thoracic movements and to maintain a constant angle between the angle-pin and the fly's body. Note that exact body pitch angle (χ) was calculated from DLTdv6 (Hedrick, 2008) instead of angle-pin angle (γ) as subtle differences existed among flies and each angle-pin placement (Fig. 1B). With the permanent magnet attached on the dorsal side, different magnetic angle-pins can be easily switched to provide different prescribed pitch angles. On the other end of the shaft, a damper was attached to create additional drag that the fly needed to overcome in steady forward flight. In other words, the damper introduced augmented aerodynamic drag to manipulate the amount of thrust that the flies needed to produce while flying at a fixed speed, therefore to better reveal their thrust production mechanism and to increase the richness of the data for identifying the force-vectoring models. Two dampers of different size (D1: 12.4×12.7×2.5 mm and D2: 15.2×16.2×2.5 mm, width×length×thickness), together with the no-damper case (D0, where the aerodynamic drag only came from the shaft and the insect body), were used in the experiments, leading to a total of three augmented aerodynamic damping conditions.

Flies mainly rely on visual feedback to regulate their flight speed. To provide a consistent visual environment and to enhance the visual cues they receive, an inner cylinder wall (diameter 203.2 mm) and an outer cylinder wall (diameter 304.8 mm) with identical square-wave grating patterns (50.8 mm interval) were used to enclose the flight mill. As a result, the flies flew in the circular corridor (width 50.8 mm) between the two walls (Fig. 1A).

To record the body and wing movements of the flies, three synchronized high-speed cameras were placed on the top, bottom and sides of an enclosure region spanning approximately 50×50×50 mm of the circular flight corridor (Fig. 1A). A circular hole was cut on the outer wall for the side camera to enable a view through the corridor. We illuminated the enclosure region with three 100 W LEDs (MonoBright LED Bi-color 750, Genaray, Brooklyn, NY, USA). The video resolution was set to 1280×1024 pixels with a frame rate of 4000 frames s^{−1} and a shutter rate of 1/8000 s. Cameras were calibrated using direct linear transformation for 3D body and wing kinematics extraction (Hedrick, 2008).

### Animal preparation

We used 4–7 day old bluebottle flies, *Calliphora vomitoria* (Linnaeus 1758) (*N*=5, 42.0±8.9 mg; Table S1), hatched from pupae purchased commercially (Mantisplace, Olmsted Falls, OH, USA) and cultured in the laboratory. For each experiment, we first cold anesthetized the flies in a refrigerator for 10 min (Duistermars and Frye, 2008) and then transferred them to a tethering stage on an oval notch plate, dorsal side up. We used UV cure glue (4305, Loctite Corp.) to attach the micro-permanent magnet on the dorsal side of the thorax. Flies were allowed to recover from anesthesia for 1 h. We then attached the flies to the rotating shaft of the flight mill and started the experiments.

### Experimental procedure

Each fly was attached to one end of the rotating shaft with angle-pins held at 0, 22.5 and 45 deg (Fig. 1B). The actual body pitch angles measured from the experiments were 5.5±3.9, 25.2±3.3 and 41.2±6.2 deg, respectively. The slight differences between the angle-pin angle and the body pitch angle were mainly due to the slight misalignment of the angle-pins to the normal of the flies' thorax. For each angle-pin, three aerodynamic damping conditions (described above) were tested. In total, nine different experimental conditions (three angle-pin angles and three damping conditions) were tested for each fly.

To initiate flight, a gentle puff of wind was applied to the fly. The fly reached steady flight as the flight mill shaft rotated at approximately constant angular velocity. At steady flight, it can be assumed that wing thrust was balanced by the total aerodynamic drag acting on the damper, shaft and fly's body. We tested the flight performance of the flies on the flight mill and only those that could complete at least five laps of flight with the 45 deg angle-pin and the largest damper (D2) were used for the experiments. We recorded the flight kinematics using high-speed cameras in the period of time following at least five initial laps of flight and before the final five laps of flight (sample recordings are available in Movies 1–3), during which the flies were flying at approximately constant flight speed with negligible acceleration or deceleration. For each individual at each experimental condition, at least four repeated trials were performed; excluding the trials with body or wing occlusion, this yielded a total of 14–19 trials for each experimental condition among all individuals (*N*=5; a total of 154 trials for all nine experimental conditions). For each trial, at least four wingbeat cycles were recorded and used in the wing kinematics analysis. After completing the trials within one condition, the flies were removed from the flight mill, allowed to rest for at least 10 min, and fed with sugar water before being used for the next condition.

### Damping calibration

The damping coefficients of the combined damper, rotating shaft and fly's body for nine different conditions were calibrated using free responses of the rotating shaft. For each calibration, a dead fly with its wings removed was attached to the shaft using one of the three angle-pins. To initiate the free response of the rotating shaft, a wind gust was applied to the damper. We recorded the time instants of the rotating shaft passing the laser beams between laser–photodiode pairs spaced with known distance, and then calculated the angular velocity (ω) of the shaft accordingly. Using blade-element analysis (Leishman, 2006) to model the aerodynamic drag, which is assumed to be quadratic, it can be shown that the equation of motion of the shaft is:
(1)where *I* is the moment of inertia of the shaft, damper and the fly's body, is the calibrated average damping coefficient, ρ is air density, *M* is the total number of objects that contribute to drag force (e.g. the shaft, dampers and insect body) and *i* is the index of an object. For the *i*th object, *A _{i}*(

*r*) is the cross-sectional area of a blade element at radial distance

*r*from the shaft center of rotation and

*C*

_{D,i}(

*r*) is the corresponding drag coefficient. Integrating Eqn 1 yields the theoretical speed profile of the shaft: (2)where

*k*is the constant of integration. We then performed least-squares curve fitting to obtain for each experimental condition; the mean, standard deviation and coefficient of determination

*R*

^{2}are reported in Table 1 (see Fig. S1 for sample calibration results). Note that the fitting error was relatively small within the range of the flies' forward velocity in the flight mill (≳0.09 m s

^{−1}, Fig. 2A, or shaft angular velocity ≳0.7 rad s

^{−1}), although the fitting error became larger at lower velocities outside this range (Fig. S1).

### Kinematics extraction

We used DLTdv6 (Hedrick, 2008) to digitize anatomical landmarks on the body and wings of the bluebottle flies (Fig. 1C), which were then used to calculate the body velocity, pitch angle and wing angles (Fig. 1C,D). We defined the body roll axis (*X*_{b}) as a unit vector from the thorax–abdomen junction to the head–thorax junction, body pitch axis (*Y*_{b}) from the left wing base to the right wing base, and body yaw axis (*Z*_{b}) using the cross-product of *X*_{b} and *Y*_{b}. We defined wing spanwise axis (*Y*_{w}) as a vector from wing base to wing tip. The cross-product of *Y*_{w} with a vector from wing base to the trailing edge location of vein CuA_{1} determined the wing normal axis *X*_{w}. Next, the cross-product of *X*_{w} and *Y*_{w} determined wing chordwise axis *Z*_{w}. Body and wing rotation matrices were then calculated based on the corresponding body and wing principal axes (Murray et al., 1994), respectively. A stroke plane coordinate frame (*X*_{s},*Y*_{s},*Z*_{s}) was defined by rotating the body frame about *Y*_{b} where *X*_{s} intersected with the maximum and minimum sweep positions formed by the wing base–wing tip vectors.

Wing Euler angles (stroke position φ, stroke deviation θ and wing rotation ψ; Fig. 1D) were calculated from the wing rotation matrices (Murray et al., 1994) (from the stroke plane frame to the wing frame) and body pitch angles were calculated from the body rotation matrices (Fig. 1D). Body translational velocity about each principal axis was calculated by taking the derivatives of the head–thorax junction positional vector. Time series of body and wing kinematics were calculated for four complete wingbeat cycles. For each trial, we then calculated time-averaged wing kinematics from the four wingbeat cycles, while also averaging the left and right wing kinematics (mirrored with respect to the *X*_{b}–*Z*_{b} plane). Euler angles of each trial were parameterized using a fifth-order Fourier series prior to kinematics analysis:
(3)
(4)
(5)where is the dimensionless time of a wingbeat cycle (ranging from 0 to 1); φ_{0}, θ_{0} and ψ_{0} are constant terms; φ_{s}* _{i}*, φ

_{c}

*, θ*

_{i}_{s}

*, θ*

_{i}_{c}

*, ψ*

_{i}_{s}

*and ψ*

_{i}_{c}

*are Fourier series coefficients; and*

_{i}*i*is the order of the Fourier series. Fourier series coefficients for each experimental condition were obtained through parametric fitting using the measured wing kinematics.

As flies change their continuous wingbeat trajectories to modulate aerodynamic forces and moments, the key changes can be captured by a finite number of wing kinematic variables that represent certain cycle-averaged features (Faruque and Sean Humbert, 2010; Sun, 2014; Taylor, 2001). Here, we selected nine distinct variables that were potentially involved in speed control and tested their contribution in the force-vectoring models (see ‘Force-vectoring model’, below). The nine wing kinematic variables were: (1) mean wingbeat frequency (*n*), (2) ratio of downstroke to upstroke duration (*T*_{d}/*T*_{u}), (3) stroke amplitude (Φ), (4) rotation amplitude (Ψ), (5) deviation amplitude (Θ), (6) mean stroke angle , (7) mean rotation angle , (8) mean deviation angle and (9) stroke plane angle (β+χ).

### Force-vectoring model

Here, we developed two force-vectoring models for the speed control of flies flying steadily in the MAGLEV flight mill: (1) a variable force-vectoring model and (2) a constant force-vectoring model. The variable force-vectoring model assumes that both body pitch and changes of wing kinematics contribute to force modulation, while the constant force-vectoring model assumes the contribution from body pitch only. By testing and comparing these two models, we can identify the roles of wing kinematics in modulating the force magnitude and direction, which will also allow us to examine the helicopter model, via identifying the degree of force magnitude and direction change due to the changes of wing kinematics. The following describes the details of model derivation.

During steady flight, the torque acting on the shaft of the flight mill was zero, which meant that the torque due to the thrust created by the flapping wings (τ_{T}) was equal to that due to the aerodynamic drag of the shaft, damper and insect body combined (τ_{D}), the latter being proportional to the linear speed of the flies or the angular velocity of the shaft. Assuming the flapping wings create a cycle-averaged aerodynamic force with a magnitude of *F* and an angle χ* _{F}* from the body longitudinal axis (Fig. 1D), it can be shown that:
(6)where

*l*is the radius of the pivot of the flight mill shaft and

*v*is the linear velocity of the fly. We can rearrange Eqn 6 into: (7)where

*v*can be predicted once

*F*and χ

*are known. It was assumed that both*

_{F}*F*and χ

*also depended on a collection of wing kinematic variables, e.g. stroke amplitude, frequency and mean rotation angle. Therefore, it was assumed that*

_{F}*F*=

*F*

_{0}+Δ

*F*(

*X*,

_{jF,i}*K*) and χ

_{jF}*=χ*

_{F}_{0}+Δχ(

*X*

_{j}_{χ,i},

*K*

_{j}_{χ}), where

*X*and

_{jF,i}*X*

_{j}_{χ,i}are vectors of wing kinematic variables of the

*i*th trial, the values of which were known from the wing kinematic data;

*K*and

_{jF}*K*

_{j}_{χ}are the regression coefficients for the wing kinematic variables, and

*F*

_{0}and χ

_{0}are the regression constant terms. Δ

*F*and Δχ are the changes of force magnitude and direction due to wing kinematic variables, which were assumed to be linear functions. To estimate the regression coefficients

*F*

_{0}, χ

_{0},

*K*and

_{jF}*K*

_{j}_{χ}from the kinematic data, we performed non-linear least-squares regression using MATLAB Statistics and Machine Learning Toolbox (MATLAB, The MathWorks, Inc., Natick, MA, USA), where the residual sum of squares (RSS) is minimized: (8)where

*m*is the number of trials and

*y*is the square of forward velocity of the

_{i}*i*th trial. Note that in this regression process, we standardized each variable by subtracting its mean and then dividing by its standard deviation. This standardization rendered all variables on the same metric so that the regression coefficients were not influenced by the variables' standard deviations (O'Rourke et al., 2005). We assumed that the changes of force magnitude (Δ

*F*) depend on eight kinematic variables (out of the nine variables mentioned above), i.e.

*X*=[

_{jF}*n*,

*T*

_{d}/

*T*

_{u}, Φ, Ψ, Θ, , , ]

*; and the changes of force direction (Δχ) also depend on eight kinematic variables, i.e.*

^{T}*X*

_{j}_{χ}=[β+χ,

*T*

_{d}/

*T*

_{u}, Φ, Ψ, Θ, , , ]

*. Note that seven out of the nine variables are shared between Δ*

^{T}*F*and Δχ, while the stroke plane angle (β+χ) exclusively affects the force direction, and wingbeat frequency (

*n*) exclusively affects the force magnitude. As a result, a total of 16 variables were included in the variable force-vectoring model, which led to 65,536 model combinations: and one with constant terms only as stated in Eqn 9.

Removing the contribution of wing kinematics, the variable force-vectoring model reduces to the constant force-vectoring model where only body pitch χ contributes to the force vectoring; this can be written as: (9)where both the total aerodynamic force and the angle between the body and the aerodynamic force become constant.

### Model selection and variable importance

The complexity of the variable force-vectoring model depended on the number of wing kinematic variables used. It is well known that models with overly large numbers of parameters suffer from overfitting that could lead to overinterpretation of the data (Burnham and Anderson, 2003). Therefore, model selection using Akaike’s information criterion (AIC) (Akaike, 1998) was performed to evaluate the trade-off between the goodness of fit and model complexity. AIC is defined as:
(10)where *K* is the number of parameters in a candidate model. As a rule of thumb (Burnham and Anderson, 2003), the small sample size-corrected version of the Akaike information criterion (AICc) is preferred if *m*/*K*<40, which was used here. AICc is defined as:
(11)Then, we calculated the AICc difference (Δ* _{i}*) between the AICc

*of the*

_{i}*i*th candidate model and the minimum AICc of all models (AICc

_{min}): (12)The relative likelihood of the

*i*th candidate model (

*g*

*), given the wing kinematic data*

_{i}*X*and

_{jF}*X*

_{j}_{χ}for the

*i*th model, can be computed as: (13)Next, Akaike weights (

*w*) for all model combinations were calculated to quantify the importance of each wing kinematic variable. The Akaike weight of the

*i*th candidate model (

*w*) is defined as: (14)We then summed the Akaike weights over the subset of models that included the

_{i}*X*variable and ranked the variable importance based on the summation of Akaike weights (

_{j}*w*

_{+}).

In addition to AICc, the Bayesian information criterion (BIC) was also calculated for evaluating the trade-off between the goodness of fit and model complexity (Burnham and Anderson, 2003), which is defined as:
(15)Note that with the natural logarithm of the trial number (*m*), BIC applies a larger penalty compared with AICc to the model complexity when *m* increases, which tends to result in simpler models.

With calibrated damping coefficients , forward velocity (*v*) and the best-approximating model, cycle-averaged thrust (*F*_{t}_{hrust}) and lift (*F*_{l}_{ift}) can be estimated according to:
(16)and
(17)

## RESULTS

### Forward flight speed and its dependency on body pitch angle and aerodynamic damping

Using three angle-pin angles and three dampers (D0, D1 and D2; Table 1), body pitch and aerodynamic damping of the flies were systematically varied. The results showed that for all individuals, forward velocity decreased sharply with body pitch (Fig. 2A), but only decreased slightly with damping coefficients (Fig. 2B), except for the medium damper condition (D1) when χ=22.5 deg. Note that the damping coefficients of the medium (D1) and large (D2) damper cases were increased by 54% and 101% compared with that of the no-damper case (D0) (Table 1).

Although the dependency of forward velocity on body pitch angle and damping coefficients was consistent among individuals, there was also considerable variance of flight speed among individuals. For example, the slowest individual (BBF#1) cruised at a mean speed of 0.59 m s^{−1} in the D0 case at χ=0 deg, while the fastest individual (BBF#3) flew at 1.25 m s^{−1} under the same condition. It is also worth noting that all individuals performed smooth steady forward flight at lower body pitch angles (χ=0 and 22.5 deg). However, at χ=45 deg or above (not reported), the forward velocity was reduced significantly to 0.15±0.06 m s^{−1} and occasionally some flies produced vertical oscillations of the rotating shaft at the beginning of the trials. This was possibly caused by the interaction between the wing lift force that tilted the MAGLEV pivot joint and the magnetic restoring torque due to the misalignment of the pivot permanent magnets and the parallel magnetic field generated by the electromagnets (Hsu et al., 2016). The oscillation usually diminished once steady-state flight had been reached.

### Wing kinematic variables during forward flight

As the body pitch angle, damping coefficient and resulting flight speed changed, there also existed considerable changes in wing kinematic patterns (Figs 3 and 4). We characterized the changes using a collection of nine wing kinematic variables representing the cycle-averaged features. The wing kinematic changes were more strongly correlated with the body pitch angle than with damping coefficients, as only wing deviation had noticeable correlation with the damping coefficients (columns in Fig. 3). The contributions of these kinematic variables to speed control were tested according to the variable force-vectoring model (see next section).

Mean wingbeat frequency (*n*) in each damping condition increased with body pitch angle (rows in Fig. 3), with approximately 11% increase from χ=0 deg to χ=45 deg. Therefore, wingbeat frequency had a clear decreasing trend with forward velocity. Mean wingbeat frequency averaged over all trials was 158.9±16.6 Hz, ranging from 172.4±12.5 Hz (χ=45 deg and D2 damping) to 147.2±18.7 Hz (χ=0 deg and D0 damping). The ratio of downstroke to upstroke duration (*T*_{d}/*T*_{u}) generally increased with body pitch angle (Fig. 4A–C), peaking at 1.084±0.077 at χ=45 deg for D1 damping and reaching a minimum at 0.996±0.065 at χ=0 deg for D2 damping.

Wing stroke amplitude (Φ), the changes of which mainly resulted from the extended excursion of the wing stroke towards the end of downstroke (forward excursion), was largest at χ=22.5 deg (123.3±7.4 deg) and was the lowest at χ=45 deg (111.1±10.6 deg) (Fig. 4D–F). The mean wingtip velocity 2Φ*nR* (where *R* is wing length) dropped from 5.08±0.59 m s^{−1} and 5.01±0.62 m s^{−1} at χ=22.5 and 45 deg, respectively, to 4.64±0.76 m s^{−1} at χ=0 deg.

Wing rotation amplitude (Ψ) increased with body pitch angle (except χ=22.5 deg and D0 damper case). Maximum rotation angle occurred at the end of wing pronation, increasing from 47.5±10.4 deg at χ=0 deg to 59.9±12.8 deg at χ=45 deg (Fig. 4G–I). Minimum rotation angle occurred shortly after the end of wing supination, bottomed at χ=22.5 deg and rose approximately 5 deg for both χ=0 deg and χ=45 deg cases (Fig. 4G–I). Wing rotation angle at down-to-upstroke reversal (supination) tended to have a 10–14 deg delay relative to stroke at χ=45 deg, while it was near symmetric (i.e. in phase with stroke) or slightly advanced at χ=0 and 22.5 deg (except with D0 damping at χ=0 deg) (Fig. 4G–I). Rotation angle at up-to-downstroke reversal (pronation) was advanced at χ=45 deg and was delayed at χ=0 deg and χ=22.5 deg.

Wingtip trajectories at χ=0 deg and χ=22.5 deg took oval shapes and those at χ=45 deg were flatter (Fig. 3A–F for χ=0 deg and χ=22.5 deg; and Fig. 3G–I for χ=45 deg). Deviation amplitude (Θ) decreased with increasing body pitch angle (Fig. 4A–C): 19.6±8.9 deg at χ=0 deg and 22.5 deg, and 12.8±5.8 deg at χ=45 deg. Deviation amplitude (Θ), which is the only kinematic variable that has noticeable correlation with damping coefficient, increased slightly with increasing damping coefficient, for example from 15.6±4.7 deg with D0 to 24.0±10.8 deg with D2 at χ=0 deg. The increasing trend was less significant at χ=22.5 deg and χ=45 deg (Fig. 4A–C). A subtle decrease of wing stroke plane angle (β+χ) was observed from χ=0 deg to χ=45 deg (rows in Fig. 3). The changes were limited, remaining within 44.5±5 deg for all trials.

### Force-vectoring models for speed control and variable importance of wing kinematic variables

Two force-vectoring models for predicting the flight speed of flies were tested through non-linear regression based on the estimated thrust, measured flight speed and a collection of wing kinematic variables (described above). The non-linear regression on the constant force-vectoring model yielded total aerodynamic force magnitude *F*_{0}=2.19×10^{−4} N with 95% confidence interval [1.99×10^{−4} N, 2.40×10^{−4} N] (or 53.2% [48.2%, 58.1%] of mean body weight) and χ_{0}=47.8 deg [45.5 deg, 50.0 deg]. The root mean squared error (RMSE) of the prediction based on the constant force-vectoring model was 0.131 m^{2} s^{−2} with *R*^{2} of 0.71. It can be seen that the residual errors of the constant force-vectoring model were relatively large as forward velocity increased (Fig. 5A).

The contributions of the 16 wing kinematic variables selected were tested in the variable force-vectoring model. Note that mean wingbeat frequency (*n*) and wing stroke plane angle (β+χ) were incorporated exclusively in force magnitude and force direction, respectively; and the other seven wing kinematic variables, i.e. *T*_{d}/*T*_{u}, Φ, , Θ,, , , were included in both force magnitude and force direction. We also confirmed that the multicollinearity bias in the regression process was low, by calculating the variance inflation factor (VIF, where VIF<5; Table S4) (Sheather, 2009). Non-linear regressions were performed on all possible variable combinations. AICc/BIC and Akaike weights (*w*) were computed for model selection and variable importance, respectively. From Fig. 5B, the AICc best-approximating model included nine wing kinematic variables that contributed to the changes in the aerodynamic force magnitude and direction (Table 2); the BIC applied a larger penalty on model complexity, which reduced the variable number to six (Fig. 5B). The AICc best-approximating model gave *F*_{0}=3.04×10^{−4} N with 95% confidence intervals [2.85×10^{−4} N, 3.22×10^{−4} N] (or 73.7% [69.2%, 78.3%] of mean body weight), χ_{0}=51.7 deg [50.6 deg, 52.9 deg], RMSE of 0.097 m^{2} s^{−2} and *R*^{2} of 0.922 (Table S3). The BIC best-approximating model gave *F*_{0}=3.06×10^{−4} N with 95% confidence intervals [2.86×10^{−4} N, 3.25×10^{−4} N] (or 74.2% [69.5%, 78.9%] of mean body weight), χ_{0}=51.8 deg [50.7 deg, 52.9 deg], RMSE of 0.104 m^{2} s^{−2} and *R*^{2} of 0.904 (Table S2). Both the AICc and BIC best-approximating model reduced the residual error compared with the constant force-vectoring model, particularly in the higher velocity range (Fig. 5A). Here, we chose the AICc model as our best-approximating model, which was biased slightly towards the goodness of fit than the BIC and included three additional wing kinematic variables (Fig. 6).

Next, from the summations of Akaike weights (*w*_{+}), the relative importance for the wing kinematic variables is indicated (Fig. 6). Wing kinematic variables that contributed to the force magnitude were: (1) wingbeat frequency (*n*), (2) stroke amplitude (Φ), (3) mean deviation angle , (4) ratio of downstroke to upstroke duration (*T*_{d}/*T*_{u}) and (5) mean stroke angle . Among these variables, force magnitude Δ*F* (in Eqn 8) depended negatively on *T*_{d}/*T*_{u}, meaning a shorter duration of the downstroke period increases the Δ*F*, while a higher wingbeat frequency and amplitude, upward shift of the mean deviation angle (so that wing trajectory becomes more oval) and a dorsal (backward) shift of mean stroke angle all led to a higher Δ*F* (see *K _{jF}* sign in Table 2). Kinematic variables that contributed to force direction were: (1) stroke plane angle (β+χ), (2) mean rotation angle , (3) mean stroke angle and (4) mean deviation angle . Among these variables, force direction (Δχ in Eqn 8) depended positively only on , meaning an upward shift of mean deviation angle (or a more oval wing trajectory) results in a backward tilt of force direction (Δχ increases), while the increases in β+χ (i.e. forward tilt of stroke plane), (i.e. increased pronation/decreased supination) and (i.e. backward shift of wing stroke angle) all result in a forward tilt of the force direction (Δχ decreases) (see

*K*

_{j}_{χ}sign in Table 2).

As the body pitch angle χ increased and flight speed decreased, changes in *n* (increase) and (increase) led to an increase in force magnitude (see effect on force magnitude as χ increases in Table 2), while changes in φ (decrease), (decrease) and *T*_{d}/*T*_{u} (increase) led to a decrease in force magnitude (Table 2). In total, the force magnitude increased slightly with increasing pitch angle as the collective result of all wing kinematic changes. In addition, as the body pitch angle increased, all changes in wing kinematic variables resulted in a forward tilt of force direction (β+χ increases, increases, increases and decreases; see Table 2), thereby compensating for the thrust loss due to the backward force tilt.

## DICUSSION

### Flies control forward velocity using their equivalent of helicopter control

Not surprisingly, the inverse dependency between the forward velocity (*v*) and the body pitch angle (χ) (Fig. 2A) in bluebottle flies is consistent with those observed in other insect species (Azuma and Watanabe, 1988; David, 1978; Dudley and Ellington, 1990; Meng and Sun, 2016; Willmott and Ellington, 1997) and birds (Brown, 1963; Pennycuick, 1968). This suggests that bluebottle flies mainly rely on body pitch adjustment to vector the wing aerodynamic forces to produce thrust and regulate flight speed. However, the current study reveals more intricacies in the force-vectoring and speed control of flies, which show close resemblance to those of helicopters or to the ‘helicopter model’. Helicopters create thrust and pitch moment using cyclic and collective pitch, in conjunction with throttle (Leishman, 2006). Collective pitch modulates the force magnitude via lifting/lowering the swashplate, as it symmetrically varies the blade angle of attack (AoA; i.e. the blade AoA is varied by the same amount over the entire rotor disc). Cyclic pitch modulates the force direction via tilting the swashplate, as it first alters the blade AoA asymmetrically, i.e. the blade AoA is altered by different amounts between the halves of the rotor disc, therefore creating an aerodynamic moment due to the lift imbalance; this moment then tilts the entire rotor disc in its orthogonal direction through precession effect and blade flapping (Leishman, 2006). The tilt of the rotor disc then alters the aerodynamic force direction relative to the helicopter body and finally creates a moment that tilts the helicopter body. In this process, the tilt of the rotor disc is relatively small compared with the tilt of the helicopter body; as a result, the angle between the aerodynamic force and the helicopter body is only modulated within a limited range, and the total vectoring of the aerodynamic force is determined mainly by the body pitch. This gives rise to the so-called helicopter model, where the thrust and speed are mainly determined by body pitch and throttle.

Through examining the variable force-vectoring model, here we show that bluebottle flies flying in the flight mill closely follow the helicopter model. The variable force-vectoring model (Eqn 8) reveals the importance of wing kinematic control in predicting the flight speed, similar to the role of collective pitch, cyclic pitch and throttle of helicopters. Specifically, the magnitude of the aerodynamic force is controlled by mean deviation angle, ratio of downstroke and upstroke duration, and mean stroke angle, which can be seen as the flies' equivalent of collective pitch (Fig. 6). In addition, stroke amplitude and wingbeat frequency also control force magnitude, resembling the function of the throttle or the ‘engine speed’ of the helicopters. The force direction is controlled primarily by stroke plane angle, mean rotation angle, mean stroke angle and mean deviation angle (Fig. 6), which can be seen as the flies' equivalent of cyclic pitch. The results also show that these kinematic variables collectively lead to moderate modulation of force magnitude (Δ*F*=8.5×10^{−5}±7.2×10^{−5} N, or 20.5±17.5% of mean body weight in Eqn 8), but only a minor change in force vector direction (Δχ=3.76±2.77 deg in Eqn 8), which resembles closely the speed control of a helicopter, and therefore confirms the helicopter model.

In summary, our results show that although flies use flapping wings instead of rotary wings, and are capable of large modulation of wing kinematics (Deora et al., 2015), their thrust generation mechanism and flight speed control still conform to the helicopter model with limited change of wing kinematics, at least within the range of speeds achieved while flying in the MAGLEV flight mill. Large modulation of wing kinematics presumably only occurs during short-period transient flight such as landing on the ceiling (P. Liu and B. Cheng, unpublished data) and recovery from extreme perturbations (Beatus et al., 2015). Finally, note that the contribution of body pitch to speed control may be subject to saturation at higher speeds, where wing kinematic modulation becomes the primary mechanism. For example, a recent study (Meng and Sun, 2016) showed that drone-flies can fly at a wide range of speed (3.1 to 8.4 m s^{−1}) for a brief amount of time prior to landing at almost the same body pitch (close to 0 deg), and relatively large changes in wing kinematics are employed by the flies to regulate speed.

### Physical significance of wing kinematic variables identified in the AICc best-approximating model

In the variable force-vectoring model, we identified a collection of wing kinematic variables (Fig. 6), which represent either symmetric or asymmetric changes of wing motion between half-strokes. In our experiments, these kinematic variables changed in response to the changes of body pitch angle; together, they control the flight speed of the flies through altering the aerodynamic force magnitude and direction. The specific roles of each wing kinematic variable in modulating the force magnitude and direction are quantified by the regression coefficients *K _{jF}* and

*K*

_{j}_{χ}in Table 2. The trend of their changes in response to body pitch angle can be quantified through Pearson's bivariate correlation between each wing kinematic variable and the body pitch (Fig. S2), and the signs of the regression coefficients are summarized in Table 2, together with the magnitude of their changes quantified by their standard derivation. With these results, below we discuss the physical significance of each wing kinematic variable in force modulation and speed control.

The magnitude of the aerodynamic force is mainly correlated with the wingbeat frequency, stroke amplitude, mean deviation and stroke angle, and the ratio of downstroke to upstroke duration (Fig. 6 and Table 2). Wingbeat frequency increased by 11% on average from χ=0 deg to χ=45 deg (rows in Fig. 3), indicating that flies increase force magnitude to compensate for thrust loss when the force vector is tilted backward with increasing pitch angle. This relatively small increase of wingbeat frequency is expected for flies with asynchronous power muscles, as the wingbeat frequency is primarily determined by the mechanical properties of the coupled wing–thoracic oscillator, which only permits slight alteration of wingbeat frequency (Bartussek et al., 2013). In addition, flies decrease the stroke amplitude Φ as χ increases, which is accompanied by a backward shift of mean stroke angle , through reducing the forward excursion (Fig. 4D–F). As the increases of both variables increase force magnitude (Taylor, 2001), the decreasing trend of stroke amplitude and increasing trend of mean stroke angle result in opposite effects on force magnitude when body pitch increases (Table 2). The backward shift of mean stroke angle also tilts the force vector forward (Zanker, 1988) (Table 2) and creates a pitch-down torque at high pitch angle (indicating the flies lowered their body pitch to compensate for thrust loss). Flies also increased the duration of upstroke (wings sweep backward), during which the thrust is mainly generated, and reduced the duration of downstroke (wings sweep forward), during which drag is mainly generated; this trend has also been found in bumblebees (Dudley and Ellington, 1990) and hawkmoths (Willmott and Ellington, 1997).

The change in mean deviation angle is also a strong contributor to force magnitude, while also being a contributor to force direction. At higher flight speed (or lower pitch), there was an increase of , mainly resulting from the increase of deviation during downstroke (Fig. 4A–C), which rendered the shape of the wingtip trajectory more oval. The oval shape introduces a velocity component perpendicular to the mean stroke plane: upward during downstroke and downward during upstroke. As suggested by Sane and Dickinson (2001), upward velocity reduces AoA and drag force during downstroke and downward velocity results in an increase in AoA and thrust during upstroke; together, they increase force magnitude.

The direction of the aerodynamic force is mainly correlated with the stroke plane angle, and mean stroke, rotation and deviation angles (Fig. 6 and Table 2). At higher body pitch angles, bluebottle flies increase stroke plane angle (β+χ) and decrease mean deviation angle (rows in Fig. 3), although with small variations (44.5±5.0 deg and 3.4±4.2 deg), to tilt the stroke plane and the aerodynamic forces more forward as body pitch angle increases, thereby compensating for the loss of thrust. At large pitch angles, they also increase the mean rotation angle (shifted forward from 8 deg to 16 deg); this increases and decreases the AoA during upstroke and downstroke, respectively, and tilts the force vector forward (similar mechanisms have previously been reported in hawkmoths and hummingbirds: Cheng et al., 2011; Cheng et al., 2016). In summary, we found all four wing kinematic variables contributing to the force direction tend to compensate for the loss of thrust as the body pitches up.

### Left and right wing asymmetry

In the analysis of wing kinematic variables responsible for forward flight, we averaged left and right wing kinematics. However, asymmetry was observed between left and right wing motion due to the rotational nature of the flight mill. For example, mean right (inner) wing stroke amplitude was 13 deg higher than mean left (outer) wing stroke amplitude. A likely explanation is that the wing asymmetry could be due to the bluebottle flies' tendency to perform corrective yaw turns in response to haltere mechanosensory feedback (Dickinson, 1999; Taylor and Krapp, 2007) or to balance the optical flow experienced by the left and right compound eyes. Note that, because of the use of identical spatial frequency of the grating patterns on the two walls, the outer wall had higher temporal frequency because of the larger radius. It has been shown that honeybees use a ‘centering response’ to mediate the unbalanced optical flow (Srinivasan and Zhang, 2004). Subtle stroke amplitude differences between left and right wings can produce roll and yaw moment (Fry et al., 2003); the 13 deg difference observed here is higher than those observed in free-flight saccades in fruit flies, which is likely to result from either continuous yawing or saturated saccadic responses (because of the tether, the error in yaw or roll control cannot be compensated, and any integral controller tends to saturate the response; Muijres et al., 2015). Nevertheless, the existence of the asymmetry does not prevent us from analyzing the forward flight by averaging the left and right wing kinematics, and interestingly it also shows the potential of using the flight mill to investigate a fly's speed control and yaw responses to (bilaterally asymmetric) visual stimuli together with mechanosensory inputs in steady forward flight.

### The advantages and limitations of experiments using MAGLEV flight mill

In this study, we have demonstrated a novel design and use of a MAGLEV flight mill in recording and studying steady forward flight of bluebottle flies. The flight kinematics of the flies were recorded in a fixed region of the flight mill enclosed by two cylinder walls displaying grating patterns, when the flies flew continuously with a speed resulting from its flight control and biomechanics. The pitch angle and therefore the forward-flight biomechanics of the flies can also be manipulated via magnetic angle-pins that tethered the flies to the rotating shaft, which can be potentially replaced by a micro-ball bearing to enable free body pitch angle. Further, the magnetically levitated pivot joint eliminated the mechanical friction compared with traditional flight mills with mechanical pivots, therefore allowing easier and more accurate calibration and manipulation of aerodynamic damping that the flies had to overcome. Together, this MAGLEV flight mill system enabled semi-free flight experiments of insects in a controlled environment, and it can be potentially exploited further to investigate insects' aerodynamics, dynamics, sensing and control in forward flight.

Compared with the free-flight experiments using wind tunnels, there are also a few limitations of the MAGLEV flight mill. First, because of the existence of aerodynamic damping acting on the rotating rod, which cannot be fully eliminated, the fastest steady-flight speed of the bluebottle flies observed in the experiments was 1.25 m s^{−1}, which is slower than that in free flight. However, flies in free flight hardly reach a steady state. For example, flying in a 1.6 m^{3} chamber, bluebottle flies showed a top speed of 2.5 m s^{−1} and an average speed of 1.3 m s^{−1} with constant acceleration and deceleration (Bomphrey et al., 2009). Second, studies showed that tethering may significantly reduce wingbeat frequency (Baker et al., 1981; Betts and Wootton, 1988; Kutsch and Stevenson, 1981). However, no significant difference in wingbeat frequency was observed between flies in the current study (158.9 Hz on average) and typical blowflies (150 Hz) (Dickinson, 1990). Likewise, the same conclusion was made in beetles' forward flight in another flight mill study (Ribak et al., 2017). Third, the rotational nature of the flight mill causes noticeable bilateral wing asymmetries as described above. While here we considered the effect of these asymmetries to be negligible on the forward speed, they need to be limited by using a sufficiently large rod radius. Nevertheless, as described above, they can potentially be exploited to study the role of visual and mechanosensory feedback. Fourth, estimated by Eqn 17, the average lift (*F*_{lift}) created by the flies in the flight mill was 69.5% of the average body weight, and therefore lower than that required in free flight (a similar result was reported in Ribak et al., 2017). One way to more accurately estimate the aerodynamic forces generated by the fly's wings is to use computational fluid dynamics simulation with detailed reconstruction of wing kinematics and deformation (Liu et al., 2016). Nonetheless, holding the effects of these limitations in check, this device provides an alternative approach to a wind tunnel to study insect forward flight in controlled conditions and has a large potential to be further exploited in the future as a common tool in insect flight research.

Lastly, there are certainly differences in speed control between the free and flight mill-tethered flight, mainly in the flies' force-vectoring capability with or without body pitch. In free flight, the flies can regulate their speed by a combination of body pitch and wing kinematic variables, while in the current study, we constrained the body pitch and only allowed the wing kinematic changes to modulate the speed at each experimental condition. However, it is through constraining the body pitch that we were able to examine the force-vectoring capability of the flies using only wing kinematic changes. Also, the observed speed of the flies in the flight mill was likely a result of visual control under the locomotor constraint described above, and the differences between the free and flight mill-tethered flight may also depend on the visual control process, which we did not attempt to address in this study. For example, one needs to test whether flies attempt to fly at a desired flight speed in a flight mill, and how this desired flight speed depends on the flight conditions. To address the above questions, we are currently examining the visual control process of flies in the flight mill, by manipulating the visual inputs through controlling the rotation of the two (motorized) cylindrical walls. In addition, attaching a micro-ball bearing to the back of the flies can potentially provide a free body pitch when they are tethered to the flight mill. With both the visual input manipulation and the free body pitch, more knowledge can potentially be gained on the speed control of flies.

## Acknowledgements

We thank Carl Delacato for assisting with the design of the MAGLEV flight mill.

## FOOTNOTES

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

**Author contributions**Conceptualization: B.C.; Methodology: S.-J.H., N.T., B.C.; Validation: N.T.; Formal analysis: S.-J.H., B.C.; Investigation: S.-J.H., N.T., B.C.; Resources: B.C.; Data curation: S.-J.H.; Writing - original draft: S.-J.H.; Writing - review & editing: N.T., B.C.; Visualization: S.-J.H.; Supervision: B.C.; Project administration: B.C.; Funding acquisition: B.C.

**Funding**This research was supported by the National Science Foundation (CMMI 1554429 to B.C.).

**Data availability**All data are available from figshare: https://doi.org/10.6084/m9.figshare.7130603

**Supplementary information**Supplementary information available online at http://jeb.biologists.org/lookup/doi/10.1242/jeb.187211.supplemental

- Received June 19, 2018.
- Accepted December 13, 2018.

- © 2019. Published by The Company of Biologists Ltd