Optimal movement in the prey strikes of weakly electric fish: a case study of the interplay of body plan and movement capability

Claire M Postlethwaite, Tiffany M Psemeneki, Jangir Selimkhanov, Mary Silber, Malcolm A MacIver

Abstract

Animal behaviour arises through a complex mixture of biomechanical, neuronal, sensory and control constraints. By focusing on a simple, stereotyped movement, the prey capture strike of a weakly electric fish, we show that the trajectory of a strike is one which minimizes effort. Specifically, we model the fish as a rigid ellipsoid moving through a fluid with no viscosity, governed by Kirchhoff's equations. This formulation allows us to exploit methods of discrete mechanics and optimal control to compute idealized fish trajectories that minimize a cost function. We compare these with the measured prey capture strikes of weakly electric fish from a previous study. The fish has certain movement limitations that are not incorporated in the mathematical model, such as not being able to move sideways. Nonetheless, we show quantitatively that the computed least-cost trajectories are remarkably similar to the measured trajectories. Since, in this simplified model, the basic geometry of the idealized fish determines the favourable modes of movement, this suggests a high degree of influence between body shape and movement capability. Simplified minimal models and optimization methods can give significant insight into how body morphology and movement capability are closely attuned in fish locomotion.

Keywords:

1. Introduction

Understanding the dynamics of behaviour is challenging because movement integrates multiple simultaneous processes and constraints, such as the mechanical and physical features of the animal and its habitat. The interplay of these aspects of the problem is complex and recommends quantitative methods for their understanding. Under the assumption that specific behavioural epochs can be understood as a solution of a constrained optimization problem (Kern & Koumoutsakos 2006; Srinivasan & Ruina 2006; Tam & Hosoi 2007), the impact of an array of neuromechanical constraints that shape the behaviour can be investigated. For example, a movement may be hypothesized to minimize the required force or time. The actual performance of the system is then compared against a model that generates optimal behaviour. If the real and optimal behaviours are similar, this is consistent with the biological system being configured to be optimal in the hypothesized respect. For this approach to be successful, the behaviour should be simple, stereotyped and well characterized. The prey capture behaviour of the weakly electric black ghost knifefish (Apteronotus albifrons) is one such example.

Black ghost knifefish continually emit a rapidly oscillating (approx. 1 kHz) weak electric field (approx. 1 mV cm−1 near the body). Perturbations of this field are sensed by over 10 000 sensory receptors scattered over the entire body surface. These perturbations are created whenever something that differs from the electrical properties of the surrounding water enters the field. This unique mode of sensing, termed active electrosense, allows weakly electric fish to hunt at night, in the muddy rivers of the Amazon where vision is rendered useless. Active electrosense in these animals has become a leading system for investigations into how vertebrates process sensory information (reviews: Bullock & Heiligenberg 1986; Turner et al. 1999), and, more recently, how sensory processing relates to mechanics (Cowan & Fortune 2007; Snyder et al. 2007). Figure 1 illustrates the fish body plan.

Figure 1

Body plan of fish and ellipsoidal model. The figure shows the outline of the ellipsoid model and a transparent body plan of the black ghost knifefish Apteronotus albifrons. Fish lengths varied from 12 to 15 cm. The length of the ellipsoid (2l1) is scaled such that the volumes of the fish and the ellipsoid match. Note that the ellipsoid is narrower than the fish at the front of the body, and wider than the fish at the rear of the body, with the consequent changes in the shading of the ellipsoid as it is covered by the surface of the fish body. (a) The unbent fish and fitted volume-matched ellipsoid. (b) How the ellipsoid is oriented if the fish body is bent. The tip of the ellipsoid coincides with the location of the fish nose, and is oriented to coincide with the unbent rostral third of the fish.

Prior work has shown that the volume in which prey are detected (approx. 3 cm out in all directions from the body surface) is similar in size and shape to the volume which the animal needs to come to a halt to capture the prey (Snyder et al. 2007). A rapid and precise strike must therefore be initiated almost immediately upon detection—otherwise, the prey will rapidly pass out of the sensory volume. Given this, and under the assumption that energy cost is proportional to muscle work, we expect that strikes will be strongly constrained by mechanical factors and that optimization of an appropriately chosen mechanical utility function may be predictive of the behaviour.

The kinematics of the weakly electric fish's prey strike has been systematically examined empirically using an approach that gives the full three-dimensional position of the body over time (MacIver et al. 2001). The results show the behaviour to be highly stereotyped. The fish generally swims forward at approximately 10 cm s−1, in the dark, using its electrosensory system to detect prey in all directions around the body up to approximately 3 cm away. Only a small fraction of the prey are detected ahead of the body. In the majority of cases, the fish rapidly reverses its body (in ≈100 ms) to bring its mouth to the detected prey.

In our highly abstracted approach, we do not model the way the fish generates propulsion. Propulsion in aquatic animals generally occurs through flapping of fins or body undulation. These movements impart momentum to the fluid, and are generally associated with the development of characteristic vortices in the wake of the swimmer. This has been observed in both experiments and numerical simulations. A few examples include vortices shed at the trailing edges of anguilliform swimming (i.e. eels; Tytell & Lauder 2004), vortex rings shed by jellyfish (Dabiri et al. 2005a,b) and knifefish (Shirgaonkar et al. in press), and vortex structures generated by the pectoral fins of sunfish (Drucker & Lauder 1999). Weakly electric fish swim using a unique adaptation, an elongated ribbon-like anal fin along the ventral midline of the body. Rather than bend their body for propulsion as many fishes do, these fish generate travelling waves along their ribbon fin, often keeping their body straight (Blake 1983; Lighthill & Blake 1990). By changing the direction, frequency, amplitude and number of waves on the fin, they can precisely control the magnitude and direction of thrust (Shirgaonkar et al. in press). During exploration of novel objects or in confined spaces, body bends can occur to facilitate sensing or for turning, and this bending is decoupled from propulsion (Assad et al. 1999). The prey strikes we are modelling are high-speed transient movements. During these high-speed movements, just as in steady swimming, the body is kept straight—or if it was bent prior to prey detection, it is rapidly straightened (MacIver et al. 2001), presumably to maximize the surge force arising from the streamwise jet generated by the ribbon fin (Shirgaonkar et al. in press).

Given these considerations, as well as significant mathematical advantages, we model the fish body as a rigid ellipsoid with no fins—which is clearly unable to generate the types of forces described earlier, which are required for locomotion. In our model, the fish is propelled by forces and torques which act at the centre of mass of the body, but we do not describe how these forces originate. Our model does not consider the complicated fluid mechanics underlying the propulsion methods discussed above, but assumes only that the fish has some mechanism of generating thrust, and that the generation of this thrust has some cost to the fish.

As with most fishes, the motion capabilities of the knifefish are limited. More precisely, if a rigid body can actively move in all six degrees of freedom (surge, heave and sway translational motions, and yaw, pitch and roll rotational motions—see figure 2), it is said to be ‘fully actuated’, while if it cannot it is ‘underactuated’. As we will provide evidence for below, the knifefish can move forward, backward, heave upward, roll and pitch—but cannot yaw, sway or heave downward. Thus, it is underactuated.

Figure 2

Velocity vector representation. We show the representation for translational and rotational velocities simultaneously, with labels related to the rotational velocity representation shown in parentheses. (a) Representing the directions of the linear (v) and rotational (ω) velocity vectors as points on the surface of a sphere. The six shaded circles are regions close to a coordinate axis. If a vector intersects the sphere in one of the circles, it can be associated with a simple movement in the direction of or about that axis. A linear velocity vector v corresponds to a translation in the direction of that vector. A rotational velocity vector ω corresponds to a rotation of the body about that axis. We show an example vector, denoted by v (ω). (b) The spherical angles Θ and Φ (defined in §2.1.1) as rectangular coordinates. The shaded regions correspond to those in (a) and are labelled according to the type of translational (or, in parentheses, rotational) movement that is associated with the velocity vectors within those regions. The indicated point v (ω) illustrates how a vector is mapped from the sphere shown in (a).

