Controlling locomotion with CPGs#

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 (ie., 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 (ie. 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).

 1import numpy as np
 2
 3
 4def calculate_ddt(theta, r, w, phi, nu, R, alpha):
 5    """Given the current state variables theta, r and network parameters
 6    w, phi, nu, R, alpha, calculate the time derivatives of theta and r."""
 7    intrinsic_term = 2 * np.pi * nu
 8    phase_diff = theta[np.newaxis, :] - theta[:, np.newaxis]
 9    coupling_term = (r * w * np.sin(phase_diff - phi)).sum(axis=1)
10    dtheta_dt = intrinsic_term + coupling_term
11    dr_dt = alpha * (R - r)
12    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.

 1# Fix random seed for reproducibility
 2np.random.seed(0)
 3
 4
 5class CPGNetwork:
 6    def __init__(
 7        self,
 8        timestep,
 9        intrinsic_freqs,
10        intrinsic_amps,
11        coupling_weights,
12        phase_biases,
13        convergence_coefs,
14        init_phases=None,
15        init_magnitudes=None,
16    ) -> None:
17        """Initialize a CPG network consisting of N oscillators.
18
19        Parameters
20        ----------
21        timestep : float
22            The timestep of the simulation.
23        intrinsic_frequencies : np.ndarray
24            The intrinsic frequencies of the oscillators, shape (N,).
25        intrinsic_amplitudes : np.ndarray
26            The intrinsic amplitude of the oscillators, shape (N,).
27        coupling_weights : np.ndarray
28            The coupling weights between the oscillators, shape (N, N).
29        phase_biases : np.ndarray
30            The phase biases between the oscillators, shape (N, N).
31        convergence_coefs : np.ndarray
32            Coefficients describing the rate of convergence to oscillator
33            intrinsic amplitudes, shape (N,).
34        init_phases : np.ndarray, optional
35            Initial phases of the oscillators, shape (N,). The phases are
36            randomly initialized if not provided.
37        init_magnitudes : np.ndarray, optional
38            Initial magnitudes of the oscillators, shape (N,). The
39            magnitudes are randomly initialized if not provided.
40        """
41        self.timestep = timestep
42        self.num_cpgs = intrinsic_freqs.size
43        self.intrinsic_freqs = intrinsic_freqs
44        self.intrinsic_amps = intrinsic_amps
45        self.coupling_weights = coupling_weights
46        self.phase_biases = phase_biases
47        self.convergence_coefs = convergence_coefs
48
49        self.reset(init_phases, init_magnitudes)
50
51        # Check if the parameters have the right shape
52        assert intrinsic_freqs.shape == (self.num_cpgs,)
53        assert coupling_weights.shape == (self.num_cpgs, self.num_cpgs)
54        assert phase_biases.shape == (self.num_cpgs, self.num_cpgs)
55        assert convergence_coefs.shape == (self.num_cpgs,)
56        assert self.curr_phases.shape == (self.num_cpgs,)
57        assert self.curr_magnitudes.shape == (self.num_cpgs,)
58
59    def step(self):
60        """Integrate the ODEs using Euler's method."""
61        dtheta_dt, dr_dt = calculate_ddt(
62            theta=self.curr_phases,
63            r=self.curr_magnitudes,
64            w=self.coupling_weights,
65            phi=self.phase_biases,
66            nu=self.intrinsic_freqs,
67            R=self.intrinsic_amps,
68            alpha=self.convergence_coefs,
69        )
70        self.curr_phases += dtheta_dt * self.timestep
71        self.curr_magnitudes += dr_dt * self.timestep
72
73    def reset(self, init_phases=None, init_magnitudes=None):
74        if init_phases is None:
75            self.curr_phases = np.random.random(self.num_cpgs) * 2 * np.pi
76        else:
77            self.curr_phases = init_phases
78
79        if init_magnitudes is None:
80            self.curr_magnitudes = np.random.random(self.num_cpgs) * self.intrinsic_amps
81        else:
82            self.curr_magnitudes = init_magnitudes

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

https://raw.githubusercontent.com/NeLy-EPFL/_media/main/flygym/simple_cpg.png

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.

 1intrinsic_freqs = np.ones(3)
 2intrinsic_amps = np.array([1.0, 1.1, 1.2])
 3coupling_weights = np.array(
 4    [
 5        [0, 1, 0],
 6        [1, 0, 1],
 7        [0, 1, 0],
 8    ]
 9)
10phase_biases = np.deg2rad(
11    np.array(
12        [
13            [0, 120, 0],
14            [-120, 0, 120],
15            [0, -120, 0],
16        ]
17    )
18)
19convergence_coefs = np.ones(3)
20
21network = CPGNetwork(
22    timestep=1e-3,
23    intrinsic_freqs=intrinsic_freqs,
24    intrinsic_amps=intrinsic_amps,
25    coupling_weights=coupling_weights,
26    phase_biases=phase_biases,
27    convergence_coefs=convergence_coefs,
28)
29
30num_steps = int(10 / network.timestep)
31phase_hist = np.empty((num_steps, 3))
32magnitude_hist = np.empty((num_steps, 3))
33
34# Simulate the network
35for i in range(num_steps):
36    network.step()
37    phase_hist[i, :] = network.curr_phases
38    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.

 1import matplotlib.pyplot as plt
 2
 3
 4fig, axs = plt.subplots(2, 1, figsize=(5, 5), sharex=True)
 5t = np.arange(num_steps) * network.timestep
 6axs[0].plot(t, phase_hist % (2 * np.pi), linewidth=1)
 7axs[0].set_yticks([0, np.pi, 2 * np.pi])
 8axs[0].set_yticklabels(["0", "$\pi$", "$2\pi$"])
 9axs[0].set_ylabel("Phase")
10axs[1].plot(t, magnitude_hist, linewidth=1)
11axs[1].set_ylabel("Magnitude")
12axs[1].set_xlabel("Time (s)")
13fig.savefig("./outputs/simple_cpg_rollout.png")
https://raw.githubusercontent.com/NeLy-EPFL/_media/main/flygym/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 (ie. 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 (ie. how far into the step the leg is), while the magnitude of the CPG represents the magnitude of the step (ie. 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 2.0 paper.

https://raw.githubusercontent.com/NeLy-EPFL/_media/main/flygym/tripod_cpg.png

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

 1intrinsic_freqs = np.ones(6) * 12
 2intrinsic_amps = np.ones(6) * 1
 3phase_biases = np.pi * np.array(
 4    [
 5        [0, 1, 0, 1, 0, 1],
 6        [1, 0, 1, 0, 1, 0],
 7        [0, 1, 0, 1, 0, 1],
 8        [1, 0, 1, 0, 1, 0],
 9        [0, 1, 0, 1, 0, 1],
10        [1, 0, 1, 0, 1, 0],
11    ]
12)
13coupling_weights = (phase_biases > 0) * 10
14convergence_coefs = np.ones(6) * 20
15
16network = CPGNetwork(
17    timestep=1e-4,
18    intrinsic_freqs=intrinsic_freqs,
19    intrinsic_amps=intrinsic_amps,
20    coupling_weights=coupling_weights,
21    phase_biases=phase_biases,
22    convergence_coefs=convergence_coefs,
23)
24
25# Simulate the network
26num_steps = int(1 / network.timestep)
27phase_hist = np.empty((num_steps, 6))
28magnitude_hist = np.empty((num_steps, 6))
29for i in range(num_steps):
30    network.step()
31    phase_hist[i, :] = network.curr_phases
32    magnitude_hist[i, :] = network.curr_magnitudes
33
34# Visualize
35fig, axs = plt.subplots(2, 1, figsize=(5, 5), sharex=True)
36t = np.arange(num_steps) * network.timestep
37axs[0].plot(t, phase_hist % (2 * np.pi), linewidth=1)
38axs[0].set_yticks([0, np.pi, 2 * np.pi])
39axs[0].set_yticklabels(["0", "$\pi$", "$2\pi$"])
40axs[0].set_ylabel("Phase")
41axs[1].plot(t, magnitude_hist, linewidth=1)
42axs[1].set_ylabel("Magnitude")
43axs[1].set_xlabel("Time (s)")
44fig.savefig("./outputs/tripod_cpg_rollout.png")
https://raw.githubusercontent.com/NeLy-EPFL/_media/main/flygym/tripod_cpg_rollout.png

Now, let’s load the behavior kinematics data:

1import pickle
2from flygym.common import get_data_path
3
4
5single_steps_path = get_data_path("flygym", "data") / "behavior/single_steps.pkl"
6with open(single_steps_path, "rb") as f:
7    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:

 1preprogrammed_steps_length = len(single_steps_data["joint_LFCoxa"])
 2preprogrammed_steps_timestep = single_steps_data["meta"]["timestep"]
 3print(
 4    f"Preprogrammed steps have a length of {preprogrammed_steps_length} steps "
 5    f"at dt={preprogrammed_steps_timestep}s."
 6)
 7for k, v in single_steps_data.items():
 8    if k.startswith("joint_"):
 9        assert len(v) == preprogrammed_steps_length
10        assert v[0] == v[-1]
Preprogrammed steps have a length of 1278 steps at dt=0.0001s.

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)\):

 1from scipy.interpolate import CubicSpline
 2
 3legs = [f"{side}{pos}" for side in "LR" for pos in "FMH"]
 4dofs_per_leg = [
 5    "Coxa",
 6    "Coxa_roll",
 7    "Coxa_yaw",
 8    "Femur",
 9    "Femur_roll",
10    "Tibia",
11    "Tarsus1",
12]
13phase_grid = np.linspace(0, 2 * np.pi, preprogrammed_steps_length)
14psi_funcs = {}
15for leg in legs:
16    joint_angles = np.array(
17        [single_steps_data[f"joint_{leg}{dof}"] for dof in dofs_per_leg]
18    )
19    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:

 1theta_ts = np.linspace(0, 3 * 2 * np.pi, 10000)
 2
 3joint_angles_by_leg = {}
 4for leg, psi_func in psi_funcs.items():
 5    joint_angles_by_leg[leg] = psi_func(theta_ts)
 6
 7fig, axs = plt.subplots(3, 2, figsize=(7, 5), sharex=True, sharey=True)
 8for i_side, side in enumerate("LR"):
 9    for i_pos, pos in enumerate("FMH"):
