Characterizing control between interacting subsystems with deep Jacobian estimation
Authors: Adam Joseph Eisen
Mitchell Ostrow
Sarthak Chandra
Leo Kozachkov
Earl K Miller
Ila R Fiete
Venue: NeurIPS 2025 spotlight ()
Keywords: control theory
dynamical systems
neuroscience
deep learning
modular systems
Primary area: neuroscience_and_cognitive_science
TL;DR: A data-driven nonlinear control-theoretic framework to characterize subsystem interactions, leveraging a deep-learning method to learn dynamical system Jacobians.

Abstract: Biological function arises through the dynamical interactions of multiple subsystems, including those between brain areas, within gene regulatory networks, and more. A common approach to understanding these systems is to model the dynamics of each subsystem and characterize communication between them. An alternative approach is through the lens of control theory: how the subsystems control one another. This approach involves inferring the directionality, strength, and contextual modulation of control between subsystems. However, methods for understanding subsystem control are typically linear and cannot adequately describe the rich contextual effects enabled by nonlinear complex systems. To bridge this gap, we devise a data-driven nonlinear control-theoretic framework to characterize subsystem interactions via the Jacobian of the dynamics. We address the challenge of learning Jacobians from time-series data by proposing the JacobianODE, a deep learning method that leverages properties of the Jacobian to directly estimate it for arbitrary dynamical systems from data alone. We show that JacobianODE models outperform existing Jacobian estimation methods on challenging systems, including high-dimensional chaos. Applying our approach to a multi-area recurrent neural network (RNN) trained on a working memory selection task, we show that the “sensory” area gains greater control over the “cognitive” area over learning. Furthermore, we leverage the JacobianODE to directly control the trained RNN, enabling precise manipulation of its behavior. Our work lays the foundation for a theoretically grounded and data-driven understanding of interactions among biological subsystems.

