Geometric Methods for Efficient Planar Swimming of Copepod Nauplii

Copepod nauplii are larval crustaceans with important ecological functions. Due to their small size, they experience an environment of low Reynolds number within their aquatic habitat. Here we provide a mathematical model of a swimming copepod nauplius with two legs moving in a plane. This model allows for both rotation and two-dimensional displacement by the periodic deformation of the swimmer’s body. The system is studied from the framework of optimal control theory, with a simple cost function designed to approximate the mechanical energy expended by the copepod. We find that this model is sufficiently realistic to recreate behavior similar to those of observed copepod nauplii, yet much of the mathematical analysis is tractable. In particular, we show that the system is controllable, but there exist singular configurations where the degree of non-holonomy is non-generic. We also partially characterize the abnormal extremals and provide explicit examples of families of abnormal curves. Finally, we numerically simulate normal extremals and observe some interesting and surprising phenomena.


Introduction
Microcrustaceans known as copepods are one of the most abundant animals on Earth. They are a type of zooplankton that serve as an important link in the marine food web. As they are prey for many larger aquatic creatures, they must adapt strategies that help them maximize their survivability. Observations have shown that different copepods have adapted different types of movement to efficiently forage for food and evade predators [9,10,20,28]. Figure 1 shows a nauplius of the copepod species Bestiolina similis.
In the world of microorganisms, water becomes a very viscous fluid in which movement produces negligible inertia. This is known as a low Reynolds number environment, as the Reynolds number R represents the ratio between the inertial force due to momentum and the viscous force experienced from the resistance of the liquid. An object swimming in some fluid experiences the Reynolds number R = avρ η , where a is the characteristic dimension of the object, v is the velocity of the object, ρ is the density of the fluid, and η is the fluid viscosity. In water, organisms as small as bacteria have a Reynolds number of approximately 10 −6 to 10 −4 , whereas humans have a Reynolds number of approximately 10 6 . Nauplii of the paracalanid copepod Bestiolina similis, as shown in Figure 1, have lengths 70-200 µm and swim at Reynolds numbers of 10 −1 − 10 1 (see [16] and the references therein). To put swimming at low Reynolds number in perspective, humans have a Reynolds number around 10 2 when swimming in molasses.
Edward Mills Purcell's talk from 1976, Life at Low Reynolds Number ( [22]), first popularized the concept of swimming in an environment with low Reynolds number and was foundational in the study of microswimmers. He stated the scallop theorem, which says that complete reciprocal motion cannot produce any displacement when swimming at low Reynolds number. Here a reciprocal motion involves two sequences of motion where the second motion, called the recovery stroke, is the reverse of the first motion, which is called the power stroke. For humans swimming in water, the time it takes to do a stroke plays a role in the induced displacement due to the inertia terms in the Navier-Stokes equations. So humans are able to swim forward by completing the second recovery stroke faster than the first. But for microswimmers, this inertia is negligible, so any displacement induced by the first stroke is reversed by the second recovery stroke. As a consequence, microorganisms must move in other ways, like utilizing a flagella or moving pairs of legs in an asynchronous manner. For example, for copepods that are bilaterally symmetrical, the symmetric pairs of legs move simultaneously, but adjacent legs move out of sync in order this range, we have used one of mmers available, the nauplii of lina similis (length 70-200 mm) m at Re of 0.1-10 [18], which is nd intermediate Re. Simplificaon predictions can allow direct gical and kinematic parameters e are free. A relatively simple n applied that can be confined hout sacrificing predictive capaine how well such a simplified for observed swimming behavtransition zone, deviations are new insights into swimming at nd inertial forces are important. employed is based on slenderpted from one that was recently differs from previous models in ertia for propulsion. By accountcally measured dimensions and appendages, our model was f the body over time and comct observations to assess the addition, the vetted model was n to displacement of each appennd appendage stroke phase in role in naupliar propulsion. s y of naupliar swimming f angular position of individual ment were made for nauplii of from cultures maintained in the r under standard conditions as nz [20]. Briefly, B. similis adults kton collections from Kaneohe cultured at ambient temperature ime, and fed ad libitum with live ). Experimental nauplii were isofied to stage using morphological idth measurements [17]. re placed into small Petri dishes od levels. Experimental nauplii corresponding to developmenus fast swims were recorded at eo system (Olympus Industrial rted microscope (Olympus IX70) he video files were converted into at) and analysed using IMAGEJ web.nih.gov/ij/). Six swim epidage angles and location over cles at 0.2 ms intervals. The angle using the main axis of the nauplius nning electron micrograph of an cation was determined by tracking erior medial margin of the head in swim sequence. Five additional r location during rapid swims to d net displacements. Swims were 1b), which was characterized by pendage: first antenna (A1) point-

