Network Architecture Producing Swing to Stance Transitions in an Insect Walking System

The walking system of the stick insect is one of the most thoroughly described invertebrate systems. We know a lot about the role of sensory input in the control of stepping of a single leg. However, the neuronal organization and connectivity of the central neural networks underlying the rhythmic activation and coordination of leg muscles still remain elusive. It is assumed that these networks can couple in the absence of phasic sensory input due to the observation of spontaneous recurrent patterns (SRPs) of coordinated motor activity equivalent to fictive stepping-phase transitions. Here we sought to quantify the phase of motor activity within SRPs in the isolated and interconnected meso- and meta-thoracic ganglia. We show that SRPs occur not only in the meso-, but also in the metathoracic ganglia of the stick insect, discovering a qualitative difference between them. We construct a network based on neurophysiological data capable of reproducing the measured SRP phases to investigate this difference. By comparing network output to the biological measurements we confirm the plausibility of the architecture and provide a hypothesis to account for these qualitative differences. The neural architecture we present couples individual central pattern generators to reproduce the fictive stepping-phase transitions observed in deafferented stick insect preparations after pharmacological activation, providing insights into the neural architecture underlying coordinated locomotion.

1. Introduction

Walking is based on cyclic patterns of coordinated movement among legs. A step consists of the stance phase, during which the leg has ground contact, supports, and propels the animal, and the swing phase, when the leg is lifted and moves back to its initial position to complete the stepping cycle. We know a lot about the role of sensory organs in initiating or terminating the stepping phases in a single leg (1). However, the exact organization of the central neural networks underlying stepping-phase transitions in multi-segmented locomotor organs still remains elusive.

So far there have been different approaches considered in modeling studies focusing on either vertebrates or invertebrates. According to one approach, stepping is based on motor neuron synergies controlled by central rhythm and pattern generating networks, which then interact with sensory feedback to generate stepping (2, 3). However, in other studies sensory feedback is a larger focus, Ekeberg et al. (4) reconstruct stepping in a cat's hindlimb based on sensory signals and the mechanical coupling of the legs. Similarly, in models studying the stick insect walking system, sensory signals play a large role as the position of the leg is the main factor for ending a step phase and initiating the next one (5, 6). Lastly, there are reduced models which define a single network for each leg and do not consider the individual networks controlling each leg segment (5, 7).

For our modeling approach we have decided to use biological data from the stick insect because it is known that there are dedicated central pattern generators (CPGs) to move each of the main leg joints (8). Walking in the stick insect is based on intra- and intersegmental interaction among CPG networks and peripheral sensory input (9). The relative contribution of the central vs. peripheral neuronal mechanisms for walking has been often debated. Pharmacological activation of the central networks in deafferented insect preparations has contributed significantly to our understanding regarding the possible role of sensory input in adapting the default centrally-generated patterns, potentially giving rise to behaviorally-relevant coordination. Such experiments have not resulted in generation of motor patterns similar to the coordination patterns observed in vivo. Thus, there has been no indication of fictive walking in adult insect preparations (10–12). However, previous experiments in the deafferented stick insect mesothoracic ganglion, after application of the muscarinic agonist pilocarpine, have indeed revealed coordinated activity of the motor neuron pools within a hemisegment in the absence of phasic sensory input. This activity resembles stepping-phase transitions in vivo (8) and are therefore referred to as fictive stepping-phase transitions. There have been distinct patterns of coordinated motor activity described which spontaneously emerge in a non-cycle-to-cycle fashion throughout the recording, and are therefore called spontaneous recurrent patterns (SRPs). The SRP1 is characterized by a switch from protractor to retractor motor neuron activity during a depressor burst consisting of both the slow and fast depressor motor neurons. This switch resembles a transition from swing to stance stepping phase during forward walking. The SRP2 is characterized by the opposite switch from retractor to protractor activity during a combined slow and fast depressor burst, resembling a transition from swing to stance stepping phase during backward walking (8, 13). These patterns have only been described for the mesothoracic ganglion of the stick insect so far, with the SRP1 showing a higher occurrence throughout the recording compared to SRP2. Although there have been neurons identified in the ventral nerve cord which influence the occurrence of SRPs, the underlying mechanisms resulting in SRP generation are still unknown (12, 13).

SRP-like activity patterns have also been reported in neurophysiological recordings of motor activity in the fourth thoracic ganglion of the crayfish after bath application of cholinergic agonists (14). In their study, 90% of the recorded intervals showed coordinated activity consisting of the levator motor neurons bursting in phase with the retractor motor neurons. This pattern is indicative of fictive backward walking. In the stick insect mesothoracic ganglion, according to the only study focusing on intrasegmental coordination in insects, 78% of the patterns of coordinated activity belonged to the SRP1 type, representing a fictive transition from swing to stance phase during forward stepping (8). Here, we hypothesize that SRPs can also be found in other ganglia, such as the metathoracic ganglion, and that an underlying network architecture exists to generate these SRPs in the absence of sensory feedback. We have chosen to investigate potential architectures through simulation, following the sentiment by (15) that “… many interactive processes on the subcellular, cellular and network levels are dynamic and complex. Computational methods are therefore required to test whether tentative explanations derived by intuition can account for experimental findings.”

The majority of current engineering research develops coordination through feedback and subsequent rules based on the returned values. Popular controllers such as WalkNet (5) and neuroWalknet (16) use behavioral data from biology to create an artificial neural network capable of controlling 18 degrees of freedom. The individual leg controllers use coordination rules, network architecture, and sensory feedback to produce walking behaviors on a six-legged robot. WalkNet can perform forwards and backwards walking on smooth and uneven terrain as well as curved walking. neuroWalknet extends the walking ability to produce different footfall patterns including tripod, tetrapod, and pentapod patterns as well as other stable intermediate patterns observed in stick insects (17) and Drosophila (18). Ekeberg et al. (19) realize a neuronal control network sufficient for controlling stick insect legs in each segment—front, middle, and hind. They use a simulation to send feedback to the control network as the primary driver for joint coordination. von Twickel et al. (20) use a similar strategy by decentralizing the joint controllers and relying on feedback from a simulated leg to develop a single-leg controller. All of the mentioned studies have proposed different control models for joint (intralimb) and leg (interlimb) coordination to reproduce biological behaviors and confirm the importance of sensory feedback in coordinated walking. However, there have not been any studies investigating the underlying neural network architecture before sensory feedback is applied in order to reproduce the fictive stepping-phase transitions observed in deafferented stick insect preparations. We address this by using biologically-plausible spiking and non-spiking neuron models to create a single joint architecture similar to the one suggested by (21). The setup connects spiking rhythm-generating populations (RGPs) to non-spiking interneurons (NSIs) to inhibit the spiking motor neuron populations (MNPs). Coordination NSIs (cNSIs) are also used to connect the individual joints to create inter-joint coordination similar to the biological interneurons found by (13). The addition of NSIs is able to decouple the neuronal dynamics of the spiking populations, removing the issue of competing dynamics when communicating between spiking neurons. In this way, we find that the integration of NSI's into spiking neural networks reduces complexity when building the network while maintaining biological plausibility.

Our study shows that SRPs occur in both the meso- and meta-thoracic ganglia of the stick insect and highlight an interesting difference between the two thoracic ganglia. We propose a neural architecture that is able to reproduce the prominent SRP in each ganglia confirming that network architecture can account for the difference in SRP type.

2. Methods

2.1. Neurophysiology

2.1.1. Experimental Animals

The animals used in this study are adult female Indian stick insects of the species Carausius morosus bred in our colony at the Biocenter, University of Cologne. Animals are kept at 20–26°C with 45–60% humidity, under a 12 h light/12 h dark cycle. The experimental procedures described below comply with the German National and State Regulations for Animal Welfare and Animal Experiments.