Characterizing control between interacting subsystems with deep Jacobian estimation
Adam Eisen Brain and Cognitive Sciences MIT Cambridge, MA 02139
Mitchell Ostrow Brain and Cognitive Sciences MIT Cambridge, MA 02139 ostrowêmit. edu
jisena @mit edu
Sarthak Chandra"
Leo Kozachkoy
Brain and Cognitive Sciences MIT
IBM Thomas J. Watson Research Center
IBM Research
Cambridge, MA 02139 sarthak. chandra@icts. res. in
Yorktown Heights, NY 10598
1eokoz80brom. edu
Earl K. Miller Brain and Cognitive Sciences
lla R. Fiete Brain and Cognitive Sciences
MIT Cambridge, MA 02139 edu
MIT
Cambridge, MA 02139 fiete@mit.cm edu
Abstract
Biological function arises through the dynamical interactions of multiple subsys- tems, including those between brain areas, within gene regulatory networks, and more. A common approach understanding these systems is model the dy- namics of each subsystem and characterize communication between them. An alternative approach is through the lens of control theory: how the subsystems control one another. This approach involves inferring the directionality, strength, and contextual modulation of control between subsystems. However, methods for understanding subsystem control are typically linear and cannot adequately describe the rich contextual effects enabled by nonlinear complex systems. To bridge this gap, we devise data-driven nonlinear -hhtrri-t framework characterize subsystem interactions via the Jacobian of the dynamics. We ad- dress the challenge of learning Jacobians from time-series data by proposing the JaCobianODE, deep learning method that leverages properties of the Jacobian to directly estimate it for arbitrary dynamical systems from data alone. We show that JacobianODE models outperform existing Jacobian estimation methods on challenging systems, including nigh-dimensional chaos. Applying our approach to multi-area recurrent neural network (RNN) trained on a working memory selection task, show that the "sensory" area gains greater control over the "cog- nitive" area over learning. Furthermore, we leverage the JacobianODE to directly control the trained RNN, enabling precise manipulation of its behavior. Our work lays the foundation for a theoretically grounded and data-driven understanding interactions among biological subsystems.
"This author now at The International Centre Theoretical Sciences. This author is now at Brown University.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Introduction
Complex systems are ubiquitous in nature. These systems exhibit a wide range of behavior and function, in large part through the dynamic interaction of multiple component subsystems within them. One approach to understanding such complex systems is build detailed models of their underlying dynamics. An alternative and simpler yet powerful approach is offered by control theory, focusing instead on how subsystems influence and regulate one another, and how they can be controlled.
Control theory thus offers complementary approach to both understanding and manipulating biological systems. The theory describes how inputs must be coordinated with system dynamics to achieve desired behaviors, and can be applied across domains ranging from robotics to biology (Figure [A). The brain coordinates neural activity across multiple interconnected brain areas, dynamically modulating which regions receive information from which others depending on need and context
2 Interareal interactions play central roles in cognition and consciousness [377], in selective
attention 8 16 decision making [7 [8] working memory 19 20). feature binding (21 motor
control 23 27], and learning and memory (8 32].
common approach to characterizing interareal interactions is quantify communication between them, using methods such as reduced-rank regression to define "communication subspaces". These subspaces determine low-dimensional projections that maximally align high-dimensional states of the input area with high-dimensional states of the target area 35). However, effective control not only involves alignment of high-variance input states with high-variane target states, but also appropriate alignment of the inputs with the dynamics of the target area. Given connected subsystems and B, an identical signal from B will have dramatically different control effects on A, depending on whether the signal aligns with stable unstable directions of A's dynamics: projections onto more unstable eigenvectors can much more readily drive the system to novel states. (For more detail Appendix C.I
Accordingly, recent work in neuroscience has espoused control-theoretic perspectives on interareal interactions (Figure (B) [24 42 The dominant approach has involved linear control 43 43
However, the inherently nonlinear dynamies of the brain enable richer contextual control than possible to fully model with linear systems, necessitating nonlinear control approach.
A
B
contral theory
DOC
area A
area
Figure I: Schematic overview of control-theorett framework applied to neural interactions. (A) Control theory generalizes across diverse systems. (B) Illustration of interareal control, highlighting how neural activity in one area directly influences dynamics in another.
One approach to extend control-theoretic analyses to nonlinear systems is by linearizing the nonlinear dynamics through Taylor expansion, which involves the Jacobian. This converts the nonlinear system into a linear state- and time-dependent one 40 (4)5))) allowing for simple control. Jacobian linearization for control is straightforward with access analytical expressions for the non-linear system, but it becomes non-trivial in purely data-driven scenarios. Estimating the Jacobian involves conjunctively inferring both a function and its derivative, yet good function approximation need not yield good approximations of its derivatives (see Section 4 and Appendix
This paper introduces several key contributions:
Robust data-driven Jacobian estimation. We present JacobianODE, a deep learning- based method for estimating Jacobians from noisy trajectories in high- -igh-iimensoonal dynamical systems. We demonstrate the validity of this approach in several sample systems that include high-dimensional chaos.
A framework characterize control between interacting subsystems via data-driven Jacobian estimation. Harnessing Jacobian estimation method, we devise rigorous,

---

data-driven approach for nonlinear control-theorede analysis of how paired interacting systems, including brain areas, drive and regulate each other across different contexts. Data-driven inference of control dynamics in trained recurrent neural networks. We apply our data-driven framework to a recurrent neural network (RNN) trained on working memory task. We show that, purely from data, we can identify key interactions between the areas, and that these interactions crucially evolve over the course of learning to produce the desired behavior.
Demonstration of accurate control of rich interacting high-dimensional coupled dy- namical subsystems. We demonstrate high aceuracy in challenging nigh-dimensional data-driven nonlinear control task, enabling precise control of the behavior of the RNN.
Overall. our work lays the foundation for data-driven control therrettt analyses in complex high- dimensional nonlinear systems, including the brain.
Related work
Interareal communication A wide range of tools have been developed and harnessed to study
interareal communication in neural data 34 57]. This includes, but is not limited to, methods based
reduced rank regression [33 35 recurrent neural network models of neural dynamics (40 588 63),
Gaussian process factor analysis 64
65], canonical correlation analysis [66-68], convergent cross
mapping (19 switching dynamical systems [72 [73,, granger causality [74 75). dynamic causal
mapping 76.. point process models 777 and machine learning methods 78 79
Nonlinear controllability Classical results assess the controllability of nonlinear control systems via the Lie theory [80)8]] 84 Another approach to nonlinear network controllability is based on attractor strength [85] Liu et al. [54] and Parkes et al. (55) note that the large literature of linear controllability analyses could be extended locally to nonlinear systems with an appropriate lincarization method.
Neural network controllability Linear structural network controllability has been implicated across a wide range of contexts, tasks, and neuropsychiatrie conditions [36
42 44
50). This
approach has been extended to functional brain networks [31553] 53 Recent work has characterized
the subspaces of neural activity that are most feedback controllable as opposed to feedforward controllable using linear methods 43). Other approaches to neural control analyses consider the
identification of structural driver nodes 86 and input novelty (87]
Data-driven Jacobian estimation Besides approaches estimating Jacobians through weighted [88 (89). some methods have used direct parameterization via neural networks to learn Jacobians of general functions (90] 9I]. These approaches inform our method but do not explicitly address dynamical systems. Applying path-integral-bar
linear regressions
Jacobian estimation dynamical systems is challenging, as the target function (the system's 's time derivative) is typically unobserved. Siik-Mohammadi et al. [92] utilized this idea in dynamical systems to learn contracting latent dynamics from demonstrations.
3
JacobianODE: learning Jacobians from data
3.1 Jacobian linearization
We consider nonlinear dynamical systems in IR., defined by x(() f(x(t)). The Jacobian of the dynamics a matrix-valued function Jr IR" + Roxn (henceforth, J) given by
8f
0f
8/1
dag
0x.
0
Df2
0f
a 0x
(x(())
f(x(t))
Of .0m
af
afm
dag
Or...
At each time I, the Jacobian defines a linear subspace relating input and output changes, capturing how perturbations to the system will propagate. This recasts nonlinear dynamics as linear time-varying dynamics the tangent space locally along trajectories (formally, (x = Jr(x(1))xx(I), see Figure 2 left, also see Lohmiller and Slotine (93] for discussion in the context of virtual displacements).

---

(xx(t)
I(x(t))
6x(1')
area A
R
X
6x
game
g
8x
x
6x
area B
Jacobian linearization
controi systems
+
5*
J
6x
xx
Time in trajectory
with within-area dynamiics
control
(x(() = ((x())xxx()
Figure 2: Analytical framework pairwise interacting subsystem control. Trajectory dynamics (left, top) locally linearized via Jacobians (left, bottom), explicitly separating within-area (diagonal blocks) and interareal (off- ooff-diagonal blocks) dynamics (right). These separated dynamics can be used to construct interaction-speeifii control systems.
path integrate around loop
B
C
At
t
path intogrtte
j
CDE estimate
gencrated predictions toacher forood prodictions
path integrate
(x((t))
K(ty
Jds
to estimate
Gwop
loop closure loss d'
D
trajectory data
neural network Jacobian estimate
path integration Rictt)) f(x66))
integration Ar)
training
:r(t)
trajectory loss 4xc++ an) Ar))
teacher forcing R((+) Ar)
aj8i A0) oox(
Figure 3: Jacobian estimation with JacobianODE models. (A) Path integration of the Jacobian predicts future states. (B) Generalized teacher forcing stabilizes trajectory predictions during training. (C) Loop-closure constraints enforce consistency of Jacobian estimates. (D) Training pipeline, combining neural Jacobian estimation, path integration, teacher forcing, and self-supervised loop- closure loss.
3.2 Parameterizing differential equations via the Jacobian
We now turn to the problem of how to estimate the Jacobian from data. We assume that we have access only to observed trajectories of the system, of the form (x) + k&t),, == 1,2,... j 1, 2, ... where At is a fixed sampling interval, j indexes the trajectory, and ((5) is a trajectory-speciff start time. Crucially, directly via neural network. To do this, we parameterize neural network function J" with learnable parameters 9 that then trained to approximate J.
do not assume access to the function f. Our method estimates the Jacobian
Path integration Following Lorraine and Hossain 91 and Beik-Mohammadi et al. 92 we exploit the fact that the path integral of the Jacobian path independent in the construction of the loss function. This is because the rows of the Jacobian are conservative vector fields (i.c., they the

---

gradients of scalar functions). For an intuitive picture, consider the work done by the force of gravity as you climb mountain. Regardless of the path you take to elimb, the resulting work is dependent only on the start and end points of the path. Formally, for the time derivative function f we have that
[c Jds :-==/ S.a J(c(r)) dr,
f(x(t,)) - f(x(t)))
(2)
where C is a piecewise smooth curve in IR"" and c [4i,ty] + C is a parameterization of C' with c(4)) = x(4+) and c(ly) x(ly) (Figure 3A). Given estimates of f(x((:)) and the Jacobian J, any time f For these integrals, we use line between the endpoints as a simple choice of path (Figure [3A]. Then, to generate an estimate x(ty + At) of the next step, we can use a standard ordinary differential equation (ODE) integrator
we can then use Equation
2 to approximate (xx())))
(e.g., Euler, fourth-order Runge- Kutta, etc.) to integrate the estimated time derivative f(x(t)) (see Appendix for more detail). To avoid the need to represent f directly (thereby enabling all gradients to backpropagate through the Jacobian network) we note that we parameterize an estimate of (x(()) in terms of the Jacobian (see Appendix A.1
in Equation 2
3.3 Loss functions
Trajectory reconstruction loss Given an observed trajectory x(to + length T. we can compute the trajectory reconstruction loss, Lay (0; between the true trajectory and the estimated trajectory using an appropriate distance measure d (e.g.. mean squared error, see Appendix A.2 for more detail on generating predictions and the trajectory prediction loss).
of
Generalized Teacher Forcing To avoid trajectory divergence in chaotic noisy systems, we employ Generalized Teacher Forcing, generating recursive predictions partially guided by true states (Figure BB and Appendix D
Loop closure loss The Jacobian captures how perturbations to the system will propagate along any direction in state space. Estimating it purely from dynamics constrains only the direction of the flow, leaving the full solution underdetermined. To address this, we again exploit the fact that each row of the Jacobian is conservative vector field. Specifically, we note that for any piecewise
smooth loop Choup we have | Jds 2 = 0 (see Figure 3 C). Thus, by integrating along loops that
contain directions orthogonal to the system's 's dynamics (and penalizing the deviation from zero), encourage the estimated J Jacobians capture information about other directions in state space for full technical details). To ensure broad coverage of tangent space directions, form loops from concatenations of line integrals between randomly selected data points. This strategy samples diverse directions from the tangent space while remaining easy to compute. The resulting self-supervised loss term, Liep x). builds on the loss introduced by lyer et al. [95]. It constrains j to satisfy both the dynamics and conservativity. This improves Jacobian estimation accuracy significantly (see Appendix C.4 for ablation studies).
(see Dp.endices/AAJJandDDD
Training loss We therefore minimize the following loss function with respect to the parameters A:
[(0; Lung (0; x) + Atoop x0
(3)
where Atoop controls the relative weighting of the loop closure loss Cloop (0;x) compared to the trajectory prediction, and is hyperparameter of the learning procedure (Figure 3D).
4 Jacobian estimation in dynamical systems
Data To evaluate the quality of the Jacobian estimation procedure, we apply our approach to several example systems for which the dynamics are known. For this analysis, we used the Van der Pol oscillator (96), Lorenz system [97], and the Lorenz 96 system across three different system sizes (12,
32, and 64 dimensional) 198] All systems were simulated using the dysts package, which samples dynamical systems with respect to the characteristic timescale T of their Fourier spectrum [99] [00]. All training data consisted of 26 trajectories of 12 periods, sampled 100 time steps per All models were trained on a 10 time-step prediction task with teacher forcing.
To evaluate the performance of the methods in the presence of noise, we trained the models on data with 1 %,, 5%, and 10% Gaussian observation noise added i.i.d over time, where the percentage is defined via the ratio of the cuclidean norm of the noise to the mean euclidean norm of the data.
5

---

JacobianODE model For the JacobianODE framework, loop closure loss weights were chosen via line search (Appendix D.8.8). where the neural network J was taken to be four-layer multilayer perceptron (MIP) with hidden layer sizes 256, 1024, 2048 and 2048. Path integration was performed using the trapezoid method from the torchquad package, with each integral discretized into 20 steps [101]. ODE integration was performed using the fourth-order Sunge-Kutte (RK4) method from the torchdiffed package [102]. The JacobianODE models used 15 observed points to generate the initial estimate of f (see Section 3.2 and Appendices A.1 and A.2). All models were built in Py PyTorch
103.. Full implementation details are provided appendices A and
Baselines We chose two different Jacobian estimation procedures for comparison. The first was neural ordinary differential equation (NeuralODE) model trained to reproduce the dynamics [102]. The NeuralODE was implemented as four-layer MLP with hidden layers of the same size as the one used for the JacobianODE model. Jacobians were computed via automatic differentiation. NeuralODEs were regularized via a penalty on the Frobenius norm of the estimated Jacobians to prevent the model from learning unnecessarily large negative eigenvalues (see Appendix (D3) [104 [06]. We also employed baseline that estimates Jacobian via weighted linear regression, which computes locally linear models at each point in the space 88 89 (see Appendix D.4).
Table 1:
Mean Frobenius norm error on Jacobian estimation, [IJJ - )) for each system and noise level. Errors are reported as mean standard deviation, with mean and standard deviation computed over five random initializations of the model architectures.
Project
Training noise
JacobianODE
NeuralODE
Weighted Linear
VanDerPol (2 dim)
1%
0.7 = 0.1 0.72 # 0.05
1.0 = 0.3 0.72 # 0.08
6.11
5%
6.09
10%
1,35 + 0.06
2,2 # 0.2
6.08
Lorenz (3 dim)
1 %
3.3 + 0.2
8.7 0.3
21.94
5%
5.1 # 0.9
26.0 1.5 26.7 # 0.9
21.90
10%
6.4 + 0.1
21.84
Lorenz 96 (12 dim)
%
1.2 + 0.2
4.8 0.2
28.67
5%
2.7 = 0.2
5.9 = 0.2
28.64
10%
4.6 + 0.1
6.1 = 0.1
28.56
Lorenz 96 (32 dim)
1 %
8.7 + 0.2
16.8 + 0.4
47.13
5%
7.8 # 0.2
17.7 = 0.6
46,96
10%
13.45 & 0.09
19.5 0.4
47.03
Lorenz 96 (64 dim)
1
30.9 + 0.5
45.5 + 0.1
66.39
5%
34.0 0.2
45.7 # 0.1
66.26
10%
36.0 # 0.2
46.0 + 0.2
66.29
task-trained RNN
1%
188.5 7.1
294.02 + 0.03
301.46
5%
166.8 3.6
294.357 + 0.003
297.63
10%
180.4 0.5
294.328 # 0.004
296.63
Performance We tested the approaches on Jacobian estimation on held-out trajectories without noise. The JacobianODE method outperforms the baseline methods in terms of mean Frobenius norm error for virtually all systems (Table T). This was also true when considering the spectral matrix 2-norm (see Table S3 in Appendix
C.3
We plot performance in Figure 4 While both the JacobianODEs and the NeuralODEs reproduce the observed dynamics (Figure 4A-D), the JacobianODE learns more accurate estimate of the Jacobian (Figure 4 E,F). In particular, looking at Lyapunov spectra learned by the models, we note that the JacobianODE exceeds the other methods in estimating the full Lyapunov exponent spectrum, indicating a better overall representation of how perturbations along different directions will evolve (Figure 4G.H).

---

state space
10 time-step trajectory prediction
Jacobian estimation
Lyapunov spectrum
E
MSE
MSE
0.25
0.00
Expenent
-20
R²
-00
-30
10%
1%
5%
10%
Observation noise
Observation noise
Lyapunov exponent A
B
D
H
MSE
.**
MSE
Expoment
+**
0.5
0.5
4
0.0
1%
5%
10%
0.0
1%
5%
10%
50
[yapunow exponent
Observation noise
Observation noise
JacobianODE
NeuralODE
Weighted linear Jacoblans
Ground truth system
Figure 4: JacobianODE surpasses benchmark models on chaotic dynamical systems. Error bars indicate standard deviation, with statistics computed over five different model initializations. (A.B) State-space trajectories for (A) Lorenz and (B) 64-dimensional Lorenz 96 systems, with time-step predictions. Spectral embedding was used to reduce the Lorenz 96 data to three dimensions. (C,D) Accuracy 10 time-step trajectory predictions at varying noise levels. (E,F) Comparison of Jacobian estimation accuracy, quantified by mean squared error (MSE) and R (G,H) Estimated Lyapunov spectra averaged over initializations.
5
Control-theoretic analyses in a task-trained RNN
5.1 Characterizing control between subsystems with Jacobian lincarization
Consider neural data recorded from two areas, A and B. Concatenating their data into a state vector RR composed of x E TR and xB E RR and assuming dynamics governed by f, linearize around a reference trajectory (óx(() = (x(()) 5x(()). Splitting the Jacobian into block matrices yields:
(x)(1) = (x1+A(((1)) 8x +J B- (x(1))xx
5x 5x [A-B((((1) 8x + J (x(1))
where diagonal blocks J A RRAAXA and J B- R represent within-ares dynamics, and off-diagonal blocks B R"axRE, JA R represent interareal interactions (Figure right). Explicit separation of each area's control dynamics quantifies the direct influence each exerts on the other, and readily generalizes beyond two areas (see Appendix B.3
Since Jacobians are time-dependent, we obtain linear time-varying representation of control dynamics along the trajectory. This enables computation of time-varying reachability ease, capturing how readily each area drives the other toward novel states [09]. Reachability is quantified via the reachability Gramian, matrix defining a local metric in tangent space. For the above control system capturing the influence of area on area A, the time-varying reachability Gramian on the interval [to, t1] is defined as
Wr I" (t... B() (T) T),
(5)
where & (computed from A- denotes the state-transition matrix of the intrinsic dynamics of subsystem A without any input (i.e., 5x (t) = (t,, (t)oxx(tt))) and B(r) (x(+))((T)) (see Appendix B.1 [07]. The Gramian W, W,(to, (1) is symmetric and positive semidefinite for every > to [107]. Each eigenvalue of the reachability Gramian quantifies how easily the target system can be driven along its corresponding eigenvector. Thus, the trace of the reachability Gramian reflects average ease of reaching new states, and its minimum eigenvalue reflects the ease along the most challenging direction of control.
7

---

B
pre-stimulus
stimuli
delay 750 ms
cue
response
cognitive area
visual area
100 ms
750 ms
250 ms
Input Matrix W
Output Matrix W
Weight Matrix Wnn
D
ugper sttr,
Outgring Neurnn
0.2
Weight Stre-Streggg
Ourgoing Neuron
0.2
0.4
Weight Strongth
0.2
Weiggtt Stree gtth
64
0.0
Incominn
64
0.0
-0.2
-0.4
-0.4
128
128
00C088CAY Incoming Input
128
Incoming Neuron
Outgoing Outppt
G
1% Noise
10%
MSE
MSE
0.005
0.000
Expoment
-12
-10
-2
-20
R2
0.0
-30
1%
5%
10%
1%
5%
10%
Observation noise
Observation noise
Lyapunow exponent
JacobianODE
NeuralODE
Weighted Inear Jacotians
Ground truth system
Figure 5: Jacobian accurately infers trained RNN Jacobians. Error bars indicate standard deviation, with statistics computed over five different model initializations. (A) Task schematic, involving stimuli presentation, delays, cueing, and response. (B) RNN architecture with distinet visual and cognitive areas interacting. (C-E) Structure trained weight matrices: (C) recurrent weights within RNN, (D) input connectivity pattern, and (E) output connectivity. (F) State-space trajectories from the trained RNN, with 10 time-step predictions. Spectral embedding was used reduce the data to three dimensions. (G,H) JacobianODE performance evaluated against other models for trajectory prediction (G) and Jacobian estimation accuracy (H) at multiple noise levels. (I) Estimated Lyapunov spectra averaged over initializations.
5.2 Estimating the Jacobian of a task-trained RNN
Task To demonstrate how JacobianODE models could be used in neuroscience, we performed control-thcoretic analysis of a task-trained RNN. We used a working memory selection task from Panichello and Buschman 110 (Figure 5 On each trial, the network presented with two of four possible "colors", denoted by one-hot vectors. After a delay, the network is presented with cue indicating which of the colors select. The network is then asked to reach state of sustained activation that corresponds to the selected color.
RNN model To perform this task, trained I 128-dimensional continuous-time RNN. The RNN had hidden dynamics defined by
rh(t) = h (mo(h(()) + Wnau(t) + o(t) W hoth(t)
(6)
where War defines the internal dynamics, Wh: maps the input into the hidden state, Woh. maps the hidden state to a niu--dimnnsioal output (t),, b is a bias term, and a is the exponential linear unit activation. The RNN had two 64-neuron areas: "visual" area (which received sensory input) and "cognitive" (which output the RNN's ' color choice) (Figure BB). To encourage multi-area structure, we initialized the within-area weights with greater connectivity strength than the across-area weights (Figure 5C). Since input comes only to the visual area and output only from the cognitive area, the two areas are forced to interact to solve the task Piggure5p.F. The RNN solves this task with 100% accuracy.
8

---

Jacobian reconstruction quality We trained JacobianODEs on RNN trajectories from the post-cue delay and response portion of the trials. We used the same baselines as in Section JacobianODEs are not only robust to noise, but can also benefit from it (for multiple systems in Table I 5% training noise improves estimation). Noise encourages the model to explore how perturbations around the observed trajectory evolve, which is crucial for learning accurate Jacobians in nigh-dimensional systems. Thus (for both the JacobianODE and NeuralODE) we add a small amount of additional noise to the data during learning. Although the models perform similarly on trajectory reconstruction (Figure 5FG), on Jacobian estimation, JacobianODEs drastically outperform both baseline models
(Figure H.I, and Table I
A
B
Reecchatilty trape)
Ground truth
JaCebienODE
Ground fnuth reahasiit
-225
Ground truth
JaD
Gramion
exjusio
+0
100
Reachabiity
33
2
Cramian
0.01
vioual
Time from della start (4)
oagriive
0.00
cogritio
Ground
Ground tuth
Visual
cognitive
oagriiv
vocol " oogretive
JacobianDDE
cogn Live visual
JoobbinC0E
Early
Early
Wassl oogritive
cogn lvw veuus!
D
F
cognitive area
visual area
Resporse Accuracy
*** 0.92
Moan Sqaadd
M 0
1.77
0,20
0.11
Graand JusoteerODE NeuralOCE
Liree
Grann
Lineer
Figure 6: JacobianODE reveals differential interareal reachability. Error bars indicate standard deviation for (A) and (C), and standard eтror for (B), (E), and (F), with statistics computed over trajectories. (A) Comparison of interareal reachability between ground truth and JacobianODE estimates. (B) Temporal evolution of reachability (Gramian trace) throughout the delay period. (C) Comparison of reachability in early and late learning.
(D) Schematic illustrating targeted control of cognitive-area via visual-area stimulation. (E) Mean accuracy of targeted responses using JacobianODE versus alternative methods. (F) MSE for multiple ILQR- -asee methods.
5.3 Reachability in the task-trained RNN across contexts
Next, we used the Jacobians learned by the JacobianODE (trained on 5% noise) to evaluate reachability control in the RNN, and compared it to evaluations using ground truth Jacobians. All analyses were performed using the 10 time-step reachability Gramian. We first found that was on average easier for the visual area drive the cognitive area, both when considering overall ease (Gramian trace) and worst-case ease (Gramian minimum eigenvalue) (Figure 6/A, larger values indicate greater reachability). We next considered how reachability ease varied throughout the delay period. We found that the visual area could drive the cognitive area more easily at the beginning of the delay, with ease decreasing into the middle of the delay period (Figure bottom). The cognitive area was able to drive the visual area more easily slightly later in the delay period (Figure top). Finally, we considered whether reachability changes over the course of learning. We found that both directions of reachability increased after learning (Figure 6 (C). The JacobianODE reproduced all results accurately (Figure A-CC, comparison with ground truth). Our analysis reveals that reachability is crucial in the RNN's s ability to perform the working memory task, with the visual area's re's ability to drive the cognitive area shortly after the cue being especially important.
5.4 Controlling the task-trained RNN
To further validate our approach, we used JacobianODE Jacobians to control the RNN (Figure()) 6)). Given a trial in which the network was cued to respond with a particular color, we tested if we could induce a specific incorrect response by input to the visual area alone. To do so, we implemented Iterative Linear Quadratic Regulator (ILQR) control, which relies on knowledge of the Jacobian |ОII 222. The controller guided the network towards the mean hidden state of

---

training trials corresponding to the desired incorrect color. We defined accuracy as the percentage of time points during which the RNN output the desired color. We found that the JacobianODE widely outperformed the baseline models on this task, achieving accuracy nearing that of the ground truth system (Figure ( The JacobianODE was furthermore able to achieve the lowest mean squared ervor on the desired control trajectory (Figure This illustrates that while both the JacobianODF and the NeuralODE can learn the dynamics, only the JacobianODE learns a representation of the Jacobian that is sufficient for control.
6 Discussion
Extended Jacobian estimation A natural extension of JacobianODE models would add inductive biases for particular classes of dynamics. For example, one could parameterize a negative definite matrix and thus ensure contracting dynamics [92]. Other extensions could include an L1 penalty to encourage sparsity, as well as the inclusion of known dynamic structure (c.g., hierarchical structure, low-rank structure, etc.). In neuroscience, connectomic constraints could be incorporated to capture interactions within a neural circuit.
Limitations A future challenge for JacobianODE models is partially observed dynamics. Recent work has identified that is possible to learn latent embeddings that approximately recover the true state from partial observation [ Jacobian-basee babbed dynamics learning has been performed in latent space [92] however is unclear whether this translates to accurate Jacobian estimation, given the challenges related to automatic differentiation and Jacobian estimation presented here. We also note that reachability estimates depend sensitively on several factors: the alignment of cross-subsystem interactions and within-subsystem dynamics, the eigenvectors and eigenvalues of within-subsystem Jacobians, and the way activity propagates within each subsystem (see Appendix C.I.. While our method reliably captures broad trends in reachability over time, fine ane-grained, time- point-specific comparisons should be interpreted with caution. Finally, although JacobianODE models scale well to moderately
dinensinnnii systems, their performance in systems that are orders of magnitude larger than those considered here (e.g., recordings of thousands of neurons via calcium imaging) remains be tested. In many practical settings, this may not be a limitation: neural representations during tasks often exhibit intrinsic dimensionalities comparable to those studied here, enabling JaC models to operate in reduced dimensionality (see Appendix C.7). Notably, most models converged in under 45 minutes on single GPU, suggesting favorable scaling (see Table S7) Future work should explore the viability of the method in higher dimensions, and assess whether dimensionality reduction strategies (such as latent state models or low-rank Jacobian approximations) can further improve scalability and accuracy.
Broader impact on dynamical systems We demonstrated in this work that it is possible to learn high accuracy Jacobians of arbitrary dynamical systems from data. While in this work we focus on control- theoretic analyses, we emphasize that JacobianODE models can be trained in any dynamical setting. Jacobians are widely applicable, and are essential in stability analysis 80 119, contraction
theory [40] 80 93 126) and nonlinear control [56, 112 127 129 among many other settings.
Broader impact on control-theored analyses The control-theoretic analysis presented here is performed with respect to a particular reference trajectory. It therefore describes the local con- trollability properties along the reference trajectory, rather than global properties for the overall system [56]. This locality is desirable in biological settings due the context-sensitive nature of the relevant processes. Applications of this approach could involve a comparison of the functioning of the autonomic nervous system, gene regulatory networks, between control and disease states, given the association of these processes with disease [30 Future work should investigate whether a more global Jacobian-based control analysis (e.g., contraction theory) is achievable in data-driven manner via JacobianODE models. This would be useful especially in engineering settings where global controllability desirable, such robotics and prosthetic limbs.
Acknowledgments and Disclosure of Funding
The authors would like to thank Andrew Kirjner, Laureline Logiaco, Federico Claudi, Lakshmi Narasimhan Govindarajan, Jaedong Hwang, and Akhilan Boopathy for providing stimulating discus- sion about this work. The authors also acknowledge the MIT Office of Research Computing and Data
10

---

for providing high-performance computing resources that have contributed to the research results reported within this paper. I.R.F. acknowledges funding from The National Science Foundation Computer and Information Science and Engineering Directorate, The Simons Collaboration on the Global Brain, and The McGovern Institute MIT. E.K.M. acknowledges funding from ONR MURI 100014-23-1-27666 NIMH 1R13177510000 NIMH RO1 MH1 1559, The Simons Center for the Social Brain, The Freedom Together Foundation, and The Picower Institute for Learning and Memory, L.K. acknowledges funding from the Goldstine Postdoctoral Fellowship at IBM Research. M.O. acknowledges funding from NSF GREP 2141064.
References
[1] Laura Lee Colgin, Tobias Denninger, Marianne Fyhn, Torkel Hafting, Tora Bonnevie, Ole Jensen, May-Britt Moser, and Edvard I. Moser. Frequency of gamma oscillations routes flow of information in the hippocampus. Nature, 462(7271):353-337, November 2009. ISSN 1476-4687, doi: 10. 10.108//ature00553. URL /ttts:s::///doi /ttp::///oooooo///////o wttp://..///...///tiuiiiii. 1008/nature085//
[2] Shencong Ni, Brendan Harris, and Pulin Gong. Distributed and dynamical communication: a mechanism for flexible cortico-cortical interactions and its functional roles in visual attention. Commun. Biol., 7(1):550, 8 May 2024.
[3] F Crick and C Koch. Towards a neurobiological theory of consciousness. Seminars in the Neurosciences, 2:263-275, 1990.
[4] Andreas K Engel and Pascal Fries. Chapter 3 neuronal oscillations, coherence, and con- sciousness. Steven Laureys, Olivia Gosseries, and Giulio Tononi, editors, The Neurology of Conciousness (Second Edition), pages 49-60. Academic Press, San Diego January 2016.
[5] Evan Thompson and Francisco J Varela. Radical embodiment: neural dynamics and conscious- ness. Trends Cogn. Sci., 5(10):118 425, October 2001.
[6] Lawrence M Ward. Synchronous neural oscillations and cognitive processes. Trends Cogn. Sci., 7(12):553-559 December 2003.
[7] Markus Siegel, Tobias H Donner, and Andreas K Engel. Spectral fingerprints of large-scale neuronal interactions. Nat. Rev. Neurosci, 13(2):121-134, 11 January 2012.
[8] Satu Palva and Matias Palva. New vistas for Sipha-frequency band oscillations. Trends Neurosci., 30(4):150-155 April 2007.
[9] Timothy Buschman and Earl Miller. Top-down versus bottom-up control attention in the prefrontal and posterior parietal cortices. Science, 3 1866580))1880-18666 30 March 2007.
[10] Georgia G Gregoriou, Stephen J Gotts, Huihui Zhou, and Robert Desimone. High-frequency. long-range coupling between prefrontal and visual cortex during attention. Science, 324(5931): 1207-1210, 29 May 2009.
[11] Georgia G Gregoriou, Stephen Gotts, and Robert Desimone. Cell-tyee-specific synchroniza- tion of neural activity in FEF with V4 during attention. Neuron, 73(9):511554, 9 February 2012.
[12] E Salinas and TJ Sejnowski. Correlated neuronal activity and the flow of neural information. Nat. Rev. Neurosci, 2(8):539-550, August 2001.
[13] P Fries, J H Reynolds, A E Rorie, and R Desimone. Modulation of oscillatory neuronal synchronization by selective visual attention. Science, 29165508):1500-1553 23 February 2001.
[14] Markus Siegel, Tobias H Donner, Robert Oostenveld, Pascal Fries, and Andreas K Engel. Neuronal synchronization along the dorsal visual pathway reflects the focus of spatial attention. Neuron, 60(4):709-719, 26 November 2008.
[15] Jochem van Kempen, Marc A Gieselmann, Michael Boyd, Nicholas A Steinmetz, Tirin Moore, Tatiana Engel, and Alexander Thiele. Top-down coordination of local cortical state during selective attention. Newron, (09(5):894-9444 -98..90...8,.. 3 March 1021.
11

---

[16] C M Gray, P König, A K Engel, and W Singer. Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties. Nature, 338 (6213):334-33) 337, 23 March 1989.
[17] Markus Siegel, Andreas K Engel, and Tobias H Donner. Cortical network dynamics of perceptual decision-making in the human brain. Front. Hum. Neurosci., 5:21, 28 February 2011.
[18] Markus Siegel, Timothy Buschman, and Earl K Miller. Cortical information flow during flexible sensorimotor decisions. Science, 348(6241):1352-1355 19 June 2015.
[19] Scott L Brincat, Jacob A Donoghue, Meredith K Mahnke, Simon Kornblith, Mikael Lundqvist, and Earl K Miller. Interhemispheric transfer of working memories. Neuron, 109(6); 1055- 1066.e4, March 2021.
[20] R F Salazar, N M Dotson, S L Bressler, and C M Gray. Content-specific fronto- canto-pariital synchronization during visual working memory. Science, 33861110) 23 November 2012.
[21] A K Engel and W Singer. Temporal binding and the neural correlates of sensory awareness. Trends Cogn. Sci., 161): January 2001.
[22] W Singer and C Gray, Visual feature integration and the temporal correlation hypothesis. Annu. Rev. Neurosci., 18:555-586 5866,558,111 1995.
[23] V N Murthy and E E Fetz. Oscillatory activity in sensorimotor cortex of awake monkeys: synchronization of local field potentials and relation to behavior. J. Neurophysiol, 76(6): 3949-3967, December 1996.
[24] Laureline Logiaco, F Abbott, and Sean Escola. Thalamic control of cortical dynamics in model of flexible motor sequencing. Cell Rep., 35(9): June 2021.
[25] Fritzie I Arce-McShane, Callum F Ross, Kazutaka Takahashi, Barry Sessle, and Nicholas G Hatsopoulos. Primary motor and sensory cortical areas communicate via spatiotemporally coordinated networks at multiple frequencies. Proc. Natl. Acad. Sei. U. S. A., 113(18):5083- 5088, 3 May 2016.
[26] Matthew G Perich, Juan A Gallego, and Lee E Miller. A neural population mechanism for rapid learning. Neuron, 100(4)964-997..7, 21 November 2018.
[27] Matthew T Kaufman, Mark M Churchland, Stephen I Ryu, and Krishna V Shenoy. Cortical activity in the null space: permitting preparation without movement. Nat. Neurosci., 17(3): 440 448, March 2014.
[28] Scott L Brincat and Earl K Miller. Frequency-specific caappocampal-prffroont interactions during associative learning. Nat. Neurosci., 18(4):57-5881 April 2015.
[29] Matthew W Jones and Matthew A Wilson. Theta rhythms coordinate hippocampal-preet preffronal interactions in a spatial memory task. PLoS Biol., 3(12)xe402, December 2005.
[30] A G Siapas and M A Wilson. Coordinated interactions between hippocampal ripples and cortical spindles during slow-wave sleep. Neuron, 21(5):1123-112 November 1998.
[31] Antonio Fernández-Ruiz, Azahara Oliva, Marisol Soula, Florbela Rocha-Almein Gergo A Nagy, Gonzalo Martin-Vazquez, and György Buzsáki. Gamma rhythm communication between entorhinal cortex and dentate gyrus neuronal assemblies. Science, 372(6537), 2 April 2021.
[32] Evan G Antzoulatos and Earl K Miller. Increases in functional connectivity between prefrontal cortex and striatum during category learning. Newron, 83(1):216-2255, 2 July 2014.
[33] João D Semedo, Amin Zandvakili, Christian K Machens, Byron M Yu, and Adam Kohn. Cortical areas interact through a communication subspace. Neuron, 102(1):249-253)) e4,
3 April 2019.
12

---

[34] João D Semedo, Evren Gokcen, Christian K Machens, Adam Kohn, and Byron M Yu. Statisti- cal methods for dissecting interactions between brain areas. Curr. Opin. Neurobiol., 65:59-69, December 2020.
[35] Camden MacDowell, Alexandra Libby, Caroline I Jahn, Sina Tafazoli, and Timothy Buschman. Multiplexed subspaces route neural activity across brain-wied wide networks. bioRxiv, 12 February 2023.
[36] Shi Gu, Fabio Pasqualetti, Matthew Cieslak, Qawi K Telesford, Alfred B Yu, Ari E Kahn, John D Medaglia, Jean M Vettel, Michael B Miller, Scott T Grafton, and Danielle S Bassett. Controllability of structural brain networks. Nat. Commun., 6:8414, 1 October 2015.
[37] Ta-Chu Kao and Guillaume Hennequin. Neuroscience out of control: -heoreic perspectives on neural circuit dynamics. Curr. Opin. Neurobiol., 12-112,, October 2019,
[38]
Danielle S Bassett and Olaf Sporns. Network neuroscience. Nat. Neurosci, 20(3):353-3644 23 February 2017.
[39] Ta-Chu Kao, Mahdieh S Sadabadi, and Guillaume Hennequin. Optimal anticipatory control as a theory of motor preparation: A thalamo-cortical circuit model. Neuron, 109(9): 1581.e12, 5 May 2021.
[40] L Kozachkov, M Ennis, and J Slotine. RNNs of RNNs: Recursive construction of stable assemblies of recurrent neural networks. Adv. Neurul Inf. Prucess. Syst., 2022.
[41] Marine Schimel, Ta-Chu Kao, Kristopher T Jensen, and Guillaume Hennequin. iLQR VAE control-based learning of input-driven dynamics with applications to neural data. In Interna- tional Conference Learning Representations, October 2021.
[42] Sarah Feldt Muldoon. Fabio Pasqualetti, Shi Gu, Matthew Cieslak, Scott T. Grafton, Jean M. Vettel, and Danielle Bassett. Stimulation- control of dynamic brain networks, PLOS Computational Biology, 12(9):1-23, 09 2016. doi: 10. 10.137// httpppcci.....1.. 1005076. URL
https://aoi.org/ 10. 1 journal. pcbi 1005076
[43] Kristofer Bouchard and Ankit Kumar. Feedback controllability is normative theory of neural population dynamics. Research Square, 29 March 2024.
[44] Urs Braun, Anais Hareiit, Giulio Pergola, Tommaso Menara, Axel Schäfer, Richard F Betzel, Zhenxiang Zang. Janina Schweiger, Xiaolong Zhang, Kristina Schwarz, Junfang Chen, Giuseppe Blasi, Alessandro Bertolino, Daniel Durstewitz, Fabio Pasqualetti, Emanuel Schwarz, Andreas Meyer-Lindenberg Danielle S Bassett, and Heike Tost. Brain network dynamics during working memory are modulated by dopamine and diminished in schizophrenia. Nat. Commun., 12(1):3478, 9 June 2021.
[45] Dale Zhou, Yoona Kang, Danielle Cosme, Mia Jovanova, Xiaosong He, Arun Mahadevan, Jeesung Ahn, Ovidia Stanoi, Julia K Brynildsen, Nicole Cooper, Eli J Comnblath, Linden Parkes, Peter Mucha, Kevin N Ochsner, David M Lydon-Staley, Emily Falk, and Dani Bassett. Mindful attention promotes control of brain network dynamics for self-regulation and discontinues the past from the present. Proc. Natl, Acad. Sei. U. S. A., 1200)):2201074119, 10 January 2023.
[46] Daniela Zöller, Corrado Sandini, Marie Schaer, Stephan Eliez, Danielle S Bassett, and Dimitri Van De Ville. Structural control energy of resting-state functional brain states reveals less cost-effective brain dynamics in psychosis vulnerability. Hum. Brain Mapp., 42(7):2181--200, May 2021.
[47] Xiaosong He, Lorenzo Caciagli, Linden Parkes, Jennifer Stiso, Teresa M Karrer, Jason Z Kim, Zhixin Lu, Tommaso Menara, Fabio Pasqualetti, Michael R Sperling, Joseph Tracy, and Dani S Bassett. Uncovering the biological basis of control energy: Structural and metabolic correlates of energy inefficiency in temporal lobe epilepsy. Sei. Adv, 8(45)eab, 11 November 2022.
13

---

[48] S Parker Singleton, Andrea I Luppi, Robin L Carhart-Harriss Josephine Cruzat, Leor Roseman, David Nutt, Gustavo Deco, Morten L Kringelbach, Emmanuel Stamatakis, and Amy Kuceyeski. Receptor-informed network control theory links LSD and psilocybin to a flattening of the brain's 's control energy landscape. Nat. Commun., 13(1):1-13,3 October 2022.
[49] John D Medaglia, Denise Y Harvey, Nicole White, Apoorva Kelkar, Jared Zimmerman, Danielle Bassett, and Roy H Hamilton. Network controllability in the inferior frontal gyrus relates to controlled language variability and susceptibility to TMS. J. Neurosci., 38(28): 6399-6410, 11 July 2018.
[50] Zhoukang Wu, Liangjiccheng Huang, Min Wang, and Xiaosong Hc. Development of the brain network control theory and its implications. Psychoradiology, 4:kkae028, 14 December 2024.
[51] Weidong Cai, Srikanth Ryali, Ramkrishna Pasumarthy, Viswanath Talasila, and Vinod Menon. Dynamic causal brain circuits during working memory and their functional controllability. Nat. Commun., 12(1):3314, 29 June 2021.
[52] Ali Moradi Amani, Amirhessam Tahmassebi, Andreas Stadlbauer, Uwe Meyer- Meyer-Baese, Vincent Noblet, Frederic Blanc, Hagen Malberg, and Anke Meyer-Baese Controllability of functional and structural brain networks. Complexity, 2024(1):7402894, January 2024.
[53] Qian Li, Li Yao, Wanfang You, Jiang Liu, Shikuang Deng, Bin Li, Lekai Luo, Youjin Zhao, Yuxia Wang. Yaxuan Wang, Qian Zhang, Fenghua Long, John Sweeney, Shi Gu, Fei Li, and Qiyong Gong, Controllability of functional brain networks and its clinical significance in first-episode schizophrenia. Schizophr. Bull., 49(3)665--668, 3 May 2023.
[54] Yang- Yu Liu, Jean-Jacques Slotine, and Albert-Lászile Barabási. Controllability of complex networks. Nature, 47377346)1177-173 12 May 2011.
[55] Linden Parkes, Jason Z. Kim, Jennifer Stiso, Julia K Brynildsen, Matthew Cieslak, Sydney Covitz, Raquel E Gur, Ruben C Gur, Fabio Pasqualetti, Russell T Shinohara, Dale Zhou, Theodore D Satterthwaite, and Dani Bassett. A network control theory pipeline for studying the dynamics of the structural connectome. Nat. Protoc., 19(12):321 29 December 2024.
[56] David R. Tyner and Andrew D. Lewis. Geometric jacobian linearization and LQR theory. J. Geom, Mech., 2(4):397-440, 2010.
[57] Robert E Kass, Heejong Bong, Motolani Olarinre, Qi Xin, and Konrad N Urban. Identification interacting neural populations: methods and statistical considerations. Neurophysiol, 130 (3):455-966, September 2023.
[58] Matthew G Perich, Charlotte Arlt, Sofia Soares, Megan E Young, Clayton P Mosher, Juri Minxha, Eugene Carter, Ueli Rutishauser, Peter H Rudebeck, Christopher D Harvey, and Kanaka Rajan. Inferring brain-wide interactions using data-constrained recurrent neural network models. bioRxiv, page 1020.12.18.423333...
11 March 2021.
[59] Matthew G Perich and Kanaka Rajan. Rethinking brain-wide interactions through multi-region 'network of networks' models. Curr. Opin. Neurobiol., 15:146-151, December 2020.
[60] Aaron S Andalman, Vanessa M Burns, Matthew Lovett-Barrow Michael Broxton, Ben Poole, Samuel J Yang, Logan Grosenick, Talia N Lemner, Ritchie Chen, Tyler Benster, Philippe Mourrain, Mare Levoy, Kanaka Rajan, and Karl Deisseroth. Neuronal dynamics regulating brain and behavioral state transitions. Cell, 177(4):970-988 0077009985222... 2 May 2019.
[61] Lucas Pinto, Kanaka Rajan, Brian DePasquale, Stephan Y Thiberge, David Tank, and Carlos D Brody. Task-dependent changes in the large-scale dynamics and necessity of cortical regions. Neuron, (04(4):810884400 20 November 2019.
[62] Michael Kleinman, Chandramouli Chandrasekaran, and Jonathan Kao. A mechanistic multi- area recurrent network model of decision-making. Adv. Neural Inf. Process. Syst. 34:23152-
23 23165, 6 December 2021.
14

---

[63] Joao Barbosa, Rémi Proville, Chris C Rodgers, Michael R DeWeese, Srdjan Ostojic, and Yves Boubenee. Early selection of task-relevant features through population gating. Nat. Commun. 14(1): 27 October 2023.
[64] Evren Gokcen, Anna I Jasper. João D Semedo, Amin Zandvakili, Adam Kohn, Christian K Machens, and Byron M Yu. Disentangling the flow of signals between populations of neurons. Nature Computational Science, 2(8):512- 525, 18 August 2022.
[65] Evren Gokcen, Anna Ivic Jasper, Alison Xu, Adam Kohn, Christian K Machens, and Byron M Yu. Uncovering motifs of concurrent signaling across multiple neuronal populations. In Thirty-seventh Conference on Neural Information Processing Systems, 2 November 2023.
[66] Sadegh Ebrahimi, Jéröme Lecoq, Oleg Rumyantsev, Tugce Tasci, Yanping Zhang, Cristina Irimia, Jane Li, Surya Ganguli, and Mark J Schnitzer. Emergent reliability in sensory cortical coding and inter- area communication. Nature, 605(7911):713-72,, May 2022.
[67] Jordan Rodu, Natalie Klein, Scott L Brincat, Earl K Miller, and Robert E Kass. Detecting multivariate cross-correlation between brain regions. J. Neurophysiol. 12044):162--1972, October 2018
[68] João D Semedo, Anna Jasper, Amin Zandvakili. Aravind Krishna, Amir Aschner, Christian K Machens, Adam Kohn, and Byron M Yu. Feedforward and feedback interactions between visual cortical areas use different population activity patterns. Nat. Commun., 13(1):1099, March 2022.
[69] Hao Ye, Ethan Deyle, Luis J Gilarranz, and George Sugihara. Distinguishing time-delayed causal interactions using convergent cross mapping. Sci. Rep., 5:14750,5 October 2015.
[70] Satohiro Tajima, Toru Yanagawa, Naotaka Fujii, and Taro Toyoizumi. Untangling wide dynamics in consciousness by cross-embedding. PLoS Comput. Biol., 11014): November 2015.
[71] George Sugihara, Robert May, Hao Ye, Chih-Hao Hsieh, Ethan Deyle, Michael Fogarty, and Stephan Munch. Detecting causality in complex ecosystems. Science, 33866166) 26 October 2012.
[72] Joshua I Glaser, Matthew Whiteway, John P Cunningham, Liam Paninski, and Scott W Linderman. Recurrent switching dynamical systems models for multiple interacting neural populations. In Advances in Neural Information Processing Systems 33 (NeurIPS 2020), page 2020. 1e..1......99... 22 October 2020.
[73] Orren Tamboor,, David M Zoltowski, E Mika Diamanti, Lucas Pinto, David W Tank, Carlos D Brody, and Jonathan W Pillow. Modeling communication and switching nonlinear dynamics in multi-region neural activity. bioRxiv, page 2022.
15 September
[74] Steven L Bressler and Anil K Seth. Wiener Wiener-granger causality: a well established methodology. Neuroimage, 58(2)323-322, 15 September 2011.
[75] Anil K Seth, Adam B Barrett, and Lionel Barnett. Granger causality analysis in neuroscience and neuroimaging. J. Neurosci, 3583)3323-3327, 25 February 2015.
[76] K J Friston, L Harrison, and W Penny. Dynamic causal modelling. Neuroimage, 19(4): 1273-1302, August 2003.
[77] Yu Chen, Hannah Douglas, Bryan J Medina, Motolani Olarinre, Joshua H Siegle, and Robert E Kass. Population burst propagation across interacting areas of the brain. J. Neurophysiol. 128 (6):157-11992, December 2022.
[78] Josuan Calderon and Gordon J Berman. Inferring the time-varying coupling of dynamical systems with temporal convolutional autoencoders. eLife, 12 November 2024.
[79] Ziyu Lu, Anika Tabassum, Shruti Kulkarni, Lu Mi, J Nathan Kutz, Eric Shea-Brown and Seung-Hwan Lim. Attention for causal relationship discovery from biological neural dynamics. In NeurIPS 2023 Workshop on Causal Representation Learning, 12 November 2023.
15

---

[80] Jean-Jacques E Slotine and Weiping Li. Applied Nonlinear Control. Prentice-Hii Hall, 1991.
[81] R W Brockett. Nonlinear systems and differential geometry. Proc. IEEE Inst. Electr. Electron. Eng., 64(1):61-72, 1976.
[82] G W Haynes and Hermes. Nonlinear controllability via lie theory. SIAM J. Control, 8(4):
450 460, November 1970.
[83] Cesar O Aguilar and Bahman Gharesifard. Necessary conditions for controllability of nonlinear networked control systems. In 2014 American Control Conference, pages 33733 IEEE, June 2014.
[84] Andrew J Whalen, Sean N Brennan, Timothy D Sauer, and Steven J Schiff. Observability and controllability of nonlinear networks: The role of symmetry. Phys. Rev. X., 5(1):01 23 January 2015.
[85] Le-Zhi Wang. Ri-Qi Su, Zi-Gang Huang, Xiao Wang, Wen-Xu Wang, Celso Grebogi, and Ying-Cheng Lai. A geometrical approach to control and controllability of nonlinear dynamical networks. Nat. Commun., 7:1 1323, 14 April 2016.
[86] Yang Tang, Huijun Gao, Wei Zou, and Jürgen Kurths. Identifying controlling nodes in neuronal networks in different scales. PLoS One, 7(7):411375, 27 July 2012.
[87] Gautam Kumar, Delsin Menolascino, and Shinung Ching. Input novelty as a control metric for time varying linear systems. arXiv [math. OC], 21 November 2014.
[88] Ethan R Deyle, Robert M May. Stephan B Munch, and George Sugihara. Tracking and forecasting ecosystem interactions in real time. Proc. Biol. Sci., 283(18 13 January 2016.
[89] Tosif Ahamed, Antonio C Costa, and Greg Stephens. Capturing the continuous complexity of behaviour in caenorhabditis elegans. Nat. Phys., 17(2):275 283, 5 October 2020.
[90] Frédéric Latrémolière, Sadananda Narayanappa, and Petr Vojtěchovský. Estimating the jacobian matrix of an unknown multivariate function from sample values by means of a neural network. arXiv LC), April 2022.
[91] Jonathan Lorraine and Safwan Hossain. Jacnet: Learning functions with structured jacobians. arXiv preprint 24ti::2200..3.3337 2024.
[92] Hadi Beik-Mohammadi, Søren Hauberg, Georgios Arvanitidis, Nadia Figueroa, Gerhard Neumann, and Leonel Rozo. Neural contractive dynamical systems. arXiv [cs. 17 January 2024.
[93] Winfried Lohmiller and Jean-Jacques E Slotine. On contraction analysis for non non-linear systems. Automatica, 34(6))683-6)) 696, June 1998.
[94] Florian Hess, Zahra Monfared, Manuel Brenner, and Daniel Durstewitz. Generalized teacher forcing for learning chaotic dynamics. arXiv [cs.LG], 7 June 2023.
[95] Abhiram Iyer, Sarthak Chandra, Sugandha Sharma, and I Fiete. Flexible mapping of abstract domains by grid cells via self-supervised extraction and projection of generalized velocity signals. Neural Inf Process Syst, 37:85441-85444 2024.
[96] Balth van der Pol. LXXXVIII. on dellaxation-oscllaaiionn Dublin Philosophical Magazine and Journal of Science, 2(11):978-99) 9989,, November 1926.
The London. Edinburgh, and
[97] Edward N Lorenz. Deterministic nonperiodic flow. J. Atmos. Sci., 30(2):130-14) 1 March 1963.
[98] E N Lorenz. Predictability: problem partly solved. In ECMWF Seminar on Predictability 4-8 September 1995. European Centre for Medium Medium-Range Weather Forecasts, 1996.
[99] William Gilpin. Chaos an interpretable benchmark for forecasting and data-driven modelling. In Vanschoren and S Yeung, editors, Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks, volume 1, 2021.
16

---

[100] William Gilpin. Model scale versus domain knowledge in statistical forecasting of chaotic systems, Phys. Rev. Res., 5(4):043252, 15 December 2023.
[101] Pablo Gómez, Håvard Hem Toftevaag, and Gabriele Meoni. torchquad: Numerical Integration in Arbitrary Dimensions with PyTorch. Journal of Open Source Software, 64(6), 2021. doi: 10.21 1i.21105//oos.03433
[102] Ricky Q Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In S Bengio, H Wallach, H Larochelle, K Grauman, N R Garnett, editors, Advances Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
and
[103] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Gimelshein, L Antiga, Alban Desmaison, Andreas Köpf, E Yang. Zach Vito,o Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An imperative style, high-performance deep learning library. Adv. Neural Inf. Process. Syst., abs/1912.01 3 December 2019.
[104] Judy Hoffman, Daniel A Roberts, and Sho Yaida. Robust learning with jacobian regularization. arXiv [stat.MLJ, 7 August 2019.
[105] Alexander Wikner, Joseph Harvey, Michelle Girvan, Brian R Hunt, Andrew Pomerance, Thomas Antonsen, and Edward Ott. Stabilizing machine learning prediction of dynamies: Novel noise-inspired regularization tested with reservoir computing. Neural Netw, 170: 94-110, 1 February 2024.
[106] Steffen Schneider, Rodrigo González Laiz, Anastasiaa Filippova, Markus Frey, and Macken- zic Weygandt Mathis. Time-series attribution maps with regularized contrastive learning. arXiv [s], 17 February 2025.
[107] Panos J Antsaklis and Anthony N Michel. Linear Systems. Birkhäuser Boston, 2006.
[108] Yu Kawano and Jacquelien M Scherpen. Empirical differential gramians for nonlinear model reduction. Automatica (Oxf.), 127(109534)1109534 May 2021.
[109] Gustav Lindmark and Claudio Altafini. Minimum energy control for complex networks. Sci. Rep., 8(1):3188 16 February 2018.
[110] Matthew F Panichello and Timothy J Buschman. Shared mechanisms underlie the control of working memory and attention. Nature, March 2021.
[111] Weiwei Li and Emanuel Todorov. Iterative linear quadratic regulator design for nonlinear bio- logical movement systems. ICINCO 2004, Proceedings of the First International Conference on Informatics Control, Automation and Robotics, Setábal, Portugal, August 25-28, 2004, 222-229 1 January 2004.
[112] Yuval Tassa, Tom Erez, and Emanuel Todorov. Synthesis and stabilization of complex behaviors through online trajectory optimization. In 2012 IEEE/RSJ International Conference Intelligent Robots and Systems, pages 4906-4913 IEEE, October 2012.
[113] S Ouala, D Nguyen, L Drumetz, B Chapron, A Pascual, F Collard, L Gaultier, and R Fablet. Learning latent dynamics for partially observed chaotic systems. Chaos, 30(10):00312, 20 October 2020.
[114] William Gilpin. Deep reconstruction of strange attractors from time series. In Proceedings of the 34th International Conference on Neural Information Processing Systems, number Article 18 in NIPS'20, pages 204-216, Red Hook, NY, USA, December 2020. Curran Associates Inc.
[115] Peter Y Lu, Joan Ariño Bernad, and Marin Soljačić. Discovering sparse interpretable dynamics from partial observations. Commun. Phys., 5(1):1-7, 12 August 2022.
[116] Elise Özalp, Georgios Margazoglou, and Luca Magri. Reconstruction, forecasting, and stability of chaotic dynamics from partial data. Chaos, 3339), September 2023.
17

---

[117] Charles D Young and Michael D Graham. Deep learning delay coordinate dynamics for chaotic attractors from partial observable data. Phys. Rev. E., 107(3-1):032215 30 March 2023.
[118] Sean Jaffe, Alexander Davydov, Deniz Lapsekili, Ambuj Singh, and Francesco Bullo. Learning neural contracting dynamics: Extended linearization and global guarantees. arXiv [cs.L.G) 12 February 2024.
[119] Steven H Strogatz. Nonlinear Dynamics and Chaos: With Applications Physics, Biology, Chemistry, and Engineering. Avalon Publishing, 29 July 2014.
[120] Rudiger U Seydel. Practical bifurcation and stability analysis. Interdisciplinary applied mathematics. Springer, New York, NY, edition, 10 December 2009.
[121] Luca Dieci, Robert D Russell, and Erik S Van Vleck. On the computation of Yуарuпo exponents for continuous dynamical systems. SIAM J. Numer. Anal., 34(1):402-423, 1997.
[122] F Christiansen and H H Rugh. Computing lyapunov spectra with continuous gram schmidt orthonormalization. Nonlinearity, 10(5):1063, September 1997.
[123] Vincent Duchaine and Clement M Gosselin. Investigation of human-robot interaction stability using дуарunov theory. In 2008 IEEE International Conference on Robotics and Automation, pages 2189-2194. IEEE, May 2008.
[124] Ian R Manchester and Jean-Jacques E Slotine Control contraction metries: Convex and intrinsic criteria for nonlinear feedback design. IEEE Trans. Automat. Contr., 62(6))3046- 3053, June 2017.
[125] Patrick M Wensing and Jean-Jacques Slotine. Beyond convexity-contraction and global convergence of gradient descent. PLoS One, 158)):0036661, 4 August 2020.
[126] Leo Kozachkov, Mikael Lundqvist, Jean Jacques Slotine, and Earl Miller. Achieving stable dynamics in neural circuits. PLoS Comput. Biol., 16(8):1-15, 2020.
[127] N M. Hootsmans and Dubowsky. Large motion control of mobile manipulators including vehicle suspension characteristics. In Proceedings. 1991 IEEE International Conference Robotics and Automation, pages 2336-2341 vol.3. IEEE Comput. Soc. Press, 2002.
[128] S Ali A Moosavian and Evangelos Papadopoulos. Modified transpose jacobian control of robotic systems. Automatica (Oxf.), 43(7):1226-1233, July 2007.
[129] Krzysztof Tchoń, Adam Ratajezak, and Ida Góral. Lagrangian jacobian inverse for nonholo- nomic robotic systems. Nonlinear Dyn., 82(4):1923-19) 1932, December 2015.
[130] Chevelle Brudey, Jeanie Park, Jan Wiaderkiewicz. Ihori Kobayashi, Thomas A Mellman, and Paul J Marvar. Autonomic and inflammatory consequences of posttraumatic stress disorder and the link cardiovaseular disease. Am. J. Physiol. Regul. Inlegr. Comp. Physiol. 309(4): R315-21, 15 August 2015.
[131] Gail Alvares, Daniel S Quintana, Ian B Hickie, and Adam Guastella. Autonomic nervous system dysfunction in psychiatric disorders and the impact of psychotropic medications: systematic review and meta anta--naalysi. J. Psychiatry Neurosci., 41(2):89-104, March 2016.
[132] Jeffrey J Goldberger, Rishi Arora, Una Buckley, and Kalyanam Shivkumar. Autonomic nervous system dysfunction: JACC focus seminar. Am. Coll. Cardiol.. 73(10):1189-13) 19 March 2019.
[133] Li Xiong and Thomas W H Leung. Autonomic dysfunction in neurological disorders. Aging (Albany NY), 19044 9 April 2019.
[134] Albert-L.ázz16 Barabási, Natali Gulbahce, and Joseph Loscalzo. Network medicine: a network- based approach to human discase. Nat. Rev. Genet., 12(1):56-68, January 2011.
18

---

[135] Matthew T Maurano, Richard Humbert, Eric Rynes, Robert E Thurman, Eric Haugen, Hao Wang, Alex P Reynolds. Richard Sandstrom, Hongzhu Qu. Jennifer Brody, Anthony Shafer, Fidencio Neri, Kristen Lee, Tanya Kutyavin, Sandra Stehling-Sun, Audra K Johnson, Theresa K Canfield, Erika Giste, Morgan Diegel, Daniel Bates, R Scott Hansen, Shane Neph, Peter Sabo, Shelly Heimfeld, Antony Raubitschek, Steven Ziegler, Chris Cotsapas, Nona Sotoodehnia, Ian Glass, Shamil R Sunyaev, Rajinder Kaul, and John Stamatoyannopoulos. Systematic localization of common discase-associatee variation in regulatory DNA. Science, 337(6099); 1190-1195, September 2012.
[136] Piyush B Madhamshettiwar, Stefan R Maetschke, Melissa J Davis, Antonio Reverter, and Mark A Ragan. Gene regulatory network inference: evaluation and application ovarian cancer allows the prioritization of drug targets. Genome Med., 4(5):41, , May 2012.
[137] Tong Ihn Lee and Richard A Young. Transcriptional regulation and its misregulation in disease. Cell, 153(6))1377-125) 14 March 2013.
[138] Paula Unger Avila, Tsimafei Padvitski, Ana Carolina Leote, He Chen, Julio Saez-Rodriguez, Martin Kann, and Andreas Beyer. Gene regulatory networks in disease and ageing. Nat. Rev. Nephrol., 20(9):616 616-633, September 2024.
[139] Jaideep Pathak, Zhixin Lu, Brian R Hunt, Michelle Girvan, and Edward Ott. Using machine learning replicate chaotic attractors and calculate lyapunov exponents from data. Chaos, 27 (12):121102, December 2017.
[140] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott. Model-free predic- tion of large spatiotemporally chaotic systems from data: reservoir computing approach. Phys. Rev. Lett. 120(2):024102, 12 January 2018.
[141] Zhixin Lu, Jaideep Pathak, Brian Hunt, Michelle Girvan, Roger Brockett, and Edward Ott. Reservoir observers: Model-free inference of unmeasured variables in chaotic systems. Chaos, 77(4):04112,, 5 April 2017.
[142] Amitava Banerjee, Jaideep Pathak, Rajarshi Roy, Juan G Restrepo, and Edward Ott. Using machine learning to assess short term causal dependence and infer network links. Chaos, 29 (12):121104, 26 December 2019.
[143] Mattia Rigotti, Omri Barak, Melissa R Warden, Xiao-Jing Wang, Nathaniel D) Daw, Earl K Miller, and Stefano Fusi. The importance of mixed selectivity in complex cognitive tasks. Nature, (97(4441)7555555 (55555555)5 30 May 2013.
19

---

A JacobianODE technical details
A.1 Parameterizing the time derivative via the Jacobian
Recall that, for the time derivative function f, we have that the path integral is independent of the choice of path:
f(x(tt)) f(x(t()) = 1 J ds = S. f. J(c(r)) dr,
(7)
where C a piecewise smooth curve in TK" and c: [1,1,] + C parameterization of with c(t+) x(4)) and c(ty) x(ty) (also see equation 2). Path integration notably provides only differences between distinct points. Thus, knowing the value of one point needed for trajectory reconstruction, which we do not assume. To address this, we note it is possible to represent solely through the Jacobian at a point x(t) as
f(x(t)) G(to. t; J.c) + x(t) x(to)
(8)
1 4-40
where
((0,,4JJ,0) [
(9)
and we have abbreviated Cr,(()) (,,;,,,,x())))) x()),, piecewise smooth curve on [s, (] parame- terized by r and beginning and ending at x(s) and x(t) respectively. Intuitively, we can circumvent the need to know f by recognizing that integrating f between time points will produce the difference in the system states at those time points, which are known to us. With this formalism, we avoid the need to represent directly, thereby enabling all gradients backpropagate through the Jacobian network. The proof of this formalism is presented below.
Proposition 1 caaoobiin-pprammttriieee ODEs). Let xt) = f(x(t)) and let Jx(x(t)) = J(x()) =
8x fxx(t)). Then given times to, t we can express f parameterized by the Jacobian as
f(x(1)) (( x(1) (()))
t t- to
where
G(,o = f(' J(ca:(rr)),
and we have abbreviated (T) c(r; s,,, x(t)), piecewise smooth curve on [s, t) parame- terized by and beginning and ending at (() and x(t) respectively.
Proof. For given times I, lo. and from the fundamental theorem of calculus, we have that
x(t) x(to) / f(x(s))ds
and
f(x(t)) - f(x(s)) = [ (1
Letting
H(s, = S: J(C f(x(I) ((x()))((x(()
We then have that
20

---

S (ssttdds [ f(x(1)) f(x(s))ds
f f(x(1))ds /'' f(x(s))ds
= (1-11))(((()) - -x((1) -((1))
Letting now
(,) S H(s,))ds -f" [. ((c.,(t))
We can see that
G(lo, (;J,c) -(1-1))((x(0)) (x(1) x(1n))
and thus
f(x(t)) G(to, t; x(t) - x(fo)
t to
A.2 Path integrating to generate predictions
Given an initial observed trajectory x(to - kdtt,, k = o,.... of length at least two (6 2 1), we compute an estimate of f(x(4o + bat)) of f(x(1o + ba()) by replacing with jé in equation 8 Specifically, we compute an estimate of f(x((t)) of f(x(ta)
f(x(tt)) G( to: Jo c) + x(to) x(to)
(10)
tb to
In practice, compute this estimate by constructing cubic spline on the initial observed trajectory using 15 points (i.e., = 14). Given that the computation G involves a double integral one over states and over time- using spline is computationally advantageous. This is because the path integral (and all intermediate steps) can be quickly computed along the full spline, with the result of each intermediate step along the path then being summed to approximate the time integral. The integral is computed by interpolating points for every gap between observed points, resulting in discretization of 58 points along the spline. Integrals are computed using the trapezoid method from torchquad [01].
Once we have constructed our estimate of f(x(In) we can estimate f at any other point x(!) as
H(to, + f(x(t))), f(x(te)) tt,,
if
f(x(t))
if if =t
f(x((t))
where by convention we integrate forwards in time and H is the path integral defined above, with J in place of J. In practice, for the integration path c(r; 8, x(s), x()), we construct line from x(s) x(t) as
c(r; ((s),, -1- x(
-- --xx(0x x(!)
to maintain the interpretability of having r in the range [s, t however can be easily seen that setting r' we recover the familiar line
c(r') = (1 (r-r)()(()) +r'x(t)
21

---

with r' taking values on [0. 1]. Line integrals are computed with 20 discretization steps using the trapezoid method from torchquad [TOT]
Using f(x(t)), we can then generate predictions
x(t + At) x(t)
f(x(F))dr
(11)
where the integral can be computed by standard ODE solver (we used the RK4 method from torchdiffeg with default values for the relative and absolute tolerance).
We can then compute the trajectory reconstruction loss as
Luai (0: x)
Σ
d (x(1o + &DI), x(oo + kat))
(12)
k=b+1
where d is a distance measure between trajectories (e.g. mean squared error).
A.3 Loop closure loss
While the space Jacobian functions J that solve the dynamics problem may be quite large, we are interested only in functions J" that both solve the dynamics problem and generate matrices with that conservative vector fields. To effectively restrict optimization this space, employ a self-supervised loss building on the loss introduced by lyer et al. (95]. This loss imposes a conservative constraint the of the Jacobian matrix. Specifically, for any piecewise smooth
curve Choop starting and ending at the same state Xo. we have that
[Gm. Jds
0.
Thus, given Thoups sets of any L (not necessarily sequential) points (x (f1) x (t2). x x(' (t2))
we can form loop Gioop ( consisting of the sequence lines from (t.) to x (l++1), 1), 1, L,.., 1-1,
followed by a line from x (tz) to x(4) (t1). We then define the following regularization loss to
enforce this conservation constraint:
Thons
Cioop (hop(0::3 x) =
js
(13)
Moooos
( keep
We note that while other choices for loops would have been possible (for instance, forward and backward passes along observed trajectories, or arbitrary circles beginning and ending at the same point), the approach presented here has several advantages (as discussed in Section 3.3). Namely, balances computational tractability with uniform sampling of the directions the tangent space. Enforcing that all the points that comprise the loop are on the data manifold ensures that the loops are integrated in directions that are informative as possible for estimating the Jacobian of the given dynamical system.
B
Controc- analysis details
B.1
Gramian computation
We begin with equation 4 from Section 5) in which we separate the locally-linearized dynamics in the tangent space into separate pairwise inter-subsystem control interactions. These pairwise control interactions are in general of the form
8x (x (t)
A(t)xx B(t)6x²
where A is the within-subsystem Jacobian, B is the across subsystem Jacobian and x xB are the states of subsystems A and B, respectively. For a given control system, the vireinarying reachability Gramian on the interval (o. 4] is defined as
2 (T)) 1" t0
Wr (2) e
22

---

where o denotes the state-transition matrix of the intrinsic dynamics of subsystem A without any imput (i.e., (t) (tt, to (t)) [07]] Differentiating with respect t, we obtain 5x (t) (to). Noting also that, the absence of input from subsystem B, we have
(x) (t) (tt)xx (o,) (fo)
Thus setting the equations equal to each other and canceling 8x (to) from both sides we obtain
0 (t. (o) -A (t)
Note that ((tu to) = I. Now, letting (t,,r) 7 (T) 7t and differentiating the Gramian expression with respect to the second argument (and using the Leibniz integral rule) yields
(1o,,, = 0
T(t, T)dT
Ot
0t
0 -I(,,))
a (1)) + S
a 7(,
(1()) r(l,
0t
(B(()B((((1) 0+ l"d 01 r(,,))d
and observing
( T
T(,,T) = (01
0t
= A(t) (1)
= +())
can continue
i Ot
W, Wr(1o, B(I)B(() (()) / T(t, + (f((()))ds (1) = =B( + A(t) (to, t) + W, (t)
which illustrates that the reachability Gramian can be solved for using an ODE integrator (with initial condition W.(4o, (n) 0) [56]. In practice, compute the reachability Gramians using the trained JacobianODE models, we fit a cubic spline c(t) to the reference trajectory x(t) and compute J(t) (c(t). We can then parse the Jacobian matrix into its component submatrices and compute the Gramians accordingly.
The reachability Gramian is symmetric positive semidefinite matrix [107] The optimal cost of driving the system from state Jxy to state 5x1 on time interval [to, t1] can be computed as
(6x1 ((1,40))x0 T w (lo, 41) (6x1 < (o)6x))
Thus each eigenvalue W, reflects the case of control along the corresponding eigenvector [107 [09]. Eigenvectors with larger corresponding eigenvalues will have smaller inverse eigenvalues, and thus scale the cost down along those directions when computing the cost as above.
B.2 Extension to non-autonomous dynamics
While we deal with autonomous systems in this work, we note that this construction can easily be
extended to include non-autonomous dynamical systems, of the form x(f) = f(x, E) by constructing an augmented state E g=+1 which simply the concatenation of x with t. This yields autonomous dynamics x = f(x), where f IR=+1 - Ru+1 is the concatenation of with a function that maps all states to 1.
23

---

B.3 Extension to more than two subsystems
Consider a nonlinear dynamical system in IR,, defined by x(t) f(x(t)). Suppose now that the system is composed of K subsystems, where the time evolution of subsystem & is given by x (t) "R x(t) is then comprised of concatenation of the x (t), with Σ Th Given a particular reference trajectory, x(t), with as the Jacobian of f. then the tangent space dynamics around the reference trajectory for subsystem {I,..., K} are given by
(x
(1) Σ J (x))))xx(x(()
where (x=(((((t) E Reax Xn, is the submatrix of J(x(t)) in which the columns correspond to subsystem i and the rows correspond subsystem a. Now, for a given subsystem 3 E {l, K} with 3 * a, we wish to analyze the ease with which 3 can control a locally around the reference trajectory, without intervention from other subsystems. Discounting the interventions from other subsystems equates to setting ox (t) = 0 for k # a, B, leaving the expression
(x(() (1) = -x (x(1)) ( (1) + + (x())xxxxx
which quantifies the influence / exert subsystem. Using this representation, the reachability Gramian can be computed as described above. Performing this procedure for all pairs of subsystems a, B E {I, K} thus characterizes all pairwise control relationships between subsystems along the reference trajectory.
in the absence of perturbations from any other
C
Supplementary results
C.1 Control and communication capture different phenomena
In the Introduction (Section I we note that measuring communication between two systems is different than measuring control. To illustrate this, consider two brain areas, and B, whose dynamics are internally linear. The areas are connected only through a linear feedforward interaction term from area area A (Figure STA). Concretely, we consider the dynamics:
x (t) = J A (t) J
x (1) = J (()
Here, the J matrices are time- inverinvariant. When considering communication between brain areas, one might aim to find the subspaces in which communication between B and A occurs, as well as the messages passed 33 73). In this construction, the subspace in which area B communicates with area A is explicitly given by J B-+A, and thus the message, or input, from area B to area
any time t is simply given by J - (t(). We pick the dimensions of each region to be nA = B = and let the eigenvectors of the (negative definite matrix) JA-AA be given by v,, i 1,2,3, 4. The eigenvectors are numbered in order of decreasing real part of their corresponding eigenvalue. Now, suppose we construct the interaction matrix J in different ways: If we let (1) Jf v, 1 (Figure SIB. left), then the signal from B is projected onto the most stable mode of region A, whereas if we let (2) B BAA v,1" (Figure SIB. right), the signal projected onto the least stable mode. While interaction 1 and interaction 2 communicated messages with identical magnitudes, interaction 2 led much lower cost reachability control (Figure STC). This illustrates that control depends not only on the directions along which the areas can communicate, but also on how aligned the communication is with the target area's 's dynamics.
C.2 Derivative estimation not implied by function estimation
An alternative approach to directly estimating the Jacobians would be to learn an approximation f of the function f, and then approximate the Jacobian via automatic differentiation (i.e., an esti- mate j 8 f, as with the NeuralODEss) While this can be effective in certain scenarios, it is not generally the case that approximating function well will yield a good approximation of its derivative. To illustrate this, recall example from Latrémolière et al. [90), in which functions
24

---

A
B
Interaction
Interaction
X
A
JA-A
-
x
A
1
J
X
vA-A
Vipa
1,2
Interaction 105 "
Interastio (f
D
[n(x) 1 cos(per)
(2(2) sin(
10.0
bi
coosinx)
sin(mx)
Reachate by esse
ricattriccioonn strenghh
(0
-20
20
-20
20
Trm 30
Figure S1: Communication versus control in linear systems and the challenge Jacobian estimation. (A) Setup of two linearly connected brain areas, and B. (B) Interaction matrices projecting signals from area onto either the most stable (Interaction 1) or least stable (Interaction 2) eigenvectors of area A. (C) Although both interactions communicate identical signal magnitudes, Interaction 2 provides significantly enhanced reachability due to alignment with unstable modes of area dynamies. (D) An illustrative example demonstrating that accurate approximation of function (left panel) does not guarantee accurate approximation of its derivative (right panel), emphasizing the necessity of directly estimating the Jacobian.
fa = cos(nr) nx) are used to approximate the function f(x) = 0 (Figure SID). While these approx- imations improve with increasing (i.e., limm-+oe fn = f). this is not the case for the derivative (lim,- f/ sin (nz) # f'). This demonstrates the necessity of learning J directly (rather than first approximating f), which we also demonstrate empirically. this setting could be interpreted as overfitting 190]. As long as function estimates match at the specific points in the training set, how the function fluctuates between these points not constrained match the true function. For this reason, included the Frobenius norm Jacobian regularization in our implementation of the NeuralODEs (Appendix D.3).
the context of machine learning,
C.3 Full benchmark dynamical systems results
We here present the full results for all considered example dynamical systems. 10 time-step trajectory predictions along with Jacobian estimation and estimated Lyapumov spectra are displayed in Figure along with MSE and R in Table SSI (trajectory prediction) and Table S2 (Jacobian estimation). Jacobian estimation errors using the 2-norm are presented in Table S3
S2
C.4 Ablation studies
To determine the value of the different components of the JacobianODE learning framework, performed several ablation studies. We chose to evaluate on the Lorenz system and the task-trained RNNs, as these together provide two common settings of chaos and stability, well low- and high-dimensional dynamics.
Ablating the Jacobian- parameterized ODEs The JacobianODE framework constructs an estimate of the initial time derivative via double integral of the Jacobian as described in Section 3.2 To briefly recall, Jacobian path integration is described as
f(x(L)) ((x(t)) = [ J ((((((()))( dr,
(14)
When we do not have access to the initial derivative and base point f(x(1.)), x(1:), we use the following formulation
25

---

10 time-step
System
n
State space
trajectory prediction
Jacobian estimation
Lyapunov spectrum
MSE
0.005
MBR
10
0.000
Expement
Van der Pol
1
1%
5%
10%
Observation noise
Observation
Iyspunov
MSE
0.25
0.00
MISE
50
Lorenz
4
1%
5%
10%
1%
-40
Observation noise
Observation noise
C
MSE
0.00
Lorenz 96
12
1%
5%
10%
115
6%
10%
Obsarvation noise
Obeervation nolee
MSE
0.25
MSE
2
0.00
Lorenz 96
32
1%
5%
10%
10%
Obsarvation noise
Oeservation
Lyspunov exponent
MSE
Expoient
Lorenz 96
64
0.5
0.0
1%
10%
5%
Observation noise
Observation naise
Lyspunov exponent
JacobianODE
NeuralODE
Weighted linear Jacobians
Ground truth system
Figure S2: Full dynamical systems prediction results. State space representations (using spectral embedding for systems that have more than dimensions). 10 time-step trajectory predictions, Jacobian estimation, and Lyapunov spectrum estimates for each of (A) the Van der Pol oscillator, (B) the Lorenz system, and the Lorenz 96 system with (C) 12, (D) 32, and (E) 64 dimensions.
26

---

Table S1: Trajectory prediction metrics (MSE and Rta) for each system and training noise. Errors
mean + standard deviation across five random initializations of the model architectures.
MSE
R²
Project
Training noise
JacobianODE
NeuralODE
JacobianODE
NeuralODE
VanDerPol (2 dim)
1%
0.00017 = 0.00007
0.0004 + 0.0002
0.99987 + 0.00005
0.9996 # 0.0002
= 0.0006
0.0009 =
0.9995 = 0,99949 0.00009
10%
0.004 0.001
0.0041 0.0009
0.9979 0.0005
0.9969 0.0003
Lorenz (3 dim)
1%
0.07 = 0.02
0,045 = 0,004
0.99979 = 0.00004
0,99986 + 0.00001
5%
0,12 = 0.03
0.10 = 0.01
0,9996 + 0.0001
0,99958 + 0.00004
10%
0.16 + 0.02
0.41 + 0.05
0,99928 # 0.00008
0.9981 1: 0.0003
Lorenz (12 dim)
1%
0.003 0.0007
0.033 0.003
0,99969 0.00007
0.9969 # 0.0003
5%
0.014 = 0.001
0.051 0.003
0.9985 0.0001
0.9955 + 0.0003
10%
0.038 = 0.002
0.063 + 0.005
0.9960 # 0.0002
0.9940 + 0.0007
Lorenz (32 dim)
1%
0.079 = 0.004
0.23 + 0.01
0.9939 + 0.0003
0.977 1 0.002
0.079 = 0.004
0.20
0.9939 + 0.0003
0.974 0.002
10%
0.205 0.001
0.34 0.02
0.9837 0.0001
0.968 0.002
Lorenz 96 (64 dim)
1%
0.72 = 0.03
1.194 + 0.005
0.946 + 0.002
0.8997 + 0.0007
5%
0.79 = 0.01
1, 1,205 + 0.007
0.9400 = 0.0010
0.8986 + 0.0010
10%
0.87 = 0.01
1.23 + 0.01
0.9332 + 0.0009
0.896 +: 0.001
Task-trained RNN
1%
0.0042 0.0006
0.00558 0.00002
0.9981 0.0003
0.99762 0.00001
5%
0.0043 + 0,0006
0.00596 + 0.0000S
0.9981 + 0.0002
0.99747 : 0.00002
10%
0.0089 = 0.0001
0.00655 = 0.00003
0.99591 0.00004
0,99721 0.00001
Table S2:
Jacobian estimation metrics (MSE and 12²) for each system and training noise. Errors are mean standard deviation across five random initializations of the model architectures.
MSE
R
Project
Training noise
JacobianODE
NeuralODE
JacobianODE
NeuralODE
VanDerPol (2 dim)
1%
0,18 = 0.06
0.4 + 0.1
0.985 = 0.005
0.97 0.01
5%
0.16 0.02
0.17 +: 0.03
0.987 = 0.002
0.985 0.002
10%
0.50 0.04
1.3
0.958 0.003
0.86 + 0.03
Lorenz (3 dim)
1%
2.4 0.1
1.70
0.954 = 0.002
0.77 0.01
5%
4.6 1.2
82.0
0.91 0.02
-0.6 0.2
10%
5.7 0.2
83.0
0.891 0.004
-0.58 0.10
Lorenz 96 (12 dim)
1%
0.012 0.003
0.18 0.01
0.9979 0.0006
0,969 0.002
5%
0.053 0.006
0.27 0.02
# 0.001
0,953 0.003
0.153 0.008
0.27 0.01
+ 0.001
0.951 0.002
Lorenz (32 dim)
1%
0.075 0.004
0.28 0.01
0.965 + 0.002
0.866
0.060 0.003
0.02
0.972 + 0.001
0.850
10%
0.180 0.003
0.38 + 0.02
0.915 = 0.001
0.820 0.008
Lorenz 96 (64 dim)
1%
0.238 0.007
0,515 + 0.002
0.773 = 0.007
0.509 0.002
5%
0,286 0.003
0.519 + 0.003
0.727 = 0.003
0.505 0.003
10%
0.321 0.004
0.527 & 0.005
0.694 + 0.003
0.497 0.005
Task-trained RNN
1%
2.2
5.2777 & 0.0009
0.59 + 0.03
-0.0008 + 0.0002
5%
1.70 + 0.08
5,2900 + 0.0001
0.68 = 0.01
-0.00315 = 0.00002
10%
1.99 + 0.01
5,2889 + 0.0001
0.623 = 0.002
-0.00295 = 0.00002
f(x(t)) (( x(1) x(1o) t-to
(15)
where
f. J(c..(t))d
(16)
Alternatively, one could use Equation 14d learn f(x((+)) and x(.)) (4:) as learnable parameters (x(4:) needed the first point of the path, i.c. c(t,) x(t.) inside the integral) [92]. Here, the path integral
27

---

Table S3: Mean 2-norm error on Jacobian estimation, [I] - - 12.)). for each system and noise level. Errors reported + standard deviation, with mean and standard deviation computed over five random initializations of the model architectures.
Project
Training noise
JacobianODE
NeuralODE
Weighted Linear
VanDerPol (2 dim)
1%
0.7 0.1
1.0 # 0.3 0.71 + 0.08
6.05
5 5%
0.71 = 0.05
6.03
10%
1.31 = 0.05
2.2 + 0.2
6.02
Lorenz (3 dim)
1
3.2 + 0.2
8.6
17.06
5 % 10%
4,9 + 0.9
25,9 1.5
17.02
6.0 + 0.1
26,4
16.95
Lorenz 96 (12 dim)
1
0.8 + 0.1
0.2
14.94
5 5% 10%
1.75 + 0.09
4.1
14.93
2.91 + 0.09
4.2
14.92
Lorenz 96 (32 dim)
%
3.75 + 0.08
7,8 0.1
16.29
5%
3.10 + 0.09
+11
16.26
10%
4.89 + 0.05
8.6 0.1
16.28
Lorenz 96 (64 dim)
1
8,4 + 0.1
12.63 # 0.02
16.84
9.04 + 0.07
12.66 0.03
16.81
10%
9.42 + 0.06
12.72 # 0.05
16.82
Task bask-trained RNN
1%
30.4 + 0.4
38.565 = 0.006
39.29
5%
29.4 + 0.9
38.5611 0.0004
38.91
10%
36.0 + 0.1
38.5600 = 0.0004
38.70
between x(4) and x(tr) is a linear interpolation, before. These ablated models were trained on 10 time-step prediction, as with the original JacobianODEs. For these models, we completed full sweep over the loop closure loss weight Atoop in order to determine the best hyperparameter.
Ablating teacher forcing We trained models without any able to generate only one-step predictions, without any recursive predictions. Again we did a full hyperparameter sweep pick Atoop.
That is, models were
Ablating loop closure loss We ablated the loop closure loss in two ways. The first was to set Auoop 0 to illustrate what would happen if there were constraints placed on the learned Jacobians. The second was to instead use the Jacobian Frobenius norm regularization that was used for the NeuralODEs (details are in Appendix D.3). We did a full sweep pick Ajaes the Frobenius norm regularization weight.
Table S4: Mean Frobenius norm enror (||J - j(e)) for different model ablations on the Lorenz system with 10% observation noise. Errors are reported as mean + standard deviation, with statistics computed over 8 test trajectories, each consisting of 1200 points.
Model variant
Frobenius norm error
JacobianODE (original)
6.46 # 1.50
With learned base derivative point
7.87 # 1.14
No teacher forcing
11.59 + 1.30
No loop closure
58.22 # 2.00
With Jacobian penalty instead of loop closure
9.59 + 2.45
Ablation results The performance of the ablated models on Jacobian estimation in the Lorenz system presented in Table S4 The original JacobianODE outperforms all ablated models, indicating that all components of the JacobianODE training framework improve the model's 's performance in this setting. Ablating the aacobian-parametter decrease in the estimation loss. This is potentially because the network could offload some of the
initial derivative estimate resulted in slight
28

---

responsibility for generating correct trajectory predictions onto the estimated base point x xf+a) and derivative estimate f(x(t))), slightly reducing the necessity of estimating correct Jacobians. Ablating the teacher forcing annealing predictably led to worse Jacobian estimation, as the network no longer has to consider how errors will propagate along the trajectory. The most dramatic increase in error was with the ablation of the loop closure loss. Without this important regularization, the learned Jacobians reproduced the dynamics but were not constrained be conservative, resulting in poor Jacobian estimation. The inclusion of the Frobenius penalty on the Jacobians mitigated this, although it did not encourage accurate Jacobian estimation to the same degree as the loop closure loss.
Table S5: Mean Frobenius norm error [IJ J (|e) for different model ablations on the task-trained RNN with 10% observation noise. Errors reported as mean + standard deviation, with statistics computed over 409 test trajectories, each consisting of the 49 points from the second delay and response epochs.
Model variant
Frobenius norm error
JacobianODE (original)
180.13 # 1.55
With learned base derivative point
186.28 # 1.17
No teacher forcing
187.72 # 1.63
No loop closure
313.31 + 29.88
With Jacobian penalty instead of loop closure
163.69 + 4.22
We then tested the ablated the models on Jacobian estimation in the task-trained RNN, with results presented in Table [S5] Again, ablating the Jacobian-pareeee paramltereaa derivative estimates, teacher forcing, and loop closure resulted in worse Jacobian estimation. Interestingly, in this setting, the inclusion of a penalty on the Frobenius norm of the Jacobians outperformed the use of the loop closure loss. This could potentially be because the loop closure loss is more difficult to drive to zero in high dimensional systems, or because the loop closure loss is more important in chaotic systems like the Lorenz system considered above. Future work should consider in what contexts each kind of regularization is most beneficial to JacobianODE models.
C.5 NeuralODEs achieve improved performance the cost of increased inference time
In the main paper, implemented both the JacobianODEs and the NeuralODEs as four-layer MIPs, with the four layers having sizes of 256, 1024, 2048, and 2048 respectively. This was done for the fairest architectural comparison between the models, to ensure that both models had the same representational capacity when generating their respective outputs. However, there are many architectural changes that we could make to this setup that impact performance. We hypothesized based on the discussion in Appendix C.2 that increasing the hidden layer size of the NeuralODEs would improve Jacobian estimation, as larger models have been known to learn smoother representa- tions. Furthermore, we wondered whether including residual blocks in place of the standard MLP implementation would improve Jacobian estimation.
To test this, we implemented the NeuralODEs as four-layer residual networks and tested three different sizes of hidden layer: 1024, 2048, and 4096. Results are in Table [S6] For nearly all models, these changes yielded only marginal improvements over the original NeuralODE model. Only the model with 409e-dimelll hidden layers under 10% training noise achieves a performance near the original JacobianODEs. However, the performance limit is still below that of the JacobianODEs, even with large increase in the representational capacity of the model. It is furthermore of note that the NeuralODEs were able to significantly improve performance only in the high noise setting. This suggests that high noise is necessary for the model be forced learn the response of the system to perturbation. In contrast, JacobianODEs perform similarly across all noise levels, indicating a stronger inductive bias to learn the response of the system to perturbation.
While the increased hidden layer size seems to induce smoother hidden representations (as expected), it comes with increasing computational cost. Recall that the NeuralODEmodels directly parameterize the dynamics as x (x). The Jacobian at state is computed by directly backpropagating through the Neural ODE Ñ thereby significantly increasing compute. On the other hand, the JacobianODEs only require forward pass. We evaluated the Jacobian inference time of JacobianODE
29

---

Table S6: Mean Frobenius norm error [|| - F] for different model types on the task-trained RNN data across observation noise levels. Errors are reported mean + standard deviation, with statisties computed over 409 test trajectories.
Training noise 5%
Training noise 10%
Learnable purameters
Model type
JacobianODE (original)
174.1 = 2,4
180.1 = 1.5
4.02c+07
NeuralODE (original)
294.0 =
293.9 1
6.85c406
NeuralODE (1024 dim. hidden layens)
291.6
289.8
3.41c406
NeuralODE (2048 dim. hidden layers)
288.5 2.3
275.0 2.7
1.31c407
NeuralODE (4096 dim. hidden layers)
275.8 2.4
184.2 4 5.6
5.14c+07
and NeuralODE models with four hidden layers of the same size H100 GPU. Each model was timed on inferring the Jacobians of 100 batches of the 128-dimensional task-trained RNN, with each batch containing 16 sequences of length 25. Timings repeated ten times for each model. Results are plotted in Figure S3
12.5
10
100 batch wall-clock time (s)
10.0
7.5
10
5.0
2.5
10 -1
0.0
2000
4000
10
10
Hidden dimension size
Hidden dimension size NeuralODE
JacobianODE
Figure S3: JacobianODEs achieve highly efficient Jacobian inference. Jacobian inference times computed over repetitions of 100 batches (error bars indicate mean + standard deviation). Red circles indicate the inference time corresponding to the largest hidden layer dimension of the highest performing NeuralODE model in Table|S6] Green eircles indicate the inference time corresponding to the largest hidden layer dimension of the highest performing JacobianODE model in Table S6 Each plot illustrates the same data but with different and y scaling.
The JacobianODE achieves much faster inference times than the NeuralODE approximately two orders of magnitude faster at large hidden dimension sizes. Furthermore, as shown in Table S6 the JacobianODE with a maximum hidden layer size of 2048 outperforms the NeuralODE with a maximum hidden layer size of 4096 Jacobian estimation, and does so with orders of magnitude faster inference (Figure S3] green and red circles, Table S6). This suggests that while both architectures appear to have inference times that scale approximately exponentially, the JacobianODE achieves more favorable scaling across every hidden layer size we tested. Our analysis therefore illustrates that it is possible to improve the NeuralODEs' Jacobian estimation with larger models, but the inference time scaling renders these models ill-equipped for important settings such as real-time control and the analysis of hogu-volume neural data.
C.6 Comparison echo state networks
Echo state networks (ESNs, or reservoir networks) are tools for time-series forecasting and analysis. ESNs are typically implemented as RNNs with fixed, sparse, hidden connectivity. The output weights
30

---

of the ESNs can then be trained via linear regression to reproduce specific temporal patterns. ESNs have been shown to be able to forecast the Lorenz 63 system for a particular choice hyperparameters 139]. We note that Pathak et al. [39 implements the same method as Pathak et al. [140] We implemented the ESNs exact hyperparameters reported in Pathak et. al. 2017 for the Lorenz system, which involved a reservoir of 300 nodes. Our implementation of the Lorenz system was simulated with sampling interval of 0.015 s. very similar to the 0.02 s used in the paper. Thus, for consistency of comparison with Jacobian we fit ESN models to data sampled every 0.015 s.
To enable Jacobian estimation, we changed one small implementation detail from Pathak et al. 139 Specifically, we trained the network output weights to predict the third variable using only the network state r(4), as opposed to the vector F(()) for which the node state is squared (t) = (1)
for half of the nodes i. This is because squaring the reservoir state complicates Jacobian computation, and was stated by Pathak et al. [139] and Lu et al. [111] to only be necessary for symmetry breaking if only the variable of the Lorenz system were provided. We trained the ESNs on the same volume of data as the JacobianODEs and NeuralODEs, with 1% observation noise. We swept the regularization parameter 3 to ensure good fit. We furthermore generated predictions using 150 time-step lead-up to ensure the ESNs converged to good internal representation. The ESNs predicted the Lorenz system well, achieving an MSE of 0.03 on the 10-step prediction task from our paper (on test data). This was very similar to the reported 10-step prediction MSEs of 0.07 and 0.04 for the JacobianODEs and NeuralODEs, respectively.
In Pathak et al. [39], the Lyapunoy exponents of the system were computed using the Jacobian of the interior dynamics of the ESN, as the recursion necessary to compute the Jacobian can be formulated uniquely only in regard to this state representation. Thus, each of these Jacobians would be a 300 x 300 matrix. While this suffices to compute the Lyapunov exponents, it does not suffice for our purposes of obtaining Jacobian estimates of the dynamics of the original system. For this purpose, it has been shown that it is possible obtain an estimate of the Jacobian of the original dynamical system (i.e., not the echo state space) using ESNs [142]. However, as Banerjee et al. [142] noted, is necessary to invert the non-square ESN output matrix in order to obtain these Jacobians, thus they are non-uniquely determined, and are used in the paper "only for causality estimation purposes". Nevertheless, we computed ESS-esttmated Jacobians on the 1% observation noise test data. We used a large (7,200 time-step) lead-up and ensured coherence to the trajectory by feeding in the true state as input. The Frobenius norm error over five random seeds was 85.8 + 3.1, much larger than the reported errors of 3.3 0.2 and 8.7 = 0.3 for the JacobianODEs and NeuralODEs respectively. Thus, similar to the NeuralODE, while ESNs can achieve comparable prediction on chaotic systems, they do so without enabling an accurate reconstruction of the Jacobian of the original dynamical system.
C.7 Analyzing the dimensionality of Jacobian manifolds
RNN activity is often confined to low-dimensional manifolds, potentially simplifying Jacobian estimation. To address the extent to which this influences our results, we performed the following analyses.
Jacobian rank: We first note that both the RNN hidden weight matrix and nearly all Jacobians had full rank (128; ~0.1% 1% had rank 127). To assess the effective rank, we computed the participation ratio (PR), singular malue-based measure of intrinsic dimensionality. The RNN weight matrix had a PR of 123, and the Jacobians had a mean PR of - 106 (std 0.7).
Intrinsic- -ata dimensionality: Another possibility is that the data itself may be low-dimensional. Since the Jacobian is a function of the trajectory data, the Jacobian manifold's intrinsic dimension is at most that of the trajectory manifold. To estimate the intrinsic dimensions, from a large sample of trajectories. For RNN state vectors, we computed PR directly; for Jacobians, we flattened them into n² -ddimensional vectors before computing PR. We found the RNN trajectory PR to be ~13.5, and the Jacobian PR to be -6. Given the 128D extrinsic space, this confirms both RNN states and Jacobians lie domer-dimensional manifolds.
pooled time points
Notably, prior work estimated the dimensionality of prefrontal cortex delay activity at using ~4000 neurons, falling to -6 when subsampled to 100 neurons [43]. These values were interpreted as high-dimensional. Thus, the dimensionalities we observe are consistent with biological systems performing complex tasks.
31

---

Furthermore, in the 64 04-dimensional Lorenz 96 system, the trajectory PR was -47.5, and the Jacobians had PR of -50, which we believe be nigh-dimensional enough demonstrate the effectiveness of our method on complex systems.
We note that in larger systems, approximating the Jacobian as low rank may be more efficient. Rather than computing each of the outputs directly, we can constrain the Jacobian to be low-rank by designing the oulput as a low-rank combination of vectors (e.g., it learns to output a low-rank singular value decomposition).
D Experimental details
D.1 Dynamical systems data
We used the following dynamical systems for the testing and validation of Jacobian estimation.
Van der Pol oscillator We implement the classic Van der Pol oscillator [96]. by the following equations
The system is governed
á =
111 Jy
We pick = 2 in our implementation.
Lorenz We implement the Lorenz system introduced by Lorenz 997
a a(y x)
ú = x(p 2) y
1 .y 3z
with the typical choices for parameters (a 10, p= 28, 3 8/3).
Lorenz 96 We implement the Lorenz 96 system introduced by Lorenz 98 and defined by
d; (I(++1) (mod N) T(i-2) (mod ()))(((1) (mod N) I2 (mod N) + F
with F' = and N (12, 64).
Simulation We use the dysts package to simulate all dynamical systems 99 100].
The charac-
teristic timescale of their Fourier spectrum 7 is selected and the systems are sampled with respect to т. For all systems, the training data consisted of 26 trajectories of 12 periods, sampled at 100 time steps per T. The validation data consisted of 6 trajectories of 12 periods sampled at 100 time steps per r. The test data consisted of trajectories of 12 periods sampled at 100 time steps per т. Trajectories were initialized using random normal distribution with standard deviation 0.2. The simulation algorithm used was Radau. Batches were constructed by moving a sliding window along the signal. The sequence length was selected such that the generated predictions would generate 10 novel time points (i.e., 11 time steps for the NeuralODE, and 25 time steps for the JacobianODE, due the 15 time steps used estimate the initial time derivative f).
Noise We define P% observation noise with P = 100p in the following way. Let Asignal
[1 be the expected squared norm of the signal with x(t) E R", Then consider a noise signal (t) E I2" where each component (t) ~ NoO, P Assgnal. Then
1 -p² Asigmal p² Asignal (22(1) = Σ
E
MON = E
(1)
=
E
1
1=1
1=1
and thus the noise percent is
E l
pa Asignal p
E |
Asigual
32

---

D.2 Task-trained RNN
The task used to train the RNN was exactly as defined in Section [ The hidden dimensionality of the RNN was 128, and the imput dimensionality was 10, where the first four dimensions represented the one-hot encoded "upper" color, the second four dimensions represented the one-hot encoded "lower" color, and the last dimensions represented the one-hot encoded cue. The RNN used for the task had hidden dynamics defined by
rh(1) - W (aar(h()) + Waau(() + ((l) = Woxh(()
with 50 ms, which for the purposes of training was discretized with Euler integration with time step of At 20 ms. Wih is the 128x dimensional matrix that defines the internal dynamics, Wh: the 128 10 dimensional matrix that maps the input into the hidden state, W oh is the 4x 48 dimensional output matrix that maps the hidden state to a four-dimensional output o(t), and b is static bias term. a was taken be an exponential linear unit activation with a 1. The RNN hidden state was split into two "areas" each with 64 dimensions. The input matrix Wh hi masked during training so that inputs could only flow into the first 64 dimensions the "visual" area. The same procedure was performed for the output matrix Wo.. except the mask was such that outputs could stem only from the second group of dimensions the "cognitive" area. The within-area subblocks of the matrix Wh were first initialized such that the real part of the eigenvalues were randomly distributed on the interval [-0.1, 0J and the imaginary part of the eigenvalues were randomly distributed on the interval [0, 2n... The eigenvectors were random orthonormal matrices. We then computed the matrix exponential be random normal, then divided by the 2-norm of the resulting matrix and multiplied by 0.05. The input and output matrices were initialized as random normal and then scaled by the 2-norm of the resulting matrix. The static bias b was initialized O. After initialization, all weights could be altered unimpeded (except for the masks). Notably, the inputs u(l) were present only during the stimulus and cue presentation epochs otherwise the network evolved autonomously. The loss was computed via cross entropy loss on the RNN outputs during the response period (the final 250 ms of the trial).
this matrix. The across-area weights were first initialized to
For the training data, we generated 4096 random trials, and used 80% for training and the remainder for validation. The batch size used was 32. Training performed for 40 epochs. The learning rate was 0.0005. For use with the Jacobian estimation models, data was batched and used for training exactly as done with the other dynamical systems data (see Appendix D.I). Observation noise was also computed in the same way.
D.3 NeuralODE details
NeuralODE models directly estimate the time derivative with a neural-network parameterized function &.. Then the Jacobians be computed as J f"
The NeuralODEs were implemented as described in Section and ODE integration was done exactly as for the JacobianODE using the torchdiffee ffeq package with the RK4 method [102] To regularize the NeuralODE we implemented Frobenius norm penalty on the estimated Jacobians, i.e.
Line
Aino ((x(t))
where j the estimated Jacobian computed via automatic differentation and Ajac a hyperparameter that controls the relative weighting of the Jacobian penalty [06] As mentioned in the main text, this penalty prevents the model from learning unnecessarily large eigenvalues and encourages better Jacobian estimation.
D.4 Weighted linear Jacobian details
We implemented baseline Jacobian estimation method using weighted linear regression models as described in Deyle et al. [88]. Given reference point x(t") which the (discrete) Jacobian will be
33

---

computed, all other points are weighted according to
(1x(((x) x(()
W/ = exp
d
where
d = Σ x((t)
the average distance from x(t)) to all other points. We then perform a linear regression using the weighted points (and bias term), the result of which is estimate of the discrete Jacobian at x(t"), which can be converted to continuous time by subtracting the identity matrix and dividing by the sampling time step (i.e., J = Jasus= where J disccrte is the discrete Jacobian). The parameter 0 tunes how strongly the regression weighted towards local points. To pick 6, we sweep over values range from 0 to 10, and pick the value that yields the best one-step prediction according
x(t" + 2 = x(t" + +t) + (x + At) - x(t"))
where J is the estimated Jacobian. This form of prediction has been previously been used to learn Jacobians in machine learning settings (90]. To test the method, we pick 0 based on data with observation noise particular noise level, then add in the denoised data to the data pool in order to
compute regressions and estimate Jacobians at the true points.
D.5
Model details
All models were implemented as four-layer MLPS, with the four layers having sizes of 256, 1024, 2048, and 2048 respectively. All models used a sigmoid linear unit activation. JacobianODE models output to the dimension n² which was then reshaped into the matrix of the appropriate dimension. NeuralODEs output to the dimension n.
D.6 Lyapunov spectrum computation
We discretize the Jacobians using the matrix exponential (i.e., J Jiscrette = eفب and then propagate bundle of small vectors through the Jacobians, using QR to ensure the perturbations remain bounded.
To compute the Lyapunov spectrum, we employ a QR based algorithm 121 122
D.7 Iterative Linear Quadratic Regulator (ILQR)
We implement the standard algorithm for ILQR, the details of which can be found in Li and Todoroy (II and Tassa al. [12]. In brief, the ILQR algorithm linearizes the system dynamics around a nominal trajectory using the Jacobian, and then iteratively optimizes the control sequence using forward and backward passes to minimize the total control cost. The state cost matrix Q was diagonal matrix with 1.0 along the diagonal. The final state cost matrix Q was a diagonal matrix with 1.0 along the diagonal. The control cost matrix R was a diagonal matrix with 0.01 along the diagonal. The control matrix was a 128 x 128 matrix in which the 64 64 block corresponding to the first 64 neurons (the "visual" area) was the 64-dimensional identity matrix. The control algorithm was seeded with only the initial state of the test trajectory with 5% noise. The control sequence was initialized random normal with standard deviation 0.001 and mean The ILQR algorithm was run for max of 100 iterations. The regularization was initialized at 1.0, with minimum of 1 x 10-6 and a maximum of 1 x 1000, Ao was set to 2, as Tassa et al. [T22]. If the backward pass failed 20 times in row, the optimization was stopped. The list of values for the line search parameter was 1.1 - for & E 0,..., 9 (see Tassa al. (122 The linear model used for the linear baseline was computed via linear regression.
D.8 Training details
All models were implemented in PyTorch, The batch size used was 16. Gradients were accumulated for batches. Training epochs were limited 500 shuffled batches. Validation epochs were limited to 100 randomly chosen batches. Testing used all testing data. Training was run for a maximum of 1000 epochs, hours, or until the early stopping was activated (see Appendix (.8.6), whichever came first.
34

---

D.8.1 Generalized Teacher Forcing
The Jacobian can be best learned when training predictions are generated recursively (i.e., replacing x(t) by x(t)). However, chaotic systems, and/or systems with measurement noise (as considered here), this could lead to catastrophic divergence of the predicted trajectory from the true trajectory during training. We therefore employ Generalized Teacher Forcing when training all models [94]. Generalized Teacher Forcing prevents catastrophic divergence by forcing the generated predictions along the line from the prediction to the true state. Specifically, for given predicted state x(1) and true state x()), the teacher forced state is
x(t) (1 a)x(t) + ax((t)
with a [0, 1]. This effectively forces the predictions along line from the prediction the true state, by an amount with proportion a. a= corresponds to fully-forced prediction (i.e., one-step prediction) and a 0 corresponds completely unforced prediction (i.e., autonomous prediction). Hess et al. [94] suggested that good estimate of a is
max
max P
1
(17)
|
(P) r:2
where 1 the Jacobians of the modeled dynamics computed at data-constrained states, p indicates the batch or sequence index, and
T-2 [[ Jr-
(
effectively computes the discrete maximum Lyapunov exponent. our implementation, we compute [[ Il using QR-dccomposition based Lyapunov computation algorithm [121]. As the Jacobian of the dynamics is necessary to compute this quantity, the JacobianODEs enjoy an advantage over other models in that the Jacobians directly output by the model, and do not have to be computed via differentiating the model itself.
We furthemore employ slightly modified version of the suggested annealing process in Hess et al. 94 which sets ag and updates Ox as
On = 70n-1 + (1 7)a
where a is computed according to equation 17 Following the suggested hyperparameters, we set y 0.999 and update an every 5 batches. Once the teacher forced state x(t) computed, it can simply replace x(1) in equation to generate predictions.
D.8.2 Loop closure loss
We implemented a loop closure loss discussed in Section 3.3 and Appendix A.3 For each loop, we used 20 randomly chosen points from the batch. For each batch, we constructed the same number of loops there were batches. Path integrals were discretized in 20 steps and computed using the trapezoid method from torchquad [TOI].
D.8.3 Validation loss
All models were validated on 10 time-step prediction task with teacher forcing parameter 0 (i.e., autonomous prediction).
D.8.4 Learning rate scheduling
For all models, the learning rate was annealed in accordance with teacher forcing annealing. Given an initial and final learning rates 77 and fy we compute the effective learning rate as
? = ? + o (a,,) m7)
where ar the current value of the teacher forcing parameter and
(aa))=
On
(1 0₂ ) kan
35

---

(a,,) is scaling function with (11) = 1 and (0)) 0 and for which the shape of the scaling is controlled by the parameter k. For positive values of k. the scaling super-linear, and for negative values of l it is sub-linear. We use k = 1, ensuring that the learning rate does not decrease too quickly at the start of learning. We set 'i = 10-4 and fy = 10-6 for all models.
D.8.5 Optimizer and weight decay
All models were trained with PyTorch's Adam W optimizer with the learning rate as descried above, and weight decay parameter 10-1. All other parameters were default (31 = 0.9, 82 0.999, c = 10 $0_$).
We also used gradient norm clipping, with a clipping value of 1.0.
D.8.6 Early stopping
For all models, we implemented an early stopping scheme that halted the training if the validation loss improved by less than 1 1% for two epochs in row.
D.8.7 Added noise during learning
For the models trained on the task-trained RNN dynamics, we added 5% Gaussian i.i.d. noise (defined relative to the norm of the training data with observation noise already added). Noise was sampled for each batch and added prior to the trajectory generation step of the learning process. Additional noise was not added for the loop closure computation.
D.8.8 Hyperparameter selection
For the JacobianODFs, the primary hyperparameter to select is the loop closure loss Aloop- To select this hyperparameter, we trained JacobianODE models with Atoop e [0, 10-1 10-5 10-4, 10-3 10- 10-1, 1, 10). For each run, the epoch with the lowest trajectory validation loss (Lmmj) is kept. Then, for this model, we compute the one-step prediction error on validation data, the validation loop closure loss (Lirop) and the percentage of Jacobian eigenvalues on all validation data that have a decay rate faster than the sampling rate At' We exclude any models that meet any of the following criteria:
1. One-step prediction error greater than the persistence baseline. The persistence baseline is computed as the mean error between each time step &åt and the subsequent time step (k: 1) across the dataset, and constitutes sanity check for whether model is capturing meaningful information about the dynamics.
2. Loop closure loss greater than /nt, where n the system dimension (see Appendix D.11 for the derivation of this bound). As discussed in the main text, we are interested in Jacobians that not only solve the trajectory prediction problem, but that also are constructed so that the rows of the matrix are approximately conservative vector fields.
3. More than 0.1 % of the Jacobian eigenvalues have a decay rate faster than the sampling rate Ar. Since large negative eigenvalues do not impact trajectory prediction, the models may erroneously learn Jacobians with large negative eigenvalues. If the decay rate of these eigenvalues faster than the sampling rate, we can infer that the eigenvalues are not aligned with the observed data.
If none of the models that meet criterion (2) meet criterion (1), we discount criterion (2), as this suggests that loop closure loss below /n bound is too strict to obtain good prediction on this system. Additionally, if none of the models that meet criterion (3) meet criterion (I), we discount criterion (1), as this suggests that noise is very high in the data, which leads to both high one-step prediction error, and large negative eigenvalues to compensate for the perturbations introduced by the noise. Of the remaining models, we select the one with the lowest trajectory validation loss Lunj.
For the NeuralODEs, we needed to select the hyperparameter Ajac. which regularized the mean frobe- nius norm of the Jacobians computed through automatic differentiation. Again, to select this hyperpa- rameter, we trained JacobianODE models with Aja E [0, 10-6, 10 10-1, 10-3 10-2, 10-1, 1, 10). We followed exactly the above procedure with the exception of criterion (2), which was deemed unnecessary, as computing the Jacobians implicitly via a gradient of the model (using automatic differentiation) ensures that the rows of the matrix are conservative.
36

---

All other hyperparameters (model size, leamnigg rate, length of initial trajectory, number of discretiza- tion steps, etc.) were fixed for all systems. Given the wide range of systems and behaviors and dimensionalities that the JacobianODfss are capable of capturing, this indicates that the method is robust given a reasonable choice of these hyperparameters.
D.8.9 Model hyperparameters and training details
Below are presented the details of all Jacobian estimation models considered in the main paper.
Table S7: Hyperparameters used and model details for each system and noise level. Training time is reported in seconds.
Syetem
Model
Leep wiight
Jacchese parilt
Traimn
spacs
Loamming
Weight droy
Loomohie pirmeiton
(7.dmm)
1020-00
9480
21
010001
1.00a-05
(1204)
5%
00001
(02 (2)
11
0/0001
1,00s
0.5001
6.508c (06
5%
NeauaDDE
249.43
0.0001
1,00c 06
0.0001
6.564; (D6
10%
678-194
10
000001
1:00-166
(12001
6:508c440
Nearak300
0400100
109.11
18
010011
0.2001
6.505=40
5%
Dao
0.10
1.147e008
18
00001
1.00s 06
0.0001
6:578c100
5%
NeaaaDEE
0
0,0000
(00165)
15
00001
1,00 06
0.0001
6.566c
10%
200
000001
1:00-1h
012001
5.7/8:400
Nearak3D
1.000w03
30
0.0001
6:507
5%
Dao
0000
716 ass
12
0.0001
1.00
0.0001
6:877x005
Nearaka0g
0
0.0030
10
GU0001
1:00-10
(00001
6.50/2400
10%
18
000001
1:00-166
6.57c440
Near
000000
0/0001
6:637--09
5%
0000
23
0.0001
1.60x-06
0.0001
5%
Neara0D0
0
0,0030
30
0U0001
1,00
000001
0.633:400
10%
0.001
14
000001
3,605:100
dr))
Nerr
30
0.0001
6.206e-49
5%
0,001
27
00001
1.00x:06
0.0001
L.9722007
5%
NeaaaODE
0
0,000
4879
30
010001
1,00 16
0.0001
6,706c (06
10%
1700lw00
010001
1/90744177
0.10
9:77wwot
16
010001
05
4.554r-05
5%
DaotiuuDD 0.001
2350ar08
20
0.0001
1,00s 06
0.0001
4.016c 107
5%
NeaaODDE
1,0c
7.370e108
27
0/0001
1,00 06
0.0001
6354; (06
10%
DasimonciD:: 0m0
2.172w008
27
000001
1:000-10
0.0001
8o4477
10%
Noui
15
3300e10
0.001
0.0001
6.542:00
D.9 Information about computing resources and efficiency
All models were able to be trained on a single H100 GPU, with 80 GB of memory.
Jacobian inference times Jacobian inference times for the JacobianODE and NeuralODE models are discussed in Appendix C.3 As discussed, models were implemented with four hidden layers of the same size, and tested on 100 batches of the 128- 228-dimensional tusketrained RNN data, with each batch consisting of 16 sequences of length 25. Timings were repeated ten times for cach model (see Figure S3 for details).
Training time Total training times for each of the chosen models are presented in Table S7 Furthermore, we include the training time (including backward pass) for 100 batches (with 16 sequences per batch), using 10 time-step prediction, in Table
D.10 Statistical details
All statistics were computed using scipy. For the comparison between Jacobian and NeuralODE trajectory and Jacobian predictions, as well as the comparison of Gramian traces and minimum eigenvalues, we used a two-sample t-test. For the comparison of II.QR control accuracies and errors, we used Wilcoxon signed-rank test.
37

---

Table S8: Trajectory training time (seconds) for each system and noise level. Training used 100 batches, with 16 sequences per batch, as well as 10 time-step prediction.
Model
Lorenz (3 dim)
VanDerPol (2 dian)
Lorenz 96 (12 dim)
Lorenz 96 (32 dim)
Lorenz 96 (64 dim)
Task-trined RNN
1% noise
JacobianODR:
15.772
13.737
606
16.469
16.214 18.192
23.855
NeuralODE
8.215
8,403
11.852
12.075
14.955
5% noise
JacobianODE NeuralODE
12,422 6.191
16.772
10.655
13.121
12,949 21.307
23.197
5.841
13.501
8.021
30.209
10% noise
JacobianODD NeuralODE
12.590
15.719
18.447
10.506 8.008
18,900
22.482
D.11
Derivation of loop closure loss bound
We consider the loop closure loss as defined in 3.3 error
We are interested in estimating a bound on the
IE
Jds
71
where the system dimension. While in theory this quantity should be in practice due to numerical estimation error, it will not be. First recall, that
px((ee (mad z))
Jds
J(c(r))
x(L.
(mod mj)
where the number of loop points and c is a line from x(1, (mod 1)) to x(1(+11) (mod L)). We assume that
| X (mod [((())) =
( (mod L))
(x (mod 1)))
E
J(c(r))
x(1,
(mod L))
Vi,j E [I,.., LJ. which is justified as the numerical error accrued will likely be similar along different lines for the same system. Thus
1
L
x(ll (1). (mod b)
2
Jds
[(c(r))c((rdrr
Caap
71
x(ti (mod 2s)
pxit(413 (mod m.)
27
J(c(r))c((rdd
12
xits (mod 14)
x((6)41) (med 1))
27
J(c(r))
17
x(t, (mod 2))
x(t))
E
J(c(r))d(rrdh
71
x(to)
L
x(t1)
E
((cr)))c(()dd
n
x(to)
where we have used the fact that the errors are assumed be equivalent for each the line segments comprising the overall loop path. Recall now that for the trapezoid integration rule, the error E in
38

---

integrating ."' (x))dr be computed E (b-a)² emor bounded as
["(u) for some E [a, b). Thus the squared 12M
a)
max aEla.b]
12M f" 11
In our case, when path integrating lines, we are effectively integrating from a = to b = 1. Furthermore, let g(r) J(c(r))c((r), the integrand of our line integrations. We assume that E ||max, [|m]x,""()]l
v where the first n comes from the fact that the norm involves summing over the n components of the vector g"(r) and the second vi involves an assumption that loop closure loss will be more difficult to compute accurately in higher dimensions, though this will be more pronounced as dimensionality initially starts increasing. Thus E [|lmax, g((r) for some IR+, Now, continuing on,
kn vn
x(()) J(c(r))c(rrdtt
Jds
E
17
x(ta)
13
E
max g(())
11
12M
2
L
max g (r)
1
2
Lkn vn
12
Lkyn
12
Finally, assuming that the number of discretization steps M was chosen to large enough such that M 12z Lk we finally obtain
Jds
S vn
7i
chip
39

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The abstract and introduction claim the the paper contributes robust data-driven Jacobian estimation, a framework to characterize control between interacting subsystems, inference of control dynamics in trained recurrent neural networks, and accurate control the trained network. All claims are substantiated in detail through a description of the presented method and analytical framework, well as comprehensive testing. See Sections 3 and
Justification:
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 in the paper and important assumptions and limitations. A No or
NA answer to this question will not 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 long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The primary limitation of the present approach is that it was not tested on partially observed system dynamics, limitation which discussed in the Discussion (S The method was tested on a wide range of datasets (Section 4 and computational efficiency discussed in Appendix D.9
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 on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low 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 of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be 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 an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
40

---

3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer:
[Yes]
Justification: complete and correct proof of the presented theoretical result, that the time derivative f can be represented in terms of its Jacobian, is provided in Appendix A.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented
by formal proof's provided in appendix or supplemental material. Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: The paper provides full information about the presented methods and experi- ments, with sufficient detail that all results could be reproduced.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, 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 to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.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 reproducibility may be tricky in some cases, in which case authors are welcome 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 possible for other researchers to have some path to reproducing verifying the results.
41

---

reproducible, they 5. Open access to data and code
Does the paper provide open to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: As part of the submission, we provide an anonymized link to download the analyzed data, well as all code implementing the models, baselines, control analyses, and data generation. the paper is accepted, the full code will be open sourced.
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, so "No" an acceptable answer. Papers cannot be rejected simply for not 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 NeurlPS code and data submission guidelines (https
/ //nips. httppsscccpbllcccc/cc somblic/guuides guides/
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments 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 much information possible in supplemental material (appended the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question:
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The paper specifies all training and test details, including data splits, hyperpa- rameter selection, and training details such optimizers.
Guidelines:
The answer NA means that the paper does not include experiments. The experimental setting should be presented in the core 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, 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: The paper includes error bars and statistical information. Additionally, the factors of variability that the errors are capturing are clearly stated.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
42

---

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 given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We provide information about compute resources needed to reproduce the experiments in Appendix D.9
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of individual experimental runs 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 failed experiments that didn't make into the paper).
9. Code of ethies
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: Ethical considerations have been taken into account and the research conforms fully to the NeurIPS code of ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is 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: [Yes]
Justification: Yes, the paper discusses societal impacts, including potential applications to modeling biological dysfunction in disease. There expected negative impacts.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
43

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition 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 a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
This paper contributes an analytical framework for characterizing nonlinear control in interacting subsystems, and poses no risk of misuse. We do not make use of any language models or image generators, and do not use any scraped datasets. We therefore do not describe safeguards.
Justification:
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model. for example by requiring that users adhere 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
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: Yes, we implement existing methods our paper and cite the original creators appropriately.
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 or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
44

---

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
batasets.paaersswt
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 is 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 well-documented code for generating the models and implementing the datasets used in this paper.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing 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: This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in 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 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
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 an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor 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.
45

---

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, 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, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does not involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does involve LLMs as any important, original, non standaadd components.
Please refer to our LLM policy http:///eerrrsssccooofeeereeeeee for what should or should not be described.
46


Reviews:
[{"review_id": "kb7wzDr3Ir", "reviewer": "Reviewer_ya7v", "summary": "Authors are modeling dynamical systems by learning from observational data of trajectories. Instead of learning the vector field of the dynamical system directly, they reformulate the problem to learn the system's Jacobian matrix.\n\nThe authors recognize an important fact: each row of the Jacobian corresponds to the gradient of a scalar function, and therefore represents a conservative vector field. This allows them to impose structure in the learning process — specifically, that line integrals over these fields are path-independent.\n\nAlthough the full vector field f(x) does not need to be conservative, its Jacobian rows must be, since they are gradients by definition.\n\nBy applying a generalized form of the fundamental theorem of calculus, they express the function f(x(t)) — which is not modeled directly — as a path integral over the Jacobian.\n\nThis approach enables them to reuse standard ODE solvers to generate trajectories of the full nonlinear system — not just a linearized approximation.\n\nIn the loss function, they combine a trajectory reconstruction loss (to ensure consistency with observed dynamics) with a loop closure loss that enforces path independence, reflecting the conservative structure of Jacobian rows.", "questions": "For the case of chaotic dynamical systems, I would like to see comparison how this approach compares to echo state networks:\n\nsee: Doan, Nguyen Anh Khoa, Wolfgang Polifke, and Luca Magri. \"Physics-informed echo state networks for chaotic systems forecasting.\" Computational Science–ICCS 2019: 19th International Conference, Faro, Portugal, June 12–14, 2019, Proceedings, Part IV 19. Springer International Publishing, 2019.\n\nPathak, Jaideep, et al. \"Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach.\" Physical review letters 120.2 (2018): 024102.\n\nPlease state limitations more explicitly.", "limitations": "Please state limitations more explicitly.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- elegant theoretical framework that is quite applicable for range of dynamical systems\n- casting the problem from learning vector field of dynamics to jacobian learning\n- efficient way of learning, without explicit need to model directly vector field\n\nNo major weaknesses.", "was_revised": true, "final_justification": "Author have performed additional experiments with echo networks, and as expected echo state networks worked quite good, but not good in estimating Jacobian. They have explained what are limitations to their paper is. I recommend acceptance of paper.", "rebuttal": "We thank the reviewer for their thoughtful feedback, and appreciate that the reviewer found our theoretical framework to be both elegant and broadly applicable. We address here the questions posed by the reviewer.\n### Comparison to echo state networks\n\n_**TL;DR: Echo State Networks (ESNs) achieve comparable short-term prediction on chaotic systems, but fail to recover accurate Jacobians. This stems from structural limitations in mapping ESN dynamics back to the source system.**_\n\nThe reviewer aptly suggested a comparison to echo state networks (ESNs, or reservoir networks) in the case of chaotic dynamical systems. ESNs have been shown to be able to forecast the Lorenz 63 system for a particular choice of hyperparameters (Pathak, J. *et. al.* Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data. _Chaos_, 2017). We note that Pathak et. al. (2017) implements the same method as Pathak et. al. (2018) referenced by the reviewer. We implemented the ESN architecture and hyperparameters from Pathak et al. (2017) for the Lorenz 63 system, which involved a reservoir of 300 nodes. Our implementation of the Lorenz system was simulated with a sampling interval of \\~0.015 s, very similar to the 0.02 s used in the paper. Thus, for consistency of comparison with JacobianODEs, we fit ESN models to data sampled every 0.015 s. \n\nWe trained the ESNs on the same volume of data (with 1\\% observation noise) as the other models, and swept the regularization parameter $\\beta$ to ensure a good fit. To ensure stable internal dynamics, we used a 150-step warm-up before generating predictions. The resulting ESNs predicted the Lorenz system well, achieving 10-step test MSE of 0.03—comparable to JacobianODEs (0.07) and NeuralODEs (0.04).\n\nHowever, obtaining the Jacobian of the original Lorenz system from the ESN presents challenges. While the Jacobian of the internal reservoir dynamics (a 300×300 matrix) can be uniquely computed and was used in prior work to estimate Lyapunov exponents (Pathak et. al. 2017), recovering the Jacobian of the _original_ system requires inverting the non-square readout matrix (Banerjee, A. *et. al.* Using machine learning to assess short term causal dependence and infer network links. _Chaos_, 2019). As shown in Banerjee et al. (2019), this inversion yields a non-unique solution, making such Jacobians ill-suited for system identification, and are used in the paper \"only for causality estimation purposes\".\n\nWe estimated these Jacobians nonetheless, using a long (7200-step) warm-up and true state input for stability. The resulting Jacobians had a Frobenius norm error (averaged over five random seeds) of $250 \\pm 0.3$, compared to $3.3 \\pm 0.2$ for JacobianODEs and $8.7 \\pm 0.3$ for NeuralODEs. Thus, while ESNs can match predictive accuracy, they are significantly less accurate in capturing the system’s Jacobian, due to architectural limitations.\n\nWe have added a discussion of this comparison into the appendix of the paper, and we thank the reviewer for suggesting this valuable benchmark.\n### Stating limitations more explicitly\n\nWe have extended the limitations section to address limitations surrounding sensitivity of reachability and scalability, in addition to partial observability. The revised section is as follows: \n\n> A future challenge for JacobianODEs is partially observed dynamics. Recent work has identified that it is possible to learn latent embeddings that approximately recover the true state from partial observation [108-112].  Jacobian-based dynamics learning has been performed in a latent space [61, 113], yet it is unclear whether this translates to accurate Jacobian estimation, given the challenges related to automatic differentiation and Jacobian estimation presented here. \n>\n>We also note that reachability estimates depend sensitively on several factors: the alignment of cross-subsystem interactions and within-subsystem dynamics, the eigenvectors and eigenvalues of within-subsystem Jacobians, and the way activity propagates within each subsystem (see appendix C.1). While our method reliably captures broad trends in reachability over time, fine-grained, time point-specific comparisons should be interpreted with caution.\n>\n>Finally, although JacobianODEs scale well to moderately high-dimensional systems, their performance in systems that are orders of magnitude larger than those considered here (e.g., recordings of thousands of neurons via calcium imaging) remains to be tested. In many practical settings, this may not be a limitation: neural representations during tasks often exhibit intrinsic dimensionalities comparable to those studied here, enabling JacobianODEs to operate in a reduced dimensionality. Notably, most models converged in under 45 minutes on a single GPU, suggesting favorable scaling. Future work should explore the viability of the method in higher dimensions, and assess whether dimensionality reduction strategies (such as latent state models or low-rank Jacobian approximations) can further improve scalability and accuracy."}, {"review_id": "ioRxZl4zVX", "reviewer": "Reviewer_DEPk", "summary": "The authors provide a novel method of local deep Jacobian estimation for non-linear dynamical systems, involving optimising a neural network jointly on trajectory prediction and adherence to conservativity. The authors used this method to model the interarea controllability of a two-area task-trained RNN, showing that their method outperforms other Jacobian estimation methods in predicting and exerting control on local dynamics.", "questions": "- Why was this method of loop curation chosen? Do you expect other design choices to significantly change performance, beyond a simple binary ablation?\n- Figure 6B suggests that the model is systematically delayed in predicting changes to area reachability in opposite directions. Is there an intuitive explanation for why this is? Did you consider any investigations into this offset?\n- Did the authors consider training on other laboratory tasks where control theory has previously been applied, e.g. motor control? Comparison of insights provided by JacobianODE compared to previous analyses would provide powerful justification for the method", "limitations": "There are no relevant negative societal impacts associated with this work.\n\nThe authors sufficiently discuss limitations to their work in the final section of the paper. In particular, they highlight the need for an extension to partially observable dynamical systems. This would be a necessary step towards applying this method to real neural dynamics", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths\n- The authors provide a well-motivated Jacobian estimation technique, providing extensive results on synthetic data. \n- The authors are careful to provide details throughout the paper about the technical and intuitive details of their method e.g. appendix C.1, C.2 provide welcome information and results to reinforce the method motivation, and ablation studies justify each component of the model\n\nWeaknesses\nThis work presents an interesting method and theoretically strong results. The primary weaknesses of this paper are around the communication and format of information:\n- Section 2.1 seems out of place. While the authors do emphasise the importance of modelling interareal/intersystem communication in the introduction, this does not play a role in the training of JacobianODE, or understanding of the method as a whole. Introducing these concepts in any technical detail may be better placed near section 5, where interarea systems’ data are first introduced.\n- Overall, there is a poor separation of information in the main text and in the appendix. Section 3.1 provides much information that is not strictly necessary for understanding the training procedure in Figure 3, or at least requires an understanding of the appendix material to take advantage of this level of explanation. This could be moved to the appendix to create a clearer narrative in the main text. This would mirror other methods that were key to appreciating the results of this paper, e.g. ILQR, being delegated to the appendix.\n- Despite being an innovation for this work, section D.8.2 seems like insufficient exposition on the formulation and importance of the loop closure loss. Table 3 shows that this is critical for performance, but no investigation or context is provided on how one could curate optimal loops for training.\n\nLess importantly:\n- Images of trajectories provide little information in their current form. Dimensionality reduction, e.g. by tSNE, could provided a better qualitative visualisation of the quality of predictions by each model.\n- A final minor weakness is that spacing is limited throughout the paper, table captions are difficult to discern, and information on figures are often too small", "was_revised": true, "final_justification": "I am continuing to recommend this paper for acceptance, as per my original review and subsequent reponse.", "rebuttal": "We thank the reviewer for their helpful comments and feedback, and are glad that they found the method to be well-motivated, interesting, and theoretically strong. We now address the questions and comments posed by the reviewer.\n\n### Section 2.1 is out of place\n\nWe have moved section 2.1 to the analysis of the task-trained RNN (section 5), which we agree is a more conceptually aligned location.\n\n### Main text/appendix split \\& loop construction explanation\n\nTo clarify the narrative flow, we moved “Constructing the time derivative” and technical loop closure details to the appendix. In their place, we provide a more intuitive discussion of loop closure in the main text (see next point).\n\n### Formulation of the loop closure loss and methods of constructing loops\n\n_**TL;DR: The loop closure loss encourages the Jacobian to learn about directions in the tangent space that are orthogonal to the flow of the dynamics. We construct loops from concatenated lines between random points on the data manifold to encourage uniform sampling of these directions.**_\n\nWe agree that the loop closure loss deserves more explanation. We have revised the main text as follows:\n\n> The Jacobian captures how perturbations to the system will propagate along _any_ direction in state space. Estimating it purely from dynamics constrains only the direction of the flow, leaving the full solution underdetermined. To address this, we again exploit the fact that each row of the Jacobian is a conservative vector field. Specifically, we note that for any piecewise smooth loop $\\mathcal{C}\\_{\\text{loop}}$, we have $\\left \\Vert \\oint\\_{\\mathcal{C}\\_{loop}} \\mathbf{J} ds \\right \\Vert_2 = 0$. Thus, by integrating along loops that contain directions _orthogonal_ to the system's dynamics (and penalizing the deviation from zero), we encourage the estimated Jacobians to capture information about other directions in state space (see appendix A.3 for full technical details). To ensure broad coverage of tangent space directions, we form loops from concatenations of line integrals between randomly selected data points. This strategy samples diverse directions from the tangent space while remaining easy to compute. This self-supervised loss term builds on the loss introduced by Iyer et. al. [63]. It constrains $\\hat{\\mathbf{J}}^{\\theta}$ to satisfy both the dynamics and conservativity. This improves Jacobian estimation accuracy significantly (see appendix C.4 for ablation studies).\n\nAs mentioned, we have also moved the technical details to the appendix.\n\nOther choices (e.g., forward and backward loops along observed trajectories or arbitrary circles) were considered, but our method offers a strong tradeoff between tractability and uniformity in sampling. We have added a discussion of these other potential choices for loops to the appendix as well. \n\nWe appreciate the reviewer’s suggestion, which helped clarify a core component of our method.\n### Images of trajectories from sample systems\n\nWe agree that the images of the trajectories in their current forms may not be optimally informative. As suggested, we will run tSNE, as well as UMAP, on the data and add the most visually salient / illustrative representations to Figures 4, 5, and S2. We have also added time-series plots of the top Principal Components in the appendix to highlight the temporal unfolding of each system.\n### Limited spacing, table captions, and small figure information\n\nRegarding spacing, as mentioned, we have moved many details to the appendix, which should create more space. If accepted, we will further space out the text and figures to fill out the extra allotted page. \n\nWe have updated the following table captions and figure information to make it easier to read:\n- More detail has been added to the caption for Table 1, Table S1, and Table S6\n- Table S2 and S3 have been made larger\n- The text in Figure 6 has been made larger\n### Delayed reachability curve in Figure 6B\n\n_**TL;DR: The reachability Gramian is highly sensitive, and takes into account the alignment between cross- and within-area dynamics, as well as how within-area dynamics propagate over time. Small estimation errors can thus cause temporal offsets. We conjecture that small estimation errors in the within-area eigenvalues are the source of the observed delay, and now include a discussion of this phenomenon in our limitations.**_\n\nAs the reviewer astutely notes, in Figure 6B, cognitive-to-visual reachability peaks slightly late, while visual-to-cognitive hits a minimum slightly early. This likely stems from the sensitivity of the reachability Gramian, which arises from a complex interplay of components. These include the alignment of the cross-area interaction terms with the eigenvectors of the within-area dynamics, the corresponding eigenvalues of each of the eigenvectors, and how activity propagates within each area over successive time points. This phenomenon can be observed in the results in appendix C.1. Mathematically, consider the definition of the reachability Gramian currently in appendix B.1, in which the $\\mathbf{\\Phi}(t_1, \\tau)$ term captures how within-subsystem dynamics propagate between times $\\tau$ and $t_1$, and the multiplicative term $\\mathbf{\\Phi}(t_1, \\tau)\\mathbf{B}(\\tau)$ captures how the directions of interaction between the subsystems aligns with the within-subsystem dynamics. The Jacobian estimation for the task-trained RNN, while robust considering the Jacobian consisted of a $128 \\times 128$ matrix at each time point, was not perfect (we achieved an $R^2$ of 0.68 with 5\\% observation noise, as reported). Thus, for such a sensitive measure, even small Jacobian estimation errors can lead to timing deviations. \n\nWe can therefore conjecture that the source of this delay might have to do with small estimation errors in the eigenvalues of the estimated within-area dynamics matrices over time. Broadly speaking, when considering control over short time windows, (moderately) unstable eigenvalues in the target subsystem yield higher reachability, as activity can more easily be driven outward along their corresponding eigenvectors. As mentioned, the exact reachability values will depend sensitively on many factors, but the eigenvalues should provide a clue. For instance, the average predicted maximum eigenvalue in the visual area peaks later than the ground truth, which may explain the cognitive-to-visual delay. In the cognitive area, a transient dip in the average predicted maximum eigenvalue could cause the early minimum in visual-to-cognitive reachability.\n\nCrucially, the delays in the reachability curves do not obscure the core result of the figure. On average, the cognitive-to-visual reachability tended to peak later in the delay than the visual-to-cognitive reachability. We furthermore note that this bias does not extend to the actual Jacobian estimation, for which the average residual is very close to zero, indicating that the JacobianODEs constitute an unbiased estimator.\n\nThis analysis highlights a limitation of our approach that we did not previously address, which relates to the interpretability of the time-resolved reachability curves. As demonstrated, reachability is a highly sensitive measure. While our method reliably captures broad trends in reachability over time, fine-grained, time point-specific comparisons should be interpreted with caution. We have added this to the discussion of limitations in the paper, and we thank the reviewer for encouraging this clarification.\n### Training on other laboratory tasks (e.g. motor control)\n\nWe agree that applying our method to laboratory tasks involving motor control, where control-theoretic frameworks are well established, would be highly compelling. This is a promising direction for future work, and we appreciate the suggestion."}, {"review_id": "7cKmEVPrs8", "reviewer": "Reviewer_aLtD", "summary": "Authors introduce a modeling approach for predicting next states in an observed dynamical system. They specifically utilize the linearized dynamics and learn the Jacobian to both achieve an interpretable modeling of the observed dynamics (though locally) that later enables control of the system and be able to predict next states for a given initial condition and inferred set of time-varying jacobians. They design a novel learning algorithm to train these Jacobians and demonstrate these two specific claims with simulations.", "questions": "1) It is still not clear to me how Jacobian at a given time point can be approximated without repetition of the same trials/time windows/states? One would ideally need to sample around a particular state to be correctly estimate the jacobian? There is some coarse-graining here that I am not seeing/not being discussed explicitly. Can you comment on this?\n\n2) The model of the RNN should be shown in the main text. In my opinion, this work would significantly benefit from using more complex RNN architectures. Specifically, the RNN architecture used here is known to have neural activities confined to low-dimensional hyperplanes, hence the jacobian estimation is actually quite simpler and maybe effectively low-dimensional. You can check this by computing the PC dimensionality in learned RNNs, which would be important to note here. This may be ok, since most behavioral tasks animals perform in the labs are low-dimensional anyways, but this point should be discussed and raised in the limitations. In my opinion, using LSTMs trained on slightly more complex tasks (and validated to have high PC dimensionality) may be a good control for this. \n\n3) The discussion on Gramian is too quick. Could you please expand on this? The reader cannot really appreciate what is exactly being computed, how it relates to the experiments etc. \n\n\nHere are some simple extra suggestions:\n\n- Line 150: waiting -> weighting?\n\n- Title was a bit misleading for me. This seems to be a neuroscience-relevant work, though only looking at the title, I would have likely not read this work. (and it would have been a mistake on my part, since this is quite interesting!)\n\n- The transition to Eq. (1) is weird. At first, I thought that was part of the figure.\n\n- Section 6 should likely come way earlier, before the results section.", "limitations": "See above.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths\n\n- The focus on predicting not just dynamics, but also the jacobian is novel and quite interesting. This is particularly relevant since the multiplication J(t) x is under constrained by just observed dynamics, but Jacobian is needed for defining optimal control dimensions.\n\n- The authors demonstrate both of their claims with simulations and RNN training, which is almost sufficient (see below).\n\n## Weaknesses\n\n- The writing is quite dense, which leads to many questions in the mind of (this) reader. In several cases, I felt that the results were too good to be true (which can also become a strength once authors clearly set the boundaries of their method).\n\n- I believe a different RNN architecture may be more appropriate to test these claims (see below).", "was_revised": true, "final_justification": "I think this work is a nice addition to Neurips, as agreed by all reviewers.", "rebuttal": "We thank the reviewer for their constructive feedback and questions about the paper, and address them here.\n\n### Dense writing and boundaries of the method\nWe rewrote section 3 to clarify key intuitions and moved technical details to the appendix for improved narrative flow. If accepted, we will use the extra page to space out text and figures.\n\nTo better define the method’s boundaries and clarify our claims, we substantially revised the limitations section (see also reviewer 1):\n\n> A future challenge for JacobianODEs is partially observed dynamics. Recent work has identified that it is possible to learn latent embeddings that approximately recover the true state from partial observation [108-112].  Jacobian-based dynamics learning has been performed in a latent space [61, 113], yet it is unclear whether this translates to accurate Jacobian estimation, given the challenges related to automatic differentiation and Jacobian estimation presented here. \n>\n>We also note that reachability estimates depend sensitively on several factors: the alignment of cross-subsystem interactions and within-subsystem dynamics, the eigenvectors and eigenvalues of within-subsystem Jacobians, and the way activity propagates within each subsystem (see appendix C.1). While our method reliably captures broad trends in reachability over time, fine-grained, time point-specific comparisons should be interpreted with caution.\n>\n>Finally, although JacobianODEs scale well to moderately high-dimensional systems, their performance in systems that are orders of magnitude larger than those considered here (e.g., recordings of thousands of neurons via calcium imaging) remains to be tested. In many practical settings, this may not be a limitation: neural representations during tasks often exhibit intrinsic dimensionalities comparable to those studied here, enabling JacobianODEs to operate in a reduced dimensionality. Notably, most models converged in under 45 minutes on a single GPU, suggesting favorable scaling. Future work should explore the viability of the method in higher dimensions, and assess whether dimensionality reduction strategies (such as latent state models or low-rank Jacobian approximations) can further improve scalability and accuracy.\n\n### The model of the RNN should be shown in the main text\nWe agree that this was an oversight. We now introduce the RNN architecture directly in section 5, as suggested.\n### Suggested analysis of more complex RNN architectures, due to potentially low-dimensional Jacobian estimation\n\n_**TL;DR: The Jacobians were not low-rank. While both RNN activity and Jacobians lie on low-dimensional manifolds relative to the 128D state space, they match dimensionality observed in neural data, as the reviewer noted. The 64D Lorenz 96 system had much higher intrinsic dimensionality, supporting our method’s effectiveness on complex systems.**_\n\nThe reviewer astutely notes that RNN activity is often confined to low-dimensional manifolds, potentially simplifying Jacobian estimation. We address this point through the following analyses.\n\n**Jacobian rank**: We first note that both the RNN hidden weight matrix and nearly all Jacobians had full rank (128; \\~0.1\\% had rank 127). To assess the _effective_ rank, we computed the participation ratio (PR), a singular value-based measure of intrinsic dimensionality. The RNN weight matrix had a PR of 123, and the Jacobians had a mean PR of ~106 (std 0.7).\n\n**Intrinsic-data dimensionality**: As the reviewer notes, the data itself may be low-dimensional. Since the Jacobian is a function of the trajectory data, the Jacobian manifold’s intrinsic dimension is at most that of the trajectory manifold. To estimate the intrinsic dimensions, we pooled time points from a large sample of trajectories. For RNN state vectors, we computed PR directly; for Jacobians, we flattened them into $n^2$-dimensional vectors before computing PR. We found the RNN trajectory PR to be ~13.5, and the Jacobian PR to be ~6. Given the 128D extrinsic space, this confirms both RNN states and Jacobians lie on lower-dimensional manifolds.\n\nNotably, prior work estimated the dimensionality of prefrontal cortex delay activity at 24 using ~4000 neurons, falling to ~6 when subsampled to 100 neurons (Rigotti _et al._, The importance of mixed selectivity in complex cognitive tasks. _Nature_, 2013). These values were interpreted as high-dimensional. Thus, the dimensionalities we observe are consistent with biological systems performing complex tasks.\n\nFurthermore, in the 64-dimensional Lorenz 96 system, the trajectory PR was \\~47.5, and the Jacobians had a PR of \\~50, which we believe to be high-dimensional enough to demonstrate the effectiveness of our method on complex systems. \n\nWe note that in larger systems, approximating the Jacobian as low rank may be more efficient. Rather than computing each of the $n^2$ outputs directly, we can constrain the Jacobian to be low-rank by designing the output as a low-rank combination of vectors (e.g., it learns to output a low-rank singular value decomposition).\n\nWe agree that applying our method to more complex architectures (e.g., LSTMs) and tasks is an important direction for future work. For now, we believe our results sufficiently demonstrate the method’s ability to capture rich, high-dimensional dynamics.\n\nAs suggested, we added a discussion of dimensionality to the limitations section, as well as to the appendix. We thank the reviewer for prompting this clarification of an important aspect of the method.\n\n### How the Jacobian at a given time point can be approximated without repetition\nWe appreciate this question and recognize it was not fully addressed. The loop closure loss provides the primary mechanism as it constrains the Jacobian’s form beyond fitting observed trajectories. We have rewritten the loop closure section as follows (and see reviewer 2 discussion):\n\n> The Jacobian captures how perturbations to the system will propagate along _any_ direction in state space. Estimating it purely from dynamics constrains only the direction of the flow, leaving the full solution underdetermined. To address this, we again exploit the fact that each row of the Jacobian is a conservative vector field. Specifically, we note that for any piecewise smooth loop $\\mathcal{C}\\_{\\text{loop}}$, we have $\\left \\Vert \\oint\\_{\\mathcal{C}\\_{loop}} \\mathbf{J} ds \\right \\Vert_2 = 0$. Thus, by integrating along loops that contain directions _orthogonal_ to the system's dynamics (and penalizing the deviation from zero), we encourage the estimated Jacobians to capture information about other directions in state space (see appendix A.3 for full technical details). To ensure broad coverage of tangent space directions, we form loops from concatenations of line integrals between randomly selected data points. This strategy samples diverse directions from the tangent space while remaining easy to compute. This self-supervised loss term builds on the loss introduced by Iyer et. al. [63]. It constrains $\\hat{\\mathbf{J}}^{\\theta}$ to satisfy both the dynamics and conservativity. This improves Jacobian estimation accuracy significantly (see appendix C.4 for ablation studies).\n\nObservation noise also helps the model learn responses in directions beyond the primary flow. We added small observation noise to each batch during training on the task-trained RNNs. As noted in the “Jacobian reconstruction quality” section: “Noise encourages the model to explore how perturbations around the observed trajectory evolve, which is crucial for learning accurate Jacobians in high-dimensional systems.”\n\nWe hope this addresses the reviewer’s question and are happy to clarify further.\n### Gramian discussion is too quick\n\nWe added the Gramian equation into the main paper, and agree it is a key quantity of interest. We also expanded the discussion of this quantity as follows:\n\n> Since Jacobians are time-dependent, we obtain a linear time-varying representation of control dynamics along the trajectory. This enables computation of time-varying reachability ease, capturing how readily each area drives the other toward novel states [56-58]. Reachability is quantified via the reachability Gramian, a matrix defining a local metric in the tangent space. For the above control system capturing the influence of area B on area A, the time-varying reachability Gramian on the interval $[t\\_0,t\\_1]$ is defined as $\\mathbf{W}\\_r(t\\_0, t\\_1) \\triangleq \\int\\_{t\\_0}^{t\\_1}\\mathbf{\\Phi}(t\\_1, \\tau)\\mathbf{B} (\\tau)\\mathbf{B}^T(\\tau)\\mathbf{\\Phi}^T(t\\_1, \\tau) d\\tau$ where $\\mathbf{\\Phi}$ (computed from $\\mathbf{J}^{A \\to A}$) denotes the state-transition matrix of the intrinsic dynamics of subsystem A without any input (i.e., $\\delta\\mathbf{x}^A(t) = \\mathbf{\\Phi}(t, t_0)\\delta\\mathbf{x}^A(t\\_0)$), and $\\mathbf{B}(\\tau) = \\mathbf{J}^{B \\to A}(\\mathbf{x}(\\tau))$ [56]. The Gramian $\\mathbf{W}\\_r(t\\_0, t\\_1)$ is symmetric and positive semidefinite for every $t\\_1 > t\\_0$ [56]. Each eigenvalue of the reachability Gramian quantifies how easily the target subsystem can be driven along its corresponding eigenvector. Thus, the trace of the reachability Gramian reflects average ease of reaching new states. The minimum eigenvalue indicates the lowest ease of control (i.e., the highest cost), and its corresponding eigenvector denotes the most difficult direction to control. For further details see appendix B.1.\n### Typo on line 150\nThis has been fixed.\n### Misleading title\nThank you for the comment. We’ve considered adding more neuroscientific context, and would appreciate your thoughts on this alternative: _Deep Jacobian estimation for characterizing control between interacting subsystems with applications to neuroscience_.\n\n### Transition to equation (1)\nEquation 1 has been moved to better separate it from Figure 1.\n### Relocating the related work section\nWe moved the related work section (section 6) to follow immediately after the introduction."}]

Author rebuttal:

Meta-review: In this paper, the authors propose JacobianODE, a method for directly estimating Jacobians from time series data using a path integral loss and a loop closure constraint. The framework enables nonlinear control-theoretic analysis of subsystem interactions, demonstrated on synthetic dynamical systems (Van der Pol, Lorenz, Lorenz-96) and a task-trained RNN. Results show JacobianODE outperforms Neural ODEs and weighted regression on Jacobian estimation error, Lyapunov spectrum recovery, and control tasks.

Reviews were overall positive. R-ya7v acknowledged the elegant theoretical formulation and applicability across nonlinear dynamical systems, and noted that the shift from vector field to Jacobian learning is well-motivated. R-DEPk highlighted the strong empirical evaluation and clear ablations, though suggested that section structure and exposition could be improved, and requested clearer discussion of the loop closure loss and more informative visualizations. Reviewers also asked for comparisons to echo state networks (ESNs) and clearer articulation of limitations. The rebuttal sufficiently addressed a lot of these points. The authors added ESN comparisons, showing that ESNs perform well on short-term prediction but fail to recover accurate Jacobians, underscoring the main improvement with their approach. They also expanded the limitations section, noting sensitivity of reachability analyses and challenges for scalability and partial observability. Clarifications were made around loop construction, delayed reachability estimates, and figure presentation. Reviewers were satisfied with the rebuttal, maintaining support and recommending acceptance after revisions.

Overall, the paper introduces a clear and well-supported methodological advance for data-driven Jacobian estimation and its application to control-theoretic analysis. I recommend acceptance, with the expectation that the authors integrate the added ESN comparisons, expanded limitations, and clarified presentation into the final version.
Decision: Accept (spotlight)