100 [] [] []


# Rigid-Body Sound Synthesis with Differentiable Modal Resonators

###### Abstract

Physical models of rigid bodies are used for sound synthesis in applications
from virtual environments to music production. Traditional methods such as
modal synthesis often rely on computationally expensive numerical solvers,
while recent deep learning approaches are limited by post-processing of their
results. In this work we present a novel end-to-end framework for training a
deep neural network to generate modal resonators for a given 2D shape and
material, using a bank of differentiable IIR filters. We demonstrate our
method on a dataset of synthetic objects, but train our model using an audio-
domain objective, paving the way for physically-informed synthesisers to be
learned directly from recordings of real-world objects.

Index Terms—  differentiable signal processing, machine learning, sound
synthesis, physical modelling

##  1 Introduction

Fig. 1: A schematic diagram of our model and its training pipeline. A CNN
generates shape embeddings from a synthetic dataset of 2D convex shapes, which
are concatenated with material parameters and excitation co-ordinates and
passed to an MLP to generate parameters for a differentiable IIR filterbank.

The synthesis of contact sounds from rigid bodies and materials has been of
continual interest in applications from music production, through sound
design, to the rendering of object sounds in virtual environments [1, 2]. For
this reason, the problem has been extensively studied, and physically-based
numerical methods are often selected for its solution. However, such
approaches typically incur significant computational cost, as well as storage
cost for cached solutions if real-time interaction is required. These
limitations inhibit flexibility in such applications, where adapting to new
object shapes and materials requires prohibitively time-consuming computation
of the new solution.

Of these numerical approaches, which ultimately rely on solving the wave
equation, the finite difference time domain and the finite element method
(FEM) are most commonly used due to their adaptability. FEM solvers are used
to precompute the vibrational modes of arbitrarily shaped rigid objects,
allowing contact sounds for these objects to be synthesized using modal
synthesis.

The computation of the modes is usually posed as a generalized eigenvalue
problem using the mass and stiffness matrices of an object. Using the
solution, sound can be rendered with an oscillator bank by projecting an
impulse onto the modes at each discrete location within the object. This
approach to computing an object’s modes can, however, be computationally
expensive. Moreover, every time the shape or material characteristics of the
object change, the system must be solved again. Thus, numerous approaches have
been proposed to accelerate this process [3, 4, 5].

###  1.1 Data-driven methods

A recent generative method, proposed by Traer et al. [6], used a perceptually
derived statistical model to approximate object impulse responses by
approximating the modes within some degree of perceptual tolerance. While this
method circumvents the need for a numerical solver and uses filters to account
for different object interactions, it relies on abstract intermediate
representations of an object’s physical characteristics, and does not take
into account object geometry.

Jin et al. [7, 8, 9] trained a deep neural network to produce an object’s
modes, circumventing the need for a FEM solver. Their model predicted
eigenvalues and eigenvectors from sparsely voxelized objects, using as
supervision the corresponding modes obtained from the solver. Additionally, in
order to generalize to different sizes and materials, the authors suggest the
adaptive scaling of the modes based on material parameters (except for
Poisson’s ratio) and the size of the object. Similarly, our method uses a
discretized shape representation as input to a neural network. However, unlike
Jin et al.’s model, no post-processing is required to account for different
material parameters. Further, our model is not limited to querying modes at
discrete positions in the shape, allowing for arbitrary co-ordinate input.

###  1.2 Differentiable resonators

While it is possible to synthesize modes using an oscillator bank where each
damped sinusoid corresponds to an impulsed mode, resonant filter banks are
commonly used for synthesis due to their flexibility, as they can be excited
with a variety of signals [10]. Realising such a filter bank, however,
requires the corresponding eigenvalue and eigenvector problems to be solved in
order to derive the filter coefficients. In this work, we propose an end-to-
end learning method for generating resonant filter banks without explicitly
solving the system.

Our method relies on the use of differentiable infinite impulse response (IIR)
filters, which allow a neural network to directly control filter coefficients
by propagating loss gradients through them. Whilst it is possible to directly
specify such a filter in its recurrent form [11, 12], this requires
backpropagation of loss gradients through time which can introduce instability
in the form of exploding or vanishing gradients, and is poorly suited to
parallel computation on a GPU. For these reasons, other recent approaches [13,
14] have used the frequency sampling method during optimisation. This approach
allows a finite impulse response (FIR) approximation of an IIR filter’s
frequency response to be computed at arbitrary resolution, whilst retaining
full differentiability.

###  1.3 Our contribution

In this work, we set out to address the shortcomings with existing neural
network approaches using an end-to-end learning paradigm with differentiable
signal processors. Our proposed method can synthesise contact sounds from
arbitrary 2D shapes and material parameters, without assuming any particular
damping model, in a fraction of the time required for a numerical method. Our
network is coupled with a differentiable IIR filter bank which can be
configured with varying arrangements of cascaded and parallel biquads.

##  2 Method

Fig. 1 provides an overview of our method, illustrating how we sample
synthetic data, prepare inputs, and train our model.

###  2.1 Learnable modules

To allow efficient synthesis of contact sounds at multiple object positions
for a given shape, our system contains two separate neural networks which are
jointly optimised. The first is a convolutional neural network (CNN), using
the EfficientNet-B0 architecture [15], which generates an embedding from two-
dimensional shapes. The network takes as input a 2D occupancy grid describing
the shape of an object in discretized space and outputs an embedding of 1000
dimensions. For each distinct shape, only a single forward pass through this
module is required.

The second learnable module is a multi-layer perceptron (MLP) which takes as
input the shape embedding concatenated with material parameters and the co-
ordinates of the vertex at which the shape will be excited. The cached output
of the CNN can be reused to evaluate different excitation positions and
materials for the same shape, eschewing the need to recompute a shape
embedding. The material parameters
$\boldsymbol{\phi}=(\rho,E,\nu,\alpha,\beta)$ consist of mass density
$(\rho)$, Young’s modulus ($E$), Poisson’s ratio ($\nu$), and Rayleigh damping
coefficients ($\alpha,\beta$). Both excitation co-ordinates $\mathbf{x}$ and
material parameters $\boldsymbol{\phi}$ are normalized in the range [0, 1].
The MLP’s output is used to parameterise a differentiable IIR filter bank.

###  2.2 Differentiable filter bank

Each biquad in our differentiable filter bank is parameterized by a single
pole $p_{l,m}$ and zero $q_{l,m}$, where $l$ and $m$ correspond to the
filter’s index in the parallel and cascade dimensions, respectively. We
constrain the remaining pole and zero of each filter to the complex conjugates
of these parameters. Each biquad also has a gain parameter $k_{l,m}$ applied
to its numerator. In order to encourage filters with strong resonant peaks, we
specify a fixed bias for pole and zero parameters, initialising them near the
edge and centre, respectively, of the unit circle.

To ensure filter stability, we constrain the complex poles to be within the
unit circle following the method described in [16]. Specifically, we apply an
activation function $h:\mathbb{C}\to\mathbb{C}$ defined as:

$h(p)=\frac{\tanh|p|}{|p|}\cdot p$ (1)

We do not apply such a constraint to the zeros as we are unconcerned with
preserving minimum phase. The transfer function of each filter is thus given
by:

$H_{l,m}(z)=\frac{1-2\mathfrak{Re}(q_{l,m})z^{-1}+\left|q_{l,m}\right|^{2}z^{-2}}{1-2\mathfrak{Re}(h\left(p_{l,m}\right))z^{-1}+\left|h\left(p_{l,m}\right)\right|^{2}z^{-2}}$
(2)

While training the model, we do not evaluate the filter’s exact recursive
implementation, and instead opt, as in recent work on differentiable IIR
modelling [13, 14], to approximate the filter response by sampling it in the
frequency domain. That is, we evaluate the filter’s transfer function at
discrete points $e^{j2\pi\frac{k}{N}}$ on the unit circle, where $N$ is the
signal length and $k\in\\{0,\dots,N/2+1\\}$. Whilst this method incurs a loss
in precision, it allows the filter computation to be performed in parallel on
a GPU, whilst also sidestepping the aforementioned issues with training
stability and exploding gradients.

Our filter bank consists of a flexible number of parallel filters, each
consisting of a cascade of some number of second order (biquadratic) sections.
The cascade depth $M$ and number of parallel filters $L$ are thus tunable
hyperparameters of our method, with the overall transfer function of the
filter bank given by:

$H(z)=\sum_{l=1}^{L}\prod_{m=1}^{M}k_{l,m}H_{l,m}(z)$ (3)

###  2.3 Loss

We train our model to minimize the following loss function:

$\mathcal{L}=\lambda\left\|X_{\text{mel}}-H_{\text{mel}}\right\|_{2}^{2}+\gamma\left\|\log
X_{\text{mel}}-\log H_{\text{mel}}\right\|_{2}^{2},$ (4)

where $\lambda$ and $\gamma$ are loss scaling hyperparameters, and
$X_{\text{mel}}$ and $H_{\text{mel}}$ are the mel-scaled magnitude spectra of
the discrete Fourier transform (DFT) and filter frequency response,
respectively. In our experiments we used $\lambda=1.0$, and $\gamma=0.1$.

###  2.4 Dataset

To train our model, we generated a synthetic dataset consisting of 500 convex
shapes [17] along with their material parameters and audio synthesised using
modal synthesis.

In particular, each convex shape was represented by a point set $P$ which was
triangulated using Delaunay triangulation. Subsequently, each triangular
surface mesh $\mathcal{M}=(V,E,F)$ with a set of $V$ vertex positions, $E$
edges and $F$ faces was rasterized into a $64\times 64$ occupancy grid.
Materials were homogeneous and isotropic, and their parameters were sampled
uniformly in the intervals $\rho\in[500,15000]$, $E\in[8\times 10^{9},5\times
10^{10}]$, $\nu\in[0.1,0.4]$, $\alpha\in[1,10]$ and $\beta\in[3\times
10^{-7},2\times 10^{-6}]$. For each shape, 500 random materials were sampled.

