I Introduction
Many important applications require robots to traverse cluttered obstacles such as search and rescue in rubble, environmental monitoring on the forest floor, and extraterrestrial exploration through Martian rocks. The dominant approach of robotic locomotion in complex environments is to avoid obstacles [4, 3, 17], which could be challenged in such cluttered terrain. By contrast, insects such as the discoid cockroach are excellent at traversing cluttered large obstacles, often using their bodies and appendages to make physical interaction with obstacles [11, 20, 10, 6, 7, 9] (for a review, see [16]). Understanding how the physical interaction between a selfpropelled locomotor and obstacles leads to or hinders traversal can help inform robot locomotion beyond avoiding obstacles and increase the terrain that they can access.
Recent animal and robot studies from our lab established a potential energy landscape approach for quantifying locomotorobstacle interaction, which provided insights into the mechanical principles of how different modes of locomotion emerges from the interaction, which can be controlled to achieve traversal [16, 15]. However, there still remains a lack of understanding of the dynamics of cluttered obstacle traversal, which is not modeled by the potential energy landscape approach. Here, we take the next step towards this by developing a dynamic model of locomotorobstacle interaction using a minimalistic model system and studying how the properties of locomotor selfpropulsion and obstacle affect traversal.
Our study is inspired by and builds upon findings from our previous work [15]. The potential energy landscape approach revealed that, when interacting with large obstacles, a selfpropelled locomotor has highly dissipated dynamics [16] and is strongly attracted to basins of a potential energy landscape resulting from locomotorterrain interaction [16, 15]. To escape from attraction to certain basins that lead to unfavorable modes (e.g., become trapped against the obstacles [9, 15]), the locomotor can use both active sensing and control and passive kinetic energy fluctuation from selfpropelled oscillation [15]. As a first step towards understanding traversal dynamics, we focus on stochastic dynamics from selfpropelled oscillation and do not consider active sensing and control. In addition to adding dynamics, our work further explores the interplay of selfpropulsion and random oscillation and the effect of asymmetry in terrain physics common in nature [12, 19] and how locomotorobstacle interaction leads to traversal over large spatiotemporal scales, which the previous work did not address.
Our minimalistic model system consists of a circular locomotor body selfpropelled forward with only two translational degrees of freedom moving on a level plane interacting with a pair of grasslike beam obstacles that rotate within the plane. A lateral random force is added to introduce body oscillation, leading to stochastic dynamics. We choose to use a random force considering that animals often have random motion during locomotion
[5, 2, 21, 13]. In addition, we introduce asymmetry in the environment by allowing the two beam obstacles to have different torsional stiffness. Moreover, we used simulation results with systematically varied propulsive and random forces and asymmetric beam stiffness to test how well the potential energy landscape can qualitatively predict traversal dynamics without solving equations of motion. Furthermore, we extended the single locomotorobstacle interaction model to a larger terrain with multiple obstacles in a lattice configuration test how well it can statistically predict how the body moves across a complex terrain on a larger scale.We introduce how the model was defined and the dynamics modeled (Sec. IIA, IIB), how the potential energy landscape was calculated (Sec. IIC), and how we generated a multiobstacle field (IID). We then describe results using the stochastic dynamics simulation to study traversal of a pair of beams and testing how well the landscape model informed traversal (Sec. III) and how well the traversal of a multiobstacle field can be predicted by a Markov chain Monte Carlo method (Sec. IIIC). Finally, we summarize our findings and their implications (Sec. IV).
Ii Methods
Iia Model Definition
Our simplistic 2D model of bodyobstacle interaction (Fig.1) consists of three rigid bodies on a horizontal plane: a circular disc representing the locomotor body (mass = 1 kg, radius = 10 m), and two plates with viscoelastic revolute joints at their far ends (mass = 0.1 kg, length
= 25 m, moment of inertia about joint
= 20.83 kgm). Hereafter, we refer to them as beams. The beams align with each other when unloaded, forming a closed “gate”. For simplicity, we assumed that the body surface is frictionless. Thus, bodybeam contact results only in forces normal to the body surface, which generates no torque. Thus, body orientation remains unchanged. In addition, we assumed that 2dimensional bodybeam collisions are inelastic with partial energy loss. We used coefficient of restitution to measure the degree of elasticity, defined as the ratio of final to initial relative normal velocities between the body and beams during each collision. Here we tuned by observing different trials and chose a value of 0.8 to conserve a large portion of mechanical energy so that traversal still occurs even after a large number of collisions. Finally, for simplicity, we assumed that there is no friction against the ground (although it can be added–see Sec. IIIC, which should only quantitatively but not qualitatively change our results).Besides forward propulsion, animals with a sprawled leg posture like cockroaches can generate substantial lateral forces during locomotion [18]. In addition, animals can randomly change their movement direction during locomotion [13]. Considering these and for simplicity, selfpropulsion of the locomotor body was modeled by a constant forward propulsive force and an oscillating random lateral force (Fig. 1
red and green arrows) with a standard Gaussian distribution:
(1) 
where is a constant characterizing the magnitudes of the random forces and
, a normal distribution with a mean of 0 and a standard deviation of 1.
Because the body does not rotate, forward and lateral forces are always in the  and  directions, respectively. Applying Newton’s second law in these two directions, the equations of motion of the body are:
(2) 
where is body mass, is body acceleration, and is the resistive force from the beam, with for the left and right beams, respectively.
Applying Newton’s second law to the beam:
(3) 
where is beam moment of inertia about its joint, are the orientation, angular velocity, and angular acceleration of the beam, and and are the torsional stiffness and damping coefficient of the beam joint.
IiB Interaction Dynamics
We used the Euler method to numerically integrate forward in time to calculate the dynamics of the system. A time step of 0.004 s was chosen to achieve good numerical accuracy while maintaining computation efficiency. The runtime of each trial varied from 5 to 20 min on a lab workstation. The key part of dynamics is to determine the interaction between the body and beams. We developed two complementary methods to model interaction dynamics for two different scenarios: a collision method and a constraint method. The collision method models repeated collisions between body and beams after initial contact, which occur due to the competing forward propulsive force and backward beam resistive forces. However, as these collisions dissipate energy, the duration of single collision decreases and eventually becomes smaller than the time step, leading to a significant increase in numerical error [1]. In this case, we assumed that the interacting body and beam have no relative motion in the direction normal to the contact point but can move relatively in the tangential direction and developed the constraint method to describe interaction dynamics.
IiB1 Collision method
When calculating collision occurring over an infinitesimal time, we considered the body and beam as a system with the center of rotation at the beam joint. The system is subject to a finite external beam joint torque and external forces . After multiplying the finite values with an infinitesimal time, the external momentum impulse is negligible. In addition, there is a reaction force acting on the fixed beam joint, which has a zero moment arm. Therefore, the angular momentum of the system about the instantaneous center of rotation is conserved during a collision.
The calculation of collisions depends on the geometric relationship between body and beam. During traversal, there are two different contact cases: tangential contact and point contact.
In the tangential contact case (Fig. 2), the beam is tangential to the body surface at the point of contact. The normal direction of the body surface at the point of contact is perpendicular to the beam. Applying conservation of angular momentum to the bodybeam system:
(4) 
In addition, applying the definition of CoR to the tangential contact case:
(5) 
where is the distance from beam joint to contact point, and are beam angular velocities before and after collision, and are local velocities of the beam at the contact point before and after collision, and are body velocities before and after collision, , , , and are their projections in the direction normal and tangential to the beam. Note that = because there is no interaction force in the tangential direction from our frictionless body assumption. Using Eqns. 4, 5, we can solve the dynamics of the collision in the tangential contact case.
In the point contact case (Fig. 3), the free end of the beam contacts the body surface. Because the normal direction of the body surface at the contact point is no longer perpendicular to the beam, it is more complex to formulate a expression in scalar form. Considering that its collision is small and only happens in the short duration before detachment, we assumed the collision to be perfectly elastic in this case for simplicity. Thus, mechanical energy is assumed to be conserved at the collision time step:
(6) 
Also the system complies with the conservation of angular momentum before and after collision:
(7) 
(8) 
(9) 
where is the change of body momentum, and are body velocities before and after the collision, is the change of angular momentum of the beam before and after the collision, is the interaction force acting on the body from the beam, is the distance from the beam joint to contact point, and is the change of angular momentum of the system, which is zero at collision. The velocities after the collision can be solved from those before using Eqns. 69.
IiB2 Constraint method
The constraint method is based on Newton’s second law and the geometric constraints that the body and beam have the same normal velocity at the contact point. It also differs between the two contact cases.
In the tangential contact case, the velocity constraint is:
(10) 
Taking time derivative of both sides of Eqn. 10, we have:
(11) 
where is the center of mass position measured from the beam joint, is body velocity, and is the body acceleration.
In the point contact case, the velocity constraint is:
(12) 
Taking time derivative of both sides of Eqn. 12, we have:
(13) 
where is beam length.
In the simulation, the collision and constraint methods are integrated following the workflow shown in Fig. 4. The simulation first runs the collision method in loops. As the change of the body’s momentum during each collision decreases below a small threshold kgm/s, i.e., collisions become small enough and exchange little momentum, the simulation switches to the constraint method in loops. During the switch, beam angular velocities are updated to satisfy the initial velocity constraint of the constraint method.
IiC Potential Energy Landscape
The potential energy landscape is a modeling approach to model the conservative forces during locomotorterrain interaction over relevant degrees of freedom. Without knowledge of nonconservative forces that are often difficult to measure, it provides a useful approach for understanding how the system may or may not move on the landscape, as long as potential energy landscape plays a major role in shaping dynamics [16]. In our system, potential energy arises from elastic torsional joints of the beams:
(14) 
where is the total beam potential energy, is the elastic stiffness of the L/R beam, and is the deflection angle of the beam, with for the left and right beams, respectively.
The landscape is calculated as a 3D surface over a mesh grid of m. At each grid point, we assumed that beams are always deflected forward and contact the body if the body is within the range of beam radius. During traversal, the body pushes forward and deflects the beams, resulting in two potential energy barriers on each side (Fig. 5AB). The two barriers overlap in the middle area, where the body can interact with both beams. Fig. 5C is an example of asymmetric landscape with different and .
Because we assumed a constant forward , the potential energy landscape can also include conservative potential energy from . The summed potential energy is:
(15) 
where is the potential energy from . The zero level of was defined at the anterior boundary = 60 m.
A basin that spans across the direction in front of the beams (around = 20 m) emerges after including the energy of of a modest magnitude in the potential energy landscape (Fig. 5D). This helps us understand how the landscape is tilted by the propulsive force. Within insufficient , the locomotor body would be trapped in this basin and not traverse. Adding kinetic energy fluctuation induced by the random force can help the body escape the horizontal basin. With sufficient , the landscape is so heavily tilted that the basin disappears, and the locomotor should always traverse.
The evolution of the potential energy of the system in the simulation traversal trials can be plotted over the landscape to visualize how the bodybeam interaction influences body motion. Typically, the trajectory of the system potential energy is above the landscape surface due to the beam inertia. However, sometimes the trajectory penetrates the landscape surface. An example is shown in Fig. 6B. In this case, the left beam lost contact with the body and returned its initial orientation (Fig. 6
A). Because the body is still within the radius of the beam, our landscape calculation that always assumes beams being deflected forward overestimates the potential energy. If we only consider the landscape from the right beam that is deflected, this artifact is removed (Fig.
6C). See multimedia material part 3 for an example video.IiD MultiObstacle Field Interaction
Using the bodybeam interaction model above with a single pair of beams (a gate), we can create a large cluttered obstacle field by composing multiple pairs of beams (gates) in a lattice arrangement. Fig. 7 shows a 9 9 gate obstacle field used for further simulation. Except for the first gate region which uses a set of initial system states (), the system state input of each gate region comes from the output of the previous gate region. An index was used to track which gate region is being activated. For this 9 9 gait obstacle field, and . Finally, the main program outputs a trajectory restricted in a single gate region, which is a series of local positions , where m, m. Later in the visualization process, according to the gate index , are converted to the global position data , where m, m. See multimedia material part 4 for a simulation example.
The pseudo code:
All the modeling calculations and simulations were performed in MATLAB R2019b.
Iii Results
Iiia Traversing symmetric beams
We first studied the probability of traversing a pair of symmetric beam obstacles under selfpropulsive and random forces. The body starts at an initial forward velocity , accelerates for a short distance before making initial contact with beams. The criterion of traversing successfully is that the body can reach = 65 m, where the body is out of the range of beams (multimedia material, part 1). Failure resulted from: (1) the body gets stuck by the beams (multimedia material, part 2), (2) the body deviates from the desired track and exits from the side boundaries ( 25 m or 25 m), and (3) the body did not reach = 65 m within maximal iterations.
From the potential energy landscape (Fig. 5A), the locomotor body needs to have sufficient kinetic energy to overcome a potential energy barrier to traverse [15]. Because higher propulsive force can tilt the landscape, facilitating to overcome barriers by lowering them, while the random force adds more kinetic energy fluctuation in the system. Thus, we hypothesized that the probability of traversing a pair of symmetric beam obstacles increases with propulsive force and random force.
To test this hypothesis, we varied propulsive force and the magnitude of random force using system parameters in Table I and ran 100 trials for each given () to obtain traversal probability (Fig. 8). For a given , traversal probability increased with propulsive force. However, for a given , traversal probability did not always increase with random force. A random force of up to 20 N increased traversal probability. However, as increased beyond 20 N, traversal probability decreased. Examination of simulation videos revealed that this was because more trials exited the side boundaries of the gate region before arriving at the beams. Thus, a sufficient but not excessive lateral random force helps traverse symmetric beam obstacles.
Parameter  Value 
Propulsive force  4,5,6,7,8,9 N 
Random force magnitude  0,10,20,30,40 N 
Beam stiffness  400 Nm/rad 
Beam damping  50 Ns/rad 
Initial forward velocity  0 m/s 
Coefficient of Restitution  0.8 
Oscillator frequency  50 Hz 
Acceleration distance  20 m 
Max iterations  3000 
IiiB Traversing asymmetric beams
Next, we studied the probability of traversing a pair of asymmetric beams with different stiffness. The potential energy landscape becomes asymmetric in this case, with a barrier lower on one side where the beam has a smaller stiffness (Fig. 5C). We hypothesized that the probability of traversing on the side of lower stiffness increases with the level of asymmetry.
To test this hypothesis, we varied the left beam stiffness while keeping the right beam stiffness constant and varied and systematically (Table II). Unlisted parameters are the same as those in Table I. Here we selected three groups with varying stiffness ratios as examples to illustrate the typical results (Fig. 9). Only successful traversal trajectories are plotted over the landscape. In the symmetric control group with both beam stiffness being 500 Nm/rad (Fig. 9A), traversal was rare with similar probability on both sides (L:13%, R:8%). This is because the propulsive force was barely sufficient to overcome the potential energy barriers. When the left beam stiffness was reduced from 500 to 300 Nm/rad, traversal was frequent (52%), with the body pushing across the left beam (Fig. 9B). As the left beam stiffness further reduced to 100 Nm/rad, traversal was almost guaranteed by pushing across the left beam (98%, Fig. 9C).
Parameter  Value 
Propulsive force  7,8N 
Random force magnitude  10,20,30,40,50,60 N 
Left beam stiffness  100:50:500 Nm/rad 
Right beam stiffness  500 Nm/rad 
For all the asymmetric tests varying and
, we classified the trajectories into “left” or “right” types, reflecting traversal tending to which barrier side, by comparing the maximal deflected angle of the left or right beam. Then, we can quantitatively evaluate the asymmetry of trajectories by calculating the ratio of the number of trials for each type to the total number of successful traversal trials (Fig.
10). For N, as the left beam stiffness decreases, traversal becomes less frequent on the right side. However, for N, this trend is still existed, but less monotonically due to the high stochasticity. In addition, traversal probability increased with (Fig. 10A vs. B). These results well supported our hypothesis and were consistent with observations in the traversal of symmetric beams.IiiC MultiObstacle Traversal
On the larger terrain, trajectories become more complicated and lack of patterns. Here we focused on how the locomotor body transits between adjacent gate regions. As a first step, we assumed that all the gates are symmetric with the same stiffness (400 Nm/s). Each gate is considered as a single system with an input and an output when crossing region boundaries. Here we discretized the boundary of a gate region into 6 segments: middle traverse (MT), left traverse (LT), right traverse (RT), left deflect (LD), right deflect (RD) and enter (EN) (Fig. 11). For the traversal problem, EN is only for input and MT is only for output, whereas the other four can be either input or output. Another possible output is that the body is trapped inside a gate region in between the two beams. Motion in each gate can be treated individually as a stochastic process independent of the history before entering. And we assume for fixed model parameters, there is a constant probability correspondence between input and output states. Thus, the Markov Chain [14] can be applied to describe the stochastic transition between gate regions.
For a discretetime Markov process, the continuous states needed to be discretized into a finite number of input and output states. Here, the state of the body when crossing the boundary is , where {EN, RD, LD, RT, LT, MT} is the boundary index, is the position on the corresponding the boundary. is measured as the distance to several benchmark points (EN: m, LT/LD: m, RT/RD: m, MT: m). is the velocity when crossing boundaries. To set up the state space for the Markov chain, we discretized the input states for boundaries EN, RD and LD (Table III). Inputs from RT or LT rarely occur (fewer than 4% of all trials) and they are not involved with the bodybeam interaction. Thus, for simplicity, we did not further discrete RT and LT like the other boundaries to improve computational efficiency. In total, there are 87 input states, , and 88 output states, , with being the state of being trapped in a gate.
EN  RD  LD  
(m) 

