Finding separatrices of dynamical flows with Deep Koopman Eigenfunctions
Authors: Kabir Vinay Dabholkar
Omri Barak
Venue: NeurIPS 2025 poster ()
Keywords: Reverse-engineering
Recurrent Neural Networks
Dynamical Systems
Koopman Operator
Primary area: neuroscience_and_cognitive_science
TL;DR: We present a neural network-based method using Koopman Eigenfunctions to efficiently identify separatrices in high-dimensional dynamical systems, validated on synthetic, ecological, and neural network models.

Abstract: Many natural systems, including neural circuits involved in decision making, are modeled as high-dimensional dynamical systems with multiple stable states. While existing analytical tools primarily describe behavior near stable equilibria, characterizing separatrices -- the manifolds that delineate boundaries between different basins of attraction -- remains challenging, particularly in high-dimensional settings. Here, we introduce a numerical framework leveraging Koopman Theory combined with Deep Neural Networks to effectively characterize separatrices. Specifically, we approximate Koopman Eigenfunctions (KEFs) associated with real positive eigenvalues, which vanish precisely at the separatrices. Utilizing these scalar KEFs, optimization methods efficiently locate separatrices even in complex systems. We demonstrate our approach on synthetic benchmarks, ecological network models, and high-dimensional recurrent neural networks trained on either neuroscience-inspired tasks or fit to real neural data. Moreover, we illustrate the practical utility of our method by designing optimal perturbations that can shift systems across separatrices, enabling predictions relevant to optogenetic stimulation experiments in neuroscience. Our code is available on [GitHub](https://github.com/KabirDabholkar/separatrixLocator), and we share an interactive description of the work and its extensions in a [UniReps blog](https://unireps.github.io/blog/2025/Separatrix-Locator/).

Finding separatrices of dynamical flows with Deep Koopman Eigenfunctions
Kabir V. Dabholkar
Faculty of Mathematics
Technion Israel Institute of Technology
Haifa, Israel 3200003
kabir@campus technion. ac.
Omri Barak
Rappaport Faculty of Medicine and Network Biology Research Laboratory
Technion Israel Institute of Technology
Haifa, Israel 3200003 omri bar.kkgmail......... com
Abstract
Many natural systems, including neural circuits involved in decision making, are modeled as high-dimensional dynamical systems with multiple stable states. While existing analytical tools primarily describe behavior near stable equilibria, charac- terizing separatrices- the manifolds that delineate boundaries between different basins of attraction remains challenging, particularly in digenniimesii settings. Here, we introduce a numerical framework leveraging Koopman Theory combined with Deep Neural Networks to effectively characterize separatrices. Specifically, we approximate Koopman Eigenfunctions (KEFs) associated with real positive eigenvalues, which vanish precisely at the separatrices. Utilizing these scalar KEFs, optimization methods efficiently locate separatrices even in complex systems. We demonstrate our approach synthetic benchmarks, ecological network models, and high-ding -igh-dimensional recurrent neural networks trained on either neuroscience- inspired tasks fit real neural data. Moreover, we illustrate the practical utility of our method by designing optimal perturbations that can shift systems across separatrices, enabling predictions relevant to optogenetic stimulation experiments in neuroscience. Our code is available on GitHubl and we share an interactive description of the work and its extensions UniReps blog
Introduction
Recurrent neural networks (RNNs) widely used in neuroscience models of computation arising from the coordinated dynamics of many neurons, motivating efforts to reverse-engineer their underlying dynamical mechanisms 2 In partieular, many cognitive tasks such decision-making 3 and associative memory [4] can be modeled as multistable dynamical systems, where distinct decisions or memories correspond different stable attractor states in phase space. Transitions between these attractors are governed by the geometry of the basins of attraction and, crucially, by the separatrix: the manifold that delineates the boundary between basins (Figure TA].
reverse-engineering method that has yielded significant insights about RNN computations involves finding approximate fixed points and linearising around them [5]. This involves minimizing scalar function- the kinetic energy q(x)) =
locate these points (Figure B). Once found, the
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

linearisation of the dynamics at the fixed point can shed light on the mechanism of computations 16
However, fixed points alone do not capture the global organization of multistable dynamics. Since inputs perturb the state in arbitrary directions, it is critical to know whether they cross the separatrix. To predict the effects of perturbations or design targeted interventions, one must characterize the separatrix itself.
Ideally, we would have a scalar function analogous to the kinetic energy-smootte yet vanishing precisely the separatrix (Figure IC). This would allow gradient-basee optimization to locate the decision boundary and enable the design of optimal decision -hanninnn perturbations (Figure [TA).
In this work, propose novel method to characterize separatrices in high-dimensional black-box dynamical systems by leveraging Koopman operator theory [7 18
Specifically, we approximate scalar-valued Koopman eigenfunctions (KEFs) with positive real eigenvalues using deep neural networks. These eigenfunctions vanish precisely on the separatrix.
We apply this framework to synthetic systems, ecological models, RNNs trained on neuroscience- inspired tasks, and trained to reproduce neural recordings. In addition, demonstrate that the learned KEFs can be used to design minimal perturbations that push the system across separatrices setting relevant experimental protocols such as optogenetic stimulation.
A
á = f(x)
B ((x) = [x]]]
C
v(x)) =
gradient descent
separatrix
gradient descent
input
optimal
choice
Figure 1: (A) Phase-portrait of a 2D bistable system. The two attractors can signify different choices, and therefore the direction between them called the choice direction. External input pushes the system across the separatrix letting it relax to the other attractor. The optimal perturbation has different direction. (B) The kinetic energy vanishes at the fixed points (green 'o' stable and '+' unstable) but does not reveal the full separatrix.
(C) We aim to learn a scalar function (xxx) that vanishes precisely on the separatrix. Gradient descent from random initial points yields a numerical method to locate the respective minima of the scalar functions.
We summarise our main contributions:
We develop tool to locate separatrices, the surfaces between basins of attraction in black- box multi-stable dynamical systems: a gap in the RNN reverse-engineering toolkit.
We demonstrate that KEFs with positive eigenvalues vanish precisely on the separatrix and can be trained using deep neural networks and a loss based on the Koopman PDE error.
We identify two degeneracies of the Koopman PDE and propose effective regularization strategies to resolve them.
We show how the learned KEFs can be used to design minimal norm perturbations that shift the system across separatrices.
We demonstrate the method on systems ranging from low-dimensional models to 668- dimensional RNN fit to mouse neural data.

---

2
Related Work
Our work builds growing body of literature at the intersection of Koopman operator theory, deep learning, and the analysis of dynamical systems, particularly in neuroscience and machine learning.
Koopman theory has recently been used to evaluate similarity between dynamical systems, both in neuroscience [he] is applied to study the temporal structure of computation- in machine learning, where it has been used to compare training dynamics across models [20]. These approaches typically analyze system-leve behavior using dynamic mode decomposition 21 23 finie-dimenniionl approximation of the Koopman operator.
In parallel, deep learning methods have emerged as powerful tools for solving partial differential equations (PDEs) in high-dimensions. Notably, the Deep Ritz Method (24] and Deep Galerkin Method (DGM) (25] which eliminate the need for meshes of points. related line of work uses physics-informed neural networks (PINNs), which incorporate known physics (often PDEs in fluid
dynamics) as part of the loss function during DNN training (26).
Koopman-based embeddings have also been proposed as a tool for analyzing the internal dynamics of RNNS. In 277] the authors show that eigenvectors of finite-dimensional approximations of the Koopman operator can uncover task-relevant latent structure in RNNs. More generally, several works explore DNN- approximations of Koopman operators for learning meaningful embeddings of nonlinear dynamics 30).
An alternative line of research for identifying Lagrangian Coherent Structures (LCS) employs the Finite- Time Lyapunov Exponent (FTLE) [31,3], which quantifies sensitivity to initial conditions by measuring the exponential rate of separation between nearby trajectories over finite time horizon. Ridges in the FTLE field reveal stable and unstable LCS, with the latter corresponding to separatrices in our terminology. These methods are most often applied to two- or three-dimensional fluid flows, where they delineate dynamically distinct regions in the velocity field.
Finally, our approach is conceptually connected to work the geometry of Koopman eigenfunctions themselves. In particular, (34) studics the level sets of KEFs and their relationship to isostables and isochrons in systems with stable fixed points. In the setting of linear systems, and nonlinear systems topologically conjugate to them, [35 36) establish theoretical links between KEF level sets and separatrices (stable manifolds). Together, these studies motivate approach of deep-learning KEFs, as a method for identifying separatrices in general high-dimensional, -tuble-ttable systems.
Results
3.1 KEF's as Scalar Separatrix Indicators
We consider autonomous dynamical systems of the form:
attractor
separatrix
(xxx)
t f=x,, x E X
(1)
where the is shorthand for the time derivative d D and f X + X' defines the dynamics on di
attractor
0
an N dimensional state space X.
Our goal is to construct a smooth scalar func- tion that vanishes precisely on the separatrix between basins of attraction. Consider two such basins (Figure 2). All we care about is the exis- tence of the separatrix manifold, with dynamics moving away from it. The simplest such dynam- ies is one miedimensional linear system gó Au, with A > 0. This motivates mapping that projects a E X' to v²()) E R, and induces these dynamics. Such mapping needs to satisfy:
a = f(() x E X
X
2 E
R
Figure 2: Mapping the high-dimensional dynamics of æ with an unstable manifold - the separatrix to the one dimensional linear dynamics of & with instability at db = O, ie.. À > 0. We approximate the mapping (xx) with a DNN.
d dt (M(a(e)) (x)()()
(2)
3

---

along any trajectory zt)) in X. See Appendix for formal link between &b and the separatrix.
This is precisely the behavior of Koopman eigenfunction (KEF) with eigenvalue À 0. Note that Koopman eigenfunctions are usually introduced a different manner, and Appendix shows the connection our description.
Equation can be re-written as:
Ve(()) Xx((a).
by employing the chain-rule of differentiation and requiring hold at as E X. This is the Koopman partial differential equation (PDE). relates to the timescale of & and is an important hyperparameter of our method (Appendix (1.
We approximate v using deep neural network (Appendix F and train it by minimizing the Koopman PDE residual. Specifically, we define the loss:
PPD pp))) (v(x) f(z) Ax)(==))
where p(ar) is sampling distribution the phase space 24 25 As with any eigenvalue problem,
this loss admits the trivial solution u = 0. To discourage such solutions, we introduce a shuffle- normalization loss where the two terms are sampled independently from the same distribution:
Shuullle
E Rm=ppaa). [Vi((z)) f(x)) {x)(2))
(5)
and optimize the ratio:
LPDE
Catio
(6)
Shullele
We train using stochastic gradient descent, where expectations are approximated by a batch of samples drawn from p(œ) and the shuffle corresponds to random permutation of the samples in the batch (see Appendix G for details).
To illustrate the method, we start with an analytically solvable system in ID (Figure 3A):
i = I 3
The system has three fixed points, corresponding minima of ((x) (Figure 33). A KEF can be derived analytically (Appendix A
(x))
(8)
x
And the zero this function corresponds to the unstable point, which serves as a separatrix this ID case. Figure 3C shows that the DNN approximates this function well, with the location of the zero (separatrix) being captured precisely.
We also upply the method to two 2D bistable systems: 2D damped Duffing oscillator (Figure 3DEF). and 2-unit GRU RNN trained on one-bit flip-flop task GGH), In both cases, the system has two stable fixed points (green circles) and one unstable saddle (green crosses). Kinetic energy functions, shown for comparison, are minimized at the fixed points. contrast, the learned À KEFs are zero on the separatrix (green contours).
3.2 Challenges and Solutions
While the examples above show cases where simple optimization leads to the separatrix, there are several crucial implementation details of our proposed methods. In particular, even v0(x) that satisfies the Koopman PDE may fail to identify the true separatrix. This arises from known degeneracies in Koopman eigenfunctions, particularly in multistable or high-dimensional systems. To enable utilization of our tool, we describe two key failure modes and our strategies to resolve them, summarized in Figure 4

---

A
dæ
B
g(x)
C
(x)
dt
1D bistable
DNN analytical
D
E
F
Damped Duffing
Oscillator
1 22
2 22 σ1
21
3
2
0.0
0.5
1.0
G
2-D GRU RNN, 1-bit flop flop
H
22
2
0.0
0.5
(
[v(x)]
Figure 3: Our method to approximate KEFs three bistable systems. (A) ID system à x² The curve shows (xx), and its fixed points in light green "o's stable and 'x's unstable. (B) The kinetic energy ((x) of this system. (C) the true KEF (8) and its DNN approximation obtained by our method. (D,E,F) Damped Duffing Oscillator in 2D (G.H.I) 2-unit GRU 37 RNN trained on 1-bit flip flop (IBFF) and our KEFs.
Degeneracy across basins. central issue stems from the compositional properties of Koopman eigenfunctions. Let 4/1 (z) and (12(x) be eigenfunctions with eigenvalues A₁ and Ag. Then, their product is also KEF:
(()((2))2)))) f(x) = (A1 + A2) (x))
In particular, consider a smooth KEF 1 with À = 1 that vanishes only on the separatrix (c.g., as in Figure 3). Now, consider piccewise-constant function 6o with X = 0 that takes constant values within each basin and may be discontinuous at the separatrix. The product yi 4o remains a valid KEF with À 1, but it can now be zero across entire basins thereby destroying the separatrix structure we aim to capture (Figure top).
We observe this behavior empirically in Appendix D where independently initialized networks converge different spurious solutions. To mitigate this, we introduce a balance regularization term that biases e to have nonzero values in opposing basins, encouraging sign changes across the

