## Abstract

Previous studies of insect flight control have been statistical in approach, simply correlating wing kinematics with body kinematics or force production. Kinematics and forces are linked by Newtonian mechanics, so adopting a dynamics-based approach is necessary if we are to place the study of insect flight on its proper physical footing. Here we develop semi-empirical models of the longitudinal flight dynamics of desert locusts *Schistocerca gregaria*. We use instantaneous force–moment measurements from individual locusts to parametrize the nonlinear rigid body equations of motion. Since the instantaneous forces are approximately periodic, we represent them using Fourier series, which are embedded in the equations of motion to give a nonlinear time-periodic (NLTP) model. This is a proper mathematical generalization of an earlier linear-time invariant (LTI) model of locust flight dynamics, developed using previously published time-averaged versions of the instantaneous force recordings. We perform various numerical simulations, within the fitted range of the model, and across the range of body angles used by free-flying locusts, to explore the likely behaviour of the locusts upon release from the tether. Solutions of the NLTP models are compared with solutions of the nonlinear time-invariant (NLTI) models to which they reduce when the periodic terms are dropped. Both sets of models are unstable and therefore fail to explain locust flight stability fully. Nevertheless, whereas the measured forces include statistically significant harmonic content up to about the eighth harmonic, the simulated flight trajectories display no harmonic content above the fundamental forcing frequency. Hence, manoeuvre control in locusts will not directly reflect subtle changes in the higher harmonics of the wing beat, but must operate on a coarser time-scale. A state-space analysis of the NLTP models reveals orbital trajectories that are impossible to capture in the LTI and NLTI models, and inspires the hypothesis that asymptotic orbital stability is the proper definition of stability in flapping flight. Manoeuvre control on the scale of more than one wing beat would then consist in exciting transients from one asymptotically stable orbit to another. We summarize these hypotheses by proposing a limit-cycle analogy for flapping flight control and suggest experiments for verification of the limit-cycle control analogy hypothesis.

## 1. Introduction

We can now claim a reasonable descriptive understanding of how insects support themselves in the air. This has been facilitated in large part by the longstanding and enthusiastic adoption of techniques and methodologies first developed in the physical sciences, from flow visualization (e.g. Magnan 1934) to computational fluid dynamics (Smith *et al*. 1996; Liu *et al*. 1998). In light of this enthusiasm, it is all the more remarkable that the rich legacy of a century of research on aircraft flight stability and control has scarcely been considered in the study of insect flight control (Taylor 2001; Taylor & Thomas 2003). Rather than paralleling the aircraft literature in using an approach founded on dynamics, the insect flight-control literature has gone its own way in adopting a largely kinematics-based approach. Whereas a dynamics approach directly links the motions of a system with the forces that produce them, a kinematics approach treats motions in isolation from the forces that produce them. This means that a kinematics approach is necessarily correlative: we know how body kinematics correlate with wing kinematics (Wakeling & Ellington 1997), we know how aerodynamic force production correlates with wing kinematics (Lehmann & Dickinson 1998) and we even know how wing and body kinematics are correlated with the firing of individual flight muscles in free flight (Kutsch *et al*. 2003)—but our understanding of the motions nevertheless rests at the level of statistical correlation, rather than at the level of a physical model of the underlying mechanics. By providing this functional link between forces and motions, a flight dynamics approach provides a new, physically proper, paradigm for understanding insect flight control. The flight dynamics approach also offers a host of new techniques for analysing experimental data. Here, we use a flight dynamics model to predict body kinematics from empirical measurements of instantaneous force production.

Placing experimental studies of insect flight control on a flight dynamics footing is ultimately essential, as it is otherwise impossible to examine how the system will respond to control inputs and external disturbances. This is important, because the end effect of a control input may be quite different from its initial effect. For example, whereas the initial effect of increasing the throttle on an aircraft is a straightforward increase in flight speed, the consequent increase in lift can cause the aircraft to climb until it eventually settles (via a lightly damped oscillation in speed and pitch) to a steeper flight path in which the flight speed is unchanged from its initial value (Etkin & Reid 1996). Such counterintuitive effects as these call for caution in predicting the function of specific control inputs without a proper model of the flight dynamics. We therefore think it completely necessary—and, indeed, completely desirable—that experimental data on insect flight control should be analysed with the aid of equations of motion, which lie at the heart of any physically realistic analysis of flight dynamics. Equations of motion are the appropriate way to analyse experimental data on insect flight control when questions of causation relating to Newtonian mechanics are being asked, just as statistical techniques are the appropriate way to analyse experimental data on insect flight control when questions of correlation are being asked.

Taylor & Thomas (2003) provided the first flight dynamics analysis of experimental data on insect flight control, using a classical linearized framework borrowed directly from the aircraft flight literature to analyse force measurements from tethered locusts. Sun & Xiong (2005) subsequently used the same framework to model the flight dynamics of hovering bumblebees, using a computational fluid dynamics approach to model the aerodynamic forces and moments. In both cases the models failed to explain the insects' flight stability fully, but a key limitation was that the models were time-invariant which meant that the periodic wing beat forces had to be time-averaged before being used in the analysis. Here, we present the previously unpublished instantaneous force–moment data from the experiments by Taylor & Thomas (2003), and analyse them using equations of motion with periodic coefficients. Having measured directly how the instantaneous forces on individual tethered locusts vary with changes in flight condition, we use the equations of motion to predict how the combined passive and active responses, which contribute to the total measured response, determine how the individual insect will respond to disturbance. The motivations for this analysis are twofold: firstly, by comparing the nonlinear models we develop with the earlier linear models of Taylor & Thomas (2003) for the same individual locusts, we are able to assess the importance of incorporating periodicity and nonlinearity into our models of insect flight dynamics; secondly, because the framework we develop is a novel one, driven by the periodic character of the new experimental data we present, we provide the most complete dynamic analysis possible with the available data. The framework we develop here is as complete as we think necessary to make major progress in analysing insect flight dynamics, and by also presenting the mathematical tools for analysis, we hope to inspire a new generation of experimental analyses of insect flight control.

We begin by presenting the Newton–Euler rigid body equations of motion on which the analysis of the experimental data is founded (§2). We then describe briefly how the instantaneous force measurements were collected and give the rationale behind the experiments (§3), before analysing these experimental data and using them to formulate a series of empirical models of force production (§4). By embedding the empirical models of force production in the equations of motion, we are then able to develop a semi-empirical model of the flight dynamics of each locust, which we analyse numerically (§5). Finally, we discuss what the models can tell us about locust flight control, before going on to present some new hypotheses about flapping flight control and suggestions for experiments to test them, inspired by the modelling approach we develop (§6) (table 1).

## 2. Equations of motion

We will assume, for the purposes of this analysis, that a locust can be treated as a bilaterally symmetric body. Pure longitudinal motions are then possible in which all of the forces, moments and motions operate parallel to the plane of symmetry. If we furthermore treat the locust as a single rigid body (which in effect means that we assume the centre of gravity to be fixed and the moment of inertia to be constant), then the Newton–Euler equations of motion for symmetric flight are(2.1)(2.2)(2.3)(2.4)where (*X*, *Z*, *M*) are the force components acting along the *x*- and *z*-axes, and the pitching moment acting about the transverse axis, respectively (see Taylor & Thomas 2003). The four state variables (*u*, *w*, *q*, *θ*) are, respectively, the components of translational velocity along the *x*- and *z*-axes, the angular velocity about the transverse axis and the pitch attitude of the body with respect to the horizontal (figure 1). Body mass *m*, pitching moment of inertia *I _{yy}*, and gravitational acceleration

*g*=9.8066 m s

