Trunk Pitch Oscillations for Energy Trade-offs in Bipedal Running Birds and Robots

Bipedal animals have diverse morphologies and advanced locomotion abilities. Terrestrial birds, in particular, display agile, efficient, and robust running motion, in which they exploit the interplay between the body segment masses and moment of inertias. On the other hand, most legged robots are not able to generate such versatile and energy-efficient motion and often disregard trunk movements as a means to enhance their locomotion capabilities. Recent research investigated how trunk motions affect the gait characteristics of humans, but there is a lack of analysis across different bipedal morphologies. To address this issue, we analyze avian running based on a spring-loaded inverted pendulum model with a pronograde (horizontal) trunk. We use a virtual point based control scheme and modify the alignment of the ground reaction forces to assess how our control strategy influences the trunk pitch oscillations and energetics of the locomotion. We derive three potential key strategies to leverage trunk pitch motions that minimize either the energy fluctuations of the center of mass or the work performed by the hip and leg. We suggest how these strategies could be used in legged robotics.


Introduction
Creating dynamic running motion for bipedal mechanisms is challenging due to the complexity in controlling an underactuated trunk-leg structure, that has nonlinear coupled dynamics and intermittent ground contacts [1]. Robotics research often focuses on the control of lower extremities for humanoids and suppresses trunk motions for simplicity [2,3]. In contrast, bipedal animals have diverse morphologies and display a wide range of motion patterns with prominent trunk movements [4]. In particular, terrestrial birds are able to generate exceptionally agile, energy efficient, and robust motion; irrespective of their vast variability in body size, posture, and habitat [4][5][6]. Birds exploit their trunk's inertia and generate trunk movements to assist the postural stability [7,8]. The concept of leveraging trunk pitch oscillations has been analyzed for humanoids, where the trunk motion assists energetics by redistributing the work between the leg and hip joint [9]. Here, we investigate whether a similar strategy exists for terrestrial birds with a pronograde (horizontal) trunk orientation. We use a spring-loaded inverted pendulum model with a controller based on a virtual point (VP) concept. Our aim is to analyze how the magnitude and direction of the oscillations depend on the VP position and running speed.
Birds show exceptional locomotor efficiency and capabilities, yet the mechanisms underlying such agile motion are not well understood [6]. The majority of the prior research builds on simplified models that address the motion of the center of mass (CoM) without the trunk [6,10,11]. These studies are able to recreate some of the essential characteristics of avian gaits. In particular, these gaits display asymmetric traits such as left-skewed ground reaction forces (GRF) [5,10,12], distinct leg lengths at leg touch-down/takeoff [10,12], and asymmetric leg protraction/retraction angles [13][14][15]. Recent studies reveal the role of the pronograde (horizontal) trunk orientation in these gait asymmetries [16][17][18][19]. Terrestrial birds show pronograde trunk posture, where the CoM is located cranial (headward) to the hip joint with an inclination of 100°-135°with respect to the vertical axis [4]. Bird's unique hip attachment necessitates continuous hip extension torques to hold the trunk in a horizontal orientation against gravity [17]. Furthermore, the trunk amounts to 70-80 % of the total body mass and has a high moment of inertia [4,15,20]. The heavy, horizontal trunk requires high hip torques, which can be realized only if the inertia is high enough to resist the trunk rotation. Hence, the trunk is an integral part of understanding bipedal locomotion.
The spring-loaded inverted pendulum (SLIP) model can predict fundamental characteristics of running motion for animals with varying morphologies. Its simplicity and its prediction capabilities make the model well suited to investigate both humans and avians [24]. The model consists of a point-mass body attached to a massless and springy leg. The SLIP can be extended with a rigid trunk (TSLIP), which is actuated by a torque at the hip [17]. One method to select the hip torque is based on the VP concept, where the ground reaction forces are redirected to intersect at a virtual point. The GRF intersecting at a point bounds the moment GRF creates around the CoM. Therefore, some research considers the VP as a pivot point that assists postural stability, and implement the VP as a control target in both avian and human TSLIP models [1,17,25].
The VP can also be regarded as a method to generate trunk oscillations with different magnitudes and directions, which is shown in [9] for the human morphology. The relative position of the VP with respect to the axis passing through the CoM and foot determines the direction of the trunk rotation. The trunk rotates backward during the stance phase if the the VP is above (VP A ) the CoM-foot axis, whereas it rotates forward if VP is below (VP B ) [9]. The distance between the VP and CoM determines the magnitude of the trunk oscillations. As a consequence of the direct mapping between VP location and trunk oscillations, there is no need to adjust the control gains when the operating point of the gait changes. Hence, the framework allows creating a variety of trunk oscillations across a range of speeds in a systematic and comparable manner.
Another advantage of the VP framework is that it generates feasible hip torque profiles. The VP concept creates a coupling between the magnitudes of the hip torque and leg force. Consequently, the hip torque profile inherits the continuity and smoothness properties of the leg force profile. It makes the VP framework well suited for investigating the role of the hip torque in generating trunk motions.
The prospect of using the VP framework for the pronograde trunk is motivated by the avian gait measurements. A VP A is observed in the running gait of chickens [25], and in walking and running gait of quails [17,19]. Inspired by these observations, the VP is implemented for the pronograde TSLIP model [17]. What is missing is the relation of the VP to the trunk oscillations for the avian morphology. It is shown that a VP A in the TSLIP framework produces backward trunk motion, which is the opposite of what is observed in human running [9]. Instead, a VP B provides the matching forwards trunk rotation in the human model. Here, we inquire what relation exists for the avian morphology.
Our overall goal is to conceptualize how trunk motion is beneficial for bipeds with different trunk orientations. We would like to transfer the gained knowledge to controller and robot design.Currently, robot studies often maintain the trunk posture at a fixed angle without any movements. There are only a few studies that address dynamic trunk stability; one of which uses a VP control scheme [26] for walking gaits, while the other tracks a sinusoidal reference pitch angle [27] within the TSLIP framework for the ATRIAS robot. However, these studies involve a single parameter set and do not extend to multiple operating points with multiple speeds and different trunk oscillation characteristics.
In this paper, we present a unified framework to analyze the effect of trunk motions for bipedal running, with a focus on avians. We implement a TSLIP model with varying VP targets for a pronograde trunk and systematically compare the resulting gaits for speeds of 4-10 m s −1 . Our work examines four hypotheses, where we question if the pronograde TSLIP model is able to Hypothesis 1. predict left-skewed GRF profiles, whose magnitude is proportional to the forward speed and matches to the avian gait data in [5,13,28]. Hypothesis 2. generate downward (dorsoventral) trunk pitch motion at stance phase, whose magnitude is proportional to the forward speed [9]. Hypothesis 3. utilize the VP to alter the gait dynamics, which gives rise to multiple solutions with different gait characteristics for a given speed. Hypothesis 4. determine the VP location in favor of energetics, in a similar manner as [9]. This work presents three novel contributions that are essential for the pronograde TSLIP model.