Model formulation
To determine the extent to which observed locomotion of a nauplius could be accounted for based on observed appendage movements and the assumptions of a low Re regime (see Introduction), we employed a model of swimming with rigid appendages adapted from one based on slender-body theory for Stokes flow [19]. The aim of the model is to predict the position of the body, as the angle of each leg changes over time. The model provides us a reasonable approximation for long and slender appendages paddling at low Re [21], which omits inertia, as explained in the Introduction. It makes several additional simplifying assumptions intrinsic to its formulation. The copepod nauplius has a compact rounded body (figure 1) that is simplified in the model as a sphere with a diameter that is the mean of the length and width of its body. Using the more accurate prolate ellipsoid shape instead made little difference in predicted displacements. Naupliar appendages are relatively rigid elongate rods, slightly tapering at both ends, again with rounded cross section. In the model, they were simplified and represented as uniform cylinders, to produce a net displacement. These observations are important as they could potentially be used to design microscopic robots that move in similar ways. One application involves using bacteria based nanoswimmers to transport drugs from a loading point to a destination such as cancer cells ( [13]).
The main microorganisms of focus in this paper are copepods. Most types of copepod are only able to move in a motion called swimming-by-jumping ( [2]). This motion is similar to the one described in the previous paragraph: it involves moving symmetric pairs of legs in a way such that adjacent pairs move asynchronously. In other words, the power stroke and recovery stroke alternate among each of the pairs of legs. One way to model this is by moving each symmetric pair of legs in a reciprocal motion while introducing a phase lag between each pair of legs ( [26]). As in such models, here we restrict the copepod's motion to a plane for simplicity, despite the fact that the actual animals live and move in three dimensions.
Here we model the copepod as a slender body in Stokes flow as in [26,4,7,8,11,6]. Other models of micro-swimmers capable of rotation appear in [12,23,15]. It is important to note that real living copepods do indeed perform rotations to both evade predators ( [24]) and to capture prey ( [10]). In [19], the authors analyze such rotational maneuvers (yaw, pitch and roll) via high-speed video observations of copepod larvae.
One dimensional translational motion for the copepod model has been well studied. It is possible to achieve positive displacement along an axis using as few as two pairs of legs moving in a reciprocal motion ( [26]). Methods from sub-Riemannian geometry and Hamiltonian dynamics have been used to find efficient optimal strokes in the translational case; efficiency was defined as the ratio between the displacement resulting from a stroke and the length of the stroke. Numerical methods were used to determine the optimal strokes maximizing this efficiency ( [4,7,8]). Here the term stroke refers to a periodic motion of the legs.
Here we generalize this prior work by analyzing planar motions. To produce orientation changes we need to break the symmetry of the pair of legs. A first attempt was made in [11] by looking at three independent legs oscillating sinusoidally; here we generalize that approach to include all strokes but using two legs. We find that rotation by strokes is indeed possible with only two legs. We also show that the two-legged system is controllable, although the difficulty in steering locally depends on the initial state; in other words, the system possesses singularities, which we classify. Taking the mechanical energy expended as our cost function, we develop the two-legged copepod movement as an optimal control problem and apply the Pontryagin maximum principle ( [21]) to study both the normal and abnormal extremals. We partially characterize the abnormal extremals, and provide some explicit examples. Finally, Figure 2. The rotating copepod with 6 legs. Note that angles θ i are associated to the body frame while the angles α i are associated to the inertial frame .
we utilize the optimal control software Bocop to simulate normal extremals ( [27]). Among our simulations we find copepod motions which produce net rotation without net displacement, we characterize the optimal motions which produce rotation with no conditions on displacement, and we discover paths in the xy-plane which appear to be Euler elastica ( [3]).

