Bridging two insect flight modes in evolution, physiology and robophysics – Nature

Ancestral state reconstruction

Muscle type labelling

We encoded the orders Odonata, Ephemeroptera, Dermaptera, Plecoptera, Orthoptera, Embioptera, Phasmatodea, Mantodea, Blattodea, Isoptera, Raphidioptera, Megaloptera, Neuroptera, Trichoptera, Lepidoptera and Mecoptera as synchronous and the orders Thysanoptera, Strepsiptera, Coleoptera and Diptera as asynchronous^3,5,6,51,52. Muscle type in Zoraptera remains unknown, and the orders Mantophasmatodea, Grylloblattodea, Siphonaptera are all wingless. The three remaining orders, Hemiptera, Psocodea and Hymenoptera are known to have both synchronous and asynchronous species^3,5,6,52,56. The muscle type for each tip species of the insect phylogeny and its associated reference are included as a raw data file (Supplementary Table 4, ‘Species muscle type with sources’). We also conducted a more detailed literature analysis on all orders.

Psocodea

Psocodea is a particularly notable clade that has the potential to strongly influence the ancestral state and single origin of asynchronous muscle. This order is historically considered to have mixed muscle types³ based on the muscle structure data for various species⁶. Of the tip species present in the phylogeny used in this study, we were able to discern muscle types based on muscle structural data or absence of wings (see Fig. 1). Muscle structure data for the longest-branched genus (Ectopsocus) are inferred from a closely related species investigated by Cullen⁶ and the wingless state of other tip species was determined through multiple sources^57,58. Additional confirmation of the states of this group was completed by cross referencing the phylogeny used in this study (Fig. 1) with a more densely sampled phylogeny of Psocodea⁵⁹ and searching for additional muscle type data for Psocodea species not represented in the phylogeny used in this study. We found that one species (Trogium pulsatorium) belonging to the most ancient suborder within the group (Trogiomorpha) is reported to have synchronous flight muscle⁶. The next most ancient subclade (suborder: Psocomorpha) within Psocodea is known to have multiple species with asynchronous flight muscle based on its structure⁶. Data from a more recently diverging clade, the Amphientometae infraorder (within the Troctomorpha suborder), are absent and the remaining species of the larger clade (Troctomorpha), where Amphientometae is nested, are wingless^57,58.

All evidence together suggests that the ancestral state of Psocodea is asynchronous. However, there remains uncertainty in this group due to poor muscle data and their large degree of winglessness. Using scaling relationships based on body mass and measurements of Psocodea wing sizes, it seems likely that most winged psocodean species fly with wingbeat frequencies well over 100 Hz (150–500 Hz), even when we allow body mass to differ by an order of magnitude^60,61,62 from 0.1 mg to 1 mg. Wingbeat frequencies above 100 Hz are strongly associated with the evolution of asynchronous flight muscle. Second, in support of our scaling argument, other authors report that all species of winged Psocodea are asynchronous based on their necessarily high wingbeat frequencies⁶³, directly conflicting with the muscle type data from Cullen⁶. Finally, the winged clade of unknown muscle type, Amphientometae, has a most recent common ancestor with the Psocomorpha clade, which does possess asynchronous muscle⁶. If Amphientometae does have asynchronous muscle as expected, this would most probably result in an asynchronous ancestral condition of the entire Psocodea order. Thus, multiple lines of evidence support asynchronous muscle type as the ancestral condition of this clade with possibly a single, independent reversion to synchrony within the Trogiomorpha clade, which is not present in the phylogeny used in this study.

Hemiptera

Although Hemiptera is another clade with interspecific variation in muscle type we are again confident in our reconstruction of the ancestral state. First, many species of Hemiptera were investigated in ref. ⁶ and the clade is relatively well sampled in the phylogeny used in this study. In addition, we mapped additional muscle type data from ref. ⁶ onto a more densely sampled Hemiptera phylogeny⁶⁴ for at least one species from most (19 out of 29) families. Here, we find that all investigated species from the Heteroptera suborder, which includes 20 (11 of which have been investigated) of the 29 families present in the phylogeny in Johnson et al.⁶⁴, use asynchronous flight muscle⁶. Second, two of the four longest branch families with the most ancient diverging Hemiptera suborder (Sternorrhyncha) are also asynchronous. The Sternorrhyncha suborder is only represented by two synchronous species in the phylogeny used in this study and thus is likely overweighted when inferring the ancestral condition of this group. Despite that, we still recover an ancestral condition of asynchronous muscle. Other authors reviewing known muscle types and flight neuromechanics also concluded that most Hemiptera species rely on asynchronous flight muscle⁶³.