2.1.2. Preparation and Experimental Setup

The experimental procedure followed in this study has been previously established by (8). Some experiments are performed on the ipsilateral hemisegment of the isolated (disconnected from neighboring ganglia) meso- and meta-thoracic ganglia. Other experiments are performed on the interconnected meso- and meta-thoracic ganglia. In all studies, the ganglia are deafferented and isolated from the rest of the nerve cord. Additionally, the first abdominal ganglion is always left interconnected to the meta-thoracic ganglion. Rhythmic activity in motor neuron pools is induced by bath application of the muscarinic acetylcholine receptor agonist pilocarpine, and is assessed by recording extracellular activity by placing extracellular electrodes on the lateral nerves nl2, nl3, nl5, and C2. This captures the motor activity innervating the leg muscles. All lateral nerves in the ganglia of interest are either crushed or cut to block afferent input. The nerves nl2 and nl5 carry the axons that innervate the antagonistic protractor and retractor coxae muscles, while nl3 and C2 innervate the extensor tibiae and depressor trocanteris, respectively (22). The protractor and retractor move the leg forwards and backwards along the horizontal plane, the extensor muscle extends the tibia and the depressor muscle allows for downward movement of the leg along the vertical plane. Collectively, these muscles allow the movement about the three main leg joints: the Thorax-Coxa (ThC), the Coxa-Trochanter (CTr), and the Femur-Tibia (FTi) joint.

The signal recorded with the extracellular electrodes is pre-amplified by 100-fold using isolated low-noise preamplifiers (model PA101, Electronics workshop, Zoological Institute, Cologne). It is further amplified by ten-fold to reach an overall gain of 1,000 and filtered (low-cut: 200 Hz, high-cut: 3 kHz) using a standard four-channel amplifier/signal conditioner (model MA102, Electronics workshop, Zoological Institute, Cologne). Finally, the signal is digitized at a sampling rate of 12 kHz, using the Micro 1401-3 acquisition unit (CED, Cambridge, UK) and it is monitored using the Spike2 software (CED, Cambridge, UK).

2.1.3. Analysis of the Neurophysiological Data

The electrophysiological data are initially assessed by using tools provided by Spike2. Time series of the action potentials (spikes) are marked by manually setting a threshold. In cases where spikes of non-interesting neurons also cross the threshold, such as those of the common inhibitory motor neuron, the respective time-series are subtracted from the data. The depressor bursts are marked, noting the timing of the burst onset because the depressor cycles are used as a reference for the analysis.

Figure 1 depicts how the phase difference between nerves is calculated. In this study, “Nerve 2” in Figure 1 represents bursts of motor activity in the depressor because it is the reference nerve. Spontaneous recurrent patterns of coordinated motor activity are manually marked throughout the recording to use for the phase calculation. The spike time-series are exported using a sampling rate of 1,000 Hz and the phase of each spike within the ongoing depressor cycle is calculated throughout the recording using MATLAB. Phase values are binned (bin size = 10°) and the mean bin value [± standard deviation (STD)] among animal preparations is plotted after normalizing data of each animal to the maximum number of spikes. Histograms are made by either considering all depressor cycles throughout the recording or selecting only those cycles during which an SRP1 or SRP2 occurred. Circular means and angular deviations, as well as the 90% confidence intervals (CI) are calculated using the Circular Statistics Toolbox in MATLAB (23).

FIGURE 1

Figure 1. Illustration of how the phase (ϕ_x) of each spike of a nerve (Nerve 1) is calculated in relation to the cycle of a reference nerve (Nerve 2). The time is recorded for all spikes reaching the set threshold within Nerve 1 (t_x). The phase is calculated for each spike and then averaged to find the average activity of a nerve as compared to the reference nerve. The cycle period is dynamic so the phase is calculated based on the corresponding cycle, the cycle period is represented by t₂ − t₁ where t₁ is the cycle onset and t₂ is the end of the cycle.

2.2. Network Simulation

The neural network designed in this study is based on observations from biological research. Individual CPGs are known to control a single antagonistic muscle pair producing alternating rhythmic activity per joint (8, 9). This means that the individual networks must be coordinated to produce the observed swing and stance phases of a single leg step cycle. Membrane potential oscillations of NSIs have been shown to correlate with SRPs (13) indicating that NSIs are most likely involved in inter-joint coordination. Our network uses this knowledge to couple three CPG networks with NSIs to produce coordinated firing.

2.2.1. Single Joint Architecture

Before coordination can be achieved, a reliable anti-phasic output must be generated corresponding to the antagonistic muscle pairs in the stick insect leg controlling each joint. Each individual joint is controlled by one CPG consisting of two RGPs mutually inhibiting each other in a half-center oscillator architecture (9). Figure 2A shows a detailed diagram of the neural network controlling a single joint.

FIGURE 2

Figure 2. (A) Network architecture of an individual joint showing the half-center oscillator comprised of two RGPs mutually inhibiting each other with a synaptic weight of −5,000 nS. The RGPs excite buffer NSIs which then inhibit the MNPs periodically, creating output oscillations. The voltage threshold potential V_th of neurons in the RGPs is updated based on the input current to the velocity NSI (vNSI) in order to manipulate output frequency. Excitatory current is sent to the RGPs to initiate bursting and to the MNPs to promote tonic spiking. The output of the MNPs is sent to a motor controller. (B,C) The network architecture for coupling individual joints through coordination NSIs (cNSIs). The rhythm-generating populations are labeled based on the MNPs they control—Protractor (Pro), Retractor (Ret), Depressor (Dep), Levator (Lev), Extensor (Ext), and Flexor (Flx). The highlighted populations are the RGPs driving coordination. The burst duration NSI (bdNSI) is used during uncoordinated firing and silenced when the cNSIs are active. The retractor RGP sends spikes to the bdNSI. When the bdNSI is active, it waits a calculated amount of time before sending excitatory current to the protractor RGP and reducing the inhibition from the retractor to protractor RGP in order to stop the Pro MNP from bursting. (B) The architecture to produce SRP1 transitions. The Dep RGP drives coordination, sending excitatory spikes to the cNSIs to coordinate bursting through excitatory current to the Flx and Ret RGPs. (C) The architecture to produce SRP2 transitions. The Lev RGP drives coordination, sending excitatory spikes to the cNSIs to coordinate bursting through excitatory current to the Ext and Ret RGPs.

The RGPs and MNPs consist of adaptive exponential integrate-and-fire (AdEx) neurons. The parameters for the RGP neurons are set to bursting whereas the parameters used for the MNP neurons produce tonic spiking. These parameters are set according to (24). The AdEx neuron dynamics are described by Equations (1) and (2).

$\begin{array}{ll}\begin{array}{l}C\frac{d{V}_m}{dt}=-g_L\left(V_m-E_L\right)+g_L\left(\Delta T\right)e^{\frac{V_m-V_{th}}{\Delta T}}-w+I_e\hfill \end{array}& (1)\end{array}$

when V_m > 0mV then V_m → V_reset

$\begin{array}{ll}{\tau }_w\frac{dw}{dt}=a(V_m-E_L)-w& (2)\end{array}$

when V_m > 0mV then w → w + b

“where C is the membrane capacitance, V_m is the membrane potential, E_L is the resting potential, g_L is the leakage conductance, I_e is the bias current plus Gaussian white noise, a is the sub-threshold adaptation conductance, b is the spike-triggered adaptation, Δ_T is the sharpness factor, τ_w is the adaptation time constant, V_th is the voltage threshold potential, V_reset is the reset potential, and w is the spike adaptation current (24). Equation (1) defines the change in membrane potential per time step, whereas, Equation (2) outlines the current adaptation” (25).