In this paper, we test the hypothesis that the trajectory of the fish body, as approximated by an appropriately sized rigid ellipsoid, minimizes mechanical effort. Importantly, we allow our idealized fish to be fully actuated. If the fully actuated idealized body follows the same trajectory as the underactuated real fish, this would be indirect evidence that the fish, even though underactuated, is sufficiently actuated to perform optimal movements. Another way to express the point is that the fish has all the movement capability it needs for effective low-cost movement. To provide more direct evidence for this claim, we also underactuate the idealized fish in some of the same ways that the real fish is underactuated, and examine whether there is any change in the generated optimal trajectories. If the optimal trajectories are essentially unchanged, then we can conclude that those additional degrees of freedom are not needed by the fish to achieve optimal movements.

This paper is organized as follows. In §2 we describe our simplified mechanical model of the fish and our optimization procedure. In §3 we apply our optimization methods to test three hypotheses. We find that our optimal trajectories match the measured trajectories with excellent agreement, suggesting that the fish indeed moves so as to minimize effort. This is discussed in §3.1. Further in §3.2 we find that optimal trajectories for a fish model that was underactuated in the same way as the real fish have little difference from the fully actuated case. The fish therefore has all the movement abilities it needs for low-cost movements. Finally in §3.3 we show that we can use the minimum-effort criterion to detect a switch from prey search behaviour to the goal-directed behaviour of a prey strike in a post hoc analysis. This is simply because, prior to prey detection, the fish is not attempting to reach a specific point in space (the one occupied by the prey) with minimal effort; we can use the onset of an optimal trajectory to that point as a method for detecting a behavioural mode switch. Section 4 concludes the paper with a discussion of the results, the limitations of our mechanical model and some suggestions of how we might further investigate the motion and sensing capabilities of the weakly electric fish using tools of mechanics and optimal control.

2. Methods

In a previously published study (MacIver et al. 2001), adult weakly electric fish (A. albifrons) were videotaped in a light-tight enclosure under infrared illumination. Individual water fleas (Daphnia magna, 2–3 mm in length) were introduced near the water surface and drifted downward; prey capture behaviour was recorded using a pair of video cameras oriented along orthogonal axes. Relative to the fish's velocity (approx. 10 cm s−1), the prey were relatively stationary (prey velocity less than 2 cm s−1). Prey capture events (from shortly before detection to capture) were subsequently digitized, and three-dimensional motion trajectories of the fish surface and prey were obtained using a model-based tracking system with a spatial resolution of 0.5 mm and a temporal resolution of 1/60 s (MacIver & Nelson 2000). The time of prey detection (tD) was defined by the onset of an abrupt longitudinal deceleration as the fish reversed swimming direction to capture the prey. These reversals are characteristic of most prey capture encounters. Initial prey sensing tends to be uniformly distributed along the entire length of the body, so a reversal of swimming direction is typically required to intercept the prey. A total of 116 prey capture manoeuvres (hereafter, ‘trials’) were recorded.

In our optimization procedure, we use an idealized model of the fish in which we model the body as a solid rigid ellipsoid. In order to compare mechanically optimal trajectories with the video-captured fish trajectories described above, we fit a trajectory of an ellipsoid to the video-captured data. Four fish of slightly different sizes were used in the prey capture study (MacIver et al. 2001). For each prey capture trial, the corresponding model ellipsoid has the same width (2l2) and height (2l3) as the fish, but the length (2l1) of the ellipsoid is chosen such that the volumes of the fish and ellipsoid are equal. This results in an ellipsoid with approximately 67 per cent of the length of the fish. We assume the fish is neutrally buoyant, so the density of the ellipsoid is equal to that of the surrounding fluid. Equal volumes of the fish and ellipsoid then also result in equal masses. For the four fish that were used in the prey capture study, the model ellipsoids have semi-axes lengths l1 ranging from 4.06 to 5.05 cm, l2 from 0.39 to 0.46 cm and l3 from 0.85 to 1.08 cm.

The data from the digitized motion capture events give the coordinates of the fish nose in a laboratory-fixed frame, and three angles that describe the orientation of the fish body, as functions of time. We create an approximation to the fish's trajectory by fitting the ellipsoid model to the trajectory of the fish. For each time step of the recorded behaviour, we position the ellipsoid so that the front of the ellipsoid and the fish nose coincide, and orient the ellipsoid so that the rostral part of the fish and the ellipsoid are aligned, as shown in figure 1. Throughout the paper, we refer to this rigid ellipsoid-fitted trajectory as the ‘motion-capture fitted (MCF) trajectory.’

In comparisons between mechanically optimal and MCF trajectories, each of the two ellipsoids is represented by a surface mesh of 256 nodes. We evaluate the distance between corresponding nodes to obtain a measure of how similar the two trajectories are. Details of this error measure calculation can be found in §2.4.

We next describe how the optimal trajectories are computed, beginning with a description of how we model the hydrodynamics of the submerged ellipsoidal approximation to the fish. We then describe the optimization procedure, the objective function that is to be minimized and the constraints on the optimization.

2.1 Equations of motion

The fish is modelled as a neutrally buoyant solid rigid ellipsoid immersed in an infinitely large volume of incompressible, inviscid fluid that is at rest at infinity. For Reynolds numbers typical of the fish motions (103–104) (Blake 1983; Shirgaonkar et al. in press), the inviscid assumption is reasonable. In some trials, owing to the small behavioural tank the original study was conducted in, the fish is bent prior to prey detection. As discussed in §1, this bending is not part of the fish's propulsion mechanism. We discuss the bending of the fish and how we exclude such trials from our comparison in §2.3.

It is well known (Lamb 1932) that the motion of such an ellipsoid through an ideal fluid can be described using the Kirchhoff equationsEmbedded Image(2.1)Embedded Image(2.2)where v and ω are three-vectors describing the velocity of the centre of mass of the ellipsoid and the angular velocity of the body, respectively. The matrices M and J are the effective mass and moment of inertia (discussed further below), and F and T are applied forces and torques.

Recall that we do not discuss the origin of the applied forces and torques, but only assume that they are generated by the fish, and producing forces in different directions costs the fish the same amount of effort. We can write Embedded Image and Embedded Image where F1 is the surge force, F2 is the sway force and F3 is the heave force, and T1 is the roll torque, T2 is the pitch torque and T3 is the yaw torque. In this formulation, the forces and torques act at the centre of mass of the body. Moreover, we allow the body to be fully actuated, i.e. F(t) and T(t) can be any three-vectors.

Since we assume the fluid is inviscid, this model does not include the effects of drag. We discuss reasons why this is not expected to significantly affect our results in §4.

All quantities in equations (2.1) and (2.2) are given in a frame of reference which rotates with the body. We refer to this frame as the ‘body frame’, and to a fixed inertial frame of reference as the ‘laboratory frame’. We choose the body frame so the x-axis is aligned with the long axis of the fish, as shown in figure 1. The matrices M and J are given byEmbedded Image(2.3)where M is the sum of the added mass matrix (due to the volume of fluid accelerated by translations of the ellipsoid) and the body mass matrix, and J is the sum of the added moment of inertia matrix (due to the volume of fluid accelerated by rotations of the ellipsoid) and the body moment of inertia matrix. Formulae to calculate mi and ji can be found in Holmes et al. (1998). Recall that we use four model ellipsoids of slightly different sizes in our calculations. For one of our four ellipsoid models, we find (m1, m2, m3)=(6.04, 17.31, 8.39) g, (j1, j2, j3)=(1.57, 27.78, 54.11) g cm2 and the mass of the ellipsoid, m=5.84 g. The other three ellipsoids have similar values. The fluid motion enters the model only through the added mass and inertia effects.