To assemble the stiffness and density matrices and solve the generalized
eigenvalue problem corresponding to each shape-material pair, we used the
scikit-fem [18] package. We computed the first 32 modes for each object and
projected a unit impulse for each node in the mesh. Since each mesh has a
distinct tessellation, the number of triangles per mesh varied between
approximately 300 and 1000. Each vertex thus represents a unique combination
yielding a total of approximately 100 million samples.

###  2.5 Training

We trained our models using the Adam optimizer with an initial learning rate
of $3\times 10^{-5}$ and an exponential decay of $0.9$ every 300 steps. Models
were trained until validation loss did not improve for 20k steps. All models
were trained with a batch size of 64 on a single NVIDIA RTX A5000 GPU.

##  3 Results

Fig. 2: Inference results for test shapes, materials and positions. Top:
Occupancy maps and indicated hit position in red. Middle: FEM results’
spectrogram using traditional modal synthesis. Bottom: Results using our
proposed method. Here the log-frequency log-magnitude spectrograms show that
the generated results closely resemble the original approach.

In Fig. 2 we illustrate our network’s predictions for samples from a test
dataset of 100 materials and shapes. We observe that our method closely
matches the frequency distribution of the target signals, while generalising
well across shapes. We provide comprehensive examples of audio rendered with
our system in the online supplement at http://rodrigodzf.github.io/iir-
modal/results.html, including audio from the objects in Fig. 2, and
interpolation of shape and material parameters.

Our model appears to produce shorter decay times than those in the target
signal for modes of higher frequency. Further, we note that the complex
activation in Eqn 1 causes the gradient magnitude to decrease as a pole
approaches the unit circle. It is thus possible that for higher frequency
modes with lower amplitudes, this diminishing gradient slowed learning of
these resonances resulting in faster decays.

###  3.1 Computation time

To ascertain our method’s computational performance, we compared its
processing time to that of a conventional FEM solver (as implemented in
scikit-fem). Table 1 lists the processing times of both methods for 96, 426,
and 1792 vertices averaged across 100 meshes of different resolutions.

Our method produces filter bank parameters for a set of 96 arbitrary co-
ordinates in a mesh in 16.75 milliseconds using a GPU. Subsequently rendering
the audio is possible at negligible computational cost using a recursive
implementation of the filter bank. In comparison, modal synthesis using the
FEM solver takes approximately 83.66 milliseconds to calculate the mode gains
and frequencies for a mesh with the same number of vertex positions. The
solutions are only available at the discrete coordinate positions of the
vertices. Increasing the number of mesh vertices, and thus obtaining solutions
at finer steps, increases the duration of processing.

To produce a mesh with 1792 vertices, the FEM solver takes slightly over 1.7
seconds, while our method requires only 0.08 seconds. Note that in these
experiments we generate a shape embedding even though the occupancy grid does
not necessarily change in the model input, and thus this can be considered a
worst case scenario bound on the actual runtime of our model. If we cache the
shape embedding, the total time is reduced to 0.01 seconds for this number of
vertices. Moreover, as our method does not need a re-discretization of the
shape, we can query filter parameters at any continuous valued co-ordinate in
the shape with the same constant processing time. Further, our method is not
restricted to simultaneously querying multiple points in the mesh to render
the sound for a single position. Our approach is therefore preferable in real-
time and interactive synthesis for arbitrary shapes.

| 96 vertices | 426 vertices | 1792 vertices
---|---|---|---
FEM | 83.66ms | 387.12ms | 1728.25ms
Ours | 16.75ms | 22.03ms | 85.86ms
Table 1: Processing time in milliseconds for the FEM solver and our method. We
average the processing time for 100 meshes of different resolutions. For the
solver, we measure the time required to solve the system for each mesh. In our
case, we add the inference time for a shape and material pair and the
inference time for a fixed number of coordinate positions.

###  3.2 Filter-bank configuration

To ascertain the optimal trade off between the cascade depth $M$ and number of
parallel filters $L$, we compared three filter bank configurations such that
the total number of poles was constant, fixed at $M\cdot L=64$. We evaluated
models trained using these configurations on the test dataset using both the
frequency domain mean-absolute and mean-squared error. In general, we observed
that filter banks with a lower cascade depth and thus a lower filter order
converged faster.

Table 2 lists the results for these experiments. Convergence was very poor for
the highest order filter bank, with dramatically higher errors than either of
the other configurations. The 4 order filter bank achieved the best
performance across both metrics, suggesting a hybrid cascade/parallel topology
is optimal.

Filters and order | MAE | MSE
---|---|---
16 filters 8 order | 4.324 | 20.854
32 filters 4 order | 0.504 | 0.574
64 filters 2 order | 0.568 | 0.697
Table 2: Average log-spectrogram MAE and MSE for different filter
configurations. We test our model using a different number of filters in
parallel of different order.

##  4 Conclusion

This paper presented a novel approach for synthesising contact sounds for
rigid surfaces of arbitrary shape and material composition. Our method
generates parameters for a resonant IIR filter bank from shape, material, and
excitation position inputs. The generated filter bank is suitable for real-
time implementation, allowing realistic contact sounds to be produced
interactively using a variety of different excitation signals.

In this work, we limited shapes to two dimensions, but we believe that
extending the proposed method to 3D objects simply requires adapting the shape
encoder to use 3D convolutional layers. Additionally, a future adaption of our
model could apply a discretization agnostic network [19] to learn from
different surface representations. The rest of the pipeline would not require
any alteration.

We further consider that, adapting a method previously applied to generative
models [20, 21], our trained model could also be used to optimise a second
network for the task of predicting materials and shapes from sound.

Finally, since the predicted filter bank can be excited with any signal, we
intend to explore the possibilities of complex interaction on the object’s
surface, along with a qualitative evaluation.

## References

  * [1] James F O’Brien, Chen Shen, and Christine M Gatchalian,  “Synthesizing sounds from rigid-body simulations,”  in Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation, 2002, pp. 175–181.
  * [2] Kees Van Den Doel, Paul G Kry, and Dinesh K Pai,  “Foleyautomatic: physically-based sound effects for interactive simulation and animation,”  in Proceedings of the 28th annual conference on Computer graphics and interactive techniques, 2001, pp. 537–544.
  * [3] Dominik Rausch,  Modal sound synthesis for interactive virtual environments,  Ph.D. thesis, Universitätsbibliothek der RWTH Aachen, 2017.
  * [4] Nicolas Bonneel, George Drettakis, Nicolas Tsingos, Isabelle Viaud-Delmon, and Doug James,  “Fast modal sounds with scalable frequency-domain synthesis,”  in ACM SIGGRAPH 2008 papers, pp. 1–9. 2008.
  * [5] Jeffrey N Chadwick, Changxi Zheng, and Doug L James,  “Faster acceleration noise for multibody animations using precomputed soundbanks,”  in Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2012, pp. 265–273.
  * [6] James Traer, Maddie Cusimano, and Josh H McDermott,  “A perceptually inspired generative model of rigid-body contact sounds,”  in The 22nd International Conference on Digital Audio Effects (DAFx-19), 2019.
  * [7] Xutong Jin, Sheng Li, Tianshu Qu, Dinesh Manocha, and Guoping Wang,  “Deep-modal: real-time impact sound synthesis for arbitrary shapes,”  in Proceedings of the 28th ACM International Conference on Multimedia, 2020, pp. 1171–1179.
  * [8] Xutong Jin, Sheng Li, Dinesh Manocha, and Guoping Wang,  “Deepeigen: Learning-based modal sound synthesis with acoustic transfer maps,”  arXiv preprint arXiv:2108.07425, 2021.
  * [9] Xutong Jin, Sheng Li, Guoping Wang, and Dinesh Manocha,  “Neuralsound: learning-based modal sound synthesis with acoustic transfer,”  ACM Transactions on Graphics (TOG), vol. 41, no. 4, pp. 1–15, 2022\.
  * [10] Simon Conan, Olivier Derrien, Mitsuko Aramaki, Sølvi Ystad, and Richard Kronland-Martinet,  “A synthesis model with intuitive control capabilities for rolling sounds,”  IEEE/ACM transactions on audio, speech, and language processing, vol. 22, no. 8, pp. 1260–1273, 2014.
  * [11] Boris Kuznetsov, Julian D Parker, and Fabián Esqueda,  “Differentiable iir filters for machine learning applications,”  in Proc. Int. Conf. Digital Audio Effects (eDAFx-20), 2020, pp. 297–303.
  * [12] Giovanni Pepe, Leonardo Gabrielli, Stefano Squartini, and Luca Cattani,  “Designing audio equalization filters by deep neural networks,”  Applied Sciences, vol. 10, no. 7, pp. 2483, 2020.
  * [13] Shahan Nercessian, Andy Sarroff, and Kurt James Werner,  “Lightweight and interpretable neural modeling of an audio distortion effect using hyperconditioned differentiable biquads,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 890–894.
  * [14] Joseph T Colonel, Christian J Steinmetz, Marcus Michelen, and Joshua D Reiss,  “Direct design of biquad filter cascades with deep learning by sampling random polynomials,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 3104–3108.
  * [15] Mingxing Tan and Quoc Le,  “Efficientnet: Rethinking model scaling for convolutional neural networks,”  in International conference on machine learning. PMLR, 2019, pp. 6105–6114.
  * [16] Shahan Nercessian,  “Neural parametric equalizer matching using differentiable biquads,”  in Proc. Int. Conf. Digital Audio Effects (eDAFx-20), 2020, pp. 265–272.
  * [17] P Valtr,  “Probability that n random points are in convex position,”  Discrete & Computational Geometry, vol. 13, no. 3-4, pp. 637–643, 1995.
  * [18] Tom Gustafsson and G. D. McBain,  “scikit-fem: A Python package for finite element assembly,”  Journal of Open Source Software, vol. 5, no. 52, pp. 2369, 2020\.
  * [19] Nicholas Sharp, Souhaib Attaiki, Keenan Crane, and Maks Ovsjanikov,  “Diffusionnet: Discretization agnostic learning on surfaces,”  ACM Transactions on Graphics (TOG), vol. 41, no. 3, pp. 1–16, 2022\.
  * [20] Han Han and Vincent Lostanlen,  “wav2shape: Hearing the shape of a drum machine,”  arXiv preprint arXiv:2007.10299, 2020.
  * [21] Zhoutong Zhang, Jiajun Wu, Qiujia Li, Zhengjia Huang, James Traer, Josh H McDermott, Joshua B Tenenbaum, and William T Freeman,  “Generative modeling of audible shapes for object perception,”  in Proceedings of the IEEE International Conference on Computer Vision, 2017, pp. 1251–1260.