Contribution 1.
We show that a bilinear leg damper [29] generates zero damping force at leg touch-down, and therefore is better suited for the pronograde TSLIP model compared to a typical linear leg damper [17,30]. Contribution 2. We investigate which frame is suitable for defining the VP location to obtain feasible gaits. The VP is typically defined with respect to the body frame and above the CoM in the TSLIP model to facilitate a self-stabilizing postural behavior [25]. There is little evidence in literature to support the choice of the VP frame. It is shown that a body or world aligned VP predicts the GRF better than a trunk aligned VP in human walking [31]. We extend this analysis to avian running using the TSLIP model. Contribution 3. The change in VP position corresponds to modifying the orientation of the GRF vector [9]. A VP A creates more vertically oriented GRF compared to VP B . A more vertically oriented GRF vector has a smaller horizontal GRF component, which accelerates and decelerates the body less (i.e., it is energetically more favorable) [32,33]. In a similar manner, we test whether a VP A in our avian TSLIP model can yield such an energetic benefit. In particular, we establish a relation between the VP location, GRF alignment and energy variations of the CoM.

Simulation Model
In this section, we describe the TSLIP model applied in this work. It consists of a trunk with mass m and moment of inertia J, which is connected to a massless leg of length l that has a parallel springdamper mechanism (see Figure 1a). The dynamics are described by a flight phase, where the CoM moves in a ballistic motion and a stance phase, where the leg force and hip torque propel the body forward. The switch between these phases occurs at the touch-down (TD) and take-off (TO) points, where the foot establishes contact with the ground and the leg extends to its rest length l 0 , respectively. Swing leg dynamics are omitted, similar to other TSLIP studies in [17,34].
The equations of motion for the CoM state (x C , y C , θ C ) during the stance phase can be written as, The linear leg spring force F sp =k (l−l 0 ) and bilinear leg damping force F dp =cl (l−l 0 ) generate the axial component of the GRF in foot frame T . The hip torque τ H generates the tangential component of the GRF T , as seen in Figure 1d.
In our formulation, k denotes the stiffness of the leg spring and c is the damping coefficient.
The leg is passively compliant, where the damper removes energy from the system by performing negative work, and the spring stores and releases elastic energy in sequence. The hip is actuated and produces net positive work to balance the energy depleted by the leg. Thus, the actuated hip torque τ H is the only element that we can actively control to induce trunk pitch oscillations. We select τ H , such that the GRF points to a VP, which is characterized by the radius r V P (i.e., distance between the hip and CoM) and angle θ V P , as shown in Figure 1d ( , ). The hip torque as a function of the VP is expressed as, The sign of the hip torque depends on whether the VP is placed above or below the the hip-foot axis (i.e., leg axis). For the pronograde trunk, all VP A within our parameter range remain above the leg axis across entire stance phase. Thus, the resulting hip torque is always negative. VP B above the leg axis yield similar torque pattern as VP A . In contrast, VP B below the leg axis yield positive hip torque at the first half of the stance phase and negative towards the second half. These region is referred to as VP BL , and is illustrated in Figure 1c.
The concept of VP control is open-loop, therefore the VP controller is highly sensitive to changes in the initial state and model parameters. This parameter sensitivity makes it challenging to find feasible gaits. In order to simplify and guide the parameter search, we use the iterative gait generation framework in [9]. The framework has an initial controller that combines the VP based torque in Equation (2) with a PID controller on the pitch angle, which yields stable gaits with semifocused GRF in Figure 5b and 5c. The solutions are fed to an intermediate control scheme, in which the PID control is disabled and the VP angle is linearly adjusted with respect to the body angle. With this scheme, simulated gaits converge to a VP trend, where the GRF gradually focus on a single point as in Figure 5d and 5e.
The morphological parameters of the TSLIP model are selected from the biomechanical literature to match an ostrich of 80 kg with 1 m leg length (see Table 1). The control gains and damping coefficient of the model are then adjusted to conform certain gait characteristics observed in avian locomotion, which is described in detail in Section 3.2 and 3.3.