The position and orientation of the body with respect to the laboratory frame are given by a three-vector b, which gives the position of the centre of mass of the body, and a rotation matrix R, which we write in terms of three Euler angles ϕ, θ and ψ (R is given in terms of these angles in appendix A). We choose the Euler angles so that ϕ is a rotation about the x-axis (roll), θ is a rotation about the y-axis (pitch) and ψ is a rotation about the z-axis (yaw). Note that the Euler angles themselves do not give useful information for visualizing the orientation of the fish, since if more than one of the Euler angles is non-zero the fish has undergone a combination of rotations. Rotations of the fish are better understood by examining the form of the angular velocity vector ω. We discuss this further along with our method of representing the rotations of the fish in §2.1.1. The choice of the Euler angles we make has singularities at θ=±90°, and so we restrict θ∈(−90°+ϵ, 90°−ϵ) for some 0<ϵ≪180°. (Typically we choose ϵ≈6°. No restrictions are required on the other angles.)

The body-frame velocity v and angular velocity ω are related to the state variables in Embedded Image and R byEmbedded Image(2.4)where if ω=(ω1, ω2, ω3), then Embedded Image is the matrixEmbedded Image(2.5)

Equations (2.1)–(2.4) with prescribed forces and torques F(t) and T(t), together with initial conditions b(tI)=b0, R(tI)=R0, v(tI)=v0 and ω(tI)=ω0, fully describe the motion of the body through the fluid for t>tI. In the following, we refer to the trajectory in terms of the state variables as a vector Embedded Image.

2.1.1 Velocity classifications

In the discussion of the motion of three-dimensional rigid objects, it is often convenient to refer to velocities in simple terms such as ‘heave’, ‘pitch’, ‘roll’, etc. Generic velocity vectors v and ω cannot, however, be classified as such—the instantaneous linear or rotational velocity of the body will be a complicated combination of several of these motions. In §3, however, we find it useful to classify the velocity vectors into such categories whenever possible, as described below.

Since the velocities v and ω are given in the body-fixed frame of reference shown in figure 1, if the vectors are close to one of the coordinate axes, they can be associated with what we term ‘simple’ translations or rotations. That is, translational velocity vectors close to one of the coordinate axes are associated with the simple translations ‘surge’ in the x-direction, ‘sway’ in the y-direction and ‘heave’ in the z-direction. Similarly, rotational velocity vectors close to one of the coordinate axes are approximately rotations about one of the primary body axes: ‘roll’ about the x-axis; ‘pitch’ about the y-axis; and ‘yaw’ about the z-axis. These simple rotations correspond to varying just one of the Euler angles.

By plotting the velocity vectors using spherical coordinates (R, Θ, Φ) instead of Cartesian coordinates, we construct a two-dimensional representation for each of the linear and angular velocity vectors. Each vector v can be written Embedded Image where −130°≤Θ<230° is the azimuthal angle and −90°≤Φ≤90° is the polar angle measured from the xy-plane. (Θ and Φ can be thought of as longitude and latitude, respectively, for a projection of the vectors onto a sphere.) The regions of the two-dimensional (Θ, Φ)-space corresponding to the simple translations and rotations are highlighted in figure 2. We use this representation of the velocity vectors throughout this paper.

2.2 Optimization procedure

For a given trajectory q(t), with t∈[tI,tF], describing the motion of the ellipsoid through the fluid, and associated forces and torques F(t) and T(t), we can associate a cost CEmbedded Image(2.6)for some function K, given below. The computation of the mechanically optimal trajectories is a constrained optimization problem: minimize the cost function C, subject to constraints on the trajectory q(t). The constraints consist of the equations of motion (equations (2.1) and (2.2)) as well as boundary conditions, and, if desired, further constraints on the applied forces and torques. The boundary conditions and further constraints are discussed in §2.2.1. Note that the duration (tFtI) of the trajectory is fixed.

The objective function we choose to be minimized can be thought of as a proxy for the metabolic cost of muscle activation for the fish to complete the trajectory. Our approach, which models the fluid as inviscid, does not lend itself to the use of an energy-based cost function. We choose an objective function similar to that used by Kanso & Marsden (2005),Embedded Image(2.7)The scaling factors αi are included so that the terms in the sum have the same dimension, and are equal to the reciprocals of the radii of gyration of the ellipsoid. For our model ellipsoids, αi range from 0.26 to 1.93 cm−1.

In order to implement the optimization, the trajectory q(t) has to be discretized. We achieve this by approximating q(t) by a piecewise linear function on an evenly spaced time grid Embedded Image. That is,Embedded Image(2.8)The equations of motion (which form constraints on the optimization) must also be discretized. The discretization of both the trajectory and the equations of motion is discussed further in appendix A.

The optimization is performed using the package Snopt (Gill et al. 2002), an implementation of sequential quadratic programming (SQP). The optimization routine finds a local minimum, subject to any number of constraints, of the objective function, given an ‘initial guess’ trajectory qinit(t). We discuss the sensitivity of the routine to choosing different qinit(t) in §2.2.2.

2.2.1 Constraints

The optimization of the cost function C is performed by varying the state variables q(t) and the applied forces and torques F(t) and T(t). Throughout the paper, the equations of motion are constraints on the optimization procedure, but, for each of the three hypotheses described in the introduction, we use a different set of additional constraints, as follows.

For the first hypothesis, that the fish trajectories are close to mechanically optimal, we set the start time of the trajectory to be the time of prey detection (defined in §2), so tI=tD, and the final time of the trajectory to be that of prey capture. (Recall that the duration of the trajectory is fixed to match that of the motion capture data.) Boundary conditions are chosen so the trajectory can be compared with the motion capture data (the MCF trajectory). Here, we fix the initial (that is, at t=tI) position, orientation and both linear and angular velocities of the trajectory to match that of the motion capture data. That is, q(tI) and Embedded Image are given. In this study, we also fix q(tF) (the final position and orientation of the ellipsoid) to that of the MCF trajectory, but Embedded Image is free to vary.

For the second hypothesis, concerning the redundant and essential motor capacities, we produce two additional sets of trajectories that we term ‘redundantly underactuated’ (RUA) and ‘essentially underactuated’ (EUA). The initial and final times of these trajectories and the boundary conditions on the position and velocity of the ellipsoid at these times are the same as for the first hypothesis, given earlier. However, we impose additional constraints on the allowed forces and torques for both the RUA and EUA trajectories. These are described in detail in §3.2.

For our third hypothesis, to determine whether a change in behaviour can be detected using optimal control results, we produce optimal trajectories for which the initial time tI is not equal to the time of detection tD. The final time of the trajectory remains equal to the time of prey capture. In this study, we additionally fix the velocity of the ellipsoid at t=tF to be that of the MCF trajectory. That is, q(tI), Embedded Image, q(tF) and Embedded Image are given.

2.2.2 Validation

The results of the optimization code described above give only approximate solutions to the Kirchhoff equations, owing to the discretization of both the equations and the trajectory. Ideally, we would run the optimization code with a very large N (the number of discretized time intervals), but this is not possible due to finite computing power. Here we describe how we choose an appropriate N.

For 10 randomly chosen trials (out of the 116 motion capture trials), we produced optimal trajectories for N∈{5,10,20,40,80} using the MCF trajectory as the initial guess qinit(t). The results from one of these trials can be seen in figure 3. Note that the CPU time required to produce the optimal trajectories increases exponentially with N since both the number of constraints and the number of variables increase with N.

