Controlling locomotion with CPGs

Note

Authors: Sibo Wang-Chen, Femke Hurtak

The code presented in this notebook has been simplified and restructured for display in a notebook format. A more complete and better structured implementation can be found in the examples folder of the FlyGym repository on GitHub.

This tutorial is available in .ipynb format in the notebooks folder of the FlyGym repository.

Summary: In this tutorial, we will introduce the concept of Central Pattern Generators (CPGs) and build a CPG-based model to control walking of the simulated fly.

Central Pattern Generators (CPGs) are neural circuits that generate rhythmic output without receiving rhythmic input (see review of CPGs in insects by Mantziaris et al, 2020). CPGs have been identified in various species as playing a crucial role in locomotion. Their principles have also been adopted in the realm of robotic motor control (see review by Ijspeert, 2008). It is hypothesized that CPGs play a more important role in animals whose locomotion is fast and therefore cannot process decentralized limb sensory signals in time to adjust their movements (i.e., like cockroaches running at ~1 m/s!).

In this tutorial, we will start by implementing an oscillator network. Then, we will use this network to control the stepping of the legs by mapping the phases of the oscillators to the phases of preprogrammed leg stepping kinematics. Next, we will plug this oscillator-based controller into NeuroMechFly to drive locomotion in the physics simulation. Finally, we will add leg adhesion to the locomotor simulation.

The CPG network

Basic CPGs can be implemented as feedforward networks of oscillators—in other words, the network behaves without taking into account sensory feedback. Similar to the formulation from Ijspeert et al (2007), the oscillator network can be described by the following ordinary differential equations (ODEs):

\[\dot\theta_i = 2\pi\nu_i + \sum_{j} r_j w_{ij} \sin(\theta_j - \theta_i - \phi_{ij})\]

\[\dot r_i = \alpha_i (R_i - r_i)\]

where \(\theta_i\) and \(r_i\) are the current phase and magnitude of the i-th oscillator. \(R_i\) is the maximum amplitude of the i-th oscillator, and \(\alpha_i\) is a constant determining the rate of convergence to synchrony. \(w_{ij}\) is the coupling weight between the i-th and the j-th oscillator, and \(\phi_{ij}\) is the phase bias between them. Intuitively, the first term of the first equation maintains an intrinsic frequency for each oscillator; the second term of the first equation keeps the oscillators synchronized (i.e., maintains the phase differences between the oscillators), and the second equation maintains the amplitudes of the oscillators.

To start, let’s write a function that computes \(\dot\theta\) and \(\dot r\). For explicitness, we will implement the CPG network using only Python and NumPy (no scipy.integrate).

import numpy as np


def calculate_ddt(theta, r, w, phi, nu, R, alpha):
    """Given the current state variables theta, r and network parameters
    w, phi, nu, R, alpha, calculate the time derivatives of theta and r."""
    intrinsic_term = 2 * np.pi * nu
    phase_diff = theta[np.newaxis, :] - theta[:, np.newaxis]
    coupling_term = (r * w * np.sin(phase_diff - phi)).sum(axis=1)
    dtheta_dt = intrinsic_term + coupling_term
    dr_dt = alpha * (R - r)
    return dtheta_dt, dr_dt

Next, let’s implement the CPG as a class. We will use Euler’s method to integrate the ODE, but you can use any higher-order methods or libraries.

