Machine learning reveals the control mechanics of an insect wing hinge

Initial prediction of wing pose with using a convolutional neural network (Flynet)

Our pose estimation strategy required constructing an anatomically accurate 3D model of a fly that we would fit, frame-by-frame, to the data in our three high speed video sequences. We based our 3D model from images of a tethered flying fly taken at different angles, from which we extracted aspect ratios and curvatures of the head, thorax, and abdomen. These 2D contours were converted into 3D surfaces via Non-Uniform Rational B-splines (NURBS)[1] Like B-splines, NURBS surfaces may be manipulated by changing the locations of control points. By collapsing control points at the edges of NURBS surfaces onto each other, we created spheres for the head, thorax, and abdomen and then manipulated the control points to reshape the objects according to the aspect ratios and curvatures derived from the fly images.

Observations from high-speed videos show that the wings exhibit chord-wise deformations during flight that are concentrated at the L3, L4, and L5 veins (Fig. 1d). Accordingly, our model wing consists of four rigid panels connected by three hinge lines. To reduce the number of parameters representing wing deformation, we assumed that the deformation angle along each hinge line was equal (Fig. 2c). We thus modeled the total deformation angle, ξ, as the sum of the local deformation angle of each hinge line (ξ = ξ / 3 + ξ / 3 + ξ / 3), an assumption that was supported by manual inspection of the raw images.

Each wing pose vector consisted of 8 values: 4 parameters for a quaternion describing the orientation of the leading edge panel of the wing relative to the world reference frame, 3 parameters for a translation vector describing the distance between the wing root and the origin of the world reference frame, and the wing deformation parameter ξ. For the head, thorax, and abdomen pose vectors we used 7 parameters: 4 for the quaternion and 3 for the translation vector. Our total fly pose vector thus required 37 parameters to account for the head (7), thorax (7), abdomen(7), left wing (8), and right wing (8). Besides this pose vector, we defined a vector that contained 5 parameters for scaling the body and wing components independently.

We developed software, called Flynet, to implement the entire wing pose estimation process, including the manual annotation required for creating a training set. A module in the Flynet GUI constructed in Python with the PyQt and PyQtGraph packages allows the user to load the 3D fly model and scale, rotate, translate, and deform all the components. The GUI uses the DLT-calibration[2] of the high-speed cameras to project a wireframe of the fly model onto the three camera views. By adjusting the pose and scaling vectors of each component until it matches the images in each view, an accurate label can be created for each frame triplet. Using this GUI, we collected a manually annotated a final dataset of 2905 frame triplets and associated pose vectors. Neural networks work best with data centered around 0 with a standard deviation of 1. Thus, the 8-bit pixels of the high speed images were rescaled by dividing by 255. To reduce computation time, the normalized frames were cropped into 224 × 224 pixel images centered around the thorax of the fly.

The CNN was implemented in Python using the Tensorflow library with Graphics Processing Unit (GPU) support. The CNN architecture started with three convolutional blocks, one for each camera view. A graphical representation of the complete Flynet code (created using the the model plotting utility in Keras 3.0) is shown in Appendices A and B; the full code is available as described in the Data Availability Section of the main manuscript.

We randomly split our set of manually annotated data into groups of 95% and 5% for for training an validation, respectively. The network was trained for 1000 epochs with a batch size of 50. After some testing, we chose a dropout rate of 0.1 as it gave the lowest validation error. Lower dropout rates reduced the training error but increased validation error; higher dropout rates increased the training error and therefore also the obtainable validation error. Mean squared error (MSE) served as the loss function. We used a learning rate of 1.0 · 10⁻⁴ and a decay of 1.0 · 10⁻⁷. The decay value was subtracted from the learning rate after each epoch, resulting in a linear decrease of the learning rate during training. Using a Nvidia Geforce GTX 1080 graphics card with 8 GB of memory, it took approximately 24 hours to train Flynet. During the training phase, the MSE for the independent components of the pose vector decayed up to 400 epochs and remained constant afterwards (Extended Data Fig. 2b). Stabilization of the validation error indicates that the network’s performance did not improve with more epochs, even though the training error continued to decline. The weights of the trained network were saved in an hdf5 file and loaded in the Flynet GUI for the pose prediction.

Refining pose with Particle Swarm Optimization

The neural network predictions of body and wing pose were close to the actual wing pose as observed in the images, but not accurate enough to study the influence of muscle activity and resulting aerodynamic forces, which are extraordinarily sensitive to small changes in kinematics[3]. A larger training set might have resolved this issue, but many additional frames would likely have been required. Instead, we opted to refine the pose estimate in a second step using Particle Swarm Optimization (PSO)[4]. This approach avoids an extensive manual annotation process while still maintaining the benefits of time-independent tracking.

To implement the refinement step, it was necessary to segment the images and create binary masks of the body and wing pixels (Extended Data Fig. 2a). Instead of using a median image for background subtraction, we determined a minimum pixel intensity image over a batch of 100 frames. Body pixels were found by a simple threshold (intensity>200) of the raw, 8-bit image data. Wing pixels were identified by subtracting the minimum image from the frame and subsequently selecting all pixels that had an intensity >10. Subtracting the minimum image removed all stationary pixels, thereby identifying those corresponding to moving body parts.

Our refinement approach used area matching as a measure of how well the 3D model aligned with the images. The advantage of area matching is that it is simple and computationally efficient; however, the disadvantages are that the cost function is not continuous and gradient-based optimization algorithms will not work well. Fortunately, PSO is a gradient-free optimization algorithm that has the ability to find the global minimum, even when the cost function contains multiple local minima. PSO is not guaranteed to converge on the global minimum in a finite number of iterations, but for complex problems such as area matching, it is one of the few optimization methods available.