Figure 3

Optimal trajectories as N is varied. Optimal trajectories obtained for one trial, for various N. Here N=5 is shown by the dotted curves, N=10 the dot-dashed curves, N=20 the dashed curves and N=40 the solid curves. As N is increased (N∈{5,10,20,40}), the optimal trajectories converge to the same solution. (a) Position of centre of mass. (b) Euler angles.

The optimal trajectory clearly converges as N is increased. We determine the N values for which the optimal solutions are sufficiently discretized by comparing the solutions with two different discretizations of the Kirchhoff equations. As N increases, we expect that the discrepancies between the solutions for different discretizations will decrease to zero.

We set a discrepancy level using the following two comparisons. First, using the functions F(t) and T(t) from the optimal trajectory, we obtained a corresponding q(t) by forward integrating equations (2.1)–(2.4) using a built-in Runge–Kutta (4,5) solver in Matlab (The Mathworks, Natick, MA, USA). This output was compared with the q(t) provided by the optimization code. Second, we implemented a different discretization of Kirchhoff's equations (from that used as constraints for the optimization): we solved equations (2.1) and (2.2) for F(t) and T(t) using (the discretized) v(t) and ω(t) from the optimal trajectory. This was compared with F(t) and T(t) provided by the optimization code.

For N ranging between 5 and 80, we found that, on average (for the 10 trials chosen above), the discrepancy measure we used decreased by an order of magnitude each time N was doubled. We set an allowed tolerance level for this discrepancy measure, and found that, for most trials, N=40 produced optimal solutions that passed this discrepancy check. Seven trials (out of 116) were excluded from the rest of our calculations because they did not meet our tolerance level on these checks.

For the same 10 trials used in our convergence tests described above, we also investigated the robustness of the optimization results to different initial guess trajectories qinit(t). This can indicate whether the code is finding a local or global optimum. In addition to the MCF trajectory, we used three other types of initial guess trajectories, as follows: a straight line between the boundary conditions for each of the state variables Embedded Image; a quadratic polynomial consistent with the boundary conditions; and a trajectory created by adding random noise (ranging up to 10%) to the MCF trajectory and then smoothing using a low-pass digital Butterworth filter. For each of the different qinit(t) used, we found that the code produced an optimal solution that was almost identical to the optimal solution produced using the MCF trajectory as the initial guess: the maximum error was 1.76×10−5 cm in position and 0.0017° in orientation. The optimization scheme (Snopt, see §2.5) only finds locally optimal solutions, but the results of this validation check provide convincing evidence that the solution obtained is the global minimum to the objective function. In the following results, we use the MCF trajectory as the initial guess qinit(t).

2.3 Body bend

In a number of the trials, the body of the fish does not remain rigid. As discussed in §1, the bending of the fish is mainly due to the small size of the behavioural tank and is not a propulsion mechanism. However, the rigidity assumption of our model is less valid in these trials; we measure the amount of bend in the body of the fish as follows, and exclude the most bent trials from our analysis.

The fish is capable of two types of body bends: lateral bend with respect to the midcoronal plane (assumed to go through the spine) and dorsoventral bend with respect to the midsagittal plane. The following calculation accounts for both types of bends. We identify six nodes on the surface of the ellipsoid. With its centre of mass at the origin, these six nodes are located at Embedded Image (recall that li are the radii of the ellipsoid). We find the corresponding locations on the unbent fish body when the ellipsoid and fish are aligned. Our bend measure is given by the distance between these six points on the fish and the ellipsoid, averaged over the number of nodes, ν=6 (note that this is similar to the error calculation given below in §2.4).

This bend measure was used to identify those trials that exhibit the most bend. Histograms of average bend (throughout the trial), maximum bend (throughout the trial) and initial bend were created. Each distribution had a well-defined ‘tail’ containing approximately 35 per cent of the trials. Any trial which fell into the tail of at least one of these histograms was classified as ‘bent’. Figure 4 shows the data for the initial bend, and the cut-off point for determining the bent trials. In total 40 trials (out of 109) were classified as bent. The MCF trajectories of the remaining 69 trials were used to obtain the below results.

Figure 4

Histogram of body bend at prey detection. Histogram of body bend values at the time of prey detection. All trials to the right of the vertical line are classified as ‘bent’. Histograms of maximum bend and average bend values over the length of the trials were similarly examined.

Figure 1b shows the body plan of the fish while bent. This amount of bend is approximately equal to the maximum amount of bend we allow in the ‘unbent’ trials.

2.4 Error calculation

In §4 that follows, we compare optimal trajectories, produced using the optimization scheme described earlier, with the MCF trajectory described in §2. The comparison between the two trajectories is calculated as given below.

We compute a ‘nodal error’ between two ellipsoids based on the nodes that make up the surface mesh of the ellipsoid. Each node is a point on the surface of the ellipsoid. The one-time nodal error at time tj is defined byEmbedded Image(2.9)where ν is the number of nodes (we use ν=256); Embedded Image are the coordinates of the ith node of the MCF ellipsoid at time tj; and Embedded Image are the coordinates of the ith node of the optimal ellipsoid at time tj. Figure 5 shows this graphically. For each trial, the nodal error, Embedded Image, is given by the average of E(tj) over the duration of the trajectoryEmbedded Image(2.10)where jI is the ‘initial frame’ (that is, Embedded Image) and jF is the final frame Embedded Image, and Embedded Image is the number of video frames of the motion capture data between tI and tF. Note that M is proportional to the duration of a trajectory because for all trials the video frames are of equal length, namely 16.7 ms. Across all trials M ranges from 22 to 66, with an average value of M=39. Using equation (2.10), we are able to compare errors between trials of different durations.

Figure 5

Nodal error. The nodal error compares locations of corresponding nodes on the surface meshes of two rigid ellipsoid models (shown here at one instance in time). Defined in equation (2.9), E(tj) is the average distance per node at time tj.

2.5 Computing

The computations of the optimal trajectories were performed on a 54 CPU (2 GHz G5, 1 GB RAM) cluster of Xserves (Apple Computer Inc., Cupertino, CA, USA) running OS X. An open source distributed computing engine (Grid Engine, Sun Microsystems, Santa Clara, CA, USA) was used to manage the computation across the nodes. Code was written in Matlab (The Mathworks, Natick, MA, USA) and compiled to portable executables for execution on the cluster.

For computation of the optimal trajectories, we used an implementation of SQP called Snopt (Gill et al. 2002). This routine finds a local optimum that may be the global optimum. We also implemented the built-in Matlab optimization routine fmincon on 12 randomly chosen trials. For these trials, we obtained results almost identical to those produced by Snopt. The maximum error between the solutions generated by Snopt and fmincon was 0.08 cm in position and 2.55° in orientation. The forces and torques produced by both optimization routines differed by less than 3 per cent. We found that Snopt generated the optimal solution at least twice as quickly as fmincon, therefore we used Snopt.

3. Results

We present results related to the three hypotheses outlined in the introduction. Firstly, we examine the extent to which the actual fish trajectories are mechanically optimal. Secondly, we investigate the extent to which yaw and sway are used in the optimal trajectories. Thirdly, we use optimal control results to detect a change in the fish's behaviour close to the time of prey detection.

The results for all the three hypotheses use optimal trajectories generated using the simplified model and optimization procedure described in §2.

3.1 Prey capture trajectories are mechanically optimal

For 69 prey capture trials, we compute mechanically optimal trajectories using the methods described in §2, to compare with the corresponding MCF trajectories. The initial time, tI, of the trajectories, at which the position and velocity of the ellipsoid are matched to the MCF ellipsoid, was equal to the prey detection time, tD (as determined by the method described in §2). The end time tF of the trajectory is the time of prey capture, and the position and orientation of the ellipsoid at this point are constrained to be equal to that of the MCF trajectory. In this section, we allow the velocity of the optimal trajectories at tF to be free to vary.