date02042015

# Distributed Recurrent Neural Forward Models with Synaptic Adaptation for
Complex Behaviors of Walking Robots

Sakyasingha Dasgupta  Correspondence<EMAIL_ADDRESS>Current
address: Riken Brain Science Institute, 2-1 Hirosawa, Wako, Saitama, Japan
Institute for Physics - Biophysics, George-Agust-University, Göttingen,
Germany Bernstein Center for Computational Neuroscience, George-Agust-
University, Göttingen, Germany  Riken Brain Science Institute, 2-1 Hirosawa,
Wako, Saitama, Japan  Dennis Goldschmidt Bernstein Center for Computational
Neuroscience, George-Agust-University, Göttingen, Germany  Florentin Wörgötter
Institute for Physics - Biophysics, George-Agust-University, Göttingen,
Germany Bernstein Center for Computational Neuroscience, George-Agust-
University, Göttingen, Germany  Poramate Manoonpong Bernstein Center for
Computational Neuroscience, George-Agust-University, Göttingen, Germany
Maersk Mc-Kinney Moller Institute, University of Southern Denmark, Odense,
Denmark

(date)

###### Abstract

Walking animals, like stick insects, cockroaches or ants, demonstrate a
fascinating range of locomotive abilities and complex behaviors. The
locomotive behaviors can consist of a variety of walking patterns along with
adaptation that allow the animals to deal with changes in environmental
conditions, like uneven terrains, gaps, obstacles etc. Biological study has
revealed that such complex behaviors are a result of a combination of
biomechanics and neural mechanism thus representing the true nature of
embodied interactions. While the biomechanics helps maintain flexibility and
sustain a variety of movements, the neural mechanisms generate movements while
making appropriate predictions crucial for achieving adaptation. Such
predictions or planning ahead can be achieved by way of internal models that
are grounded in the overall behavior of the animal. Inspired by these
findings, we present here, an artificial bio-inspired walking system which
effectively combines biomechanics (in terms of the body and leg structures)
with the underlying neural mechanisms. The neural mechanisms consist of 1)
central pattern generator based control for generating basic rhythmic patterns
and coordinated movements, 2) distributed (at each leg) recurrent neural
network based adaptive forward models with efference copies as internal models
for sensory predictions and instantaneous state estimations, and 3) searching
and elevation control for adapting the movement of an individual leg to deal
with different environmental conditions. Using simulations we show that this
bio-inspired approach with adaptive internal models allows the walking robot
to perform complex locomotive behaviors as observed in insects, including
walking on undulated terrains, crossing large gaps as well as climbing over
high obstacles. Furthermore we demonstrate that the newly developed recurrent
network based approach to online forward models outperforms the adaptive
neuron forward models, which have hitherto been the state of the art, to model
a subset of similar walking behaviors in walking robots.

##  1 Introduction

Walking animals show diverse locomotor skills to deal with a wide range of
terrains and environments. These involve intricate motor control mechanisms
with internal prediction systems and learning (Huston and Jayaraman 2011),
allowing them to effectively cross gaps (Blaesing and Cruse 2004), climb over
obstacles (Watson et al. 2002), and even walk on uneven terrain (Pearson and
Franklin 1984), (Cruse 1976). These capabilities are realized by a combination
of biomechanics of their body and neural mechanisms. The main components of
these neural mechanisms include central pattern generators (CPGs), internal
forward models, and limb-reflex control systems. The CPGs generate basic
rhythmic motor patterns for locomotion, while the reflex control employs
direct sensory feedback (Pearson and Franklin 1984). However, it is argued
that biological systems need to be able to predict the sensory consequences of
their actions in order to be capable of rapid, robust, and adaptive behavior.
As a result, similar to the observations in vertebrate brains (Kawato 1999),
insects can also employ internal forward models as a mechanism to predict
their future state (predictive feedbacks) given the current state or sensory
context (sensory feedback) and the control signals (efference copies), in
order to shape the motor patterns for adaptation (Webb 2004),(Mischiati et al.
2015). Essentially, such a forward model acts as an internal feedback loop,
that uses a copy of the motor command, in order to predict the expected
sensory input. Comparing this to the actual input, appropriate modulations of
this signal or adaptive behaviors can be carried out.

In order to make such accurate predictions of future actions to satisfy
changing environmental demands, the internal forward models require some
degree of memory of the previous sensory-motor information. However, given
that, such motor control happens on a very fast timescale, keeping track of
temporal information is integral to such very short-term memory processes.
Reservoir-based recurrent neural networks (RNNs) (Maass et al. 2002), (Jaeger
and Haas 2004), (Sussillo and Abbott 2009), with their inherent ability to
deal with temporal information and fading memory of sensory stimuli, thus
provide a suitable platform to model such internal predictive mechanisms.
Taking this perspective, here, we utilize a newly developed model of self-
adaptive reservoir networks (SARN) (Dasgupta et al. 2013), (Dasgupta 2015), to
act as the forward models for sensorimotor prediction. This works in
conjunction with other neural mechanisms for motor control and generates
complex adaptive locomotion in an artificial walking robotic system.
Specifically, by exploiting the adaptive recurrent layer of our model it is
possible to achieve complex motor transformations at different walking gaits,
which is significantly difficult to achieve by currently existing adaptive
forward models employed with walking robots (Manoonpong et al. 2013), (Dearden
and Demiris 2005), (Schröder-Schetelig et al. 2010).

We present for the first time a distributed forward model architecture using
six SARN-based forward models on a hexapod robot, each of which is for sensory
prediction and state estimation of an individual robot leg. The outputs of the
models are compared with foot contact sensory signals (actual sensory
feedback) and the differences between them are used for motor adaptation, in
an online manner. This is integrated as part of the neural mechanism framework
consisting of 1) single central pattern generator-based control for generating
basic rhythmic patterns and coordinated movements, 2) distributed reservoir
forward models and 3) searching and elevation action control for adapting the
movement of an individual leg based on the forward model predictions, in order
to deal with changing environmental conditions.

In the following section we describe the architectural setup of the neural
mechanisms used for the design of adaptive locomotion control in a walking
robot, along with a description of the simulated hexapod robot AMOS II and the
modular robot control environment used as the development platform for our
proposed control system. In section 3, we present the materials and methods
used in this study. Specifically, we introduce the setup and implementation of
the distributed reservoir-based adaptive forward model, with details of the
learning procedure. Section 4 presents experimental results of the learning
mechanism and the resulting behaviors of the simulated hexapod AMOS II on
different complex locomotion scenarios likes crossing a large gap, walking on
uneven (rough) terrains, and overcoming obstacles. The results obtained from
the reservoir based forward models are juxtaposed with the previous state of
the art adaptive neuron forward models setup. Finally, in section 5, we
discuss our results and provide an outlook of further future directions.

##  2 Neural Mechanisms for Complex Locomotion

The neural mechanisms (Figure 1 a) for locomotion control, are designed based
on a modular architecture, such that, they comprise of, i) central pattern
generator (CPG)-based control, ii) reservoir-based adaptive forward models,
and iii) searching and elevation action control. The CPG-based control and the
searching and elevation control have been previously discussed in detail in
(Manoonpong et al. 2013), thus here we will only provide a brief overview of
these mechanisms, while the reservoir-based adaptive forward models, which
forms the main topic of this work, will be presented in detail in the
following section.

The CPG-based control primarily generates a variety of rhythmic patterns and
coordinates all leg joints of a simulated hexapod robot AMOSII (Figure 1 (b)),
thereby, leading to a multitude of different behavioral patterns and insect-
like leg movements. The patterns include omnidirectional walking and insect-
like gaits (Manoonpong et al. 2013). All these patterns can be set manually,
or autonomously driven by exteroceptive sensors, like a camera (Zenker et al.
2013), a laser scanner (Kesper et al. 2013), or range sensors. While the CPG-
based control provides versatile autonomous behaviors, the searching and
elevation control at each leg uses the accumulated error signals provided by
the reservoir-based adaptive forward models in order to adapt the movement of
an individual leg of the robot and deal with changes in environmental
conditions.

Figure 1: (a) The closed-loop architectural diagram of an artificial bio-
inspired walking system consisting of the sensors (i.e., proprioceptive and
exteroceptive sensors) that receive environmental inputs and feedback, the
neural mechanisms (i,ii,iii) for adaptive locmotion control, and the
biomechanical setup of the hexapod robot AMOSII (i.e., six 3-jointed legs, a
segmented body structure with one active backbone joint (BJ), actuators, and
passive compliant components (Manoonpong et al. 2013)). (b) Modular Robot
Control Environment embedded in the LPZRobots simulation toolkit (Der and
Martius 2012), (Hesse et al. 2012). The simulation environment provides the
main testbed for developing the controller, testing it on the simulated
hexapod robot, and finally transferring it to the physical agent. Here we
evaluate our model and results primarily on the simulated robot (bottom left),
which accurately embodies the characteristics of its physical equivalent, AMOS
II robot (bottom right). Here, $FC_{1}$, $FC_{2}$, $FC_{3}$, $FC_{4}$,
$FC_{5}$, and $FC_{6}$ are foot contact sensors installed in the robot legs,
which are used as the main sensory stimuli compared against the predicted
signal from the RNN-based (reservoir) forward models. Each leg (right inset)
consists of three joints: the innermost thoraco-coxal (TC-) joint enables
forward and backward movements, the middle coxa-trochanteral (CTr-) joint
enables elevation and depression of the leg, and the outermost femur-tibia
(FTi-) joint enables extension and flexion of the tibia. The morphology of
these multi-jointed legs were designed based on a cockroach leg (Zill et al.
2004). The front and back parts of the body are connected with a backbone
joint (BJ) which primarily allows upwards and downwards tilting of the front
body segment. Thus this is used for climbing and gap crossing purposes. More
details on BJ control for climbing can be found in (Goldschmidt et al. 2014).