^{−2}are all assumed to be constants. Note that in this formulation, the flapping of the wings can manifest itself only in periodicity of the forces and moments. This means that the model takes account of wing inertial and aerodynamic forces, but does not account for changes in the moment of inertia and centre of gravity resulting from the wings' flapping. This is not unreasonable, since the wings comprise less than 4% of the total mass of the insect, and is unavoidable here, because it was not possible to collect data on wing kinematics with the experimental setup used. Equations (2.1)–(2.4) are derived theoretically, but they are so general as to be meaningless until they are parametrized empirically. We therefore term the models we develop ‘semi-empirical’ to indicate that while the overall form of the equations is derived theoretically, their parameters are measured empirically.

This deceptively simple set of four coupled nonlinear ordinary differential equations ((2.1)–(2.4)) is not amenable to analytical solution, and conceals a host of complexities since the three force–moment components *X*, *Z* and *M* are in general functionals of the state variables *u*, *w* and *q* (Tobak & Schiff 1981), and may also depend explicitly on time (in which case the system is said to be non-autonomous, in contradistinction to autonomous systems, which are time-invariant). In Taylor & Thomas (2003), equations (2.1)–(2.4) were simplified to express the flight dynamics of desert locusts *Schistocerca gregaria* (Forskål) using the linear time-invariant (LTI) model:(2.5)where the state vector comprised small deviations from the equilibrium values of the state variables (the superscript ‘T’ denotes vector transposition). The constant 4×4 system matrix * A* was given by(2.5a)where the subscript ‘where

*P*(

*t*) is a periodic function,

*ω*is the angular frequency of the wing beat and

*a*and

_{n}*b*are real coefficients. This has the effect of removing any non-harmonic content due, for example, to wind tunnel turbulence. Unfortunately,

_{n}*ω*varied slightly through each 15 s recording, causing phase shifts that resulted in an unacceptable smoothing when

*ω*was treated as constant, as in equation (4.1). Although variation in

*ω*was trivial between wingbeats, the cumulative phase shift over a 15 s recording could be quite large. Conceptually, this can be dealt with by replacing the linear function

*ωt*with a nonlinear phase function

*Φ*(

*t*) that takes account of variations in

*ω*. In practice, because the data are discrete, the continuous function

*Φ*(

*t*) was replaced by a phase vector

*, unique to each recording, and generated by assigning a phase angle to every point in the discrete time domain.*

**Φ**We began by identifying the approximate start and endpoints of successive wing beats by noting when the *Z*-component crossed its mean datum line with negative slope, having first down-sampled the data at 0.5 kHz to remove spurious mean-crossings due to higher spectral content. We then reverted to the original data sampled at 10 kHz and used an exhaustive search procedure to optimize sequentially the exact start and endpoints of successive wing beats within the 0.002 s window set by the resolution of the mean-crossing analysis, minimizing the absolute difference between the first and last *Z*-components of each successive wing beat. We then generated * Φ* as a piecewise linear series, varying linearly from 0 to 2

*π*over the course of each wing beat. Having generated

*for each recording, we next used Matlab's general linear model (GLM) ‘glmfit’ function (Model I regression, e.g. Sokal & Rohlf 1995) to fit a model of the form(4.2)where is discrete. Practically, this amounts to fitting a separate additive GLM for each force–moment component, with predictor variables cos*

**Φ***n*and sin

**Φ***n*where

**Φ***n*=1…8. Using a GLM offers the best possible fit to the data, in the sense of minimizing the residual sum of squares, and has the additional advantage of automatically providing

*p*-values and confidence limits for the Fourier coefficients. Care is needed in the interpretation of the

*p*-values because the

*a*and

_{n}*b*coefficients describe the same

_{n}*n*th harmonic; therefore, the harmonic should be considered significant if either coefficient is significant, irrespective of the significance of the other. For example, if the

*n*th harmonic of a Fourier series is an even function of

*t*—i.e. if

*P*(

*t*)=

*P*(−

*t*)—then only the cosine coefficient

*a*will be significant. We initially fitted models with

_{n}*h*=10 harmonics but since usually only the first seven, or sometimes eight, harmonics were significant (

*p*<0.05), we subsequently fitted a new set of models with only

*h*=8 harmonics, thereby condensing each force–moment component of each 15 s recording into a Fourier series with eight harmonics.

## 4.3 Alignment of Fourier series phase angles

Since the phase vectors * Φ* were generated separately for each recording, the phasing of the Fourier series coefficients is exactly consistent between the various force–moment components of a single recording, but only approximately so between recordings. In the absence of any external reference, the phase of the first harmonic of a Fourier series is arbitrary, and we therefore aligned the overall phase of the various Fourier series by adjusting their phase so that their respective first

*Z*-harmonic had zero phase angle. To do this, we first converted each series to the form(4.3)where the phase angles

*ϕ*are in radians and the coefficients

_{n}*c*are real numbers. We then subtracted

_{n}*ϕ*

_{1}for the

*Z*-component from all of the

*ϕ*for each force–moment component of each recording, so that the phasing of every Fourier series was consistent with that of every other.

_{n}## 4.4 Correction of zeroth coefficients

The locusts exhibited large temporal variation in their quasi-static force production through the course of the experiments. This was controlled for in the previous LTI analysis by pairing each recording made at a perturbed body angle with a recording made under reference conditions immediately after, and using this to correct for temporal variation in quasi-static force production (Taylor & Thomas 2003). This method of correction was found to be highly effective, so we essentially used the same method again here. The zeroth coefficients of the Fourier series correspond to the mean level of force production, and we therefore corrected each *c*_{0} term corresponding to a perturbed body angle by subtracting the *c*_{0} term of the respective paired reference and adding the mean of all the *c*_{0} reference terms. Unfortunately, no such correction was possible for the recordings in which speed was varied.

## 4.5 Reduction of multiple reference recordings to single Fourier series

In the previous LTI analysis, the multiple reference recordings were simply averaged for each locust to give a single mean level of force production under the reference conditions (Taylor & Thomas 2003). Averaging the dynamic forces in the time domain results in an unacceptable smoothing of the data, but averaging in the frequency domain is equally problematic in the issues it introduces in averaging the phase angles. To overcome these problems, we fitted a single Fourier series to all of the reference recordings, having used the phase angle adjustments of the previous section to align the phase vectors * Φ* of the various recordings. To do this, we simply subtracted

*ϕ*

_{1}for the

*Z*-component from the respective phase vector

*, thereby aligning all of the reference*

**Φ***with the first*

**Φ***Z*-harmonic. We then concatenated each of the

*,*

**Φ***Z*,

*X*and

*M*vectors for the various reference recordings, and used a GLM to fit a single Fourier series to the concatenated data for each of the force–moment components.

## 4.6 Periodic models of force production

We used the Fourier series fitted in §§4.2–4.5 to reconstruct the instantaneous forces through a single wing beat for each flight condition. We fitted 500 points per cycle, which equates approximately with the original sampling frequency. The resulting time domain reconstructions are plotted in figures 2–7 as functions of wing beat phase to show the effects of speed and body angle on each instantaneous force–moment component. These reconstructed force traces not only condense the original recordings into a manageable form, but also offer as clean a signal as possible for fitting the coefficients of a periodic model of force production. To do this, we used separate GLMs to fit the first eight harmonics of the Fourier series of the form(4.4)and(4.5)for each force–moment component and locust. Note that because the empirical models of force production in equations (4.4) and (4.5) are fitted to force traces reconstructed over one wingbeat period, they recover the strictly periodic form of equation (4.1)—cf. equations (4.2) and (4.3).

It is important to note at this juncture that, although the experiments do not allow us to distinguish between body angle *α* with respect to the relative flow and pitch attitude *θ* with respect to the horizontal, we have assumed here that the flight forces change only in response to the former. This is always the case for passive aerodynamic effects, but need not be the case for the active control responses of locusts if they can sense pitch attitude directly (e.g. Taylor & Thomas 2003). However, in the absence of any evidence to the contrary, we will use the simplest model possible, treating the flight forces as functions of body angle but not of pitch attitude.