class CPGNetwork:
    def __init__(
        self,
        timestep,
        intrinsic_freqs,
        intrinsic_amps,
        coupling_weights,
        phase_biases,
        convergence_coefs,
        init_phases=None,
        init_magnitudes=None,
        seed=0,
    ) -> None:
        """Initialize a CPG network consisting of N oscillators.

        Parameters
        ----------
        timestep : float
            The timestep of the simulation.
        intrinsic_freqs : np.ndarray
            The intrinsic frequencies of the oscillators, shape (N,).
        intrinsic_amps : np.ndarray
            The intrinsic amplitude of the oscillators, shape (N,).
        coupling_weights : np.ndarray
            The coupling weights between the oscillators, shape (N, N).
        phase_biases : np.ndarray
            The phase biases between the oscillators, shape (N, N).
        convergence_coefs : np.ndarray
            Coefficients describing the rate of convergence to oscillator
            intrinsic amplitudes, shape (N,).
        init_phases : np.ndarray, optional
            Initial phases of the oscillators, shape (N,). The phases are
            randomly initialized if not provided.
        init_magnitudes : np.ndarray, optional
            Initial magnitudes of the oscillators, shape (N,). The
            magnitudes are randomly initialized if not provided.
        seed : int, optional
            The random seed to use for initializing the phases and
            magnitudes.
        """
        self.timestep = timestep
        self.num_cpgs = intrinsic_freqs.size
        self.intrinsic_freqs = intrinsic_freqs
        self.intrinsic_amps = intrinsic_amps
        self.coupling_weights = coupling_weights
        self.phase_biases = phase_biases
        self.convergence_coefs = convergence_coefs
        self.random_state = np.random.RandomState(seed)

        self.reset(init_phases, init_magnitudes)

        # Check if the parameters have the right shape
        assert intrinsic_freqs.shape == (self.num_cpgs,)
        assert coupling_weights.shape == (self.num_cpgs, self.num_cpgs)
        assert phase_biases.shape == (self.num_cpgs, self.num_cpgs)
        assert convergence_coefs.shape == (self.num_cpgs,)
        assert self.curr_phases.shape == (self.num_cpgs,)
        assert self.curr_magnitudes.shape == (self.num_cpgs,)

    def step(self):
        """Integrate the ODEs using Euler's method."""
        dtheta_dt, dr_dt = calculate_ddt(
            theta=self.curr_phases,
            r=self.curr_magnitudes,
            w=self.coupling_weights,
            phi=self.phase_biases,
            nu=self.intrinsic_freqs,
            R=self.intrinsic_amps,
            alpha=self.convergence_coefs,
        )
        self.curr_phases += dtheta_dt * self.timestep
        self.curr_magnitudes += dr_dt * self.timestep

    def reset(self, init_phases=None, init_magnitudes=None):
        """Reset the phases and magnitudes of the oscillators.
        High magnitudes and unfortunate phases might cause physics error
        """
        if init_phases is None:
            self.curr_phases = self.random_state.random(self.num_cpgs) * 2 * np.pi
        else:
            self.curr_phases = init_phases

        if init_magnitudes is None:
            self.curr_magnitudes = np.zeros(self.num_cpgs)
        else:
            self.curr_magnitudes = init_magnitudes

To demonstrate this network, let’s simulate a network of three oscillators connected as follows:

For the sake of illustration, let’s make them oscillate at an intrinsic frequency of 1 and intrinsic amplitudes of 1.0, 1.1, 1.2. They are coupled with a weight of 1 and phase differences of 120 degrees. We will initialize the phases and magnitudes randomly.

intrinsic_freqs = np.ones(3)
intrinsic_amps = np.array([1.0, 1.1, 1.2])
coupling_weights = np.array(
    [
        [0, 1, 0],
        [1, 0, 1],
        [0, 1, 0],
    ]
)
phase_biases = np.deg2rad(
    np.array(
        [
            [0, 120, 0],
            [-120, 0, 120],
            [0, -120, 0],
        ]
    )
)
convergence_coefs = np.ones(3)

network = CPGNetwork(
    timestep=1e-3,
    intrinsic_freqs=intrinsic_freqs,
    intrinsic_amps=intrinsic_amps,
    coupling_weights=coupling_weights,
    phase_biases=phase_biases,
    convergence_coefs=convergence_coefs,
)

num_steps = int(10 / network.timestep)
phase_hist = np.empty((num_steps, 3))
magnitude_hist = np.empty((num_steps, 3))

# Simulate the network
for i in range(num_steps):
    network.step()
    phase_hist[i, :] = network.curr_phases
    magnitude_hist[i, :] = network.curr_magnitudes

We can visualize the phases (wrapped to \([0, 2\pi]\)) and the magnitudes of the oscillators over time. We observe that, after a brief period of synchronization, the oscillators converge to a state where they oscillate 1/3 of a cycle apart at their intrinsic frequencies and amplitudes.

from pathlib import Path
import matplotlib.pyplot as plt

output_dir = Path("./outputs/cpg_controller")
output_dir.mkdir(exist_ok=True, parents=True)

fig, axs = plt.subplots(2, 1, figsize=(5, 5), sharex=True)
t = np.arange(num_steps) * network.timestep
axs[0].plot(t, phase_hist % (2 * np.pi), linewidth=1)
axs[0].set_yticks([0, np.pi, 2 * np.pi])
axs[0].set_yticklabels(["0", r"$\pi$", r"$2\pi$"])
axs[0].set_ylabel("Phase")
axs[1].plot(t, magnitude_hist, linewidth=1)
axs[1].set_ylabel("Magnitude")
axs[1].set_xlabel("Time (s)")
fig.savefig(output_dir / "simple_cpg_rollout.png")