Methods
We consider a simplified copepod microswimmer in a low Reynolds number environment. The idealized copepod consists of stiff slender legs and a body of negligible radius in comparison to the length of its legs. In this section we will develop our mathematical model, derive the equations of motion, describe the copepod's motion as an optimal control problem, and develop the appropriate version of the maximum principle.
2.1. Model. We assume the copepod moves in a plane and possesses 6 independently moving legs, three on each side of the body. The position of the copepod at time t can be described by the vector (x(t), y(t), φ(t)) T , where x and y represent the usual Cartesian coordinates on the plane and φ represents the orientation of the copepod with respect to the positive x-axis. Let θ i denote the angle between the copepod's orientation and the i th leg, and let α i = θ i + φ denote the angle between the i th leg and the positive x-axis. See Figure 2 for an illustration.
We denote the state of the copepod at time t by the vector while the position and orientation coordinates alone will be written asq = (x(t), y(t), φ(t)) T . Thus our configuration space is ostensibly R 9 , however in order to prevent the legs from passing each other, we impose the constraint For the rest of this paper we will focus on a simplified copepod with two independent legs, one on each side of the body. See Figure 3. This simplification allows us to conduct a mathematical analysis and is justified by assuming the three legs on each side of the body are collapsed into one stronger leg. As will be seen in Section 4, even with this simplification we obtain swimming motions reflecting actual laboratory observations.
In most of Section 3, including all of Section 3.2, we use standard techniques from optimal control ( [21]) and sub-Riemannian geometry ( [18]). In Section 3.3, however, we utilize the optimal control software Bocop. As stated in [27], this software approximates our optimal control problem by a finite dimensional optimization problem using the direct transcription approach to time discretization. The resulting nonlinear programming problem is solved using the software package Ipopt, using sparse exact derivatives computed by ADOL-C.