PSO relies on a swarm of particles that move through the state space with a position and a velocity. The state vector consists of all variables that affect the cost function. At the start of the PSO algorithm, the particles in the swarm are given random positions and random velocity vectors. For a set number of iterations, the algorithm updates the position and velocity of each particle according to the following equations: where i is the particle index, is the particle position at iteration k, is the particle velocity, r₁ and r₂ are random weight vectors drawn from a normal distribution, is the personal best of the particle, g^best is the global best for the whole particle swarm, and w, c₁, and c₂ are weights. A particle’s personal best is the position with the lowest value for the cost function throughout the travel history. The global best is the position with minimum cost value in the travel history of all particles of the swarm. In each iteration, the cost function is evaluated for the current position of each particle in the swarm. Subsequently, the cost function values of the personal and global best are compared to the cost evaluations for the current positions. If the current cost function of a particle position is lower than the personal or global best, these values and the associated position vectors are updated. Additionally, in each iteration vectors r₁ and r₂ are updated by randomized weights drawn from a normal distribution.

The velocity update of each particle (equation 2) consists of three parts. The inertia term, w is the previous velocity multiplied with weight, w. The second term points towards the personal best of the particle with some added random noise and multiplied by the scaling factor, c₁. The third term is a vector pointing towards the global best of the swarm with added random noise multiplied by a scaling factor, c₂. This scaled randomization introduces stochasticity in the search process, thus insuring that the particles will sample the space near a minimum more often. The vectors pointing towards the personal and global best positions help the particles converge towards a minimum value. The inertia term controls how quickly the swarm converges to a global best. If the inertia term is small, the swarm will converge quickly and risks getting stuck in a local minimum. Setting the inertia term too high, however, may result in no convergence at all or convergence only after a large number of iterations.

In our particular application of PSO, the cost function was based on the projected images of the 3D model. In order to speed up evaluation, the PSO algorithm and the cost function were implemented in C++. The PSO algorithm lends itself well for parallel processing, as the cost function evaluations of all particles can be executed at the same time. Using the std-library in C++17, a separate thread for evaluating the cost function was assigned to each particle. The cost function relies on matrix operations that were implemented using the Armadillo library. To integrate the C++ code within the Flynet GUI, the Boost library was used to make the functions callable in Python.

The 37-dimensions pose vector was a large parameter space in which to search. However, we simplified the process by splitting the pose vector into 5 different components (head, thorax, abdomen, left wing, and right wing), with 7 or 8 state parameters each. Although the pose vectors of the five model components are independent, the cost function evaluation is not, as different model components can overlap in the projected view of the 3D model. To address this dependence, the PSO algorithm updated one randomly selected component at a time.

The variables in the pose vector are all linear, except for the quaternion. A unit quaternion, ∥q∥ = 1, can be seen as a point on a 4D sphere with radius 1. However, when updating the quaternion according to equations 1 and 2, the result is not likely to be a unit quaternion. Determining the quaternion difference that satisfies the constraint of a 4D sphere with radius 1 requires application of the concepts of quaternion multiplication and the quaternion conjugate[5]. The details by which we implemented the PSO algorithm using quaternion operations are detailed elsewhere[6].

The cost function for each model component j, C _j, is given by: for body components, and for the wings, where I^b, and I ^w are the dynamic bitsets of the body and wing pixels respectively, is the dynamic bitset of the model component projections in all 3 views, and Δ and ∪ are the bitwise symmetric difference and union operators, respectively. We vectorized both the images masks and the projected images to accelerate these operations.

Before the tracking started, we set the number of particles, the number of iterations, and the standard deviation of the normal distribution that is used to perturb the initial state and search parameters. The initial particle pose vectors were created by adding random noise to the predicted pose vector from the CNN. Similarly, the particle velocities were randomly sampled from a normal distribution. By specifying the standard deviation of the normal distribution used to generate the random noise vectors, we specified how close the particles searched around the initial pose. We chose a standard deviation of 0.1, which meant that most particles searched around the initial pose value.

With 300 iterations and 16 or more particles, the PSO algorithm converged on the manually annotated pose if the initial CNN prediction was sufficiently accurate. After refining wing pose using PSO, we smoothed the data using an extended Kalman filter as described in detail elsewhere[6].

Representing wing motion in the strokeplane reference frame

Although mathematically preferable, quaternions do not provide an intuitive sense of wing motion. Instead, we defined a set of Tait-Bryan angles relative to a thorax-fixed reference frame. Prior studies of free flight kinematics[3, 7, 8] have used this approach to define a reference frame at a fixed angle relative to the longitudinal body axis. In case of tethered flight, however, the animal often deflects its head and abdomen for prolonged periods as it attempts to fictively steer. This is problematic, as the body’s longitudinal axis does not necessarily align with the symmetry plane of the thorax under these conditions. Using the pose vector of only the thorax is problematic because its quaternion is not sufficiently accurate. Instead of determining a reference axis based on the body, we defined our wing reference frame using the positions of the left and right wing roots. Because the fly is tethered, the wing roots move in a sable, ’c’-shaped trajectory around the thorax (Fig. 2b). Performing a principle components analysis (PCA) on the left and right root traces provides a convenient means of defining a strokeplane reference frame (SRF). The first three principal components correspond to the y-, x-, and z-axis of the SRF, respectively. Although the axes defined in this way are orthogonal, they are not fixed consistently in space. To create a reference frame that has the x-axis pointing forward and the z-axis pointing upward, we used the thorax pose vector to establish directionality. Using the normalized axes of the SRF, it is possible to compute the quaternion of the SRF, q_SRF. The left and right wing quaternions relative to the SRF may be obtained by multiplying the left and right wing quaternions with the conjugate, . The left and right wing root traces can then be expressed relative to the SRF by subtracting the origin location of the SRF and multiplying by the rotation matrix of

The mechanics of the wing hinge are complex and the root of the wing is not realistically approximated as a single point. This is why we choose not to define any constraints between the wing and thorax in Flynet. However, to simplify the analysis of wing motion and aerodynamics, we assumed the hinge to act as a ball joint at a fixed position on the thorax with 3 degrees of freedom. From an aerodynamics perspective, this assumption is not likely to change computed or measured aerodynamic forces substantially, as the arm of the wing root is relatively small and does not increase wing velocity significantly.