Two additional phylogeny tips within Hemiptera from the genera (Xenophysella and Nilaparvata) did not have published muscle structure data from any species within the same family, warranting further investigation. First, we code the Xenophysella tip as wingless because the majority of investigated species (24 out of 25) from the Coleorrhyncha suborder that includes Xenophysella are reported to be flightless⁶⁵. Second, we code the Nilaparvata tip as asynchronous for the following reasons: Nilaparvata myofibril diameter has been reported as 1.8 μm in insects three days post-emergence⁶⁶, which is above the 1.5-μm threshold for differentiating synchronous from asynchronous muscle^6,56. From transmission electron microscopy of Nilaparvata, the sarcoplasmic reticulum appears to be sparse⁶⁶ which is a proposed hallmark of asynchronous muscle⁶. Despite this evidence, there is still some uncertainty about the Nilaparvata muscle type for the following reasons: a different species (Sogatella furcifera) from the same family (Delphacidae) is reported to have muscle with myofibril diameters ranging between 1.5 and 2.0 μm 3 days post-emergence⁶⁷, which also falls on the border of the diameter (1.5 μm) used to differentiate muscle type in this group⁶. Further, Nilaparvata belongs to the Fulgoroidea superfamily⁶⁴. The Fulgoroidea superfamily shares a most recent common ancestor with the Membracoidea superfamily based on the Hemiptera phylogeny in Johnson et al. (2018)⁶⁴, and Membracoidea contains species of both synchronous and asynchronous muscle type⁶, making the identification of Nilaparvata equivocal based on its phylogenetic position alone. However, the most direct histological evidence supports our classification of Nilaparvata as asynchronous.

The high variation of muscle type within Hemiptera makes this clade particularly interesting for future studies on the evolution of synchronous and asynchronous muscle physiology and structure. Based on our assessment of muscle type across Hemiptera, there appears to have been multiple reversions back and forth between the two types, where both types have likely evolved at least once from an ancestor of the other type. These bidirectional transitions within Hemiptera support the thesis that muscle physiology lies on a continuum rather than as two discrete types and may transition across the bridge in parameter space (Fig. 3b, e). Despite the diversity within Hemiptera, the reconstruction of the ancestral node is confidently asynchronous.

Hymenoptera

All Hymenoptera muscle types were assigned based on published muscle structures and supported by other investigations of muscle physiology. As noted above, variation in muscle myofibril diameter is directly related to muscle type, where myofibril diameters less than 1.5 μm are considered synchronous muscle⁶. Myofibril diameter was measured in 46 species distributed across the Hymenoptera phylogeny⁶⁸. All but one of the 46 species are reported⁵⁶ to have myofibril diameters greater than 1.5 μm. A second line of evidence that relies on the muscle being defined as ‘close-packed’ versus ‘fibrillar’ supports these results⁵⁶. Thus, while Hymenoptera is considered to be a group of mixed muscle type³, we find the presence of synchronous muscle to be relatively rare. In support of these conclusions, ref. ⁶³ also reports that most species of Hymenoptera rely on asynchronous flight muscle.

All other clades

All other insect orders are reported to be invariant in flight muscle type or are known to be completely wingless. Therefore, we relied on their invariant classification from published reviews^2,3,52 where authors used total evidence from both muscle structure and physiology data across species of each order to determine the muscle type classification of each order. However, we do code the single Zoraptera tip as unknown because we could not find any data on any species from this order. The coding of this tip as synchronous, which would be the most conservative classification, does not impact our results.

Implementation

We used a previously published molecular phylogeny grounded in fossil records spanning all insect orders¹⁵, which modifies the fossil calibration of the extensive insect phylogeny developed by Misof et al.⁶⁹. For ancestral state reconstruction, we assumed an equal rates model of evolution and used maximum-likelihood estimation to estimate the posterior probability of ancestral states using the Phytools R Package¹⁴. These analyses were performed in RStudio (v. 1.1.383) using R (v. 4.0.2). In Supplementary Discussion A we test other models of evolution using different character states, allowing all rates to be different rather than equal, and consider hidden rates models, using the Phytools and corHMM R Packages^16,70. The latter allows for heterogenous rates of evolution but quickly increases the number of parameters in the model^16,17,70. All models with a single rate class are consistent with the order-level reconstruction of synchronoy and asynchrony, including the single origin of asynchrony at node 200 (Thysanoptera, Hemiptera, Psocodea and holometabolous insects). Some models with multiple rate classes can produce more ambiguous reconstructions with more possible patterns but are overfit as assessed by Akaike information criterion values (see Supplementary Discussion A for further discussion).