---

Challenges:
Solutions:
À
Demand (v) 0
0
0.5
Choose p(x) to bias v
Figure 4: Top: In the presence of multiple basins, KEF can collapse to zero within single basin. This degeneracy is realised by multiplying the KEF with a piecewise constant KEF with X and invoking ( This example corresponds to ă
... We introduce regularisation term (I0) to encourage the mean value (अ)) ~ O. This encourages solutions with sign changes across basins. Bottom: In higher dimensions, degeneracy arises from directional ambiguity in solutions. We visualise the analytical solution (12) for Ø = x z²; = y². We address this by sampling from multiple local distributions around separatrix points and training ensemble of KEFs.
separatrix. Specifically, we define:
(E)( Varly
Cbal
(10)
and train using the combined loss Liatio + "hat Chat, where "hat is scalar hyperparameter.
In higher-dimensional systems, the Koopman PDE admits family of valid KEFs that differ in their directional dependence. Consider a separable 2D system:
i Si (z),
ý = -2(y)-
(11)
Solving the PDE for this system (appendix yields family of KEFs parameterised by lr E R::
(b. v) ((x))((())+>))>))
Byy)
(12)
where A(z) and (y) are KEFs the respective ID problems. For example, when μ 1. the cigenfunction depends only on and ignores therefore unable to capture y-dependent separatrices. Figure 4 (bottom) illustrates this effect: different values of μ yield KEFs aligned with different separatrices.
Even in non-separable systems, this degeneracy arise. Optimizing Lisa lead to KEF that identifies some separatrices and ignores others (Appendix DD. To address this, we train multiple KEFs (va(+)) while using their input distributions bias each one capture different separatrix. For (1(2)), we choose two points in different basins of attraction, and then use a binary search on the line connecting them to find point on the separatrix. Note that the KEF is still needed obtain the full separatrix, and not just point. Around each such point 3,, we define a local distribution A of I). using range of scales {oiy)
span both fine and global structure. For each distribution, minimize the sum Σ C We then consider the union of the separatrices obtained from each of the KEF's to complete the picture (see 2-bit flip flop demonstration below).
3.3 Demonstrations
We demonstrate the applicability of the method several qualitative examples.

---

A
B
C
Figure 5: Two-bit flip flop task in a 3-unit GRU. The system has stable fixed points (light-green points). (A.B) Two KEFs obtained by our method. They complement each other as they each captures separatrix along one direction. (C) Use of KEF design minimal perturbations that push trajectories across the separatrix.
3D GRU RNN Performing Two-Bit Flip Flop
We first demonstrate our method on low-dimensional recurrent neural network trained to perform two-bit flip flop (2BFF) task. Specifically, we use a 3-umit gated recurrent unit (GRU) network (37). The trained network exhibits four stable fixed points (Figure 5 corresponding to different memory states of the task.
To overcome the degeneracies described in Figure adopt targeted sampling strategy. We first identify points on the separatrix by interpolating between pairs of fixed points and performing binary search: each step, we simulate the dynamics to determine basin membership and refine the search. Around these discovered separatrix points, we construct concentric isotropic Gaussian distributions, and sample from them train on the loss Liotal (Appendix G).
Two resulting KEF are shown in Figure A,B). As expected, the KEFs vanish precisely along the separatrices. This result validates the ability of our method to recover boundary manifolds in neural dynamical systems, even in the presence of degeneracy. Once we know the separatrices, we can determine optimal perturbation directions (Figure Starting from a given initial condition (red star), we see that the same amplitude perturbation is sufficient to reach a different attractor when using the separatrix information, and insufficient when directed the desired attractor. A more quantitative depiction of this effect is shown below in a higher-dimensional system.
11D Ecological Dynamics
KEF Values in 2D Reduced Space
10
We next apply our method (Appendix G high-dimensional ecological model: general- ized Notka-Volterre (gLV) system fit to genus- level abundance data from mouse model of cutibittic-inducer Clostridioides difficile infec- tion (CDI) (38). The system has five stable fixed points. Following [9] we focus our analysis to two of these fixed points representing healthy and discased microbial states.
0.1
We optimize the KEF in the full 11-dimensional state space. For interpretability, we follow the projection approach of 39), visualizing the dy- namies in the 2D plane spanned by the two cho- sen stable fixed points and the origin (see Fig- ure G). Although the KEF is trained entirely in the original l-dimensional space, its zero level
0.1
0.0
0.2
6.4
0.6
0.8
1.0
Figure 6: KEF approximation in fitted 11D gLV model of CDI [8 39.. Zero level set of the KEF aligns with the separatrix 2D projection plane. set (light green curve) aligns well with the true separatrix (orange line) computed using a grid of initial conditions in the 2D slice (39].
7

---

A
B
C
bistable oscillations
(r- 2) 2 2)
= 1
2
0.0
0.5
logIo(((=))
(x)
Figure 7: Applying our method system stable and unstable limit cycles, system without any fixed points on the separatrix. (A) system equations. (B) kinetic energy, with dashed line for separatrix. (C) KEF from our method with zero level highlighted.
This result demonstrates that our technique can be applied directly to real-world fitted models, without dimensionality reduction training time.
Limit cycle separatrix
We test our method in setting where there are fixed points along the separatrix. We construct system which oscillates at fixed frequency (0 1), but converges one of two preferred amplitudes = - 2) - (r )) The system has three limit-cycles, two of them stable (r 1, 3) and one unstable (r 2). In Figure 7 visualise the flow, its kinetic energy and the limit cycles. The system has no fixed points, and thus fixed point analysis is futile. We utilize Radial basis function neural network [40] to parameterise the KEF (Appendix FD.
We show that our approximation the KEF recovers the separatrix r (Figure 7
668D RNN fit to mouse neural activity
To demonstrate our method in high-dimensional (see Appendix H for scaling results) and neurosci- entifically relevant setting, we applied it to recurrent neural network (RNN) trained to reproduce mouse neural activity from Finkelstein al. [9] The trained RNN exhibits bistability between two memory states. As in dimenssoaan systems, we first located point on the separatrix by performing a binary search along the line connecting the fixed-point attractors, simulating the dynamics at each step to determine basin membership.
In the original experiment, mice were trained to respond to optogenetic stimulation of their sensory cortices, and the RNN was fit the peristimulus time histogram of recorded neural activity. The network undergoes bifurcation from monostability to bistability as function of an external ramping input Wext- For analysis, we fixed Mext 0.9 within the bistable regime and trained the Koopman eigenfunction (KEF) network (Appendix G using samples drawn from isotropic Gaussian distributions centered at the separatrix point.
Because of the digh-dimensionality. direct visualization of the learned KEF and dynamics is not feasible. Instead, we validated the model using curve-based evaluation approach. We construct multiple Hermite polynomial curves that interpolate between the two stable fixed points. The curvature of each curve parameterized by a random vector, and each is defined by parameter
E [0, 1/, where a 0 corresponds to one attractor and a 1 to the other (see appendix E Because each curve continuously connects the fixed points, must cross the separatrix. Figure shows 2D PCA projection of several such Hermite curves. Crucially, the actual curves span the entire 668D space. We simulate dynamics from 100 points along each curve and determine their final basin to infer where each curve crosses the separatrix, forming ground-truth truth reference.

---

A
B
V
0.0
1.0
Curve parameter
C 1.0
D
Towards target fixed point
Towards random points
on separatrix
point a
Optimal direction
R² = 0.941
0.5
||
0.0
0.0
0.5
1.0
10
20 30
40
true separatrix point a
distance from base point
Figure 8: Validation of KEF approximation 668D RNN fit to mouse neural activity 9 (A)
PCA projection of Hermite curves between fixed points coloured by true basin labels and separatrix points along the curves (red crosses). (B) KEF values along three Hermite curves versus curve parameter a, well the true separatrix point along the curve. (C) Comparison between true and predicted separatrix positions along eurves. (D) Perturbation amplitudes ||], i.e., distance from Tbase perturbation targets. The KEF-guided solution yields the smallest perturbation crossing the separatrix.
Next, evaluate the learned KEF along these same curves. Figure 33 shows KEF values along sample Hermite curves as a function of a, with the zero crossing indicating our predicted separatrix. Figure [8C compares the a-locations of the ground truth and the KEF- KEF-predieted separatrix points. We observe strong agreement, indicating that the learned KEF reliably tracks the separatrix this high-dimensional system.
Finally, we demonstrate how the KEF can be used to design minimal perturbations that shift state across the separatrix (similar to Figure Figure 5C). In general, this involves input-driven dynamics & f(x u), with time-varying inputs u(t). To demonstrate the utility of the method we study a specific, simplified scenario in which the input is a strong instantaneous perturbation (t)) Aō(t), which moves the state 20+)) x(0_) A, after which the dynamics evolves according to f(x)) f(x 0).
Given a base point x (0_) = Dbsses just before perturbation, we aim to solve:
A arg min
|| subject to (a Passe + A) 0.
(13)
A
Taking advantage the differentiable nature of our method, we use the Adam optimizer to find A from random initialization. Figure shows that indeed the optimised perturbation A smaller, compared to simply aiming the perturbation towards the target fixed point, or towards random points the separatrix.

---

4
Discussion
We presented a novel framework for identifying separatrices in digh-dimensional, black-box dynami- cal systems using Koopman eigenfunctions (KEFs). This method is particularly useful for analyzing recurrent neural networks (RNNs), which are commonly used to model neural computations involving multiple stable states.
Prior efforts in reverse-engineering RNNs relied heavily on locating fixed points and linearizing dynamics locally [616]. While powerful, these methods cannot directly capture global structures or predict system responses to large perturbations that cross basin boundaries. By directly approximating scalar-valued KEFs that vanish precisely separatrices, our method complements and extends existing local linearization approaches. Practitioners can use our KEFs alongside ixed-poin analysis to achieve a comprchensive understanding of the dynamical system's 's landscape.
While FTLE methods B31 also identify separatrices as ridges of finite-time trajectory divergence, they may change sharply near separatrices while providing little gradient elsewhere. We speculate that this could limit their usefulness in nigh-dimensional systems, where gradient-basee -batd localization is needed. Moreover, differentiating FTLE requires differentiating through the dynamical function, which may be computationally expensive or even infeasible due to vanishing/exploding gradients. In contrast, our Koopman eigenfunction framework, by integrating globally, provides a smooth scalar field whose zero level set identifies the separatrix, enabling efficient gradient-baser searches without repeated forward simulations of the target system.
Our work also advances the application of Koopman operator theory dynamical systems. Previous studies primarily utilized Koopman eigenfunctions to predict or control dynamics within a single Likewise, methods comparing dynamical systems to one another use the dynamic mode decomposition which does not always discern between different basins 19 20 46). Such studies usually involve KEFs associated with negative eigenvalues (A<0), which exhibit opposite behavior to ours: they explode at separatrices and approach zero at attractors. In contrast, we specifically targeted eigenfunctions associated with positive eigenvalues, ensuring their zeros correspond exactly to separatrices.
basin of attraction 28 45
To help practitioners use our method, we highlight inherent challenges, such as degeneracy in the Koopman PDE. To overcome these, we introduced a specific regularization-a balance term ensuring eigenfunctions change sign across differrntt basins and show how to choose distributions for the Koopman PDE. These ideas build on existing work on KEF approximation 30 29 47 enabling reliable identification of separatrices in diverse and nigh-dimensional systems.
Our method provides an alternative to more direct approaches for locating separatrices, such as grid searches or bisection methods that repeatedly simulate the ODE from many initial conditions (39]. While learning a KEF involves iteratively solving a PDE over phase space, trajectory- approaches scale with simulation time and often revisit the same regions of phase space. In contrast, solving the PDE resembles dynamic programming and can be made more efficient (Appendix (H), par- ticularly near separatrices where critical slowing down makes ODE-based methods computationally demanding.
While demonstrated the applicability of method to diverse scenarios, we do not provide theoretical guarantees linking the accuracy of the KEF approximation and that of the separatrix location. Furthermore, like techniques for finding fixed points [5 48), our method requires knowing the dynamics in the entire phase space. Extending this to trajectory-based methods facilitate separatrix inference directly from neural data.
50 can
An interesting extension of our work is to stochastic dynamical systems [51 52 While one can approximate separatrices using only the deterministic component of the dynamics, full treatment of stochasticity and its effect on basin boundaries remains open. Likewise, data-driven models with uncertainty in the inferred flow [53] raise questions about how this uncertainty propagates to the separatrix. Since separatrices often occupy sparsely sampled or unstable regions of phase space, this motivates active sampling sampling-potentiallyy through targeted optogenetic stimulations to improve model accuracy along these boundaries [54].
In conclusion, we hope that by focusing separatrices, our method could inform intervention strategies in neuroscience, ecological or engineering systems, providing general-purposs tool to predict and control transitions between stable states in complex dynamical landscapes.
10

---

Acknowledgements
This work was supported by the Israel Science Foundation (grant No. 1442/21 OB) and Human Frontiers Science Program (HFSP) research grant (RGP0017/202) to QB). The funders had no role in study design, data collection and analysis, decision publish, or preparation of the manuscript.
References
[1] Omri Barak. Recurrent neural networks as versatile tools of neuroscience research. Current Opin- ion Neurrbiology, 46:1-6, October 2017. ISSN 0959-4388. doi: 10.10016f.combb2211.0..000
URL httt:t:/s///ww http:/wwww.ciincceirre....... com scissciiiie
[2] Saurabh Vyas, Matthew D. Golub, David Sussillo, and Krishna V. Shenoy. Computation
Through Neural Population Dynamics. Annual Review of Neuroscience, 43(1):249 275, July
2020. ISSN 0147-006X, 1545-4126. doi: 10.1 040144//mn---mm--000060666066600
URL
https:/ bthps::/wwwrr /twp:://wwww../. annualreviews. htt.://wwwwiii..rrro/o///oooo/// 1 1146/ -nwurer-eeeooo--022111 099111
[3] Christian K. Machens, Ranulfo Romo, and Carlos D. Brody. Flexible Control of Mutual
Inhibition: Neural Model of Two-Interval Discrimination. Science, 30)
February 2005. ISSN 0036-8075 8075, 1095-9203 doi: 10. 111///cience URL https
//www.science.. www.ncc.ic..................ee.. wttp:.ccience.rg/doo//oo///o
w1t/////e///ce..// 1104171
[4] J Hopfield. Neural networks and physical systems with emergent collective computational abilities. Proceedings of the National Academy of Sciences, (9(8):255422555 April 1982. doi: 10.1073/pnas.79.8.22
URL / /ww. www.wwww..aa..mm..mmmom... pna..paa..... org http:/doi/aaaaa/ae 10/3///0..00///m .99 8.2554 Publisher: Proceedings of the National Academy of Sciences.
[5] David Sussillo and Omri Barak. Opening the Black Box: Low-Dimensionall Dynamics in Biigh-Dimensionae Recurrent Neural Networks. Neural Computation, 25(3):626 649, March 2013. ISSN 0899-7667.
doi: 10.1162/NECO a 00009. URL httpt: / oww.oii.orgg1...11 1162/
NECO_a_ 00409
[6] Federico Carevaae,, Victor de Lafuente, Ranulfo Romo, Omri Barak, and Néstor Parga. Dynamic Control of Response Criterion in Premotor Cortex during Perceptual Detection under Temporal Uncertainty. Neuron, 86(4):1067-1077, May 2015. ISSN 0896-6273, doi: 10. hw001//j.neroon2015....... articleepa3/500566771555333
URL, httpp:://ww..ciieeeetee scwencedirect..... com/science/
[7] Niru Maheswaranathan, Alex H. Williams, Matthew D. Golub, Surya Ganguli, and David Sussillo. Reverse engineering recurrent networks for sentiment classification reveals line attractor dynamics. Advances in Neural Information Processing Systems, 32: 32:15696-15705,
December 2019. ISSN 1049-5258.
[8] Niru Maheswaranathan, Alex H. Williams, Matthew D. Golub, Surya Ganguli, and David Sussillo. Universality and individuality in neural dynamics across large populations of recurrent networks. Advances in Neural Information Processing Systems, 1199:1562911660 December 2019. ISSN 1049-5258.
[9] Arseny Finkelstein, Lorenzo Fontolan, Michael N. Economo, Nuo Li, Sandro Romani, and Karel Svoboda. Attractor dynamics gate cortical information flow during decision- making. Nature Neuroscience, 24(6))843-850, June 2021. ISSN 1097-6256, 1546-
1726. doi: 10. 00.1038/41199--0210000000 151598-021-008800
URL https: htttps:/wwwy http:///wwwn.au...eeee nature. com/articles.
[10] Valerio Mante, David Sussillo, Krishna V. Shenoy, and William T. Newsome. Context-dependent computation by recurrent dynamics in prefrontal cortex. Nature, 503(7474):78 84, November 2013. ISSN 1476-4687.
URL httpp::/www.w.... articles/nature127/0 Publisher: Nature Publishing Group.
doi: 00.1038/natur
com/
[11] Mengyu Liu, Aditya Nair, Nestor Coria, Scott W. Linderman, and David Anderson. Encoding of female mating dynamics by hypothalamic line attractor. Nature, 634(8035))901-999 October 2024. ISSN 1476-4687, doi: 10. 158-3/444118--00--799119
URI. /httpst:s/.. //www
nature. com/ Aarticlese 841584-156680 020-888---44-0000 079166
Publisher: Nature Publishing Group.

