## Abstract

Undulatory swimming animals exhibit diverse ranges of body shapes and motion patterns and are often considered as having superior locomotory performance. The extent to which morphological traits of swimming animals have evolved owing to primarily locomotion considerations is, however, not clear. To shed some light on that question, we present here the optimal shape and motion of undulatory swimming organisms obtained by optimizing locomotive performance measures within the framework of a combined hydrodynamical, structural and novel muscular model. We develop a muscular model for periodic muscle contraction which provides relevant kinematic and energetic quantities required to describe swimming. Using an evolutionary algorithm, we performed a multi-objective optimization for achieving maximum sustained swimming speed *U* and minimum cost of transport (COT)—two conflicting locomotive performance measures that have been conjectured as likely to increase fitness for survival. Starting from an initial population of random characteristics, our results show that, for a range of size scales, fish-like body shapes and motion indeed emerge when *U* and COT are optimized. Inherent boundary-layer-dependent allometric scaling between body mass and kinematic and energetic quantities of the optimal populations is observed. The trade-off between *U* and COT affects the geometry, kinematics and energetics of swimming organisms. Our results are corroborated by empirical data from swimming animals over nine orders of magnitude in size, supporting the notion that optimizing *U* and COT could be the driving force of evolution in many species.

## 1. Introduction

Undulatory swimming organisms achieve locomotory feats that in terms of maximal burst speed, acceleration or agility are unmatched by man-made aquatic vehicles. These have been the inspiration for the development of biomimetic robots [1] which were reverse-engineered based on living fishes under the assumption that their morphology is optimized for swimming. Whether fish-like organisms are indeed optimized for swimming, and whether extant morphological traits would evolve based on locomotion considerations alone, however, has not been completely established.

Despite a vast body of work on various aspects of undulatory swimming (from physiology to physics of swimming), optimization studies based on mathematical models are relatively sparse. Swimming motion for a given body shape has been optimized from a hydrodynamical perspective using theoretical [2] and numerical models [3,4]. The body shape and body stiffness for efficient hydrodynamical performance have also been studied [5]. Further insights into the relationship between morphology and swimming performance have been obtained from numerical studies on bodies of prescribed shape and motion [6–8], but without conducting optimization.

A significant drawback of all these studies is that they do not consider muscle behaviour and the associated energetics so that it is not clear whether the motions or morphologies obtained are physiologically feasible. Mathematical models of muscle behaviour during swimming have mostly been developed to study the muscle response to a given neural activation [9–12], with model parameters being often fine-tuned for particular species [12–14]. The actuation–response relationship varies widely among species [15–17], making these models unsuitable for the optimization of morphological traits of a generic organism. Furthermore, these models do not provide information about metabolic energy consumption, which is a critical component of swimming energetics.

A different approach to studying morphological effects on swimming, without resorting to detailed mechanistic models, is through a comparison of extant morphologies based on their observed performance. Qualitative studies of fish shapes [18–21] and hypothesis testing methods [22,23] have given us some intuition of the optimal body shapes and motion patterns of undulatory swimmers and have provided insight into the effects that performance trade-offs can have on morphology. These methods, however, are qualitative and generally not predictive.

In contrast to the existing studies, our objective is to *predict* optimal morphological traits, including body shapes and motion patterns, across broad ranges of size scales. We achieve this by optimizing locomotory performance measures based on a comprehensive swimming model which incorporates a novel model for periodic muscle contraction. We perform multi-objective optimization with respect to two conflicting performance measures (sustained swimming speed, *U* and cost of transport, COT) to understand the effect of the interplay between them on the morphological traits of the obtained optimal organisms. Finally, the obtained optimal morphological traits are compared with those observed in nature.

## 2. Model description

We study sustained straight-line undulatory swimming (powered by superficial red muscle [20,24,25]), where an organism passes a muscle-produced wave of curvature down its body and propels itself using the hydrodynamic forces exerted on the body as a reaction to the motion. To describe the kinematics and energetics of swimming, the main components of the swimming machine converting the energy from food into useful propulsion work have to be modelled. In addition to an effective and robust body shape and motion description, our swimming model consists of three parts: (i) hydrodynamical model describing the flow around the moving body, (ii) structural model describing the distribution of the internal forces required for swimming motion, and (iii) muscle model describing the muscle behaviour needed to achieve such forces. To facilitate optimization, these model components are sufficiently general to describe the physics for arbitrary morphologies across many scales, and highly computationally efficient to allow a large number of simulation realizations.

