Flight control in the hawkmoth Manduca sexta: the inverse problem of hovering

The inverse problem of hovering flight, that is, the range of wing movements appropriate for sustained flight at a fixed position and orientation, was examined by developing a simulation of the hawkmoth Manduca sexta. Inverse problems arise when one is seeking the parameters that are required to achieve a specified model outcome. In contrast, forward problems explore the outcomes given a specified set of input parameters. The simulation was coupled to a microgenetic algorithm that found specific sequences of wing and body motions, encoded by ten independent kinematic parameters, capable of generating the fixed body position and orientation characteristic of hovering flight. Additionally, we explored the consequences of restricting the number of free kinematic parameters and used this information to assess the importance to flight control of individual parameters and various combinations of them.

Output from the simulated moth was compared to kinematic recordings of hovering flight in real hawkmoths; the real and simulated moths performed similarly with respect to their range of variation in position and orientation. The simulated moth also used average wingbeat kinematics(amplitude, stroke plane orientation, etc) similar to those of the real moths. However, many different subsets of the available kinematic were sufficient for hovering flight and available kinematic data from real moths does not include sufficient detail to assess which, if any, of these was consistent with the real moth.

This general result, the multiplicity of possible hovering kinematics,shows that the means by which Manduca sexta actually maintains position and orientation may have considerable freedom and therefore may be influenced by many other factors beyond the physical and aerodynamic requirements of hovering flight.

Animal movement arises from the complex interplay of sensory information processing, appendage actuation and force generation. In flight, wings move in ways that may permit hovering, forward, or even maneuvering flight. The various patterns of wing motions follow from the sensory cues involved in flight control (position and orientation along with their velocities and accelerations) and lead to the forces and torques that underlie different flight behaviors. There are potentially many ways in which wings could be actuated to generate a particular flight behavior, such as maintaining a fixed position and orientation. The extent to which such multiple realizations of actuator motions may satisfy the requirements for a given flight behavior remains unexplored. In this study we treat flight control as an inverse problem and examine the set of kinematic inputs capable of generating a specific flight behavior (hovering) in the hawkmoth Manduca sextaL.

Hawkmoths, like all flying insects, generate locomotor forces by activating the flight muscles to move their wings through the air, generating aerodynamic forces and torques, which support and propel them. Variation in these inputs allows the moth to express many flight behaviors including fast forward flight, odor plume tracking, hovering in front of and feeding from flowers,decelerating upon approach and compensating for any slight perturbations to its position or orientation (Stevenson et al., 1995; Willis and Arbas,1991). Although these hovering flights are readily expressed by laboratory reared hawkmoths feeding from artificial flowers positioned in a flight arena, the kinematic and actuator variation used by the moths to establish and maintain position and orientation are unknown. Moreover, the set of possible patterns of variation is also unknown, increasing the difficulty of making specific, a priori predictions of expected kinematic variation and its likely aerodynamic consequences. In our treatment of hovering flight as an inverse problem, a system with known outputs (hovering flight) but unknown inputs (wing movements), we sought different sequences of wingbeats and types of kinematic variation that allowed simulated hovering hawkmoths to maintain their position and orientation, mimicking the feeding behavior of these animals.

Inverse problems are defined in opposition to forward problems, where forward problems are the computation of an output from known inputs. For example, computing the maximum elevation reached by a cannon ball of known mass fired with a known trajectory and initial velocity is a forward problem. The reverse of this, computing the initial velocity of the cannon ball from its maximum elevation (along with its mass and initial trajectory) is an inverse problem. These definitions for forward and inverse problems rely on the notion of a process or system, usually time related, that converts inputs to outputs. In the earlier example the process is gravity. In the inverse problem of hovering, the process is the aerodynamic forces generated by the motion of wings in the fluid. The inputs are therefore the wing motions and the output changes in position and orientation. Inverse problems, including the inverse problem of hovering flight, are often ill-posed. Ill-posed problems are defined in opposition to well-posed problems, where well-posed problems meet the following criteria: (1) a solution exists, (2) the solution is unique and (3) the solution depends continuously on the data(Hadamard, 1902). Ill-posed problems fail to meet one or more of the three conditions. For example, the inverse problem of the cannonball's initial velocity described earlier is ill-posed if the mass and initial trajectory are not known.

While inverse problems and inverse approaches are common in some fields of biology, especially metabolic and gene network analysis (e.g. Holter et al., 2001),application to biomechanics has largely been limited to inverse dynamics analyses of legged locomotion (e.g. Winter, 1990). Those studies of terrestrial locomotion seek to reconstruct joint torques, muscle forces and even neural inputs from whole limb or whole body ground reaction forces and typically reach a single solution. Models incorporating neural inputs are particularly prone to becoming ill-posed(Hatze, 2000), but any system with more degrees of freedom in the inputs than the outputs will likely result in an ill-posed inverse problem.

In some cases ill-posed problems can be made well posed viaregularization, a process for adding assumptions about the solution(Tikhonov, 1963). For example,a recent study of simulated bipedal running and walking reduced a multiplicity of solutions for each locomotor speed to the solution that minimized energy use (Srinivasan and Ruina,2006). While there are myriad constraints one could imagine that may be brought to bear on such ill-posed problems (e.g. energy, target position and its derivatives, stability to perturbations, ranges of acceptable solutions, evolutionary history), the logical first step is to examine the simplest problem of target tracking in flight control. We anticipated that the inverse problem of hovering would have multiple solutions and that there may be many different sequences that would suffice, even sequences that change different aspects of the simulated moth's wingbeat kinematics. However, the factors that have led Manduca sexta to use some limited portion of this set are unknown and cannot be considered without some knowledge of the larger realm of potential solutions.