The CPG-based control (see supplementary Figure 1 for detailed description)
itself is designed as a modular neural network that consists mainly of the
following four elements:

  1. 1.

CPG mechanism with neuromodulation for generating different rhythmic signals.
Inspired by biological findings, here the CPG circuit is designed as a two-
neuron fully connected recurrent network (Pasemann et al. 2003), such that
using different external neuromodulatory inputs different walking gaits can be
achieved.

  2. 2.

CPG post-processing units (PCPG) for shaping CPG output signals.

  3. 3.

Phase switching network (PSN) and velocity regulating networks (VRNs) for
walking directional control.

  4. 4.

Motor neurons with embedded fixed delay lines for transmitting motor commands
to all leg joints of AMOS II. These delay lines are utilized to realize the
inter-limb coordination, in which they introduce phase differences between the
transmitted signals to all leg joints. As a result, the desired walking gait
can be achieved.

The searching and elevation control at each leg, consist of single recurrent
neurons that receive the difference (instantaneous error) between the
predicted forward model signal and the actual sensory feedback. Due to the
recurrent self-connection, this error is accumulated over time. The
accumulated error can then be used to either extend specific leg joints in
order to get better foothold (searching action) during the stance phase, or
elevate further to overcome obstacles during the swing phase (see Figure 6 (e)
in section 4.1). All neurons in the CPG-based control and the searching and
elevation control are modeled as discrete-time rate-coded neurons with tan-
hyperbolic and piece-wise linear activation functions (see (Manoonpong et al.
2013) for details), respectively. They were updated with a frequency of
$\approx$ 27 Hz.

##  3 Materials & Method

###  3.1 Reservoir-based Distributed Adaptive Forward Models

We design, six identical adaptive RNN-based forward models
($RF_{1,2,3,...,6}$), one for each leg of the walking robot (Figure 2(a)).
These serve the purpose of online sensorimotor prediction as well as state
estimation. Specifically, each forward model learns to correctly transform the
efference copy of the actual motor signal for each leg joint (i.e., here the
CTr-joint motor signal111We use the CTr-joint motor signal instead of the TC-
and FTi-motor signals since this shows clear swing (off the ground) and stance
(on the ground) phases which can be qualitatively matched to the actual foot
contact signal.), into an expected or predicted sensory signal. This predicted
signal is then compared with the actual incoming sensory feedback signals
(i.e., here the foot contact signal - Figure 2 (b), of each leg) and, based on
the error accumulated over time, it triggers the appropriate action (searching
or elevation) and modulate the locomotive behavior of the robot. Each forward
model is based on a random RNN architecture of the self-adaptive reservoir
network type (Dasgupta et al. 2013), (Dasgupta 2015). Due to the presence of
rich recurrent feedback connections, the dynamic reservoir and intrinsic
homeostatic adaptations, the network exhibits a wide repertoire of nonlinear
activity and long fading memory. This can be primarily exploited for the
purpose of specific leg joint-motor signal transformation, act as motor memory
and for the prediction of sensorimotor patterns arising in the current
context.

Figure 2: (a) Neural mechanisms implemented on the bio-inspired hexapod robot
AMOSII. The yellow circle ($CPG$) represents the neural locomotion control
mechanism (see appendix. LABEL:sec:CPG). The gray circles ($RF_{1,2,3,...,6}$)
represent the reservoir-based adaptive forward models. The green circles
($SE_{1,2,3,...,6}$) represent searching and elevation control modules. The
orange circles represent leg joints where $TR_{i}$, $CR_{i}$, $FR_{i}$ are
TC-, CTr- and FTi-joints of the right front leg ($i=1$), right middle leg
($i=2$), right hind leg ($i=3$) and $TL_{i}$, $CL_{i}$, $FL_{i}$ are left
front leg ($i=1$), left middle leg ($i=2$), left hind leg ($i=3$),
respectively. $BJ$ is a backbone joint. The orange arrow lines indicate the
motor signals which are converted to joint angles for controlling motor
positions. The black arrow lines indicate error signals. The green arrow lines
indicate signals for adapting joint movements to deal with different
circumstances. b) An example of the reservoir-based adaptive forward model.
The dashed frame shows a zoomed in view of a single reservoir neuron. In this
setup, the input to each of the reservoir network comes from the CTr-joint of
the respective leg. The reservoir learns to produce the expected foot contact
signal for three different walking gaits ($z_{1}$, $z_{2}$, $z_{3}$). The
signals of the output neurons are combined and compared to the actual foot
contact sensory signal. The error from the comparison is transmitted to an
integrator unit. The unit accumulates the error over time. The accumulated
error is finally used to adapt joint movements through searching and elevation
control.

### Network Setup

The basic setup of each reservoir forward model can be divided into three
layers: input, hidden (or internal), and readout layers (Figure 2 (b)). The
internal layer consists of a large recurrent neural network driven by time-
varying stimuli (CPG motor signals). These driving signals are projected via
the input layer. The internal layer is constructed as a random RNN with fixed
randomly initialized synaptic connectivity (in this setup we only modify the
reservoir-to-readout neuron weights). Using a discrete time version of SARN,
with a step size of $\Delta t$, the discrete time state dynamics of each
reservoir neuron is given by the following equations:

${x}_{i}(t+1)=\left(1-\frac{\Delta t}{\tau_{i}}\right)x_{i}(t)+\frac{\Delta
t}{\tau_{i}}\left(g\sum_{j=1}^{N}W^{rec}_{i,j}r_{j}(t)+W^{in}_{i,1}u(t)+B_{i}\right),i=1,\dots,N.$
(1)

$r_{i}(t)=\mathit{tanh}(a_{i}x_{i}(t)+b_{i}),$ (2)

$\mathbf{z}(t)=\left[\mathbf{W}^{out}\right]^{T}\mathbf{r}(t).$ (3)

The RNN model consists of $N$ neurons, such that the membrane potential at the
soma (at time $t$) of the reservoir neurons, resulting from the incoming
excitatory and inhibitory synaptic inputs, is given by a $N$ dimensional
vector of neuron state activations. $x(t)=x_{1}(t),x_{2}(t),....,x_{N}(t)$.
The RNN here, does not explicitly model action potentials, but describes
neuronal firing rates. Where in, the variable $r_{i}(t)$ describes the
instantaneous firing rate ($N$ dimensional) of the reservoir neurons and is
calculated as a non-linear function of the state activation $x_{i}(t)$
(Equation 1). Each reservoir neuron $i$, receives inputs from other neurons in
the network with firing rates $r_{j}(t)$ via synaptic connections of strength
$W_{ij}^{rec}$ along with incoming stimuli from the input layer via synapses
of strength $W_{ij}^{in}$. Each reservoir neuron is also provided with an
auxiliary bias $B_{i}$. The parameter $g$ (Sompolinsky et al. 1988), (van
Vreeswijk and Sompolinsky 1996) acts as the scaling factor for the recurrent
connection weights allowing different dynamic regimes from stable ($g<1$) to
highly irregular chaotic ($g>1$) (Sussillo and Abbott 2009), being present in
the network.

Figure 3: (a) Plot of the change in the mean squared error for the forward
model task for one of the front legs ($R_{1}$) of the walking robot with
respect to the scaling of the recurrent layer synaptic weights $W^{rec}$ with
different $g$ values. As observed, very small values in $g$ have a negative
impact on performance compared with values closer to one being better.
Interestingly, the performance did not change significantly for $g>1.0$
(chaotic domain). This is mainly due to homeostasis introduced by intrinsic
plasticity in the network. The optimal value of $g=0.95$ selected for our
experiments is indicated with a dashed line. (b) Plot of the change in mean
squared error with respect to different reservoir sizes ($N$). $g$ was fixed
at the optimal value. Although increasing the reservoir size in general tends
to increase performance, a smaller size of $N=30$ gave the same level of
performance as $N=100$. Accordingly keeping in mind the trade off between
network size and learning performance, we set the forward model reservoir size
to 30 neurons. Results were averaged over 10 trials with different parameter
initializations on the forward model task for a single leg and a fixed walking
gait. (c) Example of the intrinsic plasticity to adjust the reservoir neuron
non-linearity parameters $a$ and $b$. Initially the the reservoir neuron fires
with an output distribution of Gaussian shape matching that of the input
distribution. However after adjustment using intrinsic plasticity mechanism
(Dasgupta et al. 2013) the reservoir neuron adapts the parameters $a$ and $b$,
such that, now for the same Gaussian input distribution the output
distribution follow a maximal entropy Exponential-like distribution. (d)
Distribution of the reservoir forward model individual neuron time constants
before and after adaptation.

The input to the reservoir $u(t)$, consists of a single CTr-joint motor
signal. This acts as an efference copy of the post-processed CPG motor output.
The readout layer consists of three neurons, with their activity being
represented by the three-dimensional vector $\mathbf{z}(t)$. Although
typically $M<N$ readout neurons can be connected to the reservoir, here we
restricted it to three neurons, as each readout here learns the predictive
signal for one of the following different walking gaits: wave ($z_{1}$),
tetrapod ($z_{2}$), and caterpillar ($z_{3}$) gaits. The wave, tetrapod, and
caterpillar gaits are used for climbing over an obstacle, walking on uneven
terrain, and crossing a large gap, respectively222These three gaits were
empirically selected among $19$ other possibilities. Previous studies have
demonstrated that the wave and tetrapod gaits are the most effective for
climbing and walking on uneven terrains, respectively. While in this
particular study we observed that the caterpillar gait was the most effective
for crossing a gap. However, without any loss of performance, additional
walking gaits can be applied easily by adding further readout neurons..
Subsequent to the supervised training of the reservoir-to-readout connections
$\mathbf{W}^{out}$, each readout neuron basically learns to predict the
expected foot contact signal associated with each of these gaits. The decay
rate for each reservoir neuron is given by $\frac{1}{\tau_{i}}$, where
$\tau_{i}$ is the individual membrane timeconstant. The input-to-reservoir
connections weights $\mathbf{W}^{in}$ and internal recurrent weights
$\mathbf{W}^{rec}$ were drawn randomly from the uniform distribution
$[-0.1,0.1]$ and a Gaussian distribution of zero mean and variance
$\frac{g^{2}}{\sqrt{p_{c}N}}$, respectively. Where, the parameter $p_{c}$
controls the probability of connections inside the recurrent layer and is set
to be 20%. In order to select the appropriate reservoir size, empirical
evaluations were carried out (Figures 3(a) and (b)), such that we achieved a
moderate network size of $N=30$, for which the minimum prediction error was
obtained at the readout layer, irrespective of the walking gait. The recurrent
weights were subsequently scaled by the factor of $g=0.95$ (see Figure 3),
such that the spontaneous network dynamics is in a stable regime and achieves
the best performance of the chosen network size. In accordance with the SARN
model, unsupervised intrinsic plasticity (Triesch 2005) and neuron timescale
adaptation (Dasgupta 2015) were carried out in order to learn the transfer
function parameters ($a_{i}$ and $b_{i}$)and the reservoir timeconstant
parameters $\tau_{i}$ for each individual neuron (Figure 3 (c) and (d)).