TSLIP Model for Avians
The model configuration in Section 2 is used in [9] for the human morphology.
In this section, we underline the modifications made in the TSLIP model Position vectors are referred to as r F H , r F V , r F C , r F H . Angles θ L , θ C , θ V P are the leg, trunk and VP angles. c) Classification of regions for placing the virtual point. There are two major axes that affects the trunk motion: CoM-foot and hip-foot axes. The CoM-foot axis sets the direction of the angular moment applied by the GRF. Virtual points above the CoM-foot axis cause mainly upward trunk motion during stance, whereas points below yield the opposite. VP needs to lie on the vertical axis passing through the CoM for the steady state motion [30]. Hence, the parameter space is divided between the VP's above the CoM (VP A ) and below the CoM (VP B ). The hip-foot axis determines the sign of the hip torque. All VP A lie above hip-foot axis for our parameter range (VP radius up to 60 cm). The VP B region is divided into points above the hipfoot axis, and below. We refer the points below both CoM-foot and hip-foot axes as VP BL . d) The leg spring and damper create the axial component of the ground reaction force F F a , whereas the hip torque creates the tangential component F F t . The vector sum of these forces passes through the virtual point VP. 1 The ostrich silhouette is adapted with the permission from https://proartinc.net/shop/4k-wildlife-relax/4k-ostrich-africa.
and control, in order to accommodate the avian morphology. First we justify our choice of using a bilinear leg damper with the avian gait data. Second we explain the changes in the control strategy of the The literature values of the model parameters are presented for an ostrich. The chosen values match an ostrich of 80 kg. 3 The trunk angle values are calculated from the mean pelvis and hip joint angles from the indicated literature. 4 The leg angle is a function of forward speed, see Figure 7b. 5 We assume the distance between the hip and CoM to be equal to the femur length.
leg angle at touch-down and the condition for the leg take-off. Lastly, we clarify the basis for defining VP position w.r.t. the body and world frames for VP A and VP B , respectively.

Bilinear Damping:
The need for having a leg damper in the our model is related to the energetic consequence of having pronograde trunk orientation. A pronograde trunk with cranial CoM requires the hip to supply power continuously, so that the trunk does not collapse into flexion. This positive power needs to be dissipated over a step cycle to maintain a constant energy level. The standard method to achieve this is by adding a parallel leg damper in the TSLIP model [17,30]. However, the linear leg damper generates non-zero leg length velocity at leg touchdown and take-off events in Figure 2b ( ), which is not consistent with the avian gait data ( ) estimated from [23]. The inaccurate leg length velocity, causes non-zero damping forces at touch-down and take-off events in the simulation. It is not realistic to have a non-zero damping force at leg touch-down, when the leg is at its rest length. In contrast, non-zero damping force is acceptable at leg take-off due to the early leg take-off observed in avian gaits. The early leg take-off can be seen in Figure 2a ( ), where the leg length at take-off is smaller than at touch-down. Consequently, we use the bilinear leg damper in [9] that combines the leg length velocity with leg deflection ( ) to obtain non-zero damping forces at touch-down.

Early Leg Take-off in the Avian Model:
The standard TSLIP model terminates the stance phase when the leg reaches to its rest length. Condition 1 causes an unexpected phenomenon for the pronograde TSLIP model, in particular at slow speeds, which is demonstrated in Figure 3  It results in non-zero damping force at leg touch-down. It is not realistic to have a damping force when the leg is at its rest length. In addition, the simulated length leg velocity does not match the avian gait data( ). Therefore, we use bilinear damper that accounts for both leg length velocityl and its deflection ∆l in the form of F d = c×(l×∆l). As a result, we obtain zero leg length velocity ( ), and therefore zero damping force at touch-down. The non-zero leg length velocity at take-off is caused by the early leg take-off, which is also observed in avian data. In the plots, leg lengths are offset to the same touch-down value, velocities are normalized. The ostrich data ( ) is estimated from the joint angles in Figures 6-10 of [23], using inverse kinematics.
while the leg length is still smaller that its rest length l 0 . Consequently, the stance phase continues and the CoM starts to accelerate downward until the leg is fully extended (Figure 3c-3d, ). In other words, the CoM behaves like an inverted pendulum (IP) towards end of the stance phase. The flight phase from leg take-off to apex diminishes, and the phase from apex to leg-touch down starts at a lower height with negative vertical velocity. The IP-like behavior is uncharacteristic for the running gaits with spring-mass dynamics. In order to prevent the undesired negative vertical CoM velocity and negative vertical GRF after the midstance (MS), we add Cond. 2-3 to the stance phase termination conditions. As a result, the leg length at take-off is shorter than the rest length (Figure 3e-3f, ) and the GRF profiles end suddenly (Figure 3g-3h, ).
In summary, the stance phase is terminated when one of the three conditions below holds. Consequently, the leg takes off earlier, before the CoM reaches to its apex, as shown in Figure 3c-3d ( ).  Figure 3: (a-b) Touch-down and take-off events of the TSLIP model are shown for VP radii of ±40 cm. The gait kinematics (c-f) and kinetics (g-h) are drawn for two case scenarios, where the stance phase terminates at conditions 1-2-3 ( ) and condition 1 only ( ). In general, the leg take-off condition 1 (l ≥l 0 ) is used in the TSLIP model to terminate the stance phase. However, condition 1 alone leads to gaits, where the CoM to reaches its apex and gains negative acceleration towards the end of the stance phase. This IP-like behavior is marked with crosses ( ) in (c-d) and is not a part of the running gait with springmass dynamics. Consequently, we extend the leg takeoff event condition to 2-3 to prevent negative vertical speed of the CoM after the midstance. The extended conditions result in early leg take-off (e-f) and a cut-off in ground reaction forces (g-h). Condition 1. The leg reaches its rest length: l ≥l 0 . Condition 2. The vertical CoM speed reaches zero after midstance. I.e.,ẏ C ≤0. Condition 3. The vertical GRF reaches zero. I.e., GRF vert. ≤0, which presents the unilateral constraint.