We examined the inverse problem of hovering by constructing a redundantly actuated (more kinematic parameters than may be needed) mathematical model of Manduca sexta, computing the aerodynamic forces and torques that arise from wing motions and applying these forces to the instantaneous centre of mass of the animal. Wing and abdominal motions were specified by ten independent parameters. We then used a microgenetic algorithm (μGA, see Materials and methods) minimization approach to find individual sequences of parameters, and thus wing motions, that allowed the simulated moth to satisfy our definition of adequate hovering by maintaining position and orientation within specified bounds. The overall flight performance of these simulations was then compared to that of several hawkmoths feeding from artificial flowers in a laboratory flight chamber. Given the redundantly actuated nature of the moth simulation, we expected that many different combinations of wingbeat kinematic parameters would allow the simulated moths to achieve hovering performance that met or exceeded that of the laboratory moths, although the number of combinations and the number of free parameters necessary to achieve performance equivalent to the real moth was not known. After this initial examination of redundantly actuated simulations, we restricted the simulation to fewer free kinematic parameters, moving from redundantly to fully actuated(a sufficient number of free kinematic parameters) and under-actuated (an apparently insufficient number of free kinematic parameters) cases to assess the relative importance of different kinematic parameters and reveal different combinations that might be employed by real moths. We also extended the simulation to level, forward flight at 3 m s^-1 to assess whether the relative importance of different kinematic parameters varies with flight behavior and find the degree of kinematic flexibility required to shift between behaviors.

We examined the inverse problem of hawkmoth flight by creating a forward dynamic differential equation model of the moth with ten different wing and body control parameters then finding parameter sequences that allowed a simulated moth to match the flight behavior of a hawkmoth flying under laboratory conditions. A forward dynamic model computes the state at a specific future time from a set of initial conditions and parameters. In this case, our model computed the simulated moth's instantaneous linear and rotational position, velocity and acceleration from the initial value of those variables and the ten kinematic parameters. We generated individual parameter sets with a microgenetic algorithm (μGA, see below), which broadly searched the available parameter space and identified local minima with sufficient precision to allow a simulated moth to maintain the specified flight behavior. We subsequently computed a large library of parameter sequences by running the simulation with different initial values for the random number generator that underlies the μGA search algorithm. Each different initial value potentially results in a different output from the search process, depending on the number of different acceptable solutions. Finally, we restricted the number of free parameters in the simulation by replacing specific parameters with their mean value from the initial results library and holding these values fixed. We then computed a second library with all possible combinations of free and restricted parameters. At all stages we compared the model's performance with that of real hawkmoths flying in laboratory conditions.

The performance of real hawkmoths was characterized by analysis of high-speed video recordings taken from three individuals hovering in front of and feeding from an artificial flower. The moths were from a colony maintained at the University of Washington and from the same population used to provide the morphological data given below. The flower was positioned in an infra-red illuminated 86 cm×53 cm×87 cm flight chamber and a lateral view of the flight recorded with a high speed video camera (Troubleshooter, Fastec Imaging, San Diego, CA, USA) operating at 250 Hz. We selected a minimum of 1.8 s of particularly steady hovering from the records for each moth, digitized the eye and the tip of the abdomen in each video frame, and then used the two points to calculate pitch. We then assumed that the centre of mass was positioned at a point on the body one third of the distance from the eye to the tip of the abdomen (see below, calculations of moments of inertia), and characterized overall flight performance by the average distance from the centre of mass to its mean location for the entire trial. This simple analysis provided sufficient information to compare the overall performance of the simulated and real moths. A kinematic study detailed enough to record flight kinematics at the same level of detail as the simulation, which uses ten independent parameters to describe each wingbeat, was well beyond the scope of this study.

The morphological parameters employed in the simulation were gathered from 10 adult individuals taken from a colony maintained in the Department of Biology at the University of Washington. The values collected from this individual were similar to those reported in prior studies (e.g. Willmott and Ellington,1997a). We assumed that thickness and density were constant throughout the wing and that moths are bilaterally symmetric.

The simulation used three coordinate systems, beginning with a wing coordinate system X_wY_wZ_w with the origin at the wing root, +X in the anterior direction, +Ylateral to the left and +Z upwards(Fig. 1). Wing positions were first computed in the wing coordinate system according to non-dimensional stroke time (see below), then translated into an anatomic coordinate system X_bY_bZ_b with the origin at the moth's centre of mass (computed from body segment and wing position, see below), +X in the anterior direction, +Y to the left and +Z upwards. Finally, the simulated moth was positioned in an inertial coordinate system XYZ, a right hand, Earth-fixed system with +Z upwards.

where Φ_t̂ is the azimuthal angle at non-dimensional time t̂,Φ _A is the azimuthal amplitude,

is the average azimuthal angle,δ _Φ is the azimuthal phase,θ _t̂ is the elevation angle at non-dimensional time t̂,θ _A is the sweep amplitude, is the mean sweep angle andδ _θ is the sweep angle phase. Note that these angles describe wing position in the wing coordinate system, an anatomic reference frame, and are not in relation to any overall stroke plane angle.

where α_t̂ is the wing rotation angle at non-dimensional time t̂, α_A is the amplitude, δ_α is the rotation angle phase offset and

is the mean rotation angle during the stroke. A rotation angle of 0° would place both the leading and trailing edges of the wing in the horizontal plane of the body coordinate system X_bY_bZ_b with the trailing edge posterior. An angle of 180° would place the trailing edge anterior to the leading edge. The span axis rotation angle matches the longitudinal rotation angle α_sp used in a prior kinematic study of freely flying Manduca(Willmott and Ellington,1997a). The π/2 term in Eqn 3 allowed the hyperbolic tangent function to more closely match the wing rotation patterns observed in that study.

