SUMMARY
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.
Introduction
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 sexta L.
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 illposed. Illposed problems are defined in opposition to wellposed problems, where wellposed 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). Illposed 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 illposed 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 illposed (Hatze, 2000), but any system with more degrees of freedom in the inputs than the outputs will likely result in an illposed inverse problem.
In some cases illposed problems can be made well posed via regularization, 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 illposed 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 underactuated (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.
Materials and methods
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.
Hovering in Manduca sexta
The performance of real hawkmoths was characterized by analysis of highspeed 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 infrared 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.
Model Manduca sexta morphological parameters
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.
Coordinate systems
The simulation used three coordinate systems, beginning with a wing coordinate system X_{w}Y_{w}Z_{w} with the origin at the wing root, +X in the anterior direction, +Y lateral to the left and +Z upwards (Fig. 1). Wing positions were first computed in the wing coordinate system according to nondimensional stroke time (see below), then translated into an anatomic coordinate system X_{b}Y_{b}Z_{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, Earthfixed system with +Z upwards.
Kinematic model
Our model specified wing and abdominal motions with respect to a nondimensional time parameter t̂ (0≤t̂≤1), which denotes a fraction of the wingbeat cycle. The flapping motion of a moth's wing was prescribed by three angles: azimuthal or elevation (Φ), sweep (θ) and rotation (α) about an axis running along the wing (span axis, Fig. 1). Azimuthal and sweep angles were specified as sinusoidal functions with a mean value, amplitude and phase: (1) (2)
where Φ_{t̂} is the azimuthal angle at nondimensional time t̂,Φ _{A} is the azimuthal amplitude, is the average azimuthal angle,δ _{Φ} is the azimuthal phase,θ _{t̂} is the elevation angle at nondimensional 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.
The span axis wing rotation angle was specified as: (3)
where α_{t̂} is the wing rotation angle at nondimensional 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_{b}Y_{b}Z_{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.
Since the values of the various control parameters may vary between successive wingbeats, it is important to insure continuity of motion between wingbeats with different kinematic parameters. Such continuity was enforced by combining the current and prior wingbeats with a hyperbolic tangent function, smoothing the transition between strokes: (4)
where P_{t̂} is the final computed value of one of the wing kinematic parameters defined in Eqn 13, P_{t̂n1} 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.
The abdomen of hawkmoths constitutes a significant fraction of the total mass of the animal and visuallymediated abdominal flexion reflexes may shift the moment generated by wing forces. Thus our model also included variation in the angle between the thorax and abdomen (Fig. 1A). Abdominal angle was specified as: (5)
where β_{t̂} is the abdominal angle at nondimensional time t̂,β _{n1} is the final abdominal angle from the previous wingbeat, and β_{n} is the final abdominal angle for the current wingbeat.
The instantaneous centre of mass of the moth was calculated from the position of the head, thorax, abdomen and wings. The head and thorax were modeled as spheres, the abdomen as an ellipsoid; wing moments were determined empirically from the mass and shape of a whole wing, assuming constant thickness and density. The moment of inertia about the pitch axis, I_{tot}, was computed as the sum of the moment of inertia of the body segments and wings: (6) (7) (8) (9) (10)
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 semimajor and semiminor 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.
Kinetic model
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 bladeelement 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 50element case; instantaneous moments differed by a maximum of 3.0% and instantaneous forces by a maximum of 2.5%.
The magnitude of the translational force acting on the wing was calculated as: (11)
where ρ is air density (1 kg m^{3}), U is the instantaneous velocity of the flow across the wing (estimated from the instantaneous velocity of the wing in the inertial frame), C_{r} is the resultant force coefficient, r̂ the nondimensional radial position along the wing (equal to the radial distance divided by wing span), R is the wing length, c̄ is the average chord length, and ĉ is the nondimensional chord length (scaled to the maximum wing chord). The resultant force coefficient C_{r} was calculated for each wing segment from: (12)
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 (1508000) 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 `quasisteady' approximations such as Eqn 12.
The magnitude of the force due to rotation of the wing about its span axis was computed from equation 12 in Sane and Dickinson (Sane and Dickinson, 2002), reprinted here for convenience: (13)
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 nondimensional chord length (see Ellington, 1984).
The magnitude of the force due to added mass was computed from Sane and Dickinson (Sane and Dickinson, 2002), equation 3 (with corrections to the sign of the final term and the shape parameter in the final integral, personal communication, S. Sane): (14)
where and are the wing's overall instantaneous angular velocity and acceleration, α is the wing's span axis angular position, is the wing's span axis rotational velocity and is the wing's span axis rotational acceleration.
All translational, rotational and added mass forces generated by the wings were assumed to act normal to the upper surface of the wing at a point¼ chord lengths behind the leading edge. Our assumption of normal forces is well supported by results from mechanical models of flapping flight at a Reynolds number characteristic of Manduca flight (∼8000) (Usherwood and Ellington, 2002) as well as the lower Reynolds number typical of Drosophila flight (∼136) (Dickinson et al., 1999). The assumption that the centre of pressure lies ¼ chord lengths behind the leading edge is longstanding (e.g. MilneThomson, 1973) and, although not directly validated in either mechanical or computational fluid dynamic studies of flapping wings has proved sufficient for estimation of torque asymmetries in a mechanical flapper (Fry et al., 2003). Wing surface orientation vectors were calculated by taking the cross product of a vector along the Xaxis and the Cartesian form of the elevation and sweep angular position of the wing, then rotating that vector about the wing axis by the span axis rotation angleα _{t̂}. The upper surface was designated as the wing surface with a normal vector more than 90° offset from the mean flow vector. The forces acting on each blade element (F_{trans}, F_{rot} and F_{acc}) were summed to calculate the total aerodynamic force (F_{tot}) acting on the wings; this force was multiplied by the orientation unit vector to create a force vector for the ith element. Instantaneous aerodynamic torques were then computed from the position of each wing element and : (15)
where τ_{wing} is the net aerodynamic torque from the wing, (x_{i}, y_{i}, z_{i}) is the position of the ith element.
Our model also included force and torque from drag acting on the moth's body. The moth's body was modeled as a cylinder of diameter 1.1 cm and length 4.6 cm with both linear and rotational motion and, like the wing, divided into a series of individual elements. The magnitude of this drag was calculated as: (16)
where C_{d} is the coefficient of drag (Willmott and Ellington, 1997b), S_{b,i} is the surface area of the ith body segment normal to the estimated incurrent flow and U_{b,i} is the magnitude of the flow past of the ith segment, estimated from the body's rotational and translational velocity. Body drag was directed opposite the velocity vector of each segment giving a drag vector ; torque from body drag was calculated from the cross product of the and the position of each body segment relative to the instantaneous centre of mass: (17)
where τ_{body} is the net torque generated by drag on the body, (x_{i},y_{i},z_{i}) is the position of the ith body segment. The wing and body torques were summed to give the net torque.
Body and wing forces were similarly summed to give .
The instantaneous forces and torques were combined with Newton's laws of motion to form a set of coupled ordinary differential equations describing the changes in rotational and translational momenta: (19) (20) (21) (22)
where is the moth's position vector in the global coordinate system XYZ, is the moth's velocity vector, is the moth's acceleration and is the gravitational acceleration vector (0,0,9.81), M is body mass, is the vector of orientation angles vector, is the orientation angular velocity and 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 Z axes 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 wingwake interaction or wingwing 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 wingwing 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 wingwake 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 wingwake 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 wingwake interaction forces, a quasisteady formulation for these effects has yet to be developed. All of these forces could be calculated with a full NavierStokes 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 quasisteady 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.
Optimization techniques
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 forwarddynamics 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 NelderMead 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 1722) were solved at each step via the MATLAB ode45 function, an implementation of an 8th order DormandPrince 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.
Simulations conducted
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^{10}1 combinations), fixing restricted parameters at their mean value in the initial library of trials with 10 free parameters. Although the illposed 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 16processor computer cluster via the MATLAB Distributed Computing Engine.
Results
Results from animals
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.
Results from simulations
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 underactuated 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 underactuated 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 underactuated 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 underactuated 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 underactuated combination and did not improve performance compared to the underactuated case.
We further investigated the validity of the apparently successful underactuated 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 Xaxis velocity appeared to be well correlated, with the simulated moth pitching downward while flying forward (Fig. 8C). A crosscorrelation analysis of the two confirmed a strong relationship with changes in the pitch angle preceding changes in the Xaxis velocity by two wingbeats (peak crosscorrelation of 0.807 at a lag of 2; crosscorrelations 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 Xaxis. 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 underactuated 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 underactuated 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 underactuated combination also differed between the two activities.
Discussion
This study revealed many redundantly, fully and underactuated 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 μGAdriven 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 underactuated 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.
Redundantly actuated flight
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 threeparameter 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 , and . In most cases these parameters were also prominent in the table of successful fully and underactuated 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.
Fully actuated flight
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 move the wing path up the Zaxis and tilt the overall stroke plane slightly forward (Fig. 10) tending to generate additional force along the Xaxis, 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 Xaxis, 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 Xaxis or pitching moment and an increase in force along the Zaxis, accelerating the moth upward. Although either or δ_{θ} occurs in every successful threeparameter 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 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.
Underactuated flight
We found a single, underactuated kinematic parameter set that allowed the simulated moth to hover adequately, both within the initial 41wingbeat test and a longer 405wingbeat 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 XZ plane. 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 handtuning the underactuated 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 underactuated flight control in actual hawkmoth hovering, we believe these results show it is an unlikely possibility.
Comparison with experimental measurements
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).
Abdominal flexion
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 nonzero 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 underactuated cases. That said, we do not know the extent to which it plays a role in the fully actuated cases.
A dilemma of delay
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.
Inverse and illposed problems
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, illposed. Hovering flight's illposed 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}
 semimajor radius of the abdomen
 b_{abdo}
 semiminor radius of the abdomen
 c̄
 mean wing chord
 ĉ
 nondimensional 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
 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
 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̂n1}
 arbitrary kinematic parameter value from the prior wingbeat at time t̂
 R
 wing length
 r̂
 nondimensional 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̂
 nondimensional 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_{b}Y_{b}Z_{b}
 body coordinate system (origin at centre of mass)
 x_{i}y_{i}z_{i}
 position of the ith wing or body section in the body coordinate system
 X_{w}Y_{w}Z_{w}
 wing coordinate system (origin at wing root)
 α
 wing span axis rotation angle
 wing span axis rotation angular velocity
 wing span axis rotation angular acceleration
 α′
 wing angle of attack
 mean span axis rotation angle
 αA
 amplitude of span axis rotation
 β
 abdominal flexion angle
 βn
 abdominal flexion angle for the current wingbeat
 βn1
 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
 mmean wing sweep angle
 θA
 amplitude of the wing sweep angle
 ρ
 air density
 torque vector for forces applied to the body
 net torque vector
 torque vector for forces applied to a wing
 Φ
 wing elevation angle
 mean wing elevation angle
 amplitude of the wing elevation angle
 Φt̂
 wing azimuthal angle at t̂
 wing instantaneous velocity
 wing instantaneous acceleration
 moth centre of mass position vector
 moth centre of mass velocity vector
 moth centre of mass acceleration vector
 ωP
 moth pitch
 moth orientation vector
 moth rotational velocity vector
 moth rotational acceleration vector
 μGA
 microgenetic algorithm
ACKNOWLEDGEMENTS
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.
FOOTNOTES

Supplementary material available online at http://jeb.biologists.org/lookup/suppl/doi:10.1242/jeb.02363//DC1
 © The Company of Biologists Limited 2006