10        leg = f"{side}{pos}"
11        ax = axs[i_pos, i_side]
12        psi_func = psi_funcs[leg]
13        joint_angles = np.rad2deg(joint_angles_by_leg[leg])
14        for i_dof, dof_name in enumerate(dofs_per_leg):
15            legend = dof_name if i_pos == 0 and i_side == 0 else None
16            ax.plot(theta_ts, joint_angles[i_dof, :], linewidth=1, label=legend)
17        if i_pos == 2:
18            ax.set_xlabel("Phase")
19            ax.set_xticks(np.pi * np.arange(7))
20            ax.set_xticklabels(["0" if x == 0 else f"{x}$\pi$" for x in np.arange(7)])
21        if i_side == 0:
22            ax.set_ylabel(r"DoF angle ($\degree$)")
23        ax.set_title(f"{leg} leg")
24        ax.set_ylim(-180, 180)
25        ax.set_yticks([-180, -90, 0, 90, 180])
26fig.legend(loc=7)
27fig.tight_layout()
28fig.subplots_adjust(right=0.8)
29fig.savefig("./outputs/three_steps_phase_only.png")
https://raw.githubusercontent.com/NeLy-EPFL/_media/main/flygym/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 preprogramed 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:

 1theta_ts = np.linspace(0, 3 * 2 * np.pi, 10000)
 2r_ts = np.linspace(0, 1, 10000)
 3
 4##### THIS SECTION HAS CHANGED #####
 5joint_angles_by_leg = {}
 6for leg, psi_func in psi_funcs.items():
 7    neutral_pos = psi_func(0)[:, np.newaxis]
 8    joint_angles_by_leg[leg] = neutral_pos + r_ts * (psi_func(theta_ts) - neutral_pos)
 9####################################