We have now built a CPG network. In the next section, we address how the states of the CPGs can be used to drive locomotion.

Controlling leg stepping with CPGs

The state variables \(\theta\) and \(r\) can be used to drive locomotion at various levels of abstraction. This is a design choice that the modeler should make depending on the scientific question being considered. For example, in Lobato-Rios et al (2022), the CPG states are used to calculate motor neuron activity \(M_i = r_i (1 + \sin(\theta_i))\), which is in turn used to drive a muscle model. By contrast, Ijspeert et al (2007) uses a more abstract control strategy — the CPG states directly control the target joint position (i.e. angle) \(x_i = r_i (1 + \cos(\theta_i))\). This target position is then provided to a proportional-derivative (PD) controller which actuates the joint.

Here, we will use an even higher-level control approach where each oscillator controls the stepping of an entire leg (as opposed to a joint). The phase of the CPG represents the phase of the step (i.e. how far into the step the leg is), while the magnitude of the CPG represents the magnitude of the step (i.e. how large the step is). We will use experimentally recorded data to execute the individual steps. In other words, we will extract the kinematics of a single step for each leg from experimental behavior recordings and modify its magnitude (modulated by \(r\)) and speed (modulated by \(\theta\)) so that the stepping of the six legs is coordinated by the CPG network.

We will set the coupling parameters for locomotion using a “tripod gait”: at each point in time, the fore and hind legs on one side and the mid leg on the other side of the body are in stance, forming a stable tripod-shaped structure; the other three legs are in swing. This is illustrated in the figure below (left, figure adapted from Emanuel et al, 2020). The tripod gait can be implemented using a CPG network shown on the right. We observe that the legs that should not swing together are coupled with a phase difference of 180 degrees, ensuring that they are out of phase once the network is synchronized. We will use other parameters from the NeuroMechFly v2 paper.

As before, we will set up the CPG network, run the simulation, and plot the time series of the state variables:

intrinsic_freqs = np.ones(6) * 12
intrinsic_amps = np.ones(6) * 1
phase_biases = np.pi * np.array(
    [
        [0, 1, 0, 1, 0, 1],
        [1, 0, 1, 0, 1, 0],
        [0, 1, 0, 1, 0, 1],
        [1, 0, 1, 0, 1, 0],
        [0, 1, 0, 1, 0, 1],
        [1, 0, 1, 0, 1, 0],
    ]
)
coupling_weights = (phase_biases > 0) * 10
convergence_coefs = np.ones(6) * 20

network = CPGNetwork(
    timestep=1e-4,
    intrinsic_freqs=intrinsic_freqs,
    intrinsic_amps=intrinsic_amps,
    coupling_weights=coupling_weights,
    phase_biases=phase_biases,
    convergence_coefs=convergence_coefs,
)

# Simulate the network
num_steps = int(1 / network.timestep)
phase_hist = np.empty((num_steps, 6))
magnitude_hist = np.empty((num_steps, 6))
for i in range(num_steps):
    network.step()
    phase_hist[i, :] = network.curr_phases
    magnitude_hist[i, :] = network.curr_magnitudes

# Visualize
fig, axs = plt.subplots(2, 1, figsize=(5, 5), sharex=True)
t = np.arange(num_steps) * network.timestep
axs[0].plot(t, phase_hist % (2 * np.pi), linewidth=1)
axs[0].set_yticks([0, np.pi, 2 * np.pi])
axs[0].set_yticklabels(["0", r"$\pi$", r"$2\pi$"])
axs[0].set_ylabel("Phase")
axs[1].plot(t, magnitude_hist, linewidth=1)
axs[1].set_ylabel("Magnitude")
axs[1].set_xlabel("Time (s)")
fig.savefig(output_dir / "tripod_cpg_rollout.png")

Now, let’s load the behavior kinematics data:

import pickle
from flygym.util import get_data_path

single_steps_path = (
    get_data_path("flygym", "data") / "behavior/single_steps_untethered.pkl"
)
with open(single_steps_path, "rb") as f:
    single_steps_data = pickle.load(f)

This gives us a dictionary containing joint angle time series for each joint. We will check if they all have the same length. The steps should be periodic, so we will also check if the first and last angles in the time series are the same:

preprogrammed_steps_length = len(single_steps_data["joint_LFCoxa"])
preprogrammed_steps_timestep = single_steps_data["meta"]["timestep"]
print(
    f"Preprogrammed steps have a length of {preprogrammed_steps_length} steps "
    f"at dt={preprogrammed_steps_timestep}s."
)
for k, v in single_steps_data.items():
    if k.startswith("joint_"):
        assert len(v) == preprogrammed_steps_length
        assert v[0] == v[-1]