### (a) Body and motion description

We consider an arbitrary three-dimensional organism of mass *m* characterized by its body length *L*, tail height *D* and body width *B* (figure 1). We assume that the body is symmetric with respect to the horizontal and vertical planes, with elliptical cross sections of area *A*(*x*) and sectional moment of inertia *I*(*x*). The lengths of axes of cross sections determine the body height and width distributions, denoted by *d*(*x*) and *b*(*x*), respectively. Wetted surface of the body is denoted by *S*. The body is assumed to be neutrally buoyant, with uniform body density *ρ*, for simplicity. Neutrally buoyant fishes often hold the fins close to the body during steady undulatory swimming [15], thus minimizing their effect on the flow around the body. In this paper, we therefore do not consider fins and other appendages.

The locomotory muscle is made of red muscle fibres arranged in a superficial longitudinal strip [20,25], located along the horizontal symmetry plane on each side of the body (figure 1*a*,*b*). The muscle cross section *A*_{m}(*x*) is a small portion *μ*_{0} (*μ*_{0}(*x*) = 2*A*_{m}(*x*)/*A*(*x*)) of the body cross section *A*(*x*) [9,25].

We express the undulatory motion of the body neutral line using a single time harmonic [15,20,24]:
2.1where *ω* is the angular frequency of tail-beat (with tail-beat period *T* = 2*π*/*ω*), *r*(*x*) is the deformation envelope and *λ*_{b} the wavelength of the body undulation.

### (b) Hydrodynamical model