By inspection, equations (4.4) and (4.5) directly yield the coefficients of a linear periodic model of force production of the form(4.6)where(4.7)

(4.8)

(4.9)

It was not possible to test for the existence of possible interactions between *P*_{α}(*t*) and *P _{U}*(

*t*) in determining the total instantaneous force, because this would have required collecting a large number of additional force recordings for each locust which was simply impractical. The additivity implied by equation (4.6) is therefore an assumption of the model, rather than an empirical conclusion. Furthermore, it is worth noting that whereas such additivity is also an assumption of the LTI model of Taylor & Thomas (2003), it is not a necessary assumption of the nonlinear models presented here and can be relaxed as necessary in the future.

Figures 2–7 show the instantaneous forces predicted by the respective periodic models (equations (4.4) and (4.5)) alongside the reconstructed empirical force traces to which they were fitted. In effect, the right-hand column of the graphs plot the regression surfaces fitted to the empirical data plotted in the left-hand column of the graphs. It is immediately clear that the periodic models fit the reconstructed force traces extremely well. The Fourier series coefficients themselves are listed in table 2, converted to (*c _{n}*,

*ϕ*) form for convenience. Coefficients are shown in bold if a harmonic was statistically significant (

_{n}*p*<0.05) in the GLM in which it was fitted. Two-thirds of the models contained significant harmonic content up to the seventh or eighth harmonic, although one (or in one case, two) of the lower harmonics could occasionally be non-significant (see table 2 for details).

The error degrees of freedom were 7466 for the models of *P*(*α*, *t*) and 3466 for the models of *P*(*U*, *t*). Histograms of the residuals and plots of residuals against fitted values showed excellent normality of error and homogeneity of variance, respectively, for the angle series data. The corresponding graphs for the speed series data showed acceptable normality of error, but displayed strong clumping of the residuals by flight condition. This heterogeneity indicates that some component of the systematic variation in the speed series data remains unexplained. The most likely reason for this difference is that whereas we corrected the angle series data for temporal variation in quasi-static force production (§4.4 above), this could not be done for the speed series data. Unfortunately, temporal variation could not be controlled for statistically without eliminating some of the variation due to speed, because speed and time are not completely orthogonal variables. Nevertheless, the clumps of residuals were reasonably evenly distributed for locusts R and B and the overall variance was sufficiently homogeneous that the assumptions of the GLM were broadly satisfied.

In the case of locust G, the residuals of the speed series data showed further patterning, indicating a strong nonlinearity in the underlying data. On inspection of the residuals, this nonlinearity was found to reflect the outlying value of mean force production at 2.70 m s^{−1}. This is clearly visible in the inflection at 2.70 m s^{−1} in the surface plot of figure 7*b* (see also figures 3*b* and 5*b*). The same nonlinearity was also found in the earlier LTI analysis (Taylor & Thomas 2003), and since the data at 2.70 m s^{−1} for locust G therefore lie outside of the range of speeds over which changes in the forces can be satisfactorily linearized, it is reasonable to drop these data from the linear periodic model of force production and to fit a new GLM. This reduced the error degrees of freedom to 2966, but nevertheless caused several harmonics that had not previously been significant to become so, which is further indication of the validity of the procedure.

The analysis provided here assumes that wing-beat frequency is invariant, taking as its value the mean wing beat frequency recorded during the angle series under reference conditions. In fact, the mean of the wing beat period through a recording (*T*) varied slightly with speed or body angle in two of the locusts: in locust R, *T* varied as a quadratic of both body angle (GLM *T*=*α*+*α*^{2}; *p*=0.018; *F*_{14}) and tunnel speed (GLM *T*=*U*+*U*^{2}; *p*=0.005; *F*_{6}); in locust G, *T* increased linearly with flight speed (GLM *T*=*U*; *p*=0.013; *F*_{14}). Nevertheless, the range of variation was relatively small, and *T* always fell within 6% of the mean (usually much less). The opposite dependence of wing beat period on flight speed has been found elsewhere for free-flying *Locusta migratoria* (L.) (Baker *et al*. 1981), and we consider the variation in wing beat period we recorded too slight and inconsistent between individuals to warrant inclusion in our models. Systematic variation in wing beat frequency will alter the system dynamics but we leave the intriguing possibilities this raises for future analyses to consider.

## 4.7 Model verification: statistical analysis

The purpose, and end result, of the above analysis is to condense all of the force recordings for each locust into a single empirical model of periodic force production:(4.10)where and **X**_{ref}(*t*), **X**_{α}(*t*) and **X**_{U}(*t*) are periodic analytic functions of time, for which Fourier series representations are given in table 2. The Fourier series coefficients for **X**_{ref}(*t*) used in the final model (equation (4.10)) are those from the angle series GLMs (i.e. fitted to data concatenated from 14 reference recordings); the equivalent coefficients from the speed series GLMs are less reliable and were therefore dropped. Before using this model in simulating locust flight dynamics, it is necessary to verify how well it predicts the original instantaneous force–moment data. To assess the predictive performance of the model, we first concatenated all of the force traces and phase vectors * Φ* from §4.2 above for each locust. Then, for consistency with the analysis above, we corrected for temporal variation in quasi-static force production by subtracting from the angle series data the mean of the corresponding paired reference and adding back the mean of all of the paired references. We then used equation (4.10) to predict the concatenated instantaneous forces and moments using the concatenated phase vectors. Finally, we fitted a least-squares linear regression of observed and predicted forces for each locust.

The *r*^{2} coefficients of these regressions are all high (see table 3 for detail), indicating that equation (4.10) explains a high proportion of the total variance in the observed forces (never less than 45% but possibly as high as 86%). By way of comparison, table 3 also includes *r*^{2} coefficients from the regression of the observed instantaneous forces on the predicted quasi-static forces, calculated by including only zeroth harmonics. The difference between the *r*^{2} coefficients calculated from this time-averaged model of force production and those already calculated using the full periodic model is the proportion of the total variance explained by the periodic components of equation (4.10) over and above any time-invariant component. This difference is very large for the *Z*-component of force (69–84%), moderate for the *X*-component (14–30%) and small for the *M*-component (4–9%). This systematic variation among force–moment components reflects differences in the scale of the periodic components of the signal relative to the variation in quasi-static force production among recordings. The general conclusion to be drawn from these regressions, however, is to verify that our periodic model of force production performs very well indeed in predicting the observed instantaneous forces, based upon the phase data calculated from the initial mean-crossing and Fourier series analysis of the *Z*-traces. Clearly, there is some circularity inherent in this procedure but this is unavoidable, and the test involves the smallest possible set of assumptions needed to test the model given its intrinsic phase dependence.

## 4.8 Model verification: conditions for quasi-static equilibrium

The statistical analysis in §4.7 verifies that the empirical models of force production we have fitted predict the experimental data well. In order to verify whether the empirical models of force production we have fitted are also consistent with the balance of forces needed for free flight, we calculated the conditions in which each of the locusts would be in exact quasi-static equilibrium. The conditions for quasi-static equilibrium are(4.11)(4.12)(4.13)where the overbar notation denotes quasi-static terms, the coefficients of which are given by the zeroth harmonics in table 2, and where the subscript ‘’ denotes the value of an independent variable at quasi-static equilibrium. Equations (4.11)–(4.13) are formulated by making the quasi-static components of aerodynamic force in the *x*-direction (equation (4.11)) and *z*-direction (equation (4.12)) equal and opposite to the corresponding weight components, and by setting the net pitching moment to zero (equation (4.13)). This set of three algebraic equations contains three independent variables: equilibrium body angle , equilibrium flight speed and equilibrium pitch attitude . Pitch attitude and body angle are only independent if we permit the locust to be in ascending or descending flight at equilibrium. This is a key difference from the LTI analysis of Taylor & Thomas (2003), in which the equations of motion were linearized about an equilibrium condition of steady level flight, leaving only two independent variables and thereby precluding a consistent solution to equations (4.11)–(4.13) or their equivalents. Instead, these had to be solved for a pseudo-equilibrium in which the equilibrium conditions were almost, but not exactly, satisfied (Taylor & Thomas 2003).

