The dynamic, three-dimensional shape of flapping insect wings may influence many aspects of flight performance. Insect wing deformations during flight are largely passive, and are controlled primarily by the architecture and material properties of the wing. Although many details of wing structure are well understood, the distribution of flexural stiffness in insect wings and its effects on wing bending are unknown. In this study, we developed a method of estimating spatial variation in flexural stiffness in both the spanwise and chordwise direction of insect wings. We measured displacement along the wing in response to a point force, and modeled flexural stiffness variation as a simple mathematical function capable of approximating this measured displacement. We used this method to estimate flexural stiffness variation in the hawkmoth Manduca sexta, and the dragonfly Aeshna multicolor. In both species, flexural stiffness declines sharply from the wing base to the tip, and from the leading edge to the trailing edge; this variation can be approximated by an exponential decline. The wings of M. sexta also display dorsal/ventral asymmetry in flexural stiffness and significant differences between males and females. Finite element models based on M. sexta forewings demonstrate that the measured spatial variation in flexural stiffness preserves rigidity in proximal regions of the wing,while transferring bending to the edges, where aerodynamic force production is most sensitive to subtle changes in shape.

During flapping flight, insect wings accelerate masses of air, generating the forces necessary to support the insect's weight and to perform complex maneuvers. At the same time, these ultralight airfoils (generally only 0.5-5%of body mass; Ellington,1984b) must withstand the forces imposed upon them by the surrounding air, as well as the inertial forces caused by accelerating and decelerating their own mass up to several hundred times per second.

Insect wings perform these roles extremely successfully, despite the fact that they are largely passive structures, with no muscular control past the wing base (Wootton, 1992). Although they are strengthened by a network of tubular veins, the wings of many species deform noticeably during flight, especially during slow flight and hovering (Willmott and Ellington,1997a). These dynamic changes in the three-dimensional shape of wings could potentially affect many aspects of force production, yet few models of insect flight have successfully incorporated passive wing flexibility or examined the effects of flexibility on force production.

In the past two decades, models of insect flight, both mathematical and physical, have contributed enormously to our understanding of the basic mechanisms of force production, despite assuming that insect wings are rigid structures [e.g. mathematical models of Ellington(1984a), Lighthill(1973), Sane and Dickinson(2002), Savage et al. (1978),Smith et al. (1996) and Wilkin and Williams (1993), and physical models of Bennett(1966, 1970), Dickinson et al.(1999), Ellington et al.(1996) and Spedding and Maxworthy (1985)]. Some aspects of wing flexibility have been mimicked in models by altering the relative positions of wing regions(Liu et al., 1998; Vest and Katz, 1996) or by modeling deformations as harmonic waves(Combes and Daniel, 2001; Daniel, 1987; Wu, 1971). These approaches provide unique insights into the mechanisms of force generation during flight,but often neglect one or more critical components of wing deflection (e.g. spanwise bending, chordwise bending or torsion), which can have large effects on aerodynamic force production (Batchelor,1967). Models of insect flight that incorporate passive wing flexibility (in which shape changes are driven by forces imposed upon the wing rather than being specified in advance) are exceedingly rare (e.g. Smith, 1996).