The three Tait-Bryan angles that specify wing orientation are: 1) the stroke angle, ϕ, that describes the angle between the y-axis of the SRF and the projection of the wingtip on the x-y plane, 2) the deviation angle, θ, that is the angle between the wingtip and the x-y plane, and 3) the wing pitch angle, η, that specifies the orientation of the leading edge panel of the wing with respect to the z-axis (Fig. 2c). Wing shape is described by ξ and remains unaltered from the Flynet definition. The wing kinematic angles of the right wing are similarly defined such that the signs are symmetric with respect to the left wing.

With the wing kinematic angles defined, we parsed the left and right wing motion into distinct wingbeats. A wingbeat was defined as the time between two subsequent dorsal stroke reversals. We computed the derivative of the stroke angle, ∂ϕ/∂t, and used the condition ϕ > 0 and ∂ϕ/∂ > 0 to find the dorsal stroke reversal times of the left and right wings in each video sequence. The dorsal reversal time of the left wing does not necessarily align with the right wing, and can be up to 3 frames apart. To remedy these small phases differences, we computed an average time point between the dorsal reversals of the two wings and subsequently rounded to the closest time point. For each high-speed video video sequence, the average dorsal stroke reversal times were determined and used to parse out individual wingbeats, which we indexed so as to preserve the instantaneous frequency and the numerical position within each sequence.

Encoding wing motion with Legendre polynomials

Depending on the activation level of the indirect power and tension muscles, wingbeat frequencies can vary between 150 Hz and 250 Hz[9]. This variation in wingbeat duration makes it difficult to compare wing kinematics between different high-speed video sequences. A wingbeat frequency of 200 Hz corresponds to about 75 high-speed video frames when captured at 15,000 frames per second. As there are four kinematic angles per wing, there are thus 300 data points representing wing motion during one wingbeat. In order to reduce the number of data points per wingbeat and allow for comparison across wingbeats, we used Legendre polynomials to parameterize wing kinematic traces. Legendre polynomials are well suited for encoding wing kinematics, as they can fit non-periodic traces and asymmetric waveforms. Fourier coefficients, which are often used for fitting time series, require periodic boundary conditions and are biased to symmetric wave forms when a low number of harmonics is used. Similar to Fourier series, Legendre polynomials form a complete and orthogonal system. An easy way to generate a Legendre polynomial is using Rodrigues’formula: where n is the order of the Legendre polynomial and x ranges between -1 and 1. A Legendre basis of order n is formed by all polynomials from order 0 up to order n. By specifying sample points in the range x = [−1, 1] and computing the values for each polynomial in the Legendre basis, a Vandermonde matrix is created. The Vandermonde matrix, X, can be used in a least-squares fit of a wing kinematic trace, Y : where β is a vector of n+1 coefficients corresponding to the order of the Legendre basis. With sufficient coefficients it is possible to fit any trace throughout a wingbeat with very little error. When using too many coefficients, however, the the polynomial fit may exhibit high frequency oscillations at the beginning and end of the waveform known as the Runge phenomenon. This may be remedied by lowering the order of the Legendre basis and imposing boundary conditions at the start and end of the wingbeat. Boundary conditions on polynomial fits can be imposed using restricted least-squares[10]. By this method, the restricted least-squares fit, β^*, makes use of the the unrestricted least squares fit, β, a restriction matrix, R, and restriction vector r:

To reduce the Runge phenomenon, we chose the restriction matrix and vectors such that the transition between subsequent wingbeats was continuous up to the fourth derivative. The j ^th derivative of a Legendre polynomial is computed from:

The 1^st, 2^nd, 3^{r d}, and 4^th derivatives of the actual wing kinematic trace were computed over a 9 frame time window centered around the start and end time points of the wingbeat. By requiring that the Legendre fit matches the actual values of the wing kinematic trace and the derivatives, the relationship between the least squares fit, β^*, the restriction matrix, R, and the restriction vector, r is:

The restriction matrix contains the Legendre basis and the derivative Legendre bases up to the 4^th derivative for x = −1 and x = 1 (start and end of the wingbeat). Matching values of the wing kinematic trace and the 1^st, 2^nd, 3^{r d}, and 4^th derivatives are contained in the restriction vector r. Inserting R and r in the restricted least-squares solution of equation 7 provides the Legendre coefficients that make the transitions between the previous and next wingbeat continuous up to the 4^th derivative.

After finding the restricted least-squares solution, we fit the waveforms using higher order Legendre polynomials, making sure to stay below the number of coefficients that induced Runge oscillations. Using higher order polynomials is preferable as it allows for a more accurate fit. For each wing kinematic angle, we tested various number of Legendre polynomials. The stroke angle, ϕ, could be fitted accurately with 16 polynomials. The deviation and deformation angles (θ and ξ) required 20 polynomials each. The wing pitch angle, η, required 24 polynomials. Thus, a we used a total of 80 Legendre coefficients to accurately describes the motion of a single wing during each wingbeat.

Convolutional neural network predicting wing motion from muscle activity

We created a CNN to perform non-linear regression between muscle activity, as measured by GCaMP fluorescence[11], and wing motion, as extracted from the high speed video data using Flynet and PSO. A graphical representation of this code (created using the the model plotting utility in Keras 3.0) is shown in Appendix C; the full code is available as described in the Data Availability Section of the main manuscript.