Equations (4.11)–(4.13) permit two exact solutions for each locust, of which only one represents forward flight. The exact solutions for forward flight are given in table 4 for each locust, and indicate that the force measurements for locusts R and G are compatible with the locusts attempting a steady ascent (*θ*_{e}>*α*_{e}), whereas the force measurements for locust B are compatible with it attempting a steady descent (*θ*_{e}<*α*_{e}). Note that this result says nothing about the mode of flight the locusts were attempting at any given moment; we have controlled for temporal variation in force production as far as possible, and, therefore, the equilibria we have solved represent the long-term average flight mode through the experiment. In any case, the results indicate, not surprisingly, that the locusts would not have remained flying in a steady level flight condition had this not been enforced by the tether. Nevertheless, it is reassuring that all of the exact quasi-static equilibria for forward flight correspond to lifelike flight conditions, consistent with the ranges of speed and body angle observed in free-flying locusts in the wild (*S. gregaria*: Waloff 1972; *L. migratoria*: Baker *et al*. 1981). There is no intrinsic reason, mathematically, why this should be the case, and the lifelike equilibria for which we have solved testify to the physical realism of the empirical models of force production embedded in the equations of motion.

## 5. Analysis: semi-empirical models of flight dynamics

### 5.1 NLTI and NLTP models

Substituting equation (4.10) into equations (2.1)–(2.3), and making use of the identities and , equations (2.1)–(2.4) expand to give the following NLTP model:(5.1)(5.2)(5.3)(5.4)where the time-variant coefficients are periodic functions represented analytically by the Fourier series in table 2. Note that since *α*_{ref} and *U*_{ref} are constants, the only dependent variables in the equations are the four state variables *u*, *w*, *q* and *θ*, and the set of four equations is therefore complete. Equations (5.1)–(5.4) can immediately be made time-invariant by including only the zeroth harmonics of the periodic terms, whence they provide an NLTI model that is a proper mathematical generalization of the LTI model developed by Taylor & Thomas (2003).

### 5.2 Numerical solutions of the NLTI and NLTP models

The models we have developed use fundamental Newtonian mechanics equations (equations (2.1)–(2.4)) to predict how each locust would have flown, had it been released momentarily from its tether, with empirical measurements of the forces generated by the tethered locusts (figures 2 and 7; table 2; equation (4.10)). The solutions to equations (5.1)–(5.4) therefore simulate the free-flight trajectory that the locust would have followed for given initial conditions. This is an initial-value problem, in the sense that once the initial conditions are specified, equations (5.1)–(5.4) define a single, unique trajectory. Since equations (5.1)–(5.4) are nonlinear differential equations, it is necessary to solve them numerically. We did this using Matlab's ‘ode45’ solver, which is based on an explicit Runge–Kutta (4,5) formula—the Dormand–Prince pair (Dormand & Prince 1980). As a check on the numerical accuracy of the solver, we first solved the NLTI equations for each locust having set the initial conditions to be in quasi-static equilibrium (, , *q*_{0}=0, ) and confirmed that this produced a steady solution. A relative error tolerance of 10^{−4} and an absolute error tolerance of 10^{−7} were quite sufficient for this.

We then solved the NLTI equations for each locust with slightly perturbed initial conditions (; , *q*_{0}=0 and ). Figure 8 compares the resulting solutions by plotting the evolution of the four state variables through time for each individual. The solutions are only plotted for as long as the speed remains within the interval 2.0 m s^{−1}≤*U*≤5.5 m s^{−1} and the body angle remains within the interval 0°≤*α*≤14°. These intervals correspond to the range of flight conditions over which the forces were measured in the experiments by Taylor & Thomas (2003), and to which the empirical models of force production are fitted. This is, nevertheless, sufficient to indicate the instability of the solutions, which diverge monotonically in response to these perturbations over the time-intervals plotted. Note that the solutions take approximately 0.5 s to diverge from the range of speed and body angle for which the forces were measured in the experiments, which is equivalent to about 10 wing beat periods. The instability of the NLTI models is consistent with the instability of the earlier LTI models by Taylor & Thomas (2003). This is a major failing of the NLTI models, which we discuss fully in §6.1 below.

We next solved the full NLTP equations for each locust, including all eight harmonics. We set the initial values of the state variables to the conditions for quasi-static equilibrium (, , *q*_{0}=0, ), which, nevertheless, results in an unsteady solution because of the forcing terms. The solutions now depend also upon the phase of *ωt*, so we solved the equations for both *t*_{0}=0 and *t*_{0}=*π*/*ω* (where *π*/*ω* is the half-period of a single wing beat). Figure 9 plots the resulting solutions through time for as long as the speed remains within the interval 2.0 m s^{−1}≤*U*≤5.5 m s^{−1} and the body angle remains within the interval 0°≤*α*≤14°. Whereas the transients for the NLTI models in figure 8 show monotonic divergence, the combination of periodic forcing and periodic response characteristics in the NLTP models makes their transients much more complicated. This complicates their interpretation, and in order to reveal the true character of the models, it is therefore necessary to let the solutions run for longer. This in turn means relaxing the constraint of plotting the solutions for only as long as they remain within the range of flight speed and body angle measured in the experiments. This is done in figure 10, which extrapolates the models for body angles within the interval −5°≤*α*≤30°. This corresponds approximately to the range of body angles used by free-flying locusts (Fischer & Kutsch 2000) and so, while there can be no guarantee of the validity of extrapolating the empirical models of force production, the solutions are nevertheless confined to a biologically relevant range of speed and body angle (but note that pitch rate *q* can be of unrealistically large magnitude).

The solutions in figure 10 show the underlying dynamics of the NLTP models better than the shorter, initial sections of those solutions plotted in figure 9. Oscillations at the wing beat frequency are clearly visible in all of the solutions, corresponding to the effects of periodic forcing. This is superimposed upon a steady drift in *u*, which causes the systems to diverge indicating that the NLTP systems—like the NLTI systems of which they are generalizations—are unstable (§6.1 below provides a thorough discussion of this issue). Unfortunately, this divergence means that even the solutions in figure 10 do not run long enough for initial transient effects to die out, with the result that the solutions for locusts R and B in particular look rather erratic. However, the solutions for locust G take longer to diverge, leaving the system long enough to approach a steady-state oscillation (albeit one superimposed on a steady divergence in *u*).

In all three individuals, the forced oscillations in dorsoventral velocity *w* and pitch attitude *θ* are of similar phase to each other (figure 10), but lag pitch rate *q* by approximately *π*/2 radians and lag forward speed *u* by just over *π*/2 radians. This phasing of the motion components is very similar to the phasing of the short period mode found in the earlier LTI analysis by Taylor & Thomas (2003). Short period mode oscillations are not visible in the solutions to the NLTI models plotted in figure 8, but can be easily induced by disturbing the dorsoventral velocity component *w* from its equilibrium value. Figure 11 overlays a solution of the NLTI model for locust G, exhibiting short period mode dynamics (solid line) with the solution of the corresponding NLTP model (dotted line) for the same initial conditions (, , *q*_{0}=0, ). As in figure 10, the solutions are only plotted for as long as the speed and body angle remain within the ranges 2.0 m s^{−1}≤*U*≤5.5 m s^{−1} and −5°≤*α*≤30°. Whereas the oscillations in the NLTI case are excited only by the parameters of the system (like the oscillations of a pendulum), oscillations in the NLTP case result from a combination of parametric excitation and forcing. It is apparent that the short period mode oscillations have a slightly higher natural frequency than the externally imposed frequency of forcing in the NLTP model. We therefore suggest that periodic forcing excites short period mode dynamics, but that the short period mode frequency locks in on the forcing frequency, as is common in forced systems (Jordan & Smith 1999).