10
11fig, axs = plt.subplots(3, 2, figsize=(7, 5), sharex=True, sharey=True)
12for i_side, side in enumerate("LR"):
13    for i_pos, pos in enumerate("FMH"):
14        leg = f"{side}{pos}"
15        ax = axs[i_pos, i_side]
16        psi_func = psi_funcs[leg]
17        joint_angles = np.rad2deg(joint_angles_by_leg[leg])
18        for i_dof, dof_name in enumerate(dofs_per_leg):
19            legend = dof_name if i_pos == 0 and i_side == 0 else None
20            ax.plot(theta_ts, joint_angles[i_dof, :], linewidth=1, label=legend)
21        if i_pos == 2:
22            ax.set_xlabel("Phase")
23            ax.set_xticks(np.pi * np.arange(7))
24            ax.set_xticklabels(["0" if x == 0 else f"{x}$\pi$" for x in np.arange(7)])
25        if i_side == 0:
26            ax.set_ylabel(r"DoF angle ($\degree$)")
27        ax.set_title(f"{leg} leg")
28        ax.set_ylim(-180, 180)
29        ax.set_yticks([-180, -90, 0, 90, 180])
30fig.legend(loc=7)
31fig.tight_layout()
32fig.subplots_adjust(right=0.8)
33fig.savefig("./outputs/three_steps_amp_modulated.png")
https://raw.githubusercontent.com/NeLy-EPFL/_media/main/flygym/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:

 1import flygym.mujoco
 2import flygym.mujoco.preprogrammed
 3
 4run_time = 1
 5sim_params = flygym.mujoco.Parameters(
 6    timestep=1e-4, render_mode="saved", render_playspeed=0.1, draw_contacts=True
 7)
 8nmf = flygym.mujoco.NeuroMechFly(
 9    sim_params=sim_params,
10    init_pose="stretch",
11    actuated_joints=flygym.mujoco.preprogrammed.all_leg_dofs,
12    control="position",
13)