Figure 6 shows, for one typical trial, snapshots of the outline of the actual fish trajectory, and the mechanically optimal solution. This figure also shows the trajectory q(t) in terms of the state variables b(t), ϕ(t), θ(t) and ψ(t). For this example the trajectories are clearly very similar. Videos of four of the 69 trajectories, showing the fish trajectory alongside the optimal ellipsoid trajectory, can be viewed in the electronic supplementary material. The remaining trajectories can be viewed at http://www.neuromech.northwestern.edu/publications/Post08a/Post08aMovies.zip.

Figure 6

Trajectory snapshots. (af) Snapshots of optimal ellipsoid and actual fish trajectories (obtained from video-captured data) for one representative trial. Snapshots are viewed approximately side-on to the tank. Here time increases from top to bottom. The MCF trajectory is given by the solid curves and the optimal trajectory is given by the dotted curves. Note that the x, y and z coordinates are in the laboratory frame and hence the origin is at an arbitrary position. Both the ellipsoid model and the fish roll about the body x-axis as they translate through the fluid, a typical prey capture manoeuvre. A video of this trajectory can be seen in the electronic supplementary material video file S1.

The similarity between the MCF and mechanically optimal trajectories is also observed in the translational and rotational velocities (v(t), ω(t)). In figure 7, we use the velocity representation described in §2.1.1 to show the velocities for both the MCF and mechanically optimal trajectories. Points are plotted at 16.7 ms intervals throughout the trajectories, for all 69 trials. Most of the points in figure 7b,c lie in the backward surge region, which indicates that the optimal trajectories exhibit the same rapid-reversal manoeuvre that is typically observed in the fish's prey capture behaviour (MacIver et al. 2001). A concentration of points near the top of the forward surge region can also be seen. This corresponds to an upward heave during forward translation. It is important to note the lack of points in both the sway regions and the heave down region. This is consistent with our understanding that the fish does not possess the means to produce force or torque in these directions. A clustering of points is clearly visible in the roll regions in both figure 7e and figure 7f. There are very few points in the regions for yaw, a body rotation that the fish is not capable of performing. We examine the yaw and sway abilities in more detail in §3.2.

Figure 7

Velocity vectors for all trials. (a–c) Translational and (d–f) rotational velocities for the MCF and the mechanically optimal trajectories, for all 69 trials. Figure 2 shows how the figures are generated. The colour bar corresponds to the magnitudes of the velocity vectors (cm s−1 for v(t) and deg s−1 for ω(t)); velocities with the smallest magnitudes (25% of the data) have been omitted here. Both (b) and (c) show clustering of points in the surge and upward heave regions and a lack of points in the sway and downward heave regions. In (e) and (f), clustering is seen in the roll regions, but points are sparse in the yaw regions. (b,e) MCF and (c,f) optimal.

We calculate the nodal error, Embedded Image (as described in §2.4), between the MCF ellipsoid and mechanically optimal trajectories over all 69 trials. The mean nodal error is 1.31 cm with a standard deviation of 0.58 cm. We note that the mean nodal error for the 40 bent trials is not significantly different: 1.38 cm with a standard deviation of 0.50 cm.

We compare the mean nodal error for the mechanically optimal results to two length scales of the fish motion. These length scales are the length of the ellipsoids, and the total (integrated) distance travelled by the centre of mass of the fish during the trajectory (from the time of prey detection to capture). The average length of the model ellipsoids is 9 cm (which is only approx. 67% of the actual fish body length), and the average distance travelled by the fish is approximately 8 cm. The mean nodal error (1.31 cm) is much smaller than both of these. We note further that the range of distances travelled is large (3.5–16.8 cm), and there is no significant correlation between the distance travelled and the nodal error across the trials. Since the error is small compared with typical length scales of the motion, we assert that the trajectories are similar, and so the MCF trajectories are therefore close to being mechanically optimal. Since the MCF trajectories are a close approximation of the actual motion of the fish, the fish trajectories are also close to being mechanically optimal.

3.2 Distinguishing essential and redundant motor capacities using optimal trajectories

As previously noted in §3.1, the prey capture trajectory motions do not contain much sway or yaw motion. This can be seen in the MCF trajectories shown in figure 7. The dominant motions seen are surge and roll. The mechanically optimal results also show similar patterns.

As discussed in §2.2.1, the structure of our optimization code makes it easy to add additional restrictions to the form of the applied forces and torques F(t) and T(t). In order to evaluate to what extent the optimal control results discussed previously make use of yaw torque and sway force, we produce optimal trajectories in which the applied force F and applied torque T are restricted. That is, we remove the ability of the ellipsoid to yaw and sway, so Embedded Image and Embedded Image. We refer to these trajectories as the redundantly underactuated (RUA) trajectories. As a comparison, we also compute optimal trajectories that have the surge and roll capabilities removed. We refer to these trajectories as the essentially underactuated (EUA) trajectories. In all the optimal trajectories, the final velocities are again allowed to be free to vary.

We compute RUA and EUA trajectories for 26 of the trials, and compute the nodal errors for each trajectory. We also compute the difference in the cost for these trajectories, relative to the fully actuated trajectories. The results are summarized in table 1. It can be seen in the table that the nodal error between the mechanically optimal and MCF trajectories increases much less when yaw and slip are removed than when the surge and roll are removed. Figure 8 shows the trajectories for the MCF ellipsoid and the fully actuated and underactuated optimal ellipsoids, for one representative trial.

View this table:
Table 1

The mean (s.d. in parentheses) nodal errors and mean costs for the comparison of optimal control trajectories to the MCF ellipsoid. (Here we consider only the 26 trials for which we have underactuated results that pass validation. All nodal errors are given in cm and all costs are given in units of force squared (i.e. (g cm s−2)2). The percentage increase in cost is compared with the mean cost of the fully actuated results.)

Figure 8

Underactuated optimal trajectories. Trajectories for one representative trial for the MCF ellipsoid (solid black), the fully actuated optimal ellipsoid (dashed blue), the RUA ellipsoid (dotted red) and the EUA ellipsoid (dot-dashed green). The RUA trajectory is clearly similar to the fully actuated trajectory, while the EUA trajectory is significantly different. (a) Position of centre of mass. (b) Euler angles.

Ideally we would produce underactuated optimal trajectories for all 69 trials (not just 26). However, almost all of the 26 trials required N≥80 in order to produce optimal trajectories (both for EUA and RUA) that satisfy our tolerance level for the validation described in §2.2.2. As explained in §2.2.2, when N is increased, the CPU time necessary to produce the optimal trajectories increases significantly. In general, computing the EUA trajectories required more computing power than the RUA trajectories. We believe that this is due to our choice of initial guess trajectory, qinit(t). In all our results, we use the MCF trajectory as the initial guess for the optimization. As shown in figure 7, most of the MCF trajectories contain a large amount of surge and roll motion. Since the EUA optimal trajectory by definition does not include any surge or roll, it will therefore be quite different from the MCF trajectory. This probably makes the MCF trajectory a bad initial guess for finding the optimum. In fact, the 26 trials that produced valid EUA trajectories exhibited significantly less surge and roll than the other 43 trials.

3.3 Using optimization to detect behavioural mode switches

In §§3.1 and 3.2, all optimal trajectories started at the time of prey detection (t=tD, as determined in MacIver et al. (2001)). In this section we investigate how the optimal trajectories vary when the initial time is some other point in the fish trajectory.

