Simulating the evolution of bipedalism and the absence of static bipedal hexapods

We adopted Simbody as our simulated physical environment, which is a toolset of SimTK, a free platform for the biomedical computation community. It provides general multibody dynamics capability solving Newton's second law in any set of generalized coordinates subject to arbitrary constraints [15]. It aims at the needs of biomedical researchers to analyze the dynamics of biological structures and make clinically meaningful recommendations, and it is also applicable to other fields of research including neuromuscular, prosthetic, and biomimetic simulations [16–18]. Simbody has a 'force toolset' that acts on objects to produce motions, which contains subsystems useful for modeling gravity, contact, molecular forces, and general mechanical elements such as springs (harmonic restraints), dampers, and motors [17]. Such features are important to allow us to implement VMC on the joints of our virtual creature. Simbody is provided as an open-source, object-oriented C++ API and maintained on GitHub, at https://github.com/simbody.

VMC is a motion control language that emulates mechanical components, like simple springs, dampers, dashpots, masses, latches, to create forces, which are applied through real joint torques [19]. These joint torques create the same effect that the virtual components would have created, thereby creating the illusion that the virtual components are connected to the robot. Due to the intuitive nature of this technique, applying VMC on a simple set of virtual components is useful for many complex tasks that are difficult to describe with traditional techniques [12].

In figure 1, we show the virtual model created in the Simbody simulation environment, and table 1 provides descriptions for the parameter names of the virtual model. First, we will present the mathematics to implement virtual components on our creature. After that, we will demonstrate the virtual components designed for the standing phase and the step phase.

Zoom In Zoom Out Reset image size

Our virtual creatures have three joints on each leg [19, 20], as shown in figure 2. The h, k and a subscripts refer to the hip, knee and ankle respectively, and l and r means left and right leg. θ (as shown in equation (1)) is a vector composed of angles of ankle, knee, and hip of both legs. The position of the hip is represented by x, as shown in equation (2). The relationship of x and θ is given by the equation (3). The Jacobian matrix is calculated in equation (4). According to the principle of virtual work [21], we perform virtual displacement on the end effector, and since the sum of the virtual work is 0 we can promptly arrive at equation (5). The joint torque can be calculated from this virtual force, as shown in equation (6):

$$\begin{equation}\boldsymbol{\theta }={\left[\begin{matrix}\hfill {\theta }_{la}\hfill & \hfill {\theta }_{lk}\hfill & \hfill {\theta }_{lh}\hfill & \hfill {\theta }_{ra}\hfill & \hfill {\theta }_{rk}\hfill & \hfill {\theta }_{rh}\hfill \end{matrix}\right]}^{\mathrm{T}}\end{equation} \tag{ 1 }$$

$$\begin{equation}\mathbf{x}=\left[\begin{matrix}\hfill {x}_{l}\hfill \\ \hfill {y}_{l}\hfill \\ \hfill {\phi }_{l}\hfill \\ \hfill {x}_{r}\hfill \\ \hfill {y}_{r}\hfill \\ \hfill {\phi }_{r}\hfill \end{matrix}\right]=\left[\begin{matrix}\hfill -{L}_{1}\enspace \mathrm{sin}({\theta }_{la})+{L}_{2}\enspace \mathrm{sin}({\theta }_{lk}-{\theta }_{la})\hfill \\ \hfill {L}_{1}\enspace \mathrm{cos}({\theta }_{la})+{L}_{2}\enspace \mathrm{cos}({\theta }_{lk}-{\theta }_{la})\hfill \\ \hfill -{\theta }_{lh}+{\theta }_{lk}-{\theta }_{la}\hfill \\ \hfill -{L}_{1}\enspace \mathrm{sin}({\theta }_{ra})+{L}_{2}\enspace \mathrm{sin}({\theta }_{rk}-{\theta }_{ra})\hfill \\ \hfill {L}_{1}\enspace \mathrm{cos}({\theta }_{ra})+{L}_{2}\enspace \mathrm{cos}({\theta }_{rk}-{\theta }_{ra})\left.\right)\hfill \\ \hfill -{\theta }_{rh}+{\theta }_{rk}-{\theta }_{ra}\hfill \end{matrix}\right]\end{equation} \tag{ 2 }$$

