Physics-based simulations of aerial attacks by peregrine falcons reveal that stooping at high speed maximizes catch success against agile prey

The peregrine falcon Falco peregrinus is renowned for attacking its prey from high altitude in a fast controlled dive called a stoop. Many other raptors employ a similar mode of attack, but the functional benefits of stooping remain obscure. Here we investigate whether, when, and why stooping promotes catch success, using a three-dimensional, agent-based modeling approach to simulate attacks of falcons on aerial prey. We simulate avian flapping and gliding flight using an analytical quasi-steady model of the aerodynamic forces and moments, parametrized by empirical measurements of flight morphology. The model-birds’ flight control inputs are commanded by their guidance system, comprising a phenomenological model of its vision, guidance, and control. To intercept its prey, model-falcons use the same guidance law as missiles (pure proportional navigation); this assumption is corroborated by empirical data on peregrine falcons hunting lures. We parametrically vary the falcon’s starting position relative to its prey, together with the feedback gain of its guidance loop, under differing assumptions regarding its errors and delay in vision and control, and for three different patterns of prey motion. We find that, when the prey maneuvers erratically, high-altitude stoops increase catch success compared to low-altitude attacks, but only if the falcon’s guidance law is appropriately tuned, and only given a high degree of precision in vision and control. Remarkably, the optimal tuning of the guidance law in our simulations coincides closely with what has been observed empirically in peregrines. High-altitude stoops are shown to be beneficial because their high airspeed enables production of higher aerodynamic forces for maneuvering, and facilitates higher roll agility as the wings are tucked, each of which is essential to catching maneuvering prey at realistic response delays.


Introduction
The stoop is a remarkable attack strategy used by peregrine falcons Falco peregrinus, and a range of other raptors [1][2][3][4]. It involves a steep, controlled dive in which the attacker strikes its prey at high-speed with a massive blow in mid-air [1]. The high momentum of the attacker places it at obvious risk of harm, especially when diving into flocks of birds [1] or when pulling out only meters from the ground [3]. Arguably, for the stoop to evolve as an habitual attack strategy, these risks must be outweighed by certain survival advantages, and stooping has therefore been proposed either to save energy [5], or to enhance catch success [6]. These hypothetical advantages remain unproven, however, because it is challenging to compare the success rates of different attack strategies empirically. Success rates are confounded by a variety of factors, including the experience [1,5] and reproductive status [7] of the attacker, the season of the attack [8], and the species of prey [6]. Even the seriousness of the attacker's behavior may be an important source of variation: falcons seem to not always focus on achieving a high success rate, and appear sometimes to be practising or playing with their prey [9]. Moreover, the outcome of the stoop is often difficult to observe due to its high speed [6].
There presumably exists a trade-off between different factors influencing catch success in a stoop. On the one hand, it has been proposed that the high speed of the attack provides an element of surprise, leaving little time for the prey to evade [5,10]. On the other hand, it is possible that the high speed of the attack decreases the precision of interception [2], and makes it harder for the attacker to follow the prey if it turns sharply [11]. Such trade-offs are difficult to investigate empirically, and we therefore turn to modeling and simulation. Because physical and physiological constraints influence catch success, we use an embodied cognition approach [12]. We investigate the success of different attack strategies by incorporating in a physicsbased simulation model the aerodynamics, flight mechanics, guidance, and control. Such detailed simulations have already proven useful in work on missile guidance: the increasing demand for better performing missiles forces the inclusion of the detailed dynamics of the missile and its target when comparing the effectiveness of different guidance systems [13]. The nonlinear nature of these dynamics restricts the use of analytic methods, such as linear-quadratic optimal control, and the effectiveness of different mechanisms must therefore be examined through parametric variation of the system between repeated simulations of interception. Here, we study the general intercept problem under the particular dynamics of flapping and gliding bird flight. Interception in this biological context differs from that of missiles in that gravity plays a pivotal role in determining the best attack strategy: in missiles, the speed and acceleration are so high that the effects of gravity are marginal, but in birds, the acceleration due to gravity dominates the dynamics.
The flight performance of the model-birds in our simulations depends on their flight morphology, and differs considerably between predator and prey. To maneuver, model-birds flap, course with its target (i.e. would hit its target if both continued flying at constant velocity thereafter). Under pure proportional navigation, the attacker turns at an angular rate proportional to the angular rate of the line-of-sight to target. Although this guidance law can be written in three-dimensional vector form, it is more intuitively explained in the two-dimensional case, for which: where γ denotes the bearing of the attacker's velocity vector and λ denotes the bearing of the line-of-sight to target, both measured in an inertial reference frame; the dot notation denotes the time derivative, and N is called the navigation constant. The numerical value of N determines the rate of convergence to a parallel navigation (CATD) course: in missiles, low values of N result in slow convergence, whilst high values can cause overshoot, leading to control instability [14]. Partly for these reasons, intermediate values of N between 3 and 5 are typical in most missile applications. The overall objective of our simulations is to identify the attack strategy that maximizes the catch success of the falcon for a given prey motion, a range of assumptions regarding the delay in the falcon's response, and the error in its vision and control (see Materials and methods).
Here, an attack strategy is defined as some particular combination of the predator's navigation constant N, and its initial vertical and horizontal distance from the prey. The optimization was conducted via parametric variation of the attack strategy, in combination with Generalized Additive Modeling (GAM), which we used to interpolate between the 10 6 randomly chosen attack strategies that we simulated in each optimization [19,20]. We hypothesise that stooping maximizes catch success, and that it does so as a direct consequence of the flight physics in our simulation model. We test this by asking whether a model-falcon's catch success is maximized by attacking from a high altitude, which couples into a high flight speed in our simulations. Remarkably, we find that the optimality of stooping depends not only on the motion of the prey, but also on the tuning of the underlying guidance law. Specifically, we show that stooping is only expected to evolve in conjunction with the same low values of the navigation constant N that have been identified empirically in peregrine falcons [15].

Starlings outmaneuver falcons in slow flight
Flight performance is expected to be a key determinant of catch success in a chase. Clearly, any prey species that can fly faster than a falcon will be able to outrun its attacker in straight flight. In practice, peregrine falcons fly much faster than starlings, and our aerodynamic model predicts that they hold a considerable speed advantage in both level flight (maximum speed: 29 versus 23 ms −1 ; Fig 3a) and vertical dives (terminal speed: 123 versus 52 ms −1 ; Fig 3b; see section K for a comparison between flight performance in the model and empirical measurements). Even so, escape is possible if the prey species can outmaneuver its attacker. For instance, if at a given flight speed the prey can produce a higher aerodynamic force relative to body weight than its attacker (i.e. produce a higher load factor), then it may escape by turning more tightly than its attacker in a smooth maneuver called a turning gambit [21,22]. Our aerodynamic model shows that starlings can indeed sustain higher load factors than peregrine falcons flying at the same speed (Fig 3c), and that although falcons can achieve even higher load factors by flying faster (Fig 3c), the net effect is such that a starling will always be able to turn on a tighter radius than a faster-flying falcon (Fig 3e). Similarly, if the prey can achieve a higher roll acceleration than the falcon, then it will be able to redirect its lift faster, and hence outmaneuver its attacker in a non-smooth jinking maneuver. Our aerodynamic model predicts that a starling can indeed produce a higher roll acceleration than a falcon flying at the same speed (Fig 3d). So great is a model-starling's advantage in this respect that a model-falcon can only be expected to match a model-starling's maximum roll acceleration by diving at close to terminal velocity (i.e. at close to its maximum speed). Hence, model-starlings may often escape model-falcons in our simulations, even though their maneuvers are not implemented as an evasive response to the falcon. In summary, a starling can always outmaneuver a falcon that is flying at a similar speed, but a falcon can always beat the load factor and roll acceleration of a starling by diving at a sufficiently higher speed. Whether this strategy enhances catch success will presumably depend on the flight pattern of the prey, and the complex ways in which the predator's flight speed, response delay, and errors in vision and control interact to affect its guidance. We explore the outcome of these complex interactions in the simulations presented below.