Preprogrammed steps have a length of 45 steps at dt=0.003s.

Now, for each leg \(i\), let’s build a function \(\Psi_i\) such that given the current stepping phase \(\theta_i\) of the leg, \(\Psi_i(\theta_i)\) provides joint angles of all DoFs on leg \(i\) based on the preprogrammed stepping kinematics. We will do this by interpolation and normalize \(\theta\) to the range \([0, 2\pi)\):

from scipy.interpolate import CubicSpline

legs = [f"{side}{pos}" for side in "LR" for pos in "FMH"]
dofs_per_leg = [
    "Coxa",
    "Coxa_roll",
    "Coxa_yaw",
    "Femur",
    "Femur_roll",
    "Tibia",
    "Tarsus1",
]
phase_grid = np.linspace(0, 2 * np.pi, preprogrammed_steps_length)
psi_funcs = {}
for leg in legs:
    joint_angles = np.array(
        [single_steps_data[f"joint_{leg}{dof}"] for dof in dofs_per_leg]
    )
    psi_funcs[leg] = CubicSpline(phase_grid, joint_angles, axis=1, bc_type="periodic")

We can then map the phase of the CPGs to the phase of the legs. Let’s visualize three stepping cycles for each leg:

theta_ts = np.linspace(0, 3 * 2 * np.pi, 10000)

joint_angles_by_leg = {}
for leg, psi_func in psi_funcs.items():
    joint_angles_by_leg[leg] = psi_func(theta_ts)

fig, axs = plt.subplots(3, 2, figsize=(7, 5), sharex=True, sharey=True)
for i_side, side in enumerate("LR"):
    for i_pos, pos in enumerate("FMH"):
        leg = f"{side}{pos}"
        ax = axs[i_pos, i_side]
        psi_func = psi_funcs[leg]
        joint_angles = np.rad2deg(joint_angles_by_leg[leg])
        for i_dof, dof_name in enumerate(dofs_per_leg):
            legend = dof_name if i_pos == 0 and i_side == 0 else None
            ax.plot(theta_ts, joint_angles[i_dof, :], linewidth=1, label=legend)
        if i_pos == 2:
            ax.set_xlabel("Phase")
            ax.set_xticks(np.pi * np.arange(7))
            ax.set_xticklabels(["0" if x == 0 else rf"{x}$\pi$" for x in np.arange(7)])
        if i_side == 0:
            ax.set_ylabel(r"DoF angle ($\degree$)")
        ax.set_title(f"{leg} leg")
        ax.set_ylim(-180, 180)
        ax.set_yticks([-180, -90, 0, 90, 180])
fig.legend(loc=7)
fig.tight_layout()
fig.subplots_adjust(right=0.8)
fig.savefig(output_dir / "three_steps_phase_only.png")

We can also modulate the amplitude of the steps using the magnitude \(r\) of the CPGs. To do this, we take the difference of the joint angles from the neutral positions and scale it by \(r\). The final joint positions are therefore \(\Psi_0 + r(\Psi - \Psi_0)\). We will use the beginnings of the preprogrammed steps (right before the start of the swing) as the neutral positions.

Let’s repeat the previous exercise, but gradually ramp up the amplitude from 0 to 1:

theta_ts = np.linspace(0, 3 * 2 * np.pi, 10000)
r_ts = np.linspace(0, 1, 10000)

##### THIS SECTION HAS CHANGED #####
joint_angles_by_leg = {}
for leg, psi_func in psi_funcs.items():
    neutral_pos = psi_func(0)[:, np.newaxis]
    joint_angles_by_leg[leg] = neutral_pos + r_ts * (psi_func(theta_ts) - neutral_pos)
####################################

fig, axs = plt.subplots(3, 2, figsize=(7, 5), sharex=True, sharey=True)
for i_side, side in enumerate("LR"):
    for i_pos, pos in enumerate("FMH"):
        leg = f"{side}{pos}"
        ax = axs[i_pos, i_side]
        psi_func = psi_funcs[leg]
        joint_angles = np.rad2deg(joint_angles_by_leg[leg])
        for i_dof, dof_name in enumerate(dofs_per_leg):
            legend = dof_name if i_pos == 0 and i_side == 0 else None
            ax.plot(theta_ts, joint_angles[i_dof, :], linewidth=1, label=legend)
        if i_pos == 2:
            ax.set_xlabel("Phase")
            ax.set_xticks(np.pi * np.arange(7))
            ax.set_xticklabels(["0" if x == 0 else rf"{x}$\pi$" for x in np.arange(7)])
        if i_side == 0:
            ax.set_ylabel(r"DoF angle ($\degree$)")
        ax.set_title(f"{leg} leg")
        ax.set_ylim(-180, 180)
        ax.set_yticks([-180, -90, 0, 90, 180])