where P_t̂ is the final computed value of one of the wing kinematic parameters defined in Eqn 1-3, P_t̂n-1 is the specified value of the parameter in the prior wingbeat and P_t̂n is the specified value for the current wingbeat. The choice of smoothing function is rather arbitrary as there are many possible forms such as such as Bezier, NURB or other splining methods. It is also possible to enforce continuity by constraining the position and derivatives to identical values for the start and end of each stroke (e.g. Deng et al.,2003), but this approach reduces the range of motion available to the simulation. The hyperbolic tangent has the advantage of an easily controlled decay about the transition point and is also differentiable. An example of the need for, and consequences of, such smoothing is seen in a plot of the temporal pattern wing elevation as the mean of the elevation parameter varies between two successive strokes (Fig. 2). The absence of smoothing leads to a significant discontinuity,making the motion physically and physiologically less reasonable.

where β_t̂ is the abdominal angle at non-dimensional time t̂,β _n-1 is the final abdominal angle from the previous wingbeat,and β_n is the final abdominal angle for the current wingbeat.

where M is the mass of the specified body segment, r is the radius of the head or thorax as indicated by subscripts, a_abdo and b_abdo are the semi-major and semi-minor radii of the abdomen, and D is the distance from the centre of mass of a segment to the whole body centre of mass. Values for the various masses and dimensions are in Table 1.

Our model included aerodynamic forces that arise from wing translation,wing span axis rotation, wing added mass and body drag. Wing forces were computed from a blade-element model using experimentally derived force coefficients. Typical simulation trials employed five elements per wing. To determine the effect of the number of blade elements on the simulation output we recomputed a typical wingbeat with 50 elements per wing rather than five. As expected, the instantaneous forces and moments were slightly greater in the 50-element case; instantaneous moments differed by a maximum of 3.0% and instantaneous forces by a maximum of 2.5%.

where α′ is the angle of attack of the wing in XYZ. Eqn 12 was derived from lift and drag polar plots for low to medium (150-8000)Reynolds number model flapping insect wings(Dickinson et al., 1999; Sane, 2003; Usherwood and Ellington,2002). The resulting large force coefficients and consequent large F_trans are the result of dynamic stall and formation of a leading edge vortex on the translating wing(Bomphrey et al., 2005; Ellington et al., 1996; Liu et al., 1998). Although these effects are typically described as `unsteady' in the sense that they require wing motion and cannot be replicated by a wing held steady in a constant flow, they do not depend on the time history of stroke or of prior strokes and can therefore be captured by `quasi-steady' approximations such as Eqn 12.

where U_T is the instantaneous velocity of the wing tip,

is the wing's instantaneous span axis rotational velocity, c̄ is mean chord length, R is the wing length and ĉ is the non-dimensional chord length (see Ellington, 1984).

where

and

\({\ddot{{\varphi}}}\)

are the wing's overall instantaneous angular velocity and acceleration, α is the wing's span axis angular position,

\({\dot{{\alpha}}}\)

is the wing's span axis rotational velocity and

\({\ddot{{\alpha}}}\)

is the wing's span axis rotational acceleration.

where τ_wing is the net aerodynamic torque from the wing,(x_i, y_i, z_i) is the position of the ith element.

Body and wing forces were similarly summed to give

⁠.

where

is the moth's position vector in the global coordinate system XYZ,

\({\dot{{\chi}}}\)

is the moth's velocity vector,

\({\ddot{{\chi}}}\)

is the moth's acceleration and

\({\vec{g}}\)

is the gravitational acceleration vector (0,0,-9.81), M is body mass,

\({\vec{{\omega}}}\)

is the vector of orientation angles vector,

\({\dot{{\omega}}}\)

is the orientation angular velocity and

\({\ddot{{\omega}}}\)

the orientation angular acceleration. In this study the simulated moth's wing kinematics were symmetric about the X_b axis, restricting it to three degrees of freedom: movement along the X and Zaxes and changes in pitch orientation. Thus we only solve a system of equations consisting of two translational and one rotational degree of freedom.

Although our simulation of hawkmoth flight makes many simplifying assumptions, especially in the calculation of aerodynamic forces, recent experiments conducted on mechanical models of flapping flight (e.g. Dickinson et al., 1999; Sane and Dickinson, 2002; Usherwood and Ellington, 2002)support the expectation that it provides a reconstruction of the aerodynamic forces adequate to the task at hand. Specifically, our simulation does not incorporate aerodynamic forces due to wing-wake interaction or wing-wing interaction (clap and fling), or consider the impact of flexible wings. As Manduca are rarely observed to execute a clap and fling type wingbeat, the absence of wing-wing interaction forces is unlikely to influence our results. While forces generated by flexible wings interacting with the fluid may be substantial (Daniel,1988; Daniel and Combes,2002), there is no direct evidence that they contribute to the forces supporting or propelling flying insects. Furthermore, there is evidence that much of the visible bending in Manduca wings is due to inertial forces and not interaction with the fluid(Combes and Daniel, 2003). Finally, forces attributable to wing-wake interactions likely occur in hovering hawkmoths as the absence of forward body motion provides greater opportunity for the wings to interact with their own wake. However,experiments with a mechanical flapper show that although the wing-wake interaction forces have a large magnitude at certain points in the wingbeat,their overall contribution to the lift and drag impulses is less than that contributed by the forces due to wing translation and rotation(Sane and Dickinson, 2002). As such, although our simulation would not doubt benefit from the inclusion of wing-wake interaction forces, a quasi-steady formulation for these effects has yet to be developed. All of these forces could be calculated with a full Navier-Stokes computational fluid dynamics model similar to several others published recently (e.g. Liu et al.,1998; Ramamurti and Sandberg,2002; Wu and Sun,2004). However, this would also require many orders of magnitude more computational time, precluding its use in a parameter search study such as this one. On the whole, we found that the quasi-steady model using experimentally derived force coefficients for forces due to wing translation and rotation, along with added mass and body drag, incorporated sufficient aerodynamic detail to generate flight behavior and kinematics quite similar to those recorded from real moths.