The role of a hydrodynamical model is to determine the relationship between the swimming speed *U* and the tail-beat frequency *ω* for steady swimming, and to provide external forces that occur during swimming. We are interested in swimming at high Reynolds numbers *Re* ≡ *UL*/*ν* (*ν* is the kinematic viscosity of water), for which potential flow models can be used. We use classic Lighthill's potential flow slender-body model for small-amplitude motion [26], which has the advantage of being three-dimensional and very simple to solve compared with other numerical models.

The hydrodynamic pressure field around a freely swimming body gives rise to a forward pointing thrust force *F*_{T} powering the forward motion, and a lateral force *F*_{L}(*x*,*t*) which causes an additional rigid-like lateral movement known as recoil. Both the imposed motion and the recoil are assumed to be small (compared with *L*), so the total deflection *h* of the body can be written as , where *y*_{0}(*t*) is the lateral and *ϕ*(*t*) the angular recoil (figure 1*c*). Equations of motion of a swimming body as a whole, relating the lateral (angular) acceleration and the total external force (moment) acting on the body, provide a way to calculate the unknown lateral (angular) recoil:
2.2where, for a slender body, *F*_{L}(*x*,*t*) = 𝒟(*m*_{a}(*x*)𝒟*h*(*x*,*t*)) [26]. Here 𝒟 ≡ ∂_{t} + *U*∂_{x} is the material derivative and *m*_{a}(*x*) the cross-sectional added mass.

To obtain the steady swimming speed *U* in the present context, we follow a standard approach [4,26,27] wherein one equates the average thrust from a potential flow model with the average drag calculated from an empirical relationship, i.e. requiring . For Lighthill's slender-body model [26], . The drag force is modelled as , using an empirical formula for the drag coefficient *C*_{D} = *C*_{D}(*Re*) (see the electronic supplementary material, equation (S.11); [27]), which exhibits a discrete jump transitioning from laminar to turbulent regime. Although there is some uncertainty about the accuracy of *C*_{D}, proper scaling with *Re* is more important for this study than its exact value.

The solution of the nonlinear system of equations ((2.2), ) determines the steady swimming condition, which can be expressed in terms of *ω*–*U* or *ω*–*Re* relationship since for a given organism *L* is known.

### (c) Structural model

The main purpose of the structural model is to obtain the internal forces acting in a swimming body so that the muscular activity required for powering the motion could be calculated. This is modelled using the standard Euler–Bernoulli beam equation [28]: 2.3

The above terms, corresponding, respectively, to forces owing to inertial, elastic, visco-elastic, hydrodynamic effects, are all balanced by the bending moment *M* produced by muscles. Aggregate Young's modulus *E* and visco-elastic coefficient *ν*_{b} include combined contribution from all the passive elements during bending: elasticity and visco-elasticity of the spine, the skin, the white muscle and the inactive part of red muscles (assuming that the morphology of the organisms is equivalent to that of fishes).

Assuming there are no muscles at the very ends of the body (*M*(*x* = 0, *L*; *t*) = 0), the boundary conditions that a feasible *h*(*x*,*t*) has to satisfy require [28]:
2.4

The sectional bending moment *M*(*x*,*t*) can then be directly obtained from equation (2.3) for a given *h*(*x*,*t*), which satisfies equation (2.4). A muscle model has to be introduced to answer the question how precisely the required bending model *M* is achieved.

### (d) Muscle model

The primary purpose of a muscle model is to determine the physiological feasibility of the prescribed motion and to determine the energy consumption by the muscle, which highly affects swimming energetics (the energy losses in real fish muscles are significant and amount to a muscle efficiency of around 20% [27]). The present model is developed for periodic swimming powered by red muscle, as is generally the case in sustained fish swimming [20,24,25]. Other modes of swimming, e.g. unsteady burst-and-glide swimming in which white muscle fibres are recruited [29], are not considered here. The model should, however, correctly describe the most important characteristics of muscle behaviour and be valid for different undulatory-swimming species and across the scales. We have focused on the facts that seem to be universally valid for swimming fishes and have built a new model based on them.

The contractive force *F*_{musc}(*x*,*t*) that the muscles at some cross section have to provide can be obtained from the calculated required bending moment *M*(*x*, *t*). For a muscle of small cross section placed from the neutral line, this corresponds to *F*_{musc}(*x*, *t*) = *M*(*x*, *t*)/0.5*b*(*x*). Since the muscle produces contractive forces only, in alternating manner from side to side at any *x* [20,30], the sign of *F*_{musc} uniquely determines the side of the active muscle fibres. According to our definition, the required contractive force *F*_{musc} is positive/negative when the muscles on the right/left side of the body are active.

The force *F*_{fib} that each muscle fibre actually produces is a function of the fibre kinematics, which is in turn dependent on body motion. Such dependence is also true for the metabolic power *P*_{fib} consumed per fibre length. During steady swimming, it can be assumed that muscle behaviour is quasi-steady [15,31] since the characteristic time for muscle fibres to adapt to a new force is typically much shorter than the characteristic tail-beat period *T*. Thus, we assume that for a contracting fibre, *F*_{fib} and *P*_{fib} are functions of instantaneous contraction velocity *v*(*x*,*t*), given by Hill's model [31] (see the electronic supplementary material, §3.2).

The contraction velocity *v*(*x*,*t*) of superficial muscle fibres (measured in lengths per second) can be determined from the time rate of change of fibre strain, which in turn can be determined from the curvature of the neutral line alone [9,15,24,32]. Based on a simple beam theory [15,24],
2.5

The sign of *F*_{musc} determines the choice of plus–minus sign in equation (2.5), where plus(minus) corresponds to the case when the fibres on the right(left) side of the body are active (the active side of the body cannot be determined from the rate of change of curvature of the spine alone).

Non-dimensional relative contraction velocity is defined as *v*_{r}(*x*,*t*) ≡ *v*(*x*,*t*)/*v*_{max}, where *v*_{max} is maximal achievable contraction velocity for given fibre characteristics.

At any cross section, the required muscle force *F*_{musc} is the sum of all the active single-fibre contractive forces *F*_{fib}. To obtain the required force *F*_{musc}(*x*,*t*) constrained by *F*_{fib}(*v*(*x*,*t*)), we assume that only a fraction *μ*(*x*,*t*) of the total muscle cross-sectional area *A*_{m}(*x*) is activated:
2.6

The condition for a physiologically feasible motion *h*(*x*,*t*) can then be stated as:
2.7The metabolic power consumption per unit length of the muscle *P*_{musc}(*x*,*t*) = *A*_{m}(*x*)|*μ*(*x*,*t*)|*P*_{fib}(*v*(*x*,*t*)) is proportional to the active muscle portion. It is always positive, corresponding to the fact that metabolic energy is being spent when the mechanical power output of the muscle *P*_{mech}(*x*,*t*) = *F*_{musc}(*x*,*t*)*v*(*x*,*t*) is positive or negative, regardlessly. With the muscle force and power consumption calculated, all relevant dynamic and energetic quantities for locomotion can be calculated. The predicted muscle efficiency matches the measured one for swimming fishes and for isolated red fibres (see the electronic supplementary material, §3.3).

## 3. Performance measures and optimization variables

The optimization problem we are trying to solve can be stated as follows: find optimal solutions for a set of conflicting objectives (locomotory performance measures) over the variables that adequately parametrize the body shape and motion, constrained by the motion feasibility, equation (2.7) and shape integrity conditions. Body shape and motion parameters are chosen as the optimization variables since they are the key mechanistic components that determine locomotory performance.

To elucidate the trade-offs between conflicting locomotion-based objectives, we focus on two performance measures of arguably great importance in the evolutionary scenario [20,33]: maximizing sustained swimming speed and minimizing energy consumption. For the latter, we use a standard non-dimensional measure called COT [20,34] (for derivation, see the electronic supplementary material, §5):
3.1where , the total metabolic power consumed by swimming at speed *U*, is the sum of the metabolic power consumed by swimming muscles and the standard metabolic rate *P*_{s} required for other physiological processes even when there is no motion at all ( denotes a length-integrated, time-averaged quantity). Note that in equation (3.1), gravity *g* is used merely for non-dimensionalization and is not related to swimming. Expressed by equation (3.1), COT is the ‘gallons-per-mile’ measure quantifying the total energy consumption per unit mass and distance, which probably governs long migrations [20].

The choice of locomotive performance measures to optimize is not unique. For example, an energetic measure can be a generic power coefficient defined as , where *P*_{0} ≡ 0.5*ρ**SU*^{3}, and is some measure of swimming power based on which *C*_{P} has different meanings and implications. In general, *C*_{P} might be more suited for studying the efficacy of hydrodynamical propulsion itself as it is normalized by the scale of hydromechanical power *P*_{0}. The ultimate justification of the present choice of *U* and COT has to be borne out on whether the consequent predictions based on it are corroborated by nature.

Optimizing conflicting objectives usually leads to an infinite number of optimal solutions. Since by the definition of conflicting objectives an organism cannot be optimal in every objective; it is considered as optimal when it is non-dominated [35], i.e. when there is no (feasible) variation of an organism's morphology that could improve every objective. We call the set of non-dominated organisms the optimal population *Π*.

To facilitate the optimization of generic swimming geometries and motions, we parametrize the body height, width and motion along the body in terms of general unbiased mathematical descriptions. We represent the body height distribution *d*(*x*)/*L* by a sum of *N*_{S} + 1 polynomial shape functions *D*_{n},
3.2where *T*_{n}(*x*) is the Chebyshev polynomial of the first kind of order *n*; shape coefficients *C*_{n} produce different shapes when varied. Without loss of generality, we assume the body width *b*(*x*) to be given by a symmetrical NACA-00 profile with relative maximum thickness *B*/*L*. We thus parametrize the body shape by *N*^{S} = *N*_{S} + 3 optimization variables (*D*/*L*, *B*/*L*, *C*_{0}, … ,*C*_{NS+1}). The body length *L* is not a parameter as it can be calculated for a given *m* once *d*(*x*)/*L* and *b*(*x*)/*L* are prescribed.

The spatial and temporal parametrization of body motion is achieved using *N*^{M} = *N* + 1 variables. The envelope *r*(*x*) is represented as a sum of *N* Chebyshev polynomials, where the coefficients of the series serve as optimization variables (see the electronic supplementary material, §1.2). To reduce the number of optimization variables and to ensure the validity of Lighthill's model (see the electronic supplementary material, §4), we set the relative body-undulation wavelength to *λ*_{b}/*L* = 1, a value characteristic for many fishes [15,20,27]. Upon parametrization, the motion is slightly corrected to satisfy motion boundary conditions (2.4).

The swimming speed *U* and the tail-beat period *T* can both be determined from *Re* using the steady swimming condition. Hence, we use *Re* as a kinematic optimization variable and the values of *U* and *T* (or *ω*) are determined as the outcome of optimization.

In the following, we use *N*^{S} = *N*_{S} + 3 = 5 and *N*^{M} = *N* + 1 = 4 as we have found that those values are sufficient to represent the extant body shapes and motion patterns to within *O*(1%). The advantage of our parametrization is that, despite *N*^{S} and *N*^{M} being small, we are capable of representing a large variety of shapes and motion patterns without introducing a particular bias.

## 4. Results

We optimize for *U* and COT using a multi-objective evolutionary algorithm [36], evolving generations of feasible populations starting from the one with random body shape and motion parameters. We perform calculations for body sizes ranging from *m* = 0.001 kg to *m* = 1 000 000 kg to obtain the optimal populations *Π*(*m*), figure 2. Given the conflicting nature of optimization objectives, *Π*(*m*) obtains a range of values for each swimming characteristic presented (*Re*, *U*, COT, *T*, relative tail amplitude *h*_{T}/*L*, *μ*_{max}). The results are compared with the empirical data, where available, for fishes and cetaceans.

For specificity, in this discussion, we focus on the values attained by organisms for which either *U* or COT is optimal. Hereafter, these predicted values are denoted as (·)_{U-opt} and (·)_{COT-opt} for those corresponding to *U*- and COT-optimal organisms, respectively. As discussed earlier, the choice of performance measure for optimization is not unique. For comparison, we provide results for the minimization of power coefficients, namely of power output-based (used in Tytell *et al.* [14]) and of power consumption-based (suggested in Yates [39]).

The Reynolds number *Re* employed by the optimal populations *Π* grows over four orders of magnitude (figure 2*a*). The prominent feature of *Re*–*m* relationship is the presence of a transition region ℛ_{T} separating otherwise allometric relationships (visible from the linear (*Re*)_{U-opt}–*m* or (*Re*)_{COT-opt}–*m* relationships in log–log plots). The transition regions (ℛ_{T})_{U-opt} and (ℛ_{T})_{COT-opt} are defined as the range of *m* for which *U*-optimal and COT-optimal organisms swim at speeds just below critical Reynolds number *Re*_{cr} to remain in the laminar regime. The ranges of (ℛ_{T})_{U-opt} and (ℛ_{T})_{COT-opt} differ, reflecting the earlier transition to turbulent flow of *U*-optimal organisms ((*Re*)_{U-opt} > (*Re*)_{COT-opt} for a given *m*). Different behaviour in ℛ_{T}, accompanied with the change of (·)_{U-opt}–*m* and (·)_{COT-opt}–*m* slopes over it, is a common feature of almost all quantities describing *Π*(*m*) (some shown in figure 2).

The optimized swimming speed *U* obtains values from *O*(0.1 − 1) m s^{−1} (corresponding to relative swimming speed *U*/*L* in body-lengths per second from *O*(1) to *O*(0.1); electronic supplementary material, figure S12*a*). As one of the performance measures being optimized, *U* is clearly maximal(minimal) for *U*-optimal(COT-optimal) organisms in *Π* of all body sizes, as expected. We find a decrease in slopes of (*U*)_{U-opt}–*m* and (*U*)_{COT-opt}–*m* over ℛ_{T}, as has been previously suggested [21,37] (figure 2*b*). The slight decrease of (*U*)_{U-opt} in the transition region (ℛ_{T})_{U-opt} is owing to the organism's inability to cross the laminar-to-turbulent transition with the available muscle. A similar, but more pronounced decrease of (*U*)_{COT-opt} in the transition region (ℛ_{T})_{COT-opt} can be explained by energetic arguments: here more muscle units could be employed but that would result in undesirably higher COT.