### Readout Weight Adaptation

Here we used a modified version of the original recursive least squares (RLS)
algorithm (Jaeger and Haas 2004),(Simon 2002) based on the FORCE learning
formulation (Sussillo and Abbott 2009), in order to learn the reservoir-to-
readout connection weights $\mathbf{W}^{out}$ at each time step, while the CPG
input ${u}(t)$ is being fed into the reservoir. The readout weights
$\mathbf{W}^{out}$ are calculated such that the overall error at the readout
neurons is minimized; thereby the network can learn to accurately transform
the CTr-motor signal to the expected foot contact signal, for each walking
gait. The instantaneous error signal ($e(t)$) at the readout layer, can be
calculated as the difference between the reservoir predicted output ($z(t)$)
and the desired output, $d(t)$ (i.e. here the expected foot contact signal).
Based on Equation 3, this can be formulated as:

$e(t)=\sum_{j=1}^{3}W_{j}^{out}(t-1)r_{j}(t)-d(t).$ (4)

Using the RLS algorithm, and minimizing this error, the readout weights
$W_{j}^{out}$ update can be defined by,

$W_{i}^{out}=W_{i}^{out}(t-1)-e(t)\sum_{j}P_{ij}(t)r_{j}(t).$ (5)

Where, $\mathbf{P}$ is a $N\times N$ square matrix proportional to the inverse
of the correlation matrix of the reservoir neuron firing rate vector
$\mathbf{r}$. $\mathbf{P}$ is initialized using the identity matrix
$\mathbf{I}$ and a small constant parameter $\delta_{c}$ as,
$\mathbf{P}(0)=\frac{\mathbf{I}}{\delta_{c}}$. $\mathbf{P}$, here, acts as the
adaptive learning rate for updating the readout weights with weight
modifications automatically slowing down as $\mathbf{P}$ decreases with time.
This allows the learning to occur stably and eventually converge to a
solution. $\mathbf{P}$ is updated as each time point as,

$\mathbf{P}(t)=\mathbf{P}(t-1)-\left(\frac{\mathbf{P}(t-1)\mathbf{r}(t)\mathbf{r}^{T}(t)\mathbf{P}(t-1)}{1+\mathbf{r}^{T}(t)\mathbf{P}(t-1)\mathbf{r}(t)}\right).$
(6)

The reservoir-to-readout neuron weights were initialized to zero at start.
Details of all the fixed parameters and initial settings for the reservoir
based forward model networks are summarized in Supplementary Table 1.

##  4 Results

###  4.1 Learning the Reservoir Forward Model (motor prediction)

In order to train the six forward models ($RF_{1}toRF_{6}$) in an online
manner, one for each leg, we let the simulated robot AMOSII walk under normal
conditions (i.e., walking on a flat terrain with the three different gaits).
Initially, we let the robot walk with a certain walking pattern, and then
every $2500$ time steps, the gait pattern was sequentially altered (this
occurs by changing the modulatory input to the CPG - see supplementary Figure
1). As a result, the robot sequentially transitions from wave gait, to
tetrapod gait, to caterpillar gait repeatedly. Using this procedure, we let
the robot walk for three complete cycles ($22500$ time steps) and collected
the corresponding CTr-motor signal and foot contact sensor readings for all
legs. Intrinsic plasticity and neuron time constant adaptations (Dasgupta et
al. 2013), (Dasgupta 2015), were then carried out using $20$ epochs of $1000$
time steps overlapping time windows. After this pre-training phase, all the
reservoir neuron non-linearity parameters and individual time constants
($\tau_{i}$) were fixed (see Figure 3 (d) for the distribution of neuronal
time constants before and after training).

Figure 4: Reservoir-to-readout weight adaptation during online learning. (a)
Changes of 30 weights projecting to the first readout neuron ($z_{1}$) of the
forward model of the right front leg ($R_{1}$) while walking with a wave gait.
During this period, weights projecting to the second ($z_{2}$) and third
($z_{3}$) output neurons remain unchanged (i.e., they are zero). (b) Changes
of the weights to $z_{2}$ while walking with a tetrapod gait. During this
period, the weights to $z_{3}$ still remain unchanged and the weights to
$z_{1}$ converge to around zero. (c) Changes of the weights to $z_{3}$ while
walking with a caterpillar gait. During this period, the weights to $z_{1}$
and $z_{2}$ converge to around zero. At the end of each gait, all weights are
stored such that they will be used for locomotion in different environments.
The grey areas represent transition phases from one gait to another gait and
the yellow areas represent convergence. The gait diagrams are shown on the
right. They are observed from the motor signals of the CTr-joints (Figure 5).
White areas indicate ground contact or stance phase and grey areas refer to no
ground contact during swing phase. As frequency increases, some legs step in
pairs (dashed enclosures). Here convergence implies no siginificant change in
the vector norm of the readout weights.

Subsequent to the pre-training phase, normal training of the reservoir-to-
readout weights $\mathbf{W}^{out}$ was carried out using the online RLS
learning algorithm with the same process of making the robot walk on a flat,
regular terrain and sequential switching between the three gait patterns every
$2500$ time steps. As such, at any given point in time only one of the readout
neurons (specific to the walking gait) are active. In this manner, synaptic
weights projecting from reservoir to the first readout neuron ($z_{1}$)
corresponding to the foot contact signal prediction for the wave gait, and
synaptic weights projecting to the second ($z_{2}$) and third ($z_{3}$)
readout neurons corresponding to the foot contact signal prediction of the
tetrapod and caterpillar gaits, are learned, respectively. Within this
experimental setup, as observed from Figures 4 (a), (b) and (c) the readout
weights corresponding to each gait converges very quickly (due to intrinsic
noise and nature of the reservoir-to-readout synaptic adaptation, the weights
still show minute fluctuations after successful learning; therefore here
convergence applies that the norm of the readout weights $|W^{out}|$ remains
constant with a small finite value (Sussillo and Abbott 2009)), in less than
the trial period of $2500$ time steps. As a result, every time the CTr-motor
signal changes due to walking gait transformations, the RF associated with
each leg learns to predict the expected foot contact signal robustly. The
training process was carried out only once under normal walking conditions.
This was subsequently used as the baseline in order to compare with the actual
foot contact signals (sensory feedback) while walking under the situations of
crossing a gap, climbing, and negotiating uneven terrains.

Figure 5 shows an example of the forward model prediction (training) during
the three different walking gaits, for the right front leg of AMOSII
($R_{1}$). Visual inspection clearly demonstrates that according to the
corresponding efference copy of CTr-motor signal at a particular gait, the
expected foot contact (FC) signal is precisely predicted at each time point.
Similarly, the foot contact signals for the other legs are also predicted
online, given the current context of CTr-signal (not shown). Note that the FC
signals of the other legs normally show slightly different periodic patterns.
Furthermore, there exists considerable lag between the expected stance phase
according to the motor signal and that observed from the FC signal (difference
between dotted green lines in Figure 5). Due to the internal memory of the
incoming motor signal in the reservoir, we see that the output neurons can
adapt to these time lags efficiently, even when the frequency of the signal
increases with a change in walking gaits. Furthermore, the reservoir-based
forward models enable the robust generation of the predicted FC signal, even
in the presence of high noise corruption or missing information in the
incoming CTr-joint motor signal (Figures 5 (j) and (k)). Due to the fact that
the CTr-motor signals are obtained after appropriate post-processing of
original CPG singals and passage through the motor neurons coupled with
different time delays. Such signal corruption can occur at various levels.
Therefore, the ability of the forward model to deal with such abrupt noise in
the motor signals in a robust manner is crucial to the adaptive mechanisms.
Furthermore, such signal corruptions can also occur, due to entrainment
mechanisms applied for the automatic tuning or adaptation of CPG outputs
(Nachstedt et al. 2013). Such online adaptation for sudden motor signal
variations, was not possible in the previous state of the art adaptive neuron
forward models (Manoonpong et al. 2013). This model inherently lacked the
ability to deal with variations in the temporal properties of the signal. As
such, a simple square wave matching the timing of the motor signal efference
copy was used, providing a limited range of behavior, as well as being
biologically implausible. However, here our reservoir-based model can
accurately estimate the spatiotemporal properties of the signal and robustly
learn the exact shape, as well as the timing of the actual FC signals.

Figure 5: (a-c) The CTr-joint motor signal of the right front leg ($R_{1}$)
for wave, tetrapod, and caterpillar gaits, respectively. This motor signal
provides the efference copy or the input to the reservoir forward models.
(d-f) The actual foot contact signal (force sensor signal under normal walking
conditions) used as the target signal of the reservoir models. (g-i) The
predicted foot contact signal or the final learned output of the forward model
for each walking gait ($RF$ output signal). The green shaded region indicates
the time interval between swing and stance phase for the CTr motor signal at
the three walking gaits. As observed the actual foot contact signal is
considerably lagged in time compared to the motor signal. Effectively, this
lag decreases with an increase in the gait frequency. The single RF adaptively
accounts for these different delay times in order to accurately predict the
expected foot contact signal. (j) above - CTr-joint motor signal demonstrated
for a single leg, with 2% Gaussian noise injected between 300-350 time steps
(yellow shaded region), below - Despite the noise corruption of the motor
signal, the reservoir forward model is able to generate the correct predicted
FC signal (blue dotted - target FC signal, red solid - predicted signal). (k)
above - The CTr-joint motor signal corrupted with missing information between
280-320 time steps. As a result, the motor signal shows a narrow spike between
310 -330 time steps (yellow shaded region), below - Reservoir forward model
predicted signal (red) as compared to the desired FC signal (dotted blue).
Although the CTr motor signal was transiently missing, the reservoir is able
to generate the desired FC signal considerably well, while at the same time
maintaining the correct temporal sequence of the signals.