We also determine the posterior probability that the Lepidoptera plus Trichoptera clade evolved from synchronous ancestry without reverting from a single asynchronous ancestor (for example, there were multiple, independent, more recent evolutions of asynchronous muscle for different synchronous orders). To do so, we first recorded the probability that the node representing the origin of asynchronous muscle (node 200) is synchronous (Supplementary Table 5, ‘Ancestral state posterior probabilities per node’). Next, we constrained that same ancient node (200) to be 100% synchronous and recorded the probability of synchrony in the next most recent node (218) in the branching path towards the origin of Lepidoptera. We iteratively continued this process, constraining each node between the origin of asynchronous muscle and the origin of the Lepidoptera clade (Node 255). We then multiplied through the probability that each node is synchronous to calculate the total probability that the Lepidoptera clade evolved from a synchronous ancestor (Fig. 1d).

Muscle physiology

Animals

M. sexta were obtained as pupae from a colony maintained at the University of Washington. Moths were kept on a 12 h:12 h light dark cycle. We used 6 female and 3 male adults, all 2–6 days post eclosion, with a mean body mass of 2.32 ± 0.46 g.

Experimental preparation

After a 30 min cold anaesthesia, we removed the head, wings, abdomen, legs and first thoracic segment from each moth to isolate the second and third thoracic segments. We then used digital callipers to measure the DLM length as the distance between the anterior and posterior phragma. These structures are the physiological attachment points of the DLMs. We measured a mean muscle length of 11.7 ± 0.5 mm.

We used a similar experimental paradigm as Tu and Daniel⁵⁵ for dynamic, whole muscle experiments on M. sexta DLMs. The key difference in our protocol is that we used a dual-mode ergometer (305 C Muscle Lever, Aurora Scientific) capable of prescribing a length trajectory while measuring the force necessary to follow that trajectory. For the anterior muscle attachment, we used cyanoacrylate glue to rigidly mount the anterior phragma to a custom three-dimensional printed ABS shaft, which was secured to our experimental table. For the posterior attachment, we attached a pair of tungsten prongs to the ergometer lever. We inserted these prongs at the invagination between the second and third thoracic segments. This ensures that the prongs adhere to the posterior face of the posterior phragma. Cyanoacrylate glue ensured a strong connection. In all preparations, we ensured that the anterior and posterior attachments were rigidly bonded following the experiment.

At this point, the intact second and third thoracic segments were rigidly mounted on our ergometer. We relieved any force buildup during this procedure by manually adjusting the ergometer length until the force reading was at zero. We then shortened the muscle by 2% because the in vivo muscle length during flight is 2% shorter than its length during rest⁴¹.

We next decoupled the anterior and posterior sections of the muscle by removing a transverse ring of exoskeleton. To minimize inertial loads on the ergometer, we removed all other muscles and any remaining exoskeleton on the anterior side of the thorax. We made sure to excise the ganglia to prevent spontaneous muscle activation.

To activate the muscles, we inserted two tungsten electrodes into the anterior end of the muscle by piercing the exoskeleton of the anterior phragma. We repeated this procedure for the posterior side by piercing the posterior end of the scutum. To ensure muscle viability, we maintained a steady drip of saline and held a constant temperature of 35 °C measured at the muscle.

Because Ca²⁺ is required for delayed stretch activation², we stimulated the muscles at 150 Hz to induce tetanus. We found that 150 Hz stimulation was the minimum stimulation frequency to establish a fused tetanus. Our experiments to measure delayed stretch activation consisted of a stretch–hold–release–hold cycle. We first maintained zero strain for 150 ms to enable a plateau in force, indicating constant activation (Fig. 2b), we stretched the muscle to in vivo strains (4.5%) at peak in vivo strain rates⁴¹ while measuring muscle force output. We calculated peak in vivo strain rate as \({\dot{\varepsilon }}_{0}=2\pi f{\varepsilon }_{0}\), where f is the wingbeat frequency (25 Hz) and \({\varepsilon }_{0}\) is peak in vivo strain. We then held the strain constant for 150 ms. We returned the muscle to rest length at the same strain rate as before.

