## INTRODUCTION

Designing a cognitive architecture that acts spontaneously in a real-world environment is one of the ultimate goals in the field of cognitive robotics (*1*). A cognitive agent is expected to have autonomy, i.e., the agent should behave independently of the designer’s control while maintaining its identity. Adaptability is another requirement for cognitive functionality. Thus, the agent must select the appropriate behavior continuously and robustly in response to the changing environment in real time. To summarize, the agent’s cognitive behavior should be implemented through the body-environment interaction while still enabling the agent to maintain its autonomy and adaptability.

In the conventional context of robotics and artificial intelligence, designers often take top-down approaches to provide an agent with a hierarchical structure corresponding to the behavioral category. This representation-based approach has a critical limitation in the design of a cognitive agent: The static, one-to-one relationship between the behavior and structure makes it difficult to adapt flexibly to the dynamically changing environment and developing body. For example, it has been considered that the motion control systems of living things, including humans, realize their high-order motion plans by combining reproducible motor patterns called motion primitives (*2*). Inspired by this viewpoint, we tend to realize an agent’s behavior control with a predetermined static hierarchical structure. However, such a hierarchical structure does not exist in living organisms from the beginning; rather, these structures are cultivated through the body’s development and dynamic interactions with the environment. Therefore, it is important to introduce a dynamical perspective to understand a hierarchical structure of behavior control generated in animals that have adaptability and flexible plasticity.

In robotics, approaches based on dynamical systems theory have been applied to analyze and control agents being modeled as sets of variables and parameters on a phase space (*3*, *4*). This dynamical systems approach can deal with both the functional hierarchy and the elementary motion in a unified form by expressing the physical constraints of the agent as the temporal development of state variables, namely, dynamics. For example, Jaeger (*4*) sketched a pioneering idea of an algorithm where the behavior of an agent is expressed as dynamics, and both the behavioral regularity (referred to as transient attractors) and the higher-order relationships among them are extracted in a bottom-up manner. Unlike the stochastic approach where the randomness of the system is realized by a probabilistic model (e.g., the Markov model), the dynamical system approach can consistently express the agent’s seemingly random behaviors via its chaoticity (*5*). Furthermore, in the stochastic approach, hierarchical structures are inevitably introduced since the mechanism of random number generation is completely independent of the system’s dynamics. Thus, the dynamical systems approach has the potential to model an agent’s spontaneous hierarchical behavior in the form of dynamic interaction without top-down structure given by the external designer. The important issue here is that the concept of the dynamical system approach itself offers no general principle for implementing these behaviors on a large number of nonlinear couplings that constitute the body-environment interaction. Therefore, it is important to study and propose the methodologies for designing spontaneous hierarchical behavior with a consistent temporal evolution rule governing the system’s dynamics.

Following this dynamical systems perspective, chaotic itinerancy (CI) (*6*–*8*) is a powerful option for modeling spontaneous behavior with the functional hierarchy realized through the dynamics. CI is a frequently observed, nonlinear phenomenon in high-dimensional dynamical systems, and it is characterized by chaotically itinerant transitions among locally contracting domains, namely, quasi-attractors (*8*). In general, a chaotic system has an initial sensitivity, and a slight difference in phase space is exponentially expanded in a certain direction with temporal development. Conversely, multiple transiently predictable dynamics can be repeatedly observed in a chaotic system, yielding CI despite the global chaoticity and initial sensitivity. This type of hierarchical dynamics frequently emerges from high-dimensional chaos even without hierarchical mechanisms, implying that explicit structure is not necessarily needed for implementing hierarchical behaviors. The chaoticity plays an important role in forming the autonomy of an agent as it is virtually impossible for the designer to completely predict and control an agent’s behavior because of the agent’s initial sensitivities, which essentially ensures the agent’s independence from the designer. Thus, CI would work as an effective tool for implementing the intellectual behavior of a cognitive agent by embedding the behavior in the form of a quasi-attractor and maintaining the autonomy of the agent with the chaoticity.

CI was first found in a model of optical turbulence (*6*). Since its discovery, similar phenomena have been numerically obtained in various setups (*7*, *9*, *10*). The properties of CI vary among previous studies, and some classes of CI present interesting features that are difficult to characterize using the stochastic processes. Tsuda *et al*. (*9*), for example, show that their asynchronous neural network model produces itinerant behavior whose transition frequency is characterized by a long-tailed distribution. The coupled map lattice proposed by Kaneko (*7*) potentially emits infinitely many states, whose transition rule should be represented by an infinite state machine. Several physiological studies have reported that CI-like dynamics have even occurred in brain activity, suggesting that CI might play an essential role in forming cognitive functions (*11*, *12*). For example, Freeman (*13*) revealed that an irregular transition among learned states was observed in the electroencephalogram pattern of a rabbit olfactory system when a novel input was given, indicating that the cognitive conditions corresponding to “I don’t know” are internally realized as CI-like dynamics. Furthermore, a recent observation of rat auditory cortex cell activity revealed the existence of a random shift among different stereotypical activities corresponding to individual external stimuli during anesthesia (*14*). On the basis of these reports, Kurikawa and Kaneko (*15*) suggested the novel idea of “memory-as-bifurcation” to understand the mechanism of the transitory phenomenon. They reproduced it in an associative memory model in which several input-output functions were embedded by Hebbian learning. An intermittent switching among metastable patterns was also observed in a recurrent neural network (RNN) by installing multiple feedback loops trained to output a specific transient dynamics corresponding to external transient inputs (*16*). CI-like dynamics arise not only in nervous systems but also in interactions between agent’s bodies and their surrounding environments (*17*–*19*).

For example, Kuniyoshi and Sangawa (*17*) developed a human fetal development model by coupling chaotic central pattern generators and a musculoskeletal system. They reported that several common behaviors, such as crawling and rolling over, emerged from the physical constraint. Therefore, CI is a nonlinear phenomenon of high-dimensional dynamical systems and is thought to play a substantial role in generating structural behavior.

Inspired by the contribution of CI to the cognitive functions and spontaneous motion generation of agents, CI has been used for motion control in the field of neurorobotics and cognitive robotics by designing the CI trajectory. For example, Namikawa *et al.* (*20*–*22*) designed stochastic motion switching among predetermined motion primitives in a humanoid robot by using a hierarchical, deterministic RNN controller. In this study, it was confirmed that lower-order RNNs with smaller time constants stably produced the trajectories of motion primitives, whereas higher-order RNNs with larger time constants realized a pseudo-stochastic transition by exploiting self-organized chaoticity. Steingrube *et al*. (*23*) designed a robot that skillfully broke a deadlock state in which the motion had completely stopped by using chaos in the RNN controller. Hence, it can be interpreted that CI-like dynamics were embedded in the coupling of the body and the surrounding environment.

While CI is an important phenomenon in high-dimensional dynamical systems, roboticists also find it a useful tool for designing an agent’s behavior structure while maintaining the agent’s autonomy. However, it has generally been difficult to embed desired quasi-attractors at will because of their nonlinearity and high dimensionality. For example, in the method of Namikawa *et al.* (*20*–*22*), the internal connections of an RNN ware trained with backpropagation through time (*24*); however, embedding a long-term input-output function in an RNN by the gradient descent method is generally unstable and requires a large number of learning epochs (*25*). Furthermore, their method required both a hierarchical structure and the same number of separated modules as the motion primitives, restricting the scalability and the range of its application. Yamashita and Tani (*26*) proposed a network model learning functional hierarchy without any modules, in which hierarchical behavior was implemented by the explicit hierarchical structure governed by predetermined multiple time scales. In addition, methods using the associated memory model (*10*, *27*–*29*) are also unsuitable for our purpose since it is difficult to embed the quasi-attractors with complicated spatiotemporal patterns.

In this study, we propose an algorithm, freely designing both the trajectories of quasi-attractors and transition rules among them in a setup of high-dimensional chaotic dynamical systems. We aim to design the properties of CI characterized by a finite state machine and finite switching time, such as those in the neurorobotics context (*19*, *21*). We prepare transition rules described by a Markov model and aim to emulate them through CI using a high-dimensional nonlinear dynamical system. Our method uses batch learning composed of the following three-step procedure (Fig. 1):

Step 1. Prepare a high-dimensional chaotic system where target quasi-attractors are embedded. We used a widely used echo state network (ESN) (*30*), one type of RNN, as a high-dimensional chaotic system. This ESN contains no hierarchical structure and modules (e.g., multiple time scales), and every network node shares the same time scale parameter. At the same time, modify the interactions (internal parameters) so that the system reproducibly generates intrinsic complex trajectories generated by an initial chaotic system (innate trajectories) corresponding to the type of the discrete inputs (named symbol). In parallel, train the linear regression model (named readout) to output the designated trajectories (output dynamics) by exploiting the embedded innate trajectory. This process can be potentially applied to the other chaotic dynamical systems not limited to RNN in silico (*31*) since neither modules nor hierarchical structures are required. In addition, this embedding process is accomplished by modifying fewer parameters using the method of reservoir computing (*32*, *33*). Therefore, our scheme is more stable and less computationally expensive than conventional methods using backpropagation to train the network parameters.

Step 2. Add a feedback classifier to the trained chaotic systems for autonomously generating specific symbolic dynamics. In the training of the feedback discriminator, the network’s internal parameters are fixed, as with the readout in step 1. Thus, by using the embedded innate trajectory, the feedback discriminator achieves multiple symbol transition rules with minimum additional computational capacity (i.e., nonlinearity and memory).

Step 3. Regulate the feedback unit added in step 2 to design designated stochastic symbol transition rules. The deterministic system is expected to imitate the stochastic process by using intrinsic chaoticity. The system repeatedly generates the quasi-attractors embedded in step 1 in synchronization with the pseudo-stochastic symbol transition, meaning that the design of the desired CI dynamics is completed.

In this study, we demonstrate that the trajectories of quasi-attractors and their transition rules can be designed using the three steps described above. In step 1, we show that the desired output dynamics can be designed with high operability by using the embedded internal dynamics reproducibly generated after the innate training. Next, in step 2, we demonstrate that various types of periodic symbolic sequences switching at a certain interval can be implemented simply by adjusting the parameters of a feedback loop attached to the system. Last, in step 3, we prepare several stochastic symbol transition rules governed by a finite state machine and show that the system can simulate these stochastic dynamics by making use of the system’s chaoticity. We also discuss the proposed method’s validity and adaptability through several numerical experiments.

## MATERIALS AND METHODS

### System architecture

In our method, we aimed to embed *M* types of quasi-attractors and the transition rules among them in an RNN. We prepared *M* discrete symbols *s* ∈ *S* (*S* ≔ {*s*_{1}, *s*_{2}, ⋯ *s _{M}*}). Each symbol corresponds to each quasi-attractor. We used an ESN as a high-dimensional chaotic system. As shown in Fig. 1A, we prepared an RNN composed of a nonchaotic input ESN (

*N*

^{in}nodes) working as an input transient generator as well as a chaotic ESN (

*N*

^{ch}nodes) yielding chaotic dynamics. The dynamics of input ESN

*x*^{in}(

*t*) ∈ ℝ

^{Nin}and chaotic ESN

*x*^{ch}(

*t*) ∈ ℝ

^{Nch}are given as the following differential equations

(1)

$$\mathrm{\tau}\frac{d{\mathit{x}}^{\mathit{ch}}}{\mathit{dt}}(t)=-{\mathit{x}}^{\text{ch}}(t)+\text{tanh}({g}^{\text{ch}}{J}^{\text{ch}}{\mathit{x}}^{\text{ch}}(t)+{J}^{\text{ic}}{\mathit{x}}^{\text{in}}(t))$$(2)where τ ∈ ℝ is a time constant, tanh is an element-wise hyperbolic tangent, *g*^{in} and *g*^{ch} ∈ ℝ are scaling parameters, *u*^{in}(*s*) ∈ ℝ^{Nin} is discrete input projected onto input ESN when symbol *s* is given, *J*^{in} ∈ ℝ^{Nin × Nin} and *J*^{ch} ∈ ℝ^{Nch × Nch} are connection matrices, and *J*^{ic} ∈ ℝ^{Nch × Nin} is a feed-forward connection matrix between input ESN and chaotic ESN. Each element of *J*^{in} is sampled from a normal distribution

. *J*^{ch} is a random sparse matrix with density *p* = 0.1 whose elements are also sampled from a normal distribution

. We used τ = 10.0, *g*^{in} = 0.9, and *g*^{ch} = 1.5 to make input ESN nonchaotic and chaotic ESN chaotic (*34*). In addition, to prevent chaotic ESN from becoming nonchaotic because of the bifurcation caused by the strong bias term, we tuned *J*^{ic} before hand to project transient dynamics converging to 0 onto the chaotic ESN when the same symbolic input continues to be given (see the Supplementary Materials for detailed information about the transient dynamics). In any case, the whole RNN dynamics ** x**(

*t*) ∈ ℝ

^{Nin + Nch}concatenating Eqs. 1 and 2 can be represented by the following single equation (⊙ represents an elementwise product)

(3)where ** x**,

**,**

*g**J*, and

**are defined by the following equations**

*u*(4)

$$\mathit{g}\u2254{[\underset{{N}^{\text{in}}}{\underset{\u23df}{{g}^{\text{in}},\cdots {g}^{\text{in}}}}\underset{{N}^{\text{ch}}}{\underset{\u23df}{{g}^{\text{ch}},\cdots {g}^{\text{ch}}}}]}^{T}$$(5)

$$J\u2254\left[\begin{array}{cc}{J}^{\text{in}}& 0\\ {J}^{\text{ic}}& {J}^{\text{ch}}\end{array}\right]$$(6)

$$\mathit{u}(s)\u2254[{\mathit{u}}^{\text{in}}(s);0]$$(7)

The output dynamics are calculated by the linear transformation of the internal dynamics ** x**(

*t*), that is, the linear readout

*w*

_{out}∈ ℝ

^{Nin + Nch}is trained to approximate the following target dynamics

*f*

_{out}(

*t*)

(8)

The symbolic dynamics *s*(*t*) itself, which is externally given in step 1, is lastly generated autonomously with a closed-loop system [Fig. 1B(2)]. In the feedback loop, the following classifier *f*_{max} : ℝ^{Nin + Nch} → *S* is attached

(9)where *w** _{s}* ∈ ℝ

^{(Nin + Nch) × M}represents the connection matrix whose elements are trained to autonomously emulate the designated symbolic dynamics

*s*(

*t*) [i.e.,

*s*(

*t*+ Δ

*t*) ≈

*f*

_{max}(

**(**

*x**t*)), where

*s*(

*t*+ Δ

*t*) is the symbolic input for the next time step and Δ

*t*is a time width for discrete temporal evolution]. To summarize, we designed the desired quasi-attractors, output dynamics, and symbolic dynamics by tuning the parameters of the RNN connections

*J*, the readout

*w*_{out}, and the classifier

*w**, respectively.*

_{s}### First order–reduced and controlled-error learning and innate training

We used two reservoir computing techniques called first order–reduced and controlled-error (FORCE) learning (*35*) and innate training (*36*). Both FORCE learning and innate training are methods that harness the chaoticity of the system. Below, we briefly describe the algorithms of both FORCE learning and innate training.

FORCE learning is a method that embeds designated dynamics in a system by harnessing the chaoticity of dynamical systems. Suppose the following ESN dynamics with a single feedback loop

$$\mathrm{\tau}\frac{d\mathit{x}}{\mathit{dt}}(t)=-\mathit{x}(t)+\text{tanh}(\mathit{gJ}\mathit{x}(t)+\mathit{u}z(t))$$(10)

$$z(t)={\mathit{w}}^{T}\mathit{x}(t)$$(11)where ** u** represents the linear feedback vector. Typically, the scaling parameter

*g*is set to be greater than 1 to make the whole system chaotic (

*34*). In FORCE learning, to embed the target dynamics

*f*(

*t*) in the system,

**is trained to optimize the following cost function**

*w**C*

_{FORCE}

(12)

Here, the bracket denotes the averaged value over several samples and trials. In particular, in the FORCE learning, ** w** is optimized online with a least-square error algorithm. It was reported from numerical experiments using ESN that better training performance was obtained when the initial RNN was in a chaotic regime (

*35*).

Innate training is also a scheme for harnessing chaotic dynamics and is accomplished by modifying the internal connection *J* using FORCE learning. The novel aspect of innate training is that the inner connection of ESN is trained in a semisupervised manner, that is, the connection matrix *J* of the ESN is modified to minimize the following cost function *C*_{innate} to reproduce the chaotic dynamics yielded by the initial chaotic RNN (*x*_{target}(*t*), innate trajectory)

(13)

Intriguingly, the innate trajectory is reproducibly generated for a certain period with the input while maintaining the chaoticity after the training. In other words, innate training is a method that allows a chaotic system to reproducibly yield the innate trajectory with complicated spatiotemporal patterns. In addition, innate training applies the FORCE learning method to the modifications of the internal connection, that is, the presynaptic connection of a node in the network is considered as the linear weight from the other nodes and trained by FORCE learning. In this study, we propose a method of designing CI by using both FORCE learning and innate training techniques.

### Recipe for designing CI

Our proposed method is a batch-learning scheme consisting of the following three-step process (Fig. 1C).

*Step 1. Designing quasi-attractor*. In step 1, the connection matrix *J*^{ch} of the chaotic ESN is adjusted by innate training to design the trajectories of quasi-attractors. First, the target trajectories

are recorded for *M* symbols under an initial connection matrix *J*^{init} and some initial states

, where

${\mathit{x}}_{\text{target}}^{s}(t)$ denotes chaotic dynamics when the symbol is switched to *s* at *t* = 0 ms (for simplification, the switching time is fixed to *t* = 0 ms in step 1; note that the symbol can be switched at any time). In step 1, *J*^{ch} is trained to optimize the following cost function *C*_{1−in}

(14)

Here, *x** ^{s}*(

*t*) represents the dynamics when the symbol is switched to

*s*at

*t*= 0 ms, and

*L*

_{innate}represents the time period of the target trajectory. We randomly choose half the network nodes (

*N*

^{ch}/2 nodes) and modify their presynaptic connections to reduce the redundancy of the training parameter. The selected elements in connection matrix

*J*

^{ch}are trained for 200 epochs for each

*s*. We lastly use

*J*

^{ch}recording the minimum

*C*

_{1 − in}(see the Supplementary Materials for the detailed algorithm used in step 1). After the innate training in step 1, the system is expected to reproduce the recorded innate trajectories

for *L*_{innate}.

Although there are no specific criteria for determining the initial states of the multiple innate trajectories, the large distances among

${\mathit{x}}_{\text{target}}^{s}(0)$are preferred since the temporal pattern of quasi-attractors is likely to differ, enhancing the separability. Moreover, the important trick of the innate training lies in its semisupervised scheme, that is, the training stability increases by guiding the offset states during training to the neighborhood of

${\mathit{x}}_{\text{target}}^{s}(0)$. Although a scheme for training the internal connections of the RNN has already been proposed in FORCE learning (*35*), the tuning of all connections to generate the same function is mostly unstable (*37*). Therefore, we randomly selected the offset state of the innate trajectory

on the phase space.

Similarly, *w*_{out} is trained to produce designated output dynamics *f ^{s}*(

*t*) corresponding to symbol

*s*. The following cost function

*C*

_{1−out}is optimized

(15)

High-dimensional nonlinear dynamical systems generally have high separability for input information, that is, it becomes easier for a linear model to solve nonlinear input-output function tasks by projecting input information into the system (*38*). In particular, the innate trajectories of chaotic systems are known to have such high expressive capability that various orbits can be designed simply by adjustment of the attached linear model (*36*). In this study, the tuned readout is also expected to stably reproduce the prepared trajectory by exploiting the high dimensionality and nonlinearity of the innate trajectories. Here, note that *L*_{innate} does not always match *L*_{out}, that is, *L*_{out} can be greater than *L*_{innate}. The training is accomplished by an offline algorithm Ridge regression based on the recorded internal dynamics *x** ^{s}*(

*t*).

*Step 2. Embedding autonomous transitions of symbol*. In step 2, we tune a feedback loop *f*_{max} to achieve the autonomous symbol transition. We especially prepare target periodic transition rules switching every *T* (ms). Suppose a target periodic symbolic time series *s*_{per}(*t*). First, the network dynamics ** x**(

*t*) of the open-loop setup [Fig. 1B(1)] is recorded with a symbolic dynamics

*s*

_{per}(

*t*) for

*T*

_{rec}≔ 500,000 ms. On the basis of the recorded dataset,

*f*

_{max}is tuned to output the next symbolic input

*s*

_{per}(

*t*+ Δ

*t*) from

**(**

*x**t*). The parameters

*w**of*

_{s}*f*

_{max}is trained to optimize the following cost function

*C*

_{2}

(16)

As the optimization algorithm, we use the limited-memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (*39*).

This optimization task is considered to be a type of timer task, a commonly used benchmark task to evaluate the temporal computational capacity, where the readout is trained to output a pulse-like wave with a certain delay after input is given. By projecting the input signal into a high-dimensional nonlinear dynamical system, the timer task can be achieved simply by adjusting the linear readout. In addition, the complex trajectory embedded by the innate training significantly increases the performance of the timer task compared with a nonchaotic random ESN (*36*). In our model, the tuned classifier *f*_{max} is expected to emulate the delayed symbolic switching by exploiting the embedded innate trajectory.

*Step 3. Embedding stochastic transitions of symbol*. In step 3, we implement a stochastic transition rule governed by a finite state machine by modifying a feedback loop *f*_{max}. As discussed in Introduction, the chaoticity of the system is expected to be used to emulate the stochastic process in the deterministic setup. We prepared the target stochastic time series *s*_{sto}(*t*) generated from a Markov model with certain switching periods. The process of the learning is same as that in step 2, that is, the pair of (** x**(

*t*),

*s*

_{sto}(

*t*)) recorded in the open-loop setup for 500,000 ms is used to train the

*f*

_{max}to emulate

*s*

_{sto}(

*t*). Here, we use the following cost function

*C*

_{3}in the training

(17)

As with the optimization of the cost function *C*_{2}, *C*_{3} is optimized with the limited-memory BFGS algorithm.

Note that the formulation of *C*_{3} is the same as that of *C*_{2}, that is, the properties of target dynamics to be embedded are not expressed in the cost function formulation. Rather, both optimization processes in steps 2 and 3 are data-driven, meaning that the properties of prepared target symbolic sequences determine whether the embedded dynamics are required to be chaotic or nonchaotic.

## RESULTS

In this section, we show the demonstration and analytic results of the numerical experiments for each step.

### Step 1. Designing quasi-attractor

As discussed in the previous section, the internal connection of the chaotic ESN *J*^{ch} is trained to reproducibly output the corresponding innate trajectories to the symbolic switching. Figure 2A demonstrates the change of the network dynamics of a 1500-node RNN (*N*^{in} = 500, *N*^{ch} = 1000) whose connection matrix is modified with innate training under the condition (*M*, *L*_{innate}) = (1,1000). The trajectory quickly spreads before *t* = *L*_{innate} in the pretrained system, whereas the target trajectory

(dotted line) is reproducibly yielded for 1000 ms (covered by the yellow rectangle) in the posttrained system. Moreover, intriguingly, the dispersion of the trajectories continues to be suppressed even after *t* = *L*_{innate}.

Next, Fig. 2B displays both the network dynamics and the output dynamics. The 1500-node RNN (*N*^{in} = 500, *N*^{ch} = 1000) trained under the condition (*M*, *L*_{innate}) = (3,1000) was used. At first, the symbolic input was absent, and then symbols were switched with random intervals from the middle. In addition, the two-dimensional readout was trained to output the Lissajous curve for symbol A, the “at” sign for symbol B, and the *xz* coordinates of the Lorenz attractor for symbol C for *L*_{out} = 1500 ms (the target trajectory for the “at” sign was made from a centerline of font data). It was observed that the desired spatiotemporal patterns were stably and reproducibly generated for a certain period in every trajectory with different initial states after the symbol transition (see movie S1). Note that the same linear model *w*_{out} was used in the demonstration, implying that the trajectory of each quasi-attractor has rich enough information to independently output the designated time-series patterns even with the single linear regressor. Our scheme for designing transient dynamics would be highly useful in the field of robotics because the process in step 1 is easily achieved by adjusting the partial elements of a high-dimensional chaotic system. For example, the system working in a real-world environment should immediately and adaptively switch its motion according to the change of environmental input like a system developed by Ijspeert *et al*. (*40*), which can be easily accomplished by our computationally cheap method. In this way, our method would work effectively in the context of robotics, where fast responsiveness and adaptability are required.

We also examined both the scalability and the validity of innate training in detail through several numerical experiments (Fig. 3). First, we examined the relationship between the number of input symbols *M* and the accuracy of innate training. To evaluate the performance of innate training, we used the normalized mean square error (NMSE) between the output and the innate trajectory

represented by the following formula

$$\text{NMSE}\u2254\frac{1}{M}{\Sigma}_{s\in S}\u27e8\frac{{\displaystyle {\Sigma}_{0}^{{L}_{\text{innate}}}}\parallel {\mathit{x}}^{s}(t)-{\mathit{x}}_{\text{target}}^{s}(t){\parallel}^{2}}{{\displaystyle {\Sigma}_{0}^{{L}_{\text{innate}}}}\parallel {\mathit{x}}_{\text{target}}^{s}(t){\parallel}^{2}}\u27e9$$(18)where the bracket represents the average over 10 trials for each symbol. We calculated the NMSE for 10 trials. Figure 3A shows the innate training performances with the different training conditions, suggesting that NMSEs are more likely to increase with a longer target trajectory and a larger number of symbols. This result implies that innate training has its limitation in the design of the quasi-attractors. We also examined the effect of network size on the capability to embed the quasi-attractors. We investigated the relationship between the number of nodes in the chaotic ESN *N _{c}* and the accuracy of innate training under the condition

*M*= 1 (Fig. 2B), suggesting that the NMSEs were less likely to increase with a larger network. To summarize, our analysis indicates that longer trajectories can be embedded in a larger network by innate training.

Next, we evaluated the effect of innate training on the capacity of the system’s information processing. We prepared a timer task and measured how long the inputted information was stored in the RNN. In the timer task, the pulse-like wave with a peak *t*_{peak} (ms) after the symbol transition was prepared as the target, and the performance was defined as the accuracy of the pulse-like wave reconstruction by a trained readout. Here, we defined the coefficient of determination value *R*^{2} between the output and the pulse-like wave as the timer task function *R*^{2}(*t*_{peak}). At the same time, we also calculated the integral value of the timer task function

and define it as the timer task capacity (see the Supplementary Materials for detailed information about the setup of the timer task). Figure 4C shows the timer task function with different innate training conditions, indicating that RNNs trained with the longer-length target trajectory *L*_{innate} perform better. It was also observed that the timer task capacity saturated around *L*_{innate} = 5000 ms in the 1000-node RNN, and the border of the saturation decreased in a smaller system (Fig. 3D). These results imply that the temporal information capacity of the system is improved by innate training with the longer target length *L*_{innate} but saturates at a certain value, which is determined by the system size.

Furthermore, we assessed the effect of innate training on the system’s chaoticity by measuring the Lyapunov exponents of the system. Since the transition among quasi-attractors is driven by the system’s chaoticity, it is necessary to keep the system chaotic. In this experiment, we measured the local Lyapunov exponent (LLE) to evaluate the degree of trajectory variation after the symbolic switching. We also measured the maximum Lyapunov exponent (MLE) without any inputs (** u**(

*t*) = 0) to estimate the global chaoticity of the system (see the Supplementary Materials for the detailed calculation algorithm of both the LLE and MLE). Figure 3E displays the LLE values of the systems with the different target trajectory length

*L*

_{innate}, suggesting that the trajectories unevenly expand after the symbol transition. In particular, it was observed from the LLE analysis that contracting regions existed (regions with negative LLEs corresponding to the lengths of the quasi-attractors) caused by the transient dynamics projected by the input ESN, and the degree of the expansion became gradual in the trained period

*t*∈ [0,

*L*

_{innate}). These results imply that innate training yields a locally contractive phase space structure, that is, a quasi-attractor. Moreover, positive MLE values were constantly obtained from the MLE analysis depicted in Fig. 3F, supporting the conjecture that the system chaoticity was maintained especially well with the larger RNNs even after the innate training. (Note that a sharp increase in MLE was observed with shorter

*L*

_{innate}, which is caused by the increase in the spectral radius of the connection matrix

*J*of the system. See the Supplementary Materials for detailed information of the analysis.)

### Step 2. Periodic symbol transition

In step 2, the system autonomously generates a symbolic sequence externally given in step 1. The additional feedback loop realizes the autonomous periodic switching of the symbols. We demonstrate that various types of periodic symbolic sequences switching at a fixed interval can be easily designed simply by tuning the parameter of the feedback loop *f*_{max}. Figure 4A demonstrates the embedding of the periodic symbolic sequence A-B-C (2000-ms interval and 6000-ms period) with a trained RNN ((*M*, *L*_{innate}) = (3,1000)). Figure 4B also exhibits the embedding of the periodic symbolic sequence A-B-C-D-E-F-G-H-I-J (500-ms interval and 5000-ms period), with the same RNN used as the demonstration in Fig. 4A (note that *f*_{max} was changed from one in Fig. 4A). In both demonstrations, the system succeeded not only in generating the desired symbol transition rules but also in stably outputting the designated output dynamics with high accuracy.

We also show that the system can solve tasks requiring higher-order memory in the same scheme. We prepared the two periodic symbolic sequences A-B-C-B and A-B-C-B-A and separately trained *f*_{max}. These two symbolic sequences are more difficult to embed because the system must change the output according to the previous output. In the symbol transition A-B-C-B, for example, the system must output the next symbol depending on the previous symbol when switching from B, though the total number of symbols is the same as in the task A-B-C. We used the same RNN and setup used in the Fig. 4A and only changed the parameters in *f*_{max} to realize the symbol transitions. Figure 4C displays the network dynamics and symbol transition of the two tasks, showing that the system successfully achieves both the periodic sequence A-B-C-B with an 8000-ms period and A-B-A-B-C with a 10,000-ms period. These results suggest that the trained RNN had the higher-order memory capacity, that is, the generated trajectories have sufficient separability to distinguish the contextual situation depending on the previous symbol sequence (see movie S2). In robotics, periodic motion control has often been implemented by an additional oscillator (e.g., a central pattern generator) to yield limit cycles (*23*, *40*–*42*). Our method in step 2 would be useful in designing limit cycles with longer periods and more complicated patterns. The analysis in fig. S2 shows that our method outperforms FORCE learning in the embedding of a long-term periodic attractor including multiple transitions in order (see the Supplementary Materials and fig. S2).

We also analyzed the effect of perturbation to investigate the stability of the embedded symbol transition. Figure 4D shows the output dynamics of both the original and perturbed trajectories, clarifying that the trajectory returned to the original one after the addition of the perturbation. We also calculated the MLE values of the system and obtained the value −1.89 × 10^{−4}, which was very close to zero. These analyses indicate that the trained feedback loop *f*_{max} made the system nonchaotic, that is, the generated internal dynamics was a limit cycle.

### Step 3. Stochastic symbol transition (CI)

In step 1, we constructed the trajectories of the quasi-attractors and the corresponding output dynamics. In step 2, we showed that periodic transitions among quasi-attractors can be freely designed by simply tuning the feedback loop *f*_{max}. In step 3, we realize a stochastic transition, that is, CI. As discussed above, the system is expected to use its chaoticity to emulate a stochastic transition in deterministic dynamical systems.

First, we demonstrate that stochastic transition can be freely designed by adjusting *f*_{max} (see Fig. 5A and movie S3). In this demonstration, we used the same RNN as in Fig. 4A. We prepared a symbol transition rule uniformly switching among symbols A, B, and C at 3000-ms intervals. Figure 5B shows the symbolic dynamics, network dynamics, and output dynamics, suggesting that the symbol transitions started to spread at around *t* = 10,000 ms and lastly settle down to completely different transition patterns. Nevertheless, the system continued to stably generate Lissajous curves. These demonstrations imply that the system constantly reproduced quasi-attractors embedded by innate training, while the quasi-stochastic transition was achieved by the global chaoticity. [Note that, although we demonstrated our approach by embedding typical stochastic processes (i.e., Markov processes) to illustrate the usability of our scheme, our method can also design a history-dependent stochastic rule that cannot be represented by a Markov model. See our demonstration represented in the Supplementary Materials and fig. S3.]

To analyze the flexibility of our method, we measured the stochastic transition matrix and the average symbolic intervals (Fig. 5B). We prepared two stochastic symbol transition rules as the targets: the transition rule governed by the uniform finite state machine (pattern 1) and the transition rule governed by the finite state machine with a limited transition (pattern 2). Note that we used the same trained RNN as in the demonstration in step 2 and embedded the transition rules simply by adjusting *f*_{max}. Figure 5B shows the results of the obtained trajectories, implying that the system successfully embedded patterns similar to the target rules, although there were some errors and variations in the transition probability and the switching time. The positive MLEs were obtained in both cases (+2.01 × 10^{−3} in pattern 1 and +1.71 × 10^{−3} in pattern 2), suggesting that the system was weakly chaotic as a whole. In addition, we analyzed the history dependence of the transition in detail, showing that transition probabilities did not differ so much according to the past symbol in both cases though the preferred routes were observed in output dynamics (see the Supplementary Materials and fig. S4). In this sense, it can be said that the system successfully expressed the random transitions in a macroscopic scale (i.e., a scale in symbol transitions) using the chaoticity.

Last, we analyzed both the structures of the obtained chaotic attractors and the symbolic dynamics in detail. Figure 6A shows the effect of small perturbations on the symbolic dynamics, implying that the patterns of symbolic dynamics varied after a certain period. To analyze the structural change of the terminal symbolic state, we measured the symbolic dynamics accompanied by the temporal development of the set of initial states on a plane constructed by the two selected dimensions (Fig. 6B), clarifying that a complex terminal symbolic structure emerges after a certain period (Fig. 6B). In particular, in the embedding of the pattern 1 rule, the entropy of the terminal symbolic pattern converges to a value close to the maximum entropy

${\text{log}}_{2}{3}^{9}\approx 14.26$(note that the entropy was measured on the basis of the probability distribution constructed by the frequency of 3 × 3 grid patterns). These results indicate that the symbol transition markedly changed even with a small perturbation and was unpredictable after a certain period, that is, the prediction of symbolic dynamics required the complete observation of the initial state value and calculation of the temporal development with infinite precision.

## DISCUSSION

In this study, we proposed a method of designing CI based on reservoir computing techniques. We also showed that the various types of output dynamics and symbol transition rules could be designed with high operability simply by adjusting the partial parameters of a chaotic system with our three-step recipe. In this section, we first discuss the scalability of our method and the mechanism of how CIs are successfully embedded by reviewing several numerical analyses that verify the validity of our methods. Next, we discuss the effectiveness and significance of our method from multiple viewpoints.

### Scalability and validity

First, the results of the innate training performances displayed in Fig. 3 (A and B) indicate that the number of RNN nodes constrains the total length of the quasi-attractors that can be embedded in the system by the innate training. However, the LLE analyses in Fig. 3E show that the system has the expanded region of the negative LLE even when the NMSE between the innate trajectory and the embedded trajectory becomes large (e.g., *L*_{innate} = 5000 ms). These results imply that, even when the innate trajectories are not successfully embedded in the system, the system stably yields high-dimensional trajectories with complicated spatiotemporal patterns for each symbol transition over *L*_{innate}, which is caused by the weakening of the system chaoticity. The same RNN trained under the condition (*M*, *L*_{innate}) = (3,1000) was repeatedly used in our series of demonstrations, the desired output dynamics (e.g., the Lissajous curves) being constantly generated for *L*_{out} = 1500-ms periods after the symbolic shift (Figs. 2, 4, and 5). We also demonstrated that the system can autonomously generate a symbol transition rule with an interval greater than *L*_{innate} (Figs. 4 and 5), suggesting that the system exploited high-dimensional reproducible trajectories longer than *L*_{innate}.

Moreover, it is assumed that the length of the quasi-attractors constrains the target stochastic transition rules that can be embedded. The system failed to imitate the stochastic transition, and the transition became periodic when the target transition had a shorter switching interval, whereas the training of *f*_{max} became unstable when it had a longer switching interval. These results suggest that the following two mechanisms should be required in the design of CI in our method: (i) The differences among the trajectories are sufficiently enlarged through the temporal development to realize the stochastic symbol transition and (ii) a similar spatiotemporal pattern should be reproducibly yielded until the switching moment to precisely discriminate the switching timing. These two mechanisms are contradictory, of course, and the desired CI can likely be embedded when both conditions are moderately satisfied.

Next, we discuss the validity of CI designed with our method through comparison with previous CI studies. We demonstrated that the embedded transition yields transition probability similar to the desired one by using the chaoticity of the network. Several CI works, however, pointed out that transition of CI is history dependent, and therefore, CI shows some preference in the transition determined by the history. Kaneko and Tsuda (*43*), for example, pointed out that the transition of CI has a specific order that is distinguished from a simple random hopping. Itoh and Kimoto (*44*) also reported that the transition presents a preferred route governed by the global phase space dynamics. Our experiment indicated that embedded CI also shows preferred trajectories in output dynamics according to the previous symbol before the switching (see the Supplementary Materials and fig. S4A). Therefore, the symbolic dynamics implemented with our methods should be history dependent.

Moreover, it should be noted that our proposed method can successfully emulate the randomness of the symbol transition despite the existence of the preferred trajectory. Our analysis revealed that the transition probabilities did not differ so much according to the difference in the previous symbol in both demonstrations (i.e., stochastic rules 1 and 2), implying that prediction of the next transition is still not easy if done only by referring to the macroscopic symbolic history (see the Supplementary Materials and fig. S4B). In this way, macroscopic symbol transitions still appear to be random, despite the trajectory preference. Thus, it can be said that our proposed method emulates the transition probability and its randomness on a high-dimensional nonlinear dynamical system by using the chaoticity of the system.

We show that our method can realize CI characterized by the random transition of a finite number of quasi-attractors as shown in (*19*, *21*). Some classes of CI, however, are difficult to design even with our method. There exists a CI whose transition frequency has a long-tailed distribution (*9*). In addition, CI can have infinitely many quasi-attractors whose dynamics should be described with an infinite state machine (*7*). These CI properties are hard to design with our method since an infinite number of auxiliary symbols and an infinite length of corresponding trajectories should be prepared, which should be solved in the future.

Last, we discuss the speed of recovery from the transient state, that is, the stability of the embedded quasi-attractor. Ahmadi and Tani (*45*), for example, demonstrated variational RNN called predictive coding-inspired variational RNN (PV-RNN) to imitate learning of stochastic transition between presented primitives, in which the stability of the attractor is determined by the coefficient of the complexity of variational evidence lower bound in Bayesian inference expressed in the cost function. Therefore, in their approach, the formulation of the cost function can explain the mechanism of how the recovery speed is determined. On this point, the stability of the quasi-attractor in our approach is considered to be determined by the spectral radius, a parameter of ESN, and the connection strength between the input ESN and the chaotic ESN. As discussed above, a transition from nonchaotic to chaotic regime occurs when the spectral radius of ESN exceeds 1.0. In addition, the chaoticity of the system is strengthened by the increased spectral radius. Hence, it is speculated that the recovery speed becomes lower as the spectral radius becomes larger since the quasi-attractor becomes more unstable. Besides, the network state is more likely to synchronize when the signals sent to the chaotic ESN are larger. Thus, the recovery speed would be also controlled by the connection strength governing the amplitude between the input ESN and the chaotic ESN. To summarize, the ESN spectral radius and the connection strength are considered to determine the recovery speed in our architecture.

### Effectiveness and significance

First, the high operability of our proposed model would be helpful to understand the underlying mechanism of brain’s information processing from a certain perspective. Previous studies have pointed out that an enormous number of nonlinear units and their interaction essentially constitute an animal’s nervous system and yields highly complicated activities characterized by nonlinear phenomena such as chaos and CI. It has been reported that chaotic behavior appears in a wide range of brain activities from the cell level (e.g., action potential) to the global measurement level (e.g., electroencephalogram) (*46*). In addition, it has been pointed out that the collective neural activities not only encode the external information but also transform it according to the history of activities (*47*, *48*), suggesting that the animal brain realizes its information processing through the high-dimensional activities. We found that high-dimensional chaos has enough rich expressive capability to design CI, implying that the high-dimensional chaotic brain activities potentially have the capability to realize various functional hierarchies. In this sense, our model would provide a clue to understand the mechanism of how high-dimensional chaos that contributes to the information processing in animal brains.

The designing method for the CI exhibited in this study offers fundamentally different benefits compared with the previous methods, which simply exploit chaotic dynamics. The dynamics of CI presents an interesting property and suggest that local coherence and global chaoticity coexist. This property of CI would effectively work especially in designing cognitive models where both autonomy and spontaneity are required. For example, a designer can implement the motion primitives through the quasi-attractor while maintaining the autonomy of the robot through the global unpredictability. Furthermore, our algorithm can design the probability distribution, that is, the global tendency of behavior, as shown in Fig. 5. This coexistence is, however, difficult to express when using only the conventional chaotic attractors. In addition, several studies pointed out that animals’ cognitive functions, such as memory recall and association, would be realized through the transition phenomenon among stereotypical activities (*13*, *14*), suggesting that CI might play an important role in animal cognition. In this sense, our method would be used for implementing cognitive models in a high-dimensional dynamical system.

Although several studies have used chaotic dynamical systems to embed desired trajectories (*35*, *36*), these conventional methods are incapable of combining multiple predetermined transient dynamics like our method. It might also be possible to merge multiple transients into one big attractor and embed it simultaneously by a learning scheme such as FORCE learning. However, the additional experiment shows that FORCE learning is more unstable than our methods when it comes to embedding long-term periodic trajectories consisting of multiple transients (see the Supplementary Materials and fig. S2). In this sense, high operability in our method is unavailable in the conventional methods.

Unlike previous methods that construct desired trajectories by tuning the entire dynamics by backpropagation algorithm, our method is accomplished by adjusting the reduced number of parameters and using the intrinsic high-dimensional chaos, which alleviates the biological implausibility and computational complexity of backpropagation algorithms. For example, recent physiological studies on the motor cortex (*49*, *50*) suggest that a large variety of behaviors can be instantaneously generated by the partial plasticity of the nervous system, supporting the biological plausibility of our learning scheme. In addition, our learning scheme is computationally cheaper than backpropagation since adjustments of entire neural circuits are not necessary. These properties would be especially helpful in the context of bioinspired robotics, where fast responsiveness and real-time processing are required.

Another advantage of our method is that it does not require the explicit structure of dynamical systems. For example, in the method proposed by Namikawa and Tani (*20*–*22*), the controller needs a fixed hierarchical structure and modularity. Therefore, the trained controller was specialized in implementing a specific behavior, making it difficult to divert it for any other purpose. In addition, it may be possible to design CI-like dynamics in an architecture where the symbolic sequence and the corresponding trajectory are separately generated, which requires an external mechanism to hold the symbol and wait until the generation of lower-order trajectory finishes. Therefore, the separation of the symbolic sequence model and trajectory encoder implicitly uses the hierarchical structure and cannot be realized by a high-dimensional chaos alone. In contrast, we proposed a method of designing CI with a setup consisting of a single chaotic ESN, auxiliary symbols, and an interface between them (input ESN) with high scalability. Moreover, the modifications of internal connections in the chaotic ESN can be realized by adding multiple linear feedback loops and training them with the FORCE learning since the presynaptic connection in the chaotic ESN can be regarded as a linear connection. Thus, our method allows us to design the various trajectories and their transition rules in a consistent high-dimensional chaotic system, thereby greatly expanding the scope of application of high-dimensional chaotic dynamical systems. Neuromorphic devices based on physical reservoir computing frameworks would be an excellent candidate for implementing our scheme (*31*). Sprintronics devices, for example, have recently been shown to exhibit chaotic dynamics (*51*, *52*) and are actively exploited as physical reservoirs (*53*–*55*). We expect that this framework would provide one of the promising application scenarios for real-world implementations of our scheme.

Our method can use a priori knowledge through the introduction of auxiliary symbols. Several neurorobotics frameworks have been proposed so far in which symbolic dynamics are self-organized on the network by end-to-end learning (*26*, *56*). Although these methods are convenient since they do not require explicit a priori knowledge, they cannot actively use the prior symbolic structure, and thus, symbolic structure only appears after the training. In contrast, we showed that the learning performance is greatly improved by auxiliary symbols (see the Supplementary Materials and fig. S3). In that sense, our method has an advantage over the conventional end-to-end scheme.

Our method is also scalable to autonomous symbol generation required in more advanced functionality. For example, in our method, *M* kinds of auxiliary symbols are given as a priori knowledge. However, in a highly autonomous system, such as humans, symbols are dynamically generated and destroyed because of developmental processes. As demonstrated by Kuniyoshi and Sangawa (*17*), these self-organizing symbolic dynamics can be realized by providing an additional automatic labeling mechanism in the system. In addition, unsupervised algorithms for extracting discrete symbols from the dynamics like the one introduced in (*57*–*59*) can be incorporated into our system. In other words, it is possible to spontaneously generate symbols by embedding an unsupervised learning algorithm in the system; this is a subject for future work.

Last, the dynamic phenomena obtained by our method are significant from the viewpoint of high-dimensional dynamical systems. As shown in Fig. 6, we demonstrated that small differences in the initial network state were expanded by the chaoticity of the system, which eventually led to drastic change in both the global symbol transition pattern *s*(*t*) and the local dynamics ** x**(

*t*). Such tight interaction between microlayer and macrolayer is a phenomenon unique to deterministic dynamical systems; that is, it cannot occur, in principle, in a system where the higher-order mechanism is completely separated from the lower-order one (e.g., independent random variables). In addition, the global characteristics of dynamical systems are often analyzed by the mean-field theory. However, the analysis by the mean-field approximation cannot capture the contribution of microscopic dynamics to the macroscopic change. Thus, our CI design method has a meaningful role in shedding light on the interaction between micro- and macrodynamics in deterministic chaotic dynamical systems.

**Acknowledgments: **This work was based on results obtained from a project commissioned by the New Energy and Industrial Technology Development Organization (NEDO). **Funding:** K.I. was supported by JSPS KAKENHI (grant number JP20J12815). K.N. was supported by JSPS KAKENHI (grant number JP18H05472) and by MEXT Quantum Leap Flagship Program (MEXT Q-LEAP) (grant number JPMXS0120319794). This work was supported by NEDO [serial numbers 15101156-0 (dated 24 June 2016) and 18101806-0 (dated 5 September 2018)] and Chair for Frontier AI Education, School of Information Science and Technology and Next Generation AI Research Center [serial number not applicable (dated 1 June 2016)]. **Author contributions:** K.I. and K.N. conceived the idea and designed the experiments. K.I. carried out the experiments and created the demonstration. K.I. and K.N. wrote the manuscript. All authors discussed and commented on the manuscript. K.N. and Y.K. directed the project. **Competing interests:** The authors declare that they have no competing interests. **Data and materials availability:** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.