$$\begin{equation}\mathbf{x}=\boldsymbol{J}\boldsymbol{\theta }\end{equation} \tag{ 3 }$$

$$\begin{align}\boldsymbol{J}& =\frac{\partial \mathbf{x}}{\partial \boldsymbol{\theta }}=\left[\begin{matrix}\hfill \frac{\partial {x}_{l}}{\partial {\theta }_{la}}\hfill & \hfill \frac{\partial {x}_{l}}{\partial {\theta }_{lk}}\hfill & \hfill \frac{\partial {x}_{l}}{\partial {\theta }_{lh}}\hfill & \hfill \frac{\partial {x}_{l}}{\partial {\theta }_{ra}}\hfill & \hfill \frac{\partial {x}_{l}}{\partial {\theta }_{rk}}\hfill & \hfill \frac{\partial {x}_{l}}{\partial {\theta }_{rh}}\hfill \\ \hfill \frac{\partial {y}_{l}}{\partial {\theta }_{la}}\hfill & \hfill \frac{\partial {y}_{l}}{\partial {\theta }_{lk}}\hfill & \hfill \frac{\partial {y}_{l}}{\partial {\theta }_{lh}}\hfill & \hfill \frac{\partial {y}_{l}}{\partial {\theta }_{ra}}\hfill & \hfill \frac{\partial {y}_{l}}{\partial {\theta }_{rk}}\hfill & \hfill \frac{\partial {y}_{l}}{\partial {\theta }_{rh}}\hfill \\ \hfill \frac{\partial {\boldsymbol{\phi }}_{l}}{\partial {\theta }_{la}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{l}}{\partial {\theta }_{lk}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{l}}{\partial {\theta }_{lh}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{l}}{\partial {\theta }_{ra}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{l}}{\partial {\theta }_{lk}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{l}}{\partial {\theta }_{rh}}\hfill \\ \hfill \frac{\partial {x}_{r}}{\partial {\theta }_{la}}\hfill & \hfill \frac{\partial {x}_{r}}{\partial {\theta }_{lk}}\hfill & \hfill \frac{\partial {x}_{r}}{\partial {\theta }_{lh}}\hfill & \hfill \frac{\partial {x}_{r}}{\partial {\theta }_{ra}}\hfill & \hfill \frac{\partial {x}_{r}}{\partial {\theta }_{rk}}\hfill & \hfill \frac{\partial {x}_{r}}{\partial {\theta }_{rh}}\hfill \\ \hfill \frac{\partial {y}_{r}}{\partial {\theta }_{la}}\hfill & \hfill \frac{\partial {y}_{r}}{\partial {\theta }_{lk}}\hfill & \hfill \frac{\partial {y}_{r}}{\partial {\theta }_{lh}}\hfill & \hfill \frac{\partial {y}_{r}}{\partial {\theta }_{ra}}\hfill & \hfill \frac{\partial {y}_{r}}{\partial {\theta }_{rk}}\hfill & \hfill \frac{\partial {y}_{r}}{\partial {\theta }_{rh}}\hfill \\ \hfill \frac{\partial {\boldsymbol{\phi }}_{r}}{\partial {\theta }_{la}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{r}}{\partial {\theta }_{lk}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{r}}{\partial {\theta }_{lh}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{r}}{\partial {\theta }_{ra}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{r}}{\partial {\theta }_{rk}}\hfill & \hfill \frac{\partial {\boldsymbol{\phi }}_{r}}{\partial {\theta }_{rh}}\hfill \end{matrix}\right]\\ & =\left[\begin{matrix}\hfill A\hfill & \hfill B\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill C\hfill & \hfill D\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill E\hfill & \hfill F\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill G\hfill & \hfill H\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill -1\hfill \end{matrix}\right]\end{align} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

where