fig.legend(loc=7)
fig.tight_layout()
fig.subplots_adjust(right=0.8)
fig.savefig(output_dir / "three_steps_amp_modulated.png")

We have now built the individual elements of the controller:

• On the level of inter-leg coordination, the CPG network controls the phase \(\theta\) of each leg of the magnitude \(r\) of its steps.

• On the level of per-leg kinematics, we find the corresponding joint states at the phase \(\theta\) based on experimentally recorded data, scaled by the amplitude \(r\).

In the next section, we will piece these components together and plug them into the physics simulation.

Plugging the controller into the simulation

We will now put everything together and control the simulated fly using the controller that we have designed. The following content assumes that you have read the tutorial “Interacting with NeuroMechFly”.

We start by initializing the simulation:

from flygym import Fly, Camera, SingleFlySimulation
from flygym.preprogrammed import all_leg_dofs

run_time = 1
fly = Fly(init_pose="stretch", actuated_joints=all_leg_dofs, control="position")
cam = Camera(fly=fly, play_speed=0.1, draw_contacts=False)
sim = SingleFlySimulation(fly=fly, cameras=[cam], timestep=1e-4)

We will also initialize a CPG network:

intrinsic_freqs = np.ones(6) * 12
intrinsic_amps = np.ones(6) * 1
phase_biases = np.pi * np.array(
    [
        [0, 1, 0, 1, 0, 1],
        [1, 0, 1, 0, 1, 0],
        [0, 1, 0, 1, 0, 1],
        [1, 0, 1, 0, 1, 0],
        [0, 1, 0, 1, 0, 1],
        [1, 0, 1, 0, 1, 0],
    ]
)
coupling_weights = (phase_biases > 0) * 10
convergence_coefs = np.ones(6) * 1000

cpg_network = CPGNetwork(
    timestep=1e-4,
    intrinsic_freqs=intrinsic_freqs,
    intrinsic_amps=intrinsic_amps,
    coupling_weights=coupling_weights,
    phase_biases=phase_biases,
    convergence_coefs=convergence_coefs,
    seed=1,
)

swing_start = np.empty(6)
swing_end = np.empty(6)
for i, leg in enumerate(legs):
    swing_start[i] = single_steps_data["swing_stance_time"]["swing"][leg]
    swing_end[i] = single_steps_data["swing_stance_time"]["stance"][leg]
swing_start /= preprogrammed_steps_length * preprogrammed_steps_timestep
swing_start *= 2 * np.pi
swing_end /= preprogrammed_steps_length * preprogrammed_steps_timestep
swing_end *= 2 * np.pi

# have the rest phase in between the swing and stance phase (as the data starts with swing initiation)
psi_rest_phases = np.ones_like(swing_start)
for i, leg in enumerate(legs):
    psi_rest_phases[i] = (swing_end[i] + 2 * np.pi) / 2

Let’s run the simulation:

from tqdm import trange

obs, info = sim.reset()
for _ in trange(int(run_time / sim.timestep)):
    cpg_network.step()
    joints_angles = {}
    for i, leg in enumerate(legs):
        psi = psi_funcs[leg](cpg_network.curr_phases[i])
        psi_base = psi_funcs[leg](psi_rest_phases[i])
        adjusted_psi = psi_base + (psi - psi_base) * cpg_network.curr_magnitudes[i]
        for dof, angle in zip(dofs_per_leg, adjusted_psi):
            joints_angles[f"joint_{leg}{dof}"] = angle
    action = {"joints": np.array([joints_angles[dof] for dof in fly.actuated_joints])}
    obs, reward, terminated, truncated, info = sim.step(action)
    sim.render()

cam.save_video(output_dir / "cpg_controller.mp4", 0)

100%|██████████| 10000/10000 [00:18<00:00, 549.74it/s]

Leg adhesion