The flight path tracking simulation functioned by comparing the moth's current location with the desired location at the end of each wingbeat. A microgenetic algorithm (μGA) was used to find wingbeat parameter sets that minimized the difference between the actual and desired locations(Krishnakumar, 1989). Genetic algorithms are a class of biologically inspired algorithms that use the familiar concepts of selection, mutation and recombination to broadly search a parameter space for minimum (or maximum) values. Micro genetic algorithms rapidly reintroduce genetic variation following a selective sweep. As shown in Fig. 3, the μGA functioned by taking an initial population of potential wingbeats (determined by the parameter sets), generated from random permutation of the prior wingbeat,computing their outcome using the forward-dynamics model outlined above, then scoring each individual wingbeat parameter set based on adherence to some predetermined general performance criteria (movement toward the target point in this case). After the individuals are scored, those that most successfully met our selection criteria were used as a basis for the next generation. This generation is created by recombination (taking individual parameters from wingbeats that scored well and combining them to create a new parameter set)and by mutation (randomly changing some of the parameters of high scoring individuals). This new population is then fed back into the loop. We ran the genetic algorithm for 50 generations with 20 individuals each generation, then took the final best available parameter set and applied the Nelder-Mead simplex search algorithm (Nelder and Mead,1965) to the μGA output to further refine it. While the μGA is capable of moving from a poor local peak to a better one, the simplex algorithm is an efficient method for finding the best local result. Appending the simplex search to the μGA output ensures that the combined algorithm reaches the local minima in the region identified by the μGA, a point that the μGA may not reach in a finite number of generations. The combinedμGA and simplex search process uses random mutation and recombination to find a good result. Therefore, running the search twice will not always result in the same answer, especially when there are many answers of similar quality available.

We implemented the μGA by modifying a portion of the MATLAB R14sp3(Mathworks, Natick, MA, USA) genetic algorithm toolbox to follow microgenetic methods, detecting and removing premature fixation of parameters in the population (Krishnakumar,1989). Parameter mutations were uniformly distributed within interval defined by the value of the parameter in the prior wingbeat ±½ the distance from minimum to maximum value for that parameter. For example, the bounds of the wing sweep angle offset were -π/4 and π/4, so if the prior value of the sweep angle offset was 0, mutations would be uniformly distributed within the interval -π/8 to π/8. In each generation, four of the 20 individuals were mutants of the best individual from the prior generation. Within each mutant, each parameter had a 50% chance of mutation. These μGA parameters reflect a broad search strategy that is relatively slow to improve upon an already good solution. However,incorporating the simplex search ensures that it will find the best local solution. We found that parameters rarely changed at the maximum allowed rate.

The core differential equations describing the simulated moth's change in position and orientation from wingbeat to wingbeat (eqn 17-22) were solved at each step via the MATLAB ode45 function, an implementation of an 8th order Dormand-Prince ordinary differential equation solver with the absolute error tolerance set to 1e^-6, approximately 0.003% of the typical value of the smallest of the motion parameters.

Constraints to the parameters prescribing the moth's wing movements were implemented at two levels. Those operating on a single parameter, such as the allowable range for one of the kinematic parameters, were enforced within theμGA mutation and population selection functions by avoiding parameter sets that violated the constraints. A further constraint preventing the intersection of the left and right wings operated within the simulation program, detecting flagging violations for removal.

After using the combined μGA and simplex search algorithms to compute a parameter set for an individual wingbeat, the simulation solved the system of differential equations describing angular and translational momenta, and moved on to the next wingbeat, recalculating the difference between actual and desired locations and searching for the next wingbeat. The final parameter set from the prior wingbeat was used as the starting parameter set for the next wingbeat.

We first used the simulation to examine steady hovering at a target location with a constant pitch angle, much like the behavior associated with nectar feeding. The initial conditions of the simulation placed the moth at the target location but with a forward velocity of 40 cm s^-1, a downward velocity of 5 cm s^-1 and an upward pitching velocity of 180° s^-1. For each trial the simulation ran until the moth either left the target region (a cube with 16 cm edges centered on the target location) or completed 41 wingbeats. We built an initial library of trials with 10 free parameters from a set of 30 simulations, each started with a different initial random seed. We then examined the effect of the number and identity of the free parameters by running three simulations for each possible combination of restricted and unrestricted parameters (2¹⁰-1 combinations), fixing restricted parameters at their mean value in the initial library of trials with 10 free parameters. Although the ill-posed nature of the inverse problem of hovering precludes traditional sensitivity analysis,these parameter restriction trials provide some of the same information.

We also simulated a 3 s long forward flight at 3 m s^-1 with an initial forward velocity of 300 cm s^-1, a downward velocity of 5 cm s^-1 and an upward pitch velocity of 180° s^-1. As before, we first computed an initial library of 30 trials with all parameters free and then the set of restricted parameter sets with three or fewer free parameters. Restricted parameters were fixed to their mean value in the unrestricted trials. Simulations were distributed over a 16-processor computer cluster via the MATLAB Distributed Computing Engine.

We found that hovering hawkmoths maintained their position within a sphere of 0.65 cm radius (and centered on their average location) during the recorded hovering sequences. During these sequences, pitch angle was 41.1±5.8° (mean ± s.d.). One of the sequences recorded in this analysis is reproduced in Fig. 4 for comparison with the simulated moth results.