All NSIs used within the network are modeled as leaky integrate-and-fire neurons with a high enough voltage threshold to avoid spiking. The neuronal dynamics are described by Equation (3).

$\begin{array}{ll}{\tau }_m\frac{d{V}_m}{dt}=-V_m+R{I}_{input}& (3)\end{array}$

“where τ_m = RC is the membrane time constant, V_m is membrane potential, I_input is input bias current plus Gaussian white noise, and R is membrane resistance” (25).

The RGPs are initialized with an excitatory current of 500 pA but this switches to a frequency regulating excitatory current. The current is designed to be low enough to induce bursting while producing the desired frequency. This current increases with frequency according to Equation (4).

$\begin{array}{l}{I}_{desc}=-21.3665·V_{th}^3-3555.1708·V_{th}^2-196971.3334·V_{th}\\ -3633546.7323& (4)\end{array}$

Where I_desc is the descending input current to the RGPs and V_th is the voltage threshold potential. The parameter calculations defined within the network architecture depend on frequency. In previous work, Strohmer et al. (26) found a linear relationship between V_th and frequency. Therefore, all linear relationship calculations use V_th as a substitute to approximate frequency instead of calculating frequency at each time step. This substitution reduces computational load.

The descending signals act as an “on/off” switch for the RGPs. Descending signals are also able to control the frequency of output oscillations from the RGPs by sending excitatory current to the velocity NSI (vNSI). This setup is based on biology, Berendes et al. (27) suggest that there could be multiple pathways for descending signals and that the signals controlling output frequency are connected to CPGs through sensory neurons. Validation that the vNSI is able to control output frequency by manipulating the V_th of the RGP neurons was previously shown by (25). The stick insect steps at a frequency of ~1–4 Hz on a slippery surface (28). These observed speeds are used to constrain the output frequency per joint during testing. Our work investigates the phase relationships observed in biological recordings of low-frequency motor activity in deafferented preparations to ultimately reproduce the phase relationships observed during live stepping. In this way, our aim is to match the phase relationship measured from deafferented samples and then visually confirm if this produces stepping-phase transitions in vivo. We can achieve this confirmation using a robot leg in simulation. We are able to directly compare network simulation and biological measurements by normalizing over degrees of a step cycle instead of timing, thereby, removing the problem of the difference in step cycle period.

In our network, the descending signals also send an excitatory current of 500 pA to the MNPs to allow the motor neurons to begin tonically spiking. The rhythm of the network's output is finally determined by phasic inhibition. From the RGPs, excitatory spikes are sent to the buffer NSIs. These NSIs act to inhibit the MNPs, mimicking the network of the stick insect (12, 21). The synaptic weight from the RGP determines the extent of membrane potential fluctuation of the NSI. We limit the fluctuation to ~15 mV to keep within a biologically-plausible range (29). The dependence of synaptic weight on frequency (V_th) is outlined in Equation (5).

$\begin{array}{ll}{w}_{nsi}=0.5125·V_{th}+31.4267& (5)\end{array}$

Where w_nsi is the synaptic weight from the RGPs to the NSIs and V_th is the voltage threshold potential of the AdEx neurons in the RGP population.

The inhibitory current from the NSIs to the MNPs is only dependent upon the momentary membrane potential of the NSI. This is described in Equation (6).

$\begin{array}{ll}\begin{array}{l}{I}_{mnp}=-\frac{1000}{15}·\Delta V_m\hfill \end{array}& (6)\end{array}$

Where I_mnp is the current from the NSI to the MNP and V_m is the membrane potential of the NSI. The selection of 15 as the denominator ensures that the maximum inhibition is −1,000 pA. The maximum is determined through trial and error.

2.2.2. Interjoint Coordination

Figures 2B,C show a high level overview of the neural network for coordinating firing between joints. Each of these sub-figures represents the same neuronal populations, it is only split so that the relevant synaptic connections can be highlighted based on the desired SRP output. Each joint is represented by a single half center oscillator to increase the clarity of the figure even though the controller per joint is as shown in Figure 2A. Network coordination is achieved through the use of cNSIs which receive excitatory spikes from the depressor or levator RGP depending on walking direction. The synaptic weight from the RGP to the cNSI is determined by Equation (7).

$\begin{array}{ll}{w}_{cnsi}=0.1462·V_{th}^2+16.5164·V_{th}+468.4534& (7)\end{array}$

Where w_cnsi is the synaptic weight from the pre-synaptic RGPs to the cNSI and V_th is the voltage threshold potential of the AdEx neurons in the RGPs. The weight is determined through trials and ensures the cNSI's membrane potential fluctuates by a maximum of 11 mV as is observed in biological NSIs found within networks generating SRPs (13). It should be noted that this change in membrane potential is smaller than the generic buffer NSIs used in the single joint architecture.

The cNSI connects to the post-synaptic RGPs with a weight determined by Equation (8). During coordinated firing, the cNSIs are selected as active by the descending signals, allowing current to pass through them and coupling the joints.

$\begin{array}{ll}{w}_{active}=33.3333·V_{th}+2091.6667& (8)\end{array}$

Where w_active is the synaptic weight from the cNSIs to the post-synaptic RGPs and V_th is the voltage threshold potential of the AdEx neurons in the RGPs. In order to send current from a cNSI to a post-synaptic RGP, the weight is multiplied by the change in membrane potential of the cNSI at each time step as shown in Equation (9).

$\begin{array}{ll}\begin{array}{l}{I}_{rgp\_act}=w_{active}·\Delta V_m\hfill \end{array}& (9)\end{array}$

Where I_{rgp_act} is the current from the cNSI to the post-synaptic RGP and V_m is the membrane potential of the cNSI.

During coordinated firing for the SRP1 swing to stance transition, the cNSIs are selected as active and Equation (9) determines the excitatory current sent to the Flx and Ret RGPs. During SRP2 coordination, the cNSIs send excitatory current to the Ext and Ret RGPs. The current injections from the cNSIs drive coordination between each of the joints but do not necessarily produce the same phase differences as recorded in deafferented samples. In order to tune the phase difference between joints, the synaptic delay from the CTr joint RGP to cNSI1 is modified (see synapse marked with “variable delay” in Figures 2B,C). The other synaptic delays are held constant and can be seen in Table 1.

TABLE 1

Table 1. Synaptic delays between spiking populations and NSIs.

Finally, Figures 2B,C show feedback current from cNSI2 to the driving CTr RGP. This feedback between the FTi and CTr joint has been observed in biological networks (30) and is added to our network to increase biological fidelity. The strength of the feedback current is set as 0.1% of the current from the cNSI (Equation 9). The strength of this current is determined through manual testing.

2.2.3. Protractor Burst Duration Modulation

Deafferented preparations of the stick insect indicate a decoupling of protraction burst duration from walking frequency during uncoordinated firing. (8)'s deafferented sample recordings show that retractor burst duration is correlated to cycle period but protractor burst duration is not. Therefore, this phenomenon must be accounted for when simulating uncoordinated firing to maintain biological plausiblity. We were able to reproduce the decoupling of protractor burst duration from cycle period using a burst duration NSI (bdNSI). The architecture is shown in Figures 2B,C. The presented solution is inspired by biological studies that find that blocking of calcium-dependent ion channels is able to lengthen burst duration of a neuron (15). We replicate this by increasing current to the protractor RGP and reducing inhibition from the retractor RGP. The prolonged excitation of the protractor RGP inhibits the protractor MNP, reducing the burst duration to the muscle. Implemented in practice, this means the retractor RGP sends excitatory spikes to a bdNSI. When the retractor RGP is spiking, the protractor MNP is also spiking so these excitatory spikes serve to inform the system that the protractor MNP has started spiking. Once the bdNSI's membrane potential is depolarized by 2 mV, the system waits for a period of time based on the frequency of stepping before sending an excitatory current of 3,000 pA to the protractor RGP and reducing inhibition from −5,000 to 0 nS from the retractor to protractor RGP. The waiting time is defined by Equation (10).