The video capture data include a portion of fish trajectory before the prey is detected. During this time, the fish is searching for prey, typically moving forward at approximately 10 cm s−1 while the body is pitched with the head down (MacIver et al. 2001). For each of the 69 trials, we produce mechanically optimal trajectories with initial times varying over the entire time of the motion capture trajectory. That is, suppose the MCF trajectory has duration tF, with m video frames captured at times Embedded Image (where ts=16.7 ms is the inter-video frame spacing). For each trial, we produce (m−1) optimal trajectories for each possible initial time tI chosen from Embedded Image. Owing to the long CPU time required, for longer trials we generate optimal trajectories only for tI equal to the 10 frames before and after detection, and then every fifth frame away from this point. This captures the important changes near tD, but ignores unnecessary information far away from tD.

The initial position and velocity of the optimal trajectory is matched to the position and velocity of the MCF ellipsoid at t=tI. The endpoint of the trajectory is fixed to be the end time of the motion capture trajectory (t=tF, as before). The final position and velocity of the optimal trajectory are set equal to that of the MCF trajectory at t=tF.

For each trial, we shift time so that the detection point tD is at t=0, i.e. tttD. We then compute the time-averaged nodal error, Embedded Image, between the mechanically optimal ellipsoid and the MCF ellipsoid, for each trajectory (that is, for each initial time tI). We plot Embedded Image against tI to examine how the error depends on the start time of the optimal trajectory.

Across all the trials, we see two distinct types of features in these plots. We classify these by eye as ‘bump’ and ‘no-bump’, and show a representative trial for each case in figure 9. A trial is classified as no-bump if its nodal error curve is approximately monotonically decreasing near (±10 video frames) tD. On the other hand, a bump trial has a local minimum and maximum in its nodal error curve near tD.

Figure 9

Nodal error versus initial time frame. The nodal error between the optimal and MCF trajectories for one representative trial (a) without and (b) with the post-detection bump.

We give a possible reason for the no-bump/bump effect in §4. Our optimal trajectories are a better model in the no-bump cases, and we now describe how we use the optimal trajectories to detect a change in the behaviour of the fish before and after prey detection.

In the no-bump trials, we see a sharp change in the gradient of the nodal error as a function of optimal trajectory start time, close to the detection point tD. This can clearly be seen in figure 9a. Figure 10a shows data averaged over all the no-bump trials, and the same change in gradient can still be seen. We attribute the sharp change in the gradient to a change in the fish's behaviour.

Figure 10

Mean nodal error versus initial time frame. Mean nodal error (shown with standard deviation bars) over all 69 trials plotted versus the initial time frame relative to the detection time frame, for (a) no-bump and (b) bump trials.

Recall that we match the initial position, orientation and velocity of the optimal trajectory to the MCF trajectory. Consider an optimal trajectory with initial time tI pre-detection. We do not expect the optimal and MCF trajectories to be similar, since before prey detection the fish has no reason to move in an optimal fashion towards the prey. However, since the optimal and MCF trajectories have initial and final positions and velocities which are equal, it is expected that the nodal error would decrease as the trajectory duration decreases (the nodal error is clearly zero when the motion-capture trajectory contains fewer than five frames as all the coordinates are completely constrained by the initial and final conditions).

However, post-detection, when we expect the fish to follow an optimal trajectory from its current position to where the (relatively stationary) prey is sensed to be located, we would not expect such a rapid decrease in nodal error as trajectory duration decreases. This is because any post-detection portion of the MCF trajectory is already close to mechanically optimal, according to the results already presented.

4. Discussion

There have been a number of prior studies of animal locomotion using optimization techniques. Two distinct methods are found in the literature. The first is parameter optimization, which is different from the method we have used. This can be seen in Tam & Hosoi (2007) (who study a three-link swimmer in low-Reynolds number flow) and Kern & Koumoutsakos (2006) (studying anguilliform swimming). In these studies, there are a number of parameters that can be varied, which may affect body shape, stroke pattern, etc. For a specific set of parameters, the resulting trajectory is computed from the equations of motion, and this results in an output efficiency. The parameters can then be varied to optimize the efficiency. A study by Lauga & Hosoi (2006) on the locomotion of gastropods uses similar ideas; the parameters over which the optimization is performed affect the rheology of the mucus produced by the gastropod. The second method, which is the one we have used, is based on optimal control. It is similar to that used by Srinivasan & Ruina (2006), who study a simplified model of a two-legged walker, and Kanso & Marsden (2005), who study a two-dimensional three-linked swimmer, and our own prior work on the knifefish (MacIver et al. 2004).

Our current approach uses the Kirchhoff equations for a solid rigid ellipsoid moving through an inviscid fluid. This approach is suitable for transient movements in which acceleration dominates drag. It will be particularly useful for application to underwater vehicle designs or to fish which keep their body relatively rigid and rely on fin movements decoupled from body movements for propulsion; however, it may also be useful to apply the approach for species of fish where these conditions do not hold. To do so, ideally one would first measure the body position of the fish from initiation of the transient movement to its termination. Then, a body-fitted ellipsoid would be fitted to the measured trajectories, corresponding to our MCF trajectories. The boundary conditions would be extracted from the initial and final positions and velocity of the fish. If such data were not available, one could instead use an approach whereby knowledge of typical fish behaviour could be used. For example, typical search velocities and orientations could form the initial conditions. In the case of a behaviour such as prey capture, any points in space which are in the prey reactive volume or sensory volume (Asaeda et al. 2002; Snyder et al. 2007) of the fish could be used as the final position condition. Typical capture velocity would be used as the final velocity condition. Orientations for the end of the transient motion would also need to be approximated. Using these boundary conditions, an optimal solution can be obtained using the equations of motion above and an optimization algorithm such as SQP. Finally, the optimal solution would be compared with the measured trajectory through the use of the error metrics presented earlier.

While this approach cannot be expected to elucidate the hydrodynamics of propulsion, it can be expected to provide insight into why the fish is able to move more effectively in some directions rather than others. In general, as discussed further below, it should predict motion capabilities consistent with the least-effort directions of movement. Failure of this prediction can be informative of other constraints not included in this highly idealized approach.

Sections 4.1–4.3 focus on each of our three hypotheses individually. We then discuss how the approach we have developed here may contribute to our understanding of the complementary nature of fish body plan and actuation capabilities. We finish with a discussion of related future work that we intend to pursue.

4.1 The optimality of prey capture trajectories

In §3.1 we generated optimal trajectories that minimize the cost function C, under the assumption that the fish is a rigid ellipsoid moving in an inviscid fluid. We show that these trajectories are similar to actual fish prey capture trajectories. It is clear that our simplified model of the fish motion suffers from a number of drawbacks. We begin this section by a discussion of these.

The most obvious simplification we make is the modelling of the fish body as (i) rigid and (ii) ellipsoidal. Concerning the rigidity assumption, the unique way in which this fish propels itself—by undulating a ventral ribbon fin while keeping the trunk semi-rigid—makes this a reasonable approach for this species. The tendency of the fish to swim with a rigid trunk has been noted in the literature for some time (Lissmann 1958, 1961; Blake 1983). One suggested reason is that the trunk carries the electrical emission system of the fish; thus, when the tail bends, there is a large modulation of the sensory receptors on the trunk surface (Bastian 1995; Chen et al. 2005). Decoupling propulsion from sensory acquisition enables the independent control of the trunk, which is populated with sensory receptors, as a sensory organ, as seen in scanning behaviours where the fish will arch their trunk around a novel object while moving forward and backward (Assad et al. 1999). During prey capture strikes, the body bends that are present prior to prey detection are probably due to the use of a relatively small tank (three body lengths in the longest direction) for the prior behavioural study, necessitating frequent turns: during the strike itself, however, any body bend that was present was observed to rapidly decrease (MacIver et al. 2001).

