We consider learning of K sequences, each of which contains M patterns, with K input patterns. We denote the μ-th targeted pattern in the α-th sequence as ξμα, and the corresponding input as η^α for μ = 1, 2, ⋯ , M over the inputs α = 1, ⋯ , K. Figure 1A illustrates the case with K = 2 and M = 3: In this case, a given sequence ξ1α,ξ2α,ξ3α (α = 1, 2) should be generated upon a given corresponding input η^α. Generally, a pattern to be generated next is determined not only by the current pattern but also by earlier patterns. Thus, a network has to retain the history of previous patterns to generate a sequence correctly.

To achieve this, we built a two-population model with different timescales, one with N fast neurons and one with N slow neurons, denoted as X and Y, respectively. X receives an external input, and Y receives the output from X and provides input to X, as shown in Figure 1A. The neural activities x_i in X and y_i in Y evolve according to the following equation:

where ui=∑j≠iNJijXxj; ri=∑jNJijXYtanh(yj). JijX is a recurrent connection from the j-th to the i-th neuron in X, and JijXY is a connection from the j-th neuron in Y to the i-th neuron in X. J^X is a fully connected network without self connections. It is modified during the learning process as described in the following subsection “Learning model” and initialized with the binary values P[JijX=±(N-1)-1/2]=1/2. The diagonal entries of J^X are kept at zero during the entire learning process. J^XY, in contrast, is a non-plastic sparse network; P(JijXY=±cN-1/2)=ρ and P(JijXY=0)=1-2ρ. X is required to generate the pattern ξμα in the presence of η^α, i.e., an attractor that matches ξμα is generated under η^α. The i-th element of a targeted pattern denoted as (ξμα)i, is assigned to the i-th neuron in X, and randomly sampled according to the probability P[(ξμα)i=±1]=1/2. The input (ηα)i is injected to the i-th neuron in X, randomly sampled according to P[(ηα)i=±1]=1/2. ξ and η are the same dimensional vectors as the fast dynamics, i.e., N-dimensional vectors. We set N = 100, β_x = 2, β_y = 20, τ_x = 1, τ_y = 100, ρ = 0.05, and c = 7. The dependence of the performance on these parameters are shown in the Supplementary Materials, while the details for the parameter setting are described.

In Equation (3), we implement the non-linear activation at the apical dendrites (Larkum et al., 2009). We assumed that the input from Y to X innervated on the apical dendrites of neurons in X, which is consistent with observations that the feedback inputs from the higher cortical areas innervate on the apical dendrites of layer 5 pyramidal neurons (Larkum, 2013), whereas the recurrent input from X to X was assumed to innervate the proximal dendrites. The synaptic inputs to the apical dendrites are integrated and evoke the calcium spike when the integrated input exceeds the threshold of spikes (refer to the detail of this type of spike in Larkum et al., 2009). To reproduce this information processing at the apical dendrites, we used two nonlinear filters by the hyperbolic tangent function for the input from Y to X. First, by adopting tanh(y_j), the activity of the j-th neuron in Y is amplified in a nonlinear way at a synapse onto the neuron x_i. Second, tanh(r_i) represents calcium spike at the branching point of the tuft dendrite. Even if these hyperbolic tangent functions are not leveraged in Equation (3), the behavior of the model is not changed qualitatively, although the performance of the model is reduced.

Only J^X changes to generate the target according to the following equation:

where τ_syn is the learning speed (set to 100). This learning rule comprises a combination of a Hebbian term between the target and the presynaptic neuron, and an anti-Hebbian term between the pre- and post-synaptic neurons with a decay term uiJijX for normalization ¹. This form satisfies locality across connections and is biologically plausible (Kurikawa et al., 2020). We previously applied this learning rule to a single network of X and demonstrated that the network learns K maps between inputs and targets, i.e., M = 1 (Kurikawa and Kaneko, 2013, 2016; Kurikawa et al., 2020). However, in that case, generating a sequence (M≥2) was not possible. In the present study, there are two inputs for X, one from input η and one from Y that stores previous information. Thus, the network can generate a pattern depending not only on the present input pattern, but also on the previous patterns.

In our model, the patterns in the sequence are learned sequentially. A learning step of a single pattern is accomplished when the neural dynamics satisfy the following two criteria: x sufficiently approaches the target pattern, i.e., mμx≡Σixi(ξμ1)i/N>0.85, and y is sufficiently close to x, i.e., Σ_ix_iy_i/N>0.5. After the completion of one learning step, a new pattern ξ21 is presented instead of ξ11 with a perturbation of fast variables x_i, by multiplying a random number uniformly sampled from zero to one. We execute these steps sequentially from μ = 1 to M to learn a sequence once, denoted as one epoch of the learning. Before finishing the learning process, this procedure is repeated 20 times (i.e., 20 epochs). The second criterion for terminating the learning step is introduced for memorizing the sequences, especially the history-dependent sequences. Further, the value 0.5 of this criterion must take an intermediate value. If this criterion is not adopted or this criterion value is small, the target pattern is switched as soon as x is close to the target during the learning process. At this time, y is far from x because y is much slower than x. In this case, y cannot store any information about x. On the other side, when the value is close to unity, y matches x and y can store only the present x. In both cases, y cannot store the history of x.

To elucidate how such a recall is possible, we analyzed the phase space of x with y quenched. In other words, y is regarded as bifurcation parameters for the fast dynamics. Specifically, we focused on the neural dynamics for 200 ≤ t ≤ 500, as shown in Figure 2A. In this period, the fast dynamics show transitions from ξ11 to ξ21 at t = 290, from ξ21 to ξ31 at t = 375, and from ξ31 to ξ11 at t = 220, 460. We sampled the slow variables every five units of time from t = 200 to 500, y_{t = 200}, y_{t = 205}, ⋯ , y_{t = 500}, along the trajectory, and analyzed the dynamics of x with the slow variables quenched at each sampled y_{t = 200, 205, ⋯ , 500}. Figure 2B shows the bifurcation diagram of x against the change in y, and Figure 2C shows the trajectories of x for specific y.

We now consider the neural dynamics for y_{t = 225}, just after the transition from ξ31 to ξ11 [Figure 2C (i)]. For this y, a single fixed point corresponding to the present pattern (ξ11) exists, leading to its stability against noise. As y is changed, the basin of ξ11 shrinks, while a fixed point corresponding to the next target ξ21 appears, and its basin expands ², as shown in Figure 2C (ii). At y_{t = 290}, the fixed point ξ11 becomes unstable. Thus, the neural state x at ξ11 goes out of there, and falls on ξ21, i.e., a transition occurs.

With a further shift of y, y_{t = 295, 300, ⋯}, a regime of coexistence of ξ21 and ξ31 with large basins appears [Figure 2C (iii)]. The basin of the attractor ξ21 shrinks and vanishes [Figure 2C (iv)], and the transition from ξ21 to ξ31 occurs at t = 375. The next transition from ξ31 to ξ11 occurs in the same manner at t = 460. These processes provide the mechanism for robust sequential recall: fixed points x of the current and successive targets coexist, and then, the current target becomes unstable when the slow variables change.

To examine the robustness of the recall, we explored trajectories from different initial conditions with Gaussian white noise with strength s (refer to Supplementary Material for details). All of these trajectories converge correctly to a target sequence after some transient period for weak noise. By increasing the noise strength, the recall performance of noisy dynamics is made equal to that of the noiseless dynamics up to noise strength s = 0.3. For stronger noise, the duration of residence at the target is decreased, because the neural state of x is kicked out of the target earlier than in the noiseless case. Even upon applying a strong and instantaneous perturbation to both x and y, the trajectory recovers the correct sequence. The sequence is represented as a limit cycle containing x and y and, thus, is recalled robustly.

Next, we test if our model flexibly infers new sequences based on the previously learned sequence. To this end, we consider the following task (Refer to Materials and methods for details). First, a network learns a sequence (ξ21,ξ31,ξ41)=(A,B,C) in response to the associated input sequences (η21,η31,η41)=(a,b,b). In addition, we provide η11=s associated with ξ11=S preceding the sequence as a fixation cue and the response to it (e. g., the subject's gaze to the fixation point), respectively. After the learning is completed, the network should generate the sequence (S, A, B, C) in response to the input sequence (s, a, b, b). Then, the network learns a new sequence (ξ12,ξ22,ξ32)=(S,A′,B), which is associated with (η12,η22,η32)=(s,a′,b). In this study, we intend to examine if the inference of C is achieved selectively from (S, A′, B) in the presence of b. For it, we need to prevent the tight and trivial association between input b and the sub-sequence (B, C). For this purpose, the network is also postulated to learn the association between η13=b and ξ13=D as a distractor. We explore if the network generates the sequence (S, A′, B, C) in response to the input sequence (s, a′, b, b) after learning the association between (S, A′, B) and (s, a′, b).

During the learning of the first sequence, the overlaps with all of the targets reach more than 0.9 after 20 epochs of learning, as shown in Supplementary Figure S4A. Actually, Figure 4A shows that the first sequence is successfully generated. Next, the network learns the new sequential patterns (S, A′, B). If the network infers (S, A′, B, C) by using the already learned sub-sequence (B, C), it generates the sequence (S, A′, B, C) after learning only (S, A′, B) without (S, A′, B, C). As expected, the average overlaps with all of the targets in the second sequence (not only A′,B,but also C) are increased through learning (Figure 4B). After the fifth epoch of learning, the overlap with C declines, whereas those with A′ and B continuously increase. As an example, we plot the recall dynamics after learning A′ and B in Figure 4C. A′ evokes B and C, although the overlap with the first target A′ is not large. Simultaneously, the network generates the first sequence (Supplementary Figure S4C). Note that the sub-sequence (B, C) is generated selectively by the input either (s, a, b, b) or (s, a′, b, b). In contrast, when a random input v is given instead of a or a′ in an input sequence, the sub-sequence (B, C)is not evoked, as shown in Supplementary Figure S4B. To examine the difference among the recall dynamics in response to (s, a, b), (s, a′, b), and (s, v, b), the success rate of generating the sub-sequence (refer to its definition in Materials and methods) is analyzed statistically in Figure 4D. We found that the input sequences (s, a, b), (s, a′, b) evoke the sub-sequence (B, C) with a significantly higher rate than the input sequence (s, v, b). Thus, our model is able to infer a new sequence based on the previously learned sequence.

We examined if the proposed model learns the history-dependent sequence (M = 6), in which the same patterns exist in a sequence such as (ξ11,ξ21,⋯,ξ61)=(A,B,C,D,B,E). The patterns succeeding B are C or E, depending on whether the previous pattern is A or D. Then, the neural dynamics have to retain the information of the target A or D, to recall the target C or E correctly. Our model succeeded in recalling this sequence, as shown in Figure 5A. Just before the target C and E are recalled, there is no clear difference in the values of fast variables x, as indicated by the circles in Figure 5B. However, the values of slow variables y are different, depending on the previous targets shown in Figure 5C, which stabilize different patterns of x. Furthermore, we demonstrate that our model succeeded in recalling more complex sequences (M = 8) such as (ξ11,ξ21,⋯,ξ81)=(A,B,C,D,E,B,C,F), as shown in Figures 5D–F. In this case, the neural dynamics have to keep three previous targets in memory to recall the target D or F after B and C. As expected, generating the sequence with M = 8 is a harder task than that with M = 6. However, some networks still can generate the sequence.

To understand the performance comprehensively, we measured the success rate in generating these sequences for different lengths. For the sequence with the length of M, the network is required to retain the information about M/2−1 preceding patterns. We examined the success rate for M = 6, 8, and 10 over 50 networks realizations in the same manner as that in the inference task. The success rate is reduced for the longer M and nearly zero for M = 10. This result indicates that our model can store three preceding patterns at a maximum, but is difficult to memorize four preceding patterns.

We, next, explored whether the model can memorize another type of the history-dependent sequence such as (ξ11,ξ21,⋯,ξ61)=(A,B,A,C,A,D), as shown in Figure 6. The network is required to discriminate three neural states in the slow dynamics just before x approaches B, C, and D, as shown by circles in Figure 6C. When the network discriminates these states successfully, it generates the sequence adequately, as shown in Figures 6A–C. We measured the success rate in generating these sequences for different numbers of the states to be discriminated (namely, M/2 states for an M-pattern sequence) in Figure 6D. For the shortest case, the success rate takes less than 0.4.

As a final example of the complex sequences, we explored learning two history-dependent sequences (Figure 7), namely, (ξ11,ξ21,ξ31)=(A, B, C) upon η¹, and (ξ12,ξ22,ξ32)=(C, B, A) upon η². In these sequences, the flow A→B→C on the state space under η¹ should be reversed under η². The learned network succeeds in generating these sequences. Although orbits of x under inputs almost overlap in the 2-dimensional space, those of y does not. This difference in y, in addition to inputs, allows the orbits of x in the reverse order of patterns. Generally, y is different depending on the history of the previous patterns and inputs even when x is the same. Different y stabilizes different fixed point of x, to generate the history-dependent sequence.

Generating these sequences is rather hard. Actually, the success rate is 0.14 for three-pattern sequences [(A, B, C) under one input, (C, B, A) under the other input] and 0.08 for four-pattern sequences. The success rate is low even for the shortest sequences. The difficulty of this task could be attributed to the request that networks have to memorize the bidirectional transitions between the target patterns (namely, the transition from A to B and its reverse transition) dependent on the external input.

Recall performance is highly dependent on the relation between τ_x, τ_y, and τ_syn. To investigate the dependence of the performance on the timescales, we trained fifty realizations of networks for various values of timescales and calculated the success rate of training as a function of τ_syn for different τ_y by fixing τ_x at 1, as are plotted after rescaling τ_syn by τ_y in Figure 8A. The ratios yield a common curve that shows an optimal value ~1 at τ_syn, approximately equal to τ_y. As an exceptional case, the success rate for τ_y = 10 yields a lower value for the optimal τ_syn because τ_y is too close to τ_x to store the information about x. The balance between τ_syn and τ_y is important to regulate the success rate when they are sufficiently smaller than τ_x.

To unveil the significance of the timescale balance, we, first, present how the recall is failed for τ_y>>τ_syn, (τ_y = 100, τ_syn = 10 in Figure 8B). Some of the targets are recalled sequentially in the wrong order, whereas other targets do not appear in the recall process. To uncover the underlying mechanism of the failed recall, we analyze the neural dynamics of fast variables with slow variables quenched in a manner similar to that shown in Figure 2 (refer to Supplementary Figure S5A). Here, all the targets are stable for certain y, although ξ21 does not appear in the recall process. We also found that fixed points corresponding to ξ11 and ξ21 do not coexist for any y: the fixed point corresponding to ξ31 has a large basin across all y. This leads to a transition from ξ11 to ξ31 by skipping ξ21, and thus, the recall is failed.

Interestingly, how a recall is failed for τ_y < < τ_syn are distinct from that for τ_y>>τ_syn. For τ_y = 100, τ_syn = 1, 000, only the most recently learned target is stable for almost all y, and thus, only this target is recalled, as shown in Figure 8C. We sampled the slow variables from the last learning step of the sequence (Supplementary Figure S5B), and analyzed the bifurcation of the fast variables against change in slow variables, in the same way as above. In this study, only the latest target (here, ξ31) is a fixed point, whereas the other targets are not. Thus, transitions between targets are missed, except the transition to the latest target.

Why does the network fail in generating the sequences when τ_y and τ_syn are not matched? The reason is as follows: the learning time for the non-optimal timescales takes a longer time than that for the optimal ones. In this study, we consider the time that is normalized by τ_syn because it characterizes the timescale regulating the synaptic plasticity and, consequently, the neural dynamics. For τ_y>>τ_syn (Supplementary Figure S5C), the trajectories of neural activities are similar to those for τ_y~τ_syn, but it takes a longer time for y to approach x after x converges to the target. As the approach to x is so slow, the learning process stabilizes the present target during the approach by modifying the connectivity, so that the target is too stable over a wide range of y. On the other hand, for τ_y < < τ_syn (Supplementary Figure S5E), the neural activities of x and y wander before x converges to the target, resulting in disruption of information about the previous targets stored in y. Thus, the networks for both cases of τ_y>>τ_syn and τ_y < < τ_syn fail to generate the sequence. These results indicate that the relative timescale τ_y to τ_syn changes the bifurcation of the fast dynamics and the memory capacity.