$$\begin{align*}A& =-{L}_{1}\enspace \mathrm{cos}({\theta }_{la})-{L}_{2}\enspace \mathrm{cos}({\theta }_{lk}-{\theta }_{la})\\ B& ={L}_{2}\enspace \mathrm{cos}({\theta }_{lk}-{\theta }_{la})\\ C& =-{L}_{1}\enspace \mathrm{sin}{\theta }_{la}+{L}_{2}\enspace \mathrm{sin}({\theta }_{lk}-{\theta }_{la})\\ D& =-{L}_{2}\enspace \mathrm{sin}({\theta }_{lk}-{\theta }_{la})\\ E& =-{L}_{1}\enspace \mathrm{cos}({\theta }_{ra})-{L}_{2}\enspace \mathrm{cos}({\theta }_{rk}-{\theta }_{ra})\\ F& ={L}_{2}\enspace \mathrm{cos}({\theta }_{rk}-{\theta }_{ra})\\ G& =-{L}_{1}\enspace \mathrm{sin}{\theta }_{ra}+{L}_{2}\enspace \mathrm{sin}({\theta }_{rk}-{\theta }_{ra})\\ H& =-{L}_{2}\enspace \mathrm{sin}({\theta }_{rk}-{\theta }_{ra})\end{align*}$$

$$\begin{equation}\left.\begin{matrix}\hfill {\boldsymbol{\tau }}^{\mathrm{T}}{\Delta}\boldsymbol{\theta }+{(-\mathbf{F})}^{\mathrm{T}}{\Delta}\mathbf{x}=0\hfill \\ \hfill {\Delta}\mathbf{x}=\boldsymbol{J}{\Delta}\boldsymbol{\theta }\hfill \end{matrix}\right\}{\Rightarrow}\boldsymbol{\tau }={\boldsymbol{J}}^{\mathrm{T}}\mathbf{F}\end{equation} \tag{ 5 }$$

$$\begin{equation}\boldsymbol{\tau }={\boldsymbol{J}}^{\mathrm{T}}\mathbf{F}=\left[\begin{matrix}\hfill A\hfill & \hfill C\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill B\hfill & \hfill D\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill E\hfill & \hfill G\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill F\hfill & \hfill H\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \end{matrix}\right]\left[\begin{matrix}\hfill {F}_{lx}\hfill \\ \hfill {F}_{ly}\hfill \\ \hfill {F}_{l\phi }\hfill \\ \hfill {F}_{rx}\hfill \\ \hfill {F}_{ry}\hfill \\ \hfill {F}_{r\phi }\hfill \end{matrix}\right].\end{equation} \tag{ 6 }$$

We take a four-legged virtual creature as an example on how we implemented virtual components for the standing and stepping phases. We adopt a similar structure to the one used in [12, 19, 22] and obtain the current virtual model shown in figure 1. As for the standing phase, the requirements are stable and remaining in an upright standing position without leaning too far forwards or backwards. The springs k₂ and damper λ₂ acting at hip height bring stability to the model, and facilitate parameter tuning. As for the step phase, the vertical wall to which k₁ and λ₁ are attached is required to move to the forward direction over time and make the creature lean in the direction of motion. The spring k_step only works in step phase, and the spring k_bodyStraight and dashpot λ_bodyStraight are only valid in standing phase, as the body is constrained in this degree of freedom.

The virtual model for the standing and step phases alone are not enough to make the virtual creature walk, and we need a high-level controller to cascade with this low-level VNC and to modulate the parameters of the virtual mechanisms [19].

We use a state machine to output discrete commands and alternate between states, which will result in a fluid motion of our virtual creature. In figure 3, we show the state machine for each leg pair, with three states in total and initiated in the 'standing on both feet' state. After the stability check, we move into the 'move right foot forwards' or 'move left foot forwards' state, which changes virtual component connections of standing phase to that of step phase as discussed in the last subsection. The choice between left and right is made by which foot is in front of the other, to then alternate at the following standing phase. The state machine cycles between step phase and standing phase to make the discrete states coherent and the creatures walk forward. Virtual creatures are defined by a set of parameters and evaluated by the walking distance. As for how to generate and select virtual creatures with different parameters, we will describe them in the next subsection.

Zoom In Zoom Out Reset image size

Genetic algorithms are a family of computational models inspired by evolution, and they have wide applications in different research fields [23, 24]. The evolution process of the genetic algorithm has three steps: selection, crossover and mutation. The following is the description of the flowchart of genetic algorithms adapted to our experiment, as shown in figure 4.