Flight kinematics from the simulated hawkmoth generally match those of a real moth; an animation of the simulated moth is visually similar to high speed video recordings of a real moth (supplementary material, movie 1). The simulated hawkmoths for which all 10 control parameters were available also performed well; their average position was 2.62 cm from the target location during the thirty 1.5 s simulated hovering bouts. Their average body angle deviated by 1.2° from the overall mean of 33.8°(Fig. 5A). The simulations also made use of all available free kinematic parameters, although some differences in the rate and magnitude of variation are apparent upon visual inspection(Fig. 5B). Moreover, the simulations maintained adequate hovering performance with a wide variety of parameter sequences and combinations (Fig. 6), showing that a large range of parameter combinations can give rise to a similar behavior. The differences between the parameter sequences shown in Fig. 6 arise from the random nature of the μGA combined with the redundantly actuated nature of the simulation with ten free parameters. Each of the simulations started with the same initial conditions and even the same set of kinematic parameters that were used to seed the first generation in the μGA. Different parameter sets arise when the μGA encounters different parameter combinations with similar outputs, leaving the exact combination chosen open to influence by the slightly random nature of the μGA search. Some effects of this are apparent; several kinematic parameters (including δ_φ,θ _A and δ_θ) follow broadly similar patterns for the first few seconds of flight where presumably there are few adequate parameter sets. However, even these similarities disappear within a few wingbeats. Because the simulation follows the μGA with a simplex minimizer, each of these parameter sets represents a local minimum, so the different sequences reached by the simulation represent different local minima, each one adequate for hovering flight.