Before the CNN could be trained, we first normalized the input and output data. To rescale muscle activity traces, the mean, μ, and standard deviation, σ, were computed for the periods when the fly was flying. A boolean value that specified whether the fly was flying (as determined from wingbeat frequency) was saved to an hdf5 file with a timestamp during each experiment, and could be subsequently used to identify flight bouts and exclude sequences when the animal stopped flying. For all steering muscles, except the b₂ and iii₁ muscles, the activity was rescaled such that a value of 0 corresponds to μ − 2σ and 1 to μ + 2σ. The b₂ and iii₁ muscles are typically quiescent during flight, and in some some experiments were not active at all. To accommodate sequences when these muscles were quiescent, the values of b₂ and iii₁ were not rescaled if the standard deviation was below 0.01. If the standard deviation was above this threshold, the b₂ and iii₁ traces were rescaled such that a value of 0 corresponded to μ and a value of 1 to μ + 3σ. Although we tracked wing kinematics for both wings, steering muscle activity was only recorded from the left side of the thorax; thus, we only used the kinematics data from the left wing in our analysis. The 80 values of the Legendre coefficients are in units of radians, which we normalized by dividing by π. Wingbeat frequency was rescaled by subtracting 150 Hz from the raw values and subsequently dividing by 100 Hz, such that a value of 0 corresponds to 150 Hz and a value of 1 to 250 Hz.

Although neural networks are capable of handling outliers, too many outlying data points could result in decreased performance. We excluded wingbeats from the dataset under any of the following conditions: 1) the normalized muscle activity was lower than − 0.5 or larger than 1.5, 2) the normalized wingbeat frequency was lower than 0 or larger than 1, and 3) if any of the normalized Legendre coefficients were not inside specified ranges: [−1 /3 <= C_θ <= 1 /3], [−2 3 <= C_η <= 2 /3], [−2 /3 <= C_ϕ <= 2 /3], and [−1 3 <= C_ξ <= 1 /3]. After removing all outliers, the final dataset consisted of 72,219 wingbeats. For the validation set, we selected the first 30 wingbeats of each high-speed video sequence (10,868 total). The remaining wingbeats in each sequence (61,351 total) were used for the training set.

In the first layer of our CNN (Fig. 3), the network processes the GCaMP7f fluorescence signal via a number of convolutional kernels over a fixed time window. After evaluating several different values, we found that a time window length of 9 wingbeats (approximately 45 ms) worked well. A shorter time window did not contain sufficient information and prediction accuracy was worse as a consequence. Longer time windows improved the training error, but the validation error tended to be higher, especially for very long windows (>50 wingbeats).

Our CNN architecture consisted of 2 convolutional layers and 2 dense layers (Fig. 3a). Input of the network consisted of a 13 × 9 matrix of normalized muscle activity and normalized frequency of 9 subsequent wingbeats. The output of the network is a prediction of the 80 normalized Legendre coefficients corresponding to the first wingbeat in a 9 wingbeat time window. Before the input data was fed into the first convolutional layer, Gaussian noise with a standard deviation of 0.05 was added to the input matrix, which helps the CNN to generalize on the data and makes the network more robust to noise in the fluorescence recordings. The first convolutional layer of the CNN consisted of 64 kernels with a 1× 9 kernel window, a1 × 9 stride and SELU activation (Fig. 3a). A total of 256 kernels were used for the second layer, with a 13 × 1 kernel window and stride, and SELU activation. The output of the second layer was a vector with 256 elements. A fully connected dense layer of 1024 virtual neurons takes in the output of the second layer and applies a SELU activation for each virtual neuron. The last layer of the CNN has 80 neurons with a linear activation function, corresponding to the 80 normalized Legendre coefficients.

The CNN was trained for a total of 1000 epochs with a batch size of 100 samples. After every epoch, the network was evaluated using the validation set. The learning rate was set as 10⁻⁴ with a decay of 10⁻⁷ per epoch with MSE as a loss function. After the first 200 epochs, the validation error stabilized at approximately 10^−3.7, whereas the training error continued to decline (below 10^−3.33) after 1000 epochs (Fig. 3d). After training, the prediction performance of the network was tested on all videos in the dataset. Examples of the CNN prediction performance are given in Fig. 3c and Extended Data Fig. 3. For most sequences, the wing motion predicted from muscle activity was within ±2^° of the tracked wing motion for each angle. This accuracy was quite remarkable, given the complex waveforms of the wing kinematic angles, as well as the sparse and filtered nature of the input data.

Virtual experiments using the trained CNN

After developing a CNN to predict wing motion from observed muscle activity, we used it as an in-silico model of the hinge to conduct virtual experiments. The wing kinematics predicted by baseline muscle activity was determined by inspecting the the dataset for sequences that contained no significant changes in muscle activity or wing motion, and subsequently averaging data over those sequences. The baseline values of activity for each muscle and wingbeat frequency, kept constant over the 9 wingbeats, were fed into the trained CNN to yield the baseline pattern of wing motion wingbeat (gray traces, Fig. 4). A naive way of using the CNN to predict the action of individual muscles would be to simply vary the activity of each muscle independently, while keeping that of all other muscles constant constant. However, we observed strong correlations among the activity patterns of many of the muscles within our dataset, determined using RANSAC[12] (Extended Data Fig. 4, Extended Data Table 1). Thus, if we performed virtual experiments by independently changing the activity of each muscle, the input patterns would likely reside outside the data subspace on which the CNN was trained, making the network predictions unreliable. To make more reliable predictions of wing motion from the trained CNN, we deliberately interrogated the CNN using input patterns that would not deviate substantially from those used to train the network by making use of the measured muscle correlations (Extended Data Fig. 4). According to this strategy, we used the trained CNN to probe the influence of each muscle by changing its input value to the maximum normalized value, while simultaneously changing the activity of the other 11 muscles according to the measured correlations.

Evaluating the function of individual sclerites using latent variable analysis

The CNN that we constructed and trained to predict wing kinematics from muscle activity deliberately mixed input from all 12 control muscles, making it difficult to draw conclusions on the mechanical function of individal wing sclerites. As a complementary approach, we trained a CNN with an encoder-decoder architecture to learn how the activity of all the muscles that attached to a particular sclerite were correlated to changes in wing motion (Fig. 6; Extended Data Fig. 8). The encoder block contained deliberate bottlenecks that forced the network to learn how to represent the input data by a small set of latent variables representing the state of each sclerite and wingbeat frequency. Training the encoder-decoder with a bottleneck is roughly equivalent to performing a non-linear PCA on the dataset[13]. A graphical representation this code (created using the the model plotting utility in Keras 3.0) is shown in Appendix D; the full code is available as described in the Data Availability Section of the main manuscript.