Zoom In Zoom Out Reset image size

First step is to initialize, we will have an empty array named 'fittest' and an array of creatures mutated from the default one. As for mutation, parameters of each virtual model include an associated probability of mutation p, an associated standard deviation σ, and a current value μ. The result of a success/failure experiment will be used to help decide whether the creature mutates or not. If the result is true, then the new parameter is chosen from the Gaussian distribution N ∼ (μ, σ²), otherwise, the parameter remains the same. After initialization, the next step is the evolutionary loop. We first obtain the fitness of all creatures through simulation and then select the creature with the higher fitness to replace inferior ones in the 'fittest' array while maintaining the specific number. Regarding selection, instead of choosing a subset directly from the current generation, we use the 'fittest' array that contains the fittest creatures ever produced to increase success rate. The algorithm then loops until reaching the specific number of evolutionary generations. There are no crossovers between different species in this experiment.

In the previous section we introduced the method and design of the simulation. In this section, we will propose our setup and specific parameter settings for the genetic algorithm. The purpose of our simulation is to explore the locomotion transition from octapods to bipods. We treated the parameter legPairs as an explorable parameter for the virtual creatures (the mutated value was rounded to the nearest integer), and the starting creatures were eight legged. The reason for starting with octapods is that we do not know in advance whether there will be six legs or four legs when multi-legged creatures evolve to bipedalism. We cannot directly start with hexapods, as it will create a natural bias toward hexapodal evolution and possibly taint final results. The generation size is sixteen and the initial fitness of the default creature was 2.3 (reaching 2.3 in its trail). The probability of legPairs mutating was 30% and the standard deviation was 2.5. We experimented with two conditions, one where the parameter legPairs presents octopods, hexapods, tetrapods and bipods as evolutionary options (W4L, standing for 'with four legs available') and one with octopods, hexapods and bipods (without four legs, or WO4L). The setting for both cases is the same except artificially preventing the emergence and evolution of four legs in WO4L case. We use Intel Core i7-8700 CPU with 12 cores and 32 GB RAM.

Simulations were performed with a maximum step size of 1/30 s, g = 9.81 ms⁻², creature material stiffness equals to 1 × 10⁻⁸ Nm⁻¹, static and dynamic friction coefficients of 0.3 and a creature material dissipation coefficient of 1 s² m⁻¹. The density of the creature material used was 1180 kg m⁻³. These coefficients were chosen to be the same as the ones used in the examples given by the authors of Simbody.