Stooping maximizes catch success if the prey maneuvers
The catch success of model-falcons was always maximized by entering a steep dive, but the optimal starting altitude varied greatly between the three different flight patterns of the prey (Table 1; see asterisked points in Fig 4). Catch success in attacks on straight-flying prey (Fig 4a)   Fig 3. Flight performance graphs in the flight simulator for the peregrine falcon (dark blue) and the common starling (light blue). The double arrows denote the direction of acceleration displayed in the graph. The starling is able to outmaneuver the falcon at a given airspeed, if there exists a region under the curve of the starling that is not overlapping with that of the falcon. (a) Level acceleration versus air speed: level flight with the requirement that lift equals weight. Dashed lines denote the speed wherein torque forces constrain the maximum acceleration (mechanical constraints). Top level flight speed is reached at the point where level acceleration is zero. (b) Vertical dive acceleration (including gravity) versus air speed. At the end of the dashed lines, flapping is substituted by gliding with retracted wings in order to maximize vertical acceleration. (c) Load factor versus air speed. The load factor is defined as lift divided by weight. The maximum load factor does not scale quadratically with forward speed due to constraints in torque forces [11]. Instead, wings are retracted optimally to increase maximum load. (d) Roll acceleration versus air speed. Roll acceleration determines the speed with which the bird can redirect its lift and is calculated by estimating the whole-body inertia around the roll-axis and the maximum net torque production [11]. (e) Turning radius is calculated as the square of air speed divided by the maximum normal acceleration.  Catch success mapped onto initial altitude and horizontal distance from the prey for non-maneuvering, smooth maneuvering and non-smooth maneuvering prey, and for 4 values of the navigation constant N: A low extreme (N = 1), the optimal value for catching nonsmooth maneuvering prey (N = 2.8), the optimal value for catching smooth maneuvering prey (N = 5.6), and a high extreme (N = 15). The yellow asterisks depict the global optima with respect to attack position and N, showing the attack strategy which uniquely maximizes catch success for a given prey motion. The yellow crosses denote local optima for a given N and prey motion. The approximate intercept speed was maximized by stooping from a low altitude (< 200m), leading to a low flight speed at the point of intercept (35 − 45ms −1 ). Optimal stoop altitude was somewhat higher (c. 350m) when prey maneuvered smoothly (Fig 4b), leading to a moderate intercept speed (50 − 55ms −1 ). Catch success with non-smoothly maneuvering prey (Fig 4c) was maximized by stooping from a very high altitude (c. 1500m), leading to a very high intercept speed (> 100ms −1 ) approaching the terminal velocity of the model-falcon (see Fig 3). Interestingly, catch success barely declined when the model-falcon attacked from a higher altitude than the optimum (Fig  4), but was greatly reduced if the model-falcon attacked from a lower starting position (Fig 4), so stooping from a high altitude is never a bad strategy provided that the guidance system is appropriately tuned (see below).

Optimization of the navigation constant for stooping
The attack strategy of a model-falcon encompasses both its initial position relative to the prey, and the setting of its navigation constant N. The global optima that we have so far discussed (asterisked points in Fig 4) assume joint optimization of the predator's initial attack position and its navigation constant N, and the optima for both parameters depend on the motion of the prey (see also Table 1). Selection on N is expected to be strongest when prey execute nonsmooth maneuvers, for which high catch success is achieved over only a narrow range of N (compare width of dark blue area denoting high catch success in Fig 5c with the equivalent areas in Fig 5a and 5b). Interestingly, for all three types of prey motion, the optimal setting of N tends to be lower the faster the stoop (see dashed lines in Fig 5 plotting the optimal setting of N conditional upon the speed at intercept). Conversely, for a given setting of N, the optimal intercept speed becomes lower the higher the value of N (see solid lines in Fig 5 plotting the optimal speed at intercept conditional upon the setting of N). Thus, for any given type of prey motion, high-speed, high-altitude stoops only maximize catch success over a small range of comparatively low values of N. At higher values of N, catch success is maximized by using a corresponding to the initial altitude is shown on the right of the graph. This is only an approximate relationship because the exact intercept speed depends on many factors within each hunt. Stooping by peregrine falcons: A physics-based simulation low-speed (Fig 5), low-altitude (Fig 4) attack, but this is generally less successful than a highspeed, high-altitude attack at a lower value of N. In summary, it turns out to be essential for our model-falcons to set their navigation constant appropriately: if a sub-optimal value of N were used, then stooping might no longer be the best attack strategy, because of poor catch success. For instance, if a model-falcon were to use the optimal value of N for smoothly maneuvering prey (N = 5.6) against prey executing non-smooth maneuvers, then a high-altitude stoop would be unlikely to result in prey capture (third panel of Fig 4c). This does not necessarily mean that a falcon must actively adjust N to match the maneuverability of its prey: the best attack strategy of a model-falcon against the best defensive flight pattern of a model-starling (i.e. non-smooth) involves entering a highspeed, high-altitude stoop at N % 3. This minimax strategy not only yields maximal catch success against non-smoothly maneuvering prey, but also yields near-maximal catch success against prey that are flying straight or maneuvering smoothly (second row of Fig 4). Hence, subject to the assumptions of our model, we expect falcons to adopt a general strategy of stooping from high-altitude at N % 3, because this strategy is effective against all of the different patterns of prey flight that we have tested here.
Some interesting flight trajectories emerge at N < 2 (see S2 Fig). In this case, the model-falcon exerts most of its acceleration towards the end of its attack (see also [14]), often diving below its prey before looping upward to intercept. This upward-curved trajectory is regularly observed in nature [4,23], and has previously been suggested to be a strategy of a falcon to fly into the blind spot of its prey's vision [9]. Our model provides a more parsimonious explanation for these flight paths, which can emerge naturally from the dynamics of the underlying feedback law.

Response delays and errors in vision and control drive the need to stoop
The most important factor that causes the reduction in catch success observed at high values of the navigation constant N is the response delay of the model-falcon. A robustness analysis (Fig 6a) shows that high values of N are no longer associated with a low catch success if the reactions of the falcon are effectively instantaneous (compare catch success as a function of N at τ = 0.1ms delay with the equivalent line for the default τ = 50ms delay used in our baseline model; see Materials and methods). Conversely, if the falcon's actual response delay is greater than the default assumed in our baseline model, then the optimal value of N is driven towards an even lower value (Fig 6a). Visual error also affects the optimal value of N in our simulations: if the falcon is subject to greater visual error than the default value assumed in our model, then the navigation constant is again driven towards an even lower value of N (Fig 6b). This reflects the fact that the propagation of this visual error into the commanded acceleration is directly proportional to N (see Fig 2). In contrast, the optimal value of N is robust to the error assumed in the control system itself (Fig 6b).
How do response delays, or errors in vision and control, impact the success of a given attack strategy? As expected, catch success declines as each of these quantities increases (Fig 2). Remarkably, however, the optimal starting altitude becomes lower when the visual error or control error is increased (Fig 6c and 6d). Thus, a high-speed, high-altitude stoop only maximizes catch success if the falcon is accurate in both vision and control. A high-speed stoop maximizes catch success for all of the response delays that we tested, noting that any much longer delay would have resulted in very low catch success (Fig 6e). On the other hand, if the falcon's response is effectively instantaneous (τ = 0.1ms), then 100% catch success is attained even in a low-altitude dive from < 200m. This implies that the falcon's flight performance is sufficient to catch a starling in a low-level stoop, but that delays in the model-falcon's response hamper its ability to catch prey. The lower catch success that results from having a slower response can be ameliorated by diving from higher altitudes at lower N.