We will also initialize a CPG network:

 1intrinsic_freqs = np.ones(6) * 12
 2intrinsic_amps = np.ones(6) * 1
 3phase_biases = np.pi * np.array(
 4    [
 5        [0, 1, 0, 1, 0, 1],
 6        [1, 0, 1, 0, 1, 0],
 7        [0, 1, 0, 1, 0, 1],
 8        [1, 0, 1, 0, 1, 0],
 9        [0, 1, 0, 1, 0, 1],
10        [1, 0, 1, 0, 1, 0],
11    ]
12)
13coupling_weights = (phase_biases > 0) * 10
14convergence_coefs = np.ones(6) * 20
15
16cpg_network = CPGNetwork(
17    timestep=1e-4,
18    intrinsic_freqs=intrinsic_freqs,
19    intrinsic_amps=intrinsic_amps,
20    coupling_weights=coupling_weights,
21    phase_biases=phase_biases,
22    convergence_coefs=convergence_coefs,
23)

Let’s run the simulation:

 1from tqdm import trange
 2
 3obs, info = nmf.reset()
 4for i in trange(int(run_time / sim_params.timestep)):
 5    cpg_network.step()
 6    joints_angles = {}
 7    for i, leg in enumerate(legs):
 8        psi = psi_funcs[leg](cpg_network.curr_phases[i])
 9        psi_0 = psi_funcs[leg](0)