$\begin{array}{ll}{t}_{wait}=-0.4051·V_{th}^3-71.1881·V_{th}^2-4136.4234·V_{th}-79425.1892& (10)\end{array}$

Where t_wait is the amount of time in milliseconds that the protractor RGP is inhibited by the retractor RGP thus allowing the protractor MNP to spike. V_th is the voltage threshold potential of the AdEx neurons in the RGPs. As the cycle period increases due to the frequency slowing down, the amount of time that the protractor MNP is allowed to spike is reduced. The length of inhibition time to the protractor RGP is found through testing.

This mechanism is added to the network architecture and allowed to modulate the burst duration during uncoordinated firing. When the cNSIs are active to generate coordinated firing between joints, the bdNSI is inhibited so that burst duration is not modulated. The plot showing protractor burst duration remaining constant over increasing cycle period during uncoordinated firing is shown in Supplementary Figure 1.

Equations (4), (5), (7), (8), and (10) are found through manual testing, plotting the values, and using the polynomial fit function in Python.

2.2.4. Testing

The network architecture is confirmed by comparing simulation output to measurements from the deafferented stick insect. The phase data is extracted from each and equated by normalizing the timing into a 360° step cycle.

Figure 3 is a visual illustration of how the phase difference and burst duration are calculated from the simulation results. The calculation of biological phase difference is outlined in the Neurophysiology section (Figure 1).

FIGURE 3

Figure 3. Example output from the protractor and retractor MNPs over two seconds of simulation. A sliding time window of 5 ms is used to count the number of spikes per step. The phase difference is calculated by subtracting the time one MNP crosses a threshold from when the compared MNP crosses the same threshold. This time is divided by the cycle period and multiplied by 360 to produce the phase difference in degrees. The burst duration is calculated by subtracting when an MNP crosses the threshold with a positive slope from when the same MNP crosses the threshold with a negative slope.

The output from each MNP is rate-coded using a sliding time window of 5 ms. At each step, the number of spikes are counted within the time window and plotted as the y-value, the time window then moves by the time resolution of the simulation, 0.1 ms, and the spikes are counted again to plot the next point.

Each simulation is run for an equivalent of 8 s using the Neural Simulation Tool (NEST) (31) [source code available on GitLab; (32)]. The rate-coded spike data is saved and analyzed in MATLAB. Figure 3 shows how the calculations are made based on the simulation output. A threshold of 35 spikes is used because it removes noise from the data. When comparing phase difference of different MNPs, the time step for each MNP as it crosses the 35 spike threshold with a positive slope is used. The time step is converted to actual time ($\frac{time\text{\_}step\text{\_}value}{10000}=time\text{\_}in\text{\_}seconds$) and the values are subtracted and normalized to find phase in degrees, as shown in Equation (11).

$\begin{array}{ll}phase\text{\_}diff=\frac{t_{mnp1}-t_{mnp2}}{cycle\text{\_}period}·360& (11)\end{array}$

Where phase_diff is the phase difference in degrees between the MNPs bursting and t_mnp1, t_mnp1 are the times in seconds when each MNP crosses the threshold. The cycle period is determined by averaging the period over a single simulation during coordinated firing. This is used to calculate the individual phase differences between MNP bursts. The single phase difference value presented in the results is the average of all these calculations (Equation 11) during coordinated firing of a single simulation. There exists biological variability between animals as high as 32% in the most extreme case when comparing the phase difference of the Ret-Dep motor activity between observed Animals 1 and 2 (Supplementary Table 1). Therefore, we have allowed a buffer of 20° to either side of the circular mean ± angular deviation of a single animal (equivalent to 11%). This is increased to 40° (equivalent to 22%) if the angular deviation was <10°. This buffer allows a variability of 11–22%, keeping our simulation within the variation observed between animals.

The burst duration is calculated by looking at a single MNP's rate-coded output and noting when it crosses the threshold as the slope is increasing and decreasing. This calculation is not normalized, it subtracts the time steps and divides by 10, 000 to formulate the burst duration time in seconds. When testing burst duration, uncoordinated firing is allowed for the complete simulation time of 8 s and the presented burst duration time is an average of all phase difference calculations from a single simulation.

Noise is added to the network through current noise to the RGP and MNP neurons. The noise is centered around 0 mV with increasing amounts of standard deviation. All tests are run at 5 noise levels, starting at 100 pA of current noise and increasing to 500 pA by increments of 100 pA. The standard deviation of noise to all NSIs in the architecture is held constant at 25 pA.

3. Results

3.1. Neurophysiology

3.1.1. Intrasegmental CPGs Can Be Coordinated in the Deafferented Preparation

Extracellular recordings of the ipsilateral motor activity in the isolated and deafferented mesothoracic ganglion after pharmacological activation with pilocarpine reveal patterns of coordinated activity similar to the SRPs previously reported by (8) (Figure 4A). SRP1 and SRP2 are denoted in Figure 4A1,A2 with a dashed magenta or orange line, respectively. The lines demarcate the onset of the respective depressor burst. All highlighted patterns on Figure 4A1 belong to the SRP1 type. Each time there is combined activity of both the fast and the slow depressor units, there is a pause in the activity of the extensor motor neurons. SRP1 is the most frequently observed pattern of coordinated motor activity in this recording (Figure 4B). This substantiates previous results regarding the mesothoracic ganglion by (8). Note that not all depressor bursts are related to an SRP1, especially not those consisting of only the low-amplitude slow depressor units ((8)). Out of a total number of 1,093 depressor bursts recorded from six mesothoracic preparations only 37.9% are accompanied by an SRP1 or SRP2. The phase analysis only displays the three preparations (n = 621) in which the extensor motor neuron activity is also recorded. However, the phase histograms corresponding to all six preparations can be seen in Supplementary Figure 2. The mean distribution of protractor, retractor, and extensor spikes within the depressor cycle show distinct peaks and troughs (Figure 4C and Supplementary Figure 2). The mean angles of spiking activity relative to the depressor cycle with a 90% confidence interval are calculated for each animal preparation and the average values for all animals are given in Table 2. Some important numbers in the table are also highlighted in the text. Figure 4C (All bursts) illustrates the overall increase in retractor activity (circular mean [90% CI]: 113.2° [111.5, 114.9°]) and decrease in protractor activity (circular mean [90% CI]: 359.8° [358.4, 1.2°]) during depressor bursts. The figure also highlights the synchronous inactivity in the extensor motor neurons (circular mean [90% CI]: 231.6° [230.4, 232.8°]) during depressor activity (circular mean [90% CI]: 55° [54.4, 55.5°]). Considering only cycles where an SRP1 occurred (n = 130) the distributions appear to be similar, as exemplified by the similar circular means. However, the histograms display sharper peaks and the switch from protractor (circular mean [90% CI]: 4.6° [3.2, 5.9°]) to retractor activity (circular mean [90% CI]: 117.4° [115.6, 119.2°]) becomes more evident (Figure 4C, SRP1). In contrast, the last group of histograms show the opposite switch from retractor (circular mean [90% CI]: 19.4° [342.4, 56.5°]) to protractor activity (circular mean [90% CI]: 104.6° [94.3, 114.9°]) during ongoing depressor activity and extensor inactivity. This is indicative of an SRP2 (Figure 4C, SRP2).

