The role of symmetry and dissipation in biolocomotion
In this paper we illustrate the potential role which relative limit cycles may play in biolocomotion. We do this by describing, in great detail, an elementary example of reduction of a lightly dissipative system modeling crawling-type locomotion in 3D. The symmetry group is the set of rigid transformations of the horizontal (ground) plane. Given a time-periodic perturbation, the system will admit a relative limit cycle whereupon each period is related to the previous by a fixed translation and rotation along the ground. This toy model identifies how symmetry reduction and dissipation can conspire to create robust behavior in crawling, and possibly walking, locomotion.
The notion of limit cycles is important in biolocomotion because simple periodic behavior is a defining characteristic of walking, running, swimming, flapping flight,… In order to construct realistic mathematical models that exhibit limit cycles, it is helpful to first identify some core mechanisms of limit cycle production. In particular, a biolocomotive gait, such as skipping or crawling, has three primary ingredients:
it is time-periodic;
with each period the body translates and rotates in space;
it is robust to noise and systemic variations.
The combination of these ingredients is known to dynamical system theorists as an -relative limit cycle. Here is the special Euclidean group for (i.e. rotations and translations of ). Such an object is a trajectory of a dynamical system with -symmetry, such that its image is a limit cycle under reduction by . In summary, an -relative limit cycle is just a limit cycle modulo rotations and translations.
Presently much of the literature on biolocomotion concerns the search for limit cycles without addressing the role of symmetry. Such limit cycles are made robust through a mixture of dissipation and the dimension reduction which occurs across the transition maps in hybrid systems. Many of these systems occur in a regime with a mixture of friction and inertial forces. In contrast, the symmetry and reduction theoretic aspects of biolocomotion are very well studied in the geometric mechanics community, but only in regimes where friction or inertial forces dominate [35, 39, 32, 33, 34, 28, 26, 23]. The mixed regime is mostly left unstudied by the geometric mechanics community. The goal of this article is to address these gaps by illustrating a simple example of an invariant system which models crawling via an -relative limit cycle. We combine techniques from geometric mechanics, hyperbolic stability and singular perturbation theory. Through this analysis, one can see how these techniques can be generalized to more sophisticated and realistic models.
1.1. Outline of the paper
We start with an overview of the background and motivation in section 2. In section 3 we review the geometric mechanics of biolocomotion in the viscous dominated and inertial regimes. We find that the use of connections in the middle regime is less natural. This lack of naturality motivates the approach of the present paper, which implements Lagrangian reduction without the explicit use of the mechanical or the Stokes connections. Then, in section 4, we introduce our model of an (unactuated) crawler: a mass-spring system resting on the ground, see Figure 1. We regularize the no-slip and no-penetration conditions imposed by the ground by ‘smearing them out’ over a small region around to smooth viscous friction (c.f. [6, 25]) and smooth potential energies (c.f. [38, 40]). What results from this regularization is a constant dimension, lightly dissipative Lagrangian system. Having described the problem as an ODE on a space of constant dimension, we apply symmetry reduction and smooth dynamical systems theory in section 5. The system is invariant under isometric transformations along the ground, and so we implement reduction by this group of transformations [29, 8]. Under mild regularity conditions (Assumption 7), we can use singular perturbation theory to find a robustly stable equilibrium for the reduced model. Next, under small time-periodic forcing (i.e. actuation of the crawler) this equilibrium persists as a limit cycle in the symmetry reduced model as a result of the persistence theorem [13, 21]. The limit cycle in the reduced space corresponds to a relative limit cycle in the original phase space, (see Figure 2). A phase reconstruction formula gives the phase shift of the lifted, relative limit cycle; this phase shift corresponds to a translation and rotation achieved upon traversing a period of the relative limit cycle. Both relative periodicity and stability are characteristics of biolocomotion, and so we can consider the relatively periodic orbit as a model of crawling in this sense.
These results are summarized by the following theorem.
Theorem 1 (main theorem).
Let a simple crawler model be described by the Lagrange–d’Alembert equations (14). Under Assumption 7 the symmetry reduced system has a stable rest state. For sufficiently small time-periodic forcings, this rest state persists as a stable limit cycle, which corresponds to a relative limit cycle in the unreduced system that models crawling.
We also show that the phase shift of the relative limit cycle depends on the magnitude of the perturbation to second order.
Finally, in section 6 we illustrate the theory with numerical experiments before stating our closing remarks in the conclusion.
2. Background & motivation
Biomechanics requires knowledge from a range of fields. This particular paper draws upon previous research in geometric mechanics, stability theory, as well as inspiration from experimental and numerical observations.
2.1. Contact problems
The regularization we are going to pursue is in contrast to the hybrid systems approach, where transitions between different types of phase spaces are given by various transition maps. The hybrid systems approach expresses the non-constant nature of the dimension in contact problems explicitly, and has yielded a number of insights and useful models. For example, a hybrid systems formulation was introduced by McGeer , where the transition maps led to a dimension reduction; it was suggested that a limit cycle was approached passively. Since the work of , the notion of walking as a limit cycle has become more common, and more sophisticated analyses have lent further support to this idea [15, 16]. The most compelling arguments are the original videos of McGeer which accompany .
2.2. Biology and Engineering
On the biological side, ‘central pattern generators’ (CPGs) have been hypothesized as fundamental neural mechanism used in biolocomotion . These CPGs are non-localized collections of neurons which produce rhythmic activity, and respond to various inputs which modulate these rhythms. Therefore the link between CPGs and limit cycle biolocomotion is one which links periodic activation of the controls to periodic motion of the body. This link is used in the creation of simple models which can be feasibly analyzed (see for example ). A similar regularization of ground contact which will be presented in this paper is used in , which studies robustness and efficiency of a simple passive dynamic walking model actuated by CPGs, although the role of symmetry was not addressed there.
The notion of a CPG is significant from the perspective of biologically inspired control theory because less demand is placed upon the control law when locomotion is achieved primarily through an open-loop control. For example, under weak assumptions, the existence of limit cycles in hybrid systems implies the existence of a reduced order model for the system as a whole .
2.3. Geometric mechanics
There is a long history of using geometric mechanics to study locomotion. Purcell’s three link swimmer  inspired Shapere and Wilczek to interpret locomotion in Stokes flow as phase shift due to the curvature of a principal connection . The simplicity of this perspective has proven useful in other dissipation dominated systems such as granular media . It was later found that a range of examples of locomotion fit within this geometric framework [33, 28]. In particular, many conservative systems could be analyzed in this way [32, 26, 23, 34].
Despite the success of the gauge theoretic picture of locomotion, the vast majority of examples of this perspective concern systems which are either conservative (i.e. Hamiltonian or Lagrangian), or friction dominated (i.e. where Newton’s law, , is replaced by ). There appear to be very few examples which invoke the gauge theoretic perspective of  in a regime which exhibits a mixture of viscous and inertial forces. The paper  by Kelly and Murray is a notable exception, where they discuss both mechanical and Stokes connections and study control of systems in the intermediate regime by viewing friction as a drift term to the system with mechanical connection. They also note that the mechanical and Stokes connections cannot simply be interpolated. This gap between the inertial and viscous regimes is one which the current paper seeks to address.
Again, the middle regime has been shown to be more than merely an interpolation
between the two extreme regimes.
For example, the scallop theorem states that a system in the viscosity dominated regime with only
one degree of freedom in shape space cannot move.
As a final point, the role of symmetry is well acknowledged within the geometric mechanics literature but this is not to say that it is absent from the biomechanics literature, for example the importance of discrete symmetries of solutions is widely acknowledged [37, 36]. However, tools such as momentum maps and connections are typically not used. This is possibly due to the fact that geometric mechanics has only addressed the extreme regimes, while many popular biomechanical models fall within the middle regime. A notable exception is , where the Noetherian momentum associated to an symmetry was used to create turning trajectories based upon the work of [3, 4]. Here we will be exploring a different, but related, application of symmetry reduction where Noether’s theorem is never invoked (nor does it apply).
3. The gap between the mechanical and Stokes connections
In this section we will explore the traditional use of connections in understanding locomotion. We will find that the use of connections is unmotivated when there is a mixture of inertial and viscous forces. This section is aimed at an audience which is familiar with these more established techniques. As the goal of the section is to illustrate why we should not use these tools, we recommend that the reader should skip this section if she is unfamiliar with the use of connections in locomotion, at least upon a first reading.
The typical setup of the configuration space in biolocomotion is that of a principal -bundle, where is the (spatial) symmetry group of the system and generates the directions in which locomotion can take place. That is, we have a configuration space and a left -action on , such that the quotient projection is a left -principal bundle. The base is conventionally called the ‘shape space’, as it describes the state of the system modulo its position.
Assuming that the system is symmetric under , we can consider the reduced dynamics on . This bundle is locally isomorphic to and is naturally a vector bundle over ; in the fibers, all velocities are retained, since the dynamics may still depend on these velocities even though it does not depend on the underlying points in the orbits of .
Given a connection, this can be decomposed in a vector bundle sum
where is the adjoint bundle and a natural vertical distribution in , while the connection is used to identify with a horizontal distribution.
Now there are two natural and useful choices of connection in the limit cases (see e.g. ): the mechanical connection, for when the dynamics is Hamiltonian, and the Stokes connection in case of a friction dominated limit, i.e. high Reynolds number in swimming-like locomotion or high Froude(-like) number in terrestrial, finite-dimensional locomotion models. We shall see that in the middle regime there is generally no natural choice of connection in order to understand locomotion. This lack of a natural choice will motivate us to avoid the use of a connection later in the paper.
The goal of this section is to illustrate this inability to address the middle regime. We shall begin by discussing the general setup of a dissipative mechanical system on a manifold before describing the role of connections in understanding locomotion.
3.1. Lagrangian mechanics with dissipation
Let us briefly return to a setup without symmetry assumption. A mechanical system with (viscous) friction can be represented using a Lagrangian to model the conservative part and a Rayleigh function to model friction. Let
denote the Lagrangian with kinetic energy metric which turns into a Riemannian manifold, and with potential . Let
denote a Rayleigh dissipation function that defines a friction force given by minus its fiber derivative, see [1, Def. 3.5.2]. Note that is assumed to be a positive definite quadratic form on that depends smoothly on , hence is a metric on , like . The parameters and will allow us to consider the Hamiltonian and dissipation dominated limits.
Now the Lagrange equations of motion are given by
to which an extra force can be added on the right-hand side. When we take the limit of , we straightforwardly converge to a conservative system. When , a more careful singular perturbation analysis shows (using the assumption that is positive definite, see Appendix A) that
is a well-defined, attractive invariant manifold for the limit dynamics. We can view as a submanifold of , but since is the graph of a section of , we can also view it as a vector field that generates first order dynamics on . This is the Stokesian limit
3.2. Mechanical connections
If is a kinetic energy metric on that is invariant under , then it defines a connection on by defining the horizontal space complementary to as its perpendicular under . This can be viewed as a sub vector bundle and descends through the quotient by to a sub vector bundle of . Note that this connection is not the Levi-Civita connection defined by , since the Levi-Civita connection induces a splitting of , at one higher level.
Now if the complete system is invariant under the (lifted) action of , then by symmetry, is an invariant submanifold for the dynamics, and it is exactly the level set of zero momentum under the (Lagrangian) momentum map induced by the action of . If we use this connection for the identification (1), then it implies that the component is constantly zero, hence we can reduce to dynamics on , and after solving that, lift solution curves to and even , using the mechanical connection and integration of the vertical component, respectively.
In this case the equations of motion (2) reduce to
and exhibit as invariant submanifold.
3.3. Stokes connections
The Stokes connection is defined in the same way as a mechanical connection, but now using the metric on . In this case we assume that friction forces dominate the inertial forces and, hence, the dynamics is only first order, and given by (4). Next, it is typically assumed, e.g. in Stokes flow swimming, that there is an external force exerted by the swimmer, say, which physically implies that , the annihilator of . Since was assumed invariant it follows that too, and using the fact that by definition, we obtain the control system
where and . In particular, given a curve , we can lift it to a curve using the Stokes connection and (6). This corresponds to the unique motion in such that no work is done in the directions of the symmetry . See also  and more details in Appendix A.
3.4. The middle regime
Let us now return to study dynamics on in the middle regime where both inertial and frictional forces are present, i.e. neither nor is negligibly small. We can rewrite the equations of motion (2) as
Since it follows that the first right-hand side term lives in , hence that part of the dynamics preserves the splitting . The mapping
however, will generally not preserve this splitting, so the mechanical connection does not yield a reduction here. If is sufficiently small, we can actually still reduce to a first order system. This can be considered the ‘perturbed Stokes regime’ where the first order ODE (4) does not accurately hold anymore, but approximations can be found using singular perturbation theory. The corrections, however, cannot be interpreted as a connection anymore, see Appendix A for the details.
The lack of a natural connection in the middle regime suggests that we should implement reduction by symmetry without the use of a connection. In the language of geometric mechanics, this means we shall derive equations of motion on the Atiyah algebroid , as in , as opposed to a Lagrange–Poincaré bundle , as in .
4. The model
The model can be broken into two distinct components: the crawler and the environment. The crawler consists of four masses connected by springs while the environment consists of the ground and a gravitational field. We will discuss the model of the crawler in empty space before we elaborate on how to include interactions with the environment.
4.1. A model of a crawler (in a vacuum)
The crawler consists of four point particles of unit mass all connected by springs of stiffness with light viscous damping , see Figure 1. We describe the crawler as a Lagrangian mechanical system with additional forces to model the spring damping. The point particles move through space with positions and velocities for . For reasons to be clarified shortly, we will exclude configurations where any of the particles overlap, and the configurations where all of the coordinates overlap. Thus the configuration space is a (dense) open set . We will use “” to denote a generic point of and to denote generalized coordinates of .
The kinetic energy is given by with the usual Euclidean metric. This endows with a flat Riemannian structure, and we will denote the Riemannian metric by or in generalized coordinates. The potential energy from the springs, , is more easily expressed in other coordinates: the spring lengths, i.e. the pair-wise distances between the points . We therefore introduce six (local) coordinate functions
for and . The potential energy of the springs is now simply given by
where are constants which denote the rest length of spring between and .
We define the viscous force of each spring by the one-form
In terms of the usual coordinates these six forces can be written as a sum of twelve force vectors describing the force exerted on particle by the viscous friction of the spring connecting it to particle . We have
where is the unit vector pointing from mass to mass .
The expression (9) constitutes the components of (8) with respect to the standard basis one-forms .
More precisely, if we denote the components of by and , then the sum is the one-form acting upon mass
Later in the paper we will make the rest lengths time dependent as a means to indirectly control the actual lengths of the springs. Upon performing the substitution by functions , one should be careful about what kind of system is modeled by the resulting equations of motion. In our case, one could imagine that the viscous damping is realized through the addition of dashpots being placed in parallel to the springs.
4.2. A regularized model of the ground
The ground is described by the plane in . Ideally, the ground is impenetrable and imposes a no-slip condition, mathematically represented by the constraints
for , where equation (10) is the no-penetration condition and equation (11) is the no-slip condition.
Both conditions present challenges of a singular nature because they abruptly ‘turn on’ at and are inactive otherwise.
It is precisely this ‘on/off’ character which we will regularize.
To do this we will repeatedly make use of the differentiable
to construct forces and potentials.
We approximate the no-penetration condition by considering a potential energy that grows rapidly for each and is zero when for . Therefore, we define the potential energy by
This penalizes particles for falling through the floor and the penetration depth for a particle at rest can be controlled with . When approaches infinity, the penetration depth goes to zero and our model approaches an exact model of a perfectly impenetrable ground. This can be viewed as modeling a one-sided holonomic constraint in the spirit of [38, 40]. A more advanced version of such an approach is used in  to model contact problems with accurate simulations without being slaved to using infinitesimal time-step sizes.
The no-slip condition is similar to the no-penetration condition in that it is only active at . However, unlike the no-penetration condition, the no-slip condition is not derivable from a potential energy but instead can be viewed as a limit of viscous friction [6, 25]. In particular, consider the viscous force given by
The force dampens the horizontal motion of particles in a region around . Moreover, we can see that is proportional to , the normal force exerted by the ground. This is consistent with standard (first-order) assumptions about the nature of slip-friction. As before, the coefficient controls the amplitude of this force and when goes to infinity we arrive at a no-slip condition.
Similarly, we dampen bouncing at the impact of a particle with the ground by including the viscous friction force
Finally, we incorporate gravity via the potential energy
which imposes the gravitational force .
4.3. The full model
Now that we have established the Lagrangian of the crawler, as well as the environmental forces imposed on it, we can finally provide the equations of motion. These equations of motion are obtained by adding the viscous forces, and , and the potential forces, and , to the equations for the crawler in a vacuum. Adding these up into the total potential energy and the total force , the equations of motion are the Lagrange–d’Alembert equations,
where . This equation implicitly determines given and . We can make this expression more explicit by writing it in the form
where denotes the cometric and the Christoffel symbols associated to the metric . Note that is linear in the velocity and can be written as for a positive semi-definite quadratic form given by
In this section we prove the existence of a robustly stable equilibrium in a symmetry reduced phase space. To begin, we review the general process of reduction by symmetry before handling the specific case at hand. We reduce our system by an symmetry to obtain a reduced vector field on the reduced phase space . Subsequently, we prove the existence of a robustly stable equilibrium which can then be periodically perturbed to obtain a limit cycle. We reconstruct from it a relative limit cycle in the unreduced system. Finally, we provide some illustrative numerical results to support our claim that the reconstructed relative limit cycle typically has a non-trivial phase shift.
5.1. Reduction by symmetry in general
The notion of reduction by symmetry in dynamical systems is conceptually very simple. If a system is invariant under a group of transformations, then it is, in some sense, more simple than a system which is not invariant.
Simply put if is a vector field on and is a Lie group which acts on , then we say that is invariant under if for all and . Here acts on by the tangent lift of the action on . For example if and acts on by rotation about the origin, then a vector field is invariant if it is of the form .
In our example, if , then is coordinatized by the radius and . We see that describes dynamics on a smaller space (), yet still captures all of the richness of . Determining from a vector field is known as reduction by symmetry. In the next section we will perform reduction by symmetry with respect to the group of rotations and translations of the plane.
Finally, if is an equilibrium of and is obtained via reduction by symmetry, then is an equilibrium of . Moreover, the linearization of about is related to the linearization of about .
Assume acts freely and properly on , and let denote the quotient projection. Let be a fixed point of . If is invariant, then is a fixed point of the reduced vector field and the linearization of about is given by , where is an arbitrary right inverse to .
Assume the setup of Proposition 2. Then the kernel of is a subset of the kernel of .
Let denote the flow of the vector field . As a consequence of [2, Prop. 4.2.4] we know that is -invariant when is -invariant. If is in the kernel of then it must be of the form for some curve which originates at . We find
Therefore, is the identity on the subspace of tangent to a -orbit. Taking the time derivative we find that must evaluate to on the subspace of tangent to a -orbit. ∎
Proof of Proposition 2.
By the commutative relation between and above we observe that is a fixed point of . As is surjective, we may define the formal inverse . By Lemma 3, , so that is a well-defined map from . In other words . We may now replace the formal inverse, , with an arbitrary right inverse, to conclude the proof. ∎
5.2. Reduction by
The group consists of all isometries of the plane, i.e. rotations and translations of . Elements of are given by an angle, , and a translation vector . We consider the action of on that translates and rotates the coordinates of each of the masses. That is, we define the action
Note that does not act freely on all of : if all masses
have the same coordinates, then is fixed by the isometries which rotate the plane about .
and we find:
The action of on is free and proper.
The action is free if implies that . Since the collection of coordinates of the points are prohibited from completely overlapping, it must be the case that fixes a non-degenerate line segment in the plane. The only such isometry which satisfies this constraint is the identity, .
To prove that the action is proper, we have to show that
is a proper continuous map, see [11, p. 53]. We shall do so by proving that has a continuous inverse, defined on its image. Let and without loss of generality assume that .
Then where and . Firstly, the angle depends continuously on the arguments vectors, since these have non-zero lengths. Secondly, the translation depends continuously on and the other arguments, where denotes the matrix of rotation over an angle . ∎
As this action is free and proper we can assert that the quotient space, , is a manifold, and is a principal bundle. In order to understand the principal bundle structure of it is useful to find a coordinate system in which the map is a Cartesian projection. Let us consider the (local) coordinates where
In words, is the average of the mass positions in the plane and is the angle between the line segment from to and the -axis. In these coordinates the action of is given by
These coordinates locally trivialize as a principal bundle in that the quotient projection simply drops the last three coordinates , and the space is a nine-dimensional space with coordinates .
The action on naturally lifts to a free and proper action on the tangent bundle, , given by
As before, we find that is a smooth manifold
and we obtain a (left) principal bundle .
These are moving frame coordinates, and the coordinates are sometimes called ‘pseudo-coordinates’ since they are not induced by coordinates on . In terms of the local trivialization , we see that form standard induced coordinates on and are moving frame coordinates on induced by left-trivialization of .
In these coordinates the left action of on is naturally given by
We can immediately see that the quotient projection merely projects out the coordinates, i.e.
On the open subset of where the coordinates are valid, the map is a principal connection if we identify as an element of . However, unlike the mechanical or Stokes connections, this map is not derived from physical properties of the system. It is merely derived from a non-canonical choice of coordinates that locally trivialize the principal bundle .
Recall that is a vector bundle over . The coordinates of can be given by and fibers of are parametrized by the coordinates . Similarly, the principal bundle is a vector bundle with base coordinates and fibers coordinates . We see that is linear in the fiber coordinates (in fact it is the identity on the fibers with respect to these coordinates), and we could say that inherits the vector bundle structure of through the map . We denote by the vector bundle dual to ; this dual vector bundle will come into play shortly.
Now that we understand , we wish to assert the existence of a unique dynamical system on which is consistent with the dynamical system on given by (14).
Note that (14) is written in terms of the total potential energy and the total dissipative force . We observe that is invariant because , , and . As a result, there exists a unique reduced potential such that . It is easy to believe that the differential which appears in (14) must be invariant as well; to understand this invariance we must consider how acts upon .
In a natural sense, the left action of on induces a right action on . In standard Cartesian coordinates for the masses we may consider the fiber coordinates on , in which case the action is given by
for . We may also consider this action in terms of the coordinates where are fiber coordinates conjugate to the fiber coordinates on . The action on is expressed in these coordinates as
We say that is invariant if for any and . In other words, is invariant if the following diagram commutes
for any . In coordinates takes the form
As is only a function of and we see that and are only functions of and as well. The group acts trivially on the variables and and therefore we find
Applying to this we indeed verify that is invariant. This invariance implies the existence of a unique map , explicitly given in coordinates by
Now that we have verified the invariance of , we must do the same for the dissipative force . If is invariant under the action, then we should find that for all and . In other words, is invariant if the following diagram commutes
for any . Intuitively, it is obvious that and are invariant because they are only functions of and , upon which acts trivially. The force is more subtle to analyze. We find that
Thus is invariant and therefore the total force is invariant. An equivalent statement of being invariant would be that acts by isometries with respect to the metric . In any case, invariance of implies the existence of a unique map such that