[10,20],[20,30]  [10,20],[20,30]  
(m/s) 




(m/s)  [10,20] 



Groups  25  30  30 
For Markov Chain analysis, we need to obtain the transition matrix , where is the probability from input to the output during a gate traversal (Eqn. 16). For our discretization, is an matrix. The matrix is obtained by running 100 trials per group with random force in the single gate simulation.
(16) 
where is the probability of output state after the traversal, is probability of input state before the traversal.
With the transition matrix, the Markov chain Monte Carlo (MCMC) method [8]
can be used to predict the consecutive transitions between gate units. Monte Carlo is a statistical simulation method relying on repeated random sampling from a probability distribution. MCMC is the combination of Markov chain and Monte Carlo. In each gate, input is known, either from the initialization or from the former gate. The probability of outputs after traversing this gate is given by multiplying the Markov chain transition matrix. We randomly pick a certain output state based on the probability distribution
and use it as the input of the next gate. By repeating the Markov chain and Monte Carlo in turn, we can obtain a path across the terrain. The simulation ends after a maximal number of iterations. The pseudo code:Unlike our dynamic simulation, the Markov chain Monte Carlo algorithm cannot precisely predict the trajectory across the terrain. But we can compare the prediction for the finial lattices where the trials finally stop. By repeating MCMC trials, we can obtain a predicted distribution of the final lattice. To compare two methods, we have the same parameters of gates in dynamic simulation as the ones in MCMC. Ideally, the MCMC simulation can predict the path on an infinitely large terrain. However, if the body keeps accelerating after traversing a series of gates, the locomotor velocity will keep increasing towards infinity, and it would be impractical to discretize velocity for the Markov chain analysis. To address this issue, we increased the beam damping ( Ns/rad) and added a viscous force ( Ns/m) between body and the horizontal plane so that the velocity can remain in the tested range (Table. III). The body started at the origin of the entire terrain ( = 0 m, = 0 m) with = 0 m/s and = 15 m/s. For our gate lattice, each MCMC trial is allowed to have a maximum of 13 steps to reach the furthest gate. But if the body exits from the terrain or becomes trapped inside a gate (i.e., reaching ), the trial will stop at that gate.
To test how well the MCMC simulation predicted traversal outcome, we ran 100 trials using both dynamic and MCMC simulation and compared the distribution of the final location of the body within the obstacle field. The patterns are strikingly similar (Fig. 12). The correlation coefficient, , between the two simulations was 0.914. In addition, the root mean squared error, , was only 1.018, i.e., the error is around 1 out of 100 trials. These results demonstrated that the MCMC simulation can predict our dynamic simulation traversal result very well. Considering that the MCMC method can finish the 100trial prediction within 10 seconds whereas the dynamics simulation took around 100 hours, there is a huge advantage in saving computational time of the MCMC method.
To evaluate the generality of this result, we tested three other sets of randomly chosen initial states and found excellent match throughout (Table IV).
(m)  (m/s)  (m/s)  
0  0  15  0.914  1.018 
23  11  16  0.9097  1.133 
4  5  10  0.8929  1.252 
17  13  12  0.9026  0.9813 
Iv Discussion
In this study, we developed a stochastic dynamic model to simulate a selfpropelled, minimalistic locomotor body with random forces traversing beam obstacles in two dimensions on a horizontal plane. We found that a larger forward propulsive force and a sufficient but not excessive lateral random force facilitate traversal, in accord with the kinetic energy needed to overcome potential energy barriers. Asymmetry because of different horizontal beams stiffness can cause asymmetric distribution in trajectories, in accord with the higher probability to traverse on the lower barrier side. The model of interacting with a basic unit of obstacle (two beams) was further extended to simulate the traversal of a multiobstacle field. Finally, we applied the Markov chain Monte Carlo method to predict how the body statistically transitions between adjacent gates, which predicted stochastic paths on the larger terrain well.
Our results demonstrated that the potential energy landscape is useful for understanding traversal dynamics without solving them exactly. This would be particularly useful for modeling systems with complex interaction, which is common in animal and robot locomotion in complex terrain. In addition, we can redesign the landscape barriers and basins by adjusting the model parameters related to physical interaction, which can be applied in robotics challenges such as motion control and path planning. Furthermore, stochasticity, which is common in locomotion and usually considered a nuisance, can be beneficial for robots when they are trapped by obstacles and need to overcome potential energy barriers. Finally, the exploration of multiobstacle field traversal supported the plausibility to decompose natural terrain as a combination of multiple obstacles [16]. Our results suggested that, by statistically understanding and predicting the interaction with each type of obstacle unit, larger scale traversal processes may be rapidly predicted.
Acknowledgment
The authors would like to thank Noah Cowan, Ratan Othayoth, and Yaqing Wang for discussion and three anonymous reviewers for suggestions.
References
 [1] (2005) Sufficient conditions for the existence of Zeno behavior. In Proceedings of the 44th IEEE Conference on Decision and Control, pp. 696–701. Cited by: §IIB.
 [2] (2000) Caribou movement as a correlated random walk. Oecologia 123 (3), pp. 364–374. Cited by: §I.

[3]
(1991)
The vector field histogramfast obstacle avoidance for mobile robots
. IEEE Transactions on Robotics and Automation 7 (3), pp. 278–288. Cited by: §I.  [4] (2002) Vision for mobile robot navigation: a survey. IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (2), pp. 237–267. Cited by: §I.
 [5] (2008) Noise in the nervous system. Nature Reviews Neuroscience 9 (4), pp. 292–303. Cited by: §I.
 [6] (2018) Bodyterrain interaction affects large bump traversal of insects and legged robots. Bioinspiration & biomimetics 13 (2), pp. 026005. Cited by: §I.
 [7] (2018) Dynamic traversal of large gaps by insects and legged robots reveals a template. Bioinspiration & biomimetics 13 (2), pp. 026006. Cited by: §I.
 [8] (1995) Markov chain monte carlo in practice. CRC press. Cited by: §IIIC.
 [9] (2021) Shapeinduced obstacle attraction and repulsion during dynamic locomotion. The International Journal of Robotics Research 40 (67), pp. 939–955. Cited by: §I, §I.
 [10] (2009) Characterization of obstacle negotiation behaviors in the cockroach, Blaberus discoidalis. Journal of Experimental Biology 212 (10), pp. 1463–1476. Cited by: §I.
 [11] (2015) Terradynamically streamlined shapes in animals and robots enhance traversability through densely cluttered terrain. Bioinspiration & Biomimetics 10 (4), pp. 046003. Cited by: §I.
 [12] (2015) Analysis and design of asymmetric oscillation for caterpillarlike locomotion. Journal of Bionic Engineering 12 (2), pp. 190–203. Cited by: §I.
 [13] (2017) Unpredictability of escape trajectory explains predator evasion ability and microhabitat preference of desert rodents. Nature Communications 8 (1), pp. 1–9. Cited by: §I, §IIA.
 [14] (1998) Markov chains. Cambridge university press. Cited by: §IIIC.
 [15] (2020) An energy landscape approach to locomotor transitions in complex 3d terrain. Proceedings of the National Academy of Sciences 117 (26), pp. 14987–14995. Cited by: §I, §I, §IIIA.
 [16] (2021) Locomotor transitions in the potential energy landscapedominated regime. Proceedings of the Royal Society B: Biological Sciences 288 (1949). Cited by: §I, §I, §I, §IIC, §IV.
 [17] (1990) Exact robot navigation using artificial potential functions. Ph.D. Thesis, Yale University. Cited by: §I.
 [18] (2000) Mechanical models for insect locomotion: dynamics and stability in the horizontal plane I. Theory. Biological Cybernetics 83 (6), pp. 501–515. Cited by: §IIA.
 [19] (2017) Changes in mechanical work during neural adaptation to asymmetric locomotion. Journal of Experimental Biology 220 (16), pp. 2993–3000. Cited by: §I.
 [20] (2020) Coordinated appendages accumulate more energy to selfright on the ground. IEEE Robotics and Automation Letters 5 (4), pp. 6137–6144. Cited by: §I.
 [21] (2020) Randomness in appendage coordination facilitates strenuous ground selfrighting. Bioinspiration & Biomimetics 15 (6), pp. 065004. Cited by: §I.