FIGURE 4

Figure 4. (A1) Extracellular recording of the ipsilateral protractor (Pro), retractor (Ret), extensor (Ext), and depressor (Dep) motor neuron activity of an isolated (disconnected from other ganglia and deafferented) mesothoracic ganglion after bath application of 5 mM pilocarpine. Motor activity becomes transiently coordinated when a switch from protractor to retractor activity occurs (blue dashed line) during an ongoing depressor burst and a pause in extensor activity (SRP1). Magenta dashed lines mark the depressor burst onsets during SRP1s. (A2) Same preparation as in (A1). There is a switch from retractor to protractor activity (blue dashed line) during an ongoing depressor burst and a pause in extensor activity (SRP2). The orange dashed line marks the depressor burst onset during an SRP2. (B) SRP1s occur more often than SRP2s in the recording in (A). (C) Spike-phase histograms relative to the depressor cycle throughout the recording (all bursts) or during cycles where only an SRP1 or SRP2 occurs. In each histogram the mean (±STD) of each bin value among animal preparations is plotted. The y-axis represents average normalized number of spikes. “N” corresponds to the number of animal preparations and “n” to the number of depressor cycles.

TABLE 2

Table 2. Mean values of the circular mean and the upper/lower 90% confidence intervals (C.I.) of three isolated mesothoracic ganglia preparations.

Taken together, these results show that hemisegmental activity in the deafferented nerve cord can indeed be coordinated. Further revealing a higher occurrence of SRP1 in the isolated mesothoarcic ganglion. This represents a fictive transition from swing to stance phase during forward stepping.

3.1.2. SRP2 Occurs More Often in the Isolated Meta-Than in the Meso-Thoracic Ganglion

Extracellular recordings of the ipsilateral motor activity after pharmacological activation with pilocarpine also reveal the occurrence of SRPs in the isolated and deafferented metathoracic ganglion (Figure 5A1,A2). SRP2s are denoted in Figure 5A1 with a dashed orange line and SRP1s with a dashed magenta line (Figure 5A2), always crossing through the onset of the respective depressor burst. Similar to the mesothoracic ganglion, from a total number of 890 depressor bursts recorded from six metathoracic preparations only 35.7% are accompanied by an SRP1 or SRP2. Extensor activity pauses each time there is combined activity of both the fast and slow depressor units (data not shown in Figure 5A). SRP2 occurs more often than SRP1 in this recording (Figure 5B). The distribution of protractor and retractor spikes throughout 752 depressor cycles show distinct peaks and troughs (Figure 5C, All bursts). Again, the mean angles of spiking activity relative to the depressor cycle with a 90% confidence interval are calculated for each animal preparation and the average values for all animals are given in Table 3. Figure 5C (All bursts) illustrates the switch from retractor (circular mean [90% CI]: 239.3° [237.3, 241.4°]) to protractor activity (circular mean [90% CI]: 69.3° [68.8, 69.8°]) shortly after the onset of the depressor cycle (circular mean [90% CI]: 67.2° [66.5, 67.8°]), as denoted by the peak at ~0° in the retractor histogram. Considering only the cycles where an SRP1 occurs (n = 35), the initial switch from retractor to protractor activity is missing and only the opposite switch from protractor to retractor becomes evident at 40–45° (Figure 5C, SRP1). Here, the circular mean of the protractor (circular mean [90% CI]: 46° [45.3, 46.7°]) is smaller than the mean angle of the depressor (circular mean [90% CI]: 66.9° [65.8, 68°]). In the last group of histograms, based only on cycles where an SRP2 occurs, the initial switch from retractor to protractor activity during ongoing depressor activity reappears (Figure 5C, SRP2) and the circular mean of the protractor activity (circular mean [90% CI]: 80.8° [79.7, 81.9°]) is larger than the depressor mean (circular mean [90% CI]: 71.1° [69.1, 73.1°]). This is similar to the histograms corresponding to “All bursts.” The circular mean values point out the difference in the activity of the protractor and retractor between the deafferented preparations of the isolated meso- and meta-thoracic ganglia. Therefore, coordination of hemisegmental activity in the deafferented metathoracic ganglion is qualitatively different compared to the mesothoracic ganglion and there is a higher occurrence of SRP2s, representing a fictive transition from swing to stance phase during backward stepping.

FIGURE 5

Figure 5. (A1) Extracellular recording of the ipsilateral protractor (Pro), retractor (Ret), and depressor (Dep) motor neuron activity of an isolated metathoracic ganglion after bath application of 5 mM pilocarpine. Motor activity becomes transiently coordinated when a switch from retractor to protractor activity occurs (blue dashed line) during an ongoing depressor burst (SRP2). Orange dashed lines mark the depressor burst onsets during SRP2s. (A2) Same preparation as in (A1), magenta dashed lines mark the depressor burst onsets during SRP1s and the blue dashed line points out the switch from protractor to retractor activity. (B) SRP2s occur more often than SRP1s in the recording in (A). (C) Spike-phase histograms relative to the depressor cycle throughout the recording (All bursts) or during cycles where only an SRP1 or SRP2 occurs. In each histogram the mean (± STD) of each bin value among animal preparations is plotted. The y-axis represents average normalized number of spikes. “N” corresponds to the number of animal preparations and “n” to the number of depressor cycles. (D) Proportion of SRP1s calculated from the sum of SRP1 and SRP2 patterns observed in recordings of the isolated meso- and meta-thoracic, and interconnected metathoracic ganglia preparations. There is an overall higher occurrence of SRP1s in the isolated mesothoracic ganglion.

TABLE 3

Table 3. Mean values of the circular mean and the upper/lower 90% confidence intervals (C.I.) of five isolated metathoracic ganglia preparations.

Figure 5D is based on data from different animal preparations. This plot reinforces the qualitative difference in ipsilateral coordination between meso- and meta-thoracic ganglia. The relative proportions of SRP1s and SRP2s are measured in the isolated meso- and meta-thoracic ganglia, and in the interconnected metathoracic ganglion [(n = 6) in each case]. The proportion of SRP1s is significantly higher than 50% in the isolated mesothoracic ganglion (p = 0.016, one-sample Wilcoxon signed rank test), whereas five out of six isolated metathoracic ganglion preparations show SRP1 proportions below 50% (Figure 5D). SRP1 proportions in the isolated and interconnected metathoracic ganglion do not significantly differ from 50% (p = 0.219 and p = 0.437, respectively). However, four out of six interconnected metathoracic preparations show SRP1 proportions over 50%, implying a possible effect of intersegmental information on ipsilateral CPG coordination in the metathoracic ganglion (see Section 4). Taken together, our results indicate there is a higher frequency of SRP1s in the isolated meso- in contrast to a higher frequency of SRP2s in the isolated metathoracic ganglion.

3.2. Network Simulation

3.2.1. Network Topology Promotes Coordinated Motor Activity

The main result of the simulation is the network architecture itself. The configuration of a single joint pictured in Figure 2A is designed and confirmed to produce anti-phasic output based on the suggested architecture by the biological study from (21). This study informed our use of NSIs between the RGPs and MNPs, revealing a critical layer in the single joint architecture to allow precise control of output amplitude. The successful regulation of amplitude across frequencies is displayed in Figure 6.

FIGURE 6

Figure 6. Plot of rate-coded output from the MNPs and the RGPs. The amplitude (amount of spikes) changes according to frequency from the RGPs but remains constant for the MNPs. This is achieved by using a constant bias current to achieve tonic spiking by the MNPs and using rhythmic inhibition from NSIs driven by the RGPs to create oscillations.