Equations of Motion.
We first develop the equations of motion for n legs, then specify to the case n = 2. Our equations consist of a system of differential equations of the form Mq = K. The equations of motion for a copepod moving in two dimensions are derived in [11], which focuses on legs moving in an oscillatory motion: θ i (t) = a cos(t + k i ) + β i . Parameters are constrained to ensure that adjacent legs never overlap but possess a phase lag. The author shows that no net rotation is possible with such a motion for two legs, thus most of the analysis concerns the case of three legs. For numerical simulations, the values of a, k 1 , β 1 , β 2 , and β 3 are fixed and the total change in orientation and displacement is computed for varying values of k 2 and k 3 . The change in displacement and orientation is maximized when (k 2 , ). In addition, the total work done by the microswimmer is calculated and a notion of turning efficiency is introduced.
The system is derived from slender-body Stokes flow, using the fact that, at low Reynolds number, inertial forces are negligible and the Navier-Stokes equations can be linearized. Here and M is the resistance matrix given by By computing the mobility matrix M −1 (well defined since M is symmetric and positive definite) we obtain the equations of motionq 2.3. Optimal Control Framework. We now consider this system from the control theoretic point of view, where the angular velocities of the leg are taken as controls. That is, we set u i =θ i , and assume these are measurable functions of time. Now let Then our control vector fields are GEOMETRIC METHODS FOR EFFICIENT PLANAR SWIMMING OF COPEPOD NAUPLII 5 where the 1 appears in the i th entry after the M −1 K i entries. Then the copepod motion is described by the driftless affine control systeṁ where n refers to the number of legs. In this work we assume no bounds on the control. In reality, of course, there are limits to how quickly an actual copepod can move its legs. But some of this issue is mitigated by the fact that we will be minimizing some form of energy; see Equation (11) below. We will see in the next section that, when n = 2, this is a controllable system; that is, it is possible to find controls steering the copepod from any given initial state q initial to any given terminal state q final . he rest of this paper will concern the special case of the copepod with two legs, so here we record these equations of motion explicitly. When n = 2, we have Our control system then takes the forṁ where q = (x, y, φ, θ 1 , θ 2 ) t , our controls are u 1 =θ 1 and u 2 =θ 2 , and the control vector fields are Straightforward calculations give In the sequel we let D denote the distribution spanned by F 1 and F 2 . Note that we can obtain F 2 from F 1 by switching the roles of θ 1 and θ 2 .
. This is a general consequence of swimming at low Reynolds number. The system is also invariant under rigid body transformations. It is obvious for translation as the F i do not depend on x or y, but it can also be verified that φ → φ + τ is invariant under x → x cos τ − y sin τ, y → x sin τ + y cos τ .
In this paper we suppose that the copepod seeks to minimize the mechanical energy expended when moving from one position to another. In [8], the authors consider a two-legged copepod moving along an axis without rotation. They describe a realistic but complicated mechanical energy functional, but show that the resulting optimal trajectories are qualitatively similar to those obtained when using the simplified energy Here t 0 is a fixed initial time, while t f is associated to the control u in the following manner. Choose some terminal boundary manifold M 1 ⊆ R 5 , which is closed, and define the target set by M = [t 0 , ∞) × M 1 . Then t f is the smallest time such that (t f , q(t f )) ∈ M , where q(t) is the state trajectory associated to the control u(t).
We therefore choose our cost function to be the energy E given in (11) corresponding to the orthonormal inner product for the two control vector fields, yielding the following optimal control formulation. Provided certain boundary conditions made explicit below, we seek solutions to the dynamical system (8) which minimize the cost (11). Note that this is a sub-Riemannian problem associated to the flat metric ( [5]).
2.4. Maximum Principle. The Pontryagin maximum principle provides necessary conditions for a solution to be optimal. The general statement can be found in the literature [5,1,17,21], and we here only state it for our application.
Consider the optimal control problem stated in the previous section, defined by the dynamics (8) and cost (11). Let j, k ≤ 5 and consider the initial and terminal boundary manifolds

Define the Hamiltonian by
Theorem 2.1 (Maximum Principle). Let u * : [t 0 , t f ] → R 2 be an optimal control and let q * : [t 0 , t f ] → R 5 be the corresponding optimal state trajectory. Then there exists a function p * : [t 0 , t f ] → R 5 and a constant p * 0 ≤ 0 such that (p * 0 , p * (t)) = (0, 0) for all t ∈ [t 0 , t f ] and having the following properties: (1) q * and p * satisfy Hamilton's equations for the Hamiltonian (12) with boundary conditions Note that this version of the maximum principle closely follows Section 4.1.2 of [17], with the addition of the initial boundary conditions M 0 and the associated transversality conditions. It important, however, to recognize that Theorem 2.1 pertains to the free-time problem. In our simulations in Section 3.2 we fix the time [t 0 , t f ] = [0, 2π], which requires a minor modification of the maximum principle. As described in Section 4.3.1 of [17], one simply introduces an extra state variable to represent time, q 6 = t, and includes the fixed terminal time in the terminal manifold M 1 . Then the free-time Theorem 2.1 applies with only the following modification: H| * = −p * 6 , which is constant. In the sequel we will refer to a trajectory satisfying the conclusions of the maximum principle as an extremal.