###  4.2 Simulated Complex Environments

In order to assess the ability of the reservoir-based forward models to
generate adaptive333Forward models for motor prediction need an internal
fading memory of the motor apparatus, in order to adjust for time delays
between motor output signal and the actual sensory feedback (Kawato 1999).
complex locomotive behaviors in a neural closed-loop control system (see
Figure 1), we conducted simulation experiments under different situations
including crossing a gap, walking on uneven terrain and climbing over high
obstacles (similar to the behaviors observed in real insects). In all cases,
we used the same training procedure for the forward models by allowing the
robot to walk under normal conditions on a flat even terrain.

During testing of the learned behavior, while AMOSII walks under different
environmental conditions and a specific gait, the output of each trained
forward model (i.e., the predicted FC signal, Figure 6 (a)) is used to compare
it to the actual incoming FC signal of the leg (Figure 6 (b)). The difference
(instantaneous error signal $\Delta$) between them determines the walking
state where a positive value ($+\Delta$) indicates losing ground contact
during the stance phase and a negative value ($-\Delta$) indicates stepping on
or hitting obstacles during the swing phase.

$\Delta_{i}(t)=RF_{i}(t)-FC_{i}(t).$ (7)

where $i\in\\{1,2,...,6\\}$ represents each leg of the robot.

Thus, we use the positive value for searching control (Figure 6 (d) above).
This is then accumulated through a single recurrent neuron $S$ with a linear
transfer function and is always reset to 0.0 at the beginning of swing phase.
Similarly, the negative value is used for elevation control (Figure 6 (d)
below). The value is also accumulated through a recurrent neuron $E$ with a
linear transfer function. These accumulated errors (Figure 6 (c)) thus allow
the robot leg to be either elevated (on hitting an obstacle) or searching for
a foothold during the swing and stance phases respectively (see (Manoonpong et
al. 2013) for more details of the searching and elevation control). As
depicted in Figures 6 (a) and (b), while walking on a rough terrain (in this
case with tetrapod walking gait), the currently recorded sensory feedback or
foot contact sensor reading differs considerably from the reservoir predicted
signal. As a result, there is a high accumulation of error between each swing
or stance phase (Figure 6 (c)). It should be noted that the initial ($\approx
50$ time steps) abruptly high amplitude signal observed in the reservoir
forward model prediction, is caused due to the transient recovery time needed
by reservoir readout neurons to settle to the exact learned patterns. This is
overcome within the next few time steps and RF predicted FC signal continues
to occur in a robust manner. The accumulated error causes the corresponding
leg action control mechanism to kick in and the robot successfully navigates
out of the rough terrain (after $\approx 4000$ time steps). Once the robot
moves into the flat terrain, the reservoir predicted foot contact signal
matches almost perfectly with the actual sensory feedback. As a result, the
accumulated error becomes zero and normal walking without any additional
searching or elevation control mechanisms, can continue. In essence based on
the reservoir forward models, while traversing from the uneven terrain (Figure
6 inset 1-4) to the flat terrain (Figure 6 inset 5), the robot can adapt its
legs individually to deal with the change of terrain. That is, it depressed
its leg and extended its tibia to search for a foothold when loosing ground
contact during the stance phase. Losing ground contact information is detected
by a significant change of the accumulated errors (Figure 6 (c)). In case of
both walking on uneven terrain and climbing, this accumulated error causes
shifting of the CTr- and FTi-joints causing the respective leg to search for a
foothold. However, in the specific case of crossing a gap (Figure 7), we use
the accumulated error in order to control tilting of the backbone joint (BJ)
and shifting of the TC- and FTi-joints such that the front legs can be
extended forward continuously till the robot can find a foothold. In addition
to this leg joint control, reactive backbone joint control using the
additional ultrasonic sensors in front of the robot can also be used to learn
to lean up the BJ for climbing over obstacles (this has been previously
successfully applied using classical conditioning based learning in
(Goldschmidt et al. 2014) and as such not discussed here).

Figure 6: Successfully navigating rough terrain with reservoir forward model
(a) The reservoir forward model predicted, expected foot contact signal. After
a small initial transient the reservoir output quickly converges to the expect
signal for normal walking condition. (b) The actual sensory feedback (foot
contact signal) while walking on the rough surface (c) Accumulated error
calculated from the instantaneous error ($\Delta(t)$) after passing through
the recurrent neuron in the searching and elevation control . (d) The
searching and elevation action control system consisting of individual
recurrent neurons as signal accumulators. After 4000 time steps, the robot
successfully overcomes the rough terrain and continuous walking on a flat
surface. As a result, there is zero accumulated error since the predicted foot
contact signal almost exactly matches the actual signal. See the experiment
supplementary video 3.

We now take the example of the more complex, multiple gap crossing experiment
in order to look in detail at the learning outcome of the forward models. This
experiment was divided into two components, consisting of one larger gap
($15$cm length) and another relatively shorter gap of $11$ cm length. The two
gaps were separated by considerable distance where the robot was allowed to
walk on a regular flat terrain. In order to learn to cross a gap, we let AMOS
II walk with a caterpillar gait (see Figure 4 (c), right), such that each left
and right pair of legs moves simultaneously. Empirically this is observed to
be the most suited gait for overcoming large gaps, as well as supported by
experimental observations in stick insects (Blaesing and Cruse 2004). As shown
in Figure 7(1), at the beginning AMOS II walked forward straight towards the
initial gap. In this period, as it walks on the flat surface of the platform,
it performed regular movements similar to the training period under normal
walking conditions (training on a flat regular surface) . Eventually, it
encounters a $15$ cm wide gap ($\approx$ 44$\%$ of body length - the maximum
cross-able distance). In this situation, during the subsequent stance phase
the front legs of the robot loose ground contact (Figs. 7(d) and (e)). As a
result, the foot contact sensors from the front legs do not record any value.
However the reservoir forward model still predicts the expected foot contact
signal, causing a positive instantaneous error (Eq. 7). This leads to a
gradual ramping of the accumulated error signal between each stance phase and
swing phase, for the front legs (Figure 7 (a)).

In order to activate the BJ and adapt the leg movements due to the difference
between the reservoir predicted FC signal and the actual sensory feedback of
the FC sensors (error signals), we used the maximum accumulated error value of
the previous step (Figure 7, (a) red line) and control the BJ and leg
movements in the subsequent step. In this manner, the BJ started to lean
upwards incrementally (step like manner) at around $680-850$ time steps
(Figure 7(2)). Simultaneously, the TC- and FTi-joint movements of the left and
right front legs were also adapted accordingly in order to carry out elevation
action (this is reflected in the higher amplitude of these two signals in this
time period). Due to a predefined time-out period for tilting upwards, at
around $850$ time steps (Figure 7(3)), the backbone joint automatically moved
downwards recording a negative value. Consequently, the front legs touch the
ground of the second platform at the middle of the stance phase; thereby,
causing the accumulated error signals to decrease. Due to another time-out
period for tilting downwards at around $900$ time steps (Figure 7(4)), the BJ
automatically moved to the normal position ($-2\deg$). Since now the situation
is similar to walking on flat terrain, the RF predicted foot contact signal
matches the one recorded by the foot sensors, with accumulated error dropping
to zero. Thereafter, the TC- and FTi-joints perform regular movements.
Subsequently left and right hind legs loose the ground contact, and AMOSII
continues to walk forward. Here the movements of the TC- and FTi-joints were
slightly adapted allowing AMOS II to successfully cross the gap and continue
walking on the second platform (Figure 7(5)). As the terrian now resembles a
regular flat surface (similar to the original training terrain) AMOSII two
continues to walk forward in normal manner with no accumulated errors being
present. However, the same procedure is repeated once again, when AMOSII re-
encounters the second gap at around $2100$ time steps. However in this case,
since the gap length is much smaller, the elevation in the BJ occurs with an
initial increment of smaller amplitude (Figure 7 (2)) as compared to the
previous case. Thereafter, a similar process is followed and AMOSII can once
again successfully overcome this gap and continue walking on the other end of
the platform (Figure 7 (9)). This clearly demonstrates the adaptive yet robust
performance of the forward model based predictions in order to successively
cross gaps of different length.

Figure 7: Real-time data of walking and crossing multiple gaps using the
forward model predictions. (a) The accumulated error (black line) and the
maximum accumulated error value at the end of each stance phase (red line) of
the right front leg ($R_{1}$). The accumulated error is reset to zero every
swing phase. (b) The backbone joint (BJ) angle during walking and gap
crossing. The BJ stays at the normal position ($-2\deg$) during normal
walking. On encountering a gap ($15$cm), it leans upwards in a step like
fashion and then finally bent downwards in order to cross the gap. This
procedure is repeated for the second gap ($11$cm), however with different
degree of elevations. (c-e) The TC-, CTr-, and FTi-joint angles of right front
leg $R_{1}$ during normal walking and gap crossing. The joint adaptation was
controlled by the maximum accumulated error value of the previous step (red
line). Below pictures show snap shots of the locomotion of AMOS II during the
experiment. Note that one time step is $\approx 0.037$ s. For further details
interested readers are recommended to see the experiment supplementary videos
1 and 2.

Figure 8 shows that the reservoir forward model in combination with the neural
locomotion control mechanisms, not only successfully generates gap crossing
behavior of AMOS II (as shown above), and learns to walk on uneven terrain,
but also allows it to climb over single and multiple obstacles (eg. up a fleet
of stairs). In all these cases, we directly used the accumulated errors for
movement adaptation via the searching and elevation control mechanisms. For
climbing, the reactive backbone joint control was also applied to the system
(see (Goldschmidt et al. 2014) for more details) and a slow wave gait walking
pattern (see Figure 4 (a), right) was used.

Figure 8: Snapshots showing the learned behavior during climbing over a high
obstacle and climbing up a fleet of stairs (a) AMOSII walked with the wave
gait and approached a 15 cm high obstacle (1). It detected the obstacle using
its range sensors installed at its front part. The low-pass filtered range
sensory signals control the BJ to tilt upwards (2) and then back to its normal
position (3). Due to the missing foot contact of the front legs, the BJ moved
downwards to ensure stability (4). During climbing, middle and hind legs
lowered downwards due to the occurrence of the accumulated errors, showing leg
extension, to support the body. Finally, it successfully surmounted the high
obstacle (5). For further details see the supplementary experiment video 4 (b)
AMOSII climbed up a fleet of stairs (1-5) using the wave gait as well as the
reactive BJ control. The climbing behavior is also similar to the one
described in the case (a). For further details see supplementary experiment
video 5.

Experimentally the wave gait was found to be the most effective for climbing,
which allows AMOSII to overcome the highest climbable obstacle (i.e., 15 cm
height which equals $\approx 86\%$ of its leg length) and to surmount a fleet
of stairs. For walking on uneven terrain, a tetrapod gait (see Figure 4 (b),
right) was used without the backbone joint control. This is the most effective
gait for walking on uneven terrain (see also (Manoonpong et al. 2013)). Recall
that in all experiments the forward models basically generate the expected
foot contact signals (i.e., sensory prediction), which are compared to the
actual incoming ones. Errors between the expected and actual signals during
locomotion serve as state estimation and are used to adapt the joint movements
accordingly. It is important to note that, the best gait for each specific
scenario was experimentally determined and fixed. However, this could be
easily extended with learning mechanisms (see (Steingrube et al. 2010)) to
switch to the desired gait when the respective behavioral scenarios are
encountered, without any additional influence on the performance of the
reservoir forward models.

Figure 9: Average time to successfully overcome uneven terrains of different
elasticity (hard, moderate, highly elastic) (a) Average success time for
reservoir-based forward model. (b) Average success time for adaptive neuron
forward model from (Manoonpong et al. 2013). Here the whiskers indicate one
standard deviation above and below the mean value. Note the difference in
scale of the y-axis in both plots. The experimental surface here consisted of
the rough terrain as presented in Fig. 6 consisting of irregular undulations,
however with varying degree of elasticity for the three cases.

In order to evaluate the performance of our adaptive reservoir forward model
in comparison to the state of the art model recently presented in (Manoonpong
et al. 2013) (single recurrent neural with low-pass filter), we carried out
simulation experiments with AMOSII walking on different types of surfaces.
Specifically, after training on a flat surface (under normal conditions) we
carried out 10 trials each with the robot walking on uneven terrains (laid
with multiple obstacles of height $8cm$), having three different elastic
properties444Here the elasticity coefficients do not strictly represent
Young’s modulus values. These were local parameter setting defined in the
simulation, with increasing values causing greater elasticity.. The surfaces
were divided into hard ($1.0$), moderately elastic ($5.0$) and highly elastic
($10.0$). A tetrapod walking gait was used in all three cases. Starting from a
fixed position, we noted the total time taken by the robot to successfully
cross the uneven terrain region and move into a flat surface region. As
observed in Figs. 9 (a) and (b), the reservoir forward model enables the robot
to traverse the uneven region considerably faster as compared to the adaptive
neuron forward model, in all three scenarios. Both the models can be seen to
overcome the hard surface much better as compared to the elastic ones. This
was expected due to the changes in surface stiffness resulting in additional
forces on the robot legs. However, the reservoir model performance was
considerably more robust with a mean difference in success time of $1.86$ mins
for the hardest surface and approximately $2$ mins for the most elastic
surface, cases. Given that the walking gait was fixed, here the success time
can be thought as an indicator of the robot’s energy efficiency. In the
absence of additional body mechanisms to deal with changing surface stiffness,
the reservoir based model outperforms the previous implementations of adaptive
forward models by $\approx 25\%$ order of magnitude on average.

##  5 Discussion

In this study, we presented adaptive forward models using the self-adaptive
reservoir network for locomotion control. The model is implemented on each leg
of a simulated bio-inspired hexapod robot. It is trained online during walking
on a flat terrain in order to transform an efference copy (motor signal) into
an expected foot contact signal (i.e., sensory prediction). Afterwards, the
learned model of each leg is used to estimate walking states by comparing the
expected foot contact signal with the actual incoming one. The difference
between the expected and actual foot contact signals is used to adapt the
robot’s leg through elevation and searching control. Each leg is adapted
independently. This enables the robot to successfully walk on uneven terrains.
Moreover, using a backbone joint, the robot can also successfully cross a
large gap and climb over a high obstacle as well as up a fleet of stairs. In
this approach, basic walking patterns are generated by CPG-based control along
with local leg control mechanisms that make use of the reservoir prediction to
adapt the robot’s behavior. The key neural mechanisms presented in this work,
namely, CPG -based neural control, internal forward models and local leg
control, are essential for robust, adaptive locomotion control. However, only
individual instances of them has been successfully realized on artificial and
bio-mimetic robotic systems (Bläsing 2004), (Lewinger and Quinn 2011),
(Schilling et al. 2012), (Ren et al. 2012), (Christensen et al. 2014),
(Pfeifer et al. 2007); thereby achieving partial solutions. Furthermore,
although a few studies have focused on a combination of these neural
mechanisms, they have largely been tailored for adaptive locomotion in
quadruped robots (Lewis and Bekey 2002), (Silva et al. 2012), without the
ability to climb obstacles or cross large gaps, as observed in real animals
and insects. Thus, this work demonstrates how the combination of these
essential components, coupled with the power of the adaptive recurrent neural
forward models can achieve very rich behavioral repertoire in bio-inspired
hexapod robots. Thus supporting the idea that such embodied neural control
(Floreano et al. 2014) is indeed a potential powerful future alternative of
more conventional control methods.

It is important to note that the usage of reservoir networks, as forward
models here, provides the crucial benefit of an inherent representation of
time and fading memory (due to the internal feedback loops and input dependent
adaptations). Such memory of the time-varying motor or sensory stimuli is
required to overcome intrinsic time lags between expected sensory signals and
motor outputs (Wolpert et al. 1998), as well as in behavioral scenarios with
considerable dependence on the history of motor output (Lonini et al. 2009).
This is very difficult in most of the previous implementations of forward
internal models using either simple single recurrent neuron implementations
(Manoonpong et al. 2013), feed-forward multi-layered neural networks
(Schröder-Schetelig et al. 2010), or Bayesian network models (Dearden and
Demiris 2005), (Sturm et al. 2008). Furthermore, in this case, online
adaptation of only the reservoir-to-readout weights (readout) makes such
networks beneficial for simple and online learning.

The concept of forward models with efference copies in conjunction with neural
control has been suggested since the mid-20th century (Holst and Mittelstaedt
1950), (Held 1961) and increasingly employed for biological investigations
(Webb 2004). This is because it can explain mechanisms which biological
systems use to predict the consequence of their action based on sensory
information, resulting in adaptive and robust behaviors in a closed-loop
scenario. This concept also forms a major motivation for robots inspired by
biological systems. Within this context, the work presented here, verifies
that a combination of CPG-based neural control, adaptive reservoir forward
models with efference copies, and searching and elevation control can be used
for robustly generating complex locomotion and adaptive behaviors in an
artificial walking system. Additionally, although in this study we
specifically focused on locomotive behaviors for walking robots, (such) SARN
based motor prediction systems can be easily generalized to a number of other
applications. Specifically for neuro-prosthetic (Ganguly and Carmena 2009),
sensor-driven orthotic control (Braun et al. 2014), (Lee and Lee 2005) or
brain-machine interface devices (Golub et al. 2012), that require the learning
of such predictive models using highly non-stationary, temporal signals,
applying SARN models can provide high performance gains with embedded memory,
as compared to the current static feed-forward neural network solutions. In
the future, we will transfer the reservoir-based adaptive forward models to
the physical hexapod robot AMOS-II (Manoonpong et al. 2013) in order to test
the adaptive behaviors in a real environment. Furthermore, although in this
work, we specifically focused on a single CPG-based control mechanism, in the
future we plan to augment the distributed forward model architecture with
multiple CPG-based control (one for each leg) (Ren et al. 2015). Thereby,
truly enabling decentralized control of the robot legs for greater degree of
adaptation as observed in biology.

## Acknowledgement

This research was supported by the Emmy Noether Program (DFG, MA4464/3-1), the
Federal Ministry of Education and Research (BMBF) by a grant to the Bernstein
Center for Computational Neuroscience II Göttingen (01GQ1005A, project D1) and
the International Max Planck Research School for Physics of Biological and
Complex Systems scholarship.

## Author contributions

S.D, F.W and P.M designed the research. S.D and P.M implemented the model,
analyzed data and carried out simulations. D.G carried out the climbing
experiments. S.D and P.M wrote the manuscript.

## References

  * Blaesing and Cruse (2004) B. Blaesing and H. Cruse.  Stick insect locomotion in a complex environment: climbing over large gaps.  _Journal of experimental biology_ , 207(8):1273–1286, 2004.
  * Bläsing (2004) B. Bläsing.  Adaptive locomotion in a complex environment: simulation of stick insect gap crossing behaviour.  _From animals to animats_ , 8:173–182, 2004.
  * Braun et al. (2014) J. M. Braun, F. Wörgötter, and P. Manoonpong.  Internal models support specific gaits in orthotic devices.  In _Mobile Service Robotics_ , number 17 in Proceedings of the International Conference on Climbing and Walking Robots, pages 539–546, 2014\.
  * Christensen et al. (2014) D. J. Christensen, J. C. Larsen, and K. Stoy.  Fault-tolerant gait learning and morphology optimization of a polymorphic walking robot.  _Evolving Systems_ , 5(1):21–32, 2014.
  * Cruse (1976) H. Cruse.  The control of body position in the stick insect (carausius morosus), when walking over uneven surfaces.  _Biological Cybernetics_ , 24(1):25–33, 1976\.
  * Dasgupta (2015) S. Dasgupta.  _Temporal information processing and memory guided behaviors with recurrent neural networks_.  PhD thesis, Georg-August University, Göttingen, 2015.
  * Dasgupta et al. (2013) S. Dasgupta, F. Wörgötter, and P. Manoonpong.  Information dynamics based self-adaptive reservoir for delay temporal memory tasks.  _Evolving Systems_ , 4(4):235–249, 2013.
  * Dearden and Demiris (2005) A. Dearden and Y. Demiris.  Learning forward models for robots.  In _International Joint Conference on Artificial Intelligence_ , volume 5, page 1440, 2005.
  * Der and Martius (2012) R. Der and G. Martius.  The lpzrobots simulator.  In _The Playful Machine_ , pages 293–308. Springer, 2012.
  * Floreano et al. (2014) D. Floreano, A. J. Ijspeert, and S. Schaal.  Robotics and neuroscience.  _Current Biology_ , 24(18):910–920, 2014.
  * Ganguly and Carmena (2009) K. Ganguly and J. M. Carmena.  Emergence of a stable cortical map for neuroprosthetic control.  _PLoS biology_ , 7(7):e1000153, 2009.
  * Goldschmidt et al. (2014) D. Goldschmidt, F. Wörgötter, and P. Manoonpong.  Biologically-inspired adaptive obstacle negotiation behavior of hexapod robots.  _Frontiers in neurorobotics_ , 8, 2014.
  * Golub et al. (2012) M. D. Golub, B. Yu, and S. M. Chase.  Internal models engaged by brain-computer interface control.  In _Engineering in Medicine and Biology Society (EMBC), 2012 Annual International Conference of the IEEE_ , pages 1327–1330. IEEE, 2012.
  * Held (1961) R. Held.  Exposure-history as a factor in maintaining stability of perception and coordination.  _The Journal of nervous and mental disease_ , 132(1):26–hyhen, 1961.
  * Hesse et al. (2012) F. Hesse, G. Martius, P. Manoonpong, M. Biehl, and F. Wörgötter.  Modular robot control environment testing neural control on simulated and real robots.  In _Frontiers in Computational Neuroscience, Conference Abstract: Bernstein Conference_ , pages 1416–1420, 2012.
  * Holst and Mittelstaedt (1950) E. Holst and H. Mittelstaedt.  Das reafferenzprinzip.  _Naturwissenschaften_ , 37(20):464–476, 1950\.
  * Huston and Jayaraman (2011) S. J. Huston and V. Jayaraman.  Studying sensorimotor integration in insects.  _Current opinion in neurobiology_ , 21(4):527–534, 2011.
  * Jaeger and Haas (2004) H. Jaeger and H. Haas.  Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication.  _Science_ , 304(5667):78–80, 2004.
  * Kawato (1999) M. Kawato.  Internal models for motor control and trajectory planning.  _Current opinion in neurobiology_ , 9(6):718–727, 1999.
  * Kesper et al. (2013) P. Kesper, E. Grinke, F. Hesse, F. Wörgötter, and P. Manoonpong.  Obstacle/gap detection and terrain classification of walking robots based on a 2d laser range finder.  _Chapter_ , 53:419–426, 2013.
  * Lee and Lee (2005) J.-W. Lee and G.-K. Lee.  Gait angle prediction for lower limb orthotics and prostheses using an emg signal and neural networks.  _International Journal of Control, Automation, and Systems_ , 3(2):152–158, 2005.
  * Lewinger and Quinn (2011) W. A. Lewinger and R. D. Quinn.  Neurobiologically-based control system for an adaptively walking hexapod.  _Industrial Robot: An International Journal_ , 38(3):258–263, 2011.
  * Lewis and Bekey (2002) M. A. Lewis and G. A. Bekey.  Gait adaptation in a quadruped robot.  _Autonomous robots_ , 12(3):301–312, 2002.
  * Lonini et al. (2009) L. Lonini, L. Dipietro, L. Zollo, E. Guglielmelli, and H. I. Krebs.  An internal model for acquisition and retention of motor learning during arm reaching.  _Neural computation_ , 21(7):2009–2027, 2009\.
  * Maass et al. (2002) W. Maass, T. Natschlaeger, and H. Markram.  Real-time computing without stable states: A new framework for neural computation based on perturbations.  _Neural computation_ , 14(11):2531–2560, 2002\.
  * Manoonpong et al. (2013) P. Manoonpong, U. Parlitz, and F. Wörgötter.  Neural control and adaptive neural forward models for insect-like, energy-efficient, and adaptable locomotion of walking machines.  _Frontiers in neural circuits_ , 7, 2013.
  * Mischiati et al. (2015) M. Mischiati, H.-T. Lin, P. Herold, E. Imler, R. Olberg, and A. Leonardo.  Internal models direct dragonfly interception steering.  _Nature_ , 517(7534):333–338, 2015.
  * Nachstedt et al. (2013) T. Nachstedt, F. Worgotter, P. Manoonpong, R. Ariizumi, Y. Ambe, and F. Matsuno.  Adaptive neural oscillators with synaptic plasticity for locomotion control of a snake-like robot with screw-drive mechanism.  In _Robotics and Automation (ICRA), 2013 IEEE International Conference on_ , pages 3389–3395. IEEE, 2013.
  * Pasemann et al. (2003) F. Pasemann, M. Hild, and K. Zahedi.  So (2)-networks as neural oscillators.  In _Computational methods in neural modeling_ , pages 144–151. Springer, 2003.
  * Pearson and Franklin (1984) K. Pearson and R. Franklin.  Characteristics of leg movements and patterns of coordination in locusts walking on rough terrain.  _The International Journal of Robotics Research_ , 3(2):101–112, 1984.
  * Pfeifer et al. (2007) R. Pfeifer, M. Lungarella, and F. Iida.  Self-organization, embodiment, and biologically inspired robotics.  _science_ , 318(5853):1088–1093, 2007.
  * Ren et al. (2012) G. Ren, W. Chen, C. Kolodziejski, F. Worgotter, S. Dasgupta, and P. Manoonpong.  Multiple chaotic central pattern generators for locomotion generation and leg damage compensation in a hexapod robot.  In _Intelligent Robots and Systems (IROS), 2012 IEEE/RSJ International Conference on_ , pages 2756–2761. IEEE, 2012.
  * Ren et al. (2015) G. Ren, W. Chen, S. Dasgupta, C. Kolodziejski, F. Wörgötter, and P. Manoonpong.  Multiple chaotic central pattern generators with learning for legged locomotion and malfunction compensation.  _Information Sciences_ , 294:666–682, 2015.
  * Schilling et al. (2012) M. Schilling, J. Paskarbeit, J. Schmitz, A. Schneider, and H. Cruse.  Grounding an internal body model of a hexapod walker control of curve walking in a biologically inspired robot.  In _Intelligent Robots and Systems (IROS), 2012 IEEE/RSJ International Conference on_ , pages 2762–2768. IEEE, 2012.
  * Schröder-Schetelig et al. (2010) J. Schröder-Schetelig, P. Manoonpong, and F. Wörgötter.  Using efference copy and a forward internal model for adaptive biped walking.  _Autonomous Robots_ , 29(3-4):357–366, 2010.
  * Silva et al. (2012) P. Silva, V. Matos, and C. P. Santos.  Adaptive quadruped locomotion: learning to detect and avoid an obstacle.  In _From Animals to Animats 12_ , pages 361–370. Springer, 2012.
  * Simon (2002) H. Simon.  Adaptive filter theory.  _Prentice Hall_ , 2:478–481, 2002.
  * Sompolinsky et al. (1988) H. Sompolinsky, A. Crisanti, and H. Sommers.  Chaos in random neural networks.  _Physical Review Letters_ , 61(3):259, 1988.
  * Steingrube et al. (2010) S. Steingrube, M. Timme, F. Wörgötter, and P. Manoonpong.  Self-organized adaptation of a simple neural circuit enables complex robot behaviour.  _Nature Physics_ , 6(3):224–230, 2010.
  * Sturm et al. (2008) J. Sturm, C. Plagemann, and W. Burgard.  Adaptive body scheme models for robust robotic manipulation.  In _Robotics: Science and systems_ , 2008.
  * Sussillo and Abbott (2009) D. Sussillo and L. F. Abbott.  Generating coherent patterns of activity from chaotic neural networks.  _Neuron_ , 63(4):544–557, 2009.
  * Triesch (2005) J. Triesch.  A gradient rule for the plasticity of a neuron’s intrinsic excitability.  In _Artificial Neural Networks: Biological Inspirations–ICANN 2005_ , pages 65–70. Springer, 2005.
  * van Vreeswijk and Sompolinsky (1996) C. van Vreeswijk and H. Sompolinsky.  Chaos in neuronal networks with balanced excitatory and inhibitory activity.  _Science_ , 274(5293):1724–1726, 1996.
  * Watson et al. (2002) J. T. Watson, R. E. Ritzmann, S. N. Zill, and A. J. Pollack.  Control of obstacle climbing in the cockroach, blaberus discoidalis. i. kinematics.  _Journal of Comparative Physiology A_ , 188(1):39–53, 2002.
  * Webb (2004) B. Webb.  Neural mechanisms for prediction: do insects have forward models?  _Trends in neurosciences_ , 27(5):278–282, 2004\.
  * Wolpert et al. (1998) D. M. Wolpert, R. C. Miall, and M. Kawato.  Internal models in the cerebellum.  _Trends in cognitive sciences_ , 2(9):338–347, 1998.
  * Zenker et al. (2013) S. Zenker, E. E. Aksoy, D. Goldschmidt, F. Wörgötter, and P. Manoonpong.  Visual terrain classification for selecting energy efficient gaits of a hexapod robot.  In _Advanced Intelligent Mechatronics (AIM), 2013 IEEE/ASME International Conference on_ , pages 577–584. IEEE, 2013.
  * Zill et al. (2004) S. Zill, J. Schmitz, and A. Büschges.  Load sensing and control of posture and locomotion.  _Arthropod structure & development_, 33(3):273–286, 2004.

See pages - of SupplementaryData.pdf