Furthermore, we find NSIs to be an effective method to loosely couple joints because they separate network dynamics. The addition of NSIs for intra-leg communication is inspired by the observations of (13) indicating the role of NSIs during coordinated firing. In this way, our study provides corroborating evidence that NSIs may be important for coordination. Furthermore, our results indicate the insect may switch between different network architectures (communication pathways) based on walking direction. The architecture for the individual networks producing either forwards or backwards fictive step transitions is presented in Figures 2B,C.

3.2.2. Simulated Coordination Is Comparable to Biological Measurements

In order to confirm that the developed network architecture is biologically-plausible, its output is compared to biological measurements from a deafferented stick insect preparation. The comparison of network output to biological measurements starts with visual inspection to ensure the MNPs are bursting in the correct order to produce a swing to stance transition. Figure 7 shows the biological measurements on the left and the simulation results on the right for SRP1 and SRP2. As the measurements are at different timescales, the plots are normalized over the cycle period of each test, respectively.

FIGURE 7

Figure 7. Side by side comparison of biological measurements (left) and output from the network simulation (right) during SRP1 (top) and SRP2 (bottom). The x-axis is the phase in degrees for all plots. The y-axis is normalized smoothed spiking activity for the biological output or number of spikes for the simulation output. The biological output for SRP1 is from the same animal as pictured in Figure 4A, also referred to as “Animal 1.” The simulation output also pictures “Simulated Animal 1.” The biological output for SRP2 is from the same animal as pictured in Figure 5A, also referred to as “Animal 4.” The simulation output also pictures “Simulated Animal 4,” which is the same as “Simulated Animal 3” because the two animals' phase difference measurements for SRP2 are similar. During the depressor burst there is a switch from protraction to retraction in the SRP1 plots. Conversely, there is a switch from retraction to protraction during the depressor burst in the SRP2 plots. These transitions are highlighted within the blue boxes. Data was not available from the extensor during biological SRP2 recordings so this plot is empty but still pictured to keep consistency. The recordings from Animal 1 are from the isolated mesothoracic ganglion and Animal 4's recordings are from the isolated metathoracic ganglion.

The figure confirms that during SRP1 there is switch from protraction to retraction during the depressor burst. Similarly, during SRP2 retraction ends and protraction begins during the depressor burst. This pattern is seen in both the biological and simulated results.

All tests of the network architecture are completed by switching from uncoordinated to coordinated firing after 2.5 s, allowing the network to initialize before switching to coordination. Four different phase difference ranges are evaluated, mimicking measurements from four individual animals (Animals 1–4) producing either SRP1 or SRP2 fictive step transitions. Each phase range is achieved at four different frequencies from 1 to 4 Hz to ensure results are validated across observed walking speeds of the stick insect (28). The results confirm the developed architecture is able to produce a biological phase range when switching from uncoordinated to coordinated firing across the frequencies tested.

Figure 8 plots the calculated phase differences from the network simulation on top of the measured biological phase differences in motor activity from a single animal using the angular deviation to create the minimum and maximum values.

FIGURE 8

Figure 8. Phase difference between motor activity measured within deafferented samples (calculation method described by Figure 1) and in simulation (calculation method described by Figure 3 and Equation 11). The filled circles and lines indicate the measured biological phase differences using the angular deviation from the mean to create a minimum and maximum. The empty circles and dashed lines indicate the buffer as described in Section 2.2.4. The “X” marker indicates the phase difference recorded from simulation. There are four X's per motor activity comparison reflecting each of the four frequency values tested. Each plot represents a single animal for a total of four different animals represented. The plots are labeled accordingly from Animal 1 to 4. Animal 1 is the same animal as pictured in Figure 4A, Animal 4 is the same animal presented with results in Figure 5A. The recordings from Animals 1 to 2 are from the isolated mesothoracic ganglion and from the isolated metathoracic ganglion for Animals 3 and 4.

Testing reveals that the network architecture must be altered to produce SRP1 and SRP2. The depressor population drives coordination during SRP1 whereas the levator population is used during SRP2. This difference is a result of manually testing multiple architectures, finding that the topology in Figure 2C is the only one capable of meeting biological phase difference measurements. However, regardless of the driving population, tuning the single parameter of synaptic delay from the ThC joint RGP to cNSI1 (see Figures 2B,C in Section 2) is enough to match biologically measured phase differences in two different animals during both forward and backward walking. The exact synaptic delays from the RGP to cNSI1 used in this study can be found in Table 4.

TABLE 4

Table 4. Synaptic delay from the CTr RGPs to cNSI1 producing biological phase differences per animal—Animal 1 (A1), Animal 2 (A2), Animal 3 (A3), Animal 4 (A4).

The exact phase differences calculated from simulation output are also recorded in Supplementary Tables 2–4. As seen in Figure 8, the simulation is able to produce biologically-plausible phase ranges across the tested frequencies.

The system contains noise and investigation at different noise levels injected to the spiking populations indicates certain ideal levels depending on the individual animal's coordination and walking direction. A standard deviation of 400 pA is found to work for simulating Animal 1 to produce SRP1 while simulating Animal 2 needed 300 pA of current noise to produce the correct phase differences. Phase ranges exhibited by two different animals (Animals 3 and 4) producing SRP2s could be reached in simulation using 400 pA of current noise.

3.2.3. Applying the Simulated Network on a Robot Leg Shows Stepping Transitions

The output of the network is also tested on a physical robot simulator, CoppeliaSim (33), to visually confirm walking coordination on a robot leg. A link to the video can be found in the Supplementary Video 1. Figure 9 shows a sequence of stepping behavior obtained from the video during coordinated firing.

FIGURE 9

Figure 9. Stills from simulation of SRP1 using phase difference ranges measured from Animal 1 to create Simulated Animal 1. The stills reflect a full swing-stance step cycle.

The ability to control a single robot leg in simulation indicates the feasibility of eventually using this network with an additional interlimb coordination control mechanism (34) to control a physical legged robot.

4. Discussion

4.1. Neurophysiology

Our study finds that SRPs occur not only in the meso- but also in the metathoracic ganglion. We substantiate the higher occurrence of SRP1s in the mesothoracic ganglion, as previously described in (8), and we report a higher occurrence of SRP2s in the metathoracic ganglion. As SRP2 represents a fictive transition from swing to stance phase during backward stepping, it is plausible to assume there is an inherent tendency for forward walking in isolated meso- and for backward walking in isolated metathoracic ganglia. This result is in line with findings by (35), reporting that the hind legs of the stick insect perform backward stepping when they do not receive input from anterior legs, indicating an inherent backward direction of movement for the metathorax. Interestingly, this segment specificity in relation to stepping direction can also be demonstrated after sensory stimulation. Load increase at the level of the mesothorax after pilocarpine application results in retractor activation, whereas at the level of the metathorax the protractor is activated, pointing toward activation of forward and backward stepping, respectively (36). Based on the mentioned studies and our results, we conclude that the meso- and meta-thoracic networks for stepping may rely on different architectures.

We also show the proportion of SRP1s is above 50% in four out of six preparations of the interconnected metathoracic gangion. Future studies should increase the sample size and verify our observation, which implies that intersegmental information descending from the meso-thoracic ganglion may affect coordination of motor activity at the metathoracic ganglion, ultimately changing the fictive backward stepping phase transition (SRP2) to forward (SRP1). A similar change in the activity is also observed in vivo, as the inherent backward direction of the metathorax changes to forward when all other legs are stepping (35). Additionally, it has been found that intersegmental information from the middle legs is necessary to produce regular stepping in the hind legs in freely walking animals (37). Another notable phenomenon is that left-right coordination of motor activity is modified in the metathoracic ganglion when it is interconnected to the mesothoracic ganglion or the rest of the nerve cord in deafferented stick insect and locust preparations (10, 11). Considering all of the above, we present evidence for segment specificity in the walking system of the stick insect. The difference in the underlying neural mechanisms among segments and the functional significance of our results will be the focus of our future experiments.