We measured the peak twitch, delayed stretch activation and tetanic force produced by the muscle. In all cases, we normalized force production by the cross-sectional area of the muscle. To determine the rate constants of delayed stretch activation we fit the muscle force data with the equation

where K_i and r_i are the coefficients and rate constants associated with a particular phase of the delayed stretch activation response: a fast decay (r₂ ≫ wingbeat frequency), a slower rise (r₃ ≈ wingbeat frequency), and a very slow decay (r₄ ≪ wingbeat frequency). The constant c represents the passive stiffness of the muscle. Phase 1 is the immediate material response of any muscle to a transient strain and is not relevant to characterizing delayed stretch activation. We quantified delayed stretch activation force magnitude as the difference between the lowest force immediately following stretch to the peak force during the plateau phase.

Following experiments, we removed and weighed the DLMs. We measured an average muscle mass of 0.123 ± 0.012 g. This corresponds to a body mass-normalized muscle mass of 5.28 ± 1.41%, which is in rough agreement with prior measurements⁵⁵ of 5.96 ± 0.62%. From muscle length and muscle mass, we calculated a cross-sectional area of 10.7 ± 2 mm² under the assumption that muscle density is 1 g cm⁻³.

All conditions were replicated on all individual preparations. Twitches were always done first to confirm that stimulation of muscle preparations produced force and that the muscle was viable. No further randomization was necessary because all data were collected from a single continuous ramp-hold experiment. Experiments were not blinded because there were not multiple conditions.

Accounting for non-ideal strain rate

Our input stimulus to the ergometer for stretch–hold–release–hold experiments was a ramp with a speed matching the in vivo strain rate of muscle contraction in a hawkmoth. Our modelling assumes that r₃ is the rate of tension rise in response to an infinite impulse, which is not possible to implement in any real physical system. To examine the discrepancy between muscle’s response to an infinite impulse and a non-ideal finite impulse, we follow the following procedure. First, we construct a rectangular pulse with a width and height that match the width and height of the actual strain rate pulse we imposed in experiment. We then compute the empirical transfer function between the sum-of-exponentials fit to our force data and the rectangular pulse. This transfer function represents the response of hawkmoth muscle to an infinite impulse in strain rate, and is equal to our fit multiplied by a scalar. We compute this scalar to be 2.29, by dividing the IIR by our fit. We then scale our force data by this constant factor and plot it, labelled as ‘IIR’ in Fig. 2e. Data presented in Fig. 2b–d are raw and unscaled by the method described here.

Delayed stretch activation model

To study how delayed stretch activation produces wing oscillations we needed to generate a feedback model for delayed stretch activation. No current detailed muscle model can predict both neural and stretch-activated force components under general dynamic conditions, in part because of limitations in our understanding of the multiscale interactions in muscle⁷¹. Thus, to model asynchronous force, we do not try to build a detailed molecular model that can predict force from arbitrary activation and strains. Instead, we seek to capture the basic functional input–output relations for delayed stretch activation between an imposed strain and the resulting force. We first constructed a reduced order model of the delayed stretch activation that was able to capture the stretch–hold–release–hold behaviour we observed in experiment. This single-parameter model is described by the ‘time to peak’ (t₀) of the delayed stretch activation force response, which allows us to study how the relative timescales of delayed stretch activation (t₀), body mechanics (the natural resonance period, T_n), and the synchronous timescale \(\left(\frac{1}{{f}_{{\rm{s}}}}\right)\) govern the emergent wingbeats. To implement delayed stretch activation in simulation, we then generated a computational representation of delayed stretch activation and coupled it to a computational representation of body mechanics.

One-parameter model of asynchrony

Measurements of delayed stretch activation in the literature⁷² consist of imposing a step change in muscle length and fitting the force response, F_step(t), with a sum of three exponentials given by equation (5). This seven-parameter model was used to fit the delayed stretch activation response in the hawkmoth muscle (Fig. 2c). However, the initial viscoelastic drop following lengthening is unlikely to be important for generating self-excited oscillations, and the symmetry between r₂ and r₄ leaves those parameters sensitive to initial conditions of a curve fit procedure. Because the delay between stretch and force production is likely the critical feature of delayed stretch activation, we only considered the delayed tension rise (determined by r₃) and the delayed tension drop (determined by r₄). In doing so, we eliminated K₂ and r₂ from the fit equation (5). We further eliminated the passive muscle stiffness from equation (3) because it does not contribute significantly to the elastic response of the thorax in M. sexta³⁸. It is possible that active muscle stiffness contributes to body mechanics. Incorporating estimates of active muscle elasticity and body dissipation that may be present in the thorax does not affect the overall conclusions or features of the simulation (see Supplementary Discussion D). We fit the hawkmoth muscle data to the reduced convolution kernel