One difficulty in modeling passive wing flexibility is that forces applied at the wing base lead to bending and twisting that are influenced not only by overall stiffness and gross anatomical features (e.g. flexion or fold lines),but also by the distribution of flexibility throughout the wing(Wootton, 1999). Extremely detailed finite element models of insect wings (incorporating vein configuration, three-dimensional relief, and variations in vein and membrane characteristics) could potentially reproduce the distribution of flexibility in real wings by accounting for the precise structural and material properties of the wing. Unfortunately, precise data about the local properties of insect wings are often unavailable and finite element models must be constructed using simplifying assumptions. For example, many finite element models assume that the material stiffness of the wing (Young's modulus, E) is equivalent to that of insect chitin and constant throughout the wing(Smith, 1996; Kesel et al., 1998); however,recent measurements reveal that E can vary widely within a wing(Smith et al., 2000) and that other proteins, such as resilin, occupy key positions (e.g. wing vein joints)in insect wings (Gorb, 1999; Haas et al.,2000a,b). Measuring spatial variation in these material properties, as well as the details of vein and membrane structure, is a time-consuming process that would need to be repeated for each new species studied.

Although these detailed approaches can provide important information about functional wing morphology in pivotal, well-studied species (e.g. Herbert et al., 2000), a more general approach to wing stiffness measurements and modeling could facilitate comparative studies and attempts to incorporate passive flexibility into models of force production. Rather than measuring (and modeling) geometric and material properties separately, the overall bending response of the wing, or flexural stiffness (EI, the product of material stiffness Eand second moment of area I) can be determined. Performing measurements of EI averaged over the whole wing is relatively straightforward (see Combes and Daniel,2003a), but measuring spatial variation in flexural stiffness throughout a wing is more challenging. Steppan(2000) approached this issue by measuring average flexural stiffness over increasingly larger sections of dried butterfly wings, and Wootton et al.(2000) measured average flexural stiffness in three isolated sections of locust hindwing. However, no studies to date have demonstrated a method of measuring the spatial variation in local flexural stiffness of intact insect wings.

In this study, we developed a method of approximating spatial variation in flexural stiffness along two axes of the wing (in the spanwise and chordwise direction). We measured flexural stiffness variation in two insects, the hawkmoth Manduca sexta and the dragonfly Aeshna multicolor. These insects have wings of similar size and are both agile fliers, capable of hovering as well as fast, forward flight. However, they are very distantly related and display large differences in wing shape and venation pattern (see Combes and Daniel, 2003a) that may underlie differences in the distribution of flexural stiffness in their wings.

To determine the spatial pattern of flexural stiffness in wings, we first developed a method to measure displacement (due to a point force) along the wing in the spanwise and chordwise directions. Next, we proposed various alternatives for how flexural stiffness might vary along the wing (represented by simple mathematical functions) and predicted the patterns of displacement that a loaded wing with these theoretical stiffness distributions would display. We then found the flexural stiffness distribution that produced a pattern of displacement most similar to the pattern measured in a real wing. We used this method to estimate the flexural stiffness distribution of the wing in response to forces applied from both the dorsal and the ventral side,as dorsal/ventral asymmetry has been noted in previous studies(Ennos, 1988; Steppan, 2000; Wootton, 1993; Wootton et al., 2000).

We also created two different finite element models of an insect wing,based on the planform geometry of a Manduca sexta forewing. Rather than realistically reproducing the structure and behavior of Manducawings, these generalized models provide a method of assessing the consequences of flexural stiffness distributions to wing bending. We attempted to capture major elements of wing structure in the models (e.g. vein position and flexural stiffness distribution) while avoiding detailed structural features,such as variation in wing thickness and three-dimensional relief. We used these models to examine how different spatial patterns of EIdetermine bending patterns in wings subjected to both static and dynamic loads.

### Flexural stiffness measurements

Insect handling and wing preparation were performed as in Combes and Daniel(2003a). Calculations of EI distribution require both the applied force and the displacement along the length of the wing. We measured displacement of wings in the spanwise direction along a line running from the wing base to the tip, and in the chordwise direction along a line running from the leading edge to the trailing edge (approximately midway between the base and tip). We attached a fresh wing at either the wing base or leading edge, and applied a point force to the wing with a pin attached to a force beam (see Combes and Daniel, 2003a),using a small drop of cyanoacrylate glue to prevent the pin from slipping off the wing tip or trailing edge. We illuminated the wing with sheets of laser light; these sheets appear as lines on the wing, and can be adjusted so that one line illuminates a strip of the wing from the base to the tip(Fig. 1A) or from the leading to the trailing edge (intersecting the pin). When a force is applied to the wing, the laser lines follow the surface of the wing, shifting their position in proportion to the displacement of the wing(Fig. 1B). We used a Coolpix 900 digital camera (Nikon, Tokyo, Japan) angled approximately 45° from the surface of the wing to photograph unloaded and loaded wings, as well as a calibration object containing a series of steps of known height.

Fig. 1.

Unloaded (A) and loaded (B) Manduca sexta wing illuminated with laser lines. A Matlab program was used to find the center of a laser line running from the base of the wing to the tip before and after applying a point force at the tip (yellow and red lines, respectively). The change in position of the laser line y(x) was then used to find displacement of the wing δ(x).

Fig. 1.

Unloaded (A) and loaded (B) Manduca sexta wing illuminated with laser lines. A Matlab program was used to find the center of a laser line running from the base of the wing to the tip before and after applying a point force at the tip (yellow and red lines, respectively). The change in position of the laser line y(x) was then used to find displacement of the wing δ(x).

We measured displacement in response to loads of varying magnitude,returning to the unloaded position to photograph the wing before applying each new load. We checked the repeatability of displacement measurements on the wings of a hawkmoth (in both the spanwise and chordwise direction) by applying five different loads to the dorsal surface of the wing and repeating each load three times (applying loads in random order).

For all other hawkmoth and dragonfly wings, we measured displacement on the dorsal surface of the wing at four different loads, gently removed the glued pin from the edge of the wing (without tearing the wing membrane), flipped the slide over, and measured displacement on the ventral side of the wing at four different loads. Dorsal stiffness in this study describes the stiffness of the wing in response to loading on the dorsal surface (resulting in a dorsally convex surface; see Steppan,2000), which is equivalent to ventral flexion in other studies(Ennos and Wootton, 1989; Wootton, 1981, 1992).

After completing the flexural stiffness measurements, we photographed the wings on a white background from above and measured wing area, span and maximum chord length in NIH Image. For Manduca sexta (whose wing morphology may depend on gender; Willmott and Ellington, 1997b), we also calculated aspect ratio(span2/area) of the forewings and wing loading (body mass divided by the area of both forewings).

In each hawkmoth, we measured spanwise flexural stiffness on one forewing and chordwise flexural stiffness on the other. We did not measure the flexural stiffness of the hindwings, which are much smaller and overlap with the forewings to varying degrees during flight. In dragonflies (which have independent fore- and hindwings of similar size), we measured spanwise flexural stiffness on one forewing and one hindwing, and chordwise flexural stiffness on the other forewing and hindwing.

We analyzed images of loaded and unloaded wings with a custom Matlab program (developed by A. Trimble) that finds the center of the laser line running from the wing base to the tip or from the leading to the trailing edge(Fig. 1A,B). We used a 2nd order Butterworth filter to remove noise from the line position data, and splined the unloaded and loaded data to an equal number of points for comparison. We found the difference in line position between the two data sets [y(x); Fig. 1B] and converted this difference to actual wing displacement along the wing δ(x) with a factor derived from analysis of the calibration object.

With this displacement data and the applied force, we were able to estimate local flexural stiffness by using a continuous beam equation to approximate EI variation along the wing (see Appendix). This equation defines local flexural stiffness as a function of the local curvature (second spatial derivative of displacement) and the local moment (applied force times distance from the point of force application). To avoid errors caused by differentiating displacement data, we solved this as an inverse problem; we posed several simple mathematical functions that might approximate flexural stiffness variation along the wing (constant, linear, exponential, or 2nd degree polynomial) and calculated the expected pattern of displacement due to the applied force and these possible EIdistributions. We then used a simplex minimization program to find the EI distribution that provided the best fit to the measured displacement along the wing. Finally, we calculated average flexural stiffness(

$$EI$$
⁠) by integrating the equation describing flexural stiffness variation along the wing's length (see Appendix).

To verify the method of determining flexural stiffness distribution, we measured δ(x) and calculated local flexural stiffness of a rectangular glass coverslip, which is made of homogeneous material and has a constant material stiffness (E) and second moment of area(I).

### Comparisons within and between individuals

Because wing displacement in hawkmoth and dragonfly wings does not increase linearly with applied force (see Combes,2002), our estimates of flexural stiffness varied with the relative displacement of the wing (δT/L, tip displacement divided by wing span or trailing edge displacement divided by chord length). To facilitate comparisons between individuals and species, we standardized measurements at a relative displacement of 0.05 for spanwise measurements and 0.08 for chordwise measurements (values that are within the range measured on real wings). For each individual, we used the slope of the relationship between EI variables (

$$EI$$
and coefficients describing the distribution of flexural stiffness) and relative displacement to estimate the value of these variables at a givenδ T/L.

In Manduca sexta, we collected spanwise data from nine males and four females, and chordwise data from ten males and nine females. In Aeshna multicolor, we used only males, but tested both fore and hindwings (which have different morphologies), so we examined stiffness separately in the two sets of wings. We collected spanwise and chordwise data on both fore and hindwings from eight individuals.

Within each group, we tested for differences between the dorsal and ventral side of the wing with a Wilcoxon signed rank test. We also tested for differences between male and female Manduca sexta with a Mann-Whitney U-test.

### Finite element modeling

To investigate how flexural stiffness variation affects wing bending, we created two simplified finite element models of an insect wing (based on the forewing of a male Manduca sexta) with the same geometry, but dramatically different patterns of spatial variation in flexural stiffness(Fig. 2). We used MSC Marc/Mentat 2001 to construct model wings with an accurate planform configuration of veins and membrane, but a simplified three-dimensional geometry (with no camber or three-dimensional relief; see Discussion). We increased the material stiffness of vein elements beyond that of the surrounding membrane to mimic the higher second moment of area of tubular veins and produce the measured spanwise-chordwise anisotropy in flexural stiffness (see Combes and Daniel,2003a).

Fig. 2.

Finite element models based on Manduca sexta wings. (A) Model wing with membrane elements (blue) and vein elements (pink), each of homogeneous material stiffness E. (B) Model wing in which declining material stiffness of membranes and veins results in an exponential decline in flexural stiffness. Each color represents a different value of material stiffness.

Fig. 2.

Finite element models based on Manduca sexta wings. (A) Model wing with membrane elements (blue) and vein elements (pink), each of homogeneous material stiffness E. (B) Model wing in which declining material stiffness of membranes and veins results in an exponential decline in flexural stiffness. Each color represents a different value of material stiffness.

The models were composed of thin shell elements with a density of 1200 kg m-3 (as measured in insect wings; Wainwright et al., 1982) and a thickness of 45 μm. We used a Poisson's ratio of 0.49, as measured in some biological materials (Wainwright et al.,1982); because the Poisson's ratio of insect wings is unknown, we tested the effects of using a Poisson's ratio of 0.3 and found that the difference in model behavior was negligible. The uniform mass distribution of the models may cause dynamic bending in the distal regions to be overestimated, but qualitative differences in wing bending between the two models should not be affected. To determine the minimum number of elements necessary to capture the bending behavior of wings, we performed a sensitivity analysis with models composed of 200, 350, 865 and 2300 total elements, and found that 865 elements were sufficient to ensure asymptotic performance of the model.

We adjusted the material stiffness E of vein and membrane elements in the two models to produce different spatial patterns of flexural stiffness,but the same overall bending performance (so that tip and trailing edge displacement in response to a point load were the same as displacements measured in real wings). In the first model, we assigned all vein elements a single (homogeneous) Young's modulus of E=1.5×108 N m-2 (similar to values measured in locust hindwing; Smith et al., 2000), and all membrane elements a Young's modulus of E=2.1×1012 N m-2 (representing both increased material stiffness and the increased second moment of area of tubular veins; Fig. 2A). In the second model,we applied declining values of material stiffness to the model wing in 12 strips oriented diagonally (Fig. 2B); these strips are perpendicular to most of the wing veins,which decrease in diameter towards the wing edge and thus are likely to decrease in stiffness along this axis. We adjusted the values of material stiffness in these strips to approximate patterns of overall wing flexural stiffness measured in real Manduca wings (E in the model varies from 4.7×107 N m-2 to 4.5×109 N m-2 in membrane elements, and from 1.9×1011 N m-2 to 1.8×1013 N m-2 in vein elements).

To determine the resulting pattern of flexural stiffness variation in the spanwise direction of the model wings, we fixed each wing at its base (with no displacement or rotation) and applied a point force at the tip. For chordwise measurements, we fixed the model at its leading edge (from the base to 2/3 span) and applied a point force at the trailing edge. We recorded the displacement at 22 nodes (junctions between elements) aligned between the point of attachment and the point of force application, and used this information to estimate spanwise and chordwise flexural stiffness distribution with the Matlab simplex minimization program, as in real wings.

We compared static bending performance of the two models by fixing each wing at its base and applying either a point load of 0.003 N at the tip(within the range of loads applied to real wings) or a pressure load of -14.43 Pa to the wing surface (equivalent to the average lift force that the wing would experience during steady flight). We compared dynamic bending by flapping each model wing at a realistic wing beat frequency and stroke amplitude. We applied boundary conditions to the nodes at the wing base so that they could not translate in any direction and could rotate only in the dorsal-ventral direction (around the y-axis, see Fig. 2). We began the simulation with initial conditions of zero displacement and zero velocity at all nodes, and gradually increased the rotation at the wing hinge to a sinusoidal motion with the following function:
$\ {\theta}(t)=(1-\mathrm{exp}^{-\mathrm{t}{/}{\tau}}){\cdot}\mathrm{sin}({\omega}t),$
1
where θ is rotation at the base nodes, t is time in seconds,τ is the time constant and ω is the angular frequency(2πf, where f is the flapping frequency). We flapped the wings at 26 Hz, and found that a time constant of 1/20 s avoids transient artefacts of rapid initial acceleration and allows the wings to reach their full stroke amplitude of 108° after 6.5 flaps. We ran each model for 19200 time steps, simulating 12.5 flaps in 0.48 s.

Because the finite element program does not calculate aerodynamic forces acting on the wing, a damping factor had to be applied to stabilize the models. We compared the motions of the model wing with declining material stiffness to those of a real Manduca wing attached to a motor and rotated in the same way (Combes and Daniel,2003b; Daniel and Combes,2002). We found that a mass damping factor of 10 reduced high frequency vibrations and provided the closest match to motions of the real wing, and applied this damping factor to both model wings.

### Flexural stiffness measurements

The simplex minimization analyses revealed that proposed EIdistributions based on either exponential or 2nd-degree polynomial equations can provide relatively good predictions of measured wing displacement, while constant or linear distributions of EI cannot(see Combes, 2002). However,the proposed exponential EI distributions were more consistent within and between individuals, and effectively predicted displacement in the spanwise and chordwise directions of both Manduca sexta and Aeshna multicolor. Because these more consistent results allow us to make comparisons within and between individuals, we used an exponential equation to approximate flexural stiffness patterns in the species tested.

Our control measurements validated the procedure for estimating local flexural stiffness and revealed that the method is repeatable. Bending tests on a glass coverslip showed that a constant distribution of flexural stiffness provided the best approximation of measured displacement (as expected for a homogeneous beam). In addition, measurements of Manduca flexural stiffness in which each load was repeated several times provided consistent estimates of the EI distribution each time that a given load was applied (see Combes, 2002).

Morphological measurements showed that female Manduca sexta are significantly heavier and have larger wings than males (body mass, P=0.004; wing area, P=0.003; wing span, P=0.002,chord length, P=0.004; wing mass, P=0.012), but aspect ratio and wing loading were not significantly different in the individuals sampled(aspect ratio, P=0.946; wing loading, P=0.262). Although overall flexural stiffness scales with wing size across a broad range of species (see Combes and Daniel,2003a), spanwise

$$EI$$
was not significantly different in male and female Manduca (dorsal, P=0.165; ventral, P=0.308; Fig. 3A,B). In the chordwise direction,
$$EI$$
was higher in females than in males on the dorsal side of the wing (P=0.021), but not on the ventral side(P=0.477). The exponents of the stiffness distribution in the chordwise direction were significantly higher in females than in males(dorsal, P=0.006; ventral, P=0.008), but exponents in the spanwise direction were the same in males and females (dorsal, P=0.165; ventral, P=0.089).

Fig. 3.

Flexural stiffness distribution in wings of Manduca sexta and Aeshna multicolor, and in finite element models of M. sextawings. In each graph, spanwise flexural stiffness is shown above (longer lines) and chordwise flexural stiffness below; dorsal measurements are in black and ventral measurements in gray. Each line within these groups represents the flexural stiffness distribution estimated from wing displacement measurements performed on a different individual. (A) Flexural stiffness distribution of male Manduca sexta forewings (N=9 spanwise, N=10 chordwise), and of finite element models with homogeneous (blue) and exponentially declining (red) vein and membrane material stiffness. (B) Flexural stiffness distribution of female Manduca sexta forewings (N=4 spanwise, N=9 chordwise). (C)Flexural stiffness distribution of male Aeshna multicolor forewings(N=8). (D) Flexural stiffness distribution of male Aeshna multicolor hindwings (N=8).

Fig. 3.

Flexural stiffness distribution in wings of Manduca sexta and Aeshna multicolor, and in finite element models of M. sextawings. In each graph, spanwise flexural stiffness is shown above (longer lines) and chordwise flexural stiffness below; dorsal measurements are in black and ventral measurements in gray. Each line within these groups represents the flexural stiffness distribution estimated from wing displacement measurements performed on a different individual. (A) Flexural stiffness distribution of male Manduca sexta forewings (N=9 spanwise, N=10 chordwise), and of finite element models with homogeneous (blue) and exponentially declining (red) vein and membrane material stiffness. (B) Flexural stiffness distribution of female Manduca sexta forewings (N=4 spanwise, N=9 chordwise). (C)Flexural stiffness distribution of male Aeshna multicolor forewings(N=8). (D) Flexural stiffness distribution of male Aeshna multicolor hindwings (N=8).

The wings of all Manduca sexta displayed dorsal/ventral asymmetry;the average flexural stiffness of their wings was higher in response to forces applied to the dorsal side of the wing, in both the spanwise and chordwise direction (spanwise, P=0.002; chordwise, P<0.001; Fig. 3A,B). In contrast, Aeshna multicolor displayed no dorsal/ventral difference in its hindwings or in the chordwise direction of its forewings (hindwing spanwise, P=0.327; hindwing chordwise, P=0.735; forewing chordwise, P=0.069; Fig. 3C,D). In the spanwise direction, the forewings were stiffer on the dorsal side than on the ventral side (P=0.018), but the difference between dorsal and ventral flexural stiffness was far smaller than that seen in Manduca.

### Finite element modeling

The finite element model with diagonal strips of declining material stiffness accurately reproduced the sharp decline in flexural stiffness measured in male Manduca wings(Fig. 3A, red lines), while the model with homogeneous vein and membrane regions displayed a vastly different pattern of flexural stiffness (Fig. 3A, blue lines).

The effects of these differences in EI distribution were apparent in static tests on the model wings. Although an applied point force resulted in the same tip displacement, most bending in the exponential wing was confined to the outer third of the wing, while the homogeneous wing bent gradually along much of its span (Fig. 4A). When subjected to a pressure load, the homogeneous wing again bent along a large portion of its length, and its maximum displacement was nearly twice as large as the displacement of the exponential wing(Fig. 4B). In the exponential wing, displacement was localized to the tip and trailing edge of the wing, and bending was apparent in both the spanwise and chordwise directions.

Fig. 4.

Results of static bending tests on finite element model (FEM) wings. Wings are fixed at the base, and displacement from original position (black outline)in the z-direction is indicated by the color bar. The FEM wing with homogeneous membrane and vein material stiffness is shown on the left and the wing with exponentially declining material stiffness is shown on the right.(A) Displacement due to a point force of 0.003 N (green arrow) at the wing tip. (B) Displacement due to a normal face load of -14.43 Pa (green arrows),the approximate pressure on a wing during steady flight.

Fig. 4.

Results of static bending tests on finite element model (FEM) wings. Wings are fixed at the base, and displacement from original position (black outline)in the z-direction is indicated by the color bar. The FEM wing with homogeneous membrane and vein material stiffness is shown on the left and the wing with exponentially declining material stiffness is shown on the right.(A) Displacement due to a point force of 0.003 N (green arrow) at the wing tip. (B) Displacement due to a normal face load of -14.43 Pa (green arrows),the approximate pressure on a wing during steady flight.

The effects of an exponential decline in flexural stiffness were further illustrated by dynamic tests on the models. The homogeneous wing showed little chordwise bending and spanwise bending was most pronounced near the base,while the distal portion of the wing remained relatively rigid(Fig. 5, blue wings). In contrast, the exponential wing bent considerably in the chordwise direction,and bending in the spanwise direction was confined mainly to the outer portion of the wing (Fig. 5, red wings). For a movie of the model wings in motion, see http://faculty.washington.edu/danielt/movies.

Fig. 5.

Sequence of images from flapping finite element model wings with homogeneous (blue) or exponentially declining (red) vein and membrane material stiffness. Wings are moving to the right, viewed from their leading edge. Models are shown near the end of the simulation at the same time steps,beginning at 0.4325 s andproceeding in 0.0029 s intervals. For a movie of the model wings in motion, see http://faculty.washington.edu/danielt/movies.

Fig. 5.

Sequence of images from flapping finite element model wings with homogeneous (blue) or exponentially declining (red) vein and membrane material stiffness. Wings are moving to the right, viewed from their leading edge. Models are shown near the end of the simulation at the same time steps,beginning at 0.4325 s andproceeding in 0.0029 s intervals. For a movie of the model wings in motion, see http://faculty.washington.edu/danielt/movies.

Modeling passive wing flexibility requires some knowledge of the spatial distribution of stiffness throughout a wing. Although few studies have attempted to quantify spatial patterns of stiffness, a number of qualitative hypotheses have been put forward concerning the relationship between the structure of insect wings and their regional mechanical behavior. Most wings have relatively stiff, supporting zones near the wing base and leading edge,and relatively deformable areas near the edges of the wing(Wootton, 1981). Wing veins taper in diameter from base to tip(Wootton, 1992), and the presumed reduction in stiffness that this provides could serve many purposes,including reducing the inertial load at wing tips (thus reducing energy expenditure and stress at the wing base) and allowing wing tips to buckle and rebound from collisions without damage(Wootton, 1992). In addition,a wing that is stiff at the base and more flexible at the tip would be well suited to withstand the bending torques (which decline sharply from base to tip) generated by flapping (Ennos,1989).

Our measurements of spatial variation in flexural stiffness agree with these qualitative observations; flexural stiffness declines sharply from wing base to tip and from the leading to the trailing edge in the wings of both Manduca sexta and Aeshna multicolor. Displacement in the spanwise and chordwise directions of these wings can be predicted fairly well when the flexural stiffness distribution is approximated by an exponential equation. The observed similarities in flexural stiffness distribution in these two species (despite large differences in wing shape and venation pattern) suggest that a sharp decline in wing flexural stiffness towards the tip and trailing edge may be a common feature of many insect wings.

Although detailed inter-specific comparisons are not feasible with only two species, this study demonstrates several interesting intra-specific differences, particularly in Manduca sexta. Male and female Manduca sexta appear to differ in the spatial patterns of flexural stiffness in their wings; chordwise flexural stiffness declines far more sharply in male moths than in female moths(Fig. 3A,B), and as a result,the wings of female moths are significantly stiffer in the chordwise direction(dorsal

$$EI$$
is higher). The functional significance of these differences is unclear, but other aspects of Manduca wing morphology (such as the physical coupling of fore- and hindwings) also display sexual dimorphism(Eaton, 1988), which may be related to higher wing loading in females carrying eggs(Willmott and Ellington,1997b). The traits measured in female hawkmoth wings could in some way compensate for this intermittent higher wing loading, for example by influencing the pattern or extent of chordwise wing deflections during flight.

The wings of both male and female Manduca sexta display a large dorsal/ventral difference in average flexural stiffness, and in the exponents of flexural stiffness distribution. Average flexural stiffness is greater in response to forces applied on the dorsal side than on the ventral side (in both the spanwise and chordwise direction). In the spanwise direction, this difference is greatest near the base (Fig. 3A,B). Steppan(2000) found a similar dorsal/ventral asymmetry near the base of butterfly wings, and higher dorsal flexural stiffness has also been measured in the leading edge of locust hindwings (Wootton et al.,2000).

Although measurements of such bending asymmetry are rare, theoretical studies have suggested that insect wings might display dorsal/ventral stiffness asymmetry due to the camber inherent in most wings(Ennos, 1995, 1997; Wootton, 1993). We explored the effects of camber on bending asymmetry in the finite element model of a Manduca wing, adding 4% spanwise and 5% chordwise camber (the maximum values measured in real wings). Adding this degree of camber to the model wings had a relatively small effect on tip and trailing edge displacement (see Combes and Daniel, 2003a), and failed to produce any dorsal/ventral bending asymmetry (tip displacement was identical whether the force was applied to the convex or the concave side; see Combes, 2002).

Camber may in fact influence the maximum load each side of the wing can resist before buckling (the buckling load), and the finite element model does not simulate elastically stable buckling. However, our measurements on real wings do not display the strong non-linear dependence on loading that would be expected if elastically stable buckling were occurring during these static tests. In a more dynamic situation (such as flapping flight), where continuous shape changes can alter flexural stiffness, camber may in fact lead to dorsal/ventral differences in wing deformation. However, wing camber does not appear to explain the large differences measured in this study on wings subjected to static point loading.

The structural source of this bending asymmetry remains unclear. Although Manduca wings have none of the gross anatomical features often associated with bending asymmetry (e.g. a ventral flexion line), their veins may contain one-way hinges or other micro-structural features that facilitate asymmetric bending in other insect wings(Wootton, 1981; Wootton, 1992). In addition,the stress-stiffening effect explored by Kesel et al.(1998) in models of dipteran wings could potentially apply to Manduca. When unloaded, the membranes of Manduca wings lie rather loosely between the veins extending to the trailing edge. A point force applied to the dorsal (convex)side of the trailing edge appears to pull the trailing edge veins further apart, which would first remove any slack in the membrane, and then possibly increase the stiffness of the wing as force is applied to the already-taut membrane (the stress-stiffening effect). In contrast, when a point force is applied to the ventral side of the trailing edge, the veins are pushed together and the membranes between them become looser.

Although the finite element models do not reproduce this dorsal/ventral asymmetry, the strikingly different bending patterns of the two models demonstrate that flexural stiffness distribution is critical in determining how insect wings bend. In the model wing with exponentially declining flexural stiffness, both static and dynamic bending are limited to the tip and trailing edge, whereas bending in the homogeneous model occurs along the entire span.

These results suggest that the sharply declining flexural stiffness measured in real wings helps maintain rigidity near the wing base (despite larger bending moments), while localizing bending to the tip and trailing edge, which are regions of particular importance in controlling aerodynamic force production. The trailing edge is a critical control surface in airplane wings (affecting both the magnitude of lift and the lift-to-drag ratio; Anderson, 1991), and passive deformations in this region are likely to play a similarly important role in controlling flight forces in insects. The integration of passive flexibility and spatial patterns of flexural stiffness into models of insect flight will be an important step in determining the functional significance of wing structure and dynamic bending to insect flight performance.

To calculate local flexural stiffness, we treated each wing as a heterogeneous, two-dimensional beam in the spanwise (or chordwise) direction,and used a continuous beam equation to estimate flexural stiffness along this axis. With the beam attached at one side and a point force applied at the other end, flexural stiffness at any point along the beam is the moment divided by the local curvature (or the second spatial derivative of displacement):
$\ EI(x)=\frac{M(x)}{\mathrm{d}^{2}{\delta}(x){/}\mathrm{d}x^{2}},$
A1
where M(x) is the moment and δ(x) the displacement at position x. If we assume that the point force is applied at an angle approximately 90° from the beam, then the moment is the product of the applied force and the distance between the point of force application and position x on the beam:
$\ EI(x)=\frac{F(L-x)}{\mathrm{d}^{2}{\delta}(x){/}\mathrm{d}x^{2}},$
A2
where F is the applied force and L is the total length of the beam.

If tip displacement is extremely large relative to beam length, the beam is bent down so far that its coordinates in the x-dimension while loaded are significantly different from its coordinates while unloaded, invalidating Equation A2. However, Equation A2 is valid for our measurements for the following reasons: (i) all measurements were standardized to a relative displacement (δT/L) of 0.05 in the spanwise direction and 0.08 in the chordwise direction, (ii) original measurements were rarely performed at relative displacements over 0.08 in the spanwise direction and 0.15 in the chordwise direction, and (iii) even in the most extreme cases(where relative displacement exceeded 0.15), the error in coordinates introduced by length changes in the x-dimension was less than 4%.

We solved Equation A2 in reverse, posing a functional distribution for EI(x), calculating the expected displacementδ(x), and using simplex minimization to find the parameters that minimize the difference between predicted and measured wing displacement. We proposed that stiffness in the spanwise or chordwise direction might be constant, or that its spatial pattern might be approximated by a linear,exponential, or 2nd-degree polynomial equation:
$\ EI(x)=\mathrm{d},$
A3
$\ EI(x)=mx+b,$
A4
$\ EI(x)=c{\cdot}\mathrm{exp}^{\mathrm{ax}},$
A5
$\ EI(x)=px^{2}+qx+r.$
A6
Each of the proposed EI distributions (Equations A3-A6) was inserted into the following equation to find the expected displacement along the span or chord given an applied moment [F(L-x)]:
$\ {\delta}(x)={{\int}_{0}^{L}}\left[{{\int}_{0}^{x}}\frac{F(L-x)}{EI(x)}\mathrm{d}x\right]\mathrm{d}x.$
A7
We used a simplex minimization program in Matlab that tests various values of the coefficients (for example, m and b in Equation A4) to find the specific EI distribution (for each equation type) that most successfully predicts wing displacement along the span or chord.
After estimating the spatial distribution of flexural stiffness, we calculated average flexural stiffness
$$EI$$
of the wing by integrating the continuous distribution EI(x):
$\ EI=\frac{1}{(x_{\mathrm{max}}-x_{0})}{{\int}_{x_{0}}^{x_{\mathrm{max}}}}EI(x)\mathrm{d}x,$
A8
where xo is the starting point (near the base or leading edge) and xmax is the end point of data collection.

List of symbols

• a

coefficient in exponential equation (A5)

•
• b

intercept in linear equation (A4)

•
• c

coefficient in exponential equation (A5)

•
• d constant

(Equation A3)

•
• E

material stiffness (Young's modulus)

•
• EI

flexural stiffness

•
• $$EI$$

average flexural stiffness

•
• f

flapping frequency

•
• F

force

•
• I

second moment of area

•
• L

beam length

•
• m

slope in linear equation (A4)

•
• M

moment

•
• p

coefficient in polynomial equation (A6)

•
• q

coefficient in polynomial equation (A6)

•
• r

coefficient in polynomial equation (A6)

•
• t

time

•
• x,y,z

spatial coordinates

•
• xo

starting point of data collection

•
• xmax

end point of data collection

•
• δ(x)

displacement at position

•
• x δT

displacement at tip or trailing edge

•
• δT/L

relative displacement of tip or trailing edge

•
• θ

rotation at wing base

•
• τ

time constant

•
• ω

circular frequency (2πf)

A. Trimble was critical in devising the method of measuring displacement by a laser ranging technique and in developing Matlab code to analyze loaded and unloaded images. M. Tu and E. Goldman generously provided useful comments and Matlab wisdom. J. Henry provided advice on collecting dragonflies and donated specimens to the study. J. Dierberger at MSC Software provided crucial troubleshooting of the FEM models, without which the dynamic simulations would not have been possible. This work was supported by NSF grant F094801 to T. Daniel, the John D. and Catherine T. MacArthur Foundation, an NSF graduate fellowship to S. Combes, and an ARCS fellowship to S. Combes.

Anderson, J. D., Jr (
1991
).
Fundamentals of Aerodynamics
. New York: McGraw-Hill,Inc.
Batchelor, G. K. (
1967
).
An Introduction to Fluid Dynamics
. Cambridge, UK: Cambridge University Press.
Bennett, L. (
1966
). Insect aerodynamics:vertical sustaining force in near-hovering flight.
Science
152
,
1263
-1266.
Bennett, L. (
1970
). Insect flight: lift and rate of change of incidence.
Science
167
,
177
-179.
Combes, S. A. (
2002
). Wing flexibility and design for animal flight.
PhD thesis
, University of Washington.
Combes, S. A. and Daniel, T. L. (
2001
). Shape,flapping and flexion: Wing and fin design for forward flight.
J. Exp. Biol.
204
,
2073
-2085.
Combes and Daniel (
2003a
). Flexural stiffness in insect wings. I. Scaling and the influence of wing venation.
J. Exp. Biol.
206
,
2979
-2987.
Combes, S. A. and Daniel, T. L. (
2003b
). Into thin air: Contributions of aerodynamic and inertial-elastic forces to wing bending in the hawkmoth Manduca sexta.
J. Exp. Biol.
206
,
2999
-3006.
Daniel, T. L. (
1987
). Forward flapping flight from flexible fins.
Can J. Zool.
66
,
630
-638.
Daniel, T. L. and Combes, S. A. (
2002
). Flexing wings and fins: bending by inertial or fluid-dynamic forces?
Int. Comp. Biol.
42
,
1044
-1049.
Dickinson, M. H., Lehman, F.-O. and Sane, S. P.(
1999
). Wing rotation and the aerodynamic basis of insect flight.
Science
284
,
1954
-1960.
Eaton, J. L. (
1988
).
Lepidopteran Anatomy
. New York: John Wiley & Sons.
Ellington, C. P. (
1984a
). The aerodynamics of hovering insect flight. I. The quasi steady analysis.
Phil. Trans. R. Soc. Lond. B
305
,
1
-15.
Ellington, C. P. (
1984b
). The aerodynamics of hovering insect flight. II. Morphological parameters.
Phil. Trans. R. Soc. Lond. B
305
,
17
-40.
Ellington, C. P., Van den Berg, C., Willmott, A. P. and Thomas,A. L. R. (
1996
). Leading-edge vortices in insect flight.
Nature
384
,
626
-630.
Ennos, A. R. (
1988
). The importance of torsion in the design of insect wings.
J. Exp. Biol.
140
,
137
-160.
Ennos, A. R. (
1989
). Inertial and aerodynamic torques on the wings of Diptera in flight.
J. Exp. Biol.
142
,
87
-95.
Ennos, A. R. (
1995
). Mechanical behavior in torsion of insect wings, blades of grass and other cambered structures.
Proc. R. Soc. Lond. B
259
,
15
-18.
Ennos, A. R. (
1997
). Flexible structures in biology.
Comments Theor. Biol.
4
,
133
-149.
Ennos, A. R. and Wootton, R. J. (
1989
). Functional wing morphology and aerodynamics of Panorpa germanica(Insecta: Mecoptera).
J. Exp. Biol.
143
,
267
-284.
Gorb, S. N. (
1999
). Serial elastic elements in the damselfly wing: mobile vein joints contain resilin.
Naturwissenschaften
86
,
552
-555.
Haas, F., Gorb, S. and Blickhan, R. (
2000a
). The function of resilin in beetle wings.
Proc. R. Soc. Lond. B
267
,
1375
-1381.
Haas, F., Gorb, S. and Wootton, R. J. (
2000b
). Elastic joints in dermapteran hind wings: materials and wing folding.
Arthropod. Struct. Dev.
29
,
137
-146.
Herbert, R. C., Young, P. G., Smith, C. W., Wootton, R. J. and Evans, K. E. (
2000
). The hind wing of the desert locust(Schistocerca gregaria Forskål). III. A finite element analysis of a deployable structure.
J. Exp. Biol.
203
,
2945
-2955.
Kesel, A. B., Philippi, U. and Nachtigall, W.(
1998
). Biomechanical aspects of the insect wing: an analysis using the finite element method.
Comp. Biol. Med.
28
,
423
-437.
Lighthill, M. J. (
1973
). On the Weis-Fogh mechanism of lift generation.
J. Fluid Mech.
60
,
1
-17.
Liu, H., Ellington, C. P., Kawachi, K., Van den Berg, C. and Willmott, A. P. (
1998
). A computational fluid dynamic study of hawkmoth hovering.
J. Exp. Biol.
201
,
461
-477.
Sane, S. P. and Dickinson, M. H. (
2002
). The aerodynamic effects of wing rotation and a revised quasi-steady model of flapping flight.
J. Exp. Biol.
205
,
1087
-1096.
Savage, S. B., Newman, B. G. and Wong, D. T.-M.(
1979
). The role of vortices and unsteady effects during the hovering flight of dragonflies.
J. Exp. Biol.
83
,
59
-77.
Smith, C. W., Herbert, R., Wootton, R. J. and Evans, K. E.(
2000
). The hind wing of the desert locust (Schistocerca gregaria Forskål). II. Mechanical properties and functioning of the membrane.
J. Exp. Biol.
203
,
2933
-2943.
Smith, M. J. C. (
1996
). Simulating moth wing aerodynamics: towards the development of flapping-wing technology.
AIAA J.
34
,
1348
-1355.
Smith, M. J. C., Wilkin, P. J., and Williams, M. H.(
1996
). The advantages of an unsteady panel method in modelling the aerodynamic forces on rigid flapping wings.
J. Exp. Biol.
199
,
1073
-1083.
Spedding, G. R. and Maxworthy, T. (
1986
). The generation of circulation and lift in a rigid two-dimensional fling.
J. Fluid. Mech.
165
,
247
-272.
Steppan, S. J. (
2000
). Flexural stiffness patterns of butterfly wings (Papilionoidea).
J. Res. Lepid.
35
,
61
-77.
Vest, M. S. and Katz, J. (
1996
). Unsteady aerodynamic model of flapping wings.
AIAA J.
34
,
1435
-1440.
Wainwright, S. A., Biggs, W. D., Currey, J. D. and Gosline, J. M. (
1982
).
Mechanical Design in Organisms.
Princeton, New Jersey: Princeton University Press.
Wilkin, P. J. and Williams, M. H. (
1993
). Comparison of the aerodynamic forces on a flying sphingid moth with those predicted by quasi-steady theory.
Physiol. Zool.
66
,
1015
-1044.
Willmott, A. P. and Ellington, C. P. (
1997a
). The mechanics of flight in the hawkmoth Manduca sexta. I. Kinematics of hovering and forward flight.
J. Exp. Biol.
200
,
2705
-2722.
Willmott, A. P. and Ellington, C. P. (
1997b
). The mechanics of flight in the hawkmoth Manduca sexta. II. Aerodynamic consequences of kinematic and morphological variation.
J. Exp. Biol.
200
,
2723
-2745.
Wootton, R. J. (
1981
). Support and deformability in insect wings.
J. Zool., Lond.
193
,
447
-468.
Wootton, R. J. (
1990
). The mechanical design of insect wings.
Sci. Am.
November,
114
-120.
Wootton, R. J. (
1992
). Functional morphology of insect wings.
Annu. Rev. Entomol
.
37
,
113
-140.
Wootton, R. J. (
1993
). Leading edge section and asymmetric twisting in the wings of flying butterflies (Insecta,Papilionoidea).
J. Exp. Biol.
180
,
105
-117.
Wootton, R. J. (
1999
). Invertebrate paraxial locomotory appendages: design, deformation and control.
J. Exp. Biol.
202
,
3333
-3345.
Wootton, R. J., Evans, K. E., Herbert, R. and Smith, C. W.(
2000
). The hind wing of the desert locust (Schistocerca gregaria Forskål). I. Functional morphology and mode of operation.
J. Exp. Biol.
203
,
2921
-2931.
Wu, T. Y. (
1971
). Hydromechanics of swimming propulsion. Part 1. Swimming of a two dimensional flexible plate at variable forward speeds in an inviscid fluid.
J. Fluid Mech.
46
,
337
-355.