Results
Most of the results in this section concern the two-legged copepod. Unless explicitly stated otherwise, we assume n = 2. The motivation for this choice is contained in the next theorem, which says that one leg is insufficient for producing rotation via periodic strokes, but two legs are sufficient. As in [6], we define a stroke to be a periodic deformation of the swimmer's body. That is, a stroke of period T is any path in configuration space satisfying θ i (0) = θ i (T ) for all i = 1, . . . n. For simulations and examples we also impose the following realistic constraint forcing each leg to stay on one side of the copepod's body: However, most of the mathematical analysis in this section is valid for the configuration space R 2 × (S 1 ) 3 . Note that [11] proves that a two-legged copepod is incapable of producing net rotation via the specific oscillating strokes considered in their work.
Theorem 3.1. A one-legged copepod moving in strokes can neither produce a net rotation nor net displacement. A two-legged copepod moving in strokes can produce net rotation.
Proof. For a copepod with one leg, we compute Then the equationq = M −1 K simplifies to Suppose the copepod moves in strokes of period T , so θ(0) = θ(T ). Without loss of generality assume φ(0) = 0. Then φ(t) = − 1 9 θ(t) + 1 9 θ(0). Because θ(t) is periodic, φ(t) is also periodic and thus no net rotation is produced after one period. Now we also have that α(t) = θ(t)+φ(t) = 8 9 θ(t)+ 1 9 θ(0). To calculate the total x displacement over the period we compute the integral Computing this integral gives us which equals zero again since θ is periodic. Thus the displacement in the x direction is 0; a similar argument gives the result for y. The proof of the second statement in the theorem is contained in the following example.
The next example demonstrates that two legs moving in strokes can produce a net change in displacement and orientation.
We consider the following motion. Both legs rotate π radians counter-clockwise in π time; then one at a time, each leg moves π radians clockwise in π 2 time. See Figure 4. Explicitly, we have We can explicitly solve the equations of motion for the orientation over time, which implies a net rotation of φ(2π) − φ(0) = − π 6 = 0. The orientation over time, along with the displacements over time obtained by numerical integration, are shown in Figure 5. The net change in position, equal to the final position, is given bŷ q (2π) = (x(2π), y(2π), φ(2π)) T = (0.0071, 0.0019, −π/6) T .  Note that the energy (11) for the motion in Example 3.2 is equal to 6π. While this motion is dynamically valid, it is likely not minimizing the cost. See Section 3.3 for further discussion.
3.1. Controllability. In this subsection we show that the two-legged copepod is indeed a controllable system. That is, given any initial and final configuration, controls exist which steer the copepod from the initial to the final configuration. The main tool here is the Chow-Rashevskii theorem; a formal statement, along with definitions of all the terminology in this subsection, can be found in [18]. The proofs here are essentially just calculations, which we performed using MATLAB and Mathematica.
Our two control vector fields F 1 and F 2 are given by (9) and (10). Denote their iterated Lie brackets (which are too complicated to display) by Note that F 1 and F 2 , and consequently their iterated brackets, only depend on θ 1 , θ 2 , and φ. A computation shows that the five vector fields F 1 , F 2 , F 3 , F 4 , and where ψ = θ 1 − θ 2 . Note that this is only a condition on two of our five variables. Solving (14) yields two sets of solutions: ψ = 2πn and ψ = 2πn ± 2 arctan(2) for n ∈ Z. We therefore have two sets of configurations which are singular for our distribution D: (2) for some n ∈ Z} S 2 = {q | θ 1 − θ 2 = 2πn for some n ∈ Z}.
Proof. The fact that the small growth vector is (2,3,5) at generic points is immediate from the fact that F 1 , F 2 , F 3 , F 4 , and F 5 are linearly independent there. These points are regular since the singular sets are closed. Points in the singular sets are analyzed by computations. Let For points in S 1 , we find that F 4 = −F 5 , but F 1 , F 2 , F 3 , F 4 , and F 6 are linearly independent. For points in S 2 , we find that F 3 = F 7 = F 8 = 0 and F 4 = −F 5 and F 6 = −F 9 , but F 1 , F 2 , F 4 , F 6 , and F 10 are linearly independent.
Corollary 3.4. The two-legged copepod system is controllable at all points. At generic points the degree of non-holonomy is 3. At points in S 1 the degree of non-holonomy is 4. At points in S 2 the degree of non-holonomy is 5.
Proof. Controllability follows from the Chow-Rashevskii theorem, as the vector fields F 1 and F 2 Lie generate the tangent bundle at every point. The degree of non-holonomy is simply the length of the small growth vector.
3.2. Abnormal extremals. Abnormal extremals are intrinsic to the dynamics; they do not depend on the cost. It is well known that they play a very important role for the optimal synthesis ( [5]). They correspond to imposing p 0 = 0 in the Hamiltonian (12). It follows that, for our application, the abnormal Hamiltonian is H a (q, p, u) = u 1 p, F 1 (q) + u 2 p, F 2 (q) .
According to the Pontryagin maximum principle, abnormal extremals are curves (q(t), p(t)) satisfying Hamilton's equations for H a as well as Differentiating these equations leads to the additional requirements p, F 3 (q) = 0 (17) The following results partially characterize the abnormal curves.
Proposition 3.5. The horizontal lifts of the singular curves in Figure 6 are projections of abnormal curves for the two-legged copepod. They are the integral curves for the vector field F 1 + F 2 restricted to the singular set S = S 1 ∪ S 2 for the distribution D, and they project to uniform circular motion in the xy-plane. More precisely, let q 0 = q i (0), and take any n ∈ Z and any c 1 , c 2 ∈ R not both zero. Then the curves are (q 1 (t), p 1 (t)) and (q 2 (t), p 2 (t)) are abnormal, where and c x = x 0 + √ 5 6 cos(± arctan 2 + φ 0 ) and c y = y 0 + √ 5 6 sin(± arctan 2 + φ 0 ).
Proof. Note that q i (t) is just the horizontal lift of the naive parametrization of the lines which constitute the connected components of S i projected to the θ 1 θ 2 -plane; these are the colored lines in Figure 6. Thus for any time t we have q i (t) ∈ S i . The curves q i are integral curves for F 1 + F 2 and therefore horizontal. In fact, span{F 1 + F 2 } is the intersection of D and the set of Cauchy characteristics for D + [D, D]. It is also the intersection of D with the tangent bundle T S. It is straightforward to check that the pairs (q 1 (t), p 1 (t)) and (q 2 (t), p 2 (t)) satisfy Hamilton's equations. Observe that p i is constant in the first two components since our control vector fields do not depend on x or y. We also have by construction thaṫ Interestingly, the abnormal Hamiltonian also satisfies ∂H a ∂q (q 2 ,p 2 ) = 0.
Note that p 1 is actually a 2-parameter family of curves, and p 2 is constant and therefore an integral of motion. It is similarly straightforward to check that (q i (t), p i (t)) satisfy abnormal equations (15 -18). Recall that F 3 | S 2 = 0 so equation (17) is satisfied for any p on S 2 . Similarly, on any continuous curve within S we have F 4 = −F 5 and u 1 = u 2 , so equation (18) is satisfied for any p.
Corollary 3.6. Abnormal strokes contained in the singular set for D can produce neither rotation nor displacement.
Proof. The only possible abnormal stroke lying in S would require the legs tracing a segment of one of the colored lines in Figure 6 first forward and then backward. The symmetry of such a stroke prevents any net rotation or displacement. Proposition 3.7. Neither control can be zero along an abnormal extremal.
Proof. Without loss of generality, assume u 2 = 0 along an abnormal extremal (q(t), p(t)), so θ 2 is constant. Assume u 1 = 0; otherwise the system is at a stationary point. Equation (18) then implies that p, F 4 (q) = 0. Differentiating this equation yields which reduces to p, F 6 (q) since u 2 = 0. A straightforward calculation shows that the vector fields F 1 , F 3 , F 4 , F 6 are linearly independent along such curve. Moreover, they all are identically zero in the fifth component, so the only p mutually orthogonal to these four vector fields is of the form (0, 0, 0, 0, p 5 ). Then using (10) and p, F 2 = 0 we also have p 5 = 0. Thus p = 0, which contradicts the maximum principle.
Proof. Note that Proposition 3.5 considers the abnormal curves within S. Assume (q(t), p(t)) is an abnormal curve not lying within S for any time interval. By equations (15 -17), we have that p is orthogonal to F 1 (q), F 2 (q), F 3 (q). But on the complement of S, we have that F 1 , F 2 , F 3 , F 4 , F 5 are linearly independent, so p cannot also be orthogonal to both F 4 and F 5 . Without loss of generality, assume p is not orthogonal to F 4 . Then we can solve equation (18) for Thus we scale H a to obtain a new Hamiltonian in which the controls do not appear at all: H a (q, p) = p, F 1 p, F 5 − p, F 2 p, F 4 . As the abnormal equations have been satisfied by construction, any solution to Hamilton's equations for this Hamiltonian will indeed be an abnormal curve.