where g₀ is a scalar that normalizes by the area under the kernel (which depends on the kernel rate constants and has units of seconds). We sought to further reduce this convolution kernel to be parameterized by a single variable, the time to reach peak tension from a step input (t₀; Fig. 1b). We first assumed a constant ratio between r₃ and r₄ such that r₄ = κr₃ (for M. sexta, κ = 0.62). We can then solve for t₀, from the reduced kernel to obtain

For hawkmoth muscle this yields a relationship between t₀ and r₃ of

The two-parameter fit model yields an r₃ of 36.39 ± 0.09 s⁻¹ and r₄ of 22.80 ± 0.04 s⁻¹, a corresponding t₀ = 0.034 ± 0.001 s, and this two-parameter fit matches the experimental data well (Extended Data Fig. 4a). While our single-parameter model is simplified from the classic seven-parameter delayed stretch activation model, the qualitative features of the Fig. 3 heat maps and associated conclusions are insensitive to the precise value of κ.

Discrete FIR filter implementation

To implement delayed stretch activation in our simulation and robot experiments, we require a model of delayed stretch activation that can produce stretch dependent forces in real-time which can provide a delayed stretch activation force to the simulation or robotics experiments. We describe the delayed stretch activation response of asynchronous muscle as a convolution of the muscle strain velocity with a kernel g such that

where F_a is the asynchronous forcing magnitude and dictates the strength of the delayed stretch activation feedback taken from the quasi-static experiments, and µ is a fitting parameter that scales the stretch activation response to flight conditions.

We implement this convolution in simulation and in the robotic models in MATLAB simulink (Mathworks) using a finite impulse response (FIR) filter, which is an instantaneous (real time at 10 kHz) evaluation of a convolution operation. We construct the filter such that the input is the muscle strain rate, and the output is the delayed stretch activation force. We convert from the angular rotational units of our wing to actuator strain through the equation \(\dot{\phi }=LT\dot{\varepsilon }\), where we divide wing velocity (\(\dot{\phi }\)) by a factor LT where L is the resting muscle or actuator length (in the robot models, L = T = 1 because all scaling can be captured in µ). This yields the following delayed stretch-activated force

The value of µ is tuned to each system (simulation, roboflapper and robobee wing; see below for details), but is not changed between experiments where the relative magnitudes of synchronous and asynchronous forcing are varied (that is, when K_r is varied as in Fig. 3).

For simulations and experiment the delayed stretch activation force is generated from an FIR filter which requires a numerical evaluation of the convolution \((\,-\,g\ast \dot{\phi })(t)\). First, we generate the response curve g(t) based on the delayed stretch activation parameters (r₃, r₄) in MATLAB (Mathworks), sampled at the system rate ∆t over the simulation/experiment duration (simulation, roboflapper and robobee wing; see below for details). The normalization parameter g₀ in equation (6) is the numerical area under the curve. We then find the value of t for which g(t) < 0.01 × max(g(t)) and truncate the response since the finite impulse response filter requires a finite kernel. Last, we multiply the output of the numerical convolution by the coefficient \(\frac{\mu {F}_{{\rm{a}}}}{LT}\). The delayed stretch activation step response is thus represented as a vector of numbers that are supplied to the ‘filter coefficients’ input of the FIR block in Simulink. When synchronous and asynchronous forces are applied together we scale F_async by (1 − K_r) according to the combined forcing equation (equation (1)).

Hawkmoth simulation

Hawkmoth body mechanics model

We used a dynamics model of the wing rotation in the stroke plane (ϕ) for the hawkmoth which has previously been derived by Gau et al.³⁸. In brief, we assumed a body mechanics model with aerodynamic drag whose magnitude depends on angular velocity squared (\(\left|\dot{\phi }\right|\dot{\phi }\)) with a constant drag coefficient (Γ), a parallel-elastic spring due to thorax elasticity (k), and rotational inertia from wing and added mass (I). These assumptions yield the equations of motion presented in equation (3). To generate equivalent torques from the linear muscle force and the linear thorax elasticity we require the transmission ratio between linear muscle displacement and the rotational wing movement. We assumed a linear transmission ratio such that T = ϕ/X, where ϕ is the wing rotation and X is the linear displacement of the muscle and thorax. We calculated the transmission ratio for M. sexta as T = ϕ₀/X₀, where ϕ₀ is the peak-to-peak wingstroke amplitude and X₀ is the peak-to-peak muscle displacement amplitude (values can be found in Supplementary Table 1 and ref. ⁴²). The equivalent torque about the wing hinge produced by the muscle force F_m is given by F_m/T. The equivalent elastic torque from the thorax linear stiffness is calculated as k/T². These two transformations can be derived using conservation of energy: work done at the rotational joint must equal work done on the linear elements (spring and muscle).