---

[12] Laura N. Driscoll, Krishna Shenoy, and David Sussillo. Flexible multitask computation in recurrent networks utiliizes shared dynamical motifs. Nature Neuroscience, 27(7): 1363. July 2024. ISSN 1544-1726. doi: 1038//41593-0000 024-01668--4444
http:///ww.nnuure.eeeeeee nature Publisher: Nature Publishing Group.
com /ttp/arrcccee///1555-242222
[13] Isabel IC Low, Lisa M Giocomo, and Alex H Williams. Remapping in recurrent neural network model of navigation and context inference. eLife, 12:RP86943, July 2023. ISSN 2050-084X. doi: 10.7554/eLife. 86943. URL httthttt:://da org 7554..5554//11. 86943 Publisher: eLife Sciences Publications, Ltd.
[14] Warasinee Chaisangmongkon, Sruthi K. Swaminathan, David J. Freedman, and Xiao-Jing Wang. Computing by Robust Transience: How the Prontt-Parietal Network Performs Se- quential, Category- -Based Decisions. Neuron,
March 2017. ISSN URL Intpss :tttps:/wwwrr //www. cell com/neuron/ abstract/S00996 6273 (17)30 X Publisher: Elsevier.
0896-6273. doi: 10. www../...o.o.oooooooooooo
[15] Marino Pagan, Vincent D. Tang, Mikio C. Aoi, Jonathan W. Pillow, Valerio Mante, David Sussillo, and Carlos D. Brody. Individual variability of neural computations underlying flex- ible decisions. Nature, 639(8054)421--42,, March 2025. ISSN 1476-4687, doi: 10. 1038/
s4 1466-4244-08 URLhhttps://// http///wwwaatt..../t www.:natueeecmmm coatureecom/areeeee http://www.../...../e//////-----
Publisher: Nature Publishing Group.
[16] Jing Wang, Devika Narain, Eghbal A. Hosseini, and Mehrdad Jazayeri. Flexible timing by temporal scaling of cortical responses. Nature Neuroscience, 21(1):102-110, January 2018 8. ISSN 1546-1726, doi: 10. (0000///1113301--000000
URL. https: //wwww wwwure.ttee..mmmm com/ Publisher: Nature Publishing Group.
articles/s41593-0000 017-
[17] B. O. Koopman. Hamiltonian Systems and Transformation in Hilbert Space. Proceedings of the National Academy of Sciences, 17(5)315--31, May 1931. doi: 10.1033/pnas.1 17.5.31..
URL/https:// h/ww::///wwww.. pnas. wwwg//og//oio /tt/////./..110/ppaa..... 1073//0a...p..... 17 515 Publisher: Proceedings of
the National Academy Sciences.
[18] Marko Budišié, Ryan Mohr and Igor Mezić. Applied Koopmanism. Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(4))04710, December 2012. ISSN 1054-1500. doi: 10.1063/ 1.4772195. URL httpht /oi. 1ttg:/0..1066/// 4772195 177199
[19] Mitchell Ostrow, Adam Eisen, Leo Kozachkov, and Ila Fiete. Beyond Geometry: Compar- ing the Temporal Structure of Computation in Neural Circuits with Dynamical Similarity
Analysis. Advances in Neural Information Processing Systems, 36:33824 -333377, Decem- ber 2023. URL
www.reeedddnnnne.eeen.e......e.. cttpppecccppppppeeepp files paper / wttpe///pppe////////aaaa 2023/hash// 8cc60799222994445577733366222223 Abstract- wow.eetteoooeeeeeeee.o html
[20] William Redman, Juan Bello-Rivas, Maria Fonoberova, Ryan Mohr, Yannis G. Kevrekidis, and Igor Mezić. Identifying Equivalent Training Dynamics. Ad- vances in Neurul Information Processing Systems. 37:23603-23629, December
cwturpps.ccpppe www.slillippeee httpees//pper//2///aaaae 2024. URL https: htttps:/prroceedingg Abstract- Conf html 2.07333.faf6ccc00aabeeeeeeeeeea. -tt.e///etett-ttttttttt-----
[21] Peter Schmid. Dynamic mode decomposition of numerical and experimen- tal data. Journal of Fluid Mechanics, 656:5-28, August 2010. ISSN 1469- 7645. 00022-112)
cambridge org/ doi: 10.1017/550221120000001
URL. https://ww/
aas./dnamic- mode ctw.ressiiiiiitoooooooiiiiiillll www.nnnaeeaa---a and 7M477688825511114421116666668113
-naeeriiaalaaddeeeeeeeemaaallll experinental
[22] Jonathan H. Tu, Clarence W. Rowley, Dirk M. Luchtenburg, Steven L. Brunton, and J. Nathan
Kutz. On dynamic mode decomposition: Theory and applications. Journal of Computa- lional Dynamics, 1(2):3914211, December 2014. ISSN 2198- doi: 10.3934/jcd. 201 4. 1.391. URL https: /wwww aimsciences. otgg/en/article//o///oo 3934/ jed. 2014.1 Publisher: Journal of Computational Dynamics.
391
12

---

