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.

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.

Every simulation ends in either the success or failure of the model-falcon to catch the model-starling. 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.

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 vector determines the position of the bird in a right-handed inertial axis system (x, y, z) in which , , are unit vectors defining the inertial axes, where 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′, y′, z′), in which , , and are unit vectors directed along the principal axes of the bird, and called the roll, pitch, and yaw axes, respectively. The roll axis 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 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: (2) (3) (4) where 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 force is composed as follows: (5) where is the net thrust minus drag force, L the lift force, and g the gravitational acceleration. The time-varying aerodynamic forces 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: (6) (7) (8) where Ω is the roll angle, ω is the roll angular velocity, I x′ is the total moment of inertia of the bird about its roll axis , and M x′ is the roll torque. The time-varying roll inertia I x′ and roll torque M x′ are outputted by the roll controller described in Section D.3. Eq 8 assumes that the body axes (x′, y′, z′) 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 in the line-of-sight to target , where is the position of the falcon and 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 calculated as: (9) where is the line-of-sight vector measured with error, and where 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 , 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 vector . This visual error vector explicitly models the effects of angular error in the estimation of the direction of the line-of-sight vector , and is therefore scaled by the target range when computing the estimate of the line-of-sight vector required by Eq 9: (10) We calculate the corresponding relative velocity as: (11) 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-of-sight rate 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: (12) 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: (13) 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 vector is perpendicular to the model-starling’s velocity vector and is directed horizontally: (14) 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). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Table 3. Model settings. https://doi.org/10.1371/journal.pcbi.1006044.t003 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: (15) where is a random unit vector that is updated by the following equation: (16) where U denotes sampling from a uniform distribution, and where 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 non-smooth 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: (17) then projected onto the plane defined by the bird’s transverse and dorsoventral axes, and : (18) The magnitude of the resulting acceleration 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 . 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]: (19) (20) (21) 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): (22) (23) where c l is the wingbeat-averaged lift coefficient, c l.max is the maximum achievable wingbeat-averaged 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, 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: (24) 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: (25) and determines the aspect ratio as: (26) PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Table 4. Morphological parameters of the bird species in the model. https://doi.org/10.1371/journal.pcbi.1006044.t004 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: (27) (see Section J for derivation), and then b* is truncated if required to fall within the closed interval [b min , b max ]. An additional upper bound on b* is the wingspan b at which the torque constraint is exceeded when exerting L: (28) c l is subsequently calculated by Eq 22. Again, c l is truncated not to exceed c l.max and if truncated, b* is recalculated as: (29) 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 : (30) 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′ and the roll moment of inertia I x′ . 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: (31) The roll moment of inertia is calculated as: (32) (see section H for justification), where I x′ is the total moment of inertia about the roll axis, I b is the body inertia, I wing the inertia of one wing around the shoulder, ϕ the proportion of wing retraction, m w the mass of the wing, m b the body mass, and J = ∑ i m i r i , where r denotes the distance from the bird’s shoulder when the wing is fully extended and where m i denotes the mass of a blade element of the wing at that point. 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 is either the positive or negative maximal angular acceleration. The angular velocity is (33) (34) where ω 0 is the angular velocity at the current time step. We calculate the time t it takes to reach γ by solving (35) 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: (36) Thus the following rule obeys Bang-Bang control: should have the same sign as γ, unless both ω has the same sign as γ and , in which case should have the opposite sign. When , 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–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–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. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: 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 The goal is to come up with a simple approximate function (37) Where is the time-averaged maximum thrust minus drag, L is the lift, v the airspeed, and 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: (38) (39) 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. is defined to be perpendicular to in the plane spanned by the blade, while is parallel to and has the opposite sign. The drag coefficient c d is composed of the body, induced and wing-friction drag coefficients: (40) 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]: (41) 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: (42) where Re denotes the Reynolds number, which is calculated as (43) with μ the dynamic viscosity of the air. Furthermore, the lift coefficient c l is related to the angle of attack α by the higher-order lifting line theory [24]: (44) 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 that maximizes under several constraints. Because the friction drag and body drag are independent of c l , we maximize , which is thrust minus induced drag. The objective function is a double integral over the wing length i and wing-beat cycle p: (45) with the equality constraint that lift equals the desired lift: (46) and the inequality constraints (47) (48) 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 (49) (50) (51) (52) 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 . To start solving this objective function, we make a few observations. First, the unconstrained optimal can be written as (53) and when U a ≫ U w (i, p), the value of the unconstrained objective function remains constant with respect to changes in airspeed and converges to zero. Also, because the wings are symmetric, the maximum of the unconstrained always corresponds to L = 0. Therefore, to get the maximum value of the unconstrained , we solve for . First, we rewrite the equation: (54) where . Removing all brackets: (55) We can note that: (56) and (57) Therefore (58) Thus (59) which is (60) 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: (61) Also, our numerical simulations show that, when L is maximized, by flapping converges to by gliding. Therefore, we know two properties of the quadratic function: (62) (63) Because when L = 0, we know that . From Eq 62 we know that: (64) Therefore a is: (65) The relation between thrust and lift becomes: (66) which can also be written as: (67) All that is left to do is to compute the double integral over the wing span and wing-beat cycle: (68) Assuming sinusoidal flapping, the position (height) of the wing section i in time is determined by (69) where t denotes the time. Differentiating this with respect to t, we get the speed of the wing section over time (70) Assuming an elliptical wing, the double integral is: (71) We first compute the inner integral: (72) and evaluated at x = 1, we get: (73) We thus need to solve (74) The solution of the integral is (75) And evaluated at 1 and 0, we get: (76) Thus the solution of the double integral is (77) Of which by far the most significant factor is the first term. Keeping only the first term, we get (78) Thus the final thrust becomes: (79) 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′ 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: (80) The scaling of the body width has the relation (81) Because these approximately preserve relationships in scaling, the inertia of the body in the roll direction scales as follows: (82) 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 m2. 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 m2. 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: (83) 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: (84) 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: (85) where J is the summation term of Eq 32. And, with a decreased wing span, we calculate: (86) where l wing is the wing length and where ϕ is the fraction of wing retraction defined as: (87) 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: (88)

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 wing span in the following way: (89) with the constraint that M ≤ M max . When b > b* the torque will not change, but the inertia will increase, so the roll acceleration is less. When b < b*, the torque declines by ϕ2, while the inertia declines less than ϕ2. Therefore, b* maximizes the roll acceleration.

J Derivation of optimal wing retraction. The object of optimizing wing retraction is to maximize the net thrust minus drag, given a desired lift L*. We assume a simple relation between wing retraction, aspect ratio and wing area, as proposed by Tucker [11]: (90) (91) Assuming this relation, the total induced drag becomes: (92) and the total thrust minus drag is (93) where β is the dive angle (i.e. the negative of the elevation angle). The optimization becomes: (94) with the equality constraint that lift equals some desired value: (95) and the inequality constraints: (96) The optimization problem turns out to have a simple solution. Namely, we calculate the partial derivative with respect to b (97) then set it to zero (98) and solve for b: (99) 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.