The COT (figure 2*c*) is one of the quantities that heretofore could not be predicted from theoretical or numerical considerations owing to the lack of a comprehensive muscle model. The results we obtain show a slight general under-prediction of the COT range which might imply that the values of *P*_{s} or *ν*_{b} we use might be lower than those in many natural organisms.

The obtained tail-beat period *T* in the laminar regime seems to be slightly greater than the measured one (but of the same order; figure 2*d*). Over the entire *m*-range, (*T*)_{COT-opt} > (*T*)_{U-opt} consistently. We find that the increase in *T* with *m* is correlated with the decrease in maximum max *v*_{r}, as has also been empirically found for cyclical muscle contractions [40] (electronic supplementary material, figure S12*c*). Note that even for the smallest organisms investigated, *T* > 0.1s (figure 2*d*) which is greater than the 30–50 ms needed for the muscle fibre to adapt to a new force [15], thus not violating the quasi-steady assumption.

The relative tail amplitude *h*_{T}/*L* shows a decreasing trend in each of the sub-regions (figure 2*e*). Generally, *h*_{T}/*L* ≪ 1, which does not violate our small-amplitude motion assumption. A non-obvious prediction is the fact that (*h*_{T}/*L*)_{U-opt} < (*h*_{T}/*L*)_{COT-opt} below ℛ_{T} but (*h*_{T}/*L*)_{U-opt} > (*h*_{T}/*L*)_{COT-opt} above ℛ_{T}.

The maximum active muscle portion *μ*_{max}, i.e. the maximum value of *μ*(*x*,*t*), exposes some of the driving constraints behind the obtained overall results (figure 2*f*). As expected, (*μ*_{max})_{U-opt} = 1 for all *m*, limiting the maximal achievable swimming speed. Generally, (*μ*_{max})_{COT-opt} < 1 indicating that only a portion of the muscles is required, as suggested [41].

Compared with these results, it appears that power output-based power coefficient *C*_{P}^{M} is not an adequate objective function as its predictions deviate from empirical data by several orders of magnitude for larger *m* (figure 2*a*–*e*). On the other hand, owing to the presence of muscle-consumed power, optimizing power consumption-based power coefficient *C*_{P}^{T} gives reasonable results (similar to optimizing *U* or COT), suggesting that other reasonable measures could be at play in living organisms.

The optimal motion envelopes *r*(*x*) converge to fish-like motion envelopes, figure 3 (cf. *r*(*x*) of initial population; electronic supplementary material, figure S11). We show here (*r*(*x*))_{U-opt} and (*r*(*x*))_{COT-opt} for select *m*; the envelopes within each optimal population *Π*(*m*) and with the change in *m* vary smoothly between those presented. Considering (*r*(*x*))_{COT-opt}, the motion is largely confined to the aft part of the body which, together with *λ*′_{b} = 1, consistently resembles the (sub)carangiform motion [15,20] (the terminology is not uniform in the literature [27,42]). Carangiform swimming has previously been associated with low energetic costs [21]. Interestingly, a computational fluid dynamics (CFD) study of mackerel and lamprey swimming [8] found that at high Reynolds numbers, the (sub)carangiform motion is faster than prescribed anguilliform motion. However, a direct comparison with our results (after matching *Re* and *λ*_{b}/*L*) is not easy since a muscle model is not considered in Borazjani & Sotiropoulos [8] and so it is not clear whether the prescribed motion is physiologically feasible (see the electronic supplementary material, §4 for details). Very small amplitudes of (*r*(*x*))_{COT-opt} in ℛ_{T} (cf. figure 2*e*) are in line with the decrease in (*U*)_{COT-opt}. It is, however, possible that Lighthill's theory together with the *C*_{D}(*Re*) model provide less accurate results in the boundary layer transition region *R*_{T}. We note that in some cases, there is significant motion of the head. This less-than-intuitive kinematics is a limitation of the present body model wherein the muscle actuation extends throughout the fish body, including the head.

The changes in kinematic and energetic quantities across the scales and among performance measures are accompanied by the shape modifications of optimal organisms (figure 4; also electronic supplementary material, figure S12*d–f*). Relative to fishes found in nature, the shapes show qualitative resemblances, for example, the emergence of the caudal peduncle that is more pronounced for COT-optimal organisms in the range *m* = 1∼100 kg. Over the transition region ℛ_{T} in the middle, optimal organisms have generally smaller *U* and *L* than the allometric expectation as they try to remain in the laminar regime. Such adaptations might be observed in nature with organisms that swim predominantly near *Re*_{cr}. Although the shapes are mostly slender, lateral dimensions *D* and *B* reach 0.4*L* in some cases (electronic supplementary material, figure S12*e,f*), where Lighthill's theory may cease to be valid. Results corresponding to these shapes should be considered with care.

## 5. Discussion

This study shows, using relatively standard hydrodynamic and structural descriptions and a novel muscular model, how optimal undulatory swimming organisms might look and move if the driving force behind evolution were locomotory performance measures, in particular, the swimming speed *U* and the COT. If submodels of different complexity or different performance objectives are used, the overall optimization framework should still be useful, although the detailed predictions would of course vary.

The body shape in nature primarily affects the hydrodynamics of swimming (in our model, it also influences the muscle performance through muscle disposition). The effect of shape on hydrodynamics in Lighthill's model is accounted for by the recoil equations (2.2), wherein the overall shape affects the total deflection *h*(*x*,*t*). Lighthill's model exhibits limitations, however. The hydrodynamics at very long motion wavelengths *λ*_{b} ≫ *L* is not correctly captured (see the electronic supplementary material, §4), therefore, a constraint on the value of *λ*_{b}/*L* is required. The model also neglects vortex shedding, lateral flow separation and viscous drag (relevant at lower *Re* numbers [13]). Despite these restrictions, Lighthill's model has been shown to provide sufficiently accurate values for the obtained lateral force [14,43]. It is important to point out, however, that our primary interest is in the correct scaling of quantities with *Re* and the proper dependence on kinematic and geometric parameters, rather than in the quantitative accuracy (requiring substantially greater computational cost). For example, we have compared the scaling of the stride length with *Re* calculated by Lighthill's model and empirical drag formula with that from a more sophisticated hydrodynamic model [7]. Over the wide range of *Re*, the slopes of the predicted scaling agree to within approximately 10 per cent.

Modelling hydrodynamics with higher accuracy might be achieved at low Reynolds numbers *Re* (*O*(10^{3}–10^{4})) where CFD models solving the viscous flow equations [3,7] are computationally feasible. However, the above *Re*-range covers only a small range of *Re* considered in this paper (which basically covers the entire range of fish and cetacean swimming). For such large *Re* numbers, potential flow models [6,26,43,44] are often the only option. The large numbers of simulation evaluations required (*O*(10^{7}) for this study) further limits the computationally feasible models to only the simplest ones. Lighthill's model provides a reasonable choice because it is valid for large Reynolds numbers *Re* and it is computationally very efficient.

Regardless of the complexity of the hydrodynamical model used, it alone cannot account for the losses that occur during the conversion of metabolic energy from food to useful mechanical work, nor can it assess the physiological feasibility of the prescribed motion, both of which are of a fundamental concern. For these reasons, the addition of the muscle behaviour model is absolutely necessary if the overall swimming physics is to be considered.

We have introduced a novel muscle model primarily because the existing models for muscle behaviour during swimming do not provide the metabolic power consumption information [9–12]. Our model of muscle behaviour considers the contraction velocity *v*(*x*,*t*) and the required contraction force *F*_{musc}(*x*,*t*) as primary quantities, which avoids relying on a still uncertain and variable relationship between *F*_{musc} and neural activity [15–17] as was done in previous studies [9–12]. The fact that the feasible combinations of the imposed motion and the required *F*_{musc}(*x*,*t*) are determined by the available muscle and the intrinsic properties of muscle fibres is often overlooked in studies which only consider the hydrodynamical aspect of swimming [3,4,7,8].

Our results compare reasonably favourably across many scales, which lends validity to the present overall model, despite the assumptions and simplifications therein. None of the quantities presented in §4 are *prescribed*; they are all outcomes of the optimization procedure, i.e. our results give the values the optimal organisms would choose to employ. As such, our results are fundamentally different from previous studies where a kinematic quantity (either *Re* [7,8] or neural activation [13,14]) that directly sets the swimming speed was prescribed. We limit motion of the wavelength to *λ*_{b}/*L* = 1, but that is a restriction on the degrees of freedom by which we describe the motion, not on a parameter that drives the motion. The value *λ*_{b}/*L* = 1 is roughly between those characteristic for the anguilliform and the carangiform swimming and is used by many fish species [15,20,27]. With such a choice, motion and geometry features of both swimming forms are found in optimal population *Π*. However, organisms with *λ*_{b}/*L* significantly different from 1, like lamprey or scup [15], or even ‘anguilliform mackerel’ [8], cannot be correctly modelled with the present model. Relaxing the constraint on *λ*_{b}/*L*, which is undoubtedly an important parameter for swimming, would further enrich this study.

The lack of artificially imposed constraints enables us to obtain the intrinsic scaling of kinematic and energetic quantities as it emerges from the optimization. Inherent allometric relationships (based on body mass *m*) are found for some quantities and they exhibit boundary layer regime dependence. Such scaling results have heretofore not been predicted from theoretical/numerical considerations alone. Discrepancies between the measured and predicted values might result from the likelihood that some measured values have not been obtained under the sustained swimming regime we assume, or that values of actual muscle and tissue properties differ from those we use. Improved predictions could presumably be achieved by tailoring the model parameters to a particular species (e.g. geometry, muscle properties and distribution); although uncertainty in measured data still remains, especially for larger *m*. Larger deviations might also indicate that other swimming or muscle behaviour not modelled here, or performance measures not presently considered, are involved.

Realistic overall results (figures 2–4) make it interesting to make a direct quantitative comparison between predicted shapes and kinematics of optimal organisms and select aquatic species over a range of *m* (figure 5). Despite the relative simplicity of the present model, including the low degrees-of-freedom in the modal representations of the shape and motion, we obtain a good match. The conflicting nature of optimizing COT and *U* contributes to the diversity of the obtained morphologies and behaviour. While parallels could be drawn between the performances of the real organisms and the theoretically predicted ones (e.g. the tuna-similar organism being close to COT-optimal—a feat for which tuna is often noted), the intent here is primarily to show that swimmers in the real world do exhibit rather similar characteristics to those predicted. In spite of a possible bias in the selection of the specific examples, the overall corroboration of the model predictions by swimming animals in nature for diverse measures and across the scales is noteworthy.

While locomotion-based performance measures studied here are not necessarily the (only) important ones in nature, the present study provides a direct evidence of their impact on morphology. Comparisons of model-predicted morphological traits and those of real organisms also provide some means for deducing possible roles that specific performance measures might have played (causation) in the organisms' adaptation. A further insight into understanding the diversity of extant morphologies could be achieved by varying the choice of performance objectives and studying the predicted morphologies, based on the present optimization framework. Understanding whether and how living morphologies are related to specific performance measures would also pave the way for improved biomimetic swimming vehicles.

## Acknowledgements

This study is financially supported by the US Office of Naval Research.

## Appendix A

#### A.1. Assumed body/muscle/fluid properties

For simplicity, in all our calculations muscle and tissue properties are taken as length and size independent, but characteristic for fishes (red fibre isometric force *F*_{0} = 150 kN m^{−2}, *v*_{max} = 5 lengths s^{−1} [48], *E* = 10^{5} N m^{−2}, *ν*_{b} = 10^{4} m^{2} s^{−1} [9,14,28], *μ*_{0} = 0.1 [25]). The standard metabolic rate used here is *P*_{s} = 0.1327 *m*^{0.80} [W] [49]. Fresh water properties are used throughout (*ρ* = 10^{3} kg m^{−3}, *ν* = 10^{−6} m^{2} s^{−1}).

#### A.2. Optimization algorithm

The optimization is conducted for organisms of mass *m* = *a*10^{b}, with log*a* = 0,1/4,1/2,3/4 and *b* = −3, … ,6. We use a multi-objective covariance matrix adaptation evolutionary strategy [36], with default parameters. For every case, an initial randomly generated feasible population of 500 individuals is evolved through 500 generations. The optimization converges in all cases, and the bounds imposed on the variables are never active in the final population.

#### A.3. Shape similarity measure

We define the shape similarity measure *𝒮* as:

It is bounded from above by 1, which marks a perfect similarity in shape. Here, *x* represents longitudinal coordinate normalized by the respective organism length *L*, such that both *d*(*x*) and the height distribution of living organisms *d*_{r}(*x*) (omitting fins and other appendages) are defined on *x*∈ [0,1].

- Received January 10, 2012.
- Accepted March 9, 2012.

- This journal is © 2012 The Royal Society