We constructed the encoder-decoder with two decoder heads, one that predicts muscle activity from the latent variable space and another that predicts wing kinematics (Extended Data Fig. 8). Input to the encoder-decoder was split into 5 streams, grouping the activity of muscles that insert on each of the 4 wing sclerites with a separate stream for wingbeat frequency. Each stream uses two convolutional layers, as in our original CNN model (Fig. 3a), but followed by a dense layer of 512 neurons and a final layer of one neuron with linear activation that corresponded to one of the five latent variables. After the latent space, the network is split into two streams: one stream uses a decoder to predict the input to the network (i.e. the muscle activity), the second stream uses two dense layers to predict the Legendre coefficients of the wing motion. We placed a back-propagation stop between the latent space and the muscle activity decoder layers, such that the latent space is only trained on correlations with wing motion. The muscle activity decoder layers are still functional, as the decoder predicts the correlation between the latent space and muscle activity, but without the noise of the input to the network.

We trained the network on the muscle activity and wing motion dataset using 1000 epochs, a batch size of 100, a learning rate of 10⁻⁴, and a decay rate of 10⁻⁷. The network predicts two vectors: the latent space (5 × 1) and Legendre coefficients (80 × 1), and one matrix: muscle activity and wingbeat frequency over 9 wingbeats (13 × 9). It is important to note that the prediction error of the wing kinematics by the sclerite function encoder-decoder is expected to be worse than our original CNN, as the bottleneck required to create the latent variable space greatly restricts the amount of information that can be used to learn the correlations. In Figure 6 and Extended Data Fig. 8, only the decoder sections of the network were used to predict muscle activity and wing motion for each latent variable, which we varied systematically between − 3σ and + 3σ in 9 steps while keeping the other latent variables constant at 0. The muscle activity and wing motion were then predicted from the latent input by the two decoders.

Dynamically-scaled flapping wing experiments

In the last three decades, our laboratory has developed several versions of a dynamically-scaled flapping wing robot[14–18]; however, none of these had the capability of creating the chord-wise deformations that we observed and tracked in our high-speed videos. The robotic wings we designed for this study were actuated by one stepper motor and two servo motors each (Extended Data Fig. 5a). Stroke angle was controlled by a stepper motor (10 kHz clock) via two gears with a 1:3 gear ratio. The position of the stepper motor was controlled via micro-stepping, which permitted fine motor control (7.8·10⁻³ degrees per microstep). Two magnets were positioned on the gear so that a Hall-effect sensor was activated when the wing moved out of bounds. Besides protecting the wing, the Hall-effect sensor was also used to home the stepper motor. During an experiment, a Teensy 3.2 microcontroller tracked the number of microsteps travelled relative to the home position. To account for motor slip, which sometimes occurred due to large torques, the steppers were homed after each experiment. During an experiment, the stepper’s position, ϕ, and velocity, , were controlled via a feed-forward process. The stepper made micro-steps on a 10 kHz clock; position and velocity set points were updated at 200 Hz. The custom instrumentation software computed the required stepping frequency and direction to reach the specified position and velocity values after each update.

The deviation and wing pitch angles were controlled by two servos (HiTec D951TW). The deviation angle servo operated via two gears (1:1 gear ratio), and the wing was directly attached by the wing pitch servo. The position of the servos was specified by a pulse-width modulation (PWM) signal at 50 Hz. The deviation servo could move between −45^° and +45^° and the wing pitch servo could move between +90^° and −90^°. During an experiment, the position of the servo was updated at 50 Hz in a feed-forward control loop.

We measured forces and torques on the wing using a 6 degree-of-freedom force and torque (FT) sensor (ATI Nano 17), mounted on the rotation axis of the wing pitch servo. Custom-machined aluminium mounts coupled the base of the FT sensor to the servo axis and the wing. During each experiment, six 16-bit unsigned integers corresponding to all forces and torques (F_x, F_y, F_z, T_x, T_y, T_z) were sampled at 200 Hz. The robot was submerged in an acrylic tank (2.4× 1× 1.2 m) filled with mineral oil (Chevron Superla white oil), with a kinematic viscosity of 115 · 10⁻⁶ m²·s⁻¹ and a density of 880 kg ·m⁻³ at 22 C^°[19]. The oil level in the robofly tank was approximately 1 m. The robot had the wingtip radius of 31 cm and was submerged at a depth of 50 cm. A previous study that systematically mapped the effects of the tank boundaries on lift an drag determined that the system approximated the forces generated within an infinite volume, provided the wing remains further than 12 cm from any boundary (top, bottom, or side)[14].

Besides the three Tait-Bryan angles describing wing orientation, Flynet tracked a fourth angle describing wing shape, the deformation angle ξ. To implement this fourth angle on the robot, the wing was composed of four panels connected by three hinge lines (Exended Data Fig. 6a). The four panels were cut out of an acrylic sheet (2.75 mm thickness). Each hinge consisted of a 2 mm steel rod core, surrounded by an acrylic tube with inner and outer diameters of 2 mm and 4 mm, respectively. The acrylic tube was cut into sections of 20 mm and these sections were glued in an alternating pattern to two adjacent panels. The rotation angle between two subsequent panels was controlled by a micro servo (HiTec HS-7115TH), which was screwed onto one panel and connected to the next panel via a 1 mm metal rod that was coupled to the servo arm. Three micro servos were used to deform the wing. Following the assumption that the wing bending angle was uniformly distributed over the three hinge lines, the deformation angle ξ was divided by 3 to obtain the rotation angle per hinge line.

Dynamic scaling was ensured by matching the Reynolds number of the robotic wing and a real fly, which requires: where n _fly is the fly’s wingbeat frequency, n_robo the wingbeat frequency of the robot, R _fly the fly’s wing length, R_robo the wing length of the robot, v_air the kinematic viscosity of air and v_robo the kinematic viscosity of the mineral oil. Rewriting the Reynolds number equality yields:

By entering the values for a typical fly (R _fly = 2.7mm, v_air = 15.7 ·10⁻⁶m²s⁻¹ at 25^°C, n _fly = 200Hz) and the RoboFly parameters (R_robo = 310mm, v_robo = 115· 10⁻⁶m²s⁻¹ at 22^°C), the required flapping frequency for the robot is 0.11Hz. The force scaling factor may be written as[6]: with relevant density values of ρ_air = 1.18kgm⁻³ and ρ_robo = 880kgm⁻³. The torque scaling factor is the product of force and moment arm:

Any experiment using the robotic wing required an input file that specified the following angles: ϕ,, θ, η and ξ at intervals of 5 ms. A C++ program read the five angles and sent them to the Teensy microcontroller via a USB-cable at 200 Hz (RawHID protocol). The Teensy microcontroller converted the angles to PWM commands for the servos and step&direction commands for the stepper motors. At the same time, the force and torque data were sent from the Teensy micro-controller to the C++ program, which logged the measurements.

At the start of an experiment, the wing was moved to the home position (ϕ = 0, θ = 0, η = 0, ξ = 0). In order to prevent rapid acceleration of the wing, the wing kinematic angles during the first wingbeat of an experiment were multiplied with the following function: where T corresponds to the wingbeat period. Similarly, the last wingbeat of the experiment ends at the home position and the wing kinematic angles are multiplied by:

When wings start flapping in the oil, it takes approximately two to three wingbeats before the wake is fully developed[15]. In order to exclude these start-up effects, a wing kinematic pattern was repeated for 9 wingbeats and only wingbeats 4 through 8 were used for analysis. After each experiment, the data were converted from 16-bit unsigned integers into actual force and torque vectors and smoothed using a linear Kalman filter. Besides aerodynamic forces, the FT-sensor also measured inertial and gravitational components. The low flapping frequency of the robot assures that the inertial forces may be ignored; however, the gravitational forces are significant. To remove the gravitational force from the analysis, every wing kinematic pattern was replayed at a 5x slower frequency, which reduces aerodynamic forces by a factor of 25. The forces and moment measured during the slow frequency, which approximate the gravity component, were subtracted from the fast frequency data after interpolation.

The aerodynamic and gravitational forces were measured by the FT sensor in the wing reference frame. For subsequent analysis it is useful to translate the forces and torques to the stroke reference frame (SRF), which was accomplished by the following Tait-Bryan operations: and where F^w and F^SRF are the forces in the two different reference frames, and F^w and R_η, R_θ, R_η, correspond to the rotation matrices. More details on thse operations are provided elsewhere[6].

Computing inertial forces via the Newton-Euler equations

Besides aerodynamic forces, wing inertia plays a significant role in fly flight. Although the wing mass is only ∼0.2% of the body mass[9], the angular velocity and acceleration of wing motion is very high. We estimated the inertial forces and torques of a rotating rigid body using the Newton-Euler equations: where m_w corresponds to the wing mass, c_w is the position of the center of gravity on the wing, I_w is the inertia tensor of the wing, a is the linear acceleration of the wing, is the angular acceleration, and ω the angular velocity, with all values relative to the wing reference frame. We ignored the effects of the linear acceleration of the wing because they are very small compared to the angular acceleration. The first right-hand term in equation 18 corresponds to the inertial forces and torques due to wing acceleration, F _acc and T _acc. Inertial effects dependent on angular velocity, such as the Coriolis and the centrifugal forces, are captured by the second right-hand term and are collectively referred to as F_angvel and T_angvel. We assumed that the inertial effects of wing deformation are small, and thus calculated inertial forces and torques based on the assumption that the wing is a rigid, flat plate.

The angular velocity and acceleration had to be computed from the Legendre polynomials describing wing motion. We computed the emporal derivatives of the Legendre polynomials using equation 8. However, the temporal derivatives of the wing kinematic angles do not correspond to the angular velocity and acceleration. Angular velocity in the wing reference frame is given by: where R_θ and R_η correspond to the rotation matrices of the Tait-Bryan operations for θ and η, respectively. Angular acceleration is given by: with being the temporal derivative of the rotation matrix.

To compute the inertial forces and torques, the mass, center of gravity, and the inertia tensor of the wing must be determined. These parameters were computed using the scaled 3D model we created for Flynet, with a wing length of 2.7mm, an estimated cuticle density of 1200kg·m⁻³, and a wing thickness of 5.4μm[20]. The inertial parameters of the (left) wing are given by:

With the wing inertia parameters, angular velocity, and acceleration defined, it is possible to compute the inertial forces and torques for any arbitrary pattern of wing motion pattern. Extended Data Fig. 5b shows the aerodynamic and inertial forces and torques in the SRF for the baseline wingbeat. The aerodynamic forces were measured on the robot and rescaled to the appropriate units for an actual fly using equations 12 and 13. To make interpretation easier, the forces and torques have been non-dimensionalized by dividing the forces by the body weight, and the torques by product of body weight and wing length. Using the estimated cuticle density and the scaled body components of the 3D fly model, the body mass was found to be 1.16 ·10⁻⁶ kg. Inspection of the traces in Extended Data Fig. 5b indicate that the inertial forces and torques in wing motion are substantial, as found in prior studies[21].