[23] Steven L. Brunton, Bingni W. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Koop- man Invariant Subspaces and Finite Linear Representations of Nonlinear Dynamical Sys- tems for Control. PLOO ONE, 11(2))e01 50171, February 2016. ISSN 1932-6203. doi:
10.1377/ouumaal.onne 015017.. URL https ///a:://jourall plos..los....o org/ wttpplosooe//aeeeeoooooooooooooe aatticel Publisher: Public Library of Science.
id=10. (d=10..371//5 1171/ pone 0150171
[24] Weinan E and Bing Yu. The Deep Ritz Method: A Deep Learning-Base Numerical Algorithm for Solving Variational Problems. Communications in Mathematics and Statistics, 6(1):1-12. March 2018. ISSN 2194-671X. doi: 10.10077440304-018-0017---- URL https:///ooii org/
10. 10.1007//40304-01 018-00000000-000000000 -0000/8800004-11-1112--
[25] Justin Sirignano and Konstantinos Spiliopoulos. DGM: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 375:1339 1364, December 2018. ISSN 0021-9991. doi: 10. 10.11166/..jo jep. 2018.001 com/science//ani htt.aa/eeeppp////000009911133330
/tnttps:/wwwwr.e tttp:/www.ccc....c.co www.sccciieeeee.....
[26] M. Raissi, P. Perdikaris, and G. E. Kariiadakis. Physics
Physics-informed informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics. 778:686-707, February 2019. ISSN 0021-9991, doi: 10.1016/j ptt.:///.jjppp.......... 2018.10.45..
URL https: wwwp://wwwwwaaao/ www. www.nceiieeccccccccmm... com/
science/ /cc/ene//arrccee/////00099999919
[27] Ilan Naiman and Omri Azencot. Operator Theoretic Approach for Analyzing Sequence Neural Networks. Proceedings of the AAAI Conference on Artificial Intelligence, 37(8))9268 9276, June 2023. ISSN 2374-3468, doi: 10. 1l.16609aai.v31177 v37.i...77....11.. 2.11.. https://///: aaai org/index, php/ /tt.////aatttteeviiieeeeeee 26111
Number: 8.
[28] Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature Communications, 91)):4950, November 2018.
ISSN 2041-1723. doi: 10. 00888/44166-011--77100 articles/s41467- -1110000000000000000000--0000000
URL https: //www. www.rentuee...m com/ Publisher: Nature Publishing Group.
[29] Enoch Yeung. Soumya Kundu, and Nathan Hodas. Learning Deep Neural Network Represen- tations for Koopman Operators of Nonlinear Dynamical Systems. In 2019 American Control Conference (ACC), pages 4832-4839, July 2019, doi:
/0/19/.2009. ISSN: 2378-5861.
URL
httpss:/jeeoexplorss ieee org/ www.oueeet//11555000
[30] Shankar A. Deka, Alonso M. Valle, and Claire J. Tomlin. Koopman-based Neural Lyapunov functions for general attractors. 2022 IEEE 61st Conference on Decision and Control (CDC), pages
December 2022. doi: 10.1 (0.1109/CDCC51059 (02.1105...2229999 2022.9992927. URL
ieee. wtt.ee/e..r//.ssttttta/oocaetete
ISSN: 2576-2370.
[31] G. Haller and G. Yuan. Lagrangian coherent structures and mixing in two-dimensional turbu- lence. Physica D: Nonlinear Phenomena, 147(3):352-370, December 2000. ISSN 0167-2789. doi: 10.10166501167-77799000000000110
http::/::://wwe http://wwww.ciincccirettt.tt wtw.e/tttteccteeeoeooooeeeeeeeee
arricile/pp1//0167777990000//
[32] Martin L. Tanaka and Shane D. Ross. Separatrices and basins of stability from time se- ries data: an application to biodynamics. Nonlinear Dynamics, 58(1):1-21, October 2009.
ISSN 1573-269X. doi:
URL https:://ai.oorrg//////
81 811071-008-9457-
[33] George Haller. Lagrangian Coherent Structures. Annual Review of Fluid Mechanics, 47(Volume 47, 13751:137-162, January 2015. ISSN 0066-4189, 1545-4479. doi: 10.
1 1146/annurv- fluid- 1ate--110111---- 141322. URL https: /www. annualreviess org/contente Publisher: Annual Reviews.
j journals/10. http://aanrooo--idd1oo1-00-----0
141322
[34] A. Mauroy, Mezić, and J. Moehlis. Isostables, isochrons, and Koopman spectrum for the action-angle representation of stable fixed point dynamics. Physica D: Nonlinear Phenomena, 261:19-30, October 2013. ISSN 0167-2789. doi: 10. po
URL fhttps
13

---

[35] Igor Mezić. On applications of the spectral theory of the Koopman operator in dynam- ical systems and control theory. 2015 54th IEEE Conference on Decision and Con- trol (CDC), pages 7034-7041, December 2015. doi: 10.1 1011109CCDCC221154433333 URL
htttps://eeexploree ieee org//oocmmett7746339
[36] Igor Mezić. Koopman Operator, Geometry. and Learning of Dynamical Systems. Notices of the American Mathematical Society, 68(07):1, August 2021. ISSN 0002-9920, 1088-9477. doi: 10. 10.1099/noti2306. URL http:: ams http/////o....c////////////oooo/
pdf
[37] Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio. Empirical Evaluation of Gated Recurrent Neural Networks on Sequence Modeling, December 2014.
URL http:
7/ 7/arxiv. org/ Berrg/bbs/4415 /412. 3555 arXiv:1412.3555 [cs].
[38] Richard R. Stein, Vanni Bucci, Nora Toussaint, Charlie G. Buffie, Gunnar Rätsch, Eric G. Pamer, Chris Sander, and João B. Xavier. Ecological Modeling from Time-Series Infer- ence:
Insight into Dynamics and Stability of Intestinal Microbiota. PLOS Computational Biology, (12))e 1003388, December 2013. ISSN 1553-7358. doi: 10. 10.1371/journal.ppp
1003388. URL
plos Publisher: Public Library of Science.
journal pcbi..pi...... 1003388
[39] Eric W. Jones and Jean M. Carlson. Steady state reduction of generalized Lotka- Volterra systems in the microbiome. Physical Review E, 99(3)x032403,, March 2019. doi: 10.1 10.3/hyyysseeE. 032403. URL https: //1////inn..////////// aps wtt.aa..p...oggooo/.o.o......... http://////....1//////eeeeeee 99. 032403 American Physical Society.
Publisher:
[40] M. D. Powell. Radial basis functions for multivariable interpolation: a review. In Algorithms for approximation, pages 143-167. Clarendon Press, USA, January 1987. ISBN 978-0-19- 853612-3.
[41] Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization, January 2017. URL http://arxiv org/ www.////....c... 6980 11Xiv:1412.6998 [cs].
[42] Ido Cohen and Guy Gilboa. Latent Modes of Nonlinear Flows: A Koopman Theory
Analysis. Elements in Non-local Data Interactions: Foundations and Applications, May 2023. doi: 10. 10.1017/9781009922866
URL //bsttp::/ww/// /wwww/wwwwww////////////.//. cambridge wwwbridge.oog/coooooooooooeooeee org/ htt///oreeeeeeeeeeeeeee elements/ /ttef//00///////0000000110000000
Latent-modest of matent-modes--f-nnoni 9781009323826 9781009323857
ISBN:
Publisher: Cambridge University Press.
[43] Igor Mezić. Spectral Properties of Dynamical Systems, Model Reduction and Decom-
positions. Nonlinear Dynamics, 41(1):309 325, August 2005. ISSN 1573-269X. doi:
10. 10.1007/111071-005--28244
/tL/htttss://ooi.rrr/// 1000000000010011-00-22222/
[44] Igor Mezic and Amit Surana. Koopman Mode Decomposition for Periodic/Quasi-peiinn Time Dependence*. PTAC-PaperSOnLine, 49(18))(1)))) 49618))690-677, January 2016. ISSN 2405-8963, doi: 10.
10.2244 URL httpt:://ww..cciieeeeeeee http://////cc..ceeecccceeeecccc/ www. http:///wwrrrtroo//c////////eraa p11/524058963663182222
[45] Yiming Meng, Ruikun Zhou, Melkior Omik, and Jun Liu. Koopman-Based Learning of Infinitesimal Generators without Operator Logarithm. In 2024 IEEE 63rd Conference on Decision and Control (CDC), pages 8302-8307, December 2024. doi: 10. 10//CD557744 2024. 108860848 URL htllltttps//ieeepplooo ieee. orgg/oocuent/1086600 ISSN: 2576- 2370.
[46] Igor Mezić and Andrzej Banaszuk. Comparison of systems with complex behavior. Physica D: Nonlinear Phenomena, 197(1):101-133 October 2004. ISSN 0167-2789. doi: 10.1016/
ipphydd.2004.6..00 URL htttt::t///// htttss//www.scceecciieeettc..c wwwencedirect.comeee httpr/dddeeec.cooccccc////rrtaat
$01672789040005077
[47] Shankar Deka and Dimos Dimarogonas. Supervised Learning of Lyapunov Functions Using Laplace Averages of Approximate Koopman Figenfunctions.. IEEE Control Systems Letters, 7:3072-3077, 2023. ISSN 2475-1456. doi: 10.10099 URL.
https: fttps::/ieeexploo// weee http///aabaaatt/oooooeetttcce
14

---

[48] Garrett E. Katz and James A. Reggia. Using Directional Fibers to Locate Fixed Points of Recurrent Neural Networks. IEEE Transactions on Neural Networks and Learning Systems, 29 (8):3636 August 2018. ISSN 2162-2388. doi: 10. AAT10//TNNN..2211227755555
URL
htttps:/Aeeeexploree ieee. httpeeerrrr/occmmmt///o//3222
[49] Elia Turnrr Kabir V Dabholkar, and Omri Barak. Charting and Navigating the Space of Solutions for Recurrent Neural Networks. In Advances in Neu- ral Information Processing Systems, volume 34, pages 25320-25333 Curran Asso-
ciates, Inc.. 2021. URL https::/ppoceeeding: www.gg...urri.... pttpr/pp..////p//////////aaaaa wttpp//ppppp//////aallll 453445447719994427eebbbbba htg.:/oow....abaaaaaaeeeaoo....o Abstract.htmi
[50]
Connor Brennan, Adeeti Aggarwal, Rui Pei, David Sussillo, and Alex Proekt. One dimen- sional approximations of neuronal dynamics reveal computational strategy. PLOS Computa- tional Biology. 1911)ce111784, January 2023. ISSN 1553-7358. doi: 10. /0.13a// 1010784. URL https://journalss. plos. http://oooccapoo/oaarreeeeeeeeee
journal pcbi.c................. 1010784
Publisher: Public Library of Science.
[51] Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian Learning and Inference Recurrent Switching Linear Dynamical Sys- tems. In Proceedings of the 20th International Conference on Artificial Inielligence and Statistics, pages 914-922. PMLR, April 2017. URL hatpps:/pproceedings mlr, press/
154/1inderman11 html
ISSN: 2640-3498.
[52] Lea Duncker, Gergo Bohner, Julien Boussard, and Maneesh Sahani. Learning interpretable continuous-time models latent stochastic dynamical systems. In Proceedings of the 36th International Conference on Machine Learning, pages 1726-1734. PMLR, May 2019. URL
https://proceeddngg. mlr. wwwps/ vttess/v97/duncker111 duuckknnnkkerre 19a. html ISSN: 2640-3498.
[53] Amber Hu, David Zoltowski, Aditya Nair, David Anderson, Lea Duncker, and Scott Linder- man. Modeling Latent Neural Dynamics with Gaussian Process Switching Linear Dynamical Systems. Advances in Neural Information Processing Systems, 37:33805-33833 Decem- ber 2024. URL https: httt:ss:/ppocceeedngge nww.ips. wwwpcpppppafiiiieeeeee files wttee/feeeeeeppeee/e///2 pttpp/pppp//aaaaa///ee 9b64416991502666444bb 04004444440110011040004444440-00
html
[54] Andrew Wagenmaker, Lu Mi, Marton Rozsa, Matthew S. Bull, Karel Svoboda, Kayvon Daie, Matthew D. Golub, and Kevin Jamieson. Active learning of neu- ral population dynamics using two-photon holographic optogenetics. Advances in Neural Information Processing Systems, 37:31659-31687. URL
December 2024.
https://proceeddingg neurips wttpccpppppp...ee http/pppr_ffless__ wateees/apprr/22/2 (al.efee.ood16666055577770000000
-wwt--bttt--ooo Conferenceeeeeeooo html
[55] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep Residual Learning for Image Recognition. In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770-778, Las Vegas, NV, USA, June 2016. IEEE. ISBN 978-1-4673-8851--1 10.1109/CVPR.2016
doi:
URL http::///eeeexpporeeeeeee www.:iieeee...eeeooooooee wtw.eeeeeeeeee...oeeoeee///////e org//dccuent//780055)
[56] Yulong Lu, Jianfeng Lu, and Min Wang. Priori Generalization Analysis of the Deep Ritz Method for Solving High Dimensional Elliptic Partial Differential Equations. In Proceedings of Thirty Fourth Conference on Learning Theory, pages 3196 3241. PMLR, July 2021. URL
https://proceedingg: mlr http://dnnnn.nn.ppprrrss//////// html
ISSN: 2640-3498.
[57] Jiequn Han, Arnulf Jentzen, and Weinan E. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115(34): $505-85100 August 2018. doi: 10. 10.1073/pnas..m 17189422111911115. URL htRh:ttss:///w .nwa org/
doi/ful1/10. 107.0pna....... 17 1718942115
[58] Zheyuan Hu, Khemraj Shukla, George Em Karniadakis, and Kenji Kawaguchi. Tackling the curse of dimensionality with physics-informed neural networks. Neural Networks, 176: 106369, August 2024. ISSN 0893-6080, doi: 10. 10.11066j.neume22224...666 URL hhtpps
15

---