Normal Extremals.
Taking p 0 = − 1 2 , our normal Hamiltonian is We analyze the normal extremals indirectly, using the optimal control software Bocop (see Section 2 and [27]). Our investigations have led to two interesting observations regarding the behavior of normal extremals. First, certain trajectories seem to show the copepod moving along a curve in the xy-plane which is a type of Euler elastica. See Figure 7 for one example. In particular, we observe this behavior when fixing the start and end positions in the plane, demanding that the net change in orientation is zero, and demanding that the copepod completes a stroke, with no other imposed boundary conditions. It can be observed that the legs follow a periodic motion, and in turn the orientation of the copepod is periodic as well. The motion in the angular phase plane (θ 1 , θ 2 ) is a perfect ellipsoid within the constraint space, reflecting the symmetry of the motion of the two legs. A possible route to proving that this phenomena holds is suggested by [3]. Our optimal control problem can be translated into a geodesic problem in sub-Riemannian geometry. Our control vector fields F 1 and F 2 have dual momenta P 1 = p, F 1 and P 2 = p, F 2 , and the sub-Riemannian Hamiltonian H sR = 1 2 (P 2 1 + P 2 2 ) generates normal geodesics corresponding to our normal optimal copepod trajectories. These geodesics parametrized by arc length correspond to solutions of Hamilton's equations for H sR (geodesic equations) with energy H = 1/2. These equations could potentially allow us to show that the curvature κ of the projection (x(t), y(t)) satisfies one of the defining differential equations for Euler elastica. The obstacles here are that the computations are unwieldy, and it is not clear how to impose the boundary conditions which seem to lead to elastica-like behavior in the copepod.
Our second interesting observation concerns the triangle T , appearing in the lower right corner of the constraint square (13), consisting of the boundary of the set {(θ 1 , θ 2 ) : 0 ≤ θ 1 ≤ θ 2 − π ≤ π}. Our simulations show that following this triangle is optimal for a copepod desiring to rotate a prescribed amount. More precisely, suppose we specify the net rotation ∆φ but impose no other boundary conditions: we do not specify the start or end points in the plane, or require the motion be a stroke. Then the optimal motion of the legs traces out the triangle T from the top right corner counterclockwise; it may go around T more than once, not necessarily an integer number of times. In fact, it will never go around an integer number of times (which would constitute a stroke). Our observations suggest the following characterization of the motion: just follow hypotenuse: 0 to .5 times around T ( 2π 3 , 5π 6 ] once around T , then hypotenuse: 1 to 1.5 times around T ( 5π 6 , π] twice around T , then hypotenuse: 2 to 2.5 times around T . In any of these cases, the hypotenuse need not be traced out completely. For example, to rotate π/3 radians one would simply traverse half the hypotenuse then stop. See Figure 8 for an example of the third case with ∆φ = π. Note the symmetry of the legs in Figure 8(d), which is implicit in the triangle T itself. Of course, to rotate more than π radians one simply reverses this process (starting at the bottom left and following T clockwise). Note that traveling along the hypotenuse induces no displacement and traversing the complete triangle induces very small net displacement (see Example 3.2). Thus these motions represent optimal swimming for a copepod attempting to rotate any amount without much net displacement. Any rotation amount less than or equal to 2π/3 can be achieved optimally with no displacement at all. Intuitively, this demonstrates the fact that traveling along the hypotenuse gives the strongest possible power stroke for inducing rotation -the legs of the triangle are simply necessary to move the copepod legs back into position for another power stroke in a way that minimizes backwards rotation. A motion which includes the legs of the triangle, as in Figure 8, does require the copepod to move around in the plane, but it returns to nearly its original position.
In Figure 9 we provide a catalog of the type of topological curves in the xy-plane obtained when varying the boundary conditions. The boundary conditions themselves appear in Table 1. Table 1. The boundary conditions corresponding to the paths in Figure 9. In all cases the initial position is (0, 0). In some simulations certain conditions were unspecified, labeled "Free" here. For some simulations we impose the stroke condition that θ i (0) = θ i (2π). The total energy of each path is not imposed, but listed here for comparison.

