The way neurotransmitters are released through the synaptic cleft may change transiently upon repeated presynaptic activity (Tsodyks and Markram, 1997), both for excitatory (Wang et al., 2006) and for inhibitory (Gupta et al., 2000) synapses. Physiologically this is, on the one side due to an increase of the Ca-concentration u ∈ [1, U_max] within the presynaptic bulge, facilitating the release of the respective neurotransmitter, and, on the other side, due to the decrease of the number φ ∈ [0, 1] of available vesicles of neurotransmitters. We use here with

a modified version of the original Tsodyks-Markram model (Tsodyks and Markram, 1997; Hennig, 2013), in which the the Ca-concentration u and the number of vesicles φ of a given synapse relax to target values U = U(y) and Φ = Φ(u, y), determined in turn by the level y ∈ [0, 1] of the presynaptic activity. A prolonged maximal presynaptic activity y ≡ 1 would lead with φ → 0 to a full depletion of the reservoir of vesicles.

The dynamics of the full depletion model (2.1) is determined by the relaxation time constants T_u and T_φ, and by the maximal level U_max of the Ca concentration. For U_max = 1 a monotone depression is present, whereas U_max > 1 initially generates facilitation by a fast calcium influx, being annulled later on by the depletion of neurotransmitters. Overall, the synaptic efficiency is proportional to uφ, viz to the number of vesicles and to the release probability (which in turn is assumed to be proportional to u). We use T_u = 300ms and T_φ = 600ms, together with either U_max = 1 or U_max = 4. These values are within the typical range of what is physiologically observed (Gupta et al., 2000; Wang et al., 2006).

The movement of robot illustrated in Figure 1 is induced by the relative gravitational pull of the three weights, together with the rolling friction and angular momentum conservation. The individual neurons i = 1, 2, 3 are modeled as rate-encoding leaky integrators,

where x_i and y(x_i) are the respective membrane potentials and firing rates. Γ is the relaxation rate, R the diameter of the robot, p ∈ [0, 1] a rescaling factor, xi(a)∈[-R,R] the sensory reading of the actual position of the weight on the rod, w₀ > 0 the weight of excitatory input and z₀ > 0 the magnitude of the inter-neural inhibitory connections. We note that the variables of the STSP, u_j and φ_j, as described by Equation (2.1), depend only on the presynaptic activity and can hence be attributed altogether to the presynaptic neuron. For the slope of the sigmoidal a = 0.4 has been selected. The weight of the excitatory input w₀ is not modulated here by short-term synaptic plasticity, corresponding to a direct sensory reading.

We selected with p = 1/2 a reduced range for the target position xi(t),

This choice allows to avoid dynamic overshooting of the weight when accelerated from its actual to the target position. The force accelerating the weight is calculated by the LPZRobots package by simulating a damped spring:

where k is the spring constant and γ the damping. Centrifugal and other induced forces, F_i, act additionally in Equation (4) on the individual weights. The complete setup of the three-neuron network is illustrated in Figure 2.

The LPZRobots simulation environment (Der and Martius, 2012) is an interactive simulator based on the ODE (Open Dynamic Engine) (Smith, 2005). LPZRobots contains rigid body dynamics in terms of a library of basic primitive objects, such as spheres and cuboids, as well as a variety of joints, sensors and surface materials.

We used roughness = 0.8, slip = 0.01, hardness = 40 and elasiticity = 0.5 for the collision and friction properties together with friction = 0.3 (the rolling friction coefficient), gravity = −9.81 (the gravitational constant) and noise = 0 (for the global noise level). All parameters are in SI units. For the stepsize of the physical simulation simstepsize = 0.001 was used (corresponding to a millisecond). With controlinterval = 1 one ensures that the controller, viz Equation (2), is updated as often as the physics of the environment.

The robot itself has a diameter of 2R = 0.5, a mass off M = 1 and a motorpowerfactor = 120. The parameters for the damped oscillator (Equation 4) are m = 1, k = m^*motorpowerfactor and γ=2k*m (critical damping). The relaxation rate for the membrane potential entering Equation (2) has been set to Γ = 20, retaining the bare excitatory and inhibitory weights, w₀ and z₀, as free simulation parameters.

In Figure 3 we present the stability regions for the various regular movement patterns found, with respective close-ups given in Figure 4. The results are for U_max = 1 (depressing short-term synaptic plasticity without Ca dynamics) and for the parameters specified in Section 2.3. They are obtained by adiabatically continuing stable states along a grid until stability is lost. Without STSP only a globally attracting fixpoint corresponding to a motionless robot is present. We note that regular motion arises for a wide range of bare excitatory (w₀) and inhibitory (z₀) synaptic weights. z₀ needs however to be larger than w₀.

All motion patterns observed are self-organized. There is no objective function (Gros, 2014), such as a maximal velocity, to be optimized. This implies that the quantitative features of the individual motion patterns change smoothly within their respective stability regions, and that one can identify the observed regular movement patters as stable limit cycles in the sensorimotor loop (Sándor et al., 2015). Fast switching between motion primitives would be possible by a putative overarching controller, since more than one limit cycle may be stable for given synaptic weights w₀ and z₀. Interactions between robots or with external obstacles might also lead to the automatic selection of another coexisting mode (see for instance Supplementary Video 1).

It is evident that the body plan of the robot examined here tends to produce meandering motion pattern. T1 and T2 are sun- and star-like movements with small (T1) and large (T2) processing angles (compare Figure 4; “T” stands for torus in phase space). There is, in addition, a (nearly pure) circular motion, C1, and three types of forward snake-like meandering motion patters, S1, S2, and S3. From these S3 partly overlaps with itself. These modes are characterized by distinct motion patterns of the three weights, as shown in Figure 5, as measured by their positions along their respective rods. The differences between the distinct modes are in part qualitative, in terms of the time sequences in which the three neurons are subsequently active, and in part only quantitative. The difference between T1 and S1 is, in this respect, that the up-times of the two active neurons are symmetric for S1, but not for T1. A spontaneous symmetry breaking can be furthermore observed in case of T1, S1, S2, S3, for which two weights always have alternating dynamics, the third one showing a qualitatively different behavior. In contrast to that, the time-series of the C1 and T2 modes reveals the symmetrical but phase shifted oscillation of the three weights. Note that the positions of the weights may overshoot the interval [−pR, pR] for the target positions xi(t), both due to inertia and due to the additional gravitational pull. Motion patterns similar to the ones shown in Figure 4 have been observed in a self-organized two-wheeled robot in the frozen mode (Der and Martius, 2013).

The dynamics of the robot takes place in a phase space combining the internal variables, of both body and controller, with the ones of the environment. The stability regions of the individual limit cycles presented in Figure 3 will therefore be bounded, generically, by a suitable bifurcation, such as a supercritical Hopf bifurcation or a fold bifurcation of limit cycles (Gros, 2015; Sándor et al., 2015). Alternatively, a transition to chaos may occur. It is on the other side also possible that chaotic attractors emerge from previously unstable manifolds and that the stability region of chaotic and stable manifolds overlap.

Close to a chaotic phase long transients may occur, which makes it difficult to study systematically the exact extend of the chaotic region. In Figure 3 we have indicated however a few representative combinations of parameters, for which stable chaos is observed both in the limit of long simulations times and for a wide range of stepsizes of the ODE simulator. No regular motion patterns can be observed in the screenshots presented in Figure 6. We have also evaluated the long-time behavior of the square of the covered real-space distance,

We found diffusive transport d~τ for the chaotic mode and ballistic transport d ~ τ for the forward meandering modes S1, S2, and S3. Both as expected.

It has been observed, that chaotic locomotion of an embodied system may be considered as a basic explorative behavior, both of the environment and of the own motor pattern (Steingrube et al., 2010; Shim and Husbands, 2012). As a test of this hypothesis we have set our three-rod robot into a restricted playground containing movable objects in the form of blocks, which can be pushed, to a certain extend, over the ground. A screenshot is presented in Figure 6. One can observe, that the robot stays for a while close to the object, bumping around, and retracting in part a trajectory having a shape similar to the one generated by a C1 limit cycle. This is possible, as the set of parameters (w₀, z₀) = (210, 400) considered is located close to but outside the C1-stability region. The C1 limit cycle is hence only weakly unstable in the chaotic phase. The active exploration of the environment, occurring here when bumping into obstacles, gives the robot hence access to otherwise unstable locomotion options. The overall behavior may be interpreted alternatively in terms of non-representational sensorimotor knowledge (Buhrmann and Di Paolo, 2014). For a longer simulation see the Supplementary Videos.

In the movie presented in the Supplementary Material one can observe, furthermore, that the robot is pushing the blocks around in a seemingly “playful” manner (see Supplementary Video 3). A remarkable behavior, in our view, considering that the sphere robot disposes of a mere total of three controlling neurons. We note, that this complex behavior results from the interplay of the autonomous dynamics, as resulting from the inter-neural short-term synaptic plasticity, with environmental feedback.

One can consider the controlling 3-neuron network in isolation by identifying the sensory reading xi(a) for the actual position of the weight along the rod with the respective target position xi(t), viz by setting xi(a)=xi(t) in Equation (2). The resulting network contains a self-excitatory coupling w₀ together with all-to-all inhibition with a bare synaptic strength z₀. The short-term synaptic plasticity then induces an autonomous activity, as illustrated in Figure 7, which is topologically equivalent to the C1 mode. This equivalence becomes even more pronounced when suspending the robot in air, which can be achieved in turn by simply removing gravity from the physics simulation (bottom time-series in Figure 7). One can hence consider the C1 mode as the driver for the observed physical motion.

The isolated 3-neuron network has, however, only a single stable limit cycle. Numerically integrating the isolated network for parameters settings (w₀, z₀) corresponding to the six modes of Figure 5, as well as for chaotic states, we find always an identical sequential activation of the three neurons illustrated in Figure 7, with only slight changes in the overall shape. It is hence clear, that the other modes T1, T2, S1, S2, and S3, as well as the chaotic behavior, do result from the closed-loop feedback of the environment. The interaction of the environment with the intrinsic dynamics then results in the emergence of alternative types of locomotion.

We present in Figure 8 an analysis of the stability of the various modes found, with respect to noise in the sensory readings, where the level of the noise is given by the relative standard deviation σ of the sensory readings xi(a). Comparing with the phase diagram, as presented in Figure 3, one notices that first modes to disappear, T1 and S3, are the ones with small stability regions in the phase diagram. Ramping up the noise level the T1 and S3 modes turn respectively, above their corresponding critical noise levels, into C1 and S1 modes. The other modes, including the chaotic phase, are in contrast very stable with respect to noise.

We present in Figure 9 the phase diagram obtained when using U_max = 4 for the maximal Ca-level entering Equation (2.1). Within the range of (w₀, z₀) scanned we find four out of the six modes observed for U_max = 1 (compare Figure 3). The range of inhibitory weights z₀ for which stable locomotion is found is rescaled down, in addition, with respect to the U_max = 1 case. Interestingly we found a chaotic state at (180, 80) which lies just inside the stability region of the C1 mode.

We did let the robot evolve within the borders of a simple maze, as shown in Figure 10 and Supplementary Video 4. Most of the time the robot is in the chaotic state, which is the dominant mode for the parameters used, namely (w₀, z₀) = (180, 80) and U_max = 4. Intermittently, after colliding with a wall, the robot switches to the coexisting C1 mode. The radius of the stable C1 limit cycle in real-world coordinates is however so large, for (w₀, z₀) = (180, 80), that it does not fit into the maze. The robot hence continues exploring. We have obtained similar results when using a U_max = 1 chaotic mode.

A screenshot of a trajectory in open space is presented in Figure 11. One notices, that the U_max = 4 and (w₀, z₀) = (180, 80) chaotic mode wanders around aimlessly in much smother manner, than the U_max = 1 chaotic mode shown in Figure 6. This is the result of topologically different attractor structures, as seen in the phase space of internal variables (see the Supplementary Materials). Different types of chaos are indeed known to exist (Wernecke et al., 2016).

The autonomous mode switching observed for the regular motion primitives can also be seen in Supplementary Video 1. For a detailed discussion of the possible switching scenarios see the Supplementary Materials.

In Figure 12 we compare for the two chaotic modes, realized for U_max = 1 and for U_max = 4 respectively, the time series for the positions of the weights along the rods. One observes, that the movements of the weight is qualitatively similar, on short time scales, to an S2 mode (compare Figure 5, see also Supplementary Video 3). It is interesting, in this context, that the S2 mode has two types of degeneracies.

For the U_max = 4 chaotic mode we did not observe discrete mode switching, in above sense, which however occurs frequently for the U_max = 1 mode (see Figure 12). The chaotic meandering observed for the U_max = 4 chaotic mode, as evident in Figure 11, is hence a consequence of a smooth diffusion of the angle of propagation on the manifold of unstable S2 limit cycles (or limit cycle relicts Linkerhand and Gros, 2013). In the phase space of the neural activity (as shown in Supplementary Figure 5), the trajectory corresponds to a chaotic phase diffusion along a limit cycle (Wernecke et al., 2016). This process is determinstic and not due to numerical errors, as we have checked by systematically reducting the stepsize used for the numerical integration. Noise is absent.

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.