Insects, including flies, have evolved highly specialized adhesive structures to facilitate locomotion over complex 3D terrain. Substantial normal forces (10–100 times body weight) and frictional forces emerge from interactions between the adhesive pads and underlying substrates. These allow insects to navigate 3D terrain with ease. Because we cannot fully represent the physics underlying real, biological adhesion, we added a more abstract leg adhesion to our model by injecting an additional normal force to the pretarsus of each leg when it is in contact with a substrate. This adhesive force increases the normal force toward the object and the frictional force.

Despite the huge forces generated by adhesive pads, insects can still lift their legs, seemingly with out effort. The mechanisms for lifting off are not well understood in Drosophila. Therefore, we abstracted the mechanisms used by other insects for lifting by turning adhesion forces on during stance and off during swing phases. In the preprogrammed stepping data, we have also indicated the start (in seconds) of the swing and stance periods:

single_steps_data["swing_stance_time"]

{'swing': {'RF': 0.0, 'RM': 0.0, 'RH': 0.0, 'LF': 0.0, 'LM': 0.0, 'LH': 0.0},
 'stance': {'RF': 0.051000000000000004,
  'RM': 0.048,
  'RH': 0.042,
  'LF': 0.051000000000000004,
  'LM': 0.048,
  'LH': 0.042}}

Let’s write a function that, given the phases of the legs, return a boolean mask indicating whether adhesion should be on (during stance) or off (during swing):

def get_adhesion_onoff(theta):
    theta = theta % (2 * np.pi)
    return ~((theta > swing_start) & (theta < swing_end)).squeeze()

To illustrate this binary signal (low = off, during swing; high = on, during stance):

onoff_signal = np.zeros((6, phase_grid.size), dtype=bool)
for i in range(phase_grid.size):
    onoff_signal[:, i] = get_adhesion_onoff(phase_grid[i])

fig, ax = plt.subplots(figsize=(4, 2), tight_layout=True)
for i in range(6):
    ax.plot(phase_grid, onoff_signal[i, :] - i * 1.5)
ax.set_yticks(-np.arange(6) * 1.5 + 0.5)
ax.set_yticklabels(legs)
ax.set_xticks(np.arange(5) * np.pi / 2)
ax.set_xticklabels(["0", r"$\pi/2$", r"$\pi$", r"3$\pi$/2", r"$2\pi$"])
ax.set_xlabel("Phase")
ax.set_ylabel("Adhesion on/off")
fig.savefig(output_dir / "adhesion_signal.png")

We can rerun the NeuroMechFly simulation with adhesion enabled. The parts of the code that have been changed are indicated with comments.

run_time = 1

fly = Fly(
    init_pose="stretch",
    actuated_joints=all_leg_dofs,
    control="position",
    enable_adhesion=True,
    draw_adhesion=True,
)
cam = Camera(fly=fly, play_speed=0.1, draw_contacts=False)
sim = SingleFlySimulation(fly=fly, cameras=[cam], timestep=1e-4)

cpg_network.reset()

ang = []

obs, info = sim.reset()
for _ in trange(int(run_time / sim.timestep)):
    cpg_network.step()
    joints_angles = {}
    for i, leg in enumerate(legs):
        psi = psi_funcs[leg](cpg_network.curr_phases[i])
        psi_base = psi_funcs[leg](psi_rest_phases[i])
        adjusted_psi = psi_base + (psi - psi_base) * cpg_network.curr_magnitudes[i]
        for dof, angle in zip(dofs_per_leg, adjusted_psi):
            joints_angles[f"joint_{leg}{dof}"] = angle
    adhesion_onoff = get_adhesion_onoff(cpg_network.curr_phases)
    ang.append([joints_angles[dof] for dof in fly.actuated_joints])
    action = {
        "joints": np.array([joints_angles[dof] for dof in fly.actuated_joints]),
        ##### THIS LINE IS NEW #####
        "adhesion": adhesion_onoff.astype(int),
        ############################
    }
    obs, reward, terminated, truncated, info = sim.step(action)
    sim.render()

cam.save_video(output_dir / "cpg_controller_with_adhesion.mp4", 0)

100%|██████████| 10000/10000 [00:26<00:00, 384.60it/s]

In summary, in this tutorial we have (1) implemented a Python class for CPG networks, (2) used it to modulate the stepping of legs using experimentally recorded data, (3) plugged this controller into the NeuroMechFly embodiment, and (4) added leg adhesion to the simulation. Note that the controller we built here is feedforward — that is, mechanosensory feedback is not used by the controller (except the position feedback in the PD controller for individual joints). In the next tutorial, we will build a rule-based controller where leg coordination is accomplished using sensory feedback in a more distributed manner.