Network architecture, short-term synaptic plasticity (STSP) mechanism and sequence working memory task paradigm.

(A) Schematic of the recurrent neural network (RNN) architecture. Sensory inputs from 30 tuned neurons are projected into a recurrent layer consisting of 80 excitatory (red) and 20 inhibitory (blue) neurons. The recurrent layer then projects to output units representing different decision states (fixation, match, non-match). (B) Microscopic dynamics of STSP. The evolution of the available neurotransmitter fraction (x, red line) and the utilization fraction (u, blue line) in response to presynaptic spikes (green) is illustrated for both short-term facilitation (left) and short-term depression (right). The transient synaptic efficacy (u ⋅ x, black line) dynamically determines the amplitude of the resulting postsynaptic current (yellow). (C) Timeline of the sequence working memory task. Following a 500 ms fixation period, a sequence of six oriented bar stimuli is presented during the sample period, followed by a 1000 ms delay. In the rhythmic condition, sample items are presented with constant inter-onset intervals (equal IOIs). In the arrhythmic condition, the sequence features temporal jitter with variable IOIs. (when stimulus duration (SD) = 100 ms, IOI ∈ {250, 300, 350, 400, 450, 500} ms; when SD = 120 ms, IOI ∈ {300, 350, 400, 450, 500} ms; when SD = 140 ms, IOI ∈ {350, 400, 450, 500} ms) The test sequence is semi-rhythmic in both conditions, prompting a final match/non-match behavioral response.

Rhythmic temporal structure improves sequence memory and organizes population-level representations.

(A) Behavioral accuracy in the working memory task under rhythmic and arrhythmic conditions (match and non-match trials). Rhythmic sequences showed significantly higher task accuracy (two-sided paired t-test, n = 100 independent networks, p < 0.001). (B) The behavioral benefit of rhythmicity was positively correlated with the base IOI (Pearson correlation, two-sided, n = 67 independent networks with SD = 100 ms, p < 0.05). (C) Distribution of the trial regularity index (R) under rhythmic and arrhythmic conditions. Larger R values indicate a higher degree of rhythmicity (see Methods: Task paradigm). (D-E) Population-level demixed principal component analysis (dPCA) of network representations under rhythmic (D1-D3) and arrhythmic (E1-E3) conditions. Rhythmic input were associated with more separable low-dimensional population trajectories across temporal, ordinal, and stimulus-related components. Unless otherwise stated, data are shown as mean ± s.e.m. Statistical tests are two-sided where applicable. Significance levels are denoted as *p < 0.05, **p < 0.01, and ***p < 0.001.

Temporal regularity selectively organizes sample-period oscillatory dynamics.

(A1) Frequency-resolved response strength T (f) of population activity during the sample epoch, plotted as a function of relative frequency (f /f0). The solid line denotes the across-trial mean response (n = 4000 trials). Blue horizontal bars indicate frequency clusters showing significant stimulus-evoked responses, assessed using a one-sample sign-flip cluster-permutation test against zero (t statistic, B = 1000 permutations). Significant stimulus-responsive clusters were observed at 0.71–1.58 f /f0 (1.78–3.96 Hz; Cohen’s dz = 1.74, cluster-level p < 0.001) and 2.05–2.85 f /f0 (5.14–7.14 Hz; Cohen’s dz = 0.97, cluster-level p < 0.001). Purple bars indicate the portions of stimulus-responsive frequencies that overlapped with significant regularity-sensitive clusters. (A2) Regularity-dependent modulation of spectral responses across frequency. The curve denotes the frequency-wise effect size ω2 from a one-way ANOVA across regularity bins. Red horizontal bars indicate significant regularity-sensitive clusters, assessed using label-permutation cluster testing on F statistics (B = 1000 permutations). Significant regularity-sensitive clusters were observed at 0.53–0.68 f /f0 (1.33–1.69 Hz; peak/mean ω2 = 0.53/0.38, cluster-level p < 0.001) and 0.75–3.50 f /f0 (1.87–8.75 Hz; peak/mean ω2 = 0.80/0.41, cluster-level p < 0.001). Purple bars indicate overlap between stimulus-responsive and regularity-sensitive frequency ranges. (B) Regularity-dependent response strength around the retained stimulation-frequency band. Each row corresponds to one regularity bin, and color indicates the bin-averaged response strength T (f) in arbitrary units. Directional modulation was quantified by Spearman’s rank correlation between regularity-bin medians and bin-averaged responses (ρ = 0.65, permutation test, p < 0.001), indicating stronger f0-band responses with increasing temporal regularity. (C) Preferred encoding phase for each sample stimulus position. Polar plots show the distribution of stimulus-arrival phases and their associated local response profiles. Grey bars denote phase-frequency histograms, colored curves denote density-regularized local response scores, and black vectors indicate the preferred phase , defined as the phase-bin center with the maximal density-regularized local response score for each stimulus position. (D) Temporal regularity reduces the deviation between stimulus-arrival phase and preferred encoding phase. The absolute circular phase deviation Dj between the stimulus-arrival phase and is plotted as a function of temporal regularity for each stimulus position. The inset shows the across-position summary of phase alignment as a function of regularity. Statistical annotations report Spearman’s rank correlations between trial-wise temporal regularity and phase deviation. (E) Temporal regularity modulates stimulus-window phase kick. Phase kick Kj was defined as the circular difference between the observed phase advance during the stimulus window and the nominal phase advance expected from f0. Positive values indicate faster-than-f0 phase advancement during stimulus presentation, whereas negative values indicate slower-than-f0 advancement. The inset shows the across-position summary of stimulus-window phase kick as a function of regularity. Statistical annotations report Spearman’s rank correlations between trial-wise temporal regularity and phase kick. (F) Phase kick as a function of stimulus-arrival phase offset. Phase offset was defined as the signed circular distance between the stimulus-arrival phase and the preferred encoding phase . Negative offsets indicate that the stimulus arrived before the preferred phase, whereas positive offsets indicate that it arrived after the preferred phase. Colors denote temporal-regularity groups (R = 1, mid-high R, mid-low R, and low R). Dashed lines indicate zero phase offset and zero phase kick. (G) Temporal regularity reduces the deviation of ISI phase progression from the dominant input frequency. The plotted metric is the absolute frequency deviation | IFISI − f0 |, which quantifies how far the ongoing phase-advance rate deviated from the dominant temporal scale of the input. Smaller values indicate that the internal phase advanced at a rate closer to f0. The first transition was not reliably modulated by regularity (O1→O2, ρ = −0.02, p = 0.69), whereas later transitions showed significant negative correlations with regularity (O2→O3, ρ = −0.15; O3→O4, ρ = −0.22; O4→O5, ρ = −0.24; O5→O6, ρ = −0.20; all p < 0.001). The inset shows the sequence-level summary across ISI windows (ρ = −0.31, p < 0.001). Unless otherwise stated, points indicate regularity-bin means and error bars denote mean ± s.e.m. Statistical tests are two-sided where applicable. Significance levels are denoted as *p < 0.05, **p < 0.01, and ***p < 0.001. For permutation-based cluster tests, significance refers to cluster-level p values.

Oscillatory organization selectively supports temporal encoding during the sample period.

(A) Representative trials illustrating the construction of external temporal phase and its relationship to internal oscillatory phase under low-regularity (R = 0) and perfectly regular (R = 1) input. Shaded regions indicate stimulus presentation windows. The external temporal phase θ(t) advances linearly from 0 to 2π within each inter-onset interval and resets at the next stimulus onset. The internal phase was obtained from the retained f0-band population signal using the same causal complex filter as in the phase-progression analyses. (B) Amplitude-weighted phase-locking value between internal neural phase and external temporal phase as a function of trial-wise temporal regularity. Points indicate regularity-bin means, and error bars denote mean ± s.e.m. across trials (n = 4000 trials). Statistical annotations report Spearman’s rank correlation between R and PLVj for each IOI window. The inset shows the average PLV across IOI windows, , as a function of regularity, with the corresponding Spearman’s rank correlation (permutation test, ρ = 0.78, P < 0.001). (C) Joint heatmaps show the joint distribution of trial regularity (x-axis) and neuronal selectivity (η2, y-axis) for oscillatory power (top row) and phase-locking value (PLV; bottom row). Columns correspond to selectivity for stimulus order, interval order, and stimulus direction. The statistics shown above and to the right of each panel summarize the monotonic trends along the regularity and selectivity dimensions, respectively. Statistics indicate the median model-level Spearman’s rank correlation; significance was evaluated using Fisher’s z transformation followed by a sign-flip test. Unless otherwise stated, data are shown as mean ± s.e.m. Statistical tests are two-sided where applicable. Significance levels are denoted as *p < 0.05, **p < 0.01, and ***p < 0.001.

Rhythmic temporal structure enhances information decodability and prolongs decoding-derived memory persistence.

(A) PLV-binned decoding accuracy analysis during the sample period. Item-level trials were grouped into fixed PLV intervals, and decoding accuracy was averaged within each bin. Higher PLV was associated with higher decoding accuracy for both neuronal activity and synaptic efficacy, consistent with an association between stronger phase locking and improved item-level information readout (Pearson correlation, two-sided, n = 6000 trials, p < 0.001). (B1-B2) Linear regression analysis linking decoding accuracy to the input regularity index (R) during the sample (Left) and delay (Right) periods. The significant positive correlations indicate that temporal regularity was associated with higher-fidelity stimulus information during both encoding and maintenance (Pearson correlation, two-sided, n = 6000 trials, p < 0.001). (C1-C2) Mean decoding accuracy for each sequence item (averaged across the sample and delay periods). At both the neuronal (C1) and synaptic (C2) levels, rhythmic input showed higher decoding accuracy across all 6 positions. This U-shaped pattern is consistent with classic serial-position effects. (D) Decay time constants (τ) of the decoding trajectories for the six items. For each independent network, decoding accuracy was estimated from 1000 trials and fitted with an exponential decay function to obtain item-specific τ values. Rhythmic input were associated with larger decay constants, indicating slower information decay. For each sequence item, rhythmic and arrhythmic conditions were compared across networks using a two-sided Wilcoxon signed-rank test with Bonferroni correction for six item-wise comparisons (n = 15 independent networks). (E) Impact of perturbing different network components on task accuracy. Violin plots show task accuracy under original conditions, shuffled neuronal activity, and shuffled synaptic efficacy. Shuffling synaptic efficacy caused a larger drop in accuracy compared to shuffling neuronal activity. Significance was assessed using a two-sided paired t-test with Bonferroni correction (n = 15 independent networks; corrected p < 0.01 and corrected p < 0.001). Unless otherwise stated, data are shown as mean ± s.e.m. Statistical tests are two-sided where applicable. Significance levels are denoted as *p < 0.05, **p < 0.01, and ***p < 0.001.

Hyperparameters used for network architecture and training.