With appropriate choice of *t*_{0}, solutions of the NLTP model for locust G can remain within tight limits for over 2 s (figure 12), staying well within the range of speed and body angle used by free-flying locusts (2.0 m s^{−1}≤*U*≤5.5 m s^{−1}; −5°≤*α*≤30°). This dependence on *t*_{0} emphasizes the time-invariant character of the NLTP model. After the initial transients have died out, the system converges on an almost steady oscillation superimposed upon a slow drift in *u*. This oscillation is qualitatively similar to the other forced oscillatory solutions already discussed (figure 10), but is steady for longer, and therefore easier to interpret. After the initial transients have died out, the flight speed oscillates sinusoidally between about 2.80 and 2.95 m s^{−1}, while the body angle oscillates sinusoidally between about 4 and 22°. These values are well within the range of static speed and body angle used by free-flying locusts (Fischer & Kutsch 2000) but the amplitude of oscillation (8.5°) is larger than occurs naturally. The slow drift in *u* eventually causes the solution to diverge (after about 2.5 s), and as slight perturbations to the initial conditions also excite divergence (dotted lines in figure 12), the system as a whole is still unstable.

### 5.3 Graphical analysis of the NLTI and NLTP equations

Plotting the state variables against time (figures 8–12) is a direct way of showing how the solutions to the equations of motion evolve, but a major limitation of this approach is the need to plot so many different solutions; even after this exercise, there is no way of knowing whether all of the different qualitative behaviours of the models have been explored. However, since the solution of equations (5.1)–(5.4) for initial conditions (*u*_{0}, *w*_{0}, *q*_{0}, *θ*_{0}, *t*_{0}) comprises the set of four state variables *u*=*u*(*t*), *w*=*w*(*t*), *q*=*q*(*t*), *θ*=*θ*(*t*), time can be eliminated completely in the time-invariant case by plotting the state variables against each other. Moreover, because the right-hand side of the equations defines a unique tangent direction at every point in the state space {*u*, *w*, *q*, *θ*}, the paths that the system follows are unique and non-intersecting. This makes it possible to interpolate between paths, so that the qualitative behaviour of a time-invariant system can, in principle, be completely described by plotting a finite number of paths in state space (Jordan & Smith 1999). A collection of such curves is called a phase portrait because it provides a good indication of the qualitative character of the system.

Paths through state space can only terminate at points where the tangent direction is singular; paths that do not begin or end at a singular point must therefore either extend to infinity or form a closed curve. The most significant salient features of a phase portrait are therefore its singular points (which here correspond to fixed points of equilibrium), together with any closed paths that may be present. In simple cases, it is often possible to infer the qualitative behaviour of a system from its fixed points alone. We therefore include both fixed points of equilibrium in the phase portrait of each NLTI model (figure 13*a*–*c*); only by including both fixed points are we able to characterize the qualitative behaviour of the equations unambiguously. One of the two fixed points is always unrealistic for locust flight, but its inclusion should not be taken to imply that we are inferring anything about locust flight behaviour outside of the fitted range of the model. Furthermore, in discussing globally extrapolated behaviours, we explain what their physical interpretations would be in order to provide the physical grounding for the otherwise abstract mathematics, but these do not reflect on the flight behaviour of the locusts. The purpose of the state-space analysis is simply to characterize the equations fully, which is important because it allows us to answer a key question: are the dynamics of the NLTI models we have developed for the three locusts qualitatively the same?

The phase portraits of the NLTI models for locusts R and B (figure 13*a*,*c*) are structurally similar, in that it is possible to transform one topology into the other while preserving the qualitative arrangement of the paths (Wiggins 2003). This structural similarity implies that the dynamics of the two systems are qualitatively the same. Note that we could not reach this conclusion without plotting all of the fixed points because a necessary condition for structural similarity is that the nature and number of the fixed points is the same. In both cases there are two fixed points of equilibrium (denoted ‘F1’ and ‘F2’, respectively), each characterized as an unstable focus because of the system's tendency to spiral out through state space when disturbed slightly. Paths spiralling out from F1 follow a conical-helical trajectory aligned with the *q*-direction. This helical trajectory corresponds physically to an exponentially accelerating spin in the nose-up direction for paths above F1 (i.e. for positive *q*) and in the nose-down direction for paths below F1 (i.e. for negative *q*)—note that this allows us to predict the qualitative behaviour of solutions in the vicinity of the equilibrium in the fitted range of the model. Superimposed upon this broad helical pattern, however, is the influence of the second fixed point ‘F2’, which imparts a spiral character upon trajectories in its vicinity. This locally spiralling trajectory corresponds physically to an unstable pitch oscillation, with oscillations in *q* and *θ* growing exponentially until the system escapes the influence of F2 and diverges into an exponentially accelerated nose-down spin. Note, however, that these behaviours occur outside of the fitted range of the model.

The topology of the phase portrait of the NLTI model for locust G (figure 13*b*) is structurally distinct from that of the models for locusts R and B, and its dynamics therefore differ qualitatively from those described above. There are still two fixed points of equilibrium but, whereas F1 is still an unstable focus, F2 is now a stable focus. Paths leading out from F1 for *u*_{0}<*u*_{e} eventually enter a trajectory aligned with the *q*-axis, corresponding physically to an exponentially accelerating spin, which begins only after the trajectory has followed a divergent path curving asymptotically towards the *u*-direction. This steady divergence in *u* completely dominates the system dynamics in the vicinity of the equilibrium in the fitted range of the model. The transition from a path aligned with the *u*-direction to a path aligned with the *q*-direction is very sudden, and is characterized by an inflection marked ‘I’ in figure 13*b*. Paths leading out from F1 for *u*_{0}>*u*_{e} also display a steady divergence in *u*, but the trajectory curves asymptotically towards the *w*-direction and eventually enters F2. Although paths leaving F1 in either direction undergo monotonic divergence from the outset, paths passing close to F1 spiral in on these monotonically divergent trajectories, corresponding physically to a damped pitch oscillation superimposed on an overall divergence. These dynamics are apparent in the time plots of figure 11 (solid lines), in which short period mode pitch oscillations are superimposed upon an overall divergence.

Graphical analysis of differential equations is more complicated for time-variant systems like the NLTP models we have developed here. Although equation (5.1)–(5.4) still defines unique solutions for given initial conditions, it is possible for the resulting paths to intersect in state space because the right-hand side of the equations define a tangent direction to the curve at every point * p*=(

*u*,

*w*,

*q*,

*θ*), but involves the time-varying quantities

*X*=

*X*(

*t*),

*Z*=

*Z*(

*t*) and

*M*=

*M*(

*t*). Merely specifying the coordinates of

*is therefore insufficient to determine the tangent*

**p***at that point: the current time*

**t***t*must also be given to evaluate

*X*(

*t*),

*Z*(

*t*) and

*M*(

*t*). It follows that at two different time instants,

*t*

_{1}and

*t*

_{2}, we may have

*(*

**p***t*

_{1})=

*(*

**p***t*

_{2}) but

*(*

**t***t*

_{1})≠

*(*

**t***t*

_{2}) because, for example,

*X*(

*t*

_{1}) differs from

*X*(

*t*

_{2}). Hence, we would see two different curves passing through the same point

*. One partial solution to this problem is to fix*

**p***t*

_{0}, so that each point

*=(*

**p***u*,

*w*,

*q*,

*θ*) has one unique tangent , though this does not permit interpretation of the plot as a phase portrait. We have arbitrarily set

*t*