Discussion and Conclusions
Here we have provided a mathematical model of a swimming copepod nauplius with two legs moving in a plane. This model allows for both rotation and two-dimensional displacement by periodic deformation of the swimmer's body. The system was studied from the framework of optimal control theory, with a simple cost function designed to approximate the mechanical energy expended by the copepod. We have found that this model is sufficiently realistic to recreate behavior similar to those of observed copepod nauplii, yet much of the mathematical analysis is tractable. In particular, we have shown that the system is controllable, but there exist singular configurations where the degree of non-holonomy is non-generic. We have also partially characterized the abnormal extremals and provided explicit examples of families of abnormal curves. Finally, we have numerically simulated normal extremals and observed some interesting and surprising phenomena.
This work suggests a plethora of interesting open problems and directions for future research. First, there are a number of potential generalizations and modifications to our model which may lead to even more realistic behaviors. For example, one can study the model with four or six legs coupled with the appropriate constraints. Real copepods have six legs. Further, this model has the potential to design soft small-scale synthetic robots [8,14,25]. Alternatively, or perhaps additionally, one could work in a three-dimensional environment, which is obviously more realistic. Moreover, there are other reasonable cost functions to consider, including more complicated versions of mechanical energy. Instead, it may be that copepods seek to minimize the time needed to perform a given motion, or the total distance traveled, either to evade predators or capture prey more effectively.
Without generalization, our current model already offers ideas for future research. In particular, it would be quite interesting to find a mathematical, physical, or biological explanation for the observed elastica-like paths, such as the one shown in Figure 7. Even numerical verification that these paths are indeed forms of elastica would be worth pursuing. A potential approach is described in Section 3.3, which in turn leads to other questions of a differential geometric flavor. Our optimal control problem can indeed be cast as the geodesic problem for a particular sub-Riemannian geometry, which appears geometrically interesting. The vast sub-Riemannian literature may yield geometric or metric tools providing deeper insight into the copepod system.   Figure 9. Gallery of simulated normal extremals showing copepod paths in the xy-plane. In all cases the copepod begins at the origin (0, 0) and t ∈ [0, 2π]. Additionally, we specify boundary conditions for x(2π), y(2π), φ(0), and φ(2π); these are given in Table 1. No other boundary conditions are imposed.
While the experimental approach in Section 3.3 led to some interesting observations, the normal extremals are still largely not understood. It would be particularly interesting to explore path-planning for the copepod system. It is also important to recognize that Bocop, like any mathematical software, has limitations, some of which we encountered. In particular, this is local optimization software, and we are not working on a convex optimization problem with one global extremum. Thus, despite the symmetries of the problem, the numerical results were sensitive to transversality conditions (for example, specifying that copepod start at the origin).
Finally, we consider how well our results approximate actual observations of copepods in motion. In [9] the authors discuss copepod swimming and escape behavior, based on observations of their swimming patterns and activity. In particular, Figure 4 in that paper shows a helical pattern projecting onto the xy-plane like an ellipse. This correlates with our motion presented in Figure 8. Figure 8 in [9] depicts escape trajectories for nauplii and copepodis which follow helical patterns that project on the xy-plane as Euler elastica. In [10], the authors observe the positions of the appendages during prey capture and prey handling; in their Figure 5 we see the leg motions are oscillatory and mostly periodic, as in our simulations, during the prey handling phase. The most striking comparison comes with the observed behavior provided in [20]. Indeed, the projection of the 3D swimming motion of the nauplii and early copepodid in their observations provide a similar complexity to our simulated trajectories. Compare our Figure 9 to Figures 1 through 9 in their paper. It is quite remarkable that despite the simplified assumptions made on the number of legs and the cost, our results still capture the essence of swimming behavior for copepods.