With the methods to measure aerodynamic forces and compute inertial forces established, it is possible to evaluate the forces generated by changes in steering muscle activity. For each steering muscle, seven wing kinematic patterns were tested on the robotic wing. The seven wing kinematic patterns were found by sampling the muscle activity patterns on a line between the baseline muscle activity and the maximum muscle activity of a selected muscle (Fig. 4), and subsequently predicting the corresponding wing motion using the trained CNN. The workflow to obtain the control forces and torques of muscle activity consisted of measuring the force and torque traces of 84 wing kinematic patterns on the robot, performing the gravity subtraction, smoothing the traces with a Kalman filter, converting traces to the SRF, computing the median force and torque traces from the measured wingbeats (F _aero and T _aero), estimating the inertial components (F _acc, T _acc,F _angvel, T _angvel), and finally, summing the inertial and aerodynamic forces and torques (F_total and T _total). A time history of the individual force and torque components are plotted over a baseline wingbeat in Extended Data Fig. 5b. The corresponding instantaneous total force vectors are shown superimposed over the 3D kinematics in Extended Data Fig. 5c, with comparable plots for the kinematics in each virtual muscle activation experiment (Fig. 4) shown in Extended Data Fig 6. Further details on the construction and operation of the robot and the analysis of aerodynamic forces and torques are provided elsewhere[6]

Simulating arbitrary free flight maneuvers with Model Predictive Control

Using our CNN, the virtual muscle activation experiments, and robotic wing experiments provides us with the components necessary to construct a state-space simulation of a flying fly using Model Predictive Control (MPC)[22], in which the virtual fly regulates its flight trajectory via changes in the time course of steering muscle activation. By specifying an initial and goal state and a time period to achieve the goal state, a MPC controller finds the optimal trajectory given a cost function and constraints. The explicit discrete time-invariant state-space system is governed by the following equations: where x(k) and x (k+ 1) are the state vectors at times k and k+ 1 respectively, u(k) is the control vector, A (k) is the system matrix, B (k) is the control matrix, C (k) the output matrix, D (k) the feed-forward matrix, and y (k) the output vector. In the case of fly flight, there is no feed-forward process and the output of the state-space system is the state vector. This means that matrices C (k) and D (k) do not need to be defined (Extended Data Fig. 7a).

The state vector describes the orientation, position, velocity, and angular velocity of the fly’s body: where v and ω are the linear and angular velocity of the body in the SRF, respectively, q is the quaternion describing SRF orientation relative to the inertial reference frame, and p is the position of the SRF in the inertial reference frame. For simplicity, we assume that the head and abdomen are stationary relative to the thorax. The temporal derivative of the state vector, , is required for updating the state for each time step: where a and are the linear and angular acceleration of the body in the SRF, respectively.

At each time step, the state vector is multiplied with the system matrix, A (k), to compute the temporal derivative of the state vector. The time step, Δt, corresponds to the wingbeat period, which is assumed to be constant at 1/f = 1/200 = 0.005 s. Although flies regulate wingbeat frequency during flight, it typically takes 20 wingbeats to increase or decrease the wingbeat frequency to a new level[3]. As our model was developed for simulating rapid maneuvers of only 10 wingbeats, frequency was held constant in our simulations.

The equations of motion of the system matrix include the following flight forces and torques: body weight, body inertia, body aerodynamics, and the aerodynamic and inertial damping of the wings. The mass, center of gravity, and inertia tensor of the body were determined from the scaled 3D body model (head, thorax, abdomen) for the average fly in the dataset. For the wing mass, the cuticle density was assumed to be 1200 kg · m⁻³, and the body mass estimated as m_b = 1.16 · 10⁻⁶ kg. The center of gravity of the body in the SRF was estimated as cb=[0.04 0 −0.24] ^T mm. The inertia tensor of the body is given by:

The inertial forces and torques on the body can be computed with the Newton-Euler equations, evaluated at the center of gravity of the body:

Obtaining the rotational accelerations of the body requires solving the inverse problem:

As the state is evaluated at the center of gravity, there are no torques due to the fly’s weight. Although the aerodynamics of the wings generate larger forces, the body of the fly experiences drag during flight. The combined 3D shape of the head, thorax, abdomen and legs is complex, but we chose to model the body conservatively as sphere with a radius of 1 mm. While the actual body drag is likely to be lower, the simplicity of the drag model makes it easier to implement in the state-space system. The drag on the sphere can be calculated as: with the drag coefficient as C_D = 0.5 and sgn as the sign operator. Because the center of pressure of the sphere is assumed to be at the center of gravity of the body, there are no torques due to body drag.

Because the aerodynamic and inertial damping terms due to the motion of the flapping wings are dependent on body rotational velocity, these components are added to the system matrix. The equations of motion to compute the relevant linear and rotational accelerations can be written as: where F_dampA and T_dampA are the aerodynamic damping terms and F_dampI, and T_dampI are the inertial damping terms. F_G is the gravity component, expressed in the SRF[6].

The control inputs that are available to the fly are the muscle activity patterns of the left and right steering muscles: u^L and u^R. The sum of aerodynamic and inertial forces and torques generated via changes in these inputs can be computed as:

Adding these control forces and torques to equation 31 gives the complete equations of motion for body rotational acceleration:

The computed accelerations of equation 33 are subsequently be used to update the state vector. As the body accelerations and velocities are both in the SRF, the update of the body velocity is relatively simple:

The body quaternion and position are in the inertial reference frame and the angular velocity update must be transformed from the SRF to the inertial reference frame. Using quaternion multiplication, the quaternion update can be written as:

After computing q_k+1, the quaternion must be normalized such that ∥q∥ = 1. The position update is given by: where is computed using q_k.

The aerodynamic and inertial forces for the maximum muscle activity patterns were measured and computed in a stationary reference frame. During free flight, the body translates and rotates, which has an effect on both the inertial and aerodynamic forces on the wing[23–25]. We used a quasi-steady aerodynamic model and the Newton-Euler equations to compute the relevant aerodynamic and inertial damping coefficients for fruit fly flight. The effects of inertial damping, can be computed given the angular velocity of the wing: where the last term adds the body angular velocity, converted to the wing reference frame. By inserting the redefined angular velocities of the left and right wings into the Newton-Euler equations, it is possible to compute the inertial forces and torques given a constant body angular velocity.

The wingbeat-averaged inertial forces and torques that are generated by body rotation show a linear trend with angular velocity. Translational body velocity does not cause any changes in inertial force or torque. The slopes of the linear trends form the inertial damping matrix:

Rotational body velocity generates a torque in opposite direction for ω_x and ω_z, but a torque in the same direction for ω_y.The positive damping torque for pitch rotation means that flight is unstable around the pitch axis. There is a strong coupling between the roll rotation and yaw torque, and vice versa. Aerodynamic damping was computed using a quasi-steady model, the details of which are described elsewhere[6].

The angular velocity of the body can be added to the angular velocity of the wing using equation 37. Inclusion of the translational body velocity into the wing’s angular velocity term is difficult, however. We therefore implemented the translational and rotational quasi-steady terms using a blade-element formulation, in which the wing is partitioned into spanwise sections and the air velocity vector is computed on each section[26]. Simulating the baseline wing kinematics on the left and right wings, with the relevant body velocities yields the damping matrix for body translational velocity: and the damping matrix for body angular velocity:

The aerodynamic damping coefficients are approximately an order of magnitude larger than the inertial damping coefficients. Both dF_Ax/dv_x and dF_Az /dv_z have a negative sign, which means that there is a force opposite to the direction of body motion. In cases of dF_Ay /dv_y, dF_Az/dv_x, and dF_Ax /dv_z, the signs are positive, however; which indicates a force in the direction of body motion. Whether the damping coefficient is negative or positive depends on the effects of body velocity on parameters such as angle-of-attack, instantaneous air velocity, and wing orientation. All the aerodynamic damping torques, dT_A /dω, have a negative sign. This means that for the baseline wing motion, flight is stable to perturbations in angular velocity, a conclusion consistent with a prior study[23].Matrices dFT_I dω, dFT_A /dv, and dFT_A /dω can be converted into a single damping matrix. Multiplication of the damping matrix with the linear and angular velocity vectors of the body yields the damping forces and torques:

A central feature of our model is that the control inputs available to the fly are the muscle activity patterns of the left and right steering muscles: u_L, u_R. Using the measurements from the robotic wing, it is possible to compute the Jacobians (dFT/du_L, dFT/du_R) of the aerodynamic force and torque production for maximum muscle activity pattern of each muscle. Steering muscle activity is ordered as follows in the control vectors: and

A full table of numerical values of the calculated Jacobians for the aerodynamic force and torque production of the left and right for each steering muscles are provided elsewhere[6]. Control forces and torques can be computed by multiplying the Jacobian with the control vector:

The muscle activity correlation analysis in Extended Data Fig. 4 shows that not all combinations of muscle activity are possible. It is, therefore, necessary to impose constraints on the left and right control vectors. Specifically, this constraint is that muscle activity has to lie on a 12-D surface representing the muscle correlation matrix (Extended Data Table 1). A mathematical formulation for the 12-D surface constraint is: where C is the 12 ×12 correlation matrix (Extended Data Table 1), and u_base is the baseline muscle activity pattern.

MPC optimizes trajectories in state-space, over a finite time window, for a given cost function[22]. The state-space system is used to predict a state trajectory over a finite time horizon. At each time step, the controller determines the optimal control input that generates one or more trajectories that minimize the cost function, while being subject to system constraints. This process is repeated for every time step, shifting the horizon forward. MPC is thus also known as receding horizon control.

We used the the Python-package do-mpc[27] to implement our MPC simulation. For the free flight simulations, the controller uses discrete time and a non-linear model. The duration of all free flight simulations was 10 wingbeats, and the MPC horizon was set to 10 wingbeats as well. This means that the goal state is always visible to the controller. At the center of MPC is the objective function, which is defined in do-mpc as: with N being the number of time steps of the problem, l (x_k, u_k)is the Lagrange term, is the r-term, and m (x_N)is the Meyer term. The Lagrange term provides the option to provide a reference-state trajectory and add a penalty for any deviation. Rapid changes in control input are penalized using the r-term: Δu_k = u_k−u_k−1, and R is a weighting matrix. The Meyer term evaluates how close the specified goal state is to the final state, x_N.

In case of the free flight simulations, the MPC objective function is defined as follows. First, the initial state, x₀, and the goal state, x_n, are specified. Subsequently, the Lagrange and Meyer terms are both defined as: with S_x as a scaling matrix. By setting the diagonal values of the S_x and R matrix, the user can tune the importance of certain aspects of the objective function. The Lagrange and Meyer terms we used were always the same, as we did not want to specify any reference trajectory for the Lagrange term. In this way, the objective function allowed the MPC controller to explore any trajectory that ends at the goal state. By setting a diagonal element to zero in the S_x matrix, the MPC controller will not optimize the trajectories for this parameter. If a diagonal element in S_x is set to a high value, the MPC controller will prioritize this element. For example, by setting the weight for the goal, v_y, to 1, the controller will punish any deviation from this goal heavily, and allows the user to enforce straight flight. In a similar way, the user can set the diagonal values of the R matrix, and determine the behavior of the steering muscles during the simulation. For the free flight simulations, we set all the diagonal values of R to 1.

Besides tuning the cost function, the bounds of the state and control vectors must be specified. In case of the control vectors, all steering muscle activity was bounded between 0 and 1. In the state vectors, linear velocity was bounded by − 10⁴ and 10⁴ mms⁻¹, angular velocity was bounded by -1000 and 1000 rads⁻¹, quaternion values were bounded by -1 and 1, position was bounded by −10⁴ and 10⁴ mm, and the body quaternion was constrained to be a unit quaternion, ∥q∥= 1. Do-mpc allows the user to specify the tolerance by which a state or control vector can deviate before a constraint will be enforced. For the free flight simulations, we set the constraint tolerance on the control vector, c_tol, to 0.001.

Extended Data Fig. 7a shows the workflow of the MPC controller. The user only needs to specify the initial and goal states, time in which the fly should reach the goal state, and the weight matrix of the objective function, S_x. With these parameters, the MPC controller will try to find the optimal state trajectory and required control inputs, that minimizes the objective function while satisfying the constraints. Further details on the design and implementation of our MPC controller may be found elsewhere[6]. Our software is available in the form of a Jupyter notebook as described in the Code Availability Statement associated with our manuscript.