The wing inertia, I, includes the added mass effects from aerodynamics. The parameter values can be found in Supplementary Table 1. Following the derivations of Ellington⁷³, wing inertia (I) is the sum of inertia due to wing mass (I_w) and added mass (I_a),

where R₂(m) is the radius of the second moment of wing mass, R₂(v) is the radius of the second moment of wing added mass, and L_w is the wing length. Note that in the aerodynamics literature, the second moments are often denoted r₂ and the wing length denoted R, but we change the convention here to avoid confusion with the rate constants r, used in the delayed stretch activation experiments. Dimensional added mass (v) is defined as \(v=\frac{2\rho \pi \widehat{v}{L}_{{\rm{w}}}^{2}}{{({\rm{AR}})}^{2}}\), where \(\widehat{v}\) is the non-dimensional added mass of the wing pair and AR is the aspect ratio of the wings. Parameter values are in Supplementary Table 1.

The lumped aerodynamic parameter Γ was calculated by following the work of Whitney and Wood⁷⁴. The quasi-steady aerodynamic drag force (F_aero) on insect wings over a single wingstroke can be modelled as

where R₂(s) is the radius of the second moment of wing shape. Setting the drag torque τ_aero = F_aerol_cp and \({\tau }_{{\rm{aero}}}=\Gamma \left|\dot{\phi }\right|\dot{\phi }\), where l_cp is the centre of pressure⁷⁵, yields the velocity-squared aerodynamic damping coefficient (Γ) as

Simulation details

We used MATLAB and Simulink (Mathworks) to run simulations of combined synchronous and asynchronous forcing on a mechanical model of the hawkmoth. Extended Data Fig. 4b,c presents a representation of the Simulink model. The system dynamics block implements the equation of motion (equation (3)) using hawkmoth parameters. It takes the combined muscle forcing as an input and generates the wing angle and angular velocity as outputs (Extended Data Fig. 4b). The wing rotational velocity is then an input into the delayed stretch activation simulation described in the previous section and calculated by equation (11).

As insect flight is driven by pairs of antagonistic muscles, we represent the upstroke and downstroke muscles separately in our Simulink simulation. The antagonistic configuration means that the sign on both strain velocity and output force is different for each muscle, as shown in Extended Data Fig. 4c. Additionally, a sine wave generator is used to produce synchronous forcing based on the amplitude F_s and frequency f_s. The output force is a weighted sum of synchronous and asynchronous forces defined by K_r.

Both the synchronous and asynchronous forces in the muscle block are saturated so that they only output tension forces. Additionally, the sine wave generators are operated with different initial phase θ₀ = 0 for the upstroke muscle and θ₀ = π for the downstroke muscle. The overall effect is that all of the negative torque is produced by the upstroke muscle, and all of the positive torque is produced by the downstroke muscle. The phase shift in the sine wave generator blocks also enforces that the fully synchronous output is identical to a single sinusoidal torque source.

Iterative force tuning procedure

The parameter K_r describes the relative amounts of synchronous and delayed stretch-activated forcing. To study how an insect that is actuated purely through delayed stretch activation (K_r = 0) can transition to being purely actuated through synchronous forcing (K_r = 1) we need to establish values of µF_a and F_s that produce feasible wingbeat motions in both of these regimes. In the hawkmoth simulation we determined that a sinusoidal forcing amplitude of F_s = 2,720 mN generates wingbeat kinematics that match in vivo observation of 117 degrees peak-to-peak. This value was previously used to synchronously drive an identical simulation to physiological wingbeat amplitudes⁴².