10        adjusted_psi = psi_0 + cpg_network.curr_magnitudes[i] * (psi - psi_0)
11        for dof, angle in zip(dofs_per_leg, adjusted_psi):
12            joints_angles[f"joint_{leg}{dof}"] = angle
13    action = {"joints": np.array([joints_angles[dof] for dof in nmf.actuated_joints])}
14    obs, reward, terminated, truncated, info = nmf.step(action)
15    nmf.render()
16
17nmf.save_video("./outputs/cpg_controller.mp4")
100%|██████████| 10000/10000 [00:18<00:00, 534.28it/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:

1single_steps_data["swing_stance_time"]
{'swing': {'RF': 0.0098,
  'LF': 0.0098,
  'RM': 0.0012000000000000001,
  'LM': 0.0012000000000000001,
  'RH': 0.0012000000000000001,
  'LH': 0.0012000000000000001},
 'stance': {'RF': 0.0408,
  'LF': 0.0408,
  'RM': 0.0318,
  'LM': 0.0318,
  'RH': 0.027200000000000002,
  'LH': 0.027200000000000002}}

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):

 1swing_start = np.empty(6)
 2swing_end = np.empty(6)
 3for i, leg in enumerate(legs):
 4    swing_start[i] = single_steps_data["swing_stance_time"]["swing"][leg]
 5    swing_end[i] = single_steps_data["swing_stance_time"]["stance"][leg]
 6swing_start /= preprogrammed_steps_length * preprogrammed_steps_timestep
 7swing_start *= 2 * np.pi
 8swing_end /= preprogrammed_steps_length * preprogrammed_steps_timestep
 9swing_end *= 2 * np.pi
10
11
12def get_adhesion_onoff(theta):
13    theta = theta % (2 * np.pi)
14    return ~((theta > swing_start) & (theta < swing_end)).squeeze()

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

 1onoff_signal = np.zeros((6, phase_grid.size), dtype=bool)
 2for i in range(phase_grid.size):
 3    onoff_signal[:, i] = get_adhesion_onoff(phase_grid[i])
 4
 5fig, ax = plt.subplots(figsize=(4, 2))
 6for i in range(6):
 7    ax.plot(phase_grid, onoff_signal[i, :] - i * 1.5)
 8ax.set_yticks(-np.arange(6) * 1.5 + 0.5)
 9ax.set_yticklabels(legs)
10ax.set_xticks(np.arange(5) * np.pi / 2)
11ax.set_xticklabels(["0", "$\pi/2$", "$\pi$", "3$\pi$/2", "$2\pi$"])
12ax.set_xlabel("Phase")
13ax.set_ylabel("Adhesion on/off")
14fig.savefig("./outputs/adhesion_signal.png")
https://raw.githubusercontent.com/NeLy-EPFL/_media/main/flygym/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.

 1run_time = 1
 2sim_params = flygym.mujoco.Parameters(
 3    timestep=1e-4,
 4    render_mode="saved",
 5    render_playspeed=0.1,
 6    enable_adhesion=True,  # THIS HAS CHANGED
 7    draw_adhesion=True,  # THIS HAS CHANGED (tarsus color indicates adhesion on/off)
 8)
 9nmf = flygym.mujoco.NeuroMechFly(
10    sim_params=sim_params,
11    init_pose="stretch",
12    actuated_joints=flygym.mujoco.preprogrammed.all_leg_dofs,
13    control="position",
14)
15
16cpg_network.reset()
17
18obs, info = nmf.reset()
19for i in trange(int(run_time / sim_params.timestep)):
20    cpg_network.step()
21    joints_angles = {}
22    for i, leg in enumerate(legs):
23        psi = psi_funcs[leg](cpg_network.curr_phases[i])
24        psi_0 = psi_funcs[leg](0)
25        adjusted_psi = psi_0 + cpg_network.curr_magnitudes[i] * (psi - psi_0)
26        for dof, angle in zip(dofs_per_leg, adjusted_psi):
27            joints_angles[f"joint_{leg}{dof}"] = angle
28    adhesion_onoff = get_adhesion_onoff(cpg_network.curr_phases)
29    action = {
30        "joints": np.array([joints_angles[dof] for dof in nmf.actuated_joints]),
31        ##### THIS LINE IS NEW #####
32        "adhesion": adhesion_onoff.astype(int),
33        ############################
34    }
35    obs, reward, terminated, truncated, info = nmf.step(action)
36    nmf.render()
37
38nmf.save_video("./outputs/cpg_controller_with_adhesion.mp4")
100%|██████████| 10000/10000 [00:38<00:00, 260.59it/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.