What mechanisms underlie the increased catch success in a stoop?
When a falcon stoops from high altitude, its attack is characterized by both a very high flight speed, and a very steep descent angle-either of which could promote catch success. To investigate the effect of steepness of the descent, we altered the initial conditions of the simulations so as to model a horizontal attack at very high initial speed (112 ms −1 ). This effectively simulates the final approach of a falcon that stoops from a very high altitude to gain speed before levelling off to intercept. Remarkably, the maximum catch success of these model-falcons is only 3% lower than for those intercepting their prey at the same speed in a steep dive (61 vs 64%). This implies that the steep descent angle is not directly responsible for the overwhelming success of a stoop, and hence that the key reason for starting from a very high altitude is to gain airspeed by converting potential energy to kinetic energy.
The very high airspeed attained in a stoop enables model-falcons to exceed the modelstarlings' maximal load factor and roll acceleration (Fig 3). To test which of these two dimensions of flight performance causes an increase in catch success in a stoop, we artificially capped the maximal load factor or maximal roll acceleration of our model-falcons. We thus investigated the catch success of a bird flying at the same high speed achieved in a high-altitude stoop (> 100ms −1 ), but with the lower maximal acceleration associated with sustained level flight (30ms −1 ). Limiting either component of the model-falcon's flight performance resulted in a substantial drop in catch success (51% when limiting roll acceleration and 42% when limiting load factor). This suggests that the high speed that a falcon attains in a stoop is important partly because of the higher load factors and higher roll accelerations that can be achieved in high-speed flight. Interestingly though, model-falcons flying at high speed still performed considerably better than model-falcons in sustained level flight (31% vs 26%) even though the maximal load factor and roll acceleration was made the same, and even when these fast falcons levelled off before interception. This implies that a high flight speed is beneficial in and of itself, independent of the higher acceleration performance that is usually also associated with fast flight.
This result might seem surprising, because model-falcons fly faster than model-starlings even in sustained level flight (Fig 3a), and the most obvious consequence of increasing the falcon's flight speed further is to increase its turning radius, potentially causing it to overshoot when attacking sharply turning prey. However, the flight speed of a falcon varies continuously in our model, on account of its varying acceleration demand, and work on missile guidance and control has shown that the accelerations commanded in response to variations in speed are lower when the angle between the current line-of-sight and the expected point of intercept is smaller. The faster the falcon, the smaller this angle, which reduces the risk of control saturation, and thereby decreases the probability of missing the target.

Discussion
Previous authors have suggested that falcons stoop at high-speed to add an element of surprise to their attack, thereby preventing escape maneuvers by the prey. Our simulations suggest that there are many other, previously unrecognized advantages to stooping: for example, a highspeed stoop still provides a clear advantage over an attack from lower altitudes even if the prey individuals fly erratically, which models how real prey behave when alerted to the presence of a predator. This is because fast-flying falcons can obtain a higher load factor, can roll faster into a turn, and can slow down less when increasing their load factor to maneuver-each of which increases the falcon's catch success. Moreover, the steep angle of the stoop, and the details of the attack geometry of a fast-flying falcon further enhance catch success. The functional reasons for stooping are therefore far richer than considered previously, and are closely related to the physical constraints upon the problem (see [24] for a wider discussion of the operation of natural selection in relation to biomechanical constraint).
In order to intercept prey successfully at high speed, our model-falcons required a suitably optimized guidance law. Informed by a recent empirical study [15], our model-falcons used a pure proportional navigation guidance law, but we identified the optimal value of the navigation constant N post hoc through Monte Carlo simulation. A range of different values of N can be used successfully in low-speed flight, or when attacking non-maneuvering and smoothlymaneuvering prey, so the tightly-defined optimum of N % 3 that applies in a high-speed stoop works well under all of the conditions that we tested (Fig 4). This most broadly effective value of the navigation constant coincides closely with the median value of N = 2.6 that has been found empirically in captive-bred peregrine falcons attacking artificial targets [15]. It also coincides with a classical result of linear-quadratic optimal guidance theory, which shows that proportional navigation with an effective navigation constant N 0 = 3 minimizes the control effort needed to intercept a non-maneuvering target [14] (the effective navigation constant is defined as N 0 = N(v cos δ)/v c for pure proportional navigation, where v is the attacker's speed relative to the ground, v c its closing speed relative to the target, and δ the deviation angle between the attacker's velocity vector and its line-of-sight to target [14]). The fact that a navigation constant of N % 3 is found to be optimal or near-optimal in so wide a range of circumstances-and in so wide a range of systems, from birds to missiles-strongly indicates the robustness of our analysis and conclusions.
The results of our simulations also offer insight into the variability that is intrinsic to the attack behavior of real falcons. Captive-bred falcons have been shown to use a range of navigation constants around the median value of N = 2.6 that is comparable to the range of values of N that are optimal under the various conditions simulated in our model [15]. Whereas these real falcons seem to maintain an approximately constant value of N during a single interception attempt, they also vary it appreciably between attacks. It is currently unknown what drives this variability in the tuning of the navigation constant. Our model suggests that the details of the prey's motion, and the details of the falcon's flight speed, response delay, and precision in vision and control are all important determinants of the optimal tuning of N, and might therefore explain any adaptive variation in N.
Surprisingly, attacking at high speed does not require a faster response from the falcon than attacking at low speed (Fig 6e). By increasing the falcon's flight performance, the stoop compensates for the decrease in catch success that would otherwise result from a slow response. In contrast, the need for accuracy in vision and control is especially acute when flying at high speed, and stooping is only optimal if the falcon has reasonably low error in both (Fig 6). Hence, we assume that selection for high visual acuity will be especially strong in species that use high-speed attacks. It follows that stooping should be considered a specialist hunting technique, because only accurate falcons with optimized guidance will be able to increase catch success by stooping. This arguably poses an exploration-exploitation dilemma for a falcon learning to catch prey: either it may seek to optimize its current catch success by adopting the easy strategy of a low-speed attack, for which the details of the parameter tuning are not critical; or, it may explore the more difficult strategy of a highspeed stoop, which could decrease catch success at first in an unskilled falcon, but can be expected to increase catch success in the long-run. The playful attacks by falcons in which they do not seriously attempt to kill their prey [6], may be necessary for acquiring sufficient skill in stooping.