The simulation does vary all of the available kinematic parameters, but it is not clear from these results which parameters are most important for hovering flight and which are superfluous. To examine this question we restricted variation in some of the kinematic parameters, fixing them at their average values from the unrestricted trials whose results are shown in Fig. 6. We systematically examined the effects of restriction by running three simulations for each possible combination of restricted and free parameters; performance under these varying conditions is summarized by the number of restricted parameters in Fig. 7. This slowly shifted the simulated moth from what we term a `redundantly actuated animal' (all parameters), to one that is fully actuated (a sufficiently large subset of the parameters for successful flight), to an under-actuated one (too few to represent successful flight). As before, performance was quantified by the average distance from the simulated moth to the target point. In some cases the simulated moths also failed to stay within the volume we used to define successful hovering and we also quantified performance by the number of wingbeats prior to any departure from this volume. As Fig. 7B shows, removing just one of the parameters was sufficient to compromise flight performance,increasing the simulated moth's average distance from the target and reducing the number of successful wingbeats. However, these results also show that the simulation is able to maintain hovering for the full 1.5 s in an under-actuated case with up to eight restricted (and two free) kinematic parameters. The minimum average distance from the target point was similar for all redundantly and fully actuated cases, those with up to seven restricted parameters (Fig. 7B). These results could be explained by the presence of a few crucial parameters, whose absence greatly impacts on performance even if all other parameters are available for control, but whose presence is sufficient for adequate flight. To explore this possibility we collected all successful simulated flight bouts for fully and under-actuated trials and from the systematic parameter restriction data set shown in Fig. 7, a set of 21 different combinations, then ranked them by their distance from the target point (Table 2). We found that several kinematic parameters were particularly important to hovering: the sweep angle phase was free in 17 of the 21 successful parameter sets, the mean azimuthal and rotation angle phase were both free in 9 of the 21. Additionally, we found that, as suggested by the change in minimum mean distance between 7 and 8 restricted parameters in Fig. 7B, several of the fully actuated combinations perform better than the only successful under-actuated combination. However, we also found that some of the parameters were completely ineffective: the rotation angle phase and abdominal angle were only valid when matched with the pair that made up the under-actuated combination and did not improve performance compared to the under-actuated case.

We further investigated the validity of the apparently successful under-actuated combination (δ_θ and α_A)by running 20 hovering simulations of 15 s each with these two free variables and found that these two parameters were sufficient for long duration hovering flight (Fig. 8). The position and orientation of the moth in these trials varies more widely than that of the unconstrained simulations (Fig. 8A): mean distance from the target region was 3.16 cm while the mean pitch angle was an average of 6.0° distant from its overall mean of 33.7°. Both free parameters varied at a high frequency throughout the simulation (Fig. 8B), while variation in velocity and orientation vary at a lower frequency. Pitch angle and X-axis velocity appeared to be well correlated, with the simulated moth pitching downward while flying forward(Fig. 8C). A cross-correlation analysis of the two confirmed a strong relationship with changes in the pitch angle preceding changes in the X-axis velocity by two wingbeats (peak cross-correlation of -0.807 at a lag of -2; cross-correlations were calculated from mean removed signals standardized to an autocorrelation of 1.0 at a lag of 0).

Finally, while these results to this point show that there are many parameter sequences and redundantly or fully actuated parameter combinations that lead to similar outcomes and adequate hovering performance; this is not necessarily the case for other possible flight behaviors. To test the generality of our results we ran a series of 3 s long forward flight simulations, targeting the moth to maintain a flight speed of 3 m s^-1 along the X-axis. The general pattern of many solutions persists (Fig. 9),although as might be expected the exact range of parameters adopted and their rates of variation differ somewhat from those used in forward flight(Fig. 5). We also ran three simulations for each fully and under-actuated parameter combination, analogous to the data presented in Table 2. As in the hovering case, several parameters were especially important and there was a single under-actuated combination that met our standard for adequate performance (Table 3). However, the set of important parameters for forward flight was different from those for hovering; the sweep angle mean value(⁠

⁠) increased greatly in importance while that of the sweep angle phase (δ_θ)declined somewhat. The successful under-actuated combination also differed between the two activities.

This study revealed many redundantly, fully and under-actuated sets of kinematic parameters that allowed the simulated hawkmoth to hover, maintaining position and orientation despite an initial perturbation and the small scale variation inherent in the μGA-driven simulation. In many of these cases the simulated moth's performance greatly exceeded our definition of adequate performance and in some cases it exceeded that of a real moth. This general result, the multiplicity of possible hovering kinematics, shows that the means by which Manduca sexta actually maintains position and orientation may have considerable freedom and therefore may be influenced by many other factors. These could include suitability for many modes of flight or minimum deviation from a single, basic wingbeat pattern. However, the simulation results do offer considerable insight into the modes of kinematic variation most likely to accompany hovering flight and the benefits or lack thereof to redundantly, fully or under-actuated control systems in animal flight. The modes of kinematic variation available to real moths are unlikely to precisely match any of the kinematic parameters used in the simulation as the parameters were based on a set of orthogonal axes with no special biological relevance beyond its association with the body axes. Studies matching variation in the neural inputs to the flight muscles, along with the resulting changes in wing kinematics, typically show neural modulation simultaneously influencing several kinematic parameters (e.g. Balint and Dickinson, 2004; Kammer,1971; Wolf, 1990). This mismatch does not present a problem because the aerodynamic effects of kinematic variation do not depend on their source, only on the resulting motion. Moreover, the simulation kinematics are rooted in anatomic coordinate system, so changes to a particular parameter do not depend on body orientation, flight speed or any other external factors and could conceivably be the result of modulation of a small number of neural inputs. Finally, no matter how many of the simulation parameters a single neural input influences,one input can represent only one degree of freedom from an actuation perspective.

As noted above and shown in Fig. 8 and Table 2, we found many kinematic parameter combinations adequate to the task of maintaining hovering flight with three degrees of freedom. The existence of several different three-parameter combinations adequate for hovering flight demonstrates that the simulation with ten free parameters is redundantly actuated, with more free kinematic parameters than degrees of freedom. However, the degree to which it is redundantly actuated is not clear because some of the kinematic parameters may not have any effect, or may exactly duplicate the effects of other parameters. Additionally, Fig. 8 shows that the average performance of the simulated moth dropped with the change from ten to nine free parameters. This suggests that either (1) there is a benefit in this simulation to redundant actuation, such as greater maximum changes in forces and moments or (2) there exists a single, crucial parameter that uniquely actuates one of the degrees of freedom important to hovering flight. Because the minimum mean distance to the target(Fig. 8B) does not change from ten to nine free parameters, and indeed does not begin to change until the step from three to two free parameters, we discount the first possibility mentioned above. The second possibility, that of a single, crucial parameter providing by far the best response in one of the simulation's degrees of freedom, was also poorly supported. Examination of the data underlying Fig. 8 showed that only one of the 30 trials failed before completing the standard 41 wingbeats. The restricted parameter in that case was θ_P, one of the most important parameters for maintaining hovering flight(Table 2). However, because restricting this parameter only led to failure in one of the three trials, it is clearly not the only parameter capable of influencing one of the degrees of freedom. Guaranteed failure within the 41 wingbeat window does not occur until three of the ten parameters were restricted. The sets that then always led to failure were

⁠,

\([{\bar{{\phi}}},{\delta}_{{\theta}},{\bar{{\alpha}}}]\)

and

\([{\bar{{\phi}}},{\bar{{\theta}}},{\delta}_{{\theta}}]\)

⁠. In most cases these parameters were also prominent in the table of successful fully and under-actuated parameter sets(Table 2), and may either provide the only means of influencing one of the simulation's degrees of freedom or be particularly potent actuators along two or three of the degrees of freedom.

Given the large number of independent control inputs available to insects,redundant actuation is likely typical and has been hypothesized to allow finer control (Taylor, 2001). This hypothesis was not directly supported by our results, which demonstrate equally fine control for fully actuated and redundantly actuated systems(Fig. 8B). However, it is possible that this would not be the case if our simulation enforced lower rates of change in the different kinematic parameters. Our comparison of hovering and forward flight showed that redundant actuation may also allow greater flexibility in flight mode. The average kinematic parameters used by the simulation in unconstrained hovering and 3 m s^-1 forward flight simulations differ, implying that there is a performance benefit to changing all parameters when switching between flight modes. Moreover, the fully actuated parameter sets that led to the best performance differ between the hovering and forward flight cases (Tables 2 and 3). Because the utility of individual kinematic parameters depends on flight mode, a redundantly actuated flight apparatus may allow adequate performance in a wider variety of flight behaviors. Redundant actuation may also enhance performance in flight behaviors more challenging than hovering and steady forward flight. For example, hawkmoths not only feed from still flowers, they track the motion of swaying flowers and cope with the changes in local flow environment.

We found twenty different fully actuated kinematic parameter sets that led to adequate hovering flight (Table 2). However, three of the ten kinematic parameters recur with high frequency:

⁠, θ_P andδ _θ, the mean elevation angle, sweep angle phase and rotation angle amplitude, respectively. Each successful set included either the mean elevation angle or the sweep angle phase. Increases in the mean elevation angle

\(({\bar{{\phi}}})\)

move the wing path up the Z-axis and tilt the overall stroke plane slightly forward(Fig. 10) tending to generate additional force along the X-axis, accelerating the moth forward without changing the pitching moment or more than slightly reducing the vertical force. Increases in the phase of the sweep angle(δ_θ) tilt the stroke plane without changing its centre(Fig. 10) and generate a moderate amount of force along the X-axis, accelerating the moth backward while generating a strong upward pitching moment. The influence of the rotation angle amplitude (α_A) is difficult to perceive in the wing tip path, but increases in its magnitude lead to little change in force along the X-axis or pitching moment and an increase in force along the Z-axis, accelerating the moth upward. Although either

\({\bar{{\phi}}}\)

or δ_θoccurs in every successful three-parameter hovering combination(Table 2), they do not allow generation of the same suite of aerodynamic forces and torques and are not directly interchangeable. Instead, they likely make other kinematic parameters useful. For example, δ_θ generates a strong upward pitching moment, which could counteract the effect another parameter such as

\({\bar{{\theta}}}\)

that generates a downward pitching moment in conjunction with additional upward force. Finally, we note that there are influential parameters of every type, including all three of the wing angles considered and the three different types of angle specification: amplitude, mean value and phase. All of these may be of interest in studies measuring kinematic variation or changes in response to stimuli in real animals.

We found a single, under-actuated kinematic parameter set that allowed the simulated moth to hover adequately, both within the initial 41-wingbeat test and a longer 405-wingbeat trial (Fig. 8). In this particular case the two parameters appear to only directly influence the magnitude of the aerodynamic force and torque, using changes in body pitch to orient the aerodynamic force in the X-Zplane. While this did allow the simulation to approximate hovering within the performance bounds we defined, it also led to less accurate positioning around the target point than observed in the actual hawkmoth(Fig. 4). Although hand-tuning the under-actuated control strategy discovered by the μGA might improve performance somewhat, we do not believe that it could approach the performance of the actual hawkmoth, given the very small changes in pitch shown in Fig. 4. Although our simulation approach cannot rule out under-actuated flight control in actual hawkmoth hovering, we believe these results show it is an unlikely possibility.

Several other researchers have made experimental recordings of both wing kinematics and activation potentials from a selection of the flight muscles in Manduca sexta (Frye,2001; Kammer,1971; Wendler et al.,1993). The most complete set of wing kinematic results were reported by Willmott and Ellington, who measured several of the kinematic parameters described earlier on three hawkmoths flying at a range of speeds from 0 to 5 m s^-1 (Willmott and Ellington, 1997a). They reported that hovering hawkmoths used a stroke amplitude ranging from 106.5° to 123.6° with a mean of 116.3°; the average amplitude for the simulated moth in the stroke plane coordinate system was 131.4°. Willmott and Ellington reported that stroke plane angle for hovering moths varied from 11.0° to 27.2° with a mean of 18.2°, the equivalent average for the simulated moth was 22.4°. In the real moth, stroke plane angle and body angle appear to vary inversely in the actual moth and the variation in their sum is less than that in either alone, ranging from 47.1° to 60.6° with a mean of 54.5°; the equivalent mean for the simulated moth was 56.2°. The reported mean elevation and sweep angles convert to a

and of 9.5° and -6.6° for hovering, compared to the average values of 8.7° and -9.9° adopted by the simulated moth. As a reflection of these similarities, the wing tip path of the simulated moth (Fig. 10) generally matches the path recorded from an actual moth shown in fig. 7 of Willmott and Ellington. Both the actual and simulated moths continually varied the kinematic parameters in question, but it is not clear whether the range of variation or its effects are similar between the two circumstances. Willmott and Ellington do not report tabulated data on the wing rotation angle, but visual inspection of their figures suggests an amplitude ranging from 70°to 90°. The wing rotation angles should be directly comparable between studies, the simulated moth used an average rotation amplitudeα _A of 62°. Phase relationships between the different wing angles were not reported, though figs 8 and 10 in Willmott and Ellington suggest that they vary with speed(Willmott and Ellington,1997a).

Many authors have noted insects' apparent use of changes in the position and orientation of the abdomen during flight maneuvers (e.g. Baader, 1990; Camhi, 1970; Götz et al., 1979). These abdominal motions are visually mediated and have been observed in Manduca sexta where there is a correlation between body trajectory and abdominal flexion angle (Frye, 2001). Functional explanations of these movements include changing in the position of the centre of gravity with respect to the wing attachment points and modulating drag on the body. The simulated moth included a hinge at the junction of the thorax and abdomen and was capable of both changing the position of the centre of gravity and the magnitude and moment arm of drag,given a non-zero airspeed. However, the simulated hawkmoth made little functional use of its abdomen (Tables 2 and 3). The aerodynamic forces and moments that result from shifts in abdominal orientation are small in comparison with those brought about by shifts in wing kinematics; this may explain why the model did not make use of the abdomen, but does not explain why actual moths appear to. It is possible that use of the abdomen represents part of a hierarchical control system that cannot be well captured by theμGA used by the simulation to select adequate kinematic parameter sequences. Thus, the abdomen might be used to impart a slight bias to the net torque while variations in wing kinematics act to counter higher frequency changes in orientation. This could allow the abdomen to generate purposeful shifts in position while rapid changes in wing kinematics help maintain position and orientation. The abdomen may also be useful in reducing the rate of pitch generated by wings when there are large vertical excursions of the center of mass. In this case, ventral flexion of the abdomen might be viewed more as an inertial rotational brake than an instigator of pitch motions. Nevertheless, it is clear that the abdomen is not among the few parameters in the under-actuated cases. That said, we do not know the extent to which it plays a role in the fully actuated cases.

Our μGA/simplex algorithm used information about the magnitude, first and second derivatives of body position and angle of a prior wingbeat to guide parameter selections for each successive wingbeat. This assumes that sensory information processing is on the order of a singe wingbeat. Indeed, in work elsewhere we suggest that such delays are problematic if all of the sensory information is delayed by an additional wingbeat(Nishikawa et al., 2006). Thus the potentially long delays in visual motion sensing would lead to a compromised flight performance. However, if pitch angular velocity is not delayed, even though information about all other body kinematic parameters is subjected to delay, the model recovers successful hovering flight. Preliminary data on a potential gyroscopic sensor in moths(Sane et al., 2004) suggest that rotational motions may be encoded with very low (<10 ms) delays. The extent to which sensory information delay interacts with reductions in the number of control parameters remains unexplored.

Treatment of hawkmoth hovering as an inverse problem demonstrated that,like many inverse problems in biology, it is without a unique solution and by definition, ill-posed. Hovering flight's ill-posed nature appeared in our results at two levels. Firstly, different sets of kinematic parameters were sufficient for adequate hovering flight(Table 2). Secondly, even identical sets of kinematic parameters result in slightly different parameter time histories over the course of multiple simulation trials. These arise because near identical solutions exist even within the tightly restricted arena of three (or even two) free kinematic parameters controlling the three degrees of freedom. Although the existence of multiple adequate solutions appears troubling at first glance, it should be taken as a reminder that not all measurable differences propagate. The existence of many functionally equivalent inputs, combinations of kinematic parameters in this case, should be recognized as a possible explanation of variation in experimental measurements, along with the typical interpretations of experimental error and functional difference. As demonstrated in a similar study of human kicking(Hatze, 2000) this tendency toward insensitivity to variation in the inputs becomes stronger when attempting to move from nervous activation patterns to muscle contractions and whole body kinetics.

List of symbols
 
• a_abdo

  semi-major radius of the abdomen

 
• b_abdo

  semi-minor radius of the abdomen

 
• c̄

  mean wing chord

 
• ĉ

  non-dimensional wing chord

 
• C_d

  coefficient of drag of the moth's body

 
• C_r

  coefficient of resultant force for translational motion of the moth's wings

 
• \({\vec{\mathrm{D}}}_{\mathrm{i}}\)

  drag vector for the ith body section

 
• D_abdo

  distance from the centre of the abdomen to the centre of mass

 
• D_body

  magnitude of drag on the moth's body

 
• D_head

  distance from the centre of the head to the centre of mass

 
• D_thorax

  distance from the centre of the thorax to the centre of mass

 
• D_wing,i

  distance from the centre of the ith wing section to the centre of mass

 
• \({\vec{\mathrm{F}}}_{\mathrm{i}}\)

  aerodynamic force vector for the ith wing section

 
• F_acc

  aerodynamic force on the wings due to added mass

 
• F_rot

  aerodynamic force on the wings due to span axis rotation

 
• F_tot

  total aerodynamic force on the wings

 
• F_trans

  aerodynamic force on the wings due to translation

 
• I_abdo

  moment of inertia of the abdomen

 
• I_head

  moment of inertia of the head

 
• I_thorax

  moment of inertia of the thorax

 
• I_tot

  whole body moment of inertia

 
• I_wing

  moment of inertia of a wing

 
• M

  whole body mass

 
• M_abdo

  mass of the abdomen

 
• M_head

  mass of the head

 
• M_thorax

  mass of the thorax

 
• M_wing,i

  mass of the ith wing section

 
• P _t̂

  arbitrary kinematic parameter value at timet̂

 
• P_t̂n

  arbitrary kinematic parameter value from the current wingbeat at time t̂

 
• P_t̂n-1

  arbitrary kinematic parameter value from the prior wingbeat at time t̂

 
• R

  wing length

 
• r̂

  non-dimensional position along the span axis of the wing

 
• r_head

  radius of the head

 
• r_thorax

  radius of the thorax

 
• S_b,i

  surface area of the ith body segment

 
• t̂

  non-dimensional time (fraction of a wing stroke)

 
• U

  speed of the flow past the wing (estimated from the instantaneous velocity of the wing in the inertial frame)

 
• U_t

  speed of the wing tip in the inertial frame

 
• U_b,i

  speed of flow past the ith body segment (estimated from body velocity)

 
• XYZ

  earth fixed inertial coordinate system

 
• X_bY_bZ_b

  body coordinate system (origin at centre of mass)

 
• x_iy_iz_i

  position of the ith wing or body section in the body coordinate system

 
• X_wY_wZ_w

  wing coordinate system (origin at wing root)

 
• α

  wing span axis rotation angle

 
• \({\dot{{\alpha}}}\)

  wing span axis rotation angular velocity

 
• \({\ddot{{\alpha}}}\)

  wing span axis rotation angular acceleration

 
• α′

  wing angle of attack

 
• \({\bar{{\alpha}}}\)

  mean span axis rotation angle

 
• αA

  amplitude of span axis rotation

 
• β

  abdominal flexion angle

 
• βn

  abdominal flexion angle for the current wingbeat

 
• βn-1

  abdominal flexion angle for the prior wingbeat

 
• βt̂

  abdominal flexion angle at timet̂

 
• δα

  wing span axis rotation phase delay

 
• δθ

  wing sweep angle phase delay

 
• δΦ

  wing elevation angle phase delay

 
• θ

  wing sweep angle

 
• \({\bar{{\theta}}}\)

  mmean wing sweep angle

 
• θA

  amplitude of the wing sweep angle

 
• ρ

  air density

 
• \({\vec{{\tau}}}_{\mathrm{bodc}}\)

  torque vector for forces applied to the body

 
• \({\vec{{\tau}}}_{\mathrm{tot}}\)

  net torque vector

 
• \({\vec{{\tau}}}_{\mathrm{wing}}\)

  torque vector for forces applied to a wing

 
• Φ

  wing elevation angle

 
• \({\bar{{\phi}}}\)

  mean wing elevation angle

 
• \({\phi}_{\mathrm{A}}\)

  amplitude of the wing elevation angle

 
• Φt̂

  wing azimuthal angle at t̂

 
• \({\dot{{\varphi}}}\)

  wing instantaneous velocity

 
• \({\ddot{{\varphi}}}\)

  wing instantaneous acceleration

 
• \({\vec{{\chi}}}\)

  moth centre of mass position vector

 
• \({\dot{{\chi}}}\)

  moth centre of mass velocity vector

 
• \({\ddot{{\chi}}}\)

  moth centre of mass acceleration vector

 
• ωP

  moth pitch

 
• \({\vec{{\omega}}}\)

  moth orientation vector

 
• \({\dot{{\omega}}}\)

  moth rotational velocity vector

 
• \({\ddot{{\omega}}}\)

  moth rotational acceleration vector

 
• μGA

  microgenetic algorithm

We wish to thank the Daniel lab members for thoughtful comments on earlier versions of this manuscript and Andrew Mountecastle for providing the high speed video sequences of hawkmoth hovering. The manuscript was greatly improved by comments from two anonymous referees. Funding was provided by an NSF Bioinformatics Postdoctoral Fellowship to T. Hedrick, an Office of Naval Research MURI to T. Daniel, and funds from the Joan and Richard Komen endowment.