_{0}=0 in the models for locusts R and B, which is reasonable because the choice of

*t*

_{0}does not qualitatively affect their behaviour, but we have set

*t*

_{0}=−0.003 s in the model for locust B, which is the value of

*t*

_{0}giving the most persistent oscillations (see §5.2 above).

The topologies of the NLTP system plots differ from the NLTI phase portraits in one crucial respect: they contain no fixed points of equilibrium. Hence, whereas the NLTI systems need to be disturbed from equilibrium in order to set them on a trajectory through state space, the NLTP systems are always transitioning between states. In other respects, the trajectories plotted for the NLTP models for locusts R and B are broadly similar to those of the corresponding NLTI models (compare figure 13*a*,*d* to 13*c*,*f*, respectively), but this superficial similarity is another example of why it is important to plot any singular points that do exist in order to characterize the dynamics properly. Like the corresponding NLTI systems, the NLTP systems for locusts R and B have a tendency to enter a helical trajectory, corresponding physically to an exponentially accelerating spin. Periodic forcing introduces some new oscillatory character to the NLTP systems over the corresponding NLTI systems (this is visible in the initial spiralling of the trajectories in figure 13*c*,*f*; see also figure 11 for locust G), but essentially the same pitch instabilities dominate the system dynamics as in the NLTI case. This overall similarity in the trajectories plotted for the NLTI and NLTP models for locusts R and B makes the contrast between the trajectories plotted for the NLTI and NLTP models for locust G all the more marked.

Once again, the state space topology of the NLTP model for locust G (figure 13*b*) is structurally distinct from that of the models for locusts R and B, and its dynamics therefore differ qualitatively from those described above. Figure 13*e* plots a single trajectory for the NLTP model for locust G. The trajectory begins at the point corresponding to the position of F1 in the NLTI model and the system quickly converges upon an almost steady state oscillation in pitch (cf. figure 12), which in the state space representation is manifested in the orbital character of the trajectory. While we have not found a closed curve trajectory in this region of the state space, it is possible that one exists in the vicinity, in which case it would be an unstable orbit. Superimposed upon this almost steady pitch oscillation is a steady drift in *u*, which eventually causes the trajectory to diverge. Nevertheless, the drift is slow enough for the pitch oscillation to persist in an almost steady state for more than 2 s (see figure 12). However, once the trajectory diverges, it does so quickly, spiralling round to the right-hand side of the plot where the trajectory again takes on an orbital character. A true closed curve trajectory does appear to exist in this region of state space, characterized as an asymptotically stable orbit. This orbit falls well outside of the fitted range of the model, and corresponds physically to a fast inverted descent accompanied by a steady pitch oscillation. It is clear enough that the stable orbit in the NLTP model corresponds to the stable focus F2 in the NLTI model—formally, as the amplitude of forced oscillation tends to zero a stable orbit will tend to degenerate to a stable focus. We can therefore make a reasonably strong case that had the foci F1 been stable in the NLTI models, we might have expected the same region of state space in the NLTP models to be occupied by an asymptotically stable orbit. This is an intriguing possibility that we explore in more general terms in §6.

In summary, the techniques of graphical analysis that we have introduced in this section allow us to formalize and generalize our descriptions of the dynamics of the various models we have developed. This is only possible because the modelling framework that we use is nonlinear: although the empirical models of force production are local, the NLTI and NLTP equations are global, and the graphical analysis is therefore a mathematically valid technique for exploring the model dynamics. The graphical analysis allows us to show formally that the dynamics of the NLTI and NLTP models for locust G differ qualitatively from those of the models for locusts R and B, in spite of the broad similarity in their respective empirical periodic models of force production (figures 2–7). This shows that relatively small differences in the model parameters can have large qualitative effects on the models' dynamics, and further emphasizes the overriding importance of adequate sampling at the level of the individual (see §3). It follows that similarly small differences in the control responses of insects could have similarly pronounced effects on their flight dynamics. We therefore conclude that a proper dynamics analysis will always be essential to the correct interpretation of experimental results on insect flight control.

## 6. Discussion

### 6.1 Why do the models fail to explain locust flight stability fully?

Free-flying locusts are stable but, in common with the LTI models of locust flight dynamics developed by Taylor & Thomas (2003), all of the models developed here are unstable and therefore fail to explain locust flight stability fully. Taylor & Thomas (2003) identified five reasons why their LTI model could have failed, which we briefly repeat here: (i) the experimental data may have been unrepresentative of natural free flight; (ii) the model assumed that the forces and moments were functions of body angle with respect to the relative flow, but locusts may also be able to sense and respond to pitch attitude with respect to the horizontal; (iii) the model assumed that the forces and moments were not functions of pitch rate, which in fact they must be because of viscous damping; locusts may also be able to sense and respond to pitch rate; (iv) the model was fully linearized, but a linear framework may be inadequate to describe locust flight dynamics; and (v) the model was time-invariant, but most of the system parameters are strongly periodic and the system is characterized by periodic forcing. A key motivation for developing the NLTP framework was to eliminate these last two possibilities.

The nonlinear flight dynamics model remains linear periodic in the forces, but we have shown empirically that this is reasonable over the range of values of the state variables measured. We also continue to treat the centre of gravity and moment of inertia as time-invariant, which is a reasonable first approximation as far as the small oscillating mass of the wings is concerned (the wings comprise less than 4% of total body mass; Taylor & Thomas 2003). On the other hand, our measurements do not take account of possible changes in the centre of gravity resulting from abdominal movements (Taylor 2001). While these are almost certainly too slow to be important on the time-scale of the periodic flight forces, it is possible that systematic changes in abdominal position might affect the quasi-static balance of the forces (tethered locusts curl the last few abdominal segments upward when the imposed body angle is low). Nevertheless, with the introduction of nonlinear dynamics and periodic forcing, we are now confident that the models' instability does not result from deficiencies in the modelling framework. This leaves two distinct groups of possible explanations for the failure of the models: (i) failure to measure certain important parameters and include them in the model and (ii) issues related to how the experimental data were collected. We deal directly with the first group of explanations below. The second group of explanations pertains generally to the question of whether tethered flight data can be expected to produce reasonable models of free flight dynamics. This is such an important issue, and subject to so much misunderstanding, that we deal with it separately in §6.2 below.

As was also the case for the LTI model developed by Taylor & Thomas (2003), the NLTI and NLTP models developed here do not include the effects of pitch rate damping and pitch attitude damping. The neglect of pitch rate damping is especially significant, because it allows pitch rate to reach unrealistically high levels in our simulations. The solutions we have plotted are restricted by placing specified limits on speed and body angle, but we have placed no limit on pitch rate. Since pitch rate was zero in the experiments in which the empirical model of force production was derived, this rapid growth in pitch rate clearly places further limits on the validity of the simulations. Neither pitch rate damping nor pitch attitude damping can be measured using a conventional wind tunnel setup, and it will not be possible to develop an improved model of locust flight control without developing new experimental techniques and apparatus. This would need to be capable of measuring the forces and moments on a tethered locust during pitching (in order to measure pitch rate damping), and would need to be capable of rotating the airflow independently of the insect (in order to separate the effects of pitch attitude and body angle). Work is now underway to develop such a system, similar to the rotary test rigs used by aircraft engineers to analyse spin dynamics and unsteady aerodynamic effects (Ericsson 1990). This presents a considerable technical challenge but illustrates a key general point: that placing our understanding of insect flight control on a flight dynamics footing demands that future experimental work is driven by explicitly specifying the parameters that need to be measured in order to formulate a proper physical model.

### 6.2 Can tethered flight measurements be used for models of free flight dynamics?