Limitations and wider implications
Although our physics-based model is realistic enough for its intended purpose, there are obviously further constraints in nature that we have not modelled here, including the effects of unsteady aerodynamics, the dynamics of pitch and yaw instability, and the mechanics of catching or knocking the prey with the talons at intercept. There are also other complicating factors that we have not modelled, including the effects of explicit evasive maneuvers by the prey, or the impact of intra-and inter-specific variation in flight morphology and physiology, and hence variation in the flight performance of predator or prey. These factors can be studied through extensions of the model and through parametric variation of the model between simulations, and will be considered elsewhere. Nevertheless, our approach to studying the dynamics of aerial predation is unique among behavioral studies of complex systems in combining guidance and control laws inspired by missile theory [14] with a detailed simulation model of the biology and physics of animal flight. The underlying feedback laws are well-founded in the theory of optimal guidance [14], and their validity as a phenomenological model of guidance and control in peregrine falcons has already been verified in nature [15]. Furthermore, the simulation approach that we have used proves necessary because of the complexity of the flight dynamics, which precludes an analytical approach [13]. Even setting aside the aerodynamic complexities that we have handled using a blade-element model of flapping flight, the mere fact that the birds must reorient their body to redirect their lift vector generates dynamics that are known to have no analytical solution in the most-closely analogous case of bank-to-turn missiles [25]. Our modeling therefore follows an embodiment approach, which states that behavior emerges through feedback-loops between the brain, the body, and the world. Aspects of cognition, such as the guidance laws used to intercept prey, are shaped by properties of the body and therefore bodily traits need to be considered to fully understand behavior [12,26]. In summary, our agent-based simulation approach provides insights into the optimization of attack strategies by an aerial predator that could not have been reached in any other way, and thereby paves the way for a new generation of studies into the optimization of complex multiagent flight behaviors.

Summary of simulation procedure
In our simulations, model-birds fly with six degrees of freedom through an open three-dimensional space without objects or boundaries. In each simulation run, a model-falcon aims to intercept a lone model-starling in mid-air, using a pure proportional navigation guidance law (Eq 1). In model-starlings, the guidance command is a forcing function that ensures that they either fly linearly, or execute smooth or non-smooth maneuvers, always keeping within ±20m of their initial altitude. Model-birds are subject to gravitational and aerodynamic forces, and flap, glide and retract their wings to manipulate the aerodynamic forces. Model-birds maximize their forward acceleration at a given speed and orientation, subject to the constraint that they meet the normal acceleration commanded by their guidance system. A flight controller determines the changes in wing shape and motion that best meet the desired acceleration. When the commanded normal acceleration cannot be met, model-birds simply exert the maximum attainable lift force.
At the start of an attack, the model-starling is located at the origin of the global coordinate system, with its body coordinate system oriented randomly. This variation in initial orientation ensures sufficient randomization to avoid artificial results due to coupling of highly specific initial conditions of the falcon and starling. The model starling begins flying at an initial speed of 11 ms −1 , calculated as the airspeed at which the cost of transport is minimized under the model. The model-falcon initially flies at a speed of 16 ms −1 , with its longitudinal body axis pointing directly towards the starling, and its lateral body axis horizontal. We parametrically vary the falcon's initial position relative to the prey, and vary the navigation constant N (i.e. the one free parameter of the falcon's guidance law; see below) to simulate a continuum of different attack strategies. For each attack, we sample at random from a uniform distribution, sampling the navigation constant N between 1 and 20, the falcon's initial altitude above the prey between −200 and 1500 m, and the initial horizontal distance to the prey between 0 and 800 m. The simulation ends when the falcon either intercepts the starling or is unsuccessful in its attempt to intercept, according to the criteria defined below. For a visualization of the simulations, see SI videos.

Analysis
Every simulation ends in either the success or failure of the model-falcon to catch the modelstarling. A catch is defined as occurring when the model-falcon comes within 0.2m of the model-starling. Failure occurs if either the falcon has not caught the starling within 40 s, or if it experiences a near-miss from which it cannot recover. A near-miss occurs when the model-falcon comes within 5.0 m of the model-starling, but subsequently finds itself further than this from the model-starling and with the model-starling in the blind zone of the model-falcon (a cone of 45˚behind the bird) such that the falcon would effectively need to begin a new engagement in order to re-acquire its target. In order to analyze how the model parameters affect catch success, we apply Generalized Additive Modeling (GAM; [19,20]). This is a nonlinear regression method which places no assumptions on the shape of the relationship between predictor and outcome. The estimation of the smoothing functions is conducted by automated cross-validation procedures (quadratic penalized likelihood), which reduce the likelihood of over-fitting and therefore ensure that our (conditional) maxima are not spurious. We applied GAMs with a logit link function, with catch success as the outcome variable and with the navigation constant N, initial-altitude and horizontal-distance as the independent variables. We built separate models for each combination of prey motion, response delay, and error. No constraints on the effective degrees of freedom were applied.

Software
Model simulations were programmed in C++, using openGL for graphics rendering. Hildenbrandt's StarDisplay model [27] was used as the framework for graphical display. Optimization studies of the blade-element model were conducted in MATLAB 2014a, and the mgcv package [28] of R statistics [29] was used for GAM regression.