However, the wingbeat kinematics in the purely asynchronous regime (K_r = 0) are emergent and thus we need to determine an appropriate F_async that can drive our insect model to appropriate wingbeat kinematics. We used a simple iterative force tuning procedure to determine the value of µ such that asynchronous actuation (K_r = 0) can produce wingbeats with peak-to-peak amplitude of ϕ₀ = 117°. We slowly increment the value of µ until the output steady-state wingbeat amplitude is within 1% of the desired amplitude of ϕ₀. In this way we ensure that the both synchronous (K_r = 1) and asynchronous (K_r = 0) actuation can produce the same wingbeat amplitude.

Calculation of K
_r for M. sexta

Direct computation of K_r for M. sexta is challenging since realistic measures of muscle force (for example, from work loop experiments) will contain a mixture of synchronous and asynchronous effects, which are unlikely to be distinguishable under flight conditions. To estimate K_r from physiological measurements, we approximate the relative contributions of synchronous and asynchronous force using quasi-static measurements of neurogenic force and stretch activation. We use the maximum stretch-activated force above baseline from our stretch–hold–release–hold experiments, F_a, as a static representation of the asynchronous muscle force when the amplitude and rate of the stretch is equivalent to those during flight. The tetanic muscle force F_tet is a static representation of the maximum neurogenic (synchronous) muscle force at operating length. Both F_a and F_tet are generated under maximum activation, and F_tet is measured isometrically, so it will not include any stretch-activated force. We can then approximate K_r as \({\widetilde{K}}_{r}\), which is the proportion of neurogenic force (F_tet) to total force (F_a + F_tet). This gives the ratio in equation (2), which is equal to 1 when there is no stretch activation and approaches zero when the stretch activation far exceeds neurogenic force.

This approximation is based on equating the proportion of synchronous and asynchronous contributions during flight conditions to those measured at maximum activation in quasi-static conditions. While the absolute value of synchronous and asynchronous forces are likely to be very different in flight and quasi-static conditions, their relative importance is likely to be more comparable. For example, if \({\widetilde{K}}_{r}\) is related biologically to the proportion of troponin isoforms that are calcium-activated, then this proportion would likely affect quasi-static and flight conditions similarly. Nonetheless, \({\widetilde{K}}_{r}\) is still an approximation and we provide it as a way of estimating this parameter from currently available experiments. This approach also enables future experiments to characterize the relative contributions of synchronous and asynchronous force magnitudes using standardized experimental methods. The simulations and models here do not depend on M. sexta having a particular value of K_r. Future work would benefit from trying to parse the physiological contributions of synchronous and asynchronous forcing in flight conditions and resolving their specific molecular correlates.

Simulation parameter sweep and analysis

To evaluate how the presence of both synchronous and delayed stretch activation in an insect muscle influences the wing kinematics we performed simulations varying both K_r and the time to reach peak tension of the delayed stretch activation, t₀. We incorporated mechanical timescales of the system by dividing t₀ by the natural period to yield the parameter t₀/T_n. For a given insect, T_n is assumed to be constant. To sweep across delayed stretch activation timescales, we adjusted r₃ via equation (9) to sweep over a range of t₀/T_n values from 0.01 to 1. We varied K_r from 0 to 1. For each set of K_r and t₀/T_n values, we first calculated r₃ from t₀ (equation (9)). We then generated the delayed stretch activation kernel as described in the delayed stretch activation model section. With µF_a from the section above, we could now combine synchronous and asynchronous forcing (equation (1)). We initialized the wing position at 0.1 rad to initiate oscillations when there was no synchronous forcing. All simulations were performed with a fixed sample time of ∆t = 1 × 10⁻⁴ s over a duration of 5 s.

For each set of parameter values, we recorded the emergent force F_m, wing position, and wing velocity. We determined the emergent oscillation frequency by taking the Fourier transform of the last 2.5 s of position and identifying the frequency with the largest magnitude. To calculate power, we extracted five periods of oscillation after the system reached steady state. We then numerically integrated force over position and divided by the time elapsed. Lastly, we computed the variation in the peak-to-peak wing amplitude by using the findpeaks command in Matlab to locate all of the wingbeat peaks. The amplitude variation is calculated as the standard deviation of the peak-to-peak wingbeat angles.

Robophysical experiment

Robot details

Experiments were performed on a dynamically scaled robophysical model described previously in Lynch et. al.³⁹. The device consists of a silicone torsion spring with known, linear characteristics³⁹; a brushless DC motor (ODrive Robotics, D6374) under closed-loop torque control; and a rigid, fixed-pitch acrylic wing submerged in a tank of water (Extended Data Fig. 4). The wing span and chord (10 × 3.6 × 0.5 cm) were selected such that the wing, flapping in water with an amplitude between 10° and 60° and frequency between 1 and 4 Hz, has a Reynolds number between 10³ and 10⁴, which is approximately the same range as M. sexta²³. Friction is minimized via a set of radial air bearings and a thrust ball-bearing. We measured the spring stiffness and system inertia and calculated the wing drag torque coefficient (Supplementary Table 2)³⁹. We also calculated the natural period T_n of the robophysical system using equation (4) as T_n = 0.416 s.