4.2. Network Simulation

Our study reveals that the addition of 2 cNSIs is enough to coordinate joints and produce fictive step transitions. However, the synaptic connections to these cNSIs must be altered depending on the direction of walking. SRP1 requires the depressor, retractor, and flexor to fire together, therefore, coordination is created through coupling their respective RGPs directly (see Figure 2B in Section 2). SRP2 needs the levator, retractor, and extensor to fire at the same time so the architecture switches to driving coordination from the levator RGP and exciting the retractor and extensor RGPs (see Figure 2C in Section 2). This indicates that there is a distinct difference in synaptic connectivity during forwards and backwards step transitions and could explain why the metathoracic ganglion produces mostly SRP2s when disconnected from the other ganglia. This finding is similar to the biological outcome from (13) showing that exciting certain NSIs increased the likelihood of one type of SRP vs. another. The replication of this phenomenon in simulation by means of altering network connectivity, points to a partial role of network architecture in stepping patterns.

Figure 7 shows that the burst duration in simulation is much longer than in biological experiments. Therefore, there is more overlap when the MNPs are spiking. This does not affect the phase difference calculations since these only account for onset in spiking. Based on the visual observation of the robot leg in simulation, the transitions can be subjectively judged as acceptable even though there is significantly more spiking overlap. However, this must be further investigated when building upon the suggested network. It should also be noted that the phase difference calculation is from onset of spiking when evaluating the simulation data (Figure 3) as opposed to an average of all spikes in a burst when looking at biological data (Figure 1). We assert this is comparable because burst duration remains stable at any given frequency (8). Therefore, comparing the average vs. the onset will produce the same phase difference because the mean of the burst time remains equidistant to the onset of the burst time.

Our developed network architecture using both spiking and non-spiking neurons to regulate output from the MNPs should be considered for use in other research, expanding on the commonly used half-center oscillator (9, 38) architecture associated with CPGs. The addition of NSIs has several advantages including separating the dynamics of the system and increased control over input to spiking populations. As seen in (25), the output amplitude of the MNPs changes with frequency due to the use of a constant time window when rate-coding spikes. This variance is removed by adding buffer NSIs between the RGPs and MNPs (see Figure 6). The affect of the excitatory spikes from the RGPs is limited by the synaptic weight to the NSIs and neuronal dynamics of the NSIs themselves. The change in membrane potential of the NSI determines the inhibitory current sent to the MNPs acting as a way to smooth noisy spiking input. In addition to incoming spikes, current injection also affects an NSI's membrane potential. The amount of current and whether it is excitatory or inhibitory, regulates the membrane potential fluctuation of an NSI. The ability of NSI's to receive spikes or current and translate them into a predictable membrane potential fluctuation makes them useful for connecting neural network architectures or conveying sensory information to a neural network.

The rhythmicity of network output is driven by the RGP neurons. Coarse frequency modulation is achieved through changing the voltage threshold potential of the neurons in the RGPs. This relationship was found by (25) where a change of ~1 Hz is achieved by manipulating the voltage threshold potential by 1 mV. This relationship only holds until a voltage threshold minimum of V_th = −57 mV due to the other chosen neuronal characteristics. This V_th still produces an output of approximately 2Hz when using the original study's suggested excitatory current injection of 500 pA to the RGPs (25). Therefore, further tuning is required to reduce the frequency to the desired 1 Hz minimum. During testing, fine frequency modulation is achieved by adjusting the level of current injection to the RGP. Reducing the excitatory current to 20 pA is able to reduce the frequency to ~1.2 Hz which is accepted to be close enough to the desired minimum. The capacitance of the RGPs also needs to be increased to 400 pF to smooth the output at such low frequencies.

Phase difference is manipulated through synaptic delays from spiking populations to NSIs. The delay to cNSI1 (see Figures 2B,C) can be used to adjust the phase based on frequency. This tuning parameter is able to replicate the phase differences measured on animals displaying motor activity coordination as far as 114° apart (Supplementary Table 1), exhibiting the sensitivity of the system to this variable. Supplementary Video 1 of the simulation confirms that the motor activity coordination between biological animals is significantly varied. Visual observation shows the Simulated Animal 1 has a more intuitive phase difference during swing to stance transition which more closely resembles forward walking in a live animal. As noted in the results, the optimal amount of noise to produce the correct motor activity coordination in simulation for forward and backward swing to stance transitions per animal differed. This is considered acceptable since noise is by definition not constant.

The network simulation phase differences presented in Figure 8 are an average over all calculations occurring after coordinated firing is initiated. However, SRPs are spontaneous and not cyclic so further testing must be done to see if the network can produce singular instances of coordinated firing. Evaluating individual phase differences between motor activity shows that it can take between 0.14 and 2.85 s before the compared phase of the motor activity is within the accepted range as compared to biological measurements (see Supplementary Tables 2–4). This difference in “initialization time” could be caused by the significant amount of noise in the system and should also be further evaluated.

The main finding that two NSIs can be used to couple the leg joints to create coordinated firing allows for either descending signals or sensory feedback to be easily added to the system. NSIs can receive spikes or current injection, making them well-suited to connect with spiking sensory neurons or translate analog sensory information. The identification of analogous interneurons in the biological system and possible sensory organs affecting their activity should be the focus of future research.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://gitlab.sdu.dk/stroh/single-leg-coordination.

Author Contributions

BS and CM developed the main idea of the article. BS researched the biological neural network concepts, formulated a theory, and applied it in simulation. CM and DK performed the biological experiments. CM analyzed the associated data. The manuscript was written by BS and CM with support from LL, AB, and PM. AB supervised the project. LL and PM discussed the results. LL and AB provided the funding. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by the SDU Biorobotics Group at the University of Southern Denmark (BS, LL, and PM), the University of Cologne (CM, DK, and AB), and this project was supported by DFG grant Bu857/15-1 to AB.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to thank Mathias Thor for providing access to the MORF robot leg in the simulation environment.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/finsc.2022.818449/full#supplementary-material

References

1. Büschges A. Sensory control and organization of neural networks mediating coordination of multisegmental organs for locomotion. J Neurophysiol. (2005) 93:1127–35. doi: 10.1152/jn.00615.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

2. McCrea DA, Rybak IA. Organization of mammalian locomotor rhythm and pattern generation. Brain Res Rev. (2008) 57:134–46. doi: 10.1016/j.brainresrev.2007.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Markin SN, Lemay MA, Prilutsky BI, Rybak IA. Motoneuronal and muscle synergies involved in cat hindlimb control during fictive and real locomotion: a comparison study. J Neurophysiol. (2012) 107:2057–71. doi: 10.1152/jn.00865.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ekeberg O, Pearson K. Computer simulation of stepping in the hind legs of the cat: an examination of mechanisms regulating the stance-to-swing transition. J Neurophysiol. (2005) 94:4256–68. doi: 10.1152/jn.00065.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Schilling M, Hoinville T, Schmitz J, Cruse H. Walknet, a bio-inspired controller for hexapod walking. Biol Cybern. (2013) 107:397–419. doi: 10.1007/s00422-013-0563-5

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Tóth TI, Daun-Gruhn S. A three-leg model producing tetrapod and tripod coordination patterns of ipsilateral legs in the stick insect. J Neurophysiol. (2016) 115:887–906. doi: 10.1152/jn.00693.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Holmes P, Full RJ, Koditschek D, Guckenheimer J. The dynamics of legged locomotion: models, analyses, and challenges. SIAM Rev. (2006) 48:207–304. doi: 10.1137/S0036144504445133