Concerning the approximation of the body plan with an ellipsoid, the selection of appropriate dimensions for the ellipsoid is not unique. We chose to match the fish's width, depth and volume, while sacrificing matching the length. This choice affects the values for the added masses used in the equations of motion, which in turn affects the optimal trajectories. We suspect that choosing an ellipsoid of slightly different dimensions would not significantly alter the trajectories. However, the added masses themselves may not be equal to those for the actual fish body, which would be very difficult to calculate.

The simplified equations of motion do not include the effect of drag, since the fluid is assumed to be inviscid. However, since the drag forces would act in the same direction as those created by the addition of the added masses in the equations of motion, we believe that this would not significantly affect the trajectories. In addition, in a recent study (Shirgaonkar et al. in press) it has been shown that, during prey strikes, the acceleration reaction forces are significantly larger (15 mN) than those due to drag (1 mN). Hence it is reasonable to neglect the effects of drag in this model. However, one effect which might be significant is that without drag, the cost to our ellipsoid for accelerating is the same as that for decelerating. Since the fish may rely on passive drag during deceleration, while having to activate muscles for acceleration, this may not be a good approximation.

Our optimization approach allows us to vary the constraints on the ellipsoid trajectory. This invites the question—which parts of the trajectory should be fixed and which should be allowed to vary? We have fixed the initial position and velocity of the ellipsoid to that of the MCF trajectory, as this is a natural way to compare the trajectories. However, the correct constraints to use at the end of the trajectory are less obvious. We have considered two cases where the final velocity of the ellipsoid trajectory was either fixed to the final velocity of the fish trajectory or allowed to be free to vary. In both cases, the position and orientation of the body were fixed. It might be possible (and even desirable) to fix only the position of the head of the body (i.e. so it can capture the prey) and allow the orientation to vary. In fact, we intend to continue this line of research in the future, although initial investigations have indicated that this may be computationally much more difficult than the currently presented approach.

We also note that the optimal ellipsoid ‘knows’ the endpoint of the trajectory throughout, but the fish will be continually collecting sensory data and possibly adjusting its trajectory. This would be particularly noticeable in trials in which the prey moves a significant amount after its detection. In these cases, an optimal control approach that incorporates feedback (e.g. Todorov & Jordan 2002) is likely to give better results. In the majority of trials considered here, the prey move a negligible amount (less than 1 cm).

Across the trials, we see that a large proportion of the nodal error is caused by a lengthwise translational shift between the optimal ellipsoid and the MCF trajectory (see videos in the electronic supplementary material). This can be attributed to an observed ‘lag’ between the two trajectories. Examining individual trajectories in detail led us to believe that the lag is caused by ‘loitering’ of the fish in MCF trajectories near the end of the trial. That is, the fish will slow down before capturing the prey to an extent not seen in the optimal trajectories. We suspect this is an additional constraint on the fish motion not captured by our modelling approach—the fish cannot capture the prey while travelling at too high a speed. (The actual engulfment of prey occurs through negative buccal pressure induced by hyoid depression. The negative pressure causes a small volume of water to be drawn into the mouth along with the prey. Movement relative to the prey may make this a more error-prone process.)

4.2 Distinguishing essential and redundant motor capacities using optimal trajectories

In §3.2 we use optimal trajectories to distinguish between those forces and torques which are essential to fish motion and those which can be thought of as redundant. We show that yaw torque and sway force are redundant, but that surge force and roll torque are essential to the fish's prey capture motions. We can see this from table 1; the EUA trajectories have significantly larger nodal errors than either the fully actuated or the RUA trajectories.

We also examine the cost of both sets of trajectories in table 1. It is clear that the cost of the EUA trajectories is significantly higher than either the fully actuated or the RUA trajectories. We note that although RUA cost is twice that of the fully actuated trajectories, the trajectories can still be quite similar, as shown by the similarity of the nodal errors; this is because the ‘cost landscape’ of the trajectory space may be very complicated. The nature of the optimization problem (that we are minimizing the cost while satisfying constraints) means that, as constraints are added, the cost will always increase.

Computational analyses of the hydrodynamics of the ribbon fin (Shirgaonkar et al. in press), actuated with a travelling sinusoid, allows us to estimate forces. We find that this fin, the fish's primary propulsor, can generate positive and negative surges, as well as positive heave. Roll appears to be accomplished by holding the pectoral fins into the flow as steering surfaces, producing an upward jet on one side of the body and a downward jet on the other side of the body to generate a moment around the body axis (M. A. MacIver 2001, unpublished observations). Pitch also appears to be accomplished in this fashion, but the left and right pectoral fins are held at similar angles of attack so that movement in a direction parallel to the body axis creates a force pushing the head up or down, causing a pitch moment around the centre of mass of the fish which lies just posterior to the fins. In summary, the fish can generate 1.5 of 3 linear forces, surge and positive heave, and 2 of 3 rotational forces, roll and pitch.

Given these observations, we would therefore expect that, since the mechanically optimal trajectories are close to actual fish trajectories, they would not require sway or yaw thrust, since the fish cannot create these forces. This is consistent with the results we find.

We also point out that as in §4.1 we are constraining both the position and orientation of the ellipsoid at the end of the trajectory to match that of the fish. Clearly, since the fish is itself underactuated, it will only reach positions consistent with its motor capacities, which should also be positions that a similarly underactuated ellipsoid may reach. However, there may be less costly trajectories the fish could take if we did not constrain the terminal position and orientation of the ellipsoid to be the same as the fish, and these less costly trajectories may require more degrees of freedom than the fish has. In the future, we will be investigating the extent to which our minimal actuation results hold when the only constraint on the final position of the ellipsoid is the position of the head, with the final orientation left free.

4.3 Using optimization to detect behavioural mode switches

In §3.3 we investigate how the nodal error between the MCF and optimal trajectories changes as the initial time of the trajectory is varied. We find two types of behaviour which we term bump and no-bump.

We explain the reason for the bump and no-bump behaviours as described below. A typical behaviour in the prey capture trajectories is the ‘rapid-reversal’ manoeuvre; the fish moves forward until shortly after prey detection, when it changes direction and moves backwards until the prey is captured. We note a distinct difference in this motion between the MCF trajectories that produce ‘bumps’ and those which produce ‘no-bumps’, as described in §3.3. In the no-bump trials, the deceleration of the fish is approximately constant at the end of the trajectory, until the prey is captured. However, in the bump trials, both the magnitude of the deceleration and the velocity of the fish decrease towards zero for approximately 150 ms near the end of the trial, before prey capture. One possible explanation for this behaviour is that, in these trials, the fish did not accurately detect the location of the prey, or, as mentioned above, that the fish cannot capture the prey while travelling at high speeds. We find that this loitering effect is correlated with those trials that exhibit a bump as described in §3.3.

In the no-bump trials, we find a change in the slope of the graph of the nodal error Embedded Image plotted against the initial start time tI, which we relate to a change in the behaviour of the fish. However, this change occurs in an average of 50 ms after prey detection. There could be a number of reasons for this delay. We give one possible hypothesis here.

The detection point was determined in a previous study (MacIver et al. 2001) to be the time at which the fish begins decelerating due to the detection of the prey. The optimal ellipsoid, started from the time of behavioural reaction, has no uncertainty about where its trajectory must end. By contrast, the fish has considerable uncertainty about where its trajectory must end. For example, from prior work, we know that, at the time of prey detection, the prey-related signal is extremely small. The length scale of the signal change on the body at the time of detection is also roughly 10 times the diameter of the prey (Nelson & MacIver 1999; Snyder et al. 2007). Given these factors, it is highly unlikely that the fish has a clear estimate of the prey position at strike initiation. It is therefore possible that localization uncertainty causes the fish to initially pursue a trajectory that is not consistent with the actual location of the prey (and therefore not consistent with the optimal trajectory). A relevant time scale that is similar to the 50 ms delay is the total delay from sensory input to motor output, which has been estimated to be approximately 120 ms in this animal (Snyder et al. 2007). This hypothesis could be examined in more detail by varying the degree of sensory uncertainty, such as by making the water more conductive, which past work has shown decreases the distance at which prey are sensed (MacIver et al. 2001).