Simulations for all experiments were performed for 5 s, after which the simulation ended and the final position of the farthest rear foot was taken as the distance traveled, i.e. the fitness of the robot. Simulations were terminated early if the robot fell down, which means the upper body fell to below 0.5 of the total robot height or if the simulation terminated when running for more than 200 s. The supplementary material provides a video showing the simulation process (https://stacks.iop.org/BB/16/056012/mmedia).

When evolving toward bipedalism we found that the W4L case showed a higher probability of success than the WO4L case, 70% and 10% respectively, as shown in figures 5 and 6. Thus, in the evolution of bipedalism with tetrapods and hexapods, tetrapods were shown to be more conducive to the generation of bipedalism. This finding is consistent with theories developed by biologists, as in a comparative study of infants and olive baboons [25] which pointed out that the tetrapodal pattern of primates does not impede the developmental acquisition of bipedal behaviors, and a secondary locomotor mode (e.g. bipedalism) experienced during infancy as a by-product of locomotor development may lead to evolutionary novelties when under appropriate selective pressures. Although the aforementioned study focused on bipedalism in tetrapods [26] without hexapods, the complete absence in nature of a truly bipedal hexapod might also be used as an evidence to reinforce our findings. In [27], researchers present a preliminary bipedal running experiments with their Robotic Hexapod, RHex, and reach stable running at an average 1.08 m s⁻¹ with a success rate of 59%. Here, we take the bipedal hexapods walking into consideration, which might be more difficult than bipedal running when the center of gravity is taken in consideration (e.g. running cockroachers).

Figure 7 shows femur length changes of different limbs creature in the whole generation and around the transition zone. It can be seen that the femurLength parameter in W4L successful case (figure 7(a)) has more apparent declines and differences, compared to WO4L successful case (figure 7(c)), while femurLength has little changes during unsuccessful cases of W4L (figure 7(b)) and WO4L (figure 7(d)). While combining all sub-figures for analysis, we found that the femurLength parameter played an essential role in leading tetrapods to evolve into bipedalism. We believe that the initial shortening of femur length on tetrapods enabled bipedal stability, but this same shortening had little impact on the stable structure of the hexapod. As a result, right after the shift bipeds gradually reached walking speed superiority while maintaining their steady gait, and finally replaced hexapods.

We found a special case in which there was only a 10% probability of reaching bipedal creatures within the WO4L condition. In figure 9, we can see that bipeds first dominate the fitness, and then dominate the population in the number of individuals. To better understand how this superiority from bipeds is established we analyze the parameters of the virtual model in that same transition zone. We find that there are apparent fluctuations of the tibiaLength parameter, which have close connection with the change of fitness between generation 1003 (point A labeled in figure 10) and 1007 (point B labeled in figure 10). We do not find any other parametric change at the same interval with the same magnitude as this one observed in the tibial length.

Authors of [28] presented the role of the tibia length, in particular, in obtaining the optimal walking speed. Ueno et al [29] proposed that the leg bone length, especially of the tibia, is a potential morphological factor for achieving superior running performance in well-trained endurance runners. These all agree that tibia length has a high impact on obtaining the forward speed. In the work of [30], they found that modern humans have exceptionally large lumbo-sacral joints and hind limbs, which carry a greater average load than those of tetrapods or suspension animals of the same body size. These characteristics represent a major change in later hominid evolutions, but trends in hind limbs cannot indicate changes in the femur or the tibia, so it is hard to say if the femur decreases or the tibia increases in tetrapods evolved to bipeds in nature, which is still a controversial topic.

Animals in nature are vastly different, with varied geometries, arrangements, gaits and leg architectures. As it would be impossible to re-create an experiment which fully captured all constructions, from sprawled spider legs to four-linked ungulates and the complex synergistic nature of mono and biarticular muscles on legs, this work relied on a simplification of leg structure and controller. This simplification allowed us to understand how parameters (morphological and control), which were initially optimized to a stable structure (higher number of legs), could migrate to a less stable structure (two or four legs). Our choice for controller enabled an evolutionary interdependence of limbs, regardless of the number of legs, while also allowing enough locomotory independence (very often resulting in gait failure, as walking is still a difficult problem). Although a two-linked leg morphology (femur and tibia) suffers from the lack of leg elasticity which could be provided by a third joint (metatarsal joint, widely used by ostriches, kangaroos and humans to amplify their motion), we concentrated our work on static locomotion to reduce the energy and deformation required during stance phase.

In this work we adopted locomotion speed (distance over time) as the evolutionary criterion. We understand that other animals may have evolved their locomotion based on energy efficiency, or even solely on gait stability, but we believe that optimizing for speed will also, implicitly, account for stability as long as the total simulation time forces the creature to take many steps. Over the course of 5 s our bipedal creature takes more than 10 steps on a fully static gait, which points to the existence of other 'faster but unstable' gaits not enduring the full simulation time and collapsing at a shorter distance (from our observations those fade within two steps, sometimes with the creature jumping forward and landing on the ground).

An optimization based on efficiency is a possibility, as a lower energy expenditure would translate as survival in inhospitable environments. In juxtaposition, an optimization based on speed is a better representation of animals which need to escape from predators or to consume other animals. Our choice of optimization criterion solely took speed/stability into consideration, limiting our results to the development of a fast bipedal walking gait over more than ten steps, with little regards for energy efficiency.

The ratio between these two links will vary significantly between gait type and number of legs (sloths have longer tibias than femurs, but a highly arboreal life). As our work focused on a static walking the results observed here might not fit all animals: the human femur is notably longer than the tibia [31] (and the humerus longer than the ulna), but in this paper we focused on static walking and did not do a longitudinal comparison between the femur and the tibia. As the parameters analyzed are from the locomotion transition zone those will not be optimal parameters, just the ones leading to the locomotion shift. Here, we aim to compare the leg length among hexapods, tetrapods, and bipeds in their transitions.