CrossRef Full Text | Google Scholar

8. Büschges A, Schmitz J, Bässler U. Rhythmic patterns in the thoracic nerve cord of the stick insect induced by pilocarpine. J Exp Biol. (1995) 198:435–56. doi: 10.1242/jeb.198.2.435

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Bidaye SS, Bockemühl T, Büschges A. Six-legged walking in insects: how CPGs, peripheral feedback, and descending signals generate coordinated and adaptive motor rhythms. J Neurophysiol. (2018) 119:459–75. doi: 10.1152/jn.00658.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Mantziaris C, Bockemühl T, Holmes P, Borgmann A, Daun S, Büschges A. Intra-and intersegmental influences among central pattern generating networks in the walking system of the stick insect. J Neurophysiol. (2017) 118:2296–310. doi: 10.1152/jn.00321.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Knebel D, Ayali A, Pflüger HJ, Rillich J. Rigidity and flexibility: the central basis of inter-leg coordination in the locust. Front Neural Circuits. (2017) 10:112. doi: 10.3389/fncir.2016.00112

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Mantziaris C, Bockemühl T, Büschges A. Central pattern generating networks in insect locomotion. Dev Neurobiol. (2020) 80:16–30. doi: 10.1002/dneu.22738

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Büschges A. Role of local nonspiking interneurons in the generation of rhythmic motor activity in the stick insect. J Neurobiol. (1995) 27:488–512. doi: 10.1002/neu.480270405

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Chrachri A, Clarac F. Fictive locomotion in the fourth thoracic ganglion of the crayfish, Procambarus clarkii. J Neurosci. (1990) 10:707–19. doi: 10.1523/JNEUROSCI.10-03-00707.1990

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Grillner S. The motor infrastructure: from ion channels to neuronal networks. Nat Rev Neurosci. (2003) 4:573–86. doi: 10.1038/nrn1137

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Schilling M, Cruse H. Decentralized control of insect walking: a simple neural network explains a wide range of behavioral and neurophysiological results. PLoS Comput Biol. (2020) 16:e1007804. doi: 10.1371/journal.pcbi.1007804

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Graham D. A behavioural analysis of the temporal organisation of walking movements in the 1st instar and adult stick insect (Carausius morosus). J Compar Physiol. (1972) 81:23–52. doi: 10.1007/BF00693548

CrossRef Full Text | Google Scholar

18. Wosnitza A, Bockemühl T, Dübbert M, Scholz H, Büschges A. Inter-leg coordination in the control of walking speed in Drosophila. J Exp Biol. (2013) 216:480–91. doi: 10.1242/jeb.078139

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Ekeberg Ö, Blümel M, Büschges A. Dynamic simulation of insect walking. Arthrop Struct Dev. (2004) 33:287–300. doi: 10.1016/j.asd.2004.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

20. von Twickel A, Büschges A, Pasemann F. Deriving neural network controllers from neuro-biological data: implementation of a single-leg stick insect controller. Biol Cybern. (2011) 104:95–119. doi: 10.1007/s00422-011-0422-1

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Yeldesbay A, Daun S. Intra-and intersegmental neural network architectures determining rhythmic motor activity in insect locomotion. Commun Nonlin. Sci Numer Simul. (2020) 82:105078. doi: 10.1016/j.cnsns.2019.105078

CrossRef Full Text | Google Scholar

22. Goldammer BA J, Schmidt J. Motoneurons, DUM cells, and sensory neurons in an insect thoracic ganglion: a tracing study in the stick insect Carausius morosus. J Comp Neurol. (2012) 520:230–57. doi: 10.1002/cne.22676

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Berens P. CircStat: a MATLAB toolbox for circular statistics. J Stat Softw. (2009) 31:1–21. doi: 10.18637/jss.v031.i10

CrossRef Full Text | Google Scholar

24. Naud R, Marcille N, Clopath C, Gerstner W. Firing patterns in the adaptive exponential integrate-and-fire model. Biol Cybern. (2008) 99:335. doi: 10.1007/s00422-008-0264-7

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Strohmer B, Stagsted RK, Manoonpong P, Larsen LB. Integrating non-spiking interneurons in spiking neural networks. Front Neurosci. (2021) 15:184. doi: 10.3389/fnins.2021.633945

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Strohmer B, Manoonpong P, Larsen LB. Flexible spiking CPGs for online manipulation during hexapod walking. Front Neurorobot. (2020) 14:41. doi: 10.3389/fnbot.2020.00041

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Berendes V, Zill SN, Büschges A, Bockemühl T. Speed-dependent interplay between local pattern-generating activity and sensory signals during walking in Drosophila. J Exp Biol. (2016) 219:3781–93. doi: 10.1242/jeb.146720

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Graham D, Cruse H. Coordinated walking of stick insects on a mercury surface. J Exp Biol. (1981) 92:229–41. doi: 10.1242/jeb.92.1.229

CrossRef Full Text | Google Scholar

29. Burrows M, Siegler M. Graded synaptic transmission between local interneurones and motor neurones in the metathoracic ganglion of the locust. J Physiol. (1978) 285:231–55. doi: 10.1113/jphysiol.1978.sp012569

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Daun-Gruhn S, Büschges A. From neuron to behavior: dynamic equation-based prediction of biological processes in motor control. Biol Cybern. (2011) 105:71–88. doi: 10.1007/s00422-011-0446-6

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Jordan J, Mørk H, Vennemo SB, Terhorst D, Peyser A, Ippen T, et al. NEST 2.18.0. Zenodo. (2019). Available online at: https://zenodo.org/record/2605422#.YkgubBbhVPY (accessed March 11, 2022).

32. Strohmer B. Single Leg Coordination Source Code. (2021). Available online at: https://gitlab.sdu.dk/stroh/single-leg-coordination (accessed August 2, 2021).

33. Rohmer E, Singh SPN, Freese M. CoppeliaSim (formerly V-REP): a versatile and scalable robot simulation framework. In: Proceedings of the International Conference on Intelligent Robots and Systems (IROS). (2013). Available online at: www.coppeliarobotics.com

Google Scholar

34. Aoi S, Manoonpong P, Ambe Y, Matsuno F, Wörgötter F. Adaptive control strategies for interlimb coordination in legged robots: a review. Front Neurorobot. (2017) 11:39. doi: 10.3389/fnbot.2017.00039

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Bässler U, Foth E, Breutel G. The inherent walking direction differs for the prothoracic and metathoracic legs of stick insects. J Exp Biol. (1985) 116:301–11. doi: 10.1242/jeb.116.1.301

CrossRef Full Text | Google Scholar

36. Akay T, Ludwar BC, Göritz ML, Schmitz J, Büschges A. Segment specificity of load signal processing depends on walking direction in the stick insect leg muscle control system. J Neurosci. (2007) 27:3285–94. doi: 10.1523/JNEUROSCI.5202-06.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Grabowska M, Godlewska E, Schmidt J, Daun-Gruhn S. Quadrupedal gaits in hexapod animals-inter-leg coordination in free-walking adult stick insects. J Exp Biol. (2012) 215:4255–66. doi: 10.1242/jeb.073643

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Marder E, Bucher D. Central pattern generators and the control of rhythmic movements. Curr Biol. (2001) 11:R986–96. doi: 10.1016/S0960-9822(01)00581-4

PubMed Abstract | CrossRef Full Text | Google Scholar