Detailed description of the bird-flight simulator
Here we explain the detail of our simulation model, using the block structure depicted in Fig 2, and discussing each of the following four segments of the block diagram in turn: A. Kinematics, B. Vision, C. Guidance, D. Control and E. Aerodynamics. We justify each variable and mechanism by parameterization to empirical data, and justify mathematical argument in terms of physics or optimality. For symbol meanings, see Table 2.
A Kinematics. From a flight dynamics perspective, birds are implemented in the model as rigid bodies with six degrees of freedom. Model birds bank to turn, so as a simplification we take account only of their roll moment of inertia in modelling the rotational dynamics. The position vectorr determines the position of the bird in a right-handed inertial axis system (x, y, z) in whichẽ x ,ẽ y ,ẽ z are unit vectors defining the inertial axes, whereẽ z points upward and therefore defines the bird's altitude. The orientation of the bird is described by a rotating right-handed body-fixed axis system (x 0 , y 0 , z 0 ), in whichẽ x 0 ,ẽ y 0 , andẽ z 0 are unit vectors directed along the principal axes of the bird, and called the roll, pitch, and yaw axes, respectively. The roll axisẽ x 0 is assumed to be aligned instantaneously with the bird's forward velocityṽ, which amounts to assuming perfect pitch and yaw stability, and together with the yaw axisẽ z 0 is assumed to define the bird's plane of bilateral symmetry.
The kinematics of the model are governed by Newton's laws of motion. Numerical Verlet integration is used to solve the translational motion of the bird according to the following differential equations:r ðt þ dtÞ ¼rðtÞ þṽðtÞdt þ 1 2ã wherer is the position,ṽ is the velocity,ã is the acceleration, and m is the mass of the bird. The update time dt is the time step of the model at which both the numerical integration scheme and the flight forces are updated, and is set at 1 × 10 −4 s (see section L for convergence tests). The flight forceF is composed as follows: where TD is the net thrust minus drag force, L the lift force, and g the gravitational acceleration. The time-varying aerodynamic forces TD and L are outputted by the lift controller described in Section D.2. Model-birds respond to their normal acceleration command by reorienting their lift vector in a banked turn, applying lift asymmetrically so as to generate a net roll torque. To simplify the model, it is assumed that the exerted roll does not affect lift, thrust and drag.
Wing retraction is assumed to be symmetric, so that asymmetric lift forces are produced only by angle of attack asymmetries, and the roll moment of inertia depends appropriately on the retraction of the wings. The rotation of the body about its roll axis is updated as follows: where O is the roll angle, ω is the roll angular velocity, I x 0 is the total moment of inertia of the bird about its roll axisẽ x 0 , and M x 0 is the roll torque. The time-varying roll inertia I x 0 and roll torque M x 0 are outputted by the roll controller described in Section D.3. Eq 8 assumes that the body axes (x 0 , y 0 , z 0 ) are principal axes, and that no inertial coupling occurs. B Vision. The model-falcon is assumed to apply a pure proportional navigation (PPN) guidance law, such that the function of its visual system is to estimate the angular rate of change_ l in the line-of-sight to targetr d ¼r f Àr p , wherer f is the position of the falcon andr p the position of its prey. In other words, the function of the visual system is to estimate the angular rate of change in the line drawn from the model-falcon to the model-starling, which, importantly, is independent of gaze direction. Mathematically, the line-of-sight rate is wherer d is the line-of-sight vector measured with error, and whereṽ d is the corresponding velocity of the falcon relative to its prey. Although we assume that the falcon holds an explicit representation of the line-of-sight rate_ l, it is important to note that our use of Eq 9 does not imply that the bird must necessarily hold an explicit representation of its own position and velocity relative to its prey. For example, the line-of-sight rate on the lefthand side of Eq 9 could be estimated directly from the retinal coordinates of the target, provided that the falcon is able to subtract from this the apparent motion of the target due to rotation of its retinal coordinate system. Rather, Eq 9 should be thought of as providing a convenient way of computing the line-of-sight rate in the simulation, which simultaneously allows us to model the effects of error in the falcon's visual system phenomenologically.
We incorporate this visual error by defining a random error vector whose magnitude is drawn from the uniform distribution from 0 to ξ (see Section F for explanation and justification of the magnitude of ξ), and whose direction is drawn from the uniform circular distribution around the true line-of-sight vectorr d . This visual error vector explicitly models the effects of angular error in the estimation of the direction of the line-of-sight vectorr d , and is therefore scaled by the target range jr d j when computing the estimate of the line-of-sight vector required by Eq 9:r We calculate the corresponding relative velocity as: where τ is the differencing time, which we also take to represent the sampling rate from vision to guidance in the falcon. C Guidance. Changes in the model-falcon's velocity are commanded by a pure proportional navigation (PPN) guidance law. The input to this guidance law is the estimated line-ofsight rate_ l from Eq 9, whilst the output from the guidance system to the controller is the commanded accelerationã Ã . In three dimensions, the PPN guidance law is defined as: whereṽ is the falcon's velocity and N is a constant of proportionality called the navigation constant. The form of this equation is such that the accelerationã Ã commanded under PPN guidance is always perpendicular to the falcon's velocity vectorṽ.
The guidance system of the model-starlings is described by one of three different kinds of forcing-function: straight flight, smooth maneuvers, and non-smooth maneuvers (see also Fig 1).
C.1 Straight flight. In straight flight, acceleration is commanded on a constant heading drawn at random from a uniform circular distribution, and at a constant elevation drawn at random from a uniform distribution between ±2.5˚. In combination with its random initial orientation, this forcing function causes the model-starling to turn towards an almost horizontal flight direction within 1 or 2 s of the start of the simulation, after which it flies in a straight line whilst maximizing its forward acceleration. The magnitude of the commanded accelerationã Ã is such that its combination with gravitational acceleration does not exceed the starling's maximum normal acceleration as calculated in Section D.1.