Should we really expect force measurements from tethered locusts to produce a workable model of their free flight dynamics? The answer to this question requires a proper understanding of the issues raised by tethering. The objection most commonly raised to tethered flight experiments is that they are made under ‘open-loop’ conditions. Although this terminology is firmly entrenched in the biological literature (including the earlier paper by Taylor & Thomas 2003), it is misleading from the perspective of control theory. We therefore refer hereon to tethered flight measurements being made under conditions of ‘broken dynamics’, to emphasize that the feedback loops of the control system remain physiologically closed and intact, while the dynamics of the system are broken physically by the kinematic constraint of tethering. This is a crucial distinction, because it is therefore possible—in principle—to expose a tethered insect to stimuli identical to those experienced during free flight. To clarify this further, we will discuss three possible classes of tethered flight experiment. Only the first class of experiment has as yet been performed. We discuss the other two classes as illustrative thought experiments, but new experimental apparatus may soon allow us to implement functionally at least the second class.

In the first class of tethered flight experiment—of which the experiments described here and in Taylor & Thomas (2003) are an example—the insect is exposed to systematic changes in flight condition that are prescribed in advance of the experiment. It is assumed that the instantaneous values of the state variables are compatible with free flight, whether steady or accelerated, and it is further assumed that the visual stimuli can be made indistinguishable from those experienced naturally (this is not expected to be a problem for the other sensory modalities). There are then three ways in which the insect could sense the effects of tethering. Firstly, the insect could sense the tether directly via mechanoreceptors in its cuticle. Secondly, the insect could detect the broken dynamics by comparing the time-course of the stimuli it receives with the expected changes in stimuli given the control actions it effects (note that this implies that the insect possesses some internal representation of its flight dynamics). Thirdly, although the instantaneous values of the state variables may be compatible with free flight, their time-course may not be, and this may also be detectable by the insect if it is able to integrate the stimuli it receives through time. For example, if the insect tends to oscillate slightly in free flight (as may be suggested by the oscillatory trajectories of the NLTP models), then tethering the insect statically, as we have done here, will result in an unrealistic stimulus time-series. This particular issue can be controlled for directly by oscillating the insect with the appropriate phase, amplitude and frequency. This would also be important in correcting for the influence of unsteady aerodynamic effects, given that accelerations of the body will lead to some degree of aerodynamic unsteadiness over and above that caused by the unsteady flapping of the wings, which might be significant for aerodynamic force production (Ramamurti *et al*. 2002).

In the second (currently hypothetical) class of tethered flight experiment, the insect would be exposed to stimuli that are again prescribed in advance of the experiment. This time, however, the motions of the insect would be recorded in free flight, and the insect would then be moved along exactly the same trajectory and through the same visual environment when tethered. For this kind of experiment to be meaningful, it is important that a sufficiently strong stimulus is used that we can reasonably expect the insect to respond to in a repeatable fashion. For example, the insect could be blown sharply to elicit a corrective response, presented with a looming object to elicit an avoidance response or presented with a moving target to elicit a chasing response. On the assumption that the insect responds in a more or less repeatable fashion (which will be violated if the insect has any tendency to engage in spontaneous manoeuvres), this experiment ensures that the stimulus time-series is realistic and at least approximates the expected changes in stimuli given the control actions of the insect. It does not eliminate the possibility that the insect can feel the tether directly (except in the unlikely situation that the insect responds in an identical way to that recorded in free flight, so that the reaction force of the tether remains zero throughout).

In the third (also currently hypothetical) class of experiment, the stimuli to which the insect would be exposed would not be prescribed in advance of the experiment. Instead, the forces that the insect generates would be measured, and equations of motion used to calculate the instantaneous trajectory through which the insect should be moved. Provided that the time lags in the apparatus can be kept short enough that they are undetectable by the insect, the stimuli that the insect receives should then be indistinguishable from those in free flight. This is exactly the principle underlying the mechanical flight simulators used to train pilots to fly aircraft. Note, once again, that although the insect will move exactly as it would have done in free flight, the kinematic constraint of tethering breaks the direct physical link between forces and kinematics. This simple thought experiment illustrates very clearly that the effect of tethering is to break the insect's own flight dynamics, not to break its control loops. The same is true of the other two classes of tethered flight experiment, and it should therefore be apparent that if the instantaneous stimuli that the tethered insect receives are compatible with free flight, the only issue that must concern us is whether the insect responds in the same way to these stimuli as it would if it experienced them in free flight.

This is certainly an important issue, but it is one that we have no choice but to accept. There is no avoiding the need to formulate flight dynamics models using force measurements, which in turn necessitates the kinematic constraint of tethering. Measurements of free-flight kinematics cannot offer a robust substitute for tethered force measurements because the instantaneous acceleration needed to calculate the forces can only be inferred from direct measurements of position, and perhaps velocity, but as we have already shown, rapid changes in the forces may have a negligible (i.e. practically incalculable) effect on the flight trajectory. For example, the bodies of very small insects do not oscillate visibly as they beat their wings, yet the forces acting on the body are undoubtedly periodic. The only obvious resolution to this problem is to use free flight data to validate the predictions of models based on tethered flight data using a system identification approach (e.g. Klein 1989, for a review of aircraft parameter estimation from flight data).

The final question we must answer is how measurements from a system with broken dynamics can be used to model the corresponding dynamic system. The answer is simple, and forms the basis of the analysis we have presented here: there is no need to measure the system with its dynamics intact, provided that all of the relevant parameters can be measured, because the broken link between forces and kinematics is intact in the equations of motions in which the measured parameters are embedded. This is why we have called our models ‘semi-empirical’: the parameters are measured empirically in a system with broken dynamics, but the dynamical link is specified a priori by Newtonian mechanics. Tethered flight measurements are not merely suitable for formulating models of free flight dynamics—they are essential, but only when combined with equations of motion do they say anything meaningful about dynamics.

### 6.3 What do the models tell us about locust flight control?

Given that the models fail to explain locust flight stability fully, can they tell us anything useful about locust flight behaviour? The question of whether a system is stable or not is just one of many questions asked in dynamics. Indeed, it is worth pointing out at the outset, that if the models had actually resulted in being stable, then this in itself would have told us nothing we did not already know about locust flight! The more interesting questions in flight dynamics relate to the response characteristics of the system—in particular its inertial response to inputs of different frequencies. While care must be taken in interpreting a model that fails to explain the stability of the system it describes, it is nevertheless possible to draw some robust conclusions about the response of the systems to inputs of different frequencies. The response of any system is limited by its damping and inertia. We are confident that we have identified the inertial characteristics of the locusts correctly because these are mass distribution properties that can be simply and accurately measured, but we do not claim that we have measured all of the damping terms (see §6.1), and the models are therefore expected to underestimate the system damping. This means that the inertial response characteristics of the NLTP models can be reliably used to provide an upper bracket on the input frequencies to which locusts should be able to respond. Once the initial transients have decayed away, the forced oscillations of the NLTP models lack higher harmonic content (note the sinusoidal form of the forced oscillations in figure 12, and compare this with the forcing functions in figures 2–7). The mechanical system is therefore acting like a low-pass filter with a cut-off frequency below the frequency of the second harmonic.

The cut-off frequency of the inertial response of locusts is expected to be as low or lower than the frequency of the second harmonic when the missing damping terms are included. This implies that subtle changes to the higher harmonics of the periodic flight forces will not translate into subtle changes in the locust's motion, because the system will effectively filter out such effects. For example, increasing the lift generated at the two points of stroke reversal would be expected to have the same qualitative effect as increasing the lift evenly throughout the down stroke. Similarly, although the relative timing of the fore- and hind-wings may be important to the unsteady aerodynamic forces they generate, the timing of the forces resulting from their interaction is likely to be unimportant from a flight dynamics perspective. This is reassuring; it is impressive enough that four-winged insects manage to control the aerodynamic forces they generate by adjusting wing phase (Taylor 2001), but it is difficult to imagine how they would be able to use this as a workable system of flight control if the timing of the resulting forces relative to the motions of the whole insect also mattered! Manoeuvre control and flight stabilization in locusts must therefore operate on a time-scale that is relatively coarse in comparison with the period of the wing beat. This is, of course, reasonable for an insect with no great need for agility.