The robophysical experiment was designed to mimic the hawkmoth simulation, replacing the virtual hawkmoth dynamics with those of a real system. We tracked wing angle using an optical encoder (US Digital, 4096 CPR) fixed to the wing shaft and a DAQ (National Instruments, PCIe 6323) sampled at 1 kHz. The encoder angle was used as input to a Simulink Desktop Real-Time (Mathworks) model running an identical combined forcing model as described previously in the sections on the delayed stretch activation model, and the hawkmoth simulations. The velocity was calculated by taking a derivative of the encoder position and fed into the model. The model prescribed a motor torque which was sent via USB serial connection to an open source motor controller (ODrive v3.6) that converted it to a current command to the brushless DC motor. The control loop for sensing wing position and sending torque commands to the motor ran at a sample time ∆t = 1 × 10⁻³ s.

Experiment details

To study how the robophysical system transitions between delayed asynchronous and synchronous forcing modes (Fig. 3d–f) we varied K_r and t₀/T_n in experiments. The robophysical experiments used approximately the same range of actuation parameters as the simulation: K_r spanning 0 to 1 and t₀/T_n from 0.02 to 1. The synchronous gain was set manually so that oscillations did not trigger the overload-current safety features of the motor driver, and the asynchronous gain was set using the same iterative force tuning procedure described above. We ran experiments for 30 s and we measured output power and frequency over the last 15 s of the experiment.

In a separate set of experiments, we studied the frequency entrainment properties of the robophysical system under combinations of both synchronous and asynchronous forcing. We first determined a value of µF_a in experiment that yielded high-amplitude asynchronous oscillations (106 ± 3° peak-to-peak) at 3.2 Hz. Next, we performed experiments with a constant µF_a, but with varied synchronous frequency f_s = [0.815, 6.515] Hz at three levels of forcing magnitude, \({F}_{{\rm{s}}}^{* }\,\)= [0.1, 0.2, 0.3]F_s, with respect to the purely synchronous forcing magnitude of F_s. We then measured the output wingbeat angle and computed: (1) the emergent frequency using the peak frequency of the Fourier transform; and (2) the peak-to-peak variation in wingbeat amplitude. The results of this experiment are shown in Extended Data Fig. 8.

Robobee experiment

We fabricated a single-winged version of a ‘dual-actuator’ Harvard Robobee, following the smart composite microstructure (SCM) process pioneered by the Harvard Microrobotics Lab^11,12. Wing parameter values are provided in Supplementary Table 3. The carbon fibre airframe, which holds the piezoelectric bending actuator and SCM transmission, was fixed to an acrylic mount on a manual translation stage to enable displacement sensor calibration.

To implement the delayed stretch activation model, it is necessary to estimate wing velocity in real time. We achieved this via a fibre-optic displacement sensor (D21, Philtec) pointed at a small piece of reflective tape glued to the bending actuator. The sensor is able to measure actuator displacement at which are fed into a Simulink model that converts sensor voltage to displacement through a calibration curve, and then takes a numerical derivative to calculate wing rotational velocity. Wing rotational velocity is then supplied to an identical Simulink model as in the hawkmoth simulations and roboflapper experiments described above. The simulation of delayed stretch activation force was converted into an amplified voltage signal (0 V < V_sig < 200 V) and sent to the piezoelectric actuator resulting in wing oscillations. The control loop for sensing wing position and sending torque commands to the motor ran at a sample time ∆t = 1 × 10⁻⁴ s. The asynchronous gain was chosen such that flapping angles were large but the actuator did not saturate, and the synchronous gain was set using the iterative force tuning procedure above. The Robobee flapping amplitudes did not exceed 50° peak-to-peak.

Observations of the Robobee wing angle were taken via a high-framerate video camera (Phantom VEO-410) at 2,500 frames per second. Video frames were processed in MATLAB to get the wing angle. Flapping amplitude was estimated by finding oscillation peaks, and flapping period or frequency was estimated by computing the time between peaks and smoothing the resulting curve (Extended Data Fig. 9).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.