We have provided evidence that the identification of the utility function that is being minimized by the motor system during a behaviour may be useful for ascertaining when an animal switches into that behaviour. Such an approach may find application in neuroethology and other organismal-level research where whole-animal behaviour needs to be quantified and segmented. This may work particularly well when the terminal point of the behaviour is easily ascertained, since this then provides a natural starting point to work back from.

4.4 The complementarity of body plan and actuation

One of the intriguing results we found is that, even if we endow the fish model with all six rigid degrees of freedom, enabling it to move in ways that the real fish cannot, the best trajectories are the ones that use only the limited motion capabilities that the real fish possesses. A natural question is how this could arise. One possibility has already been alluded to above: the boundary conditions are those consistent with movements of an underactuated fish, so in this way we can be excluding trajectories that would require less effort for a fish with different actuation capabilities.

However, there is another way in which to think about the problem. The favourable modes of movement for our ellipsoidal model are determined by the anisotropies in the mass and moment of inertia matrices given in §2.1. For example, the mass matrix indicates that it is better for the animal to move forward and backward than to move laterally, since the effective surge mass for moving forward or backward is one-third the effective sway mass for moving laterally. In the moment of inertia matrix, the differences are even larger: to yaw the body involves over 30 times more moment than to roll the body. Indeed, turns observed prior to the initiation of a prey strike are executed through lateral body bends, not yaw rotations, as discussed above. One simple prediction, therefore, is that, to the extent that an animal's body plan is anisotropic in its mass or moment of inertia, for high acceleration transient movements at the animal's biomechanical extreme, the body will be configured so as to minimize those associated masses and inertias during the movement. In doing so, the animal can maximize the movement resulting from its limited force and torque generation capability. Thus, force generation capabilities should be so as to translate or rotate the body consistent with such minimization. It seems likely that, over the course of evolution, the shape of the body and its actuation capabilities coevolve to ensure this.

An additional level of integration to be considered is the animal's sensory capacity. Since the boundary conditions of our optimization are set by objects necessarily within the sensory volume of the fish, sensory capabilities are implicitly involved in the trajectories analysed here. In the case of the weakly electric fish, prey are detected in all directions from the body. For most visually guided fish, prey are detected only in a wedge-shaped sector of space ahead of the body (e.g. Snyder et al. 2007, figure 8). The body shape, and thus actuation abilities, must be so as to enable rapidly reaching all parts of the prey detection sensory volume.

For example, the weakly electric fish's sensor density is the highest on the dorsum, where it tends to detect prey (Snyder et al. 2007). The positive heave observed in figure 7, which is generated as a proportion of surge by the travelling waves on the ribbon fin (Shirgaonkar et al. in press), is the motor complement to this aspect of the fish's perceptual abilities. Prey that are detected behind the head and lateral to the body wall are caught by combining a rapid reversal with a roll to place the prey on the midsagittal plane, while positively heaving (MacIver et al. 2001, figure 10a). All these motions are highly favoured by the body shape.

In comparison, with the forward-directed and flattened wedge-shaped sensory volume of a visually guided fish, strikes will often be propelled by one or more rapid tail beats. This type of swimming is unstable in yaw (Weihs 2002), and consequently the fish is able to quickly manoeuvre to lateral positions of the sensory volume. The fish does not need to visit positions initially behind the head and lateral to the body wall as prey would not usually be detected there. Moreover, the very mechanism that the knifefish uses to visit these positions—roll—is largely prevented by a combination of a deep body plan and a dorsal fin (Weihs 2002).

These considerations suggest that, even in cases where fish movement is poorly modelled by the rigid-body approximation, it may still be worthwhile to apply the simple mathematical approach presented here for insight into how body shape, motor capabilities and sensory performance are integrated.

4.5 Future studies

We consider now two extensions of our current model, which we intend to pursue in future studies. Our current optimization routine considers the duration of the trajectory to be fixed. This is clearly not a realistic constraint on the fish motion—in addition to moving in a mechanically optimal manner, it is probable that the fish will try to minimize the time taken to capture the prey. It would be possible to consider a joint time/mechanical optimization in the following manner. The number of discretization points, N, would be fixed, as before, but the length of time between each point would be allowed to vary. The cost function would then be a function of both the total length (tF−tI) of the trajectory and the mechanical cost function used here.

In §3.2 we considered constraining the forces available for the model ellipsoid to be more similar to those actually found in the fish body by not allowing the fish to sway or yaw—forces and torques which the fish is not capable of generating. This procedure could be extended further to better model the forces produced by the fish as it moves. As described in §4.2, the fish can generate 1.5 of 3 linear forces and 2 of 3 rotational forces. In addition, these forces are unlikely to be independent of each other, but will be restricted by some relationships. Informed by hydrodynamic studies (Shirgaonkar et al. in press), we can include these constraints on the forces moving the model ellipsoid so that it better models the fish.

Acknowledgements

This work was supported by NSF grant IOB-0517683 to M.A.M., and in part by NSF grant DMS-0709232 to M.S. T.M.P. acknowledges support from the DoD through the NDSEG Fellowship Program. The authors would like to thank Eva Kanso for helpful discussions regarding the discrete Lagrangian approach, and James Snyder for computing assistance. We also wish to thank the anonymous referees for their helpful comments.

Appendix A. Discrete lagrangian derivation

This appendix describes the discretization procedure of the Kirchhoff equations that we use in our optimization routine. We use a discrete Lagrangian derivation of the equations of motion, following methods used in Kanso & Marsden (2005).

Recall that we use the Euler angles ϕ, θ and ψ to describe the orientation of the body frame with respect to the laboratory frame. Combining the three rotations produces a matrix R, which maps vectors in the body frame to vectors in the laboratory frame. We write the rotation matrix R in terms of the Euler angles asEmbedded Image(A1)and then the angular velocity ω is given byEmbedded Image(A2)

Since the body is neutrally buoyant, the Lagrangian is equal to the total kinetic energyEmbedded Image(A3)whereEmbedded Image(A4)Recall that Embedded Image, where bT=(x, y, z) is the position of the centre of mass of the body in the laboratory coordinate frame.

A.1 Discretization set-up

We discretize the trajectory q(t) by approximating it by a piecewise linear function on an evenly spaced time grid Embedded Image. That is,Embedded Image(A5)whereEmbedded Image(A6)The applied forces and torques are similarly discretized, and we writeEmbedded Image(A7)

We introduce ‘hatted’ variables that are evaluated at the midpoints of the grid, and write derivatives asEmbedded Image(A8)and for each component of qn, we similarly writeEmbedded Image(A9)We write Embedded Image to mean R evaluated at Embedded Image, and Embedded Image to be the vector Embedded Image. Note thatEmbedded Image(A10)

The Lagrangian is discretized by writingEmbedded Image(A11)whereEmbedded Image(A12)

The discrete Lagrange–d'Alembert principle (see Junge et al. 2005; Kanso & Marsden 2005 for details) gives the following equations:Embedded Image(A13)Embedded Image(A14)Embedded Image(A15)Here, p0 and p1 are the initial and final generalized momenta of the body and Embedded Image and Embedded Image are the left and right discrete forces,Embedded Image(A16)

Footnotes

    • Received July 4, 2008.
    • Accepted September 4, 2008.

References

View Abstract