### 6.4 A limit-cycle analogy for flapping flight stability and control

Although the graphical analysis of the NLTP models was aimed at formally describing the dynamics of the systems we have empirically defined, the state space approach to analysing flight control leads us to propose a new hypothesis for flapping flight control in general. For maximum generality, we will reason this hypothesis from first principles. By definition, the non-manoeuvring flight of insects must in some sense be steady. Given that periodically forced systems, in general, have no fixed points of equilibrium, the only steady state towards which flapping flight can evolve is an asymptotically periodic oscillation (Jordan & Smith 1999), forming a closed curve when plotted in state space. Since flight must also be stable, it follows that the proper definition of stability in flapping flight is that of asymptotic orbital stability (Hahn 1967), which means, loosely speaking, that the system should tend to return to this orbit when disturbed slightly from it. Note that if the amplitude of body oscillation is negligibly small (as in very small insects), we may approximate the stability of the system by asymptotic stability with respect to a fixed point. The dynamics of the NLTP model for locust G provide a useful illustration of these general principles: with appropriate initial conditions, the system initially tends towards a pitch oscillation that is almost periodic (and also compatible with the range of speed and body angle used by free-flying locusts), but which eventually diverges because of its instability and asymptotes towards a second oscillation that is truly periodic and asymptotically orbitally stable (although incompatible with the range of speed and body angle used by free-flying locusts).

If the proper definition of stability in flapping flight is that of asymptotic orbital stability, then we may hypothesize that the problem of actively stabilizing and controlling flapping flight amounts to stabilizing the orbit with respect to small disturbances, and moving from one orbit to another in order to effect manoeuvres over the course of more than one wing beat. This hypothesis is reminiscent of a concept commonly referred to as ‘limit-cycle control’ in the literature on bipedal walking, although the term is used loosely because most powered walking systems are forced (i.e. non-autonomous) systems and limit cycles are only strictly defined for autonomous systems. The limit-cycle control approach (Markus 1973; Somasundaran & Balachandran 1986) has been used to derive control laws for stabilizing simulated bipedal walking machines (Katoh & Mori 1984; Hmam & Lawrence 1992; Laszlo *et al*. 1996; van de Panne 2000) and has been implemented physically to stabilize the forced oscillations of an inverted pendulum (R. Q. Van der Linde, unpublished study). In principle, limit-cycle control can either enlarge the region of stability around an already stable limit cycle (Hmam & Lawrence 1992), or stabilize intrinsically unstable limit cycles (Laszlo *et al*. 1996; van de Panne 2000). Either way, the key is to ensure that each cycle is stabilized to such a degree that the next begins in an appropriate initial state. Provided that any inherent instability in the system would normally take several cycles to develop, the addition of even quite moderate levels of feedback control may then suffice to stabilize a single limit cycle to the extent that the next can also be stably maintained. For example, given that the unstable pitch oscillations in figure 12 already persist for several seconds, it is likely that adding moderate feedback to counteract the steady divergence in *u* would be sufficient to make the system converge upon a true orbital trajectory.

The limit-cycle control concept is problematic for non-autonomous (i.e. time-variant) systems because the system's trajectory will, in general, depend upon the value of *t*_{0}, and this may mean that an asymptotically stable orbit is not reached for all values of *t*_{0}. We therefore conjecture further that flapping flight cannot be stable with respect to a steady-state oscillation unless the system is asymptotically autonomous (i.e. asymptotically time-invariant), *sensu* Markus (1956), meaning that its trajectories converge asymptotically towards the same orbit irrespective of *t*_{0}. This offers a formal justification for the analogy with stable limit cycles, which are the stable orbits in the autonomous system towards which the non-autonomous system asymptotes. Manoeuvre control would then consist in moving the system from one asymptotically autonomous orbit to another, via non-autonomous transients which would need to be correctly timed in order to allow the system to transition from one orbit to another (Markus 1973). An illustrative example of this kind of transient behaviour can be seen in figure 13*e*. This shows the system initially entering an almost steady oscillation in pitch (line colour *θ* oscillating through varying shades of red), which becomes unsteady and evolves through a transient (line colour *θ* changing monotonically from red to blue) into a new steady-state oscillation (line colour *θ* oscillating through varying shades of blue). We therefore propose that the most general definition of stability in flapping flight is that of asymptotically autonomous orbital stability. Manoeuvre control on the scale of more than one wing beat would then consist in exciting non-autonomous transients from one asymptotically autonomous and stable orbit to another.

### 6.5 Future directions

Nonlinear time-periodic modelling offers the most physically realistic framework available for analysing flapping flight dynamics. The semi-empirical NLTP models that we develop here are therefore a major advance on the LTI models developed previously by Taylor & Thomas (2003). While there is still scope for improvement—for example, by allowing for the possibility that wing beat frequency itself may also vary systematically with the state variables, or by allowing for periodic variations in the centre of gravity and moment of inertia—we believe that the modelling framework for analysing experimental data on insect flight control is now mature. Where, then, should the flight dynamics paradigm take us next? In unambiguously specifying which parameters are missing and need to be measured to produce a fully realistic model, the modelling framework has leapt a step ahead of experiments on insect flight control (see §§6.1 and 6.2): the modelling framework drives the experimental data we collect. This is an excellent example of the reciprocal relationship between experiment and analysis with experiment driving analysis and analysis driving experiment—it is also a good illustration of the practical value of the flight dynamics approach.

Perhaps more importantly, in inspiring the limit-cycle control analogy, the state space method of analysing our models of insect flight dynamics poses new questions for biologists to consider. Do insects use limit-cycle control to stabilize their flight? This is a question that cannot easily be answered without first understanding the control laws written into the networks of sensors, actuators and neurons that constitute the insect's sensori-motor system. If insects do indeed use limit-cycle control of flight, then we should expect these control laws to display an underlying periodicity matched to the periodicity of the flight dynamics. We should also expect to observe switches in the underlying control laws during transients from one asymptotically stable orbit to another (§6.4). Future work will therefore need to unite external models of flight dynamics with internal models of sensory and neuromuscular physiology, by combining models of manoeuvre dynamics with recordings of individual muscles and sensory neurons responding to stimuli eliciting the manoeuvre. Only when we link the physics and physiology of flight directly in this way can we claim an integrated systems understanding of insect flight control.

## Acknowledgments

We thank Johnny Evers, Marc Jacobs, Kent Kaiser, Belinda King, Valerie Martindale, Adrian Thomas and Ric Wehling for helpful discussions and for their support of this work. This research programme is supported by the Air Force Office of Scientific Research International Research Initiative (IRI contract F61775-02-C4048). The analysis was carried out during a travelling fellowship by G.K.T. and visiting professorship by R.Ż. to the Graduate Engineering and Research Center, University of Florida, made possible by grants from the Air Force Office of Scientific Research and the European Office of Aerospace Research and Development. We thank the Graduate Engineering and Research Center and the University of Florida for their hospitality during our stay. The experiments on which the analysis is based were supported by a Christopher Welch Scholarship at the University of Oxford, and by BBSRC grants 43/S09380 and 43/S08664. G.K.T. is a Royal Society University Research Fellow, and was supported by a Weir Junior Research fellowship at University College, Oxford, and by a Royal Commission for the Exhibition of 1851 Research Fellowship for the duration of this work.

## Footnotes

The supplementary Electronic Appendix is available at http://dx.doi.org/10.1098/rsif.2005.0036 or via http://www.journals.royalsoc.ac.uk.

- Received November 29, 2004.
- Accepted February 24, 2005.

- © 2005 The Royal Society