Mechanistic Interpretability of RNNs emulating Hidden Markov Models
Authors: Elia Torre
Michele Viscione
Lucas Pompe
Benjamin F Grewe
Valerio Mante
Venue: NeurIPS 2025 poster ()
Keywords: Mechanistic Interpretability
Recurrent Neural Networks
Hidden Markov Models
Stochastic Resonance
Reverse-Engineering
Behavioral Neuroscience
Computational Neuroscience
Primary area: neuroscience_and_cognitive_science

Abstract: Recurrent neural networks (RNNs) provide a powerful approach in neuroscience to infer latent dynamics in neural populations and to generate hypotheses about the neural computations underlying behavior. However, past work has focused on relatively simple, input-driven, and largely deterministic behaviors - little is known about the mechanisms that would allow RNNs to generate the richer, spontaneous, and potentially stochastic behaviors observed in natural settings. Modeling with Hidden Markov Models (HMMs) has revealed a segmentation of natural behaviors into discrete latent states with stochastic transitions between them, a type of dynamics that may appear at odds with the continuous state spaces implemented by RNNs. Here we first show that RNNs can replicate HMM emission statistics and then reverse-engineer the trained networks to uncover the mechanisms they implement. In the absence of inputs, the activity of trained RNNs collapses towards a single fixed point. When driven by stochastic input, trajectories instead exhibit noise-sustained dynamics along closed orbits. Rotation along these orbits modulates the emission probabilities and is governed by transitions between regions of slow, noise-driven dynamics connected by fast, deterministic transitions. The trained RNNs develop highly structured connectivity, with a small set of “kick neurons” initiating transitions between these regions. This mechanism emerges during training as the network shifts into a regime of stochastic resonance, enabling it to perform probabilistic computations. Analyses across multiple HMM architectures — fully connected, cyclic, and linear-chain — reveal that this solution generalizes through the modular reuse of the same dynamical motif, suggesting a compositional principle by which RNNs can emulate complex discrete latent dynamics.

Mechanistic Interpretability of RNNs emulating Hidden Markov Models
Elia Torre
Michele Viscione
Lucas Pompe
Benjamin F. Grewe
Valerio Mante
Institute Neuroinformatics,
University of Zurich & ETH Zurich
Abstract
Recurrent neural networks (RNNs) provide powerful approach in neuroscience to infer latent dynamics neural populations and to generate hypotheses about the neural computations underlying behavior. However, past work has focused on relatively simple, input-driven, and largely deterministic behaviors little known about the mechanisms that would allow RNNs to generate the richer, spontaneous, and potentially stochastic behaviors observed natural settings. Modeling with Hidden Markov Models (HMMs) has revealed a segmentation of natural behay- iors into discrete latent states with stochastic transitions between them, a type of dynamics that may appear at odds with the continuous state spaces implemented by RNNs. Here we first show that RNNs can replicate HMM emission statistics and then reverse-engineer the trained networks to uncover the mechanisms they implement. In the absence of inputs, the activity of trained RNNs collapses towards single fixed point. When driven by stochastic input, trajectories instead exhibit noise-sustained dynamics along closed orbits. Rotation along these orbits mod- ulates the emission probabilities and is governed by transitions between regions of slow, noise-driven dynamics connected by fast, deterministic transitions. The trained RNNs develop highly structured connectivity, with small set of "kick neurons" initiating transitions between these regions. This mechanism emerges during training as the network shifts into regime of stochastic resonance, enabling it to perform probabilistic computations. Analyses across multiple HMM archi- tectures - fully connected, cyclic, and linear-chain - reveal that this solution generalizes through the modular reuse of the same dynamical motif, suggesting a compositional principle by which RNNs can emulate complex discrete latent dynamics. Code available at https: / //g thub. www/ Comm/liaaoorr.. waliaaoore/hmmmyr hmmrnn
Introduction
Modern large-scale electrophysiology and imaging techniques provide access to the simultancous activity of thousands of neurons in freely behaving animals, revealing the population-level dynamics underlying perception, cognition, and movement [39 [g
40 In parallel to advances in neural recordings. it has become possible obtain quantitative descriptions of unconstrained, natural behaviors across many timescales [45). The combination of these techniques holds great promise in unlocking the neural computations underlying behavior 7
19
Machine learning critical to interpret and link the high-dimensional neural and behavioral data produced by modern neuroscience experiments. Hidden Markov Models (HMMs) capture unconstrained, natural behaviors and decompose them into sequences of elemental motifs [45). However, their discrete state representation may oversimplify the potentially continuous nature of neural processes. On the other hand, Recurrent Neural Networks (RNNs) provide powerful approach model the dynamics of large neural populations and to generate hypotheses about the underlying Despite this potential, RNNs are primarily used to model highly constrained
computations [10 47
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

tasks that are input-driven anput-driven and largely deterministic. Hence, evidence for RNNs applicability in modelling stochastic discrete state transitions in manner similar HMMs remains limited.
RNNs and HMMs may appear incompatible, as RNNs represent latent processes as trajectories through continuous state spaces, whereas HMMs rely discrete latent states with stochastic tran- sitions between them. Understanding whether RNNs can use continuous dynamics to generate stochastic transitions between discrete states would help bridge this conceptual gap and reduce the need for strong assumptions about the structure of the latent space. To address this question, we develop training approach to directly fit RNNs to HMMs, and then reverse-engineer the trained RNNs to uncover the computations they implement. This strategy may ultimately lead to testable hypotheses about how biological neural circuits may implement discrete behavioral modes through continuous internal dynamics.
This work makes three main novel contributions:
1. training paradigm for stochastic RNN behaviors (Sections Appendix B
Figure unlike
previous work focused deterministic tasks, introduce training method combining noise- driven RNNs with Sinkhorn divergence optimization, enabling learning of the stochastic state transitions typical of HMMs.
2. A demonstration that RNNs emulate HMM statistics (Sections) Appendix C Figures 10 we show that vanilla RNNs accurately replicate the emission statistics various HMM architec- tures, matching both transition dynamics and stationary distributions.
3. A multi-level mechanistic account of how RNNs generate discrete stochastic outputs.
3.1 Global dynamics (Section 4.1 Figures 3
We reveal how training drives the RNNs into a regime where stochastic inputs sustain dynamics along closed orbits ("orbital dynamics"), while in their absence, the activity converges to a single fixed point.
3.2 Local dynamics (Section 4.2 Figure 5): We identify three functional zones: clusters, kick- zones, and transitions, each with distinct dynamical signatures.
3.3 Connectivity and single neurons (Section Figures 6 7 We uncover structured connec- tivity between "noise-integrating populations" and "kick neurons" and validate their role in triggering state transitions through causal interventions.
3.4 Computational principle (Section 6, Figure We show the inferred mechanism relies on self-induced stochastic resonance.
Our findings form a coherent picture: local dynamical features, such as clusters and kick-zones, explain how global orbital dynamics give rise to discrete outputs, while the network's connectivity reveals how individual neurons instantiate these dynamics. Together, these findings show that RNNs discover general computational motif that combines slow noise integration with fast transitions.
The mechanism we uncover offers a novel interpretation of how RNNs can generate stochastic transitions between discrete states. This mechanism (Contributions 3.1 3.44) acts as "dynamical primitive", reusable in a modular fashion to emulate more complex discrete latent structures. Multiple instances of this motif combine to govern transitions between specific pairs of states.
Our results lay the foundation for future work, both in scaling to more complex HMM structures and in investigating whether similar dynamical motifs are employed by biological neural circuits during naturalistic behavior.
Related Work
RNNs have been trained on simple perceptual, cognitive, and movement tasks, both with supervised learning 30
and reinforcement learning 43 46). These RNNs are often trained with few constraints and then reverse-engineered reveal potentially novel hypotheses about the computations performed by biological neural circuits. Fully understanding the dynamics of RNNs and characterizing the mechanisms they implement (i.e., mechanistic interpretability) remains
challenging and open problem [37 36]. [37]6]]
A common reverse-engineering approach focuses on characterizing the topology of fixed points and on obtaining linear approximations of the dynamics around them [1]. For example, RNNs trained on motion discrimination [30]
228 integrate inputs along line- attrattractors.

---

Networks trained to reproduce reach movements implement rotational dynamics similar to that in primate motor cortex 42 and RNNs performing beat continuation generate low-dimensional
oscillatory dynamics 47 Across many tasks, the topology of fixed points is invariant to the details of the RNN implementation (e.g. vanilla RNN, GRU, LSTM), suggesting "universality" in the types of solutions implemented by RNNs (29).
Recent work has pushed mechanistic interpretability beyond population dynamics, to the underlying RNN connectivity. Constraining RNNs weights to being low-rank can reveal the relation between connectivity, computations, and dynamics [31] Such relations may extend to RNNs trained without constraints, as their function may rely on low-rank component of the full connectivity [24, 20]. Highly structured weights may emerge in particular when connectivity and latent dynamics are trained concurrently while imposing biological constraints [25] or when training on many different tasks at
once 22 35
Building on this line of work, here we investigate whether the above insights can be extended to RNNs implementing different type of computation, namely internally-driven, probabilistic transitions between discrete latent states.
3 Approximating HMMs with RNNs
To understand how RNNs encode discrete, probabilistic structure in their state spaces, trained them to replicate the outputs of HMMs, which implement a process that discrete and probabilistic by construction (Fig.
B
c
D
M
At
M
MAAAAM
.5
05
Noise
RNN
Outputs
RNN
Matching
Sinkhorn
Targets
HMM
HMM
Figure 1: Training pipeline. A: At each time-step the network receives i.i.d. Gaussian input vector It ~ (0, La), where the input dimensionality is varied across experiments (d E {1, 10, 100, 200 B: The signal is processed by a vanilla recurrent neural network, h2 = ReLU(h,-1W Ad + xz WT):: the hidden-state dimension is also varied ((h) E (50, 150, 200)). ). Finally, the hidden state is mapped linearly to three-dimensional output via yt A. To mimic HMM emissions, logits are converted to categorical samples via Gumbel-Softmax. C: Parameters θ [Wer. W.., A} are optimised by minimising the Sinkhorn divergence between batches of predicted outputs Y {y+ )r and target sequences Y = { -= D: Target sequences Y* are generated by linear-chain, fhteur-chin,,ulll--coneeee fully--onnectedd and cyclic HMMs (described Section and Appendix A); the example shows a 2-state linear- chain model.
3.1 Families of HMM architectures
To investigate how RNNs encode discrete, probabilistic structure across different latent state topolo- gies, train them to replicate the outputs of HMMs with three distinet types of architectures: finea--hain.,ffllycoooeee fully--onecceeeed and cyclic structures.
Linear-chain HMMs. We first consider a family of linear-chain HMMs that span spectrum from maximally discrete to quasi-continuous representations. These models have M e (2, 3,,3,4,5, 4,5) latent
states ((M) {1,.... . .M } and a common alphabet of observations O = 2,33. Each model HM
fully specified by an emission matrix (MM [c MM and a band-diagonal transition matrix T = [1( M
3

---

For each latent state i, we obtain emission probabilities as linear interpolation between state skewed toward observation 1 = 1) and one skewed toward observation 3 (i = M), while keeping P( = 1) e constant. At each timestep, the chain either remains in the current state or transitions to a neighboring state. Given a change-rate pepee--arameter 0.05, we set QM = pM E (0, 2)
a; = (i- 1)/(A) 1), fix 0.01 and define:
{ j=i,
(1 (1 0:) E
QM:
(M)
(M)
2QM
1<i M, j=i;
(I)
qM
j=itl,
E) 04
0,
otherwise.
This parametrization maintains constant probability p of reaching the most distant state in M steps. When M = 2 the system is maximally discrete; as M increases, the same 0 0++2 continuum is partitioned into progressively finer bins, yielding quasi-continuous latent representation. This design systematically varies discreteness to reveal how RNN dynamics adapt as the latent space becomes more continuous.
Fully-connected HMMs.
To assess generalization beyond simple structures, we use a 3-state fully connected HMM where each state can transition to any other. Each state favors one of the three outputs but retains small probabilities for the others, testing RNNs on fully interconnected dynamics.
We also consider a 4-state cyclic HMM with bidirectional closed-loop transitions. For each state, emissions are biased toward characteristic pair of outputs, one dominant and one weaker component, with adjacent states sharing common output.
Cyclic HMMs.
The complete specification of all HMM architectures, including transition and emission matrices, is provided in Appendix [A] (Fig. 9)
3.2 Training Networks with Sinkhorn Loss Function and Performance Metrics
Network Architecture We employ standard, vanilla RNNs [2] of hidden-state size |/] E (50, 150, 200). At each time-step, the network receives Gaussian input 2N (N La), with d E {1, 10, 100, 200}. The hidden state is updated and projects onto the three logits via:
h ReLU(h-- + w..W.)
V: = M.
To mimic the discrete emissions of an HMM, these logits are converted to categorical samples using the Gumbel-Seftmax reparametrization trick: we add i.i.d. Gumbel noise, divide by a temperature (set to I in all experiments), and apply a soft-max. This continuous relaxation acts as a differentiable proxy for the non-differentiable argmax in the Gumbel-Mar method, letting gradients flow through the sampling step while still converging to one-hot hot vectors as T -
Sinkhorn Divergence. Because our target sequences are probabilistic. we use loss function that is appropriate for comparing distributions: the Sinkhorn divergence [15] 16 38, an optimal transport (OT) distance. Typically, OT-distances find a binary coupling matrix II linking single samples in the outputs and targets, which minimizes the total Euclidean distance between coupled samples. the euclidean distance between coupled samples is zero, the output distribution reproduces the target distribution. Finding such a coupling matrix is computationally expensive and non-differentiable. The Sinkhorn divergence overcomes these issues by finding smoothed coupling matrix allowing non-unique couplings (5). Details on the training regime are provided in Appendix
Performance Metrics. To quantify how closely the trained RNNs replicate their target HMMs, we track four statistics: (i) Euclidean distance between Sinkhorm-aligned output sequences (global re- construction error): (ii) the 3 x 3 transition matrix of successive emissions (long-range range dynamics); (iii) marginal observation frequencies (stationary distribution); and (iv) observation volatility (proportion of time steps with output changes). For all metrics, the trained RNNs replicate the emission statistics of the target HMMs. Definitions and results for all combinations of hidden size (Ih) E 50, 150, 200) and input dimensionality (d E 10, 100, 200) are shown in Appendix C (Figs. 10 TT). Higher input dimensionality improves convergence, while increasing hidden size beyond 150 yields only marginal benefits. Below, we thus focus on the four configurations with [6| = 150 and d = 100.

---

Mechanistic Interpretability: Latent Dynamics
To uncover the mechanisms implemented by the trained RNNs, we first consider their global (Sec- tion 4.1 and local (4.2) latent dynamics. These analyses reveal how several distinct "dynamical motifs" together contribute the RNNs' ability approximate HIMM-like emission statistics (Fig. 2
E
F
G
H
1
Figure 2: A compositional dynamical primitive. Schematic of the dynamical motifs uncovered through reverse-engineering of trained RNNs. A: Sample trajectory in the continuous state space. B: A single fixed point in the absence of input. C: Stochastic input displaces the system from this attractor, sustaining dynamics along closed orbits whose radius grows with input variance. D: Distinct slow regions (clusters) where trajectories linger, separated by transition regions. E: Local "kick-zones" that trigger transitions between clusters. F-I: Schematic of a richer discrete latent structure (F) expressed as a composition of three instances of the same dynamical primitive (G-I). Each instance preserves the same global and local dynamical properties described in panels (A-E).
4.1 Global Latent Dynamics: Noise-sustained Orbital Dynamics
Orbital dynamics under stochastic input. We investigate the global latent dynamics of the trained RNNs by projecting their hidden states onto the first two principal components (PCA, Figure 3 When initialized from random hidden states, and without external input, activity converged to single fixed point (Figs. 2 3), with no evidence for separate attractors that might be expected to represent the discrete latent states of the target HMMs. However, the latent dynamic changes markedly for RNNs receiving stochastic (Gaussian) inputs, to regime where trajectories exhibit orbital dynamics (Figs. A BC). The stochastic input pushes activity outward (Figs. 2C, 3D), while the recurrent component pulls back, together implementing a stable closed orbit along which activity evolves uni uni-directionaalyy Along this orbit, the RNN exhibits regions of slow dynamics (clusters) each corresponding to a distinct output class the reference HMMs, with transitions occurring between them (Fig. (P). Critically, transitions between slow zones result in large changes the output probabilities. As the number of latent states increases in linear-chain HMMs, RNNs do not form additional slow regions along the orbits (App. Fig. 12p. Instead, they capture finer emission discretizations by modulating the alignment of readout axes with the plane containing the orbital dynamics (Appendix L Fig 25). In contrast, for fully-connected and cyclic architectures, RNNs develop multiple orbits connecting distinct pairs of slow regions, reflecting the richer structure of the underlying HMMs. These dynamics reveal compositional dynamical motif in which the same fundamental unit combines to generate increasingly complex latent structures (Fig. 3G and App.
Fig. 13
Emergence of orbital dynamics during training. PCA projections reveal a clear transition in the latent dynamics across training epochs (Fig. AAA purple), marking the shift from a stable fixed point to orbital dynamics. As training progresses, unstable eigenvalues emerge, marking the destabilization of the fixed-point regime. This transition coincides with a rise in eigenvalues near the imaginary axis, indicating the onset of oscillatory activity (Fig IC) and with a characteristic double-dexcent in the loss curve (Fig T6p. Functionally, this transition enables sustained quasi-periodis oscillations that encode the temporal dynamics of the target HMM (Fig. g) The radius of the resulting orbits scales
5

---

RNN Trained on 2 State HNM
RMN Trainod on State Fully Colly-Coneeeted HMH
B
Logk
H
K
M
Figure 3: Global latent dynamics of trained RNNs. Panels A and E illustrate the HMM architectures approximated by the corresponding RNNs shown on the left and right. Hidden-state trajectories are projected onto the first two (or three) principal components of the activity and colored by the dominant logit. B, F: the absence of input, trajectories from random initial conditions converge to single fixed point (cross). C, G: under Gaussian input noise, activity evolves along stable orbits with slow regions associated to distinct outputs; arrows indicate flow. More complex dynamics in G can be decomposed into the same fundamental motif observed in C by applying a second-level PCA within the first PCA space, shown in H-J. D, K-M: contour plots (95% CI) of hidden-state densities under increasing imput variance (o² E {0.1, 1.0, 2.0, 3.0, 4.0)) reveal linear scaling between input variance and orbit radius. Results for all HMM structures are provided in Appendix Figures) 12 13
a relationship confirmed by a perturbative analysis of the RNN dynamics (Appendix[H] Hb. Under unbiased Gaussian input, first-order perturbations average out, while second-order terms scaling linearly with variance dominate after the transition, explaining how stochastic input shapes the emergent orbital regime (Fig 4D). Over training, the network converges to a stable average rate of transitions between clusters (Fig 4B... ). In Section we examine how this behavior arises from the interplay between recurrent connectivity and input noise.
linearly with input variance
PCA of RNN Dynamics across Epochs
B
Average Transition Rate
Complex Eiganvalues
Expacted di
Bare Complex Pigeruatios
Crbitul crust
Transition Rate
0.5
Unutadee
Camplec
0.4
0.3
0.2
0.2
0.0
50
150
150
200
230
50
1.00
150
200
450
50
100
150 Spacts
200
250
PC2
Epochs
Spoc
Figure 4: A: PCA projections of the latent dynamics across epochs for an RNN trained on a 2-states HMM, green and red indicate the dominant logit across epochs. The model first learns single-fixed point, then becomes unstable (purple), and finally transitions to orbital dynamics. B: The RNN learns to encode emission probabilities of the target HMM, by means of specific rate of transitions between the two clusters. C: The RNN becomes unstable just before the transition, then the complex eigenvalues appear on the imaginary axis. D: The expected second order perturbation vector E[dh(2) emerges in the proximity of the transition. This perturbation approximates how the variance of the noise affects the free recurrent dynamics, capturing the dependency the orbital dynamics on the input noise variance. Results validating these findings for all HMM structures provided in Appendix E (Figs. 15).

---

4.2 Local Latent Dynamics: Clusters, Transitions, and Kick-Zones
Having characterized the global dynamics, we zoom in on the local properties of the orbital dynamies to understand where and how the network implements transitions between slow zones. Standard linearization techniques around fixed points (Section 2 offer limited insight here, as the observed convergence towards the single fixed point does not explain the full emergent dynamics (Fig. D. More insights can be obtained by considering how activity evolves in the presence of the input noise on different, short rollouts initialized from given state-space location (Fig. ( Specifically, we segment state-space into different "zones" based on three measures computed at each location: the residency time (RT), defined the average number of steps required in rollout until change in the dominant logit occurs; the logit sign-change count, quantifying how often the dominant logit's 's gradient flips before change in the dominant logit occurs, with lower values indicating more directed flow; and the number of unstable directions, obtained via Jacobian linearization and Möbius tratuuss-taaffommed eigenvalues, which reveals local sensitivity to perturbations (Appendix G
Three dynamical regimes. Segmenting state space by residency time (RT) reveals three func- tional zones, each with distinet dynamical signatures (Fig. 53)). We find similar zones in all RNNs, irrespective of the architecture of the target HMM (Appendix Figs. 17 19.
Clusters (RT > 8): regions where trajectories linger the longest, with frequent logit- agittgrrdieet sign changes (in the range 5-20), and essentially only contracting cigenvalues. These are locally stable regions, each corresponding to different probability distribution over the outputs.
Kick-zones (2 < RT < 8): located downstream of clusters, these regions exhibit moderate logit-gradient sign changes (around 2-4) and few unstable directions locally stretch the flow, indicating a local push away from the cluster.
Transitions (RT < 2): Once trajectories cross the kick-zone, they enter short-lived corridors where the system moves nearly deterministically toward the next cluster. These regions exhibit few logit-gradient sign changes (< 1), reflecting a stable and directed flow.
Noise Sensitivity. To further validate the functional relevance of these regions we explicitly probe their sensitivity to noise (Fig. 533 for the RNN trained on a 2-state HMM; Appendix [ Fig. 20 for all other linear-chain RNNs. We sample initial conditions from both cluster and transition regions and generate trajectories under three noise conditions: identical (y = 0), partially resampled (7 = 0.5), and fully independent (7 1). Transition regions show minimal variability across noise conditions once the kick-zone is crossed, trajectories proceed quasi-deterministcaal] toward the next cluster. In contrast, cluster regions are noise-sensitive: as 7 increases, trajectories diverge, with some crossing the kick-zone and others returning to the cluster. We validated these qualitative differences with quantitative measures of divergence: (i) the trace of the covariance matrix per timestep, capturing dispersion; and (ii) the time-course of the average Euclidean distance between individual trajectories and the mean trajectory (Appendix Fig. 21).
Mechanistic Interpretability: single neuron computations and connectivity
In the previous section described computational mechanism at the level of population-levé dynamics, which relies in particular on transitions between slow regions. Here we aim to explain how these transitions emerge based on key features of the RNN connectivity and the resulting single unit properties. We focus on an RNN trained on the two-state HMM and report the analyses for the other configurations in Appendices and K
Discovering "kick-neurons". Among all RNN units, two separate triplets of neurons have pre- activation values (before ReLU) with a distinctive spatial profile: pre-activation values are strongly negative within clusters, pass through near-zero regime in the kick-zones, and become positive during transitions (Appendix) J Fig. 22). The intermediate regime places the units near the ReLU activation threshold, where small variations in input can determine the opening of the ReLU gate mechanism we describe in the next paragraph. Each triplet is linked to one transition direction, firing at the onset of movement between clusters and emerging the dominant non-zen components of the second-order perturbation vector ([dh (described in Section 4.1 which reflects the network's
7

---

Avg. Expervlluzz Rcx0
(eltee Mobius frensform)
Frequency
6.6
603
0.4
603
603
203
011
10
RT<2
25AT5S RT>8
PC1
Avg. Sign Changes Logits Gradient
[Tarsteenj(rkk/oo]]][luaa
PC1
40
2
Residency Time (RT)
Eigenvalues Re 0
Cluster (V-0)
Custe (v-0.5)
Cluster (y-1)
Transition (v-0)
Transition (y-
Transition
PCI
Figure 5: Clusters, Transitions, and Kick-Zones. A: Left to right: (i) Residency time (RT) reveals slow (gray) and fast (colored) regions; (ii) logit-gradien sign changes histogram shows a bi-modal distribution, separating stable clusters from transitions; (iii) average number of unstable directions per region (via Möbius-trans tranfoored Jacobian), peaks in kick-zones and is lowest in clusters; (iv) spatial distribution of the number of real eigenvalues with positive real part, highlighting localized instability in kick-zone contrasted with the stability of clusters. B: Noise sensitivity analysis from sampled states (black star) in clusters (left) and transitions (right) with 30 trajectories (in purple) generated under increasing noise resampling conditions (7 E 0.5, 1)... Transitions are robust across noise conditions, while clusters exhibit increasing dispersion, indicating higher noise sensitivity. Panels A, B show results for RNNs trained on 5-state HMMs. Results validating these findings for the
remaining configurations are provided in Appendix (Figs. 17 19 and 20
sensitivity input variance. Together, these properties suggest causal role in generating noise-driven "kicks" that initiate state changes, prompting us to term them kick-neurons.
Noise as the trigger. To better understand the mechanich,,,mm we examine the recurrent weight matrix (Waa) and find that kick-neurons within each triplet form mutually excitatory connections, while inhibiting the opposing triplet (Fig. (). Similar to work by [14], showing that input noise integration in RNNs can give rise independent subpopulations, we observe two larger neuronal groups (comprising ~ 70 neurons), each forming recurrent excitatory loops within themselves while projecting inhibitory connections to the other. These populations interact with the kick-neurons through structured excitatory and inhibitory connections, suggesting role as noise integrators that modulate the transition gate opening, hence we refer them as inttee-ittggriing populations.
Causal Interventions. To confirm this mechanism, performed targeted interventions (Figure 7 We modulated either the neurons' activity directly, or the roive-drive they receive through projections from the associated moixe-integrating populations. A modulation factor l controlled both types of perturbations.
Ablation 0). Silencing the kick-neurons or ablating the input noise to the corresponding integrating population traps the trajectory within its current cluster, preventing state transitions. Both manipulations mirror each other, confirming that noise-driven activation is required to open the transition gate. This intervention leads to a loss of critical eigenvalue pairs from the Jacobian spectrum, reflecting the collapse of the orbital dynamics and reversion to single fixed- -oint regime.
Enhancement (p- 2). 2). Doubling kick-neuron activity or amplifying the noise-drive from the cor- responding integrating population causes overshoots beyond the target cluster. These effects resemble increase in noise variance, and the number of critical eigenvalue pairs remains unchanged, indicating that enhancing activity preserves the rotational regime in this setting.
8

---

Pepatation Prpusations
Propuaation
Pepulazio
Topukdoon
Pepulation
Peputtton
Weight
Figure Recurrent weight matrix and Kick-cireuit. A: Left right: (i) Sub-matrix of WhA restricted to the six kick-neurons. Within-triplet weights are positive (red), whereas cross-triplet weights are negative (blue), indicating mutual excitation and reciprocal inhibition. (ii) Weights from noise-inlegrating neurons to kick-newrons (sorted). Each integrating population excites one triplet and inhibits the other. (iii) Recurrent weights within the integrating populations show within-population excitation and cross-population inhibition. B: Kick-circuil: red arrows indicate excitation, blue inhibition. The circuit comprises two self-exciting, mutually inhibiting loops that project to opposing kick-neuron triplets, forming noise-inted integration mechanism that implements stochastic transitions between slow regions. Results for Anear-chain and fully-connected HMMs in App. Figs, 23 24
Popelation
Population
hommin nol
Figure 7: Causal interventions validate the kick-cireuit. A: Latent trajectories in PCA projection. We sample initial conditions (ICs) from the two cluster regions and let trajectories evolve from the same IC with identical input noise draws unperturbed (gray), ablated (blue), and enhanced (red) regimes for 60 time-steps. Suppressing either the kick-neurons (columns: 1, 4) or their noise- drive from integrating populations (columns: 2, prevents transitions between clusters. while enhancement overshoots trajectories past the opposite cluster. Control manipulations on neurons not part of noise-ine iniise-ittegrating populations (columns: 6) preserve kick neurons noise-drive and "cluster- switching". B: Stability signature of the interventions. Mcants.d. number of critical eigenvalue pairs function of the modulation factor μr. Ablation (μ 0, blue) eliminates critical-patt -airras consistent with reversion to single fixed-point regime; enhancement (μ 2, red) maintains the number of mntical--pirss preserving the orbital dynamics. Grey points show intermediate values of μ.
6 Stochastic Resonance
The above analyses show that our RNN exhibits noise-sustained orbital dynamics (Figure 3) driven by structured interactions between two functionally distinct neural populations: slow, noise woise-integrating units, and fast-responding kick-neurons (Figure 5 6). As noise accumulates in the slow subsystem, trajectories drift along quasi-stable manifolds until they reach region the kick-zone which triggers a rapid, deterministic transition to the next attractor-like cluster. This mechanism results in robust, quasi-periodic alternation between network states, even in the absence of any external periodic input. This phenomenon bears resemblance to class of noise-induced dynamics known as self-induced stochastic resonunce (SISR) 33 In contrast to classical stochastic resonance, which requires an external periodic signal 2 SISR arises intrinsically in systems with time-scale separation. Weak noise applied to fast subsystem perturbs the system off slow manifold, initiating

---

deterministic excursions that form coherent, woise-controlie oscillations [48] operate in an analogous regime: the slow populations accumulate stochastic input over time while in the cluster zones, whereas the fast kick-neurony drive sharp transitions once a noise-modulated threshold is reached. This mechanism enables the emergence of stable oscillatory patterns, with their period governed by the interplay between noise variance and slow integration dynamics. The oscillation period, linked in our case to the emission probabilities, can be shaped through learning by adjusting the effective time-constant of the slow subsystem. This interpretation is supported by the population-level dynamics shown in Figure g where each cycle reflects a transition between -itractor-like states, coordinated through the interaction of slow integrators and fast responders whose activity rises sharply before the transition. In this way, the network effectively harnesses internal noise computational signal, leveraging SISR-like dynamies perform structured, probabilistic inference, thus allowing the RNN to emulate the HMM behavior.
Our RNN appears to
Moan Activity of Kick Neurons and Noise-Intggrating Populations
K
Mean Activity
Pepuistor
25
75
125
150 175 200
Time
Figure 8: Oscillatory dynamics. Mean activity of kick-neuron triplets (solid lines) and noise- integrating populations (dashed lines) as the RNN switches cluster regions (alternating gray bands). These populations exhibit clear anti-phase oscillations: while one integrating group dominates, it reinforces itself and its matching kick-neurons, driving the system toward a transition. This dynamic reflects a relf-induced stochastic resonance (SISR) regime [33] 48), where slow noise integration and fast resets ("kicks") produce quasi periodde switching without external periodic input.
7 Conclusion
In this work, we explored how Recurrent Neural Networks (RNNs) can implement discrete, stochastic latent structure through continuous dynamics. We introduced training pipeline that fits vanilla RNNs to various families of Hidden Markov Models (HMMs). Contrary to the expectation of an n-well landscape with one fixed point per HMM state, trained RNNs converge to single fixed point in the absence of input and, under noise, exhibit noise-sustained orbital dynamics: slow regions that encode distinct emission probabilities, separated by short, deterministic transitions. Mechanistically, it relies on two complementary sub-circuits: large noise-integrating populations and fast kick-neurons, whose interplay converts input variance into quasi-periodic transitions. This connectivity structure harnesses internal noise as computational signal, enabling RNNs to reproduce HMM-like probabilistic behavior via a reusable dynamical motif and paralleling features of neural activity observed in the brain.
Cortical activity is intrinsically noisy both micro- and macro-scales, driven by stochastic ion channel dynamics (6), probabilistic synaptic transmission ①1. and the high-dimensional, recurrent nature of cortical circuits, which produces ongoing, noisy background activity even in the absence of external input (13), potentially as a consequence of chaotic dynamics (8). Far from being nuisance, such variability can enhance information processing - increasing sensitivity weak or sub-threshold signals and facilitating the coordination of distributed brain regions in the processing of sensory
information via stochastic resonance
32 34
Our findings demonstrate that unconstrained RNNs can uncover both discrete and continuous latent structure directly from data without any imposed topological priors and reveal that these networks naturally converge toward biologically plausible circuit motifs. Taken together, these results point compositional dynamical primitive in which slow noise integration and fast kick- triggered resets cooperate to generate discrete state transitions, and multiple instances of this motif can combine produce richer latent dynamics. Overall, this work positions Recurrent Neural Networks as a powerful alternative to Hidden Markov Models for modeling the latent structure and neural mechanisms of natural behaviors.
10

---

References
[i] Christina Allen and Charles F Stevens. An evaluation of causes for unreliability of synaptic transmission. Proceedings of the National Academy of Sciences, 91 10380-13333) 1994.
[2] Roberto Benzi, Alfonso Sutera, and Angelo Vulpiani. The mechanism of stochastic resonance. Journal of Physics A: mathematical and general, 14(11):L453, 1981.
[3] Dean Buonomano and Wolfgang Maass. State-dependent computations: spatiotemporal processing in cortical networks. Nature Reviews Neuroscience, (M(2):133-11 2009.
[4] John P Cunningham and Byron M Yu. Dimensionality reduction for large-scale neural record- ings. Nature neuroscience, 17 50:1 2014.
[5] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In C.J. Burges,
L. Bottou, M. Welling, Z. Ghahramani, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.
[6] CE Dangerfield, David Kay, and Kevin Burrage. Stochastic models and simulation of ion channel dynamics. Procedia Computer Science, 1(1):1587-1596, 2010.
[7] Sandeep Robert Datta, David Anderson, Kristin Branson, Pietro Perona, and Andrew Leifer. Computational neuroethology: call to action. Neuron, 104(1):11-24 2019.
[8] Brian DePasquale, David Sussillo, LF Abbott, and Mark M Churchland. The centrality of population-level factors to network computation is demonstrated by a versatile approach for training spiking networks. Neuron, (5):31311-649 2023.
[9] Laura N Driscoll, Krishna Shenoy, and David Sussillo. Flexible multitask computation in recurrent networks utilizes shared dynamical motifs. Nature Neuroscience, 27(7):134--136; 2024.
[10] Daniel Durstewitz, Georgia Koppe, and Max Ingo Thurm. Reconstructing computational system dynamics from neural data with recurrent neural networks. Nature Reviews Neuroscience, 24(11):693-71) 710, 2023.
[11] Lukas Eisenmann, Zahra Monfared, Niclas Göring, and Daniel Durstewitz. Bifurcations and loss jumps in rnn training. Advances in Neural Information Processing Systems, 36:7051 36:70511-70547 2023.
[12] Jeffrey L Elman. Finding structure in time. Cognitive science, 14(2):179-211 1990.
[13] Aldo Faisal, Lue Selen, and Daniel M Wolpert. Noise in the nervous system. Nature reviews neuroscience, 9()):292-305 2008.
[14] Arnaud Fanthomme and Rémi Monasson. Low-dimensional manifolds support multiplexed integrations recurrent neural networks. Neurul Computation, 33(4): 2021.
[15] Jean Feydy, Thibault Séjourné, François- -aanei
Xavier Vialard, Shun-ichi Amari, Alain Trouve, and Gabriel Peyré. Interpolating between optimal transport and mmd using sinkhorn divergences. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2681-2690, 2019.
[16] Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning generative models with sinkhorn divergences. In International Conference on Artificial Intelligence and Statistics, pages 1608- 1617. PMLR, 2018.
[17] Bruce J Gluckman, Theoden Netoff, Emily J Neel, William Ditto, Mark Spano, and Steven Schiff. Stochastic resonance in neuronal network from mammalian brain. Physical Review Letters, 77(19)44098, 1996.
[18] Christine Grienberger and Arthur Konnerth. Imaging calcium in neurons. Neuron, (3)::8862-
885, 2012.
[19] Michael Häusser. Optogenetics: the age of light. Nature methods, 11(10):012210144 2014.
[20] Elizabeth Herbert and Srdjan Ostojic. The impact of sparsity in low-rank recurrent neural 2022.
networks. PLOS Computational Biology,
[21] Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax. arXiv preprint arXiv:161..0144 2016.

---

[22] W Jeffrey Johnston and Stefano Fusi. Abstract representations emerge naturally in neural networks trained to perform multiple tasks. Nature Communications, 14(1):1040, 2023.
[23] Diederik P. Kingma. Adam: A method for stochastic optimization. arXiv preprini
arXiv: 2014.
[24] Renate Krause, Matthew Cook, Sepp Kollmorgen, Valerio Mante, and Giacomo Indiveri. Operative dimensions in unconstrained connectivity of recurrent neural networks. Advances in Neural Information Processing Systems, 35:1 15:17773-17085, 2022.
[25] Christopher I Langdon and Tatiana A Engel. Latent circuit inference from heterogeneous neural responses during cognitive tasks. Nature Neuroscience, pages 1-11, 2025.
[26] Soon Hoe Lim. Understanding recurrent neural networks using nonequilibrium response theory, Journal of Machine Learning Research, 22(47):1-48 2021.
[27] Chris J Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. arXiv preprint vrXiv:/611.007/1 2016.
[28] Niru Maheswaranathan, Alex Williams, Matthew Golub, Surya Ganguli, and David Sussillo. Re- verse engineering recurrent networks for sentiment classification reveals line attractor dynamics. Advances in neural information processing systems, 32, 2019.
[29] Niru Maheswaranathan, Alex Williams, Matthew Golub, Surya Ganguli, and David Sussillo. Universality and individuality in neural dynamics across large populations of recurrent networks. Advances in neural information processing systems, 32, 2019.
[30] Valerio Mante, David Sussillo, Krishna V Shenoy, and William Newsome, Context-dependent computation by recurrent dynamics in prefrontal cortex. nature, 503(7474):78 84, 2013.
[31] Francesca Mastrogiuseppe and Srdjan Ostojic. Linking connectivity, dynamics, and computa- tions in low-rank recurrent neural networks. Neuron, 99(3):609-623, 2018.
[32] Mark D McDonnell and Derek Abbott. What is stochastic resonance? definitions, miscon- ceptions, debates, and its relevance to biology. PLoS computational biology. 555)ee 2009.
[33] Cyrill B Muratov, Eric Vanden-Eijnden, et al. Self-induced stochastic resonance in excitable systems. Physica D: Nonlinear Phenomena, 210(3-4):227-2440 2005.
[34] Alexander N Pisarchik and Alexander E Hramov. Coherence resonance in neural networks: Theory and experiments. Physics Reports, 1000:1157, 2023.
[35] Reidar Riveland and Alexandre Pouget. Natural language instructions induce compositional generalization in networks of neurons. Nature Neuroscience, 27(5):988-999, 2024.
[36] Lee Sharkey, Bilal Chughtai, Joshua Batson, Jack Lindsey, Jeff Wu, Lucius Bushnaq, Nicholas Goldowsky-Dill, Stefan Heimersheim, Alejandro Ortega, Joseph Bloom, et al. Open problems in mechanistic interpretability. arXiv preprint arXiv:2501. 2025.
[37] Jimmy Smith, Scott Linderman, and David Sussillo. Reverse engineering recurrent neural networks with jacobian switching linear dynamical systems. Advances in Neural Information Processing Systems,
2021.
[38] Christos Sourmpis, Carl Petersen, Wulfram Gerstner, and Guillaume Bellec. Trial matching: capturing variability with data-constrained spiking neural networks. Advances in Neural Information Processing Systems, 36:74787 2023.
[39] Micha E Spira and Aviad Hai. Multi-electrode array technologies for neuroscience and cardiol- ogy. Nature nanotechnology, 8(2):83-94. 2013.
[40] Nicholas A Steinmetz, Cagatay Aydin, Anna Lebedeva, Michael Okun, Marius Pachitariu, Marius Bauza, Maxime Beau, Jai Bhagat, Claudia Böhm, Martijn Broux, et al. Neuropix- els 2.0: A miniaturized high-density probe for stable, long-term brain recordings. Science, 37266539) 2021.
[41] David Sussillo and Omri Barak. Opening the black box: low-dimensional dynamics in high- dimensional recurrent neural networks. Neural computation, (5(3):626-649, 2013.
[42] David Sussillo, Mark M Churchland, Matthew T Kaufman, and Krishna V Shenoy. A neural net- work that finds naturalistic solution for the production of muscle activity. Nature neuroscience, 1827)3305-- 2015.
12

---

[43] Jane X Wang, Zeb Kurth-Nelson Dharshan Kumaran, Dhruva Tirumala, Hubert Soyer, Joel Z Leibo, Demis Hassabis, and Matthew Botvinick. Prefrontal cortex as meta-reinfoorcemet learning system. Nature neuroscience, 2166)8600-668 2018.
[44] Jing Wang, Devika Narain, Eghbal A Hosseini, and Mehrdad Jazayeri. Flexible timing by temporal scaling of cortical responses. Nature neuroscience, 21(1):102-116 2018.
[45] Caleb Weinreb, Jonah E Pearl, Sherry Lin, Mohammed Abdal Monium Osman, Libby Zhang, Sidharth Annapragada, Eli Conlin, Red Hoffmann, Sofia Makowska, Winthrop F Gillis, et al. Keypoint-moseq: parsing behavior by linking point tracking to pose dynamics. Nature Methods, (((7):329-11399 2024.
[46] Guangyu Robert Yang and Xiao-Jing Wang. Artificial neural networks for neuroscientisss: primer. Neuron, 107(6): 1048)):044-1000 1070, 2020.
[47] Klavdia Zemlianova, Amitabha Bose, and John Rinzel. Dynamical mechanisms of how an mn keeps a beat, uncovered with a low-dimensional reduced model. Scientific Reports, 14(1):26388. 2024.
[48] Jinjie Zhu and Hiroya Nakao. Stochastic periodic orbits in fast-slow systems with self-induced stochastic resonance. Physical Review Research, 3(3):033070, 2021.
13

---

Contents
Introduction
2
Related Work
3
Approximating HMMs with RNNs
Families of HMM architectures
3.2
Training Networks with Sinkhorn Loss Function Performance Metries
4
Mechanistic Interpretability; Latent Dynamics
Global Latent Dynamics: Noise-sustained Orbital Dynamics
4.2
Local Latent Dynamics; Clusters, Transitions, and Kick-Zones
Mechanistic Interpretability; single neuron computations and connectivity
6
Stochastic Resonance
7 Conclusion
10
Appendix
14
A
Hidden Markoy Model Architectures
15
Training Regime Computational Resources
15
Performance Metries
16
Bankhorn-Alignt Euclidean Distances
18
Transition Matrices Squared Differences
19
Trajectories for RNNs trained on HMMs
20
D.I Trajectories for Linear-Chain HMMS
20
D.2 Trajectories Fully-Connected and Cyclic HMMs
E
Learning Trajectories Across Training Epochs
22
Curves
24
Jacobian Linearization and Stability Analysis via Möbius Transformation
25
Second Order Perturbation
26
Residency times and noise sensitivity
28
Kick-Neurons Activity
33
Recurrent weight matrices and Kick-circuits)
34
L Alignment of Output Dimensions with Orbits Plane
36
Limitations and Broader Impact
37
M Discussion of Limitations
37
Broader Impact
37
14

---

A Hidden Markov Model Architectures
37
8
D
9
5
22
47
A0
25
AB5
A05 685
35
Figure 9: Hidden Markov Model architectures used for training. Panels A-D illustrate spectrum of linear-chain HMMs with M E (2,3, 4, 5] latent states (from left to right, top to bottom). Arrows indicate state transitions, with labels showing transition probabilities; the hidden states of the HMM (larger circles) have been colored according to the most likely outcome associated with a specific state, or by a mixture of them. Smaller circles represent emission distributions for each state over the observation alphabet O = {1, 2,3). 3). with green, blue, and red corresponding to emissions o1, O2, and 03, respectively. All models are constructed to preserve constant probability p 0.05 of reaching the most distant state in M-1 steps. As M increases, the latent representation becomes progressively finer, transitioning from highly discrete (A) to quasi-continuous (D). E: Filly-connected HMM with three latent states, each capable of transitioning to any other (including self self-transitions) Emission probabilities are biased toward one dominant output per state, while retaining smaller weights for the others, producing symmetric, all-to-all transition structure. F: Cyclic HMM with four latent states arranged in a bidirectional loop. Each state emits a characteristic pair of outputs, with one output dominant and the other weaker, while adjacent states share one common output.
B Training Regime and Computational Resources
We conducted a systematic sweep across key architectural and task parameters to evaluate how RNNs replicate the emission statistics of different HMM families. Specifically, we trained networks on four lineur-chain HMMs (M- = 2, 2, 4, 5 latent states) as well as on the fully-connected and cyclic archi- tectures described in Section 3.1 For the linear-chain HMMs, we train networks with three different
hidden state sizes (I4) = 50, 150, 200) and four input noise dimensionalities (d = 1, 10, 100, 200), yielding 48 unique configurations. Each configuration is trained with three independent random seeds, resulting in total of 144 models. For the fully-connected and cyclic architectures, we fixed the hidden size to //| = 150 and input dimensionality to d = 100, training three independent seeds per architecture, resulting in six additional models. Only seeds that successfully converged were retained for analysis.
The training regime is kept consistent across all configurations. Each RNN is trained on 30,000
sequences of fixed length, sampled from its corresponding HMM: 100 for M = 2, 30 for M = 3, 4, 40 for M- 5, 30 for the fully-connected, and 40 for the cyclic architectures. Optimization is performed in batches of 4096 using the Adam optimizer (23 with learning rate of 0.001. Hidden states are initialized to zero at the start of training, and all weights are drawn from uniform
15

---

distribution u(-vk, Vk), where k hidden sige" To stabilize learning and mitigate exploding gradients, we apply gradient clipping with maximum norm of 0.9 for the linear-chain models and 0.3 for the fully-connected and cyclic architectures. Training proceeded until convergence, typically reached within 1000 epochs. However, convergence becomes less consistent as the number of hidden states in the target HMM increases, in which case shorter training sequences were used to facilitate convergence.
The complete pipcline-including HMM data generation, RNN training, computation of evaluation metrics, PCA analyses, and visualization was optimized for efficient execution on modern hardware. On an NVIDIA RTX 4090 GPU, each model completed training in approximately 5-20 minutes, depending on sequence length and network size.
C Performance Metrics
To evaluate how well task-optimized RNNs replicate the behavior of the reference HMMs, we developed four performance metrics that capture both global sequence similarity and fine-grained statistical properties. These include: Euclidean distances between matched sequences, squared differences in transition matrices, observation frequencies, and observation volatility. Below detail the definition, motivation, and implementation of each metric.
Sinkhorn-aligned euclidean distance. To quantify the global discrepancy between HMM and RNN sequences, we compute the Euclidean distance after aligning predicted and reference sequences via the Sinkhorn divergence. This procedure ensures principled pairing between sequences by solving a soft optimal transport problem, where matched sequences minimize expected pairwise transport plan.
distances under antropy-regulane
Given two sequences of categorical outputs. ymu and Y/mm- represented as one-hot vectors in TxC (where 7 is sequence length and C' the number of output categories), we flatten them into vectors of dimension T C and compute the Fuclidean distance:
d( yran: Yumm)
Σ
Ymnn
y (4) 9
Distances are computed across batches of = 5000 sequences, averaged over random seeds and configurations. As baseline, the same metric is computed between pairs of HMM-generated sequences. These distributions are visualized in Fig. 10, demonstrating that RNN outputs closely match the HMM baseline.
Transition matrix squared differences. To assess how well the RNN captures the temporal depen- dencies between successive outputs, we compute empirical transition matrices from the sequences. For each model, we extract the most likely output at each time step and count the empirical frequency of transitions from output i to
# transitions from output i to j Tú
# total transitions from
We compute T and 7hhmm and compare them via element-wise squared differences:
Aij
hhmm
The resulting matrices averaged across models and shown in Fig. providing detailed account of how accurately each RNN reproduces the internal transition structure of its target HMM.
Observation frequencies. This metric assesses the RNN's ability to reproduce the stationary distribution HMM outputs. We count the frequency of each output class across all time steps and sequences, yielding probability vector / E C We compare this to the ground-truth distribution ham and, although not shown in figures, this analysis confirms that most trained networks match long-term HMM statistics closely.
16

---

Observation Volatility. To quantify short-term dynamics, we measure how frequently the RNN changes its output across time. Volatility is then averaged across sequences and compared with that of the HMM. This metric complements observation frequency by detecting under- or over-smoothing in the RNN's emission process.
17

---

C.1 Sinkhorn-Aligned Euclidean Distances
1.40
3.8
46
1.40
4.0
1.38
1.5
14
30-
42
1.25
1.5
1.20
2.0
RMN State
aNH on State
Fufly-Cennacted H99
1.75
10
144
1.70
2.0
74
233
1.00
1.50
7.5
Erguct Dimennsion
Sxput Diasonsiion
Figure 10: Sinkhorn-aligned Euclidean distances. Mean s.d. Euclidean distances between Sinkhorn-aligned output sequences generated by trained RNNs and their corresponding HMM references, across all combinations of hidden-state size (b) E 50, 150, 200) and input-noise dimen- sionality (d E 1, 10, 100, 200). Panels A-D show results for linear-chain HMMs with M 2, 4,5 latent states, while panels E-F correspond to the cyclic and fully-connected architectures, respec- tively. Error bars indicate variability across three random seeds per configuration. Horizontal green bands represent the baseline distance obtained by comparing HMM-ggncrated sequences against themselves. Red-cireled points mark the selected RNN configuration (hh) 150, 100) used in the mechanistic analyses presented throughout the paper. For linear-chain HMMs, increasing input dimensionality systematically improves alignment with the reference emissions. approaching the HMM baseline across hidden-state sizes. For fully-connected and cyclic architectures, RNNs trained with the same configuration (Ih) = 150, 100) achieve similarly close correspondence to their HMM targets.
18

---

C.2
Transition Matrices Squared Differences
on
Hidden Dimension
25
Fully-Connected
17
306
14
15
28
25
28
15
15.
25
200
1
18
12
30
30
26
26
RNN State HNM
RNN on 5 State HMM
2G
RNN on Cyclic HHM
350
Hidden
200
355
20
18
25
Input
Input
0.001
0,030
0.015
0.020
Euclidean Distance
Figure 11: Transition matrices squared differences. Mean + s.d. of the element-wise squared differ- ences between empirical transition matrices computed from RNN outputs and those of their reference HMMs, across all combinations of hidden-state size (h) 50, 150, 200) and input-noise dimension- ality (d 1, 10, 100, 200). Panels A-D correspond linear-chain HMMs with M 2, 3, 4. 5.5 latent states, while panels E-F report results for the fully-connectedd and cyclic architectures, respectively. Errors are averaged over three random seeds per configuration. Darker colors indicate larger devia- tions from the reference transition matrices, while lighter tones reflect closer alignment. Red-outlined boxes highlight the selected RNN configurations (hh 150, d 100) used in the mechanistic analyses throughout the paper. For the linear-chain architectures, transition-matrix alignment im- proves systematically with increasing input dimensionality, approaching near-perfect correspondence at d 2 100. The filly-connected and cyclic models, trained with the same configuration, achieve comparably low deviation levels.
19

---

D Trajectories for RNNs trained on HMMs
D.1 Trajectories for Linear-Chain HMMs
RNN 2 State HMM
RNN on State HMM
RNN State HMM
RNN State HMM
PEL (36A
PGL (52.99)
PCR (96.39)
KC (96,99)
PCi (96.26)
FCI (5.2%)
PCI (48338)
(44.99)
K
PCI
X
Variance
Logit
Varianas 1.0
Varience 2.0
Locit
Variance 3.0
Logh
Varience
Figure 12: Trajectories for linear-chain HMMs. Each column shows results from RNNs (1h 150) trained to reproduce the emission statistics of linear-chain HMMs with M {2, 4, 5} latent states. Hidden-state trajectories are projected onto the first two principal components (same axes across rows). A-D: trajectories from random initial conditions converge to single fixed point (cross) in the absence of input. E-H: under Gaussian input noise, trajectories evolve along stable orbits with distinct regions corresponding to different dominant logits (colors). Arrows indicate average flow direction. I-L: hidden-state density contours (95% CI) under increasing input variance (o² E 10.1, 1.0, 3.0, 4.0}) reveal linear scaling between orbit radius and input variance, while preserving overall shape.
20

---

D.2 Trajectories for Fully-Connected and Cyclic HMMss
RNN on Fully-Connected HMM
(7766))
FOT (00.36)
(85.69)
RNN Cyclic HMM
X
Food Point
Logs:
Logh
Legit
Veienco
M
N
Verenna
Verang
Figure 13: Trajectories for connecco and cyclic HMMs. Hidden-state dynamics of RNNs
(/h) 150) trained to reproduce the emission statistics fully-connected (top) and cyclic (bottom) HMMs. Hidden-state trajectories are projected onto the first three principal components in A-B and I-J. To visualize orbital dynamics linking any given pair of clusters, within the subspace in A-B and I-J we computed second set of PC component separately for activity restricted to each pair, shown C-H and K-N. B, J: Under Gaussian input noise, trajectories evolve along multiple stable orbits connecting slow regions dominated by distinct logits (colors). These orbital trajectories can be decomposed into the same fundamental dynamical primitives described in the text, shown in C-E and K-L, respectively. A, I: In the absence of input, trajectories from random initial conditions converge to a single fixed point (cross). F-H and M-N: Hidden-state density contours (95 % CI) under increasing input variance (o² (0.1, 1.0, 2.0, 3.0, 4.0)) reveal that the linear scaling between input variance and orbit radius preserved.
21

---

E
Learning Trajectories Across Training Epochs
PCA of RNN Dynamics scress Epochs
Average Transition Rate
Complex Bigenvalues
Expected di
0.
Corsegee
Orhice onner
Camplex Equmaalue
Tharettoon Rate
0.5
Complex
0.4
0.3
0.3
0.3
0.2
0.24
RNN on 3 State HMM
1.0
0.5
Compler Eganvelas
0.0
Transkion Rate
0.5
0.4
0.3
0.4
100
200
750
200
250
250
3
Camples Eipenvaluue
0.5
Transitinn Ratt
0.5
0.4
Hoten Liits
0.4
0.3
0.3
State HMM
M
Comples Blammalloss
6.8
Trerettio Bate
0.5
0.4
E Llm
0.4
200
400
1.00
200
300
400
300 400 500
Figure 14: Learning trajectories across training epochs for linear-chain HMMs. Each row corresponds to an RNN trained on a linear-chain HMM with M {2, 4,5 latent states (top to bottom). A, E, I, M: PCA projections of hidden-state dynamics across training epochs reveal the emergence of structured orbital dynamics: networks initially converge to single fixed point, undergo a transient unstable regime (purple), and ultimately form stable, noise-sustained orbital dynamics (colored by dominant logit). B, F. J, N: the average transition rate between output clusters, showing clear shift after the transition epoch (purple dashed line) aligning with the converged transition rate (red dotted line). C, G. K, O: the fraction of unstable and complex dibbius-transfomme eigenvalues of the Jacobian rises sharply around the transition, reflecting, respectively, the destabilization of the fixed point and the onset of oscillatory dynamics. D, H, L,, P: the expected second-order perturbation vector E[dh(2) (2) emerges after this transition, capturing how input noise variance drives the recurrent dynamics and sustains the orbital regime. Together, these results demonstrate a consistent learning trajectory across architectures, in which the network transitions from stable to oscillatory dynamics.
22

---

RNN on Fully-Connected HMN
A
Transition
Eigenvalues
Expected dhift
Complex Expmmaives
Transtonn Rate
0.6
0.
Meddam Unbb
0.5
Eport
0.2
0.3
0.2
300
RNN on Cyclic HMH
H
0.7-
1.0
0.6
0.8
Therstoon Rate
0.4
Complex
0.6
Meddan JIMEN
0.3
0.4
0.1
0.0
0.0
100
300
400
500
200
200 400
500
300
200
300
400
300
Epacha
Epochs
Ebochs
Figure 15: Learning trajectories across training epochs for fully-connected and cyclic HMMs. Each row corresponds to RNN (hhl 150) trained reproduce the emission statisties of fully- connected (top) or cyclic (bottom) HMM. A, E: PCA projections of hidden-state dynamics across training epochs show the progressive emergence of structured orbital dynamics: networks initially converge to a single fixed point, undergo a transient unstable regime (purple), and ultimately form stable, noise-sustained orbital dynamics (colored by dominant logit). B, F: the average transition rate between output clusters, showing a clear shift after the transition epoch (purple dashed line) aligning with the converged transition rate (red dotted line). C, G: the fraction of unstable and complex Möbius- transformed eigenvalues of the Jacobian rises sharply around the transition, reflecting, respectively, the destabilization of the fixed point and the onset of oscillatory dynamics. H: the expected second-on aecoed-order perturbation vector E[dh (2) emerges after this transition, capturing how input noise variance drives the recurrent dynamics and sustains the orbital regime. These results demonstrate that
RNNs trained more complex HMMs follow the qualitative learning trajectory observed for linear-chain models, transitioning from a single attractor to oscillatory dynamics.
23

---

Loss Curves
A
RNN 2 State HMM
RNN 3 State HMM
Training
Training Loss
Validat on Loss
Vo idation Loss
Sinkhorn Loss
RNN 4 State HMM
RNN 5 State HMM
Tha Training Loss
Training
Validation Loss
Validation Loss
Sinkhorn Loss
RNN on Cyclic
on Fully-Connected HMM
Training
Validation Loss
Volication Loss
Sinkham Loss
100
200
Epochs
Epochs
Figure 16: Training and validation loss curves. Sinkhorn loss over training epochs for the six RNN models used the main analyses (h) = 150, d 100), each trained reproduce the emission statistics of HMM with different transition structures: linear-chain M E {2, 4, 5), fully-connected, and cyclic. Training and validation losses (blue and orange) closely overlap across all configurations, indicating signs of overfitting, All models converge within 500 epochs, exhibiting a characteristic couble-descent profile. The second loss drop coincides with the transient unstable regime identified in Figures 14 15 marking the emergence of sustainnd orbital dynamics and the transition from fixed-poin to rotational behavior. This correspondence confirms that the onset of orbital dynamics constitutes the dynamical solution through which the networks minimize the Sinkhorn loss.
24

---

G Jacobian Linearization and Stability Analysis via Möbius Transformation
In vanilla RNN with ReLU non-linearity the hidden state evolves
h (xx Wi 14-1 W
((z) = max(0, 2).
(A.1)
Because the derivative of the ReLU either 0 or the Jacobian that propagates infinitesimal perturbation 0 to 8h, takes a particularly simple form:
J. W AA Dr,
Dr = diag 1 2.00]
2, T. W + 14-1 W Ah-
(A.2)
The diagonal gating matrix Dr switches rows of Wa on or off depending on whether the correspond- ing pre-activation is positive.
Discrete- time stability and the need for a spectral map. In discrete time a fixed point is lo- cally stable iff all eigenvalues satisfy |yx < To interpret this spectrum with the intuition of continuous-time systems (where the relevant boundary is Re X; = 0) we apply the classical Möbius transformation
1 + À
1 1.
(A.4)
1-X'
This bijection maps the unit circle the imaginary axis and preserves complex conjugacy, giving the correspondence
Condition in À
After mapping
Interpretation
A < 1
Re Re ()) < 0
contraction
A
Re
neutral Hopf boundary
A > 1
Re > 0
expansion
25

---

H
Second Order Perturbation
Inspired by the perturbative approach on RNN of [26], here we compute simple terms estimate the impact of the first and second order effects of the input noise on the stable dynamics. Without loss of generality, we reformulate the RNN dynamics as:
h ((h,) 2,W
With representing the ReLU activation, h, IR" the pre-activation recurrent state, 2 RR with Ie ~ (O, In), and 2, the i - th pre-activation neuron. We recover the formulation of equation (2) in Section 3.2 with h; @(ha). Let's define the unperturbed dynamics
Such that Ag. = - for some te > 0. We consider the first order perturbation around hn hg She
By substituting we obtain the following:
h1+1 + 8lle41 = ((h + Sht) + txWWi'e
8h,+1 <(h Sh, (hh) + a.W
Given the ReLU activation, we approximate for small 8h₃
ohh + 5h,) ~ ((h)) + 5h, (x())
where (o(h is diagonal, with entries if (he): > 0, 0 otherwise. Thus:
8h4+1 ~ (b hh + W.W.T. c.
Naming 4, WT and M,,, A = I1 Ar. with A = I,. By iterating obtain:
8h+
Σ ==0
x. W
II As
Tg W WTTMa M
k=st1
==0
hr hg +
3.
W in
Mr.s
==0
We observe that, given the unbiased Gaussian noise, E[5h] 0. Therefore, the RNN integration mechanism must rely higher order terms. Henceforth, we herby derive the second order term as well. To make the higher-order term explicit, Instead of 5h oh()) O(a2) we express δh, into orders of noise:
8h, 5h(") + 5h + O
where 6/ (1) is the first-order term and (1) is the second-order term. For an element-wise activation 8 de 0 for * 1), we get the following recurrent form
dz? So. (i9) e; WTT + x.W 0(1 (8hx)
8h+41 Sh. +
Given that (1 captures only the second order effects, it must satisfy the following relation
26

---

po
sh xe ōh, 2 Do(h?)
Σ
(ôhi )²
(hg)
Es W
Ozp
with 6hg (2) obtain:
This is linear recursion, driven also by quadratic term in oht Solving it
0o. (h:)) ei Mr
8hg (2)
(on!)))
2
dzk
==0
i=1
Which leads to the following expectation, depending only on the noise variance.
2
t
2
80 (h2)
Blöh (2)
d²
Σ w M
e, Mr
2
dz?
i=1
Unfortunately, the ReLU activation does not possess second derivative. Henceforth, approximate the impact of the second order perturbation with the square of the first order perturbation over several trials whenever the pre-activation neurons activity is below O:
(mn 2 + Sk Mc
ahe
Where @ is the Hadamard product and
(2x); <0 0 otherwise
(Sa);
for 10 and averaging with respect 100 different trajectories. This vector emerges after the orbital dynamics onset (Figure 4A.D) and scales linearly with the variance, analogous to the scaling of the orbits observed in Figure B., K-M. Furthermore, empirically, we observe that kick neurons comprise some of the top non-zero components of this vector, suggesting deeper connection between variance and the firing mechanism.
27

---

I Residency times and noise sensitivity
A
RNN on 2 State HMM
500
Avd. Espervaluus Res Ceffer Nicbls Transform)
Frcquenty
sco
0.
100
Sign Changes in Logits Gradient
55
Residency Time
Eigenvalues
RNN on 3 State HMM
1000
Avryy Eigernattes Reso
Cahte Noblus Trandform)
1200
1000
600
0.13
400
0.03
4
5
10
Gradient
12
258Ts88
31
Residency Time (RT)
Eigenva aiwaluess
1200
Avgi Elgenvalues Re-0
(ohtrr Mooius Ticnsform)
1000
a00
200
00
6.0 2.5 S6. 755 10.0 15 17. 20.
RTx2
2:RTTs8
RT
PC1
Avg. Sign Changes Loges
PC
Residency Time (RT)
Eigenvalues Re
D
RNN on 5 State HMM
Ang. Elgervalues Re-0 (ofter Mobius Tronstorm)
Rc
Frequency
1913
003
8
002
g
200
2
Re
RT=8 (Tanettonn/odoeooe)))CCo
RT<2
25RTsB
Changes Gradient
40
Residency Time (RT)
Eigenvalues Re
Figure 17: Residency analyses for linear-chain HMMs. Each row corresponds to an RNN trained on linear-chain HMM with M (2, 4,5) latent states (panels A-D). Left: PCA projections of the latent trajectories colored by residency time (RT) reveal distinct dynamical regimes, with slow regions (clusters, dark gray) where trajectories linger and fast regions (Iransitions, orange) where rapid output switching occurs. Center-left: distributions of average sign changes in the logit gradient exhibit robust bimodal structure separating stable clusters from directed transition flows. Center-right: the average number of unstable directions (via brbiuustrranffrmedd Jacobian eigenvalues) peaks in intermediate RTs, identifying kick-zones that mediate transitions. Right: spatial maps of the number of eigenvalues with positive real part confirm local instability confined to these zones. The coherent structure across all HMM configurations demonstrates that the RNNs converge on common solution: orbital dynamics divided into slow clusters, unstable kick-zones, and rapid transitions, suggesting generic mechanism by which RNNs can emulate discrete stochastic HMMs emissions through continuous dynamics.
28

---

RNN on Fully-Connected HMM
Transform)
1/5
1230
1900
Eigematucs
1.27
190
0.73
500
0.80
Avg. Sign Changes Logit Gradhant
10
(Transitio Conteee
44
Residency (RT)
Eigenvalues
B
(ofter Transform)
1.50
Frequettyy
Eleenalles
2500
1.00
1000
6.79
6.25
RT<2
25RT:88
PCV
Aug. Sign Changes Leghs Cradient
PCI
45
Residency Time (RT)
Elgenvalues Re
Avg. Elgenralies Reci hobius Transform)
2500
1.2
2000
1500
1000
64
RT
25AT
RT=8
Avg. Sign Changes Logits Gradient
Residancy Time (RT)
Biganvalues
Aug. Eigervalues RaSO
(effer Mibius Tre-dfon Trs-sform)
2.00
1.79
1500
1250
1.b
1.00
rog
0.79
500
0.29
14
0.00
Avg. Sign Charges In Logits Gredlent
Residency
Elgenvalues
Figure 18: Residency analyses for fully-connected HMMs. Each corresponds to the principal space or one of the three principal subspaces identified in the RNN trained to reproduce the emission statistics of a fully-connected HMM (panels A-D). Left: PCA projections of the latent trajectories colored by residency time (RT) reveal distinct dynamical regimes, with slow regions (clusters, dark gray) where trajectories linger and fast regions (transitions, orange) where rapid output switching occurs. Center-left: distributions of average sign changes in the logit gradient exhibit robust bimodal structure separating stable clusters from directed transition flows. Center-right: the average number of unstable directions, computed from intermediate RTs, identifying locally unstable kick-zones that mediate transitions. Right: spatial maps of the number of eigenvalues with positive real part confirm that instability is confined to kick-zones while clusters remain locally stable. Together, these results show that the RNN decomposes the fully-connected HMM into three coupled dynamical subspaces, each expressing the same tripartite organization of clusters, kick-zones, and transitions described for the linear-chain architectures supporting the compositional reuse of the same dynamical primitive across HMM families.
Jacobian eigenvalues, peaks in
29

---

RNN on Cyclic HMM
Egenralees Rey (efter Mobiux Transform)
600
25 Sign Changes Logits Gredient
Residency Time (RT)
Elgenvalues
Avr Eigevvalue Rooo
(afurr Wobus Trandorn) (o
3.75
2.30
1.25
1.80
0.75
0.50
0.29
KT<2
25RT28
RT=9
PCL
Avg. Sign Changes in Logits Gradient
PCL
2 Elgenvalues Re 0
Residency Time (RT)
1790
Avo. Eigervatucs Re-0 (aftrr Mobus Transform)
0.8
1:20
750
0.
00
20 Avg. Sign Changes Logits Gradhent
Residency Time
Eigenvalues
Figure 19: Residency analyses. Each row corresponds to the principal space one of the two principal subspaces identified in the RNN trained reproduce the emission statistics of cyclic HMM (panels A-C). Left: PCA projections of the latent trajectories colored by residency time (RT) reveal distinct dynamical regimes, with slow regions (clusters, dark gray) where trajectories linger and fast regions (transitions, orange) where rapid output switching occurs. Center-left: distributions of average sign changes in the logit gradient exhibit robust bimodal structure separating stable clusters from directed transition flows. Center-right: the average number of unstable directions, computed from crobius-transforme locally unstable kick-zones that mediate transitions. Right: spatial maps of the number of eigenvalues with positive real part confirm that instability is confined to kick-zones while clusters remain locally stable. Together, these results show that cyclic architectures preserve the same dynamical primitive observed in simpler HMMs.
Jacobian eigenvalues, peaks in intermediate RTs, identifying
30

---

Cluster (y=0)
Cluster (y=0.5)
Cluster (y=1)
Transition (y=0)
Transition (y=0.5)
Transition (y=1)
RNN State HMM (y=1)
Cluster (y=0)
Cluster (y=0.5)
Transition (y=0)
Transition (y=0.5)
Transition
PC1
RNN State HMM (y=1)
(y=0)
(y=0.5)
(y=0)
(y=0.8)
Transition (y=1)
PC1
Cluster (y=0)
Cluster (y=0.5)
Cluster (y=1)
Transition (y=0)
Transition (y=0.5)
Transition (y=1)
PC1
Figure 20: Noise sensitivity analyses for linear-chain HMMs. Each row corresponds an RNN trained to reproduce the emission statistics of a linear-chain HMM with M E {2, 3, 4, 5} latent states (panels A-D). Left to right: Sampled initial conditions (black stars) taken from representative cluster (left triplets) and transition (right triplets) regions the latent space. From each location, 30 trajectories (magenta)
simulated under three levels of input-noise resampling: identical partially resampled (y = 0.5), and fully independent (? = 1). In all configurations, transitions remain compact and exhibit robust, detrrmitteriiiiiiicc flow once the kick-zone is crossed, whereas clusters show increasing divergence with noise resampling, reflecting high noise sensitivity.
0),
31

---

RNN on State HMM
RNN on State HMM
RNN on State HMM
State HMM
Trace Covariance
30
43
201
30-
15
20
10
10
Dlsar (ym(.))
Ulsar (ym1)
Fansttion (pe0.5)
Avg. Distance to Mean Trajectory
Pansition (pw2)
74
6-
Mean Distance
4-
0 2
0 2
3
Figure 21: Quantification of noise sensitivity for linear-chain HMMs. Each column corresponds an RNN trained to reproduce the emission statistics of a linear-chain HMM with M (2, 4,53 latent states (panels A-D). For each configuration, we quantify the dispersion of latent trajectories initialized from representative cluster (purple and blue) and transition (orange and red) regions under two levels of noise resampling (y = 0.5, 1). A: the trace of the trajectory covariance matrix captures the overall spread trajectories in latent space. B: the average Euclidean distance to the mean trajectory reflects deviations across trials. In all configurations, transitions remain compact and exhibit robust, quasi-deterministic flow once the kick-zone crossed, whereas clusters show progressively larger divergence with increasing noise resampling, reflecting high noise sensitivity.
32

---

Kick-Neurons Activity
Kick Group 1
8
4
Pre-Activatton Value
2
Kick Group 2
3
Figure 22: Pre-activation values of kick kick-neurons, Latent trajectories projected onto the first two principal components (PCI and PC2) of the RNN trained on 2-state HMM. Each point along trajectory corresponds to one time step and is colored by the pre-activation value of specific neuron. The first row shows the three neurons forming Kick Group I, while the second row displays the three neurons of Kick Group 2, which mediate transitions in the opposite direction. In both groups, pre-activations are strongly negative within one of the two clusters, increase to near-zero values in the kick-zones, and become positive during transitions and in the subsequent cluster. These dynamics reveal clear functional role: within clusters, the neurons are fully suppressed (ReLU outputs zero), while in the kick-zones, pre-activation values fluctuate near the ReLU threshold. This intermediate regime places the units near the ReLU activation threshold, where small variations in input can determine the opening of the ReLU gate.
33

---

K
Recurrent weight matrices and Kick-circuits
RNN on 2 State HMM
Kick
Populations
Populations
Populations
2 Populdonn
1ook
KOOk
Kt2k
Mick
Kock
Kick
Population
Population
Population
Population
RNN State HMM
Kick Neurons ant Neurons
Populations
Neurons
Populations
Populations
NOCK
0.1
Kick
Kick
onio omppr.uutio
Kick
Kick
Population
Population
Population
Population 2
RNN
State HMM
Kick Neurons a-2 Kick Neurons
Populations
Kick Neurons
Populations
Populations
Pepulationn
Poputtetio
Population
Population
Popuiation
State
Kick Neurons 2
kick Neurons
Populations
Kick Naurons
Populations
Populations
Pupulation
KOO
KOK
K00k 2
Ktck 2
Population
Population
Figure 23: Recurrent weight matrices and Kick-cireuits for linear-chain HMMs. Each row (A-D) corresponds
RNN trained to reproduce the emission statistics of linear-chain HMM with M E (2,3,4,5, latent states. Left: Sub-matrices of the recurrent weights WAA restricted to the identified kick-neurons. Within-triplet connections cross-triplet connections are inhibitory (blue), indicating mutual excitation and reciprocal inhibition between opposing kick groups. Center: Weights from the noise-inlegrating populations to the kick-neurons (sorted). Each population excites one kick-neuron triplet while inhibiting the other, implementing selective gating of transition direction. Right: Recurrent weights within and across the two noise-interd integratien populations, showing strong within-populan inhibition. Across all architectures. this structured connectivity forms kick-cireuit: two self-exciting. mutually inhibiting loops that project to opposing tick-neurons, enabling noise-driven direction- selective transitions between cluster regions.
predominantly excitatory (red), while
excitation and cross cross- poosspoopulttin
34

---

RNN on Fully-Connected HMM
Kick Neurons
Populations
1OOK
KOK
Kikk
Kick
0.2
Kick
Kick
Population
Population
Population
Population
Kick Nourons - Kick Neurons
Populations
Kick Nouroms
Populations
Populations
Population Pepultion
ICC
Ktck
Kick
Population
Population
Population
Population
Kick Neurons =
Kick Neurons
Populations
Kick Neurons
Populations
Populations
--02
Pooul otion 3" Populaation
Popultatio
Population 1"
Population
Population
Kick Neurons = Kick Ncurons
Kook
Kick
Figure 24: Recurrent weight matrices and Kick-cireuits for fully-connected HMMs. Each row (A-C) corresponds to one of the three principal subspaces identified in the RNN trained to reproduce the emission statistics a fully-connected HMM. Within each subspace, we isolate the kick-neurons and associated noise-ine iaiseeintggraiing populations that together implement self-contained instance of the dynamical primitive described for linear-chain architectures. Left: sub-matrices of the recurrent weights Wa, restricted to the identified kick-neurons. Within-group connections are predominantly excitatory (red), while cross-group connections are inhibitory (blue), indicating mutual excitation and reciprocal inhibition between opposing kick groups. Center: weights from the two noise- integrating populations to the corresponding kick-neurons (sorted). Each population excites one kick group while inhibiting the other, implementing selective gating of transition direction. Right: recurrent weights within and across noise-integrating populations, showing strong within-population excitation and cross poosspoppllttion inhibition, mirroring the architecture observed in linear-chain models. Together, these subspaces (A-C) represent three instances the same kick-circuit motif, confirming that the RNN decomposes the fully-con mully-connected HMM into compositional dynamical primitives that each generate noise-driven, direction-selectivee transitions. The notation population ion (panel B) and population" (panel C) indicates that population groups re-identified independently for each subspace, whereas the the recurrent weights among the three kick-groups show the same organization self-excitation within each group and cross-inhibition between groups revealing that the same circuit principle recurs hierarchically across levels, from single-neuron motifs multi-orbit compositions.
remain fixed across panels. D: at the higher compositional level,
35

---

Alignment of Output Dimensions with Orbits Plane
Output Dimensions Allgnment with Orbits Plane
0.50
Projection on PCA plane
0.45
0.40
0.35
0.30
0.25
States
States
4 States
5 States
Figure 25: Alignment of output dimensions with orbit plane. Mean + s.d. projection of the three output readout axes onto the principal plane spanned by the orbital dynamics, computed across three independently trained RNNs (Ih) 150, d 100) for each linear-chain HMM with M E (2, 3, 4,5] 5] latent states. As the number of latent states increases, the alignment between the readout and the orbit plane systematically decreases, indicating that the networks progressively encode output probabilities through subtler modulations within a shared low-dimensional dynamical manifold rather than by forming additional slow regions distinct attractors. This geometric reorientation provides a mechanistic account of how RNNs approximate increasingly fine latent discretizations transitioning from discrete quasi-continuous regimes by reusing the same dynamical primitive while generating smoother, less skewed output distributions.
36

---

M Limitations and Broader Impact
M.1 Discussion of Limitations
This work focuses on RNNs trained to emulate relatively simple HMM structures. Despite the diversity of HMM families examined (linaar-chain, fully-connected and cyclic), the target models share symmetric transition graphs and emission probabilities constrained to two or three output classes. While this design allows for clear interpretation and mechanistic analysis, natural behaviors are often governed by richer latent structures featuring asymmetric transitions and more diverse emission profiles. Extending the present framework capture these more complex dynamics remains an important direction for future work. Nevertheless, preliminary analyses naturalistic animal behavior suggest that the core dynamical primitive and its underlying mechanisms identified here generalize beyond the simplified settings considered, providing potential foundation for modeling spontancous behavioral sequences biological systems.
M.2 Broader Impact
This work aims to advance our understanding of how Recurrent Neural Networks (RNNs) can represent discrete, stochastic latent structure through continuous dynamics. By reverse-engineering trained RNNs, we uncover mechanistic account of how such representations can emerge in the absence of explicit structural constraints. This contributes to a growing body of work exploring how neural architectures can serve as data-driven models of behavior.
key implication is that RNNs may offer powerful alternative traditional models like Hidden Markov Models (HMMs) for inferring latent structure in naturalistic behaviors. Unlike HMMs, which impose strong assumptions about the topology of the latent state space, RNNs allow for both discrete and continuous structures to emerge directly from the data. This flexibility has the potential to yield more biologically plausible representations of the computations underlying behavior.
More broadly, the dynamical mechanisms we identify particularly how structured connectivity and noise interactions give rise stochastic, quasi-periodic transitions may inspire novel hypotheses about the neural representations of spontaneous, sequential behaviors in biological systems.
We do not foresee any direct negative societal impacts arising from this work.
37

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction in this paper accurately reflect the paper's contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We discuss limitations in the Appendix M.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
38

---

Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: All the results in this paper can be reproduced.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
39

---

Answer: [Yes]
Justification: We provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips... cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run reproduce the results. See the NeurlPS code and data submission guidelines (https /mips.cc//pullc//uuui guwde//ddee/ee/ Cot//ooooeeesiseeeeeeeeee wwwpobmmsssioooooooooooo for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The paper specify all the training and test details necessary to understand the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We report eтror bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments in this paper.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results accompanied by enror bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
40

---

The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of exccution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: The paper provide sufficient information on the computer resources needed to reproduce the experiments in
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics httpp://neurips httpp//s.../ppb.ic/cccccsseeteee
Answer: [Yes]
Justification: The research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there is a special consid- eration due to laws or regulations their jurisdiction).
Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative socictal impacts of the work performed?
Answer: [Yes]
Justification; The paper discuss both potential positive societal impacts and negative societal impacts of the work performed in Appendix M.2
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
41

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The creators original owners of assets used in the paper are properly credited and the license and terms of use are explicitly mentioned and properly respected.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
42

---

If assets are released, the license, copyright information, and terms of use in the thcode has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
latasets.paperssi
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper doess not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
43

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https: httpps////eurpps.... /ctps////ooeeeeecc///2///////// for what should or should be described.
44


Reviews:
[{"review_id": "8HS7w37soH", "reviewer": "Reviewer_qHFi", "summary": "The authors describe their efforts to train and then reverse engineer vanilla RNNs to reproduce the probabilistic dynamics of an HMM.  The authors utilize a range of mathematical techniques from nonlinear dynamics, perturbation theory, dimensionality reduction, machine learning, and computational neuroscience to describe in detail the mechanisms that the trained RNNs deploy and the training dynamics that they display.", "questions": "None.", "limitations": "yes", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 4, "strengths_and_weaknesses": "This paper is exceptional in every regard. It is surprisingly easy to read despite being extremely dense. The authors were both careful and complete in both their experiments and their analysis. The paper addresses an interesting and important question regarding the learning of a stochastic process by an RNN, and the analyses in Sections 4-6 were a feast of methods and details that were seamlessly integrated into a mechanistic description of the RNN dynamics from the population, connectivity, and single-neuron levels. The ease and simplicity of paper belies a mountain of work both in the analysis and the writing.", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time they invested reviewing our manuscript and the very positive feedback."}, {"review_id": "WHvyB2j4DT", "reviewer": "Reviewer_FrtZ", "summary": "In this work, the authors train Gaussian noise-driven RNNs with deterministic and continuous dynamics to behave like stochastic discrete HMMs of a specific class. RNNs are trained directly on the outputs of the HMMs with varying number of hidden states with a distributional loss function. Analyses of trained RNNs show a transition from input-free fixed point to noise-sustained limit cycle globally, with distinct clusters corresponding to the output states of the HMM. Locally, there appears to be 3 distinct types of “zones” that are biased to be within vs. move across the clusters. Finally, the authors identify and causally test single neuron mechanisms that, coupled with the noise input, appear to be responsible for kicking the system state into one of the clusters corresponding to the discrete HMM outputs.", "questions": "1. Is it not relevant to discuss HMMs, stochastic RNNs, and “hybrids” like switching LDS in the related works section? Similarly, perhaps a brief mention of the relationship between HMMs and linear Gaussian systems (stochastic linear RNNs) may be instructive?\n2. What’s the rational for using the specific class of HMMs defined in eq 1? The overall motivation for this is not clear, as well as specific choices, e.g., why reduce the generality of the HMM model, e.g., by limiting the transition matrix to being band-diagonal? (line 76) What happens when these assumptions are violated?\n3. what happens if the middle output class (O2) is removed?\n4. observed state labels in Fig1d shows O1,2,3, whereas everywhere else in text is labeled as 0-2?\n5. line 82: I’m confused what is meant by “the same 0 to 2 continuum is partitioned…” Does this refer to the latent state, which increases with M and is hence further partitioned, or the observation space, which seems to be always 3 discrete states?\n6. Based on the schematic shown in Figure 2, how is the system different from the normal bistable / double well potential? \n7. could you please clarify how Fig 3 row iii shows less skewed distribution as HMM state increase?\n8. I’m not sure about the generality of the SISR claim: does this not depend on the exact structure of the transition and emission matrices? The current setup, even in the 5-hidden state HMM, has a strong preference for the two extreme output classes, both through high-probability transitions into the 1st and last hidden state, and the fact that the middle output state O2 always has very low probability. Do all these findings hold for flatter distributions?\n9. there is some discussion of stochastic dynamics in the brain in the conclusion section, but it’s not at all clear to me how that relates to the current study. Is the idea that brain circuits can be well-modeled by continuous dynamics RNNs, but behavior has been observed through methods such as MoSeq to be composed of discrete set of syllables, so we need to find the equivalence? This is briefly described as the motivation of the study in the introduction, which is informative, but I think much more can be said throughout the paper to concretely relate the current findings to that overarching idea / goal, though I understand that this is difficult to do with a purely computational modeling study.", "limitations": "In my opinion, there is an inadequate discussion of the setup and assumptions underlying the current study, in particular the narrow range of HMMs used and therefore possible observed behaviors.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 1, "originality": 3, "strengths_and_weaknesses": "Strengths\n- it’s in my opinion an interesting investigation into the equivalence of two popular classes of dynamics, RNNs and HMMs, and the mechanisms that enable the former to behave like the latter, i.e., stochastic discrete dynamics\n- the analyses are thorough, covering multiple “levels” of investigation from population dynamics to single neuron mechanisms, as well as causal perturbations\n- the writing and presentation are in general clear and crisp\n\nWeaknesses\n- the findings, while interesting, are based on a very narrow family of HMMs with a specific transition and emission probability structure, and thus unclear how generalizable the insights are (e.g., results in 4.1)\n- the writing is very clear at the sentence level. But structurally, it was difficult to understand and follow the whole picture. Given the density of results and a lack of high-level signpost (perhaps at the end of the intro), it read to me like a list of disparate findings.\n- I’m ultimately not sure what the implication or significance of these results are. I suppose it’s not a given that RNNs can be trained to behave like HMMs. But even accepting that this is a novel finding, in addition to the in-depth investigation into its mechanisms, the fact that it’s a purely modeling study employing a very narrow class of HMMs leaves me wondering what this means for brain circuits.", "was_revised": true, "final_justification": "The work is technically solid and showcases interesting and in-depth mechanistic analyses of the learned RNN. Originally restricted to a narrow class of HMMs, but later expanded during revision, with apparently similar conclusions.", "rebuttal": "Thank you for the helpful feedback. We have substantially revised the manuscript and added new experiments to address your comments.\n\n**Summary:**\n1. Revised *Introduction* and *Conclusion* and added *Contributions*.  \n2. Fitted additional RNNs to more complex HMM structures (*App B.3)* showing a wider applicability of our findings.  \n3. Expanded the *Related Work* section to include suggested references.  \n4. Added a paragraph to *Conclusion* clarifying how the RNN dynamics differs from the N-well potential landscape.  \n5. Clarified the rationale for our HMM design choices in *Section 3.1*.\n6. Added further analyses in response to the other reviewer's feedback (other responses)\n\n**Weaknesses:**   \nWe would like to clarify that the goal of our work is not to provide a comprehensive account of all possible HMM variants. Instead, we focus on a more general question: *Can RNNs, in principle, reproduce the behavior of HMMs, and if so, through what mechanisms?* This question is motivated by work in computational neuroscience, where training and reverse-engineering RNNs has proven to be an effective approach for uncovering how neural circuits, both artificial and biological, implement computations. In several cases, this methodology has led to mechanistic hypotheses that were subsequently tested through experiments (e.g., Sussillo et al. '15 or Hennequin et al. '14).\n\nIn our case, the computations of interest involve *stochastic transitions between discrete states*, a common motif in models of naturalistic behavior (e.g., Calhoun et al. '19). Studying these computations required the development of new training and analysis pipelines, since existing methods primarily target deterministic tasks and do not naturally extend to modeling stochastic, discrete dynamics. To make this motivation more explicit, we revised the *Introduction*, *Conclusion*, and added a *Contributions* section that outlines our goals and coherent narrative.\n\n***Introduction and Contributions***:  \nMachine […] experiments. `Hidden Markov Models (HMMs) can capture natural, unconstrained behaviors as sequences of elemental motifs [46], but their discrete state representation may oversimplify the continuous complexity of neural processes. On the other hand, Recurrent Neural Networks (RNNs) provide a powerful approach to model the dynamics of large neural populations and to generate hypotheses about the underlying computations [10, 47]. Despite their potential, RNNs are mostly applied to deterministic, highly-constrained and input-driven tasks, and evidence for their ability to model stochastic, discrete state transitions like HMMs remains limited`.\n\nRNNs […] them. `Understanding whether RNNs can use continuous dynamics to generate stochastic transitions between discrete states would help bridge this conceptual gap and reduce the need for strong assumptions about the structure of the latent space. To address this question, we develop a training approach to directly fit RNNs to HMMs, and then reverse-engineer the trained RNNs to uncover the computations they implement. This strategy may ultimately lead to testable hypotheses about how biological neural circuits may implement discrete behavioral modes through continuous internal dynamics.`\n\n`This work makes three novel contributions to understand how RNNs can implement discrete, probabilistic computations:`\n1. `A training paradigm for stochastic RNNs (Sec 3.1-2, App A.2): unlike previous work focused on deterministic tasks, we introduce a training method combining noise-driven RNNs with Sinkhorn divergence optimization, enabling learning of the stochastic behaviors typical of HMMs.`  \n2. `Emulation of HMM statistics (Sec 3.2, App A.3, Fig 10-11): we show that RNNs accurately emulate HMM emission statistics across the discrete - continuous spectrum.` \n3. `A novel mechanism for RNNs with discrete stochastic outputs`  \n\t1. `Global dynamics (Sec 4.1, Fig 3-4): We reveal how training induces rotational dynamics and orbital trajectories.` \n\t2. `Local dynamics (Sec 4.2, Fig 5): We identify three functional zones: clusters, kick-zones, and transitions, each with distinct dynamical signatures.`\n\t3. `Connectivity and single neurons (Sec 5, Fig 6-7): We uncover structured connectivity between “noise-integrating populations” and \"kick neurons\" that trigger state transitions.`\n\t4. `Computational principle (Sec 6, Fig 8): We show how the mechanism resembles self-induced stochastic resonance.`\n\n`Our findings reveal a coherent mechanism across levels of description: local dynamical features like clusters and kick-zones give rise to global rotational dynamics, producing discrete outputs. The network’s connectivity reveals how individual neurons instantiate these dynamics. Together, these contributions show that RNNs discover a general computational motif that combines slow noise integration with fast transitions. Preliminary results on more complex HMMs (App B.3) suggest that the identified mechanism may act as a reusable “dynamical primitive”, which RNNs combine in a modular fashion (Contributions 3.1–4) to govern transitions between specific state pairs and approximate more complex discrete latent structures`.\n\n`Overall, this work offers a novel interpretation of how RNNs can generate stochastic transitions between discrete states, contributing to the broader goal of reverse-engineering neural dynamics. These results lay the groundwork for scaling to more complex HMM structures and investigating whether similar dynamical motifs underlie biological computations in naturalistic behavior.`\n\nTo show a wider applicability of our findings, we have added analyses of RNNs trained on HMMs with different graph structures (*App B.3*), which demonstrate that the same mechanisms emerge. In particular, we consider 3 new classes of HMMs:\n1. *Fully-connected*: 3-state HMM, where each state is mainly associated with a different output.\n2. *Cyclic*: 4-state HMMs, with closed-chain transitions and  2 or 3 outputs emitted.  \n3. *Linear-chain*: 3- to 5-state HMMs where each hidden state emits with a high probability one of the 3 outputs or a mixture of them.\n\nUnfortunately, we are not allowed to show the new related additional plots, which demonstrate that all RNNs exhibit the *global and local dynamics* described in the manuscript. Preliminary analyses also reveal analogous *connectivity structure* at the single neuronal level. In models 1, 2, 3, which have 3 outputs, we observe 3 *cluster regions* in the RNN state space and two sets of rotational trajectories, each implementing the transitions between *pairs of clusters*, analogous to those described in the manuscript. We thus hypothesize that the identified mechanism acts as a *“dynamical primitive”*, reusable in a modular fashion to emulate more complex HMM structures. \n\n**Questions**:\n1. Thank you for these relevant references, we have expanded the *Related Work* section to include a discussion on stochastic RNNs, sLDS and linear Gaussian systems.  \n2. We used band-diagonal HMMs to systematically vary discreteness while preserving interpretability (clarified in *Sec 3.1*). To test generality, we added experiments (*App B.3*) with *fully connected*, *cyclic*, and *linear-chain HMMs*. Despite differing structures, all trained RNNs show the same *global and local dynamics* described in the main text.  \n3. We trained additional RNNs without O2 and, given its very low emission probability in the original HMMs, found that its removal has no effect on the latent states or dynamics.  \n4. We have corrected inconsistent labeling (O1,2,3 vs 0,1,2) to O1,O2,O3.  \n5. Here we refer to the fact that as the number of hidden states increases, the discrete jump between states 1 and 2 (in emission distributions of the 2-states) gets divided by intermediate states with interpolating emission probabilities between O1 and O3 (clarified in *Sec 3.1*).  \n6. This is a key question that we now address in the *Conclusion*. We initially expected the RNNs to implement a solution via a double well potential (or N-well potential for N-state HMM), with transitions between wells driven by the noisy input. However, the identified mechanism differs in several ways from such a solution. First, trained RNNs consistently have only a fixed point (Fig 3(i)). Nonetheless, they generate “clusters” of activity characterized by specific emission distributions (Fig 3(ii)). Second, transitions between these clusters are “asymmetric” (different trajectories from cluster \"A to B\" vs \"B to A\"), while in an N-well potential, transitions could occur across the same path in both directions.  \n7. In Fig. 3(iii), trajectories are colored by sampled output (while Fig 3(ii) uses the dominant logit), illustrating how output labels intermix more as HMM state count increases. However, we agree with the reviewer that this can be better evidenced. Hence, we now quantify in the *Appendix* how the alignment of output dimensions with the principal component plane of the orbits systematically degrades as the number of HMM states increases, which leads to a more uniform output distribution.  \n8. As described above, we now also include HMMs with different graph structures. Preliminary analyses of these RNNs show analogous dynamics as those described in the initial submission. Further, we clarified how RNNs implement the “continuous” transition between two “extreme” states in the 5-state HMMs (analyses on alignment of pc-plane/output dimensions, *Question 7*). Together, these basic mechanisms seem to be suited to capture the dynamics of a broad class of HMMs, beyond the examples we present in the manuscript.   \n9. The reviewer correctly points out the motivation for our manuscript, which we now discuss more clearly in the *Introduction* and *Conclusion*.\n\nWe would like to thank you for your review. We have addressed your concerns, please let us know if you have any others in order for your score to change."}, {"review_id": "xFiXONxPCs", "reviewer": "Reviewer_2JTJ", "summary": "This work shows how RNNs can learn the stochastic dynamics of 4 different HMMs where the transition graph between hidden states is a line (the number of hidden states ranges from 2 to 5). The authors show that trained RNNs can reproduce to stochastic state transitions of these HMMs. Examining the projections onto the first two PCs of network activity, they find that the latent dynamics of the RNNs always produce “noise-sustained limit cycles” that are reminiscent of self-sustained stochastic resonance. Then, they examine single-neuron activity and perform targeted perturbation experiments identifying groups of neurons with different causal roles in these “noise-sustained limit cycles”.", "questions": "Weakness 1.1 above is the most important weakness of this work. Addressing this point satisfactorily, either by adding/updating current explanations or by completely removing all reference to the “Hopf bifurcation” explanation throughout the text, could change my rating from \"reject\" to \"weak accept\".\n\nRegarding weakness 1.2, can the authors fix Appendix B and explain, in a self-contained manner, what $\\delta h^{(1)}$ and $\\delta h^{(2)}$ are? I also do not understand what allows the authors to write \"without loss of generality\" on l. 479. Shouldn't the input be inside the nonlinearity $\\phi$? Clarifying these points would also positively affect my final rating.\n\nCan the authors address weaknesses 2, 3, and 4?\n\nRegarding weakness 5, here are some suggestions:\n- Fig 3: add legends for each rows of panels\n- Fig. 14: add titles, and legend to each rows\n- Fig 15: The labeling of the panels is confusing. In addition, a better choice of colors for the lines in the bottom panels would improve readability.", "limitations": "yes.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "### Strengths\nThe modeling setup the authors consider in this work is simple and clearly presented; its rationale is convincing. The text is well organized and the logic of the arguments are easy to follow. Overall, the text is well written. \n\nI find the insights provided by this work on how RNNs can produce dynamics that are seemingly far away from standard RNN dynamics valuable. \n\nOn an algorithmic level, I find it impressive that they can effectively train RNNs using a loss function that compares probability distributions estimated through 30,000 randomly generated sequences.\n\n### Weaknesses\n\n1. The theoretical analysis of the latent dynamics of the trained RNNs appears to contain an important self-contradiction, and the theoretical arguments the authors present are sometimes incomprehensible.\n    1. There is a contradiction between the first and the second paragraph of Section 4.1. In the first paragraph the author say that, in the absence of stochastic inputs, the latent dynamics of the RNNs suggest that there is a unique stable fixed point and no oscillatory dynamics; oscillatory dynamics (what the authors call “limit cycles”) only appear when stochastic inputs are added. This is coherent with Section 6, where the authors claim that this behavior is reminiscent of self-induced stochastic resonance. However, in the second paragraph of Section 4.1, the authors claim that oscillatory dynamics appear because the RNNs go through a Hopf bifurcation. This seems to contradict the explanation of the first paragraph. There are two other issues with the authors’ “Hopf bifurcation” explanation. First, in Muratov et al. (2005), the paper the authors refer to when introducing the notion of noise-induced stochastic resonance (SISR), Muratov et al. write explicitly than SISR is a phenomenon that appears *before* the Hopf bifurcation. Second, I do not understand the analysis presented in Appendix A.6. I have never heard of input-dependent stability analysis (the input being stochastic in this case) and I am not convinced by this approach. (Related to this point, I do not understand what is the y-axis on Fig. 4(c).)\n    2. I could not follow Appendix B, despite spending a fair amount of time trying to understand the arguments. This could simply be due to poor writing. For example, $\\delta h^{(1)}$ appears on line 484 without having been previously defined. It seems that a clear definition of what $\\delta h^{(1)}$ and $\\delta h^{(2)}$ is missing.\n2. On the methodological and theoretical sides, I cannot tell whether this work is original or not. The authors do not say explicitly in what aspects their methods or theoretical analyses are original. \n3. Some words with a precise meaning in dynamical systems are used with a much weaker meaning in the text. For example, the use of the term “limit cycle” when describing the latent dynamics of the RNNs is abusive. The dynamics of the latent variable only *resemble* that of a noisy limit cycle.\n4. The authors mentioned some preliminary results on HMMs with more diverse transition graphs but do not show them.\n5. The presentations of some figures could be improved.", "was_revised": true, "final_justification": "Following the discussion period with the authors, I have changed my rating from 2 to 4. The reason my rating was 2, originally, was because the author used mathematical terms in a way that was imprecise and misleading, which I found inacceptable. The authors have now corrected and clarified the terminology they use and their mathematical arguments; to me, these changes are satisfactory and, in line with what I wrote in my original review, I am updating my rating to 4.", "rebuttal": "Thank you for the careful assessment of our manuscript. Your feedback proved invaluable in refining our work. Below, we provide a summary of the revisions we’ve made to address each of your points.\n\n**Summary**:\n- Removed all references to the Hopf bifurcation and rephrased the relevant sections to describe the RNN dynamics in terms of rotational dynamics\n- Clarified the distinction between the analysis of the trained RNN and the one of the RNN during training\n- Revised the Appendix to clarify and align with the manuscript's notation\n- Added a novel contributions section to the introduction.\n- Corrected the use of the term limit cycle throughout the manuscript\n- Trained additional RNNs to show the generality of our findings with more complex HMM structures\n- Improved formatting of several figures following reviewer’s comments\n- Added further analyses in response to the other reviewers' feedback.\n\n**Weaknesses:**\n- **1.1**: We apologize for the confusion in our original explanation regarding the dynamics of the trained RNNs. The confusion stems from an unclear distinction between the dynamics *during training* and those observed *after training* with fixed weights. At fixed weights the removal of the input noise results in a loss of oscillatory/rotational dynamics and leads to a convergence of trajectories to the stable fixed point. During training, when input noise is always provided, over training steps we observe a transition which leads the RNN to produce closed orbits. In the manuscript, we referred to this transition at training time as a Hopf bifurcation. However we agree with the reviewer that making this analogy leads to a lack of descriptive clarity. In order to fully pursue this analogy we would have to go through an extensive formal proof, which is outside the scope of the current manuscript. We have thus decided to remove references to a Hopf-bifurcation from the paper. In addition, instead of referring to “limit cycles”  we use the more accurate term \"rotational/orbital dynamics\" to describe the onset of these oscillations and report the emergence of the associated purely imaginary eigenvalues (Fig. 4(c)) during training. To clear up the difference between the analysis of the trained, fixed weights RNNs (1st paragraph sec 4.1) and the training time analysis (2nd paragraph sec 4.1), we now divide section 4.1. into two subsections: \n\n\n\t- _`Emergence of Orbital Dynamics under Stochastic Input:`_\n\n\t\t_We investigate [...]. When initialized from random hidden states, and without external input, activity converged to a single stable fixed point, with no evidence for separate attractors that might be expected to represent the discrete latent states of the target HMMs. However, the latent dynamics change markedly for RNNs receiving stochastic (Gaussian) inputs, to `a regime where the fixed point becomes the center of noise-sustained rotational dynamics.` The stochastic input pushes [...] As the number of HMM hidden states increases, the RNNs do not develop additional slow regions along the closed orbits. Instead, they approximate the [...] without requiring additional dynamical structure._\n\t\n\t- _`Onset of Orbital Dynamics during Training:`_\n\n\t\t_`Across training epochs, we observe a qualitative shift in the RNN latent dynamics (see Fig. 4(a)). Initially, all solutions converge to a single fixed point. As the RNN is incentivized to reach both output states, it passes through a transitional phase, leading to a closed orbits solution that cycles between the two output states.` Further evidence for the onset of this rotational behavior is provided by `the emergence of the purely imaginary` eigenvalues of the RNN Jacobian matrix after Möbius-transformation (Appendix A.6). As training [...]. Figure 4(c) shows a marked increase in both the number of eigenvalues near the imaginary axis and the count of unstable eigenvalues, indicating the precise epoch at which the `transition` occurs, which also corresponds to a double-descent of the loss curve [Eisenmann et al. ('23)] (Appendix A.5). `This behavior is functionally significant, as the rotational trajectories enable quasi-periodic behavior that encodes the temporal dynamics of the target HMM (depicted in Fig. 6(c))`. As shown in Figure 3, the radius of the orbit depends linearly on input variance. This relationship [...]. Under unbiased Gaussian noise, the first-order perturbation averages out, while the second-order term — scaling linearly with input variance — dominates `post training transition`, revealing how variance impacts the RNN dynamics (Fig. 4(d)). The network leverages this variance to move `along the orbit` and converges over training [...] input noise._\n\t\n\tWe have also revised the *Abstract* and *Conclusion* accordingly. Regarding Appendix A.6 and Figure 4(c), the plot's eigenvalues were already computed using the Jacobian of the noise-free dynamics, as such we fixed the appendix by removing the input dependency to reflect this fact. For what concerns the axis of Figure 4(c) we now specify that the y-axis refers to the percentage of eigenvalues.\n\n- **1.2**: We appreciate the reviewer’s comments on Appendix B, and we agree that the notation and explanations were insufficiently clear in some areas. We originally adopted the RNN formulation of Soon Hoe Lim ('21). In light of your feedback, we now align it more closely with the formulation and algorithm used in the manuscript, $h'_{t+1} = \\phi(h'_t)W\\_{hh}^T + x_tW\\_{in}^T$ where $W\\_{hh}$ is now outside the activation function, and $h'\\_{t}$ represents the pre-activation hidden state. This should resolve the concern about the loss of generality, as we can recover the formulation used in other sections by setting $h_t = \\phi(h’_t)$. In spite of this revised form, most of the calculations remain the same. The difference is that now, $A_t = D\\_\\phi(h_t^{'o})W\\_{hh}^T$ and the derivative of the activation is computed on the pre-activations. As the actual estimation algorithm was not derived from this solution, the empirical results remain unchanged. We have also made explicit that the $\\delta h$ of the first part of the calculations is referred to as $\\delta h^{(1)}$ in the second part, and we have reported the steps to derive $\\delta h^{(2)}$ more clearly.\n\n- **2**: To clarify the originality and contributions of our work, we have revised the Introduction by adding a Contributions paragraph:\n\t`This work makes three novel contributions to understand how RNNs can implement discrete, probabilistic computations:`\n\t - `A training paradigm for stochastic RNNs (Sec 3.1-3.2, App A.2): unlike previous work focused on deterministic tasks, we introduce a training method combining noise-driven RNNs with Sinkhorn divergence optimization, enabling learning of the stochastic behaviors typical of HMMs.`\n\t- `Emulation of HMM statistics (Sec 4.1, Fig 3-4): we show that vanilla RNNs accurately emulate HMM emission statistics across the discrete - continuous spectrum.`\n\t- `A novel mechanism for RNNs with discrete stochastic outputs`\n\t\t- `Global dynamics (Sec 4.1, Fig 3-4): We reveal how training induces rotational dynamics and orbital trajectories.` \n\t\t- `Local dynamics (Sec 4.2, Fig 5): We identify three functional zones: clusters, kick-zones, and transitions, each with distinct dynamical signatures.`\n\t\t- `Connectivity and single neurons (Sec 5, Fig 6-7): We uncover structured connectivity between “noise-integrating populations” and \"kick neurons\" that trigger state transitions.` \n\t\t- `Computational principle (Sec 6, Fig 8): We show how the mechanism resembles self-induced stochastic resonance.`\n\t\n\t`Our findings reveal a coherent mechanism: local dynamical features like clusters and kick-zones give rise to global rotational dynamics, producing discrete outputs. The network’s connectivity reveals how individual neurons instantiate these dynamics. Together, they show that RNNs discover a general computational motif that combines slow noise integration with fast transitions. Preliminary results on more complex HMMs (Appendix B.3) suggest that the identified mechanism may act as a reusable “dynamical primitive,” which RNNs combine in a modular fashion to govern transitions between specific state pairs and approximate more complex discrete latent structures.`\n\t\n\tOverall, this work offers a novel interpretation of how RNNs can generate stochastic transitions between discrete states, supporting the broader goal of reverse-engineering neural dynamics.\n\t\n- **3**: We followed the reviewer’s suggestions and reworked the use of the terms “limit cycle” and “bifurcation” throughout the text.\n\n- **4**: In the revised submission, we now show a wider applicability of our findings, we have added analyses of RNNs trained on HMMs with different graph structures (Appendix B.3), which demonstrate that the same mechanisms emerge:\n\t- Fully-connected: 3-state HMM, where each state is mainly associated with a different output.\n\t- Cyclic: 4-state HMMs, with closed-chain transitions\n\t- Linear-chain: 3- to 5-state HMMs where each hidden state emits with a high probability one of the 3 outputs or a mixture of them\n\n\tUnfortunately, we are not allowed to show the new related plots, which demonstrate that all RNNs exhibit the *global* and *local dynamics* described in the manuscript. In models 1, 2, 3, which have 3 outputs, we observe 3 “cluster regions” in the RNN state space and two sets of rotational trajectories, each implementing the transitions between *pairs of clusters*, analogous to those described in the manuscript. We thus hypothesize that the identified mechanism acts as a “dynamical primitive”, reusable in a modular fashion to emulate more complex HMM structures.\n\n- **5**: We followed the suggestions of the reviewer to improve the presentation of the figures. \n\nWe would like to thank you for your review. We have addressed your concerns, please let us know if you have any others in order for your score to change."}, {"review_id": "NA8My5OGYk", "reviewer": "Reviewer_Qvyv", "summary": "The paper reverse engineers RNNs trained to match emission probabilities of HMMs. They train the RNNs with Sinkhorn divergence to match the distribution of sequences instead of individual traces. They find that the RNNs switch between two stable regions corresponding to the two outputs that vary in probability. Switches occur when enough noise accumulates in the stable region to push the RNN state to a \"kick-zone\" that pushes the state through a deterministic transition region to the other stable region. In one example RNN, they isolate these computations to particular neurons, specifically two triplets of \"kick-neurons\" that drive the kick-zone to transition dynamics.", "questions": "1. If you were to infer HMM parameters from outputs from the trained RNNs, would you recover the original HMM transition and emission probabilities? I would find this result to strengthen the evidence that the HMMs are successfully embedded in the RNNs beyond the summary statistics that are currently used.\n2. How are the transition probabilities implemented by the RNNs? Is there a way to determine the current discrete state of the RNN?\n3. For the mechanistic analysis, how consistent are the results of the existence \"kick-neurons\", etc. across the individual RNNs trained to match the two-state HMM?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strengths:** The paper is an interesting case study on reverse engineering RNNs implementing probabilistic computations, where most works in this vein have focused on deterministic tasks. The paper analyzes the population latent dynamics in-depth but also explores mechanistic implementation of these dynamics by single units.\n\n**Weaknesses:** As noted by the authors, the HMMs used here are quite simple and thus it is not entirely clear how these insights transfer to more complex, realistic models. In addition, I feel that the analysis is still lacking in explaining the emulation of switching between discrete states and how these transition probabilities are embedded in the network solution. As such, I feel like the explanation of how the RNNs implement the HMMs is pretty incomplete.", "was_revised": true, "final_justification": "The paper introduces a new way of training RNNs to match stochastic models and identifies an interesting mechanism in RNNs for implementing switches between discrete states when driven by noise. It remains unclear to me how these results generalize to more complex HMM structures, since the main text HMMs have a very specific (and arguably problematic) structure, and the new rebuttal results are only described in brief since the authors can't share plots. Still, the paper seems fairly well-done and provides some new methods and ideas that I think are worth sharing.", "rebuttal": "We thank the reviewer for the positive and constructive feedback on the generalizability and the underlying mechanisms of the RNNs. In response, we have made substantial new experiments to address these points.\n\n**Summary:** \n1. Revised *Introduction* and added *Contributions*.  \n2. Fitted additional RNNs to more complex HMM structures in *Appendix B.3* showing the general applicability of our findings.  \n3. Added fitting of HMMs in *Appendix A.3* as additional performance metrics.  \n4. Fitted additional RNNs to 2-state HMMs with different transition probabilities.   \n5. Revised *Conclusion* to include discussion on transition probabilities.  \n6. Added further analyses in response to the other reviewers’ feedback (other responses).\n\n**Weaknesses:**  \nWe would like to clarify that the goal of our work is not to provide a comprehensive account of all possible HMM variants. Instead, we focus on a more specific question, which is currently unanswered in the literature: *Can RNNs, in principle, reproduce the behavior of HMMs, and if so, through what mechanisms?* Addressing this required the development of a new training regime and several different analysis pipelines, since existing methods primarily target deterministic tasks and do not naturally extend to modeling stochastic, discrete dynamics. \n\nThe mechanism we identify, amounts to a novel solution to the general problem of generating stochastic transitions between discrete states. Notably, the RNNs do not implement an intuitive solution based on multiple attractors (energy wells) and symmetric transitions between them, which would show a direct one-to-one mapping with the discrete states of the underlying HMMs. On the contrary, the RNNs converge to a solution implementing rotational dynamics that even with a single attractor creates clusters and asymmetric transitions between them (different trajectories from cluster “A to B” vs. “B to A”). We are not aware of any work suggesting such a mechanism to generate stochastic discrete behaviors in RNNs, and we thus believe that identifying this mechanism, even in simple HMMs, amounts to an important contribution that will guide future experiments and analyses of biological networks during complex behaviors. \n\nTo show a wider applicability of our findings, we have added analyses of RNNs trained on HMMs with different graph structures (*App B.3*), which demonstrate that the same mechanisms emerge. In particular, we consider 3 new classes of HMMs:\n1. *Fully-connected*: 3-state HMM, where each state is mainly associated with a different output.\n2. *Cyclic*: 4-state HMMs, with closed-chain transitions and  2 or 3 outputs emitted.  \n3. *Linear-chain*: 3- to 5-state HMMs where each hidden state emits with a high probability one of the 3 outputs or a mixture of them.\n\nUnfortunately, we are not allowed to show the new related additional plots, which demonstrate that all RNNs exhibit the *global and local dynamics* described in the manuscript. Preliminary analyses also reveal analogous *connectivity structure* at the single neuronal level. In models 1, 2, 3, which have 3 outputs, we observe 3 *cluster regions* in the RNN state space and two sets of rotational trajectories, each implementing the transitions between *pairs of clusters*, analogous to those described in the manuscript. We thus hypothesize that the identified mechanism acts as a *“dynamical primitive”*, reusable in a modular fashion to emulate more complex HMM structures. \n\nWe now revised the *Introduction* and added a new *Contributions (Sec 1)* to clarify these objectives and their novelty (full version in answer to reviewer FrtZ): \n\n`This work [...] :`\n\n1. `*A training paradigm for stochastic RNNs (Sections 3.1-3.2, Appendix A.2): unlike previous work focused on deterministic tasks, [...].* ` \n2. `Emulation of HMM statistics (Section 3.2, Appendix A.3, Figures 10-11): we show [...].`  \n3. `A novel mechanism for RNNs with discrete stochastic outputs` \n\t1. `Global dynamics (Section 4.1, Figures 3-4): We reveal how training induces rotational dynamics and orbital trajectories. ` \n\t2. `Local dynamics (Section 4.2, Figure 5): We identify three functional zones: clusters, kick-zones, and transitions, each with distinct dynamical signatures. ` \n\t3. `Connectivity and single neurons (Section 5, Figures 6-7): We uncover structured connectivity between “noise-integrating populations” and \"kick neurons\" that trigger state transitions.  ` \n\t4. `Computational principle (Section 6, Figure 8): We show how the mechanism resembles self-induced stochastic resonance. `\n\n`Our findings reveal a coherent mechanism: local dynamical features like clusters and kick-zones give rise to global rotational dynamics, producing discrete outputs. The network’s connectivity reveals how individual neurons instantiate these dynamics. Together, they show that RNNs discover a general computational motif that combines slow noise integration with fast transitions. Preliminary results [...] approximate more complex discrete latent structures.`\n\nFor the HMMs described in the original manuscript, in our opinion the description of the mechanism is more “complete” than “incomplete”. Specifically, we do precisely describe the mechanism at various levels of description and how these relate to each other (see contribution list above). In particular, we would argue that none of the results explained in *Section 3 of Contributions* could be a priori predicted. \n\nNonetheless, the reviewer raises an important point about how RNNs implement transition probabilities. To address this point, we trained RNNs on 2-state HMMs with varying transition probabilities and an analysis (*Appendix*) on how output dimensions align with the principal component plane of the orbits. We respond to this point in detail under *Question 2*.\n\n**Questions:**\n1. *Fit quality*. This point allowed us to elaborate on additional metrics for evaluating RNN fits. Beyond the performance measures in *Section 3.2 (i.e., Euclidean distances, transition matrices, observation frequencies and volatilities), Appendix A.3* now includes a “closed-loop” analysis:   \n\t- Firstly, we generate 10k-timestep sequences from both fitted RNNs and ground-truth HMMs, then fit HMMs (with 2–10 states) to estimate the optimal number of states using *Bayesian Information Criterion (BIC)*. For the entire set of RNNs analyzed in the initial submission (fitted to different HMM ground-truths), these  fits converge to an HMM structure with 2 hidden states. This finding highlights a problem of *non-identifiability*.  \n\t- Secondly, to sidestep the issue of *non-identifiability*, we fit HMMs to simulated HMM sequences and RNNs fit to them using the ground-truth number of states, to infer the transition and emission probability matrices. Critically, we observe that the parameters retrieved from RNN- and HMM-generated sequences are very similar between them and they closely match the ground-truth.\n   \n   These findings further strengthen the argument that the statistical properties of sequences generated by RNNs match those of the target HMMs.\n\n2. *Transition probabilities & discrete states*. We address these interesting questions in the *Conclusion* of the manuscript. In simpler HMMs (e.g., 2-state), where a clear 1-to-1 mapping exists between cluster regions in RNN state space and HMM states, transition probabilities are implemented via the rate of transitions between these clusters. As shown in Fig. 4(b), the RNN learns to adjust these rates to match the ground-truth transition probability, a result confirmed across 2-state models with varying transition matrices. We hypothesize that this mechanism relies on the stochastic resonance regime (Section 6), where noise and structured connectivity induce quasi-periodic dynamics that drive switching between these clusters at specific rates. In more complex HMMs (e.g., 5-state), where ground-truth states are similar and a clear cluster-to-state mapping breaks down, the output layer plays a larger role. In the *Appendix*, we show that the alignment of output dimensions with the orbit’s principal component plane degrades as the number of HMM states increases. This misalignment likely leads to more uniform output distributions that enable the emulation of more complex HMM structures and their transition probabilities. \n\n\tRegarding the second part of your question: analyses of neuronal activation profiles (Fig. 16) and ablation studies (Fig. 5) that uncovered kick-neurons also revealed units — likely belonging to noise-integrating populations, Fig. 6 — with activity tuned to specific cluster regions. These neurons could serve as markers of specific cluster regions and that their activation values could be employed to determine the current discrete state of the RNN.  \n\n3. *Consistency across RNNs*. To answer this question, we have extended our analyses of the connectivity structure presented in the original manuscript and added the results as *Appendix A.7*. In the new analyses we show how the structured connectivity discovered in the RNN trained on 2-state HMM generalizes to: RNN trained on 2-state HMM (different random seed) and the whole range of RNNs presented in the initial submission (3, 4 and 5-states). In all cases, we identified the same mechanism of *kick-neurons* and *noise-integrating populations.* Further preliminary evidence on the newly introduced HMM classes (i.e., fully-connected, cyclic and linear-chains) suggests the presence of *kick-neurons* and *noise-integrating populations* also in these cases.  \n\nWe would like to thank you for your review. We have addressed your concerns, please let us know if you have any others in order for your score to change."}]

Author rebuttal:

Meta-review: This paper shows that a certain class of Hidden Markov Models (HMMs) can be well approximated by recurrent neural networks (RNNs).  This represents a novel and important contribution to the literature due to the popularity of HMMs for the analysis of behavior and neural responses, and the fact that the brain itself must implement such discrete latent states via a system with continuous valued units (like RNNs).  I have some concerns about the validity of the paper's claims about the existent of limit cycles, which should (at least in classical dynamical systems theory) be identified with specific a manifold of the deterministic dynamical system, and not dependent on the presence of stochastic inputs.  (That is, it isn't enough to simply examine the distribution of points in PC projections of the noisy high-dimensional trajectories to claim the existence of a limit cycle; rather, one would need to show that nearby trajectories of the (deterministic) system approach this trajectory in the limit t-> infinity.) So I would urge the authors to either make these claims more rigorous (using the dynamics of the deterministic system) or else remove them from the paper. (I would say the same for claims about a Hopf bifurcation). Nevertheless, the reviewers were unanimous in their assessment that the paper makes a worthwhile contribution to the literature and should be accepted.  Congratulations! Please be sure to address all reviewer comments and criticisms in the final manuscript.
Decision: Accept (poster)