A Analytical KEF derivation in 1D bistable system
We would like to find an analytical Koopman eigenfunction for the scalar dynamical system:
(14)
In the 1D case, the Koopman PDE 3 reduces to first-order ordinary differential equation
dv
dx f(x) (xb(x).
(15)
With f(x) = - -x and X 1 we have:
v' x2 = vx(x))
(16)
(z)z) 1
(17)
(x)
2 3
To solve this integral we first simplify the integrand.
r - 3 = xx1 - x²) = - x2)++))
(18)
So,
A
B
C
(19)
x(1 x1 +
x
1 x
Solving for A, B. C yields A 1, = 1.
Now we can integrate,
/
1
13 dx =
+
dr
(20)
2(1 x)
2(1 + x.)
1
=log
|z|
log
x
log + +
(21)
1
(x)) = log [x|
log
|1
x
log 1 + + + C
(22)
2
[=]
(1x) C².
(23)
1 x
To bring it into the form in the main text our solution by the sign(r) function which is a À = 0 eigenfunction because it remains constant in each basin (see Figure 4A). In other words, we flip the sign of our solution (xx) + (2) for < 0.
the product composition rule 99) We can multiply
((x) = C'
(24)
VII - 1
and this remains a KEF with À = 1.
B Eigenfunction Degeneracy in higher dimensions
Consider a separable 2D dynamical system:
1 = fr(x,
(25)
ý = fv),
(26)
16

---

which we write compactly as:
[f((z)
X = x=f(x,)) v) -
[v(()
(27)
We seek a Koopman eigenfunction (ax y) satisfying: Vih- f(x, v) (1p
(28)
Assume and a separable form ((r,)) = (v) Then:
O
dr (z)))()))
(29)
0 (())Y'(w)
(30)
dy
VΦ. f = XXXXXXX
(31)
Dividing both sides by (v)) gives:
(z() Y'()) ((x) fo(
fo()) =1
(32)
Y(y)
The above equation requires that the sum of the above two terms, which each depend on different variables must be for all z, y. follows that each term is also constant function.
X'(x)
X(x) f=(x) u4
(33)
Y'(y) Y(y) fy(v) -1- M.
(34)
for an arbitrary constant j E R.
Define the antiderivatives:
A(x) Fa(z) da,
(35)
B(y) / fu(e) dy.
(36)
Then the logarithms the separated components are:
X(x) ( (e)(r)) "
log X (x) = pA(x)
(37)
log Y (y) = (1 - Byy) Y(y) (eBlv)) 1-#
(38)
Thus, the general separable Koopman eigenfunction is:
(1x, y) (edlc))) ((B(w)) 1-p
(39)
C Relation of our definition to the Koopman Operator
In the main text, we introduced Koopman eigenfunctions scalar functions 3 X that evolve exponentially along trajectories s=t)) E X' of a dynamical system á = f(xx):
(x(()) = (t(((tt)).
(40)
dt
Here, we clarify the origin of this equation by defining the Koopman operator, linking our approach to the broader theory.
17

---

Let : X' - I be a real-valued function of the system state commonly referred to an observable. The collection of such observables forms an nnfnite-dimensional function space, typically a Hilbert space once equipped with an inner product (9.9") g") The Koopman operator acts linearly on this space.
Remark. When studying nonlinear dynamical systems with multiple basins of attraction, as we do, the corresponding Koopman eigenfunctions are generally not square-integrable, and therefore fall outside the Hilbert space defined by the standard inner product. We aware of this theoretical limitation and continue to employ the Koopman framework regardless. In practice, neural networks learn finite, smooth approximations to these otherwise singular structures.
For continuous-time system, the Koopman operator K,, evolves observables according to the flow map F+ : X' + X' which advances the state forward by time 7:
(K+1)()((()) (=(()())) g g( + r)).
(41)
The infinitesimal generator of the Koopman semigroup )-20 often denoted simply as K, is defined as:
(F)(()) g(x)
Kgg lim Krg
lim
(42)
When evaluated along trajectory xtt), this yields:
g g(x(t + r)) g(x(tt)
Kg(x((t) = lim
(43)
7
d
((x(()) Vo(x((t) (tt) = Vg(((()) s(e(t))
(44)
dt
This operator is also known as the Lie derivative of g along the vector field f.
Thus, an eigenfunction ub of satisfying
K Xa/
(45)
recovers the Koopman cigenfunction equation
3 used in the main text.
D KEF degeneracy in randomly initialised DNN solutions
Main text Figure illustratess challenges arising due the degeneracy of the Koopman PDE (3). In
Figure we train several DNNs 2-unit GRU trained the 2BFF. Each DNN is independently initialised and trained on a single distribution without the balance regularisation term L i.e.,
That = O. The resulting KEF approximations exhibit the same modes of degeneracy zero on certain basins as well vertical and horizontal variants.
E Curve-based validation approach
In high-dimension, we cannot visualize the entire phase space to check whether zeros of the KEF coincide with the separatrix. Instead, we generate a family of smooth curves that connect two attractors, and hence must pass through separatrix. In the 64D GRU flip flop example, the two attractors are the two stable fixed points x, E R We use cubic Hermite interpolation with randomized tangent vectors at the endpoints. Each curve is defined by:
(aa) (ao x + (1o(a) m. + hoi (oi(a)) + h11(a) my, a a [0
(46)
where the Hermite basis functions are:
hoo(a) = 2a - 3a + 1, hoi (a) = -2a3 3 + 30.
(47)
(48)
hio(a) a 2a + a,
(49)
h(((a) = a o².
(50)
18

---

Model-5
Model-6
Model-7
Model-8
Model-9
Figure 9: Many KEFs of to for 2 bit flip flop in 2D
Notice that (0) and H = The tangent vectors my and my initialized as y and perturbed with Gaussian noise:
m2 = (y- x) Exs my = (y - x) + Exy Ez, N(0,
(51)
We sample multiple such curves with independently drawn perturbations. This produces a family of curves that interpolate between and y, while varying in geometry, enabling randomized exploration of intermediate regions state space. Crucially, the curves are not limited to the manifold spanned by the attractors, but extend to all dimensions (controlled by Optional constraints (e.g., non-negativity) can be imposed by rejecting any curve that violates them. For each such curve, we both evaluate the KEF and simulate the ODE to determine the position of the separatrix.
F Neural network architectures
In most of the demonstrations we use a ResNet architecture [55 with a tanh activation function.
Let the input to the network be Xin E Rdi,, and let the hidden activations be ( Rehia, with output Nout E Raa and layers.
The network receives inputs at the first layer (0) Pad(a where Pad appends zeros to the input (we always choose dnia > din), and E R2+ is scalar chosen so that the elements of x (() are (1). The network then updates the hidden state at each layer l as
x ((+1) ( + tanh [w(0 +b(O), e 0...... L-1,
(52)
where (t)) E Ru x and ( E IRd
The output is obtained by applying final linear layer:
Xout Wont x (4) + b out
(53)
where Wrout Ro xdhi and hout E Rdos, During optimization, gradients V CoCiotal are computed for all parameters 0 (0)., b (0:--1),
Wrutt but).
Choices for each system
For the results in Figures D
use L 20, dhid 400. dout and d,n N the dimension of the dynamical system. For the gLV system in Figure we use L. 25 and dhia 1000.
Radial Basis Function (RBF) Layer
For the limit cycles example in Figure 7 we single Radial Basis Function layer [4]]
19

---

Given an input X E Rd,, the RBF layer maps it to output E R through set of M radial basis functions, each centered at c, E R with a shape parameter E; and linear combination weights dij E IR.
To compute RBF activations for i ,.... M, define the scaled radial distance:
(x) Ei ||x call
(54)
and then apply a gaussian radial basis function
fi(x)
(55)
The final output is a linear combination of the basis activations:
M
xx)))
aj (x), j=
dont
(56)
We use .M = 300, and doat = 1. During optimization, gradients VOLual are computed for all parameters 0 (Lajs), {ea)).
G
Optimisation
Our optimisation procedure described in the main text summarised in Algorithm T imple- ments the training of the neural network Koopman eigenfunction using the ratio loss and balance regularisation, for multiple sampling distributions.
We minimise the total loss:
Crotal Σ Ciaioo + That Cd
(57)
where j corresponds to the 7"h sampling distribution (see main text section 3.2p. N-dimensional points in the state space are sampled from each distribution a! ~ (x)) The ratio loss is the Koopman PDE error, normalised by sample-shuflled version:
(LHS - RHS!)²
Caaio
(58)
(LHS - RHS
perm()) )
LHS/- = V(((a)) f(x))
left-hand-side of the Koopman PDE
(59)
RHS/ = {ap()))
nighthhand-sie of the Koopman PDE
(60)
where perm (i) is random permutation of the numbers ......., B sampled during each training iteration.
The balance regularisation loss is the squared mean of the KEF values divided by their variance:
(20)
L
(61)
B
(v(x²) - )²)
List
J = B (ax)).
(62)
In general 7hat 0.05. For the limit cycles Figure 7 we Tal 0.
We compute Vv(x) using Pytorch's torch. autograd gradd specifying create_ since we differentiate through this second time to compute the gradients VoLroual with respect to the neural network parameters
We use the Adam optimiser [41 with leamiigg rate 10-4 and 12 normalisation 10 We use B 1000 and train for 1000 iterations.
Only in the case the 11D gLV, Figure 6 we use B 5000 and train for 5000 iterations.
A summary of all hyperparameters is provided in Table I
20

---

Algorithm Train Koopman Eigenfunction Network
Require: Sampling distributions {z) vector field f; eigenvalue A; neural network architec-
ture tha batch size B; iterations T; balance weight "hali learning rate yi small 10-12
1: Initialize 0
2: for 1 T do
3: Liotal f- 0
4: for
J do
5:
{ar vh. V V } SAMPLEANDEVALUUATEO), Vo, B)
6:
L Littio COMPUTERATIOOLOSS((E. vy f(z))
7:
1i + COMPUTEALLNNGGLLOSSS(
Liotll Liotal + Limi + 7bal Lia 9: end for
8:
10: Compute gradients of Lotal w.r.t. 0
11: Update weights using gradients and learning rate
12: end for
13: return trained parameters 9
14: procedure SAMPLEANDE EALAAATE(ALL ch. B)
Sample {x: Ji Pj (x)
16: Compute v/r +- vho(a;))
17: Compute V V +- V Vntis (2x)
18: return {ar dh W
19: end procedure
20: procedure COMPUTERATIOLOSSS(E VV.
21:
Compute LHS/ = V f(x))
22:
Compute RHS/ A
23:
Draw random permutation of [1....,B, B}
24:
ny Σ, (LHS) RHS )2
d, + E,HS RHS )
return L = nd,/dd, +
27:
end procedure
28: procedure COMPUTIBALANCELLOSS((f
29:
1 E. w/ E.
30:
- 22
31: return Liai = - +5)
32:
end procedure
Table I: Algorithm details and hyperparameters for various systems. System dimensionality N, Koopman eigenvalue X, balance regularisation weight %aal batch-size B, training iterations T. leaming rate ResNet depth L and width daia, number of Radial Basis Functions M.
Dynamical System
N
B
T
L
duid
M
Bistable 1D
0.05
1000
1000
10-1
10
20
400
Damped Duffing oscillator
2
0.05
1000
1000
10
20
400
IBFF, 2D GRU
0.05
1000
1000
10-4
400
2BFF, 3D GRU
3
0.2
0.05
1000
1000
10-4
20
400
IBFF, 64D
64
0.1
0.05
1000
1000
10
20
400
Two Limit Cycles
2
0
1000
1000
10
300
Ecology gLV
11
0.1
0.05
5000
5000
10-4
10
25
1000
Data-trained RNN 9
668
0.02
0.05
1000
1000
10
7
1200
21

---

G.1 Choice of Training Distributions
Choosing suitable training distributions (z) an important step in applying our method. The distributions are selected to satisfy the following criteria:
1. They are approximately bisected by the separatrix, ensuring roughly equal sampling from both basins and enabling the balance loss to be satisfied.
2. They sample sufficiently near the attractors capture the global bistable dynamics.
3. They approximately forward-invarian, i.e., trajectories initialized from the distribution remain within its support when evolved forward time. This avoids loss of mass due to transient amplification along unstable directions.
4. In systems with multiple spurious attractors, they avoid sampling from basins outside the domain of interest.
Below we list the specific distributions used for each system.
ID bistable system (Fig. 3C): (V 1), (),
0.5, 1.0, 2.0), where I2 the 2D damped Duffing oscillator (Fig. WV(,, o2 with oy 0.5,
2 x 2 identity matrix.
2D 1-hit flip-flop GRU (Fig.
3 (O, o212), with oj = [0.01, 0.1, 0.5, 1.0,
3D GRU 2-bit flip-flop (Fig. 5) N Is). with o; |0.01, 0.05, 0.2, 1.0, 5.0). The mean μ is
chosen as a point on the separatrix found by interpolating between two attractors and performing iterative binary search with ODE (T. A second μ obtained by interpolating a different attractor pair for training the second KEF.
11D ecological dynamics (Fig. :o: We similarly identify single separatrix point μ. Each coordinate di is sampled independently from a Gamma distribution æ[i] [[a[ii, B[i)), where a(i] and 3[i] are chosen such that the mode of æ[i] equals μ[ii] and the variance equals oj. with oy [0.01, 0.1, 1.0)
Data-trained RNN 19 (Fig, 8 N( 14, oyE). where we first construct a distribution which is
oblongated along the direction of the attractors and isotropic along the remaining directions:
Σ = o IN + (o - -0 uu
with
u E R unit vector along the attractor axis, i.e. u (a (a- ba-b)/lla bll. where a and are the two attractors,
> 0 the standard deviation along u.
d > the standard deviation along all orthogonal directions, and
Ly the N x N identity matrix.
As before μ is the point on the separatrix along the line joining the attractors.
We set σ to include both attractors, and choose as large as possible while avoiding spurious basins. For 300 samples drawn from from (H, E), we evolve the dynamics forward for 3.0 and estimate the covariance of the resulting approximately
distribution, denoted Σ.
H Scaling of compute with Dimensionality
We ran all experiments on a system with four GeForce GTX 1080 GPUs with 10 Gbps of memory each.
All the 2D systems take 1-5 minutes to train the KEFs. The 1ID gLV takes up to 20 minutes. The 668D data-trained RNN 9 takes 5 minutes train the KEF
Scalability to high-dimensions key strength of our approach. We evaluated the scaling behavior of our method vanilla RNNs of varying sizes, each trained on the 1-bit flip-flop task. In table 2
22

---

we report the wall-clock training time and curve-based based validation performance for learning a single Koopman eigenfunction:
Table Training time and performance across dimensionalities. Vanilla RNNs with N units were trained on the 1-bit flip-flop task. We then applied our method to each resulting dynamical system f to approximate its corresponding Koopman eigenfunction.
Dimension N
Wall-clock time (s)
Curve R
32
529
0.997
64
527
0.995
128
572
0.867
256
611
0.996
512
718
0.996
For each case, we verified good agreement with the ground truth using the curve-based validation metric (Fig. BA-C). CC. All models used fixed DNN architecture with depth and width 550. For larger N, we expect that the network width must scale with dimensionality. Since our method involves solving an N-dimensional PDE, its scaling behavior is expected to resemble that of Physics-Infomin Informed Neural Networks (PINNs) (26] and other PDE-solving neural methods (24] In particular, Lu et al.
[56] proved N N-independent generalization error bounds for the Deep Ritz Method within a class of PDEs. Addressing the curse of dimensionality remains an active area of research [57, 58 and we
anticipate that advances from this literature can integrated into our framework.
I Choice of eigenvalue for numerics
In the main text we look for approximations to the Koopman PDE 3 for real positive eigenvalue X. What should the value of be? It is known that products of KEFs are KEFs themselves with different eigenvalues. In particular, for a KEF & with eigenvalue A, we see that:
v [+(x)) anb)) f(x) = =x)
(63)
(64)
Therefore, (xx)" is also a Koopman cigenfunction, with eigenvalue at. This translates to changes in the shape of the KEF, i.e., the sharpness of the peaks, while maintaining the position of the zeroes.
In practice the choice of À affects training convergence, and it therefore an important hyperparameter in the optimisation procedure (see Figure TO). We attribute this to the time scale of interest in the system å f(x), and differences in the propagation of gradients for different A.
J Linking Koopman Eigenfunctions and separatrices: Formal Derivations
Setting. Let (X, d) be a smooth manifold with metric and flow ($*) ( generated by the autonomous ODE à: f(x).. Thus Φ0 = id, +i+= Th oh, and (1, æ)) ('(x)) is continuous.
Invariant set. An invariant set of a flow d is subset S C X such that
nh ('x)) S for all S for all e R.
Attracting set. A nonempty closed set A X is an attracting set if it is invariant and there exists an open neighbourhood U A such that
lim dist ( (2)))))))))))))))))) = 0 for all æ E U, t000
where for y X and A C X,
dist A) == inf d(y, a). aeA
23

---

Final Training Error Eigenvalue
Final Training Error vs Eigenvalue
10a
Final Loog/ Total Erorr
Final Loss Tocal Error
10%
10-1
10
109
10a
10 2
10
10
Eigenvalue
Eigenvalue
Figure 10: Training convergence as a function of eigenvalue X, evaluated by normalised PDE error fratio for systems: LEFT, ID bistable system J x x (see Figure 3) and 2BFF GRU 3D) (see Figure 5).
An attracting set A is an attractor if there is no proper subset of A that is also an attracting
Attractor:
Basins of attraction.
The basin of an attractor A is
((A) (ar E X lim dist( A) -}.
Each B(A) is forward invariant and open.
Separatrix. Let [AsheR be the set of all the attractors of oe in X. Define the separatrix as the complement of all basins:
= X U B(As)
AEK
Koopman eigenfunction.
Let u/ X (U= + R be continuous and satisfy
(((((x)) = ex (xx) Vx E X Uk As, 14 2 0,
for some eigenvalue A > 0. Then the sets &v 0}. { 0}. and { = 0} are forward invariant.
Constant sign near attractor (CS). Let v X' (U=A2) + R be continuous function. We say has constant sign near an attractor Ag if there exists an open neighbourhood Un with A C Uk C B(Ax) such that v has constant sign Ug \ Ag-
Proposition (Sign change a zero separatrix point). Ler Lergh: X (UkAu) + be a continuous Koopman eigenfunction with positive eigenvalue À > 0:
Φ (x)) = (x)),,
t t
and have constant sign near all the attractors As- Suppose i E X') (Uk Ax) satisfies vb(x)) = and, for every E 0, the ball Ba(x,E) contains points y y with (by+) > and ( < Then E E.
Proof. Assume for contradiction that an E B(Ag+). Since B(As)) is open, there exists ep > such that B(x, C ((A=)).
By (CS), v has constant sign on Up- \ Ag. for an open neighborhood Up. By attractivity, for any E B0), the trajectory &t (yy) eventually enters U2- and remains there for all large t.
Since (y)) = e (y)) the sign of v'y)) is preserved along the trajectory (forward invariance of (1b > 0} and (v/ < 0}), so ((r) must already share that constant sign. Thus all points in B(x, En) have the same sign, contradicting the assumption sign change in every neighbourhood. Therefore D
U, ((A)), i.c. E E.
24

---

Remark. The condition (b(x) = 0 with a sign change in every neighbourhood means that æ: is a boundary point of { 0} and { < 0}. Proposition shows such points lie precisely on the separatrix.
25

---

NeurIPS Paper Checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: We develop a numerical method to locate separatrices multistable dynamical systems. We utilise the fact that Koopman eigenfunctions identify separatrices and demonstrating with numerical experiments on several dynamical systems, that we find their separatrices.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made 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.
It is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Limitations mentioned in the discussion. Mainly lack of theoretical guarantees.
Guidelines:
The answer NA means that the paper has 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 separate "Limitations" section in 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 on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or 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 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 in the paper. The authors should use their best 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.
3. Theory assumptions and prool's
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
26

---

Answer: [NA]
Justification: We don't 't have any rigorous proofs, just simple examples to give an intuition of the behavior of Koopman cigenfunctions.
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 or referenced 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 provide short proof sketch to 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
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Code is provided in supplementary zip file. All the systems studied are contained in the code. Architectures and hyperparameters are provided in the Appendix
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 novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it 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 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, other means that appropriate to the research performed.
also be provided via detailed
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 it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it 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 reproducing or verifying the results.
5. Open access to data and code
27

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: Code attached in supplementary zip file. No datasets used in this work. Parameters of the 11D system included in the code. Code is also provided train the RNNs used for reverse-engincering.
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" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//apps.cc/pullic/uiiiee htt///CooSSSeiiaaeeeeeeee sisssssssiniyyyyyy for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data,
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments 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 as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permirted.
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: Yes. See Appendix G
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, 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: Many of our validations in 2D systems are qualitative. We provide in quantita- tive validation in Figure where we report the H²² coefficient of determination between predicted and true transition points.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
28

---

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.)
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 one 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 show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
error bars are reported in tables or plots, The authors should explain 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 execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: See Appendix
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 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 ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: /ate:://nuuipppe htt.///uulccttthiicccccclelliii.
Answer: [Yes]
Justification:
We reviewed the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS 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 a special consid- cration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: The paper proposes a method developed for understanding RNN models the brain and other dynamical systems. has no clear socictal impact.
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 or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness 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 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 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 or 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 attacks, mechanisms for monitoring misuse, mechanisms to monitor how 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:
No such risks.
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 to usage guidelines or 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 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 the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We cite all the code and models from other works.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URI.
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.
30

---

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.paperssit
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it 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: [Yes]
Justification:
We provide commentary in the code ensuring it's readability.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl 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: no crowd sourcing human subjects was involved.
Guidelines:
The answer NA means that the paper does 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: the paper did not involve crowdsourcing nor research with hum subjects.
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.
31

---

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: no non-standard involvement of LLMs
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 pettts://neurre http///arrpsscccooofeeeeeee////s cc//oonferences/22255/LI. for what should or should be described.
32


Reviews:
[{"review_id": "MdnsQy6rNx", "reviewer": "Reviewer_FEbm", "summary": "This paper presents an innovative deep learning framework for identifying separatrices in high-dimensional dynamical systems using Koopman eigenfunctions (KEFs). The authors develop a neural network-based approach to approximate KEFs with positive eigenvalues, which naturally vanish on separatrices, enabling efficient boundary detection through gradient-based optimization. The method demonstrates strong empirical performance across diverse systems including synthetic benchmarks, ecological networks, and recurrent neural networks, while also providing practical applications for designing optimal perturbations in neuroscience contexts.", "questions": "How does the method's performance scale with increasing system dimensionality, particularly beyond the 64D case demonstrated?\n\nCould you elaborate on the theoretical connection between the zero-level sets of KEFs and true separatrices in more complex, non-analytical systems?\n\nWhat modifications would be needed to apply this approach to stochastic dynamical systems where separatrices may not be sharply defined?\n\nHave you explored alternative neural architectures beyond RBF networks that might improve approximation accuracy or training efficiency?", "limitations": "Yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The methodology shows promising results across diverse systems (from 1D synthetic models to 64D RNNs), its theoretical foundations remain underdeveloped, lacking rigorous guarantees linking KEF approximation accuracy to separatrix identification. The approach creatively addresses degeneracy challenges through balance regularization and ensemble training, as illustrated in Figure 3, but computational scalability for very high-dimensional systems may be limited by the iterative PDE-solving process. The work makes significant contributions to dynamical systems analysis by bridging Koopman theory with practical neural network applications, though broader impacts on real-world neuroscience interventions require further validation.", "was_revised": true, "final_justification": "The method offers a fresh perspective on separatrix learning, but the lack of rigorous guarantees in high dimensions and stochastic environments, along with scalability concerns, restricts its broader applicability.", "rebuttal": "# Response to Reviewer FEbm\n\nWe thank the reviewer for the thoughtful and constructive comments. Below, we address the main concerns and questions raised.\n\n---\n\n**Broader impacts on neuroscience interventions**\n\nWe agree that experimental validation is ultimately required to assess real-world impact. However, as we discussed in our response to Reviewer dWGz, our method can be applied to *data-derived models*, including those learned from neural data. This opens the door to *in silico* experimentation and the design of optimal perturbations, which can then be validated through targeted neuroscience experiments.\n\n---\n\n**Scaling with dimensionality**\n\nThank you for raising this point. Scalability to high dimensions is an important strength of our work. We tested scalability on the 1-bit flip-flop system (with $N = 64$ in the main text, Figure 7) in increasingly higher dimensions. On a single NVIDIA GeForce GTX 1080 Ti GPU, we measured the runtime for training a single Koopman eigenfunction:\n\n- $N = 32$ → 529 seconds  \n- $N = 64$ → 527 seconds  \n- $N = 128$ → 572 seconds  \n- $N = 256$ → 611 seconds  \n- $N = 512$ → 718 seconds  \n\nFor each, we verified a good match to ground-truth using the curve based validation (Fig 7A,B,C). We used a fixed DNN with depth = 20 and width = 550. For larger dimensions, we expect that width will need to scale with $N$. Since our method involves solving an $N$-dimensional PDE, its scaling behavior is similar to that of Physics-Informed Neural Networks (PINNs) and Deep Ritz methods, and addressing the curse of dimensionality is an active area of research [2,3]. We anticipate that advances from this literature can be integrated into our framework.\n\n---\n\n**Theoretical connection between KEFs and separatrices**\n\nPrior work [5] has shown that for linear and globally stable nonlinear systems, there exist Koopman eigenfunctions whose level sets align with stable manifolds (i.e., separatrices). Our setting is more general: we make no assumptions about fixed points or global stability and learn eigenfunctions purely from the dynamics.\n\nAs such, there is no strict one-to-one correspondence between the zero level set of a learned KEF and a true separatrix. In fact, we illustrate this explicitly in Figure 3:\n1. A zero level set not aligned with a separatrix (top row, right)\n2. A separatrix without a vanishing KEF (bottom row, $\\mu = 0$ and $\\mu = 1$)\n\nThese examples are not pathologies, but expected outcomes under our relaxed assumptions. We analyze and address them using our degeneracy mitigation strategies (Figure 3, Section 3.2).\n\n---\n\n**Application to stochastic dynamical systems**\n\nFor stochastic systems of the form  $$\ndX = f(X)dt + \\sigma(X)dW,\n$$  a common approach in the literature is to analyze the deterministic flow $f(X)$ by ignoring the noise term ($\\sigma\\equiv0$). This is frequently done in phase-plane and stability analysis [1]. Our method can be applied in this manner, with no modifications needed.\n\nIn learned models, uncertainty in $f(X)$ due to limited data or inference procedures can be propagated to quantify uncertainty in the separatrix location [5]. However, this is beyond the scope of the current work.\n\nWe also consider another class of stochastic models commonly applied to neural data: Recurrent Switching Linear Dynamics Systems (rSLDS)[6]. These consist of continuous and discrete stochastic variables coupled to each other. This yields regimes of the continuous phase space that are governed by different linear dynamics, probabilistically. To apply our method to such a model we first construct a smooth, deterministic approximation to the dynamics (details in appendix below), and then apply our method to it to obtain the KEF. We tested this idea on a 2D rSLDS with a 5-state discrete variable, trained on 2D synthetic data from a stochastic bistable system, then trained a KEF on the smoothed dynamics of the rSLDS and found the separatrix reliably. We would be happy to include our results in the final version to demonstrate the applicability to stochastic models and models derived from data.\n\n\n---\n\n**Alternative architectures**\n\nThe Radial Basis Function (RBF) network was used only for the 2D limit-cycle system (Figure 6). All other experiments used deep ResNets with tanh activation, as described in Appendix F. We apologize for the lack of clarity on this point in the main text and will revise accordingly.\n\nOur (unreported) experiments showed that ResNets outperform vanilla MLPs of comparable size in terms of stability and approximation accuracy. We also found tanh activations to be more stable than ReLU in our context.\n\n---\n\n**On theoretical rigor**\n\nWe acknowledge the absence of formal guarantees linking KEF approximation error to separatrix accuracy. While our empirical results are strong and supported by theoretical intuition (e.g., zero-level sets of certain KEFs align with separatrices), establishing tighter theoretical bounds remains an important direction for future work.\n\n---\n\n**References**\n\n[1] Duncker et al., *Learning Interpretable Continuous-Time Models of Latent Stochastic Dynamical Systems*, ICML 2019  \n[2] Hu et al., *Tackling the Curse of Dimensionality with Physics-Informed Neural Networks*, Neural Networks 2024  \n[3] Lu et al., *A Priori Generalization Analysis of the Deep Ritz Method*, COLT 2021  \n[4] Hu et al., *Modeling Latent Neural Dynamics with GP-SLDS*, NeurIPS 2024  \n[5] Mezić, I., *Spectrum of the Koopman Operator, Spectral Expansions in Functional Spaces, and State-Space Geometry*, Journal of Nonlinear Science, 30(5), 2091–2145 (2020)\n[6] Linderman et al., Bayesian learning and inference in recurrent switching linear dynamical systems, AISTATS 2017\n\n---\n**Appendix: Non-linear deterministic approximation of an rSLDS**\n\nTo derive a smooth deterministic flow approximation of the rSLDS dynamics we consider the posterior $p(z|x)$ relating the discrete variable $z\\in\\{1,2,...,K\\}$ -- which chooses the recurrent dynamics matrix $A_z\\in \\mathbb R^{N\\times N}$ -- and the continuous variable $x\\in \\mathbb R^N$, whose dynamics are governed by $A_z$. One possibility for the flow is to choose the most-likely dynamics at each point,\n$$\nf(x) := A_{\\tilde z(x)} x \\quad \\text{such that} \\quad \\tilde z(x) = \\underset{z}{\\text{arg max}} \\ p(z \\mid x).\n$$ This yields a piecewise linear function that may suffer from discontinuities. Another option is to take the 'expected' dynamics:\n$$f(x):=\\langle A_z x\\rangle_{z|x}=\\sum_z p(z|x)A_{z}x.$$ This yields a smooth deterministic non-linear dynamical system. We can also choose an intermediate option by defining an inverse-temperature parameter $\\beta$ and applying a softmax:\n$$p(z|x)\\leftarrow \\frac{e^{\\beta \\log p(z|x)}}{\\sum_{z'}e^{\\beta \\log p(z'|x)}},$$ where $\\beta\\to\\infty$ recovers the argmax version. In our numerical experiment we choose $\\beta=2$, maintaining a close match to the original piece-wise linear system but smoothing the transitions between linear regimes."}, {"review_id": "gvSUYXnyd4", "reviewer": "Reviewer_dWGz", "summary": "The paper presents a new approach to identify and characterize boundaries between basins of attraction (separatrices). It uses Koopman Eigenfunctions (KEFs) with positive eigenvalues that vanish near separatrices. The method locates points in state space where these positive KEFs equal zero to identify separatrices. The authors demonstrate the approach on multiple datasets, including synthetic data, ecological models, and recurrent neural networks. They also show the method’s potential by finding optimal perturbations to shift systems between basins.", "questions": "1) Regarding “encouraging sign changes across the separatrix” (line 136): I don’t understand how the cost in Eq. 10 forces a sign change. I see how it promotes non-zero values on both sides of the separatrix, but it seems possible for the eigenfunction to have the same sign on both sides while passing through zero at the separatrix. Could you please clarify this?\n\n2) About the steps around Eqs. 13 and 14: Do you first find individual separatrices and then combine them? In other words, do you solve Eq. 14 for each $\\psi$ separately and then combine the separators afterward? Why did you choose to present Eq. 13 before Eq. 14? I think reversing this order would improve the flow.\n\n\n3) In Eqs. 15–16, why do you assume Euclidean distance? Could the relevant effort be in a different metric space instead?\n\n4) Why did you choose to use a DNN to approximate the Koopman eigenfunctions instead of directly optimizing the loss functions (e.g., via EM or variational methods) without a network?\n\n5) It seems that for most examples, the authors included some ground truth knowledge that helped with evaluation. While I agree it makes sense to evaluate the method using examples where we have some idea of what the separatrices should look like, I wonder how would you evaluate the method in the absence of ground truth data? That is, are there any properties of either the separatrix or the values of the eigenfunctions around the separatrix that could help assess the certainty of the identified boundary when ground truth is unavailable?", "limitations": "Yes, also see in previous cell under Weaknesses", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 2, "originality": 2, "strengths_and_weaknesses": "I enjoyed reading the paper and believe it addresses important challenges in dynamical systems modeling.\n\n## Strengths:\n1) The paper is very well-written, with the motivation and approach clearly explained and presented in a way that is easy to follow.\n2) Figures are visually appealing and reflect the strengths of the approach.\n3) The examples the authors chose clearly demonstrate varying aspects and strengths of the method, and encompass different dimensionalities and stability properties (e.g., fixed points vs. limit cycles) across different application domains.\n\n\n\n\n## Weaknesses:\n### **Major:**\n1) The method feels neither fully theoretical nor purely data-driven (i.e., trajectories data without knowing the dynamics in the full phase space), but more conceptual. Since it sits somewhere in between, I think it’s important to highlight how it can actually help science in practice given data. You touch on this a bit in the discussion, but I’d appreciate one or two more specific sentences explaining how it could be applied, even in future work, using data or combined with data-driven approaches. For example, maybe first applying SINDy [1] to identify basis dynamic “atoms” to better understand the flow field, or using discrete data-driven models like SLDS [2,3] / dLDS [4,5] to identify transitions that then help infer the continuous flow field, which could then be used with your method. Notably, I do not suggest applying these approaches here, as that is beyond the scope of the paper, but rather recommend including a discussion on how your method could be combined with such data-driven techniques to characterize dynamics from data.\n\n2) I do not understand how the cost in Eq. (10) forces a sign change. It seems like it promotes non-zero values on both sides of the separatrix, but the eigenfunction could still have the same sign on both sides and simply pass through zero on the separatrix. Can you clarify this?\n\n3) I miss any comparisons to existing methods. While I understand that there are not many approaches that target idnetifying separatrices, some comparison to current tools (e.g. FTLE) is important.\n\n4) While the cost terms are explained, I feel the practical step-by-step process of the method could be clearer via an algorithm or pseudocode (e.g., using the algorithm package in LaTeX).\n\n\n### **Minor:**\n1) Line 95 — please refer to the specific appendix (seems like Appendix C?).\n2) Typos: 1) Line 124:  “We” should not be capitalized. 2) Line 168: extra “)”. 3) Line 197: “it’s” -> “its”.\n3) Throughout the paper, I felt the authors tried to connect their method to neuroscience, but I think the connection should either be more clearly emphasized or removed, as right now it feels somewhat in between and unclear (parsing-wise). I suggest either focusing the paper solely on dynamical systems understanding, perhaps mentioning neuroscience as just one example among others, or, alternatively, strengthening the neuroscience link. For example, the authors could emphasize more clearly why the RNN task is actually neuroscience-relevant, or discuss the neuro significance in greater detail,like how it can be helpful in the future for BCI interventions, or  for designing brain perturbations for motor control, etc.\n4) I would be more cautious when discussing perturbations. The method seems to assume strong conditions about the effect of the perturbation, for example, that it does not alter the state space or the separatrices themselves, which may not hold in all systems. It is important to clearly state such assumptions about the nature of the perturbations.\n5) In Eqs. (15–16), I do not think the Euclidean distance is necessarily the best measure for \"minimal perturbations.\" I think it would be helpful if you could provide more details about how you assume the perturbation occurs. For example, are they operators that simply move the system to a new state directly, or do they involve pushing the system against the flow field? If the former, why should the distance matter? If the latter, if the perturbation requires pushing the system to a new location in an energy-consuming way, it might be better to consider the flow field’s density and magnitude throughout the perturbation. Geometrical distance may not correlate with the actual effort needed to push against the flow. Providing a specific example of the type of perturbation you have in mind and explaining why would help. You could also mention that energy-consumption considerations might be added in future work and justify the current use of Euclidean distance for computational simplicity.\n6) I missed a colorbar for the background in all phase portraits. Even if the exact values don’t matter, having a colorbar in at least one of the figures (and indicating the direction—e.g., blue for 0, white for high values) would be helpful for interpretation.\n\n\n\n[1] Brunton, S. L., Proctor, J. L., & Kutz, J. N. (2016). Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15), 3932-3937.\n\n[2] Linderman, S. W., Miller, A. C., Adams, R. P., Blei, D. M., Paninski, L., & Johnson, M. J. (2016). Recurrent switching linear dynamical systems. arXiv preprint arXiv:1610.08466.\n\n[3] Linderman, S., Johnson, M., Miller, A., Adams, R., Blei, D., & Paninski, L. (2017, April). Bayesian learning and inference in recurrent switching linear dynamical systems. In Artificial intelligence and statistics (pp. 914-922). PMLR.\n\n[4] Mudrik, N., Chen, Y., Yezerets, E., Rozell, C. J., & Charles, A. S. (2024). Decomposed linear dynamical systems (dlds) for learning the latent components of neural dynamics. Journal of Machine Learning Research, 25(59), 1-44.\n\n[5] Chen, Y., Mudrik, N., Johnsen, K. A., Alagapan, S., Charles, A. S., & Rozell, C. (2024). Probabilistic decomposed linear dynamical systems for robust discovery of latent neural dynamics. Advances in Neural Information Processing Systems, 37, 104443-104470.", "was_revised": true, "final_justification": "I thank the authors for their response and believe that the clarifications, additional comparisons, and discussion points will strengthen the paper.\n\nFollowing point 1, one more suggestion for the discussion: I wonder if, when applying models like rSLDS, the separatrices might be captured by the switching regions between the LDSs, or if residual errors in model fitting could slightly distort their locations, making them harder to detect. It would be interesting to compare this (discussion-wise, as I understand more experiments are beyond the scope) to what more continuous models (e.g., those I mentioned in my review) might capture.\n\nI have raised my score to 5.", "rebuttal": "# Response to Reviewer dWGz\n\nWe thank the reviewer for the thoughtful and constructive feedback. We are especially grateful for the positive assessment of the paper’s clarity, examples, and overall presentation. Below, we respond to the specific concerns raised.\n\n---\n\n## Major Weaknesses\n\n**1. Application to real neural data / data-driven models**\n\nWe agree this is a critical aspect. Our method can be applied to any model that infers a flow field $f(x)$ from data. This includes GPSDE [1], rSLDS [2], gpSLDS [3].\n\nThese are all stochastic models, but one can apply our method considering a deterministic component or approximation of the inferred dynamics. We developed this idea in the case of rSLDS and GPSDE as they have been widely applied to neural data [9]. (See our response to reviewer FEbm). We tested our method on rSLDS trained on data from a synthetic bistable stochastic system (response to reviewer FEbm). In gpSLDS [3], uncertainty in the flow is also reported; this can be propagated to estimate the uncertainty in the separatrix.\n\nWe thank you for the suggestion regarding rSLDS and believe the result strengthens our work and would be happy to include it. We will expand our discussion to explicitly describe how our approach complements these data-driven models and can be applied post hoc to real data.\n\n**2. Clarification on Eq. (10): how does it promote sign changes?**\n\nThank you for catching this. The loss term in Eq. (10) not only promotes non-zero values but also encourages the eigenfunction $\\psi$ to have both positive and negative values across a region—thus promoting sign changes. This is achieved by minimizing the squared mean relative to the variance:\n\n$$\n\\frac{\\langle \\psi \\rangle^2}{\\langle \\psi^2 \\rangle} \\ll 1\n$$\n\nWhen the mean is close to zero but the variance is high, $\\psi$ necessarily has values of opposite sign. We will clarify this in the text and add annotations to Figure 3 to highlight this behavior.\n\n**3. Comparison to existing methods like Finite Time Lyapunov Exponent (FTLE)**\n\nThank you for bringing our attention to this work. FTLE methods characterise regions of sensitivity via divergence rates of trajectories. Ridges of the FTLE correspond to separatrices.\n\nWe draw some comparisons to our method.\n\nWhile our method yields a function which is zero along the entire separatrix, the FTLE takes non-zero positive value that may vary along the separatrix. This makes it more difficult to determine if a point is a separatrix or not.\nOur KEF-based method yields a *differentiable* function whose zero set explicitly defines a boundary; this enables efficient gradient-based search for minimal perturbations. Differentiating FTLE requires differentiating through the dynamical function which may not be possible (vanishing gradients) or may be computationally expensive.\nFTLE requires forward simulation for a time $T$, while our method avoids this step and may offer improved computational efficiency once the eigenfunction is trained.\n\nWe will include a brief comparison with FTLE in the related work section.\n\n**4. Include an explicit algorithm or pseudocode**\n\nWe agree this would improve clarity. We will include pseudocode using the `algorithm` package in LaTeX to summarize our method step-by-step in the revised version.\n\n---\n\n## Minor Issues\n\n**1. Line 95 — reference to specific appendix**\n\nWe will update this to explicitly refer to Appendix C.\n\n**2. Typos**\nWe thank the reviewer for pointing them out. We will make these corrections.\n\n**3. Neuroscience framing is unclear**\n\nWe thank the reviewer for this valuable advice. We can better highlight the relevance of our method to neuroscience. While our discussions address *data-trained* dynamical models of the brain, such models may also be *task-trained*. In the past, both types have generated insights about neural computations (see [4,5] for reviews). In particular the task-trained variants serve as hypothesis generation tools for how a certain computation could be carried out by a recurrent neural circuit [5]. Such models often posses dynamical objects like fixed-points, line-attractors, oscillations, and separatrices but are embedded in a high-dimensional space making it harder to identify. This motivates methods for finding/studying these objects [6]. The flip-flop is a prototypical task yielding multi-stable dynamics in RNNs, and is often used to demonstrate dynamical tools in neuroscience [7,8]. But separatrices are more general and appear in almost any scenario involving discrete choices, memories, threshold behavior [9].\n\nWe will add this to our introduction and discussion.\n\n\n**4. Assumptions about perturbations part 1**\nThank you, we did not clearly specify whether the input perturbations affect the recurrent dynamics, i.e., $f(x)$ versus $f(x,u)$ where u is an external input, e.g., a stimulus or neural perturbation.  In the paper we assumed that it does not $f(x)$. However, our formalism can be extended to $f(x,u)$ by defining $\\psi(x,u)$, in which case the Koopman PDE becomes\n$$\\nabla_x \\psi(x,u)\\cdot f(x,u)=\\psi(x,u).$$ We leave the development of this scenario to future work.\n\nWe will specify our assumption that the flow is independent of $u$ in our paper.\n\n\n**5. Assumptions about perturbations part 2**\n\nWe apologise, we did not clearly explain our perturbation model.\nIn the paper, we consider an idealized perturbation: a powerful Dirac delta pulse that instantaneously pushes the system from its current state to a new one, i.e., $x_\\text{base} + \\delta x$. Under this model, the optimal perturbation is simply the one with minimal norm, $|\\delta x|$, subject to $\\psi(x_\\text{base} + \\delta x) = 0$ — as shown in Eq. (15).\n\nA more general version of this problem would involve finding a time-varying input signal that drives the system across the separatrix while minimizing energy or power. This still requires identifying a final state $x_\\text{final}$ such that $\\psi(x_\\text{final}) = 0$, and minimizing the cost of reaching it (e.g., via a control integral along the trajectory from $x_\\text{init}$ to $x_\\text{final}$).\n\nOur current focus is on identifying the separatrix via $\\psi(x)$ — not on solving the control problem itself. However, the use of a smooth, differentiable eigenfunction to define the boundary is compatible with broader control formulations, and could serve as a basis for future work incorporating energy-aware perturbation metrics.\n\nWe will add these details to the main text.\n\n**6. Missing colorbars in phase portraits**\n\nWe thank the reviewer for this important point. We will add colorbars to the phase plots and clarify the meaning of background colors (e.g., blue = 0, white = high) to improve interpretability.\n\n---\n\n## Responses to Specific Questions\n\n**1. Clarify Eq. (10) and sign changes**\n\nSee response under Major Weakness #2.\n\n**2. Order of Eqs. (13) and (14)**\n\nYou are right — we should present Eq. (14) before Eq. (13), since we solve for individual KEFs first and then combine them. We will revise the order and clarify in the text. For example, in Figure 3, the two KEFs are shown separately in Figure 4.\n\n**3. Why use Euclidean distance for perturbations?**\n\nSee response under Minor Issue #5. In brief: we assume a direct displacement via a delta-pulse perturbation. In that case, Euclidean norm is appropriate. Future work could incorporate flow-aware or energy-based metrics.\n\n**4. Why use DNNs for KEFs instead of EM or variational methods?**\n\nWe use DNNs as a flexible function approximator to parameterize $\\psi$ and optimize the PDE residual. DNNs are expressive and differentiable, allowing us to not only solve the PDE but also optimize over $\\psi(x)$ post-training. While one could in principle use other parameterizations, we are not aware of EM formulations that apply here, since our setup does not involve a probabilistic latent variable model. If the reviewer had a specific alternative in mind, we would be happy to consider it.\n\n**5. Evaluation when no ground truth is available**\n\nIn systems without known separatrices, we can validate a candidate point $x$ on the predicted separatrix (i.e., $\\psi(x) = 0$) by evolving nearby initial conditions. If they flow to different attractors, this supports the claim that $x$ lies on a separatrix. This is related to FTLE-type validation.\n\nWe will include a discussion of this strategy in the final version.\n\n---\n\n## References\n\n[1] Duncker et al., *Learning interpretable continuous-time models of latent stochastic dynamical systems*, ICML 2019  \n[2] Linderman et al., *Bayesian learning and inference in recurrent switching linear dynamical systems*, AISTATS 2017  \n[3] Hu et al., *Modeling latent neural dynamics with Gaussian process switching linear dynamical systems*, NeurIPS 2024  \n[4] Vyas, Saurabh, et al. \"Computation through neural population dynamics.\" Annual review of neuroscience 43.1 (2020): 249-275.\n[5] Barak, Omri. \"Recurrent neural networks as versatile tools of neuroscience research.\" Current opinion in neurobiology 46 (2017): 1-6.\n[6] Sussillo, David, and Omri Barak. \"Opening the black box: low-dimensional dynamics in high-dimensional recurrent neural networks.\" _Neural computation_ 25.3 (2013): 626-649.\n[7] Versteeg, Christopher, et al. \"Computation-through-Dynamics Benchmark: Simulated datasets and quality metrics for dynamical models of neural activity.\" _bioRxiv_ (2025)\n[8] Ostrow, Mitchell, et al. \"Beyond geometry: Comparing the temporal structure of computation in neural circuits with dynamical similarity analysis.\" _Advances in Neural Information Processing Systems_ 36 (2023): 33824-33837.\n[9] Nair, Aditya, et al. \"An approximate line attractor in the hypothalamus encodes an aggressive state.\" _Cell_ 186.1 (2023): 178-193."}, {"review_id": "DaWqwR9wut", "reviewer": "Reviewer_MwkQ", "summary": "The authors propose a deep learning approach for identifying separatrices of dynamical flows. The method first trains a model to obtain Koopman eigenfunctions with real positive eigenvalues, and subsequently detects separatrices by leveraging their precisely vanishing characteristics at these boundaries. The effectiveness of the proposed approach is validated through multiple designed experiments.", "questions": "1. Has the author provided a comprehensive review of the relevant work on the discovery of separatrices? The article lacks comparative experiments with baselines, making it unclear whether the proposed method outperforms prior approaches.\n\n2. The training of Koopman eigenfunctions requires global data near separatrices. Can the authors illustrate the training data conditions? Is it possible to directly estimate the separatrices from these trajectory data from training set? What are the advantages of the proposed method over direct estimation?\n\n3. The paper lacks theoretical and experimental studies on the robustness of the proposed method. For example, how does it perform under complex scenarios such as noisy or partially observed data?", "limitations": "Yes.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:**\n\n1. The paper is logically structured and rigorously presented.\n\n2. Extensive experimental validation has been conducted across multiple simulation systems.\n\n3. The identification of separatrices holds significant application potential for the study of complex systems.\n\n**Weaknesses:**\n\n1. The paper lacks a detailed review of related work on separatrices discovery and comparative experiments with baseline methods.\n\n2. Training Koopman eigenfunctions appears to require global data near separatrices, yet it remains unclear whether such trajectory data can be directly utilized for separatrices discovery.\n\n3. The paper lacks both theoretical and experimental investigations into the robustness of the proposed method.", "was_revised": true, "final_justification": "Although the article has certain limitations in comparative validation with baseline methods, the overall work still holds significant theoretical value in the field of nonlinear dynamics. Therefore, I am pleased to raise my score to 4.", "rebuttal": "# Response to Reviewer MwkQ\n\nWe thank the reviewer for the helpful comments and constructive critique. Below, we address the main concerns raised regarding comparisons to prior work, the role of training data, and robustness.\n\n---\n\n**1. Review of related work and comparisons with baselines**\n\nWe agree that the paper would benefit from a broader discussion of existing methods for separatrix discovery. As also noted by Reviewer dWGz, one possible comparison is with approaches based on **Finite-Time Lyapunov Exponents (FTLE)**, which identify ridges of sensitivity in the flow. We now discuss FTLE and its differences from our method in our response to Reviewer dWGz. Briefly, FTLE does not identify zero-level sets, and its computation depends on simulation time $T$, while our approach yields a smooth, differentiable function whose zero-level set can be used directly to extract separatrices and optimize perturbations. We will revise the related work section to explicitly include this discussion.\n\n---\n\n**2. Can the separatrix be directly estimated from trajectory data?**\n\nThis is a very interesting question. Our current method assumes access to a continuous or inferred vector field $f(x)$ and learns Koopman eigenfunctions from it. We agree that directly estimating the separatrix from trajectory data — without modeling the flow — is an intriguing direction, though potentially difficult in practice. Separatrices are global structures and cannot generally be identified from local trajectory features alone.\n\nThus, we view the two-step process (first infer $f$, then infer the separatrix) as a natural and principled strategy. In neuroscience, for instance, models of dynamics can be learned from neural recordings or trained task-specific networks (such as in our flip-flop examples). Both modeling paradigms — data-driven and task-driven — have yielded important insights into neural computation [2,3]. As we explain in the response to Reviewer dWGz (\"**Neuroscience framing is unclear**\"), these insights are thanks to methods in this spirit, that aid finding dynamical objects in high-dimensions [4].\n\n---\n\n**3. Robustness under noise or partial observability**\n\nWe appreciate this suggestion. While our current study focuses on fully observed deterministic systems, we believe the method can be extended to stochastic and partially observed settings. For example, the **gpSLDS** model [1] infers latent dynamics from noisy observations and provides uncertainty estimates over the learned flow field. Our method can be applied to such learned $f(x)$ functions directly, and the reported uncertainty in $f(x)$ can propagate to estimate uncertainty in the inferred separatrix.\n\nWe mention this in our response to Reviewer dWGz and will add this future direction to the discussion section.\n\n---\n\n**References**\n\n[1] Hu, Amber, et al. *Modeling latent neural dynamics with Gaussian process switching linear dynamical systems*. NeurIPS 2024  \n[2] Vyas, Saurabh, et al. *Computation through neural population dynamics*. Annual Review of Neuroscience 43.1 (2020): 249–275  \n[3] Barak, Omri. *Recurrent neural networks as versatile tools of neuroscience research*. Current Opinion in Neurobiology 46 (2017): 1–6\n[4] Sussillo, David, and Omri Barak. “Opening the black box: low-dimensional dynamics in high-dimensional recurrent neural networks.” _Neural computation_ 25.3 (2013): 626-649."}]

Author rebuttal:

Meta-review: The authors of this work devised a method for tracing out separatrices in dynamical systems, the boundaries that separate different basins of attraction, based on Koopman operators with their eigenfunctions approximated by DNNs. They illustrate their approach on several low-dimensional examples, but also on two simulated high-dimensional systems, one from ecology and one an RNN trained on a one-bit flip-flop task.

This is an important topic, of considerable theoretical and practical interest across diverse fields dealing with complex, multi-stable dynamical systems, which hasn’t received that much attention so far. The authors are mindful of different failure modes (degeneracies) of their method, and provide sound solution strategies. Initially, reviews brought up a couple of issues regarding comparisons to existing methods, the broader applicability of the method, potential robustness issues, and potential scaling problems. In my mind these were mostly addressed by the authors. There are indeed not that many methods out there for detecting these manifolds, although there is the class of numerical continuation methods which the authors may have missed (and might want to include in their discussion). An application I would have found really interesting (also discussed by the referees) is on RNNs actually trained on real-world, e.g. neuroscience, data, not just simulated or artificial tasks. This could have made a much stronger point for the methods' virtues.

Overall, however, I find this a well-crafted and important contribution, with a solid technical part and several nice ideas along the way. I support acceptance, as did all three referees following the AC-referee discussion.
Decision: Accept (poster)