C.2 Smooth maneuvers.
In smooth turning, centripetal acceleration is commanded to the bird's left or right, with oscillations in the magnitude of the commanded acceleration specified according to a harmonic forcing function. This smooth forcing function is described by the following equation:ã where t is the time in seconds, and κ is the difference between the current and initial altitude of the starling. The last term in this equation ensures that the starling always remains close to its starting altitude. The unit vectorh is perpendicular to the model-starling's velocity vector and is directed horizontally:h The terms c 1 , . . ., c 4 in Eq 13 were optimized by varying them parametrically (see Table 3 for optimized settings) so as to maximize the mean magnitude of the normal acceleration, subject to the constraints that: 1. the mean roll rate remains below 30 rad s −2 ; 2. the starling exerts accelerations equivalently in all directions (see S1 Fig); 3. the starling flies within bounds of ±20 m of its initial altitude. The starling does not exert maximal normal acceleration at each time step, because it slows down when maneuvering, thereby decreasing the maximum normal acceleration in a future time step. The optimized smooth forcing function results in high load factors (mean load factor: 3.4) and low roll accelerations (mean roll acceleration magnitude: 27 rad s −2 ).
C.3 Non-smooth maneuvers. In non-smooth turning, acceleration of approximately constant magnitude is commanded in a stepwise fashion in a randomly varying direction close to the horizontal. The non-smooth function has the following equation: whereq is a random unit vector that is updated by the following equation: where U denotes sampling from a uniform distribution, and whereQ is a unit vector pointing towards a random azimuth angle and elevation between −c 9 and c 9 degrees. The parameters c 5 , ‥, c 9 are optimized as for the smooth forcing function (see Table 3), with the modification that the non-smooth forcing function also aims to maximize roll acceleration. The nonsmooth forcing function results in similarly high load factors to the smooth forcing function (mean load factor: 3.4), but involves much higher roll accelerations (mean roll acceleration magnitude: 2012 rad s −2 ). D Control. The bird's flight controller ensures that the acceleration commanded by the guidance law is achieved or approximated by making the appropriate adjustments to the wing motion and shape. The bird's flight controller is subdivided into three subsystems determining weight support, lift control, and roll control, which we now discuss in turn.
D.1 Weight support. To achieve a change in velocity with the magnitude and direction of the commanded acceleration, gravitational acceleration needs to be considered at each time point. Specifically, the sum of the centripetal acceleration due to lift and the component of gravitational acceleration that is perpendicular to the bird's velocity should equate the commanded centripetal acceleration. Hence, to determine the required magnitude and direction of the lift force, gravitational acceleration is first subtracted from the commanded acceleration:ã steer ¼ã Ã À gẽ z ð17Þ then projected onto the plane defined by the bird's transverse and dorsoventral axes,ẽ y 0 andẽ z 0 : The magnitude of the resulting accelerationã projected is then signalled to the lift controller (described in Section D.2), whilst its direction is signalled to the roll controller (described in Section D.3).
D.2 Lift control. Lift control is governed by the acceleration objective described in section E: maximize forward acceleration given current airspeed and body orientation, subject to the constraint that the lift needed to meet the guidance commands is exerted. This objective could be achieved in multiple ways, whether by varying the wingbeat kinematics when flapping, or by retracting the wings to reduce drag when gliding. In this section, the wing-beat averaged equations that determine which of these flight modes is selected (see section E for a derivation of these equations).
The magnitude of the desired lift L Ã is mjã projected j. Constraints on the achievable lift arise physiologically due to due to constraints on the torque forces that the flight muscles can sustain, and aerodynamically due to wing stall, where c l.max is the lift coefficient at the stall limit. To calculate the maximum achievable lift for which the muscle torque constraints are not exceeded, we use the allometric scaling rule L 0 = 1.7mg, where L 0 is the maximum lift at maximum wing span [11] (see section G for a justification of the mechanical constraints in our model). In flapping flight, we assume no wing retraction, and thus L max = L 0 is the upper limit.
In gliding flight, model-birds retract their wings, and the maximum achievable lift L max is found by [11]: where b ml is the wing span that maximizes lift, b max and b min are the minimum and maximum wing span respectively, S w.max is the wing area at maximum span, and S w.ml is the wing area that maximizes lift, ρ is the air density, and v is the airspeed. The achievable lift L is the lower value of the desired lift L Ã and the maximum achievable lift L max . The following equations relate the achievable lift L to the maximum net thrust (or minimum drag): where c l is the wingbeat-averaged lift coefficient, c l.max is the maximum achievable wingbeataveraged lift coefficient, c d.fric is the friction drag coefficient on the wing, c d.body is the total drag coefficient of the body, S w is the maximum projected wing area, and S b is the frontally projected body area, A R is the wing aspect-ratio, b is the actual wing span, θ is the stroke angle from the highest point of the wing-beat to the lowest point, f is the wingbeat frequency, and l w is the wing length. The coefficient c torque accounts for the constraints in torque forces that the flight muscles can sustain. The torque forces increase linearly with airspeed v, for a given wing-beat frequency and optimal wing-twist, and exceed the sustainable threshold at a point v thresh which is determined through simulation of the blade-element model described in section E (note that this expression is almost equivalent to stating that the maximum available power for flight is constant across values of airspeed). The coefficient c torque is calculated as: It is important to note that the wingbeat-averaged lift coefficient c l can differ substantially from the local lift coefficient at a particular section of the wing at a particular time point in the wingbeat. For instance, when the downstroke delivers an upward lift force, and the upstroke a downward lift force of equal magnitude, c l will be zero. When the bird is flapping, f is set to the maximum observed wingbeat frequency for the species, jointly with a stroke amplitude derived using allometric scaling rules (see Table 4). When the bird glides, f = 0 by definition, and the wing span b is optimized to achieve L with the lowest amount of drag. The wing span b determines the wing area according to the relationship: and determines the aspect ratio as: To find the optimal wing span b Ã that achieves a given lift L with the least amount of drag, the following steps are applied. First b Ã is calculated by:  Stooping by peregrine falcons: A physics-based simulation truncated, b Ã is recalculated as: and b Ã is truncated again to fall within the closed interval [b min , b max ]. The bird chooses either to flap or to glide according to which mode of flight achieves or most closely approximates L Ã with the largest net thrust (or least drag). D.3 Roll control. The roll controller determines the roll acceleration used to reach the desired bank angle. The desired direction of roll by the bird is determined by calculating the angle between a projected andẽ z 0 : To roll in this direction, the bird must apply the appropriate roll acceleration. As Eq 8 implies, roll acceleration is determined by the the aerodynamic roll torque M x 0 and the roll moment of inertia I x 0 . The inertia is dependent on the wing span b, whilst the net torque depends on the differential lift exertion between left and right wings, together with the wing span. Therefore, the goal of the roll controller is to find the optimal wing span and lift exertion to reach the desired bank angle in the minimum amount of time. We assume that on average the lift force acts at half the wing length, thus the total torque M generated by one wing is: The roll moment of inertia is calculated as: To obey the guidance law, the bird must roll until it reaches the desired bank angle, as dictated by the guidance algorithm, and must then stop rolling. The faster the bird reaches the desired bank angle, the more accurately it obeys its guidance law. The fastest way to reach a given bank angle is to accelerate maximally until some point q, and then decelerate maximally (called Bang-Bang control). The highest roll acceleration is achieved at the wing span that maximizes lift production (see section I for a proof). Given this maximal acceleration, the bird must know this point q in the roll at which it should start to decelerate in order to reach a roll velocity of zero at the desired bank angle. We calculate q as follows. Let ω be the roll (angular) velocity and γ (see Eq 30) be the angular position relative to the desired bank angle. The angular acceleration _ o is either the positive or negative maximal angular acceleration. The angular velocity is where ω 0 is the angular velocity at the current time step. We calculate the time t it takes to reach γ by solving where the minimum positive t is the answer. If q (the point at which the bird should start to decelerate) is exceeded, there are 2 positive solutions for Eq 35, according to whether a has the opposite or same sign as γ (which is the case when turning towards the desired bank angle and decelerating). If the bird has not yet reached q, then there are no real solutions. These equations have no real solutions when: Thus the following rule obeys Bang-Bang control: _ o should have the same sign as γ, unless both ω has the same sign as γ and o 2 0 > 2 _ og, in which case _ o should have the opposite sign. When o 2 0 > 2 _ og, the equations have real and positive solutions when decelerating, which means that γ is reached within a positive time. The first moment in time when this is the case must be closest to the point q. Finally, the output of the controllers is fed to the kinematics, which completes the feedback loop. E Aerodynamics. We used a multistage modeling approach to model the aerodynamic forces on the birds. First, we constructed a blade-element model of flapping and gliding flight, assuming sinusoidal wingtip kinematics and quasi-steady flow [30][31][32][33][34][35][36]. Since our model-birds mostly flew at low Strouhal numbers (St < 0.22), this quasi-steady analysis provides sufficient accuracy for our modeling objectives [31,37]. We assumed that the relation between the lift coefficient and induced drag coefficient was that of an elliptical wing under classical lifting-line theory [38], and incorporated realistic physical constraints on the maximum aerodynamic torque, to account for the limits of sustained force production by the pectoralis and supracoracoideus [39,40]. All of the morphological parameters needed to model body, friction, and induced drag, as well as lift and thrust were derived from empirical estimates [11,30,[41][42][43][44][45], or estimated using allometric scaling laws. We then used a combination of regression statistics and analytic methods to derive a set of equations relating the maximum forward acceleration over a wingbeat to the bird's load factor (i.e. lift divided by body weight) and airspeed (Fig 7a). To determine the maximum forward acceleration for a given load factor and airspeed, we optimized the time-varying wing twist and mean angle of attack, as well as the stroke plane, and whether the bird should flap at the maximum wingbeat frequency or glide and retract its wings. To account for errors in control, we modified the load factor by a multiplicative error term (χ) proportional to the rate of change in the lift coefficient, setting χ = 0 in our baseline model. Here, we present a derivation of the time-averaged equation (Eq 23) that we used in section D.
The goal is to come up with a simple approximate function Where TD is the time-averaged maximum thrust minus drag, L is the lift, v the airspeed, andm is a vector of morphological parameters that can be gathered empirically. E.1 Assumptions. We assume the following aerodynamic and kinematic properties. For each blade-element the equations for lift and drag are: Where S(i) denotes the total area of the element i, and U(i, p) is the speed of air relative to the (moving) blade at point p of the wing-beat cycle.Lði; pÞ is defined to be perpendicular tõ U ði; pÞ in the plane spanned by the blade, whileDði; pÞ is parallel toŨ ði; pÞ and has the opposite sign. The drag coefficient c d is composed of the body, induced and wing-friction drag coefficients: where the body drag coefficient c d.body is referenced to the body's frontal projected area S b . The induced drag coefficient c d.induced is related to c l by [46]: which is the classical limit from lifting line theory for a wing with an elliptical lift distribution, and we apply Blasius' solution for laminar boundary layer of a flat plate [47] to obtain the friction drag coefficient: where Re denotes the Reynolds number, which is calculated as with μ the dynamic viscosity of the air. Furthermore, the lift coefficient c l is related to the angle Fig 7. (a) Look-up table for the accelerations due to the aerodynamic forces acting on the falcon. At each model time-step, the falcon maximizes forward acceleration (minimizes deceleration), given its forward speed and with the constraint that load factor is set to achieve the net commanded acceleration by the guidance law. If this constraint cannot be met (i.e. if it is unfeasible due to aerodynamics or high resulting torque forces), the closest approximation of the load factor is chosen. In the blade-element model, the falcon optimizes the wing twist, the wing's angle-of-attack, the wingbeat frequency and the wing retraction. Inside the trapezoidal contour, the falcon flaps at maximal wingbeat frequency, and outside it the falcon glides. Above the contour, flapping results in too high torque forces on the wing. Gravity is excluded from the accelerations in the figure. (b) The partial derivative of forward aerodynamic acceleration with respect to load factor, termed the "aerodynamic deceleration penalty". https://doi.org/10.1371/journal.pcbi.1006044.g007 Stooping by peregrine falcons: A physics-based simulation of attack α by the higher-order lifting line theory [24]: We further assume that the wing shape is elliptical, that there is no side-slip, that the airfoil is symmetric and that the wing-motion is sinusoidal.
E.2 Optimization problem. Our aim is to find the function c Ã l ði; pÞ that maximizes TD under several constraints. Because the friction drag and body drag are independent of c l , we maximize TD i , which is thrust minus induced drag. The objective function is a double integral over the wing length i and wing-beat cycle p: where c l.max is the lift coefficient at the stall limit, t max is the maximum torque that the flight muscles can sustain, and where where U w (i, p) is the speed of blade-element i in section p of the wing-beat cycle, U a is the airspeed and U 2 ¼ U w ði; pÞ 2 þ U 2 a . To start solving this objective function, we make a few observations. First, the unconstrained optimal c Ã l ði; pÞ can be written as and when U a ) U w (i, p), the value of the unconstrained objective function remains constant with respect to changes in airspeed and c Ã l ði; pÞ converges to zero. Also, because the wings are symmetric, the maximum of the unconstrained TD i always corresponds to L = 0. Therefore, to get the maximum value of the unconstrained TD i , we solve for lim U a 1 TD i . First, we rewrite the equation: where c 5 ¼ 1 2 Sr. Removing all brackets: We can note that: and lim U a !1 Therefore Thus which is where n(i) is the fraction of the total wing area. Next, numerical simulations show that, despite the non-linearities caused by the inequality constraints, we can accurately fit a quadratic model of the form: Also, our numerical simulations show that, when L is maximized, TD i by flapping converges to TD i by gliding. Therefore, we know two properties of the quadratic function: Because TD i ¼ TD i:max when L = 0, we know that b ¼ TD i:max . From Eq 62 we know that: Therefore a is: The relation between thrust and lift becomes: which can also be written as: All that is left to do is to compute the double integral over the wing span and wing-beat cycle: Assuming sinusoidal flapping, the position (height) of the wing section i in time is determined by where t denotes the time. Differentiating this with respect to t, we get the speed of the wing section over time Assuming an elliptical wing, the double integral is: We first compute the inner integral: and evaluated at x = 1, we get: We thus need to solve The solution of the integral is And evaluated at 1 and 0, we get: Thus the solution of the double integral is Of which by far the most significant factor is the first term. Keeping only the first term, we get Thus the final thrust becomes: Subsequently, the coefficient for the decrement in thrust due to torque is added, as well as the friction and body drag, to produce Eq 23.
F Justification of the bounds on visual error. In the absence of any detailed information on visual processing in falcons, we estimate the magnitude of visual error by the falcon in the following way. Hodos et al [48] found that the minimum retinal image velocity for the detection of motion in pigeons is about 8˚s −1 . We assume that this detection threshold corresponds to the point at which the apparent motion due to target motion exceeds the apparent motion due to random noise. We further assume that the differentiation time of the falcon is 50ms and that the apparent motion due to target motion must exceed the apparent motion due to random noise in this short time frame. Therefore, the retinal image must move at least 0.007 rads per 50ms to be perceived. We thus take this value as the bound ξ on the visual error (see Section B). We acknowledge that this method is somewhat arbitrary, but consider that stating it in this way helps to provided confidence that the assumed numerical bound on the visual error is likely to be of about the right order of magnitude.
G Justification of the mechanical constraints. Whether a bird is flapping or gliding, aerodynamic loads place strain on the muscles, ligaments, and bones, in particular when the bird is performing maneuvers. In our simulations, the constraints on the production of aerodynamic forces are such that the muscles, ligaments and bones would not break or tear in real birds, and that the power required to produce these forces would not exceed the muscle power available to counter the resulting torque. Here, we justify that we only need to consider the constraint on the torque that arises around the shoulder due to lift production, as this constraint is likely to be exceeded before any other constraint.
The torque around the shoulder needs to be countered by the flight musculature, including in particular the pectoralis and supracoracoideus muscles. Bayer et al. [39], provide detailed measurements of the forces and constraints on these muscles in gliding flight of the domestic pigeon Columba livia. We use their measurements to verify Pennycuick's [40] and Tucker's [11] proposal for a general allometric scaling law for maximum lift production in birds. In horizontal gliding flight, lift balances weight (mg). Therefore, each wing carries half the body weight, which we assume to be applied at half the total wing length (l w ) in producing a torque around the shoulder joint. Because the lever arm of the pectoralis is much shorter than this, it needs to act with a force of 7mg, where m is mass and g is gravitational acceleration. The maximum force that can be sustained by the pectoralis is roughly 3 to 4 times the maximum load that the muscle is required to take in gliding flight [39]. With the wings fully outstretched, this would imply a maximum acceleration of 3g to 4g, corresponding to an aerodynamic load of 1.5mg to 2mg on each wing, and hence to an aerodynamic moment of 0.75mgl w to 1mgl w . Pennycuick and Tucker propose the scaling law of a maximum lift of 1.7mg for fully stretched wings, which corresponds to a torque force of 0.85mgl w per wing, exactly in the middle of the measurements by Bayer et al. [39]. When a bird retracts its wings, the moment arm decreases. Hence the wing's maximum load increases. Tucker's model reveals that peregrine falcons, with appropriately retracted wings, may experience forces of over 18mg during a pull-out of the stoop, resulting in 9mg per wing.
Aerodynamic loads also place pressure on the ligaments and bones. Pennycuick [49] extensively examined the strength of pigeons' (Columba livia) wing bones. He investigated the breaking points in bending and twisting of the humeri and the radio-ulna. Under simple assumptions of the lift distribution on the wing during gliding flight and with fully stretched wings, the lift would need to exceed 16.8 N before torsional or bending pressures would break the bones. As the weight of his pigeons was on average 0.393 kg, this would imply that the wings could hold up to 8 or 9 times the body weight, roughly equivalent to the maximum load that the falcon's pectoralis can hold when it stoops at high speed. Therefore, the constraint on torque implicitly contains the constraint on aerodynamic loads due to breaking of bones. Further evidence that the peregrine falcon is well adapted to cope with the large aerodynamic forces that arise during maneuvers in a high-speed stoop is provided by Schmitz et al. [50]. These authors investigated the mechanical properties of the feathers of four species, including the peregrine falcon, and found that the feathers of the falcon had larger cross-sections and protrusions than the other species, resulting in a greater bending stiffness, thereby allowing greater aerodynamic loads.
H Justification of the inertia equation in the model. Here, we derive Eq 32 that determines the whole body inertia in the model for a given wing retraction. The total roll moment of inertia I x 0 is the sum of the inertia of the body I b and the inertia of the wings I wing .
H.1 Body moment of inertia. The body inertia about the roll axis is only available for one species (the rose breasted cockatoo; Table 2 of [41]). We apply several allometric scaling laws to determine the inertia of other species based on these data. According to Nudds and Rayner [42], we can adopt the body frontal area and mass relation: The scaling of the body width has the relation Because these approximately preserve relationships in scaling, the inertia of the body in the roll direction scales as follows: Where I b0 is the inertia for a width of 1m and a mass of 1kg. The rose breasted cockatoo weighs 0.293 kg with I b = 1685.5 Á 10 −7 kg m 2 . Therefore, we estimate I b0 to be 0.1346717. With I b0 , we can now estimate the body inertia for other birds. For instance, using these equations, the I b of a male Peregrine falcon is estimated to be 448.0 Á 10 −6 kg m 2 .
H.2 Wing moment of inertia. The wing inertia has been calculated for many species by Rayner and Van Den Berg [43]. Their calculations determine the moment of inertia about the attachment site of the wing to the body (i.e. the shoulder joint), whereas we are interested in calculating the moment of inertia of the wings about the center of mass of the bird. There can, of course, be a considerable difference between these two quantities (see Table 2 of [41]), but they are easily inter-converted as follows. The moment of inertia about the center of mass is given by: where m i is the mass of the ith blade element of the wing, r i its distance from the shoulder joint, and x the distance from the shoulder joint to the center of gravity. We may restate this equation in terms of the total mass of the wing m wing and its total moment of inertia about the shoulder joint I wing as follows: where the summation may be calculated using data from Rayner & Van den Berg [43] and Hedrick & Biewener [30,41] regarding the mass distribution of the wing. We assume that the distance x from the shoulder to the center of mass is half the body width. Empirically, if we do not know the wing length, but do know the body mass, then we can approximate the wing inertia as follows: where J is the summation term of Eq 32. And, with a decreased wing span, we calculate: where l wing is the wing length and where ϕ is the fraction of wing retraction defined as: Thus the total roll moment of inertia at a given point in time (assuming wing retraction is symmetrical such that the center of gravity does not change) is: I Proof that the wing span that maximizes lift also maximizes roll acceleration. Here, we explain why Eq 19 maximizes the roll acceleration with respect to wing span. Let us call this maximizing wing span b Ã . Roll acceleration is maximized when one wing produces the largest positive lift, and the other wing the largest negative lift. The roll acceleration depends on the then set it to zero À 4ðL Ã Þ 2 À 2pb 3 rv 2 À c d:fric S max 2ðb max À b min Þ rv 2 ¼ 0 ð98Þ When b is less than b min , b is set to b min . This new value of b may cause c l to exceed the constraint. Therefore, we set c l to c l.max , and solve for the corresponding b. If b is greater than b max , we set it to b max and solve for c l .
K Comparison between flight performance in the model and empirical measurements. In order for the model results to be relevant for our understanding of real falcons, the flight performance of model-birds should approximate empirical measurements in falcons and starlings, as flight performance is expected to be a key determinant of catch success in a chase. However, very little is known about the peak flight performance of birds. Most of the flight performance mentioned elsewhere is derived using modeling, and we have mentioned how our modeling relates to previous models throughout the section A-J. Here, we present the few coarse comparisons that are currently possible. First, minimum sustained flight speeds (the minimum speed of the model-starling is 4.5ms −1 and that of the falcon 7.3ms −1 ; defined as the minimum airspeed where acceleration > 0 in Fig 3a) correspond closely with the minimum speeds at which birds fly in wind tunnel experiments (starlings are not reported to fly below 6ms −2 [52, 55, 56]). Warrick [57] reports that the linear acceleration of a starling is approximately 5.8ms −2 when released to fly through a 2x2x5 tunnel, and our model-starlings linearly accelerate 4.8ms −2 . The top horizontal speed of peregrine falcons has been found to be 27.6ms −1 [58], which is very close to the speed of 28.4ms −1 in the model. Starlings have been reported to occasionally fly at 22ms −1 during migration [59], so the top speed of 24ms −1 in the model also seems appropriate. The maximum dive speed of the falcon is 122ms −1 and resembles the speed of 108ms −1 achieved by a peregrine falcon in an unpublished study by National Geographic (to date, no peer-reviewed articles include measures of the peak performance of falcons in a stoop). The maximum load factors attained in the model are considerably higher than those measured experimentally. Ponitz et al. [3] report that the load factor of their peregrine falcon during pull out was 1.15 (including gravity) at a speed of 20-22ms −1 . Model-falcons are able to generate a load factor of 2.5 at that speed (excluding gravity), but it is reasonable to assume that a real falcon will generally pull out with a load factor below maximum, to maintain stability throughout its flight. Lastly, the high roll accelerations of our model birds are of the same order of magnitude as the roll accelerations measured in pigeons (columba livia) [60]: pigeons had an average whole body angular acceleration of 601 rad s −2 flying 3-6 ms −1 . At that speed, a model-starling exerts *1200 rad s −2 and a model-falcon *400 rad s −2 .
L Tests of convergence of numerical integration. In our simulations, we discretized time and constructed a set of ordinary difference equations to numerically solve the initial value problem of the original differential equations (see Eqs 2, 3, 4, 6 and 7). Here we test the convergence of the position, speed and acceleration of predator and prey, as well as the falcon's catch success. As we do not have access to the algebraic solutions, we tested the convergence as follows. The time step of integration in the simulations was reduced iteratively in orders of magnitude until an acceptable discrepancy with next step size of integration was obtained, where we define the discrepancy as: where e dt is the discrepancy, k the number of model time steps, f(dt) the difference equation describing the value of interest, where dt is the applied time step. We simulated trajectories of stooping falcons (800m altitude and 600m horizontal distance) against erratically or circularly moving prey, where the initial seeding of random variables was the same for each applied time step dt and k was chosen such that the falcon intercepted the starling in simulations using one of the time step sizes, and thus for no dt sizes the simulation continued beyond interception. An acceptable discrepancy was achieved for dt = 1 Á 10 −4 s (see Table 5 and S4 Fig). Between dt = 1 Á 10 −4 s and dt = 1 Á 10 −5 s the catch success varied only 0.13% which is much lower than the minimum difference in main results between tested conditions (2%). Beyond dt = 1 Á 10 −6 s, e dt increases due to floating point imprecision. Additionally, we have tested the convergence and stability of the numerical integration of Eq 6, in a scenario where convergence is expected to be weak on theoretical grounds. Namely, when the birds exert large roll accelerations, oscillations around the equilibrium bank angle may arise at larger dt. S5 Fig shows the bank angle change over time of a falcon stooping at 100ms −1 , and exerting a roll acceleration of 5000 rad s −2 using Bang-Bang control to turn pi rad. The maximum absolute error of a dt = 1 Á 10 −4 s with respect to dt = 1 Á 10 −5 s is only 0.007 rad, which is unlikely to affect the model outcome. We note that, although the model equations are nonlinear, they are inherently stabilizing (i.e. robust to perturbations due to the numerical scheme), as aerodynamic drag ensures that the speed converges towards an equilibrium point, and the guidance law ensures that the falcon is attracted to the prey (see the smooth variation of catch success with respect to variation in initial conditions in Figs 4 and 6). In some curved flight paths, the peregrine falcon drops below the prey; a phenomenon also observed in nature. (EPS)