Adaptations in the Control Strategy
We use a linear controller to regulate the leg angle at touch-down, which is a function of the forward speed and apex height [2,9]. In the avian TSLIP model, we add a linear dependence of the body angle at apex on the leg angle at touch-down to bound the magnitude of the trunk oscillations. We explain the selection of the controller gains in Section 3.3.
We determine the hip torque using a VP concept, where the VP creates a passive control mechanism that guides the GRF vectors and counteracts the trunk pitch motion (see in Figure 4). If the VP does not provide countering motion , the trunk would either flip back or collapse into flexion, and the motion would fail. A major factor that determines this reaction is the coordinate system in which VP is defined. The VP frame should be selected so that the resultant GRF creates a moment around the CoM in the opposite direction of the trunk motion during at least some part of the stance phase.
In the standard TSLIP model, the VP is defined w.r.t. the body frame. When the trunk is perturbed downward as in Figure 4a-4d, the VP can generate instances with counteracting moment around the CoM for both VP A and VP B . When the trunk is perturbed upward, it is possible for VP A (Figure 4e-4f), but not feasible for VP B (Figure 4g,4h), to counteract the upward moving trunk during stance phase (i.e., in Figure 4g). For VP B , the set VP location cannot provide postural equilibrium. Since the VP angle is non-adaptive, there is no means of recovering from a upward trunk motion. We propose defining the VP B w.r.t. the world frame, which flips the direction of the VP location change w.r.t. the Figure 4: The TSLIP model is shown at touch-down and take-off events for two cases, where the trunk is perturbed downward (a-d) and upward (e-h). The VP has to react to the changes in trunk angle in a stabilizing manner. The sign indicates that the VP will counterbalance the trunk perturbation by creating a moment around the CoM ( , ) in the opposite direction of the trunk motion. denotes that it will not. When the trunk is perturbed downward, VP defined in body frame can provide counterbalancing action for both VP A and VP B . However for upward trunk perturbation the VP B in body frame have no means of preventing the trunk from flipping upward (g). The VP B is defined w.r.t. to the world frame, which flips the relative motion btw. the VP and trunk, and works for all cases. c) e) Figure 5: The simulation setup (a) and resulting GRF patterns (b-e) for the avian TSLIP model. It is difficult to find feasible gaits for the VP scheme, as the control is sensitive to the changes in the initial state and model parameters. To guide the gait search, we use the iterative control framework presented in [9] (refer to Section 2). The framework has an initial step, which combines the VP controller with a PID controller on trunk pitch angle. The resultant gaits do not have a focused GRF profile (b-c). Then, we reduce the PID gains gradually, until the GRF intersect at a single point (d-e).
trunk motion. This way the trunk can stabilize upward perturbations of trunk and can obtain steady state solutions for VP B .

Gait Generation
In our simulation setup, we sweep VP targets over r V P = ±[0, 20, 40, 60] cm 67 to create trunk oscillations with magnitude up to 10°(see Figure 5), which is the maximum value reported for avian running in literature [7,14,22,23]. We define the VP angle θ V P in body frame for VP A and in world frame for VP B , due to reasons given in Section 3.2. The VP concept requires the points to reside on the vertical axis passing through the CoM for steady state gaits [30]. Therefore, VP angles are set to θ V P =−θ C and θ V P =−π respectively. In addition, we set the desired mean body pitch angle to 100°, which is in accordance with the reported values in literature (see Table 1).
We select the model parameters and controller gains of our avian model so that the resulting gaits follow the trends observed in biomechanics. We choose to match the duty factor 8 , the GRF profile and 6 The parameter sweep for VP B at 10 m s −1 ends at r V P =52 cm due to instability caused by high angular trunk acceleration. 7 The negative r V P denotes that the VP is below the CoM. 8 Duty factor is the ratio of the leg contact time to the stride period values similar to ones reported in the biomechanical measurements in literature. Data points for the duty factor are retrieved from [37]. They cover the speeds 0.25-5 m s −1 , and show two trends: a lower valued curve for small birds ( ) and a higher valued curve for ratites ( ). We fit two curves to these data sets, which determine the lower and upper bound of the grey shaded region. We use the duty factor curve digitalized from [13] to extend this region to higher speeds up to 8 m s −1 ( ). The duty factor of our simulated gaits ( , ) lies in the grey region and decreases with the forward speed. We monitor the peak vertical GRF values in a similar fashion.
the leg angle at touch-down; which are a function of the forward speed. As forward speed increases, the duty factor gets smaller, magnitude of the GRF gets higher, and leg touch-down angle gets smaller [11,13,15,35,37,38]. We tune the damping coefficient and control gains of our model so that the resulting gait behavior is feasible and follows the trend reported for avian locomotion, in Figure 6. Simulated gaits ( , ) have a duty factor of 40-60 % in Figure 6a, and peak vertical GRF of 2-5 BW in Figure 6b. The leg touch-down angle and damping coefficient are in the range of 60-40°and 6-1.5 kN s m −1 and decrease with speed (see Figure 7a).

Simulation Results
In this section, we analyze the results of our simulation setup to investigate the effect of trunk oscillations.

Asymmetries in the Kinetics and Kinematics
The trunk inclination (θ C in Figure 1b) introduces asymmetries in the system, which are reflected in the leg dynamics and GRF profiles [17,19]. This effect is pronounced for the avian model, where the CoM is placed cranially with an inclination of 100 • .
The asymmetry is apparent for the leg angles at touch-down and take-off in Figure 7b, where the leg takes off before the leg length reaches to its resting value, l 0 . The difference between leg lengths at touchdown and take-off varies between 4-8 % of the l 0 (see max. length and ( , ) in Figure 8). The stance phase is terminated early either due to loss of foot contact caused by the vertical GRF decaying to zero, or due to the CoM reaching to its apex height (Cond. 2-3). As a result of the early leg take-off, the leg spring is unable to inject all stored energy. Effectively, the leg spring removes energy from the system, details of which we discuss in Section 4.3.1.
In our gait framework, we obtain left-skewed GRF profiles shown in Figure 9, which is consistent with the avian locomotor data in [5,17,28]. In the following, we investigate potential reasons. First, we subtract the component of the GRF created by the leg damper (F dp = F F a −F sp ) from the total GRF in Figure 9a and 9b (solid lines). We observe that the damping force skews vertical GRF to the left (Figure 9a, dotted lines, ). The effect of the damping is visible mainly in vertical GRF profile, because the leg force makes up most of the vertical GRF and only a minor part of horizontal GRF.
Next, we subtract the component of the GRF produced by the hip torque ( F F t ) in Figure 9c

Leg deflection
Leg take-off length Figure 8: The min. and max. values of the leg length are plotted as bar plot, where the leg take-off length is marked with dots ( , ). We observe an early leg takeoff with a shorter leg, which is approximately 4-8 % of the leg rest length l 0 . The leg deflection corresponds to the length of each bar and increases from 15 to 30 % of l 0 with the running speed. The GRF is the sum of the axial force created by the leg ( F F a ) and the tangential force created by the hip ( F F t ). In (a-b, solid lines), we subtract the force generated by the leg damper (F dp = F F a − F sp ) form the GRF. We observe that the damping causes a left skew in vertical GRF (a, ), and has no significant effect for the horizontal GRF. In (c-d), we subtract the tangential force generated by the hip ( F F t ) from the GRF. We observe that the hip torque causes a left skew towards left in horizontal GRF (d, ), and has no significant effect on the vertical GRF.
dotted lines, ). The integral of the horizontal GRF curve (the area) is the the fore-aft impulse, which corresponds to the forward acc/deceleration of the main body. The horizontal GRF created by the leg force causes larger negative horizontal impulse than A VP A ( ) yields higher peak vertical GRF and lower peak horizontal GRF, compared to VP B ( ). The magnitude of the peak GRF increases with speed, which is quantified in (c,e). The phase of the peak GRF (d,f) is calculated w.r.t. the step time and it increases with speed. In other words, the GRF become more symmetrical at higher speeds. the positive one, which would effectively decelerate the main body. An example is presented in Figure 9d, where the positive and negative impulses generated by the horizontal GRF (solid lines) amount to 3/−48 N s for VP A (blue) and 4/−52 N s for VP B (red) respectively. The hip torque (dashed line) produces forces to ensure sufficient forward acceleration to generate the motion. In other words, it creates equal positive and negative fore-aft impulses by shifting the zero crossing of horizontal GRF to the left.
We then investigate how the skew of the GRF profile changes with the forward speed in Figure 10. As the forward speed gets higher, the magnitude of peak horizontal and vertical GRF increase between 2.5-4.5 BW and 0.5-2 BW, respectively (Figure 10c,10e). The phase of peak GRF is calculated with respect to the gait cycle (marked with ( , ) in Figure 10a-10b). The phase of peak horizontal and vertical GRF increases between 42-84 % and 32-40 %, respectively (Figure 10d,10f). The increase in phase indicates that gaits become more symmetric at high speeds.
In terms of the VP location, VP A and VP B have similar phase values for peak GRF. VP B yields lower peak vertical GRF and higher peak horizontal GRF magnitudes shown in Figure 10a-10b, which are associated with high duty factor and high horizontal accelerations, respectively. In other words, VP B makes The trunk angular excursion (i.e., the absolute value of the difference between min. and max. trunk angles) and the peak angular pitch velocity increases with the VP radius. The stance phase is shaded in gray.
the CoM brake and accelerate more during stance phase in the horizontal direction.
Apart from the asymmetries mentioned, we detect no take-off-apex phase occurring at low speeds, around 4 m s −1 with high duty factor (see grey shaded area in Figure 11a-11b, or Figure 3c-3d). The CoM reaches its peak height at the end of the stance phase and the next step begins immediately after, which is enabled with the take-off condition Cond. 2-3 in Section 3.1.

Trunk Pitch Oscillations
In our simulations, VP B leads to -shaped, downward trunk pitch motion in during the stance phase, which are similar to the avian gait characteristics reported in [8,22]. VP A yields -shaped, upward trunk motions in Figure 11a and 11b. For equal VP radius, VP B causes larger pitch oscillations. The magnitude of the oscillations and the angular velocity increase with the VP radius for both VP A and VP B (Figure 11c-11d).
When the running speed increases from 4 to 10 m s −1 , the trunk angular excursion increases up to 2°f or VP A and 18°for VP B (see Figure 12a). The mean trunk angular velocity increases up to 13°s −1 for VP A and −91°s −1 for VP B (see Figure 12b).  Figure 12: The trunk angular excursion (i.e., the absolute value of the difference between min. and max. trunk angles) and mean angular pitch velocity increase with forward speed and increasing VP radius.

Energy Considerations
In this section, we investigate how the leg and hip contributes to the system's energy balance, and how the CoM energy evolves over step time. Moreover, we provide the mechanical cost of transport to allow comparison to the literature. To clarify, we use the term positive/negative work to address the amount of energy created/absorbed by the leg force or the hip torque. We define the energy fluctuation (∆) as the difference between maximum and minimum values of any energy type.

Work Distribution Between Leg and Hip
Temporal Analysis: We explore how the leg force, hip torque, and their respective energies evolve over the course of the stance phase in Figure 13 and 15j-15l. The leg spring deflects and stores energy during the first half of the stance phase (Figure 13g, 15j, ). It extends and returns this energy back to the main body in the second half (Figure 13g, ). However, owing to the early leg take-off, the spring is not able to recoil completely at the end of the stance phase ( Figure 13a) and return all the energy it absorbed (refer to Section 4.1 and see Figure 13g). Consequently, the spring has a net effect of removing energy from the system, which is indicated by the arrows (Figure 13g, ). Concerning the leg damper, early leg take-off interrupts the energy absorption of damper (Figure 13h,15k, ) and makes the damping force end abruptly at a non-zero value (Figure 13e, ).
The hip actuator has two purposes: to compensate for the energy losses of the leg and to provide positive net work to support the forward leaning trunk against gravity. The posterior placement of the hip w.r.t. the CoM for a pronograde trunk necessitates high positive hip work to hold the heavy trunk [12,19]. Accordingly, we see in Figure 13f that the hip torque The net work performed on the system has to remain zero to obtain steady state motion. The task of energy depletion is achieved primarily by the leg damper (h, ), and partly by the leg spring (g, ) owing to the early take-off (a). VP A causes higher positive net hip work and negative leg damper work (i,h, ), compared to VP B ( ). In contrast, VP A yields lower negative net spring work (g, ).
is always negative for VP A , and VP B with radii smaller than 30 cm. The hip solely injects energy to the system (Figure 13i,15l, ). If the VP B radius is larger than 30 cm (VP BL ), the VP is set below the leg axis at touch-down. In this case, the hip starts producing positive torque in the first half of the stance phase (Figure 13f), where the hip depletes energy (Figure 13i, ) and assists the leg force to decelerate the body. Beyond this initial negative work, the hip still has to produce net positive work to offset the energy absorbed by the leg (Figure 13i, ).
The requirement for the net positive hip work comes from the steady state condition, where the net  Figure 14: The work performed by the leg (a-f) and hip (g-l) as a function of the VP location for the speeds of 4 and 10 m s −1 . The leg produces net negative work, whereas the hip produces net positive work. At slow speeds, the positive, negative and net work performed by the leg (a-c) increases with a higher VP A radius and decreases with a higher VP B radius. The relation is inverse for the fast speeds (d-f). The positive and net hip work displays a relation similar that of the leg (see g,i and j,l) for VP A and VP B with a radius smaller than 30 cm (i.e., above the leg axis). For these VP targets, the hip produces no negative work (h,k). When the VP radius is set higher than 30 cm , the hip starts to generate negative work (h and k). The excess work is compensated with an increase in the positive hip work to maintain steady state conditions. In other words, the peak work requirement of the hip increases, while the net hip work follows the prior trend.
work done on the system must be zero 9 . Otherwise the body would accelerate or decelerate. Following this, we see that the energy injected by the hip actuator is equal to the energy absorbed by the leg in Figure 14.
The leg produces positive work via the spring and negative work through both the spring and damper, which results in a net negative work (Figure 14a-14b). The contribution of the leg spring to the overall energy removal is relatively small and amounts to less than 5 % (Figure 13g-13i). 9 Asymmetric or period-2 gaits are not considered.
Spatial Analysis: We investigate the relation between the VP location and leg-hip work. We mark that this relation changes with the forward running speed (see Figure 14). Therefore, we split our analysis into two parts involving slow and fast speeds of 4-6 m s −1 and 8-10 m s −1 . In the following, we derive statements for slow speeds and all effects described are reversed for the fast speeds. We start our analysis with the leg. As the VP radius increases from 0 to 60 cm; the positive, negative, and net leg work magnitudes (Figure 14a) increase 3.4, 2.2, 2.2 % 8 for VP A ( ), and decrease 20, 16, 13 % for VP B ( ), respectively.
The hip generates only positive work for VP A (see Figure 14g-14i, ). The lower the VP A radius, the less work the hip performs. Hip work increases about 2.2 % with the increasing VP radius. VP B up to the radius of 30 cm reduces the work requirement of the hip further about 4 % (Figure 14g). But, when the VP B radius is bigger than 30 cm (VP BL ), the hip starts to generate negative work (Figure 14g-14h) and needs to produce large amounts of positive work to compensate for both its own negative work and leg's. For example, at a VP BL radius of 60 cm the positive hip work requirement increases 30 % (Figure 14g). A larger VP BL radius yields higher positive and negative hip work independent from the speed. On the other hand, a larger VP BL reduces the net hip work at slow speeds (Figure 14i), which creates a trade-off between the peak torque demand (Figure 13f) and mean energy expenditure. Such a trade-off diminishes at fast speeds of 10 m s −1 (see Figure 14j-14l).
When the forward speed increases from 4 to 10 m s −1 , the mean values of the spring, damper and hip energies show approximately 2-fold, 3-fold and 2.5-fold increase respectively (see , Figure 15d-15f).

Energy of the CoM
In this section, we look at how kinetic and potential energies of the CoM change as a function of forward speed (on the left of Figure 15) and the normalized step time (on the right). Figure 15g-15i show that VP A causes smaller fluctuations in the linear and angular kinetic energies and higher fluctuations in potential energy, compared to VP B . We show the energy fluctuations (∆) that correspond to Figure 15a-15c in Figure 16a-16c and observe that VP A minimizes the linear kinetic and total energy fluctuations of the CoM.
When we increase the forward speed from 4 to 10 m s −1 , the mean linear and angular kinetic energies show approximately 6-fold and 2-fold increase, while the potential energy shows 1.15-fold decrease (see , Figure 15a-15c). In addition, the fluctuations within the distinct types of energies increase, which are indicated with the magnitude of the bar plot in Figure 15 and are quantified numerically in Figure 16. VP A has lower linear kinetic and higher potential energy fluctuations compared to VP B . Fluctuations in angular kinetic energy depend on the VP radius.

The Mechanical Cost of Transport
In robotics and biomechanics, the cost of transport (CoT) is a measure to compare the energy efficiency of different species and robots [39,40]. It is defined as the energy used per traveled distance and expressed as CoT= P mgẋ in dimensionless form, where P involves the overall power generated by the metabolism, muscles/actuators and so on. A smaller CoT indicates a better energy economy.
In our model, the actuator is the hip motor and the CoT reflects hip work. We consider three cases, where we use the (Figure 17a) net, (Figure 17b) positive, and (Figure 17c) absolute values of hip power to calculate the CoT ( Figure 17). As the running speed increases from 4 to 10 m s −1 , the CoT increases from 0.25 to 0.65, 1.7 and 1.6 for (a-b-c) respectively. In addition, we see that VP radii between −10 cm and 30 cm yield a small CoT. For higher VP radii, the CoT increases. We measure the power expenditure in hip joint space. Our measure involves both mechanical and partially metabolic CoT: the hip compensates for the damping in the leg and captures some of the metabolic effects indirectly.

Discussion
In this section, we suggest control design rules for determining the VP location, and discuss the benefits and drawbacks of each choice. We interpret the VP concept and our design rules in relation to the GRF alignment. Moreover, we explain if and how our simulation results conform to the biomechanical observations regarding the avian trunk oscillations and VP.

Control Design Rules
We suggest the speed dependent design rule in Table 2 for avian bipedal locomotion, which modifies trunk-leg loading through trunk pitch movements.
At slow running speeds of 4-6 m s −1 , we suggest using a VP B with a radius smaller than 30 cm to minimize the hip work. For example, assume a case where we have a robot, and its motion planner outputs a desired hip torque that exceeds the motor limit. We can make this motion feasible by integrating downward trunk oscillations to our control design, which would reduce the peak motor torque demand. Reduced torque demand allows the usage of smaller motors and hence the design of lightweight robots. Alternatively, the strategy could be used in rehabilitation for a patient with weak hip extensor strength to assist the hip. On the other hand, we suggest utilizing a VP BL to minimize the leg loading, which would also require high peak hip torques. VP BL can be beneficial for cases where the leg has to be light (e.g., in robot design), or where the leg actuation is weak (e.g., in case of injuries).
At faster running speeds, the system behavior of the pronograde trunk, and hence our disposition of the VP changes. The relation between the VP location and leg work reverses after 8 m s −1 and the hip after 10 m s −1 . In this case, VP A with a high radius minimizes both the leg and hip work. Consequently, we propose VP A for energy efficiency at high speeds.
In terms of the CoM energetics, VP A minimizes the horizontal and overall energy fluctuations of the CoM across all speeds. VP A could be useful for cases Figure 18: Four potential strategies for the GRF alignment are illustrated for the foot touch-down and take-off events of the single stance phase of running. The moment created by the GRF vector around the CoM is indicated with ( , ), whereas the foreaft acceleration/deceleration created by the GRF is indicated with ( , ) arrows. As where the horizontal accelerations are not desired, for instance when there are changes in the ground level [41], or when a biped is carrying a fragile load. VP B can be selected to minimize the vertical energy fluctuations of the CoM, for instance when the biped stumbles forward and vertical acceleration is not desired for the motion recovery.

The VP as a method for GRF Alignment
In this section, we offer a different perspective for interpreting the VP, where changing VP location accounts for modifying the ratio between horizontal and vertical GRF. We refer to this as GRF manipulation or alignment.
In biomechanics, the GRF manipulation is considered as the running technique of an animal, where the animal adjusts the forces it applies on the ground. In terms of avian locomotion, two potential GRF alignment strategies were assessed for quail in [32] (see Figure 18). The first strategy has purely vertical GRF (A) and the second strategy points the GRF towards to CoM (B). The first strategy only creates a moment around the CoM ( , ), which rotates the body in the pitch direction and causes fluctuations in the angular kinetic energy. The second strategy yields no angular motion but decelerates ( ) and accelerates ( ) the main body successively. In other words, it causes fluctuations in the translational kinetic energy. It is known that quails employ (B) in running, despite (A) being energetically more economical (i.e., (A) resulting in less kinetic energy fluctuations) [32]. This preference is motivated by the the excessive pitch angle demand of (A), which is physically not feasible due to the power limit of the hip actuation. • When stability in vertical axis is needed e.g., forward perturbations such as stumbling Table 2: Control design rules for avian bipedal locomotion with TSLIP dynamics. These rules also hold for the human morphology, with the exception of R3 [9]. The human model shows no speed dependency in energetics. Rules R1-R2 are valid for all speeds.
The VP concept lies between (A) and (B), where the GRF both induces rotation around the CoM and acc/decelerates the CoM in the horizontal direction. In our framework, GRF vectors become more vertically oriented as VP location transits from VP B to VP A , as (D) (B) (C) (A). VP A with more vertical GRF alignment yields smaller fluctuations in the linear and overall CoM energy, in accordance with [32,33].
In the light of this new perspective, we reinterpret our results in Figure 14 of Figure 4.3.1 for the avian morphology. Both the leg and hip work can be reduced, if the GRF vectors are oriented more vertically at fast speeds and more horizontally at slow speeds.

Gait Measurements vs. TSLIP Model
A set of biomechanical experiments report downward trunk motion during the stance phase of avian running, by using methods such as cineradiography or motion capture [8,14,22,23]. Another independent set of experiments estimate a VP above the CoM (VP A ) from the GRF measurements [17,19,25]. However, when we implement a VP A in the TSLIP model, the simulated gaits display upward trunk motion during the stance phase. This is in conflict with the first set of biomechanical observations (see Table 3). In the TSLIP simulation, a forward trunk motion is obtained when the VP is set below the CoM (VP B ).
In summary, there is a mismatch between the biomechanical observations and TSLIP model regarding the coupling between the trunk oscillation direction and VP location. One possible explanation is related to the disparity between the trunk and whole body dynamics. In human walking, the trunk pitching motion is reported to be in antiphase (180°out of phase) with the whole body pitching motion [31]. Given that humans have relatively heavy limbs [42], it is plausible that the trunk and whole body dynamics deviate from each other. The TSLIP model with a VP A and upward trunk motion might reflect the whole body dynamics, and not necessarily the trunk dynamics [43]. However, this argument is not convincing for the avian species with relatively light legs, where the lower extremities contribute little to the whole body dynamics [20]. The existing research is missing data to provide an antiphase correlation between the trunk and whole body pitch dynamics for human running or any kind of avian gait.
Another potential explanation is related to the data processing of the GRF measurements. GRF signals are noisy, especially at touch-down and take-off events, due to the artifacts that results from the impact, heel-strike and ankle push-off. Consequently, GRF that belong to the initial and final phases of stance phase are removed when estimating the VP [44]. In addition, GRF data is filtered to remove noise and drift. The VP is calculated as the point that minimizes the distances between GRF vectors. Truncation and modification of GRF signals might cause an error in estimating the VP. There is a need for systematic experiments that are tailored to investigate how avians and humans modify the GRF over a wide range of speeds, and how the trunk and whole body motion fits to this framework.
The TSLIP we use in our analysis is a simplified model, which can predict the dynamics of running for the CoM, close to what we observe in nature. The choice of such a simplified model for our analysis creates the question, how well the model predictions conform to more complex models and actual underlying dynamics in real-life cases. Despite these potential drawbacks, the simplification in TSLIP enables us to isolate the function of the trunk and leave out the inertial effect of other extremities. We can investigate the effect of trunk motion across various bipedal species with different leg characteristics.

Conclusion
In this paper we investigated how trunk pitch oscillations affect the dynamics and energetics of running in terrestrial birds, who possess a horizontal (pronograde) trunk. We used an avian TSLIP model with a virtual point (VP) control target to create trunk oscillations and analyzed how the VP location affects the energy economy.
We placed the VP above (VP A ) and below (VP B ) the center of mass (CoM), and performed a parameter sweep on the VP radius to assess its effect on energetics. We showed that a VP A causes upward and a VP B causes downward trunk motion during the stance phase.
Our results suggest different potential strategies to place the VP in order to shape the energy distribution. One can use a VP A to reduce the kinetic energy fluctuations of the CoM, while a VP B reduces the potential energy fluctuations. Both of these come at the expense of increased angular energy fluctuations of the CoM. At the same time, one can reduce the hip and leg work jointly by using a VP A at fast speeds of 8-10 m s −1 and a VP B at slow speeds of 4-6 m s −1 . A VP B below the leg axis (VP BL ) is useful to reduce the leg work further at slow speeds of 4-6 m s −1 at the expense of the high peak hip torques.
The aim of our study was to show how trunk motion can be leveraged to shape the energy distribution. As future work, we plan to extend our simulation analysis to other bipedal morphologies and test the validity of our control strategies with real robots.

Acknowledgement
The authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supportingÖzge Drama. This work was made possible thanks to a Max Planck Group Leader grant awarded to A. Badri-Spröwitz by the Max Planck Society.