Tractable Dendritic RNNs for Identifying Unknown Nonlinear Dynamical Systems
Authors: Manuel Brenner
Leonard Bereska
Jonas Magdy Mikhaeil
Florian Hess
Zahra Monfared
Po-Chen Kuo
Daniel Durstewitz
Venue: ICLR 2022 Submitted (2022)
Keywords: chaos
dendritic computation
piecewise linear
recurrent neural network
variational inference
interpretability
tractability
basis expansion

Abstract: In many scientific disciplines, we are interested in inferring the nonlinear dynamical system underlying a set of observed time series, a challenging task in the face of chaotic behavior and noise. Previous deep learning approaches toward this goal often suffered from a lack of interpretability and tractability. In particular, the high-dimensional latent spaces often required for a faithful embedding, even when the underlying dynamics lives on a lower-dimensional manifold, can hamper theoretical analysis. Motivated by the emerging principles of dendritic computation, we augment a dynamically interpretable and mathematically tractable piecewise-linear (PL) recurrent neural network (RNN) by a linear spline basis expansion. We show that this approach retains all the theoretically appealing properties of the simple PLRNN, yet boosts its capacity for approximating arbitrary nonlinear dynamical systems in comparatively low dimensions. We introduce two frameworks for training the system, one based on fast and scalable variational inference, and another combining BPTT with teacher forcing. We show that the dendritically expanded PLRNN achieves better reconstructions with fewer parameters and dimensions on various dynamical systems benchmarks and compares favorably to other methods, while retaining a tractable and interpretable structure.

TRACTABLE DENDRITIC RNNS FOR IDENTIFYING UN- KNOWN NONLINEAR DYNAMICAL SYSTEMS
Anonymous authors Paper under double-hlind review
ABSTRACT
In many scientific disciplines, we are interested in inferring the nonlinear dynamical system underlying set of observed time series, a challenging task in the face of chaotic behavior and noise. Previous deep learning approaches toward this goal often suffered from a lack of interpretability and tractability. In particular, the high- dimensional latent spaces often required for a faithful embedding. even when the underlying dynamics lives on a lower-dimensional manifold, can hamper theoretical analysis. Motivated by the emerging principles of dendritic computation, we augment dynamically interpretable and mathematically tractable piccewise-lines -inear (PL.) recurrent neural network (RNN) by a linear spline basis expansion. We show that this approach retains all the theoretically appealing properties of the simple PLRNN, yet boosts its capacity for approximating arbitrary nonlinear dynamical We introduce two frameworks for
systems in comparatively low dimensions. training the system, one based on fast and scalable variational inference, and another combining BPTT with teacher foreing. We show that the dendritically expanded PLRNN achieves better reconstructions with fewer parameters and dimensions on various dynamical systems benchmarks and compares favorably to other methods, while retaining a tractable and interpretable structure.
INTRODUCTION
For many complex systems in physics, biology, or the social sciences, we do not know or have only rudimentary knowledge about the dynamical system (DS) that may underlie those quantities that we can empirically observe or measure. Data-driven approaches aimed at automatically inferring the generating DS from time-series observations could therefore strongly support the scientific process, and various such methods have been proposed in recent years (Raissi et al., 2018; Zhu al., 2021; Yin et al., 2021; Norcliffe et al., 2021; Mohajerin & Waslander, 2018; Karl et al., 2017; Chen et al., 2018; Strauss, 2020). However, due to the often high-dimensional, complex, chaotic, and inherently noisy nature of real-world DS, like the brain, weather-, ecosystems, this remains formidable challenge. Moreover, although the true DS may evolve on lower-dimensional manifold its state space, the system used for approximation usually needs be of higher dimensionality achieve proper embedding (Takens, 1981; Sauer et al., 1991; Kantz, & Schreiber, 2004). This is especially true when the approximating system is of a different functional form than the one that would most naturally describe the data generation process (but is unknown), for instance, when we attempt to approximate a system of exponential or trigonometric functions by polynomials.
In this work we sought to improve the capacity and expressiveness of a specific class of recurrent neural networks (RNNs), achieving agreeable solutions with fewer dimensions and parameters while retaining set of desirable theoretical properties. Specifically, we build on piecewise-lines (PLRNNs) based on ReLU activation functions, for which fixed points, periodic orbits, and other dynamical properties can be derived analytically (Schmidt et al., 2021; Koppe et al., 2019), and for which dynamically equivalent continuous-time (ordinary differential equation, ODE) systems can be constructed (Monfared & Durstewitz, 2020b). Inspired by principles of dendritic computation in biological neurons (Fig. 1), each PLRNN unit was endowed with a set of nonlinear pre- preprocesssng subunits ("dendritic branches"), such that it effectively takes the role of an equivalent much larger network. Mathematically, this comes down, in our case, to enhancing each latent unit with linear spline basis expansion as popular in statisties (Hastie et al., 2009). Through this trick, we achieve
RNNs

---

powerful RNN which provides reconstructions of underlying nonlinear DS in lower-dimensionn latent spaces than were needed by conventional PLRNNs. At the same time, model inference can be performed within the scalable framework of sequential variational auto-encoders (SVAE) (Archer et al., 2015; Girin et al., 2020; Krishnan et al., 2017), or with classical Pack-Propaggtionn TTrotihh Time (BPTT; Rumelhart et al. (1986)) augmented by teacher-forcing (TF; Williams & Zipser (1989); Pearlmutter (1990)). We further prove that these modifications preserve the mathematical and dynamical accessibility of the resulting system, e.g., such that fixed points, cycles, and their stability, can still be computed analytically.
Besides its effectivvness in capturing complex dynamical systems in fewer dimensions within a tractable framework, our approach highlights more generally how principles of dendritic signal processing may be harvested in the design of RNNs. Strongly nonlinear local computations are known for decades to occur within dendritic trees of biological neurons (Mel, 1994; Poirazi et al., 2003), but have hardly been exploited so far for machine learning models.
2 RELATED WORK
One class of DS reconstruction models attempts to discover governing equations from the flow field estimated from data through differencing the time series. Sparse Identification of Nonlinear Dynamics (SINDy), for instance, does so by sparsely regressing on rich library of basis functions using the least absolute shrinkage and selection operator (LASSO) (Brunton et al., 2016; Rudy et al., 2017; de Silva et al., 2020). Other methods approximate the flow field using graph reconstruction via differential equations (Chen et al., 2017), sparse autoencoders (Heim et al., 2019), shallow multi-layer perceptrons reformulated as RNNs (Trischler & D' E'leuterio, 2016), or deep neural networks (Chen et al., 2018). Some works aimed at directly learning the system's 's underlying Hamiltonian (Chen et al., 2020; Greydanus et al., 2019). Generally, numerical derivatives obtained from time series tend to be more noise-prone than the time series observations themselves (Baydin et al., 2018; Chen et al., 2017; Raissi, 2018). This can be problem particularly if only comparatively short trajectories were empirically observed or when the underlying systems are very high-dimensional, in these cases the system's 's flow field may be (severely) under-sampled. Methods directly based on numerical derivatives also need be augmented by other techniques, like delay embeddings (Kantz & Schreiber, 2004) or deep auto-encoders (Champion et al., 2019), if not all the system's 's dimensions were observed.
Various RNN architectures such as Long-Short -eee--Shoor--------e-e Term-Memooe networks (LSTMs) (Zheng et al., 2017), Reservoir Computing (RC) (Pathak et al., 2018), or PLRNNs (Koppe et al., 2019; Schmidt et al., 2021) have been employed to infer DS directly from the observed time series without going through numerical derivatives. More recently, transformers (Shalova & Oseledets, 2020a;b) were used as black box approaches for DS prediction. Except for PLRNNs, however, all these systems, although optimized for DS reconstruction and prediction, rest on relatively complex model formulations that are not easy to tackle and analyze from a DS perspective (Fraccaro et al.). The ability to gain deeper insights into the specific DS properties and mechanisms of the recovered system is, however, often crucial for its applicability to science and engineering problems. Transformers, unlike RNNs, do not even constitute DS themselves (as they explicitly forgo any temporal recursions), and therefore not directly amenable to DS theory tools. Moreover, most of these models, RC in particular, need very nigh-dimensiooal latent spaces, which further adds their black-box nature.
Better interpretability and tractability is achieved by using PLRNNS (Koppe et al., 2019; Schmidt al., 2021) or by (locally) linearizing nonlinear systems through ideas from Koopman operator theory (Azencot et al., 2020; Brunton et al., 2017; Yeung et al., 2017). In such systems, certain DS properties can be analytically accessed (Schmidt et al., 2021; Monfared & Durstewitz, 2020a), or the resulting equations can be more easily interpreted by a human reader (Heim et al., 2019). On the downside, usually one needs to move to very high dimensions to represent the DS in question properly. Here aim to overcome this limitation by augmenting PLRNNs with linear basis expansions without altering their analytical accessibility.
Finally, probabilistic (generative) latent variable models such state space models have been applied to the problem of posterior inference of latent state paths 2₁ ~ (( of DS given time series observations { }or (Pandarinath et al., 2018; Ghahramani & Roweis, 1998; Durstewitz, 2017; Krishnan et al., 2017). The advantage here is that they also account for uncertainty in the model formulation or latent process itself and yield the full distribution over latent space variables (Karl
2

---

et al., 2017). For DS reconstruction, however, need to move beyond posterior inference: We require that samples drawn from the model's ' prior distribution p(z) after training exhibit the same temporal and geometric structure those produced by the unknown DS.
Here we embed PLRNNS augmented with linear spline expansion into a fully probabilistic,, varia- tional approach that scales well with system size by employing stochastic gradient variational Bayes (SGVB; (Kingma & Welling, 2014; Rezende et al., 2014)), thereby combining the advantages of the two classes of models reviewed above. On the other hand, we show that the model can also be efficiently trained by BPTT using a specific form of TF (Appx. 6.1).
3 MODEL FORMULATION AND THEORETICAL CONSIDERATIONS
3.1 PIECEWISE LINEAR RECURRENT NEURAL NETWORK (PLRNN)
Our approach builds on PLRNNs (Durstewitz, 2017; Koppe et al., 2019) because of their mathematical tractability (see Sec. 3.3). PLRNNS are defined by the M M-dimensional latent process equation
21 A + W Wo(zz-1) + + Cst + 6,,
which describes the temporal evolution of -ddimensioonal latent state vector 2; = (Zie 2M =M.)". The self-connections of the units represented by diagonal matrix A E RM whereas the connections between units are collected in off- matrix W E RMMM with the nonlinear activation function & given by the rectified linear unit (ReLU) applied element-wise:
((2--1) = max 21-1).
(2)
Additionally, the PIRNN comprises a bias term h E M potential external inputs 8, E RRKK weighted by C E R R and Gaussian noise term ~ N N0, Σ) with diagonal covariance E. The PLRNN can be interpreted as a discrete-time neural rate model (Durstewitz, 2017), where the entries of A stand for the individual neurons' time constants, W for the synaptic connection strengths between neurons, and d(2) for a (ReLU-shaped) voltage-to-spike-rate transfer function. The probabilistic latent RNN Eq. 1 is linked to the N N-dimensional observed time series (æ1))=1...7 Dy E IR,, drawn from an underlying noisy DS, by an observation function (decoder model) which, in the simplest case, may take the linear Gaussian form
32 Bz + Th:
(3)
where B E N represents factor loading matrix and m ~ NO F) is Gaussian observation noise with diagonal covariance r E IR N
3.2 DENDRITIC COMPUTATION AND SPLINE BASIS EXPANSION
Dendrites have long been known to play an active and important part in neural computation (Mel, 1994; 1999; Koch, 2004). Active, fast voltage- gated ion channels endow dendrites with strongly nonlinear behavior, giving rise for instance to dendritic Ca²+ spikes that boost synaptic inputs (Schiller et al., 2000; Häusser et al., 2000). It has been suggested previously that different dendritic branches may constitute rather independent computational sub-units whose outputs are combined at the soma, as in a 2-layer neural network (Poirazi et al., 2003; Mel, 1993; 1994). an idea that received strong empirical support especially in recent years (Poirazi & Papoutsi, 2020). Here we mimic this functional setup by modeling dendritic processing through linear combination of ReL.U-type threshold-nonlinaaride (Fig. 1), replacing Eq. 2 by
(2x-1)= = Σ og max (0, 2
-ha),
with "dendritic input/output" slopes an R and "activation" thresholds hs E RM As in real dendrites, where both ion channels and morphological structure are subject to learning (Poirazi & Papoutsi, 2020; Stemmler & Koch, 1999), we treat these as trainable parameters. We note that Eq. inserted into model Eq. takes the form of a linear spline basis expansion popular in statistics (Hastic et al., 2009) for approximating arbitrary functions (Wahba, 1990; Storace & De Feo, 2004) in regression settings. For instance, such concepts have heen frequently employed within data-analytical
3

---

models in neuroscience (Frank et al.; Huang al.; Qian al.), but never within the context of DS reconstruction enabling liver-dimensional solutions in mathematically tractable models.
To emphasize the connection to dendritic computation call the system Eqs. 1, 4, the dendPLRNN.
Q
Figure I: Inspired by principles of dendritic computation, our dendPLRNN extends each unit into a set of nonlinear branches connected soma, yielding single unit transfer functions with increased approximation capabilities. Image of dendrite from https /www, milad no no/blog/ (CC BY- SA 4.0).
3.3 MATHEMATICAL TRACTABILITY AND DYNAMICAL SYSTEMS INTERPRETATION
Sharp threshold-nonlineeaitis (like a ReLU) are a reasonable choice from neurobiological per- spective, dendrites naturally give rise this threshold-type behavior (Mel, 1999; Koch, 2004). Another important consideration in choosing this particular form, however, was that it preserves all the theoretically appealing properties of a PLRNN, we will formally establish below: For PLRNNs fixed points and cycles can be explicitly computed (Schmidt et al., 2021; Koppe et al., 2019). and they can be translated into dynamically equivalent continuous time systems (Monfared & Durstewitz, 2020b), properties which profoundly ease the analysis of trained systems from DS perspective. This crucial for application the sciences, where we are specifically interested in understanding the underlying system's ss dynamics. For PLRNNs, precise connections between the long-term behavior of the system and that of its gradients have also been established (Schmidt et al., 2021). Finally, PLRNNs belong to the class of continuous piecewise-linee (PWL) maps, for which many important types of bifurcations have been well characterized (Feigin, 1995; Hogan et al., 2007; Patra, 201 8) (ef. (Monfared & Durstewitz, 2020a) for overview). Bifurcations essential understand how geometrical and topological properties of the system's 's state space depend on its parameters or could be controlled, and hence are also important characterize improve the training process itself (Doya, 1992; Pascanu et al., 2013; Saxe et al., 2014) or to understand properties of trained systems (Maheswaranathan et al., 2019b;a).
Our first proposition, therefore, assures that by the particular form of basis expansion introduced in Eq. 4, the system will remain within the class of continuous PWL maps:
Proposition The model defined through Eg. and Eg. constitutes continuous PWI. map.
The proof essentially straightforwardly follows from the model's definition as linear spline basis expansion in each unit, but is formally provided in Appx. 6.5.4.
While Proposition is all we need to ensure we can harvest all previously established results on PLRNNs in particular, and continuous PWL maps more generally, it is revealing note that any dendPLRNN (Eqs. 1, 4) can be rewritten as a conventional PLRNN, as stated in the following theorem:
Theorem 1. Any M-dimensional dendPLRRN as defined in Eqs. can always be rewritten M x B-dimensional "conventional" PLRNN of the form
a = - max(0, 2x-1) + ho + Cs + ēt
(5)
Proof. Straightforward by construction, see Appx. 6.5.5.
4

---

This theorem highlights why the dendPLRNN will allow to reduce the dimensionality of the recon- structed system, as it suggests we may often be able to reformulate a figh-dimensional PLRNN in terms of equally powerful lower- dimenssooal dendPLRNN. In Appx. 6.5.1 we also spell out the exact computation of fixed points and k-cycles for the dendPLRNN.
Finally, the unboundedness of the PLRNN's s latent states due to the ReLU function cause divergence problems in training. The dendPLRNN, on the other hand, offers a simple and natural way contain the latent states without violating the basic model description above, as established in the following theorem:
Theorem 2. For each basis {as, h in Eq. 4 of a dendPLRNN let us add another basis {at, h with og aas and h O. Then, for max ()) 1, any orbit of this "clipped" dendPLRNN (Eq. 10) will remain bounded.
Proof. See Appx. 6.5.6.
Appx. 6.5 collects further theoretical results, assuring, for instance, that the manifold attractor regularization employed here (see next section) does not interfere with the results above (Proposition 2).
3.4 TRAINING THE DENDPLRNN
To infer the parameters 0 = (A, W. h, C,E, {ab, h2]} of the dendPLRNN, (Eq. 1, 4) from observed data, we apply two different training strategies: First, fast and scalable variational inference (VT) algorithm which maximizes the Evidence Lower Bound (ELBO) C(O, di x) == Egg [[]]()
KL (x) [[)]=)] using the reparameterization trick (Kingma & Welling, 2014), and convolutional neural networks (CNNs) for parameterizing the encoder model qa (=] (see Appx. 6.1 for details). Furthermore, as proposed in Schmidt et al. (2021), to efficiently capture DS multiple time scales, we add regularization term the ELBO that encourages the mapping of slow time constants and long-range range dependencies (so-called manifold attractor regularization, see Eq. 6, with regularization factor X). Second. we employ "classical" BPTT with variant of teacher forcing (TF) (Williams & Zipser, 1989; Pearlmutter, 1990). TF here means that the first N latent states 1,,/++1, ls S N, were replaced by observations J at times + 1, e No, where T 2 is the forcing interval (for details, see Appx. 6.1). All code used in here is made freely available at [placeholder).
EXPERIMENTS
4.1 PERFORMANCE MEASURES
In DS reconstruction, we aim capture invariant properties of the underlying DS like its geometrical and temporal structure. To evaluate the quality the reconstructions w.r.t. geometrical properties employed a Kullback-Leibler divergence (Dsss) that quantifies agreement in attractor geometries (more details in Appx. 6.2), as first suggested Koppe et al. (2019) (see also Schmidt et al. (2021)). Specifically, this measure evaluates the overlap between the observed data distribution p(a:"h)) and the distribution p 2 generated from model simulations (i.e., with za0 ~ Po after model training) across state space (not time!). Since this measure as originally defined in Koppe al. (2019) is expensive to compute, for the high-dimensional benchmark DS we used another approximation, details of which are given in Appx. 6.2. Dasp evaluated on set of 100 trajectories, pulled from the learned distribution over initial conditions, with 1000 time steps each. To assess the agreement in power spectrum correlation (PSC)
temporal structure, dimension-win wise, Gaussian-kens www.ek-eeeooteeeee between ground truth and model-generate trajectories was used (see Appx. 6.2). Finally, we also computed 20-step-ahead prediction error along test set trajectories (see Appx. 6.2), although not of primary interest in the context of DS reconstruction.
4.2 DS BENCHMARKS USED FOR EVALUATION
We evaluated our approach and the specific role the basis expansion on five different types of challenging DS benchmarks.
5

---

First, the famous 3d chaotic Lorenz attractor (Lorenz-63) originally proposed by Lorenz (1963) (formally defined in Appx. 6.4) has become a popular benchmark for DS reconstruction algorithms. Fig. 2a (l.h.s.) illustrates true (blue) and reconstructed (orange) time series from this system, while the r.h.s. illustrates the chaotic attractor's geometry in its 3d state space for both the ground truth (blue) and reconstructed (orange) systems. It is important to note that both the time and state space graphs are not merely ahead predictions from the dendPLRNN but are produced by simulating the trained dendPLRNN from some initial condition. This illustrates that the dendPLRNN has captured the temporal and geometrical structure of the original Lorenz-63 system in its own governing equations. Moreover, computing analytically (see Appx. 6.5.1) the fixed points of the reconstructed system, we see that their positions in state space agree well with those of the true system.
Second, 3d biophysical model of a bursting neuron (see Eq. 15 in Appx. 6.4; Durstewitz (2009)) highlights another aspect of DS reconstruction: Besides an equation for membrane voltage (V), the model consists of one very fast (n) and one slow (h) variable that control the gating of the model's ionic conductances. This produces fast spikes that ride on top of much slower oscillation, making this system challenging to reconstruct. One such successful dendPLRNY reconstruction is illustrated in Fig. 2b (orange) together with time graphs and state space representations of the true system (blue).
(a))
true generated
will
1
2000 2000
5000
250
1000
time steps
time steps
Figure Examples of low-dimensional model reconstructions: (a) Time series (left) and state space trajectories (right) for the original Lorenz-63 chaotic attractor and simulations produced by a dendPLRNN trained with VI (B 20, M = 15, A 1, Mreg/ /M = 0.5). Dots indicate true and reconstructed fixed points. (b) Same for the bursting neuron model, produced by dendPLRNN
trained with TF (B = 47, M = 26, t = 5). Note that the bursting is a complex limit cycle but non-chaotic.
Third, the Lorenz-96 weather model is an example of higher-dimensionall spatially organized chaotic system with local neighborhood interactions that can be extended to arbitrary dimensionality (Eq. 18 in Appx. 6.4). It has also been used more widely for benchmarking DS reconstruction algorithms. For our experiments we employed 10-dimensional spatial layout. Fig. 3a illustrates time graphs for selected dimensions (l.h.s)), the full evolving spatio-temporal pattern (center), and examples of power spectra (r.h.s.) for both the ground truth system (blue) and an example reconstruction (orange). The spatio-temporal characteristics of the true and the SoaPLLNNN--enerrtee time series tightly agree.
Fourth, as another high-dimensional example we used neural population model with structured connectivity tuned to produce coherent chaos (Landau & Sompolinsky, 2018), from which we produced 50d observations (see Appx. 6.4 for details). Fig. 3b provides example time series (I.h.s.), full spatio-tempora patterns (center), and overlaid power spectra (r.h.s.) for time series drawn from the true system (blue) and those simulated by trained dendPLRNN (orange). Again there is a tight agreement, and again we emphasize that - like all the other examples these are not mere model ahead-predica preddccioiiissonss but fully simulated from some random initial condition.
Finally, we studied real-world dataset consisting of electroenceehaloggrm (EEG) recordings from human subjects, described in more detail with results (Fig. **8) in Appx. ***
4.3 BASIS EXPANSION IMPROVES DS RECONSTRUCTION AND LOWERS PARAMETER COSTS
While, in theory, the dendPLRNN is equivalent to a larger PLRNN without basis expansion, in practice the smaller basiseex models trained more successfully. Fig. 4a summarizes observations for the VI algorithm on the impact of the basis expansion using the Lorenz-63 DS example (see Fig. S3 for further examples): Both the 20-step ahead prediction emor well

---

true
generated
true generated
time steps
time steps
dimensions
frequency
[Hx]
true
generated
time
time steps
dimensions
20
frequency [Hz]
Figure 3: Examples of high-dimensional model reconstructions: (a) Time series (left), spatio-temporal
evolution (center), and power spectra (right) for the true 10d Lorenz-96 system and for dendPLRNN
simulations (B = 50, M = 30, À 1.0, Mreg/ /M = 1.0). (b) Same for 50d neural population
model producing coherent chaos (B 5, M 12, x 1.0, Mreg/ M 0.2).
(a)
D Dstsp
20-step PE
(b)
Dstsp
20-step PE
12.5
10.0
Paramsters
206
10-1
Penemeters
530
220
14
10-
M
7.5
10
H
5.0
23
18
10
0.0
10
number bascs
number bases
number bases
number bases
Figure Effect of basis expansion for dendPLRNN trained by VI. (a) Agreement in attractor geometries (left) and 20-step ahead prediction error (right) for the Lorenz-63 system as a function of the number of bases (B) for fixed numbers of total parameters. (b) Agreement in attractor geometries (left) and 20-step ahead prediction error (right) for the Lorenz-63 system as function of the number of bases (B) for different numbers of latent states (M).
as the geometrical reconstruction quality as assessed by Daup profoundly improve with the size B of the basis expansion even for the same total number of trainable model parameters (given by (M + + B + N) + B, where is the dimensionality of the observed data). Hence, as conjectured in Sec. 3, the basis expansion yields better reconstructions at no additional costs in terms of numbers of model parameters. Fig. 4b looks at the impact of the basis expansion from the angle of dimensionality reduction by systematically varying the number of bases B and latent states M for the Lorenz-63 DS: Following the curves horizontally, it becomes clear that the basis expansion enables to reduce the model's overall dimensionality without compromising performance.
4.4 MODEL COMPARISONS
We compared our model the PLRNN without the dendritic expansion and three other algorithms purpose-tailored for DS reconstruction: First, SINDY (Brunton et al., 2016) aims to reconstruct the governing equations by approximating numerical derivatives (obtained by differencing the time series, and applying a variance regularization to reduce noise) through large library of polynomial basis functions. Sparse (LASSO) regression is used to pick out the right terms from the library (we used the PySINDy implementation (de Silva et al., 2020) with multinomials up to sixth order). Second, Vlachas et al. (2018) used hybrid of truncated LSTMs and mean-field stochastic models based on Ornstein- Ubeeenbeckheehckkkkk processes ILSTM-MSM observed time series. Third, Pathak et al. (2018) built on reservoir computing (RC) for their approach with reservoir parameters chosen satisfy the "echo state property" (Jaeger & Haas, 2004). For
approximate the true system's 's flow estimated from
7

---

Table 1: Comparison of dendPLRNN (Ours) trained by VI or BPTT+TF, RC (Pathak et al., 2018), LSTM-MSM (Vlachas et al., 2018), and SINDy (Brunton et al., 2016) on 4 DS benchmarks and one experimental dataset (top) and 3 challenging data situations (bottom). Values are mean = SEM.
Dattect/ / Setting
Meshoal
PSC
Dpp
20-eep PE
Dynamical vansills
Pertuctere
deudPLRNN VI
0.397
0,111
0,50 2: 0.25
2.1e
3 2 0.2c 3
22
10:2
deadPLRNN TP
0.957
0.102
0.13
2: 0.18
9.2
2:80
5
22
10:22
Locenz
RC
0.991
0.001
0.24
1: 0.05
1.
2
1
0.1--
345
1053
LSTM MSM
0.985
0.004
0.85
1: 0.07
1.3:-11
1
0.e--
20
1035
SINDY
0.998
0.0003
0.04
1: 0.01
6.8c-5
a: 0.2e-5
252
JasPLRNN VI
0.967
0.001
0.10
+: 0.01
3.1e-1
0.0c-1
42
4361
Lovene
96
LSTM MSM
RC
JauPPLRNN TF
0.998
0.986
0.988
0.102
0.008
0.0001
0.04
0.25
0.23 1: 0.08
2: 0.17
0.01
4.1e-2
8.2c
7.le-
2 0.1e-2
+ 0.8c-2
2 0.3c 2
62
50
440
4400
425.1
4450
SINDY
0,956
0.001
0,06
4: 0.003
6.3e
2 4 0.10
3
10
27410
dedPi
VI
0:35
0.03
7.5
6.1e-1 010--
70
20533
Nerrvon
Buntiig
dendPP RNN
LSTM-MMM
0.76
0.51
0
0.01
0.04
2.0
5.1 0.6
0.36
8.9c-2
5.0-2 0.1--2
6.1c-2
2 0.1e-2
2.6-2
70
4G
2166
21.33
20:10
SIND>
diverping
diverging
diverging
Neural
Popula
GendPLRNN
dendPLRNN
0.51
0.02
0.09
75
12
LSTM-MMM
SINDy
disrging
diveging
divering
JeuPPLRNN VI
0.01
27.9
0.046
2710%
deadPLRNN TF
0.33-
0.017
4.7
0.013
27038
BBG
RC
0.51
0.01
5.4
28672
LSTM MSM
0.005
20 0.
27728
SINDY
diverging
diverging
diverging
1351
JeuIPLRNN VI
0.967
0.007
0.10
2.
0.2c-3
10.52
Low
JendPLRNN TF
0.97 0.04
6.9
0.9c
22
10122
0.58
0.15
5.74
4.1415
1.26-
1053
0,980
0.006
2.1e
0.36
1035
SINDY
0.0003
0.03
6.8e
0.2e-
JNNN
1-4-2
Panially
JauPPLRN TF
0.1--5
22
1053
1052
LSTM-M5M
0.334
0.3c-2
29
1055
0,374
0:46-2
dendPLANN V]
0,003
0.006
4.9 0.75
3.5e-2
0.1e-2
22
1052
High
dendPLKNNN TF
RC:
0.068
0.001
2.53 0.21
0.4 0.13
3.1e-2 0.2e-2
4.0--3
0.4e-5
345
22
1053
1052
STM-MMMM
0.957
0.000
1.19 + 0:27
3-3-22 0.2c-2
79
1055
SINDY
0,981
0,006
1.01
0.05
2.3-3 0.1e-4
252
higher-dingg -ibbee-iimnniinal systems, spatially arranged set of reservoirs with local neighborhood relations employed. For all these systems, optimized hyper-pane spper-parametess were used as reported by the authors. For our own system, the dendPLRNN, we also performed a grid search for optimal hyper-parament Areg. Trr. M, and B (see Appx. 6.1 and Table $1 for details). For all four methods, to the degree possible we tried to ensure roughly the same number of trainable parameters (see Table S2).
Results for all four models on all five DS benchmarks employed here are summarized in the upper part of Table S2, using the temporal and geometrical reconstruction measures introduced in Sec. 4.2 (as well as a 20-step-ahead prediction error for comparison). To produce this table, 1 100,000 time steps for both training and testing were simulated from each ground truth system, all dimensions were standardized to have zero mean and unit variance, and process noise and observation noise (with 1 1% of the data variance) were added while simulating the (now stochastic) differential equations, and after drawing the observations, respectively (see Appx. 6.4 for further methodological details). To produce statisties, each method was run from total of 20 randomly chosen initial conditions for the parameters. We also tested all four methods on the real EEG data and on challenging data situations produced using the Lorenz-63 system (Fig. 2a), with either short time series of just 1000 time steps, only partial observations (just state variable r in Eq.14 14 in Appx. 6.4), or high process and high observation noise (drawing from Gaussian with de (O, 0.1dt x I) for the process noise as described in Appx. 6.4, and using 10% of the observation variance, respectively). SINDy cannot naturally handle missing observations, as it has no latent variables but formulates the model directly in terms of the observations. Therefore, for the partially observed system, we used delay embedding (Takens, 1981; Sauer et al., 1991) to create a 3d dataset, adding two time-lagged versions of x as coordinates.
A general observation is that indeed all four models are quite powerful for reconstructing the underlying DS. However, in most comparisons the dendPLRNN had an edge over the other methods,
We point out that this may already impose restriction for methods like SINDY as moves to very high-dimensional systems.
8

---

came out second after SINDy, especially when trained by BPTT+TF. SINDY tends to outperform the dendPLRNN on the Lorenz-63 DS, but it completely fails on the bursting-neuus neuron and population model examples, and on the real EEG data, and generally becomes comparatively slow to train on high-dimensional systems. This can he explained by the fact that SINDY already has the correct functional form for the Lorenz-63 (and also Lorenz-96) DS: Both of these have strictly polynomial form (see Eq. 14 and Eq. 18 in Appx. 6.4), and SINDY works with a set of polynomial library functions begin with. Hence, SINDy only needs to pick out the right terms from its expansion to succeed, giving it a clear advantage on these model systems by design. On the other hand, as revealed in Table S2, completely fails on systems which have a different (non-poly polyoolooomil) functional form (or when the true form, as in the EEG case, is simply not known). SINDy therefore appears less suitable as a general framework for DS reconstruction if an appropriate library of basis functions cannot be specified priori, unlike the other methods.
While our conclusion is that essentially all of the three tested models LSTM-MSM, RC, and dend- PLRNN are suitable for reconstruction of arbitrary DS even in very challenging data situations (Table S2, bottom), LSTM-MSM and RC performed worse on average and have other profound disadvantages compared to our method: First, they are quite complex in their architectures and hence not easily interpretable, i.e. much harder to track and analyze mathematically." In contrast, as summarized in Sec. 3.3, the dendPLRRN is continuous PWL map and such comes with huge bulk of already existing theoretical results (Schmidt et al., 2021; Monfared & Durstewitz, 2020b;a), as well as with mathematical tractability (see Fig. 2a and Appx. 6.5.1). On top, the dendPLRNN achieves reconstruction of all DS in (much) lower dimensions than the RC or LSTM-MSM (see Table S2), further adding its better interpretability. Second, by embedding the dendPLRNN within SVAE (Archer et al., 2015) framework we also obtain uncertainty estimates on the state trajectories and can perform posterior inference, features that the other models lack.
5 CONCLUSIONS
In this work augmented PLRNNs (Durstewitz, 2017; Koppe et al., 2019) by a linear spline basis expansion inspired by principles of dendritic computation. We show mathematically that by doing so remain within the theoretical framework of continuous PWL maps and hence can harvest huge bulk of existing DS theory (Sec. 3.3), while at the same time achieving better performance with less parameters and in lower dimensions. Another contribution of this work is transferring the PLRNN into the framework of SVAEs which allow for fast and scalable inference and training. These are two key advantages from both a scientific perspective where mechanistic insight and understanding of the system under study are sought, and prediction perspective where we are also interested in uncertainty estimates.
We close by pointing out two open issues: First, somewhat surprisingly. the BPTT+TF approach to model training clearly outperformed the more sophisticated VI approach. This could be rooted in suboptimal encoder models or in suboptimal sampling from the approximate posterior: While BPTT+TF assesses longer bits of trajectory during optimization, in VI single time-point samples are drawn and the temporal consistency is ensured only through the Kullback-Leible term in the ELBO. Other more expressive yet still fast to compute encoder models, e.g.. based on normalizing flows (Rezende & Mohamed, 2015), may boost performance. Smart initialization techniques (Talathi & Vartak, 2016) or specific annealing and curriculum training protocols (as used in Koppe et al. (2019)) are other amendments to consider. Second, we felt that quantitative measures for assessing the quality of DS reconstructions in high-dimensional, high-noise situations are an interesting research topic in their own right. It is known that "classical" DS measures like Lyapunov spectra or correlation dimensions (Kantz& & Schreiber, 2004) are very hard to robustly assess for higher- dimees-iimnsi or more noisy systems (Schreiber Kantz, 1996), and are often not even known for comparatively simple models. Yet, the geometrical and temporal measures employed here come with their own pitfalls, some of them alluded to in Appx. 6.2.
ACKNOWLEDGMENTS
"This is especially true for RC. Moreover, the fact that only the weights of the linear output layer are trainable while the recurrent connections within the reservoirs are static, may raise the question of what precisely is learnt in terms of dynamics if the reservoirs themselves cannot adapt the DS at hand.

---

Ethics statement The current work performs theoretical analysis and basic research novel type of RNN architecture, which is evaluated exclusively on a set of simulated physical and biological systems (no human or animal subjects involved, privacy concerns etc.). Although it is conceivable that the current models, training algorithms, and results may be used to develop prediction algorithms in sensitive domains (like medical time series tracking consumer data), there ethical implications from this work as far as we can see.
immediate
Reproducibility statement All theoretical results this paper were carefully and thoroughly proven, with all proofs and detailed derivations available in the Appendix. Likewise, we will make available all code used in the empirical section in a way that will allow others to easily reproduce the results from this paper. This means we will include everything, starting with the code for the benchmark models and simulations, the simulated time series data used for evaluation themselves, the code for our model and training algorithms, up to the meta-files that produce the figures in this work, on our lab github site. All of this will be clearly documented.
REFERENCES
Evan Archer, Il Memming Park, Lars Buesing, John Cunningham, and Liam Paninski. Black box variational inference for state space models. arXiv preprint arXiv:1511.07365 2015. URL
http://a /ttes//rxxiii v. org/ org/abs////eooo.oorr.ggoo/eeb /1511.111.0000
Omri Azencot, N. Benjamin Erichson, Vanessa Lin, and Michael W. Mahoney. Forecasting Se- quential Data using Consistent Koopman Autoencoders. In Proceedings of the 37th International Conference on Machine Learning, 2020. URL http: / /ttxxviivvv//////////////////av org/ abs 2000 02236.
Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey Hinton. Layer Normalization. arXiv:1607.066. 0.6.00655 [cs, stat], July 2016. URL http:://aaae org/ abs 1677. 06450. version: 1.
arXiv: 1607.06450
Atilim Gunes Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic Differentiation in Machine Learning: a Survey, Journal of Machine Learning Research, 18(153):1-43, 2018. ISSN 1533-7928. URL http: /t////mmmi jmmr. org/ papers /v18/ 8/ 17 468. html.
Justin Bayer, Maximilian Soelch, Atanas Mirchev, Baris Kayalibay, and Patrick van der Smagt. Mind the gap when conditioning amortised inference in sequential latent-variable models. URL
http:// arxiv. wttp:///g/b/bb/aaae
David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859-877 Apr 2017. ISSN 1537-274X. doi: 10. (0800016214999211.112212
URL http: //tp doi. org / 10 1 080 02144599.
2017 1285773
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc Narl Acad Sei U S A, 113 (15):3932-3937, 2016. ISSN 0027-8424. doi: 10. 10.1073/paas.15111 1517384113. URL https: http:s///www// /www. ncbi mlm. nih gov pmc/ attiiattiiio htt.sccc/l///443/3/3
Steven L. Brunton, Bingni W. Brunton, Joshua L. Proctor, Eurika Kaiser, and J. Nathan Kutz. Chaos intermittently forced linear system. Nai Commun, 8(1):19, 2017. ISSN 2041- 1723. doi: 10.1033/44446--0100000333 s41467-017-00030-
URL http:// http:://ww.nnt wwwp:/wwww.w...at..tttttwttww..w wwtp:/aatttee.ee.aaeeeeeee www.articcee///////ooooo..
Kathleen Champion, Bethany Lusch, Nathan Kutz, and Steven L. Brunton. Data-driven discovery of coordinates and governing equations. Proc Natl Acad Sci USA, 16(445):2444-2222 2019. ISSN 0027-04244 1091-6490. doi: 19wp3//aa..9069511111 16. URL http:// http:://wwwwwppp.ppp//o/ma pnas org/ fookup/doi/ /10 073/0/3///m/ 1 6.99999555 006991116.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems 31, 2018, URL http://arxtt
httt:::/axxxi..or/abbb////// www..bs/////////eooooooo /180la//////.//////// 07366.
10

---

Shizhe Chen, Ali Shojaie, and Daniela M. Witten. Network Reconstruction From High-Dimensional Ordinary Differential Equations. Journal of the American Statistical Association, 112(520): 1697-1707, 2017. ISSN 0162-1459. doi: 10.
URL https:
///oii.or//0... 1080/000001111 621 159., .0..599..11500 20..6.666000 92.197.
Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic Recurrent Neural Networks. In Proceedings of the 8th International Conference on Learning Representations, 2020.
URL http::
org/ abs/1909 13334.
Zhicheng Cui, Wenlin Chen, and Yixin Chen. Multi-scale convolutional neural networks for time series classification. Computing Research Repository, abs/1603.0699 2016. URL http:
//arxiv. wwt.o////aaa... http///bb/10333 0.803.069900 06995.
Brian M. de Silva, Kathleen Champion, Markus Quade, Jean-Christophe Loiseau, J. Nathan Kutz, and Steven L. Brunton. PySINDy: A Python package for the Sparse Identification of Nonlinear Dynamics from Data. arXiv preprint arXiv:2004 2020. URL http: / /ttpii.xivvvr...//.. org/ abs/2004 08424.
Kenji Doya. Bifurcations in the learning of recurrent neural networks. In Proceedings of the 1992 IEEE International Symposium on Circuits and Systems, 1992. ISBN 978-0-7803-0593-9 doi: 10. R..119//SCAASS 230622/.
URL http:://tteeeeee /itp:////eeeeeoreeeee//eee ieee org//oument//
Daniel Durstewitz. Implications of synaptic biophysics for recurrent network dynamics and ac- tive memory. Neural Networks, 22(8): 2009. ISSN 08936080. doi: 10.
j.neumet.2099007.01 URL https: /tt:::///nnnnnnnhhi elsevier com/ wttp///erreeeee////// $0893608009001622.
Daniel Durstewitz. state space approach for piecewise-linear recurrent neural networks for identifying computational dynamics from neural measurements. PLoS Comput. Biol., 13(6): e1005542, 2017. ISSN 1553-7358. doi: 10. Ij..ra//jou......cb. 1005542.
Mark I Feigin. The increasingly complex structure of the bifurcation tree of a piecewise-smooth system. Journal of Applied Mathematics and Mechanics, 59(6) 863, 1995. ISSN 0021- 8928. doi: 10. 001101600011-8222995000000 Seience/article//iiiie.o
URL https: www. www.::sseeeneedeeeee.m.o com/
Marco Fraccaro, Søren Kaae Sønderby, UIrich Paquet, and Ole Winther. Sequential neural models with stochastic layers. URL http:/ arriv.orgg..............o.... abs /olss/6005 01166.07777. 07571.
Loren M. Frank, Uri T. Eden, Victor Solo, Matthew A. Wilson, and Emery N. Brown. Contrasting patterns of receptive field plasticity in the hippocampus and the entorhinal cortex: An adaptive filtering approach. 22(9):38173388)) ISSN 0270-6474, 1529-2401. doi: 10. 10.1RSS/IEUUOOCI. 22-09-03817.228 2002 URL https: /www. jneurosci............ org wwontent/22/9//00
Pub-
lisher: Society for Neuroscience Section: ARTICLE.
Zoubin Ghahramani and Sam T Roweis. Learning nonlinear dynamical systems using EM algorithm. Advances in Neural Information Processing Systems /1, 1998.
Laurent Girin, Simon Leglaive, Xiaoyu Bie, Julien Diard, Thomas Hueber, and Xavier Alameda- Pineda. Dynamical Variational Autoencoders: A Comprehensive Review. arXiv preprint arXiv:2008. 12595, 2020. URL http: http://aaavvvva///a/e...aaaa..aa org/ abs 22088 12595.
Sam Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian Neural Networks. In Advances in Neural Information Processing Systems 32, 2019. URL 01563.
Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The elements of statistical learning: data mining, inference, and prediction. Springer Science & Business Media, 2009.
Niklas Heim, Václav Šmídl, and Tomás Pevný. Rodent: Relevance determination in differential equations. arXiv preprint arXiv: 2019. URL http: /ttp/v/.ivvvvvvv///////..//...vv atb/ /1912. 00656.
11

---

John R. Hershey and Peder A. Olsen. Approximating the kullback leibler divergence between gaussian mixture models. 2007 IEEE International Conference on Acoustics, Speech and Signal Processing ICASSP '07, 411--11--V---222 2007.
Stephen J. Hogan, L. Higham, and T. C. L. Griffin. Dynamics of a piecewise linear map with a gap. Proceedings of the Royal Sociery A: Mathematical, Physical and Engineering Sciences, 463(7077): 49-65, 2007. doi: 10. 10.10000/pp...20....1 1735, URL https: ptte:s/://ooallocictpobosb
org/doi/ 10. 1098/
rspa. 2006 1735.
Su- Yun Huang, Yi-Ren Yeh, and Shinto Eguchi. Robust kernel principal component analysis. 21(11): 3179-3213, ISSN 0899-7667. doi: 10.11622meeo.20000022222 02-0000000-6666 Conference Name: Neural Computation.
Michael Häusser, Nelson Spruston, and Greg J. Stuart. Diversity and Dynamics of Dendritic Signaling. Science, 290(5492):739-744, 2000. ISSN 0036-8075, 1095-9203, doi: /t.11///cience.29055544 739. URL https: //cieecc/eceeeeeeeeeeeeeee / http://ccccnnnc.ccennnnnagg////.
Herbert Jaeger and Harald Haas. Harnessing Nonlinearity: Predicting Chaotic Systems and Sav- ing Energy in Wireless Communication. Science, 304(5667):78-80, 2004. ISSN 0036-8075, 1095-9203, doi: 10. 122/sccience110912777 URL content/ /////33///666///
htw.:/enccescciieeeeeee.o... org/
Holger Kantz and Thomas Schreiher. Nonlinear time series analysis, volume 7. Cambridge university press, 2004.
Maximilian Karl, Maximilian Soelch, Justin Bayer, and Patrick van der Smagt. Deep Variational Bayes Filters: Unsupervised Learning of State Space Models from Raw Data. In Proceedings of the 5th International Conference on Learning Representations, 2017. URL org/ har//a///aa /605. 06432.
Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. In Proceedings of the 3rd International Conference on Learning Representations, 2015. URL http: /arxiv. org/ /bs//ab//ta 1412/ 6980.
Diederik P. Kingma and Max Welling. Auto-Encoding Variational Bayes. In Proceedings of the 2nd International Conference on Learning Representations, 2014. URL abs/131. 6114.
Christof Koch. Biophysics of computation: information processing in single neurons. Oxford university press, 2004.
Georgia Koppe, Hazem Toutounji, Peter Kirsch, Stefanie Lis, and Daniel Durstewitz. Identifying nonlinear dynamical systems via generative recurrent neural networks with applications to fMRI. PLOS Computational Biology, (588) 2019. ISSN 1553-7358. doi: 10. 13..or//journal. pcbi. 1007263. URL https::t///cc / plos org. http:///pw.cccpppp///aaaacaeaaa.
10 1371/ wtwpaa//oooaaaaaaa www.rnal..cbbi... ptt.ppiii.....................1. 1007263.
Rahul G. Krishnan, Uri Shalit, and David Sontag. Structured Inference Networks for Nonlincar State Space Models. In 31st AAAI Conference Artificial Intelligence, 2017. URL http:
//arxiv. www.iv.orgg.......... abs/ /10es/166...00 609. 09869.
Itamar Daniel Landau and Haim Sompolinsky, Coherent chaos in a recurrent neural network with structured connectivity. PLoS Comput Biol, 14(12)xe1003309, 2018. ISSN 1553-7358, doi: 1o..377/jummaa..cb..00o URL https:t/p///c /xx www.ppl...oooooooo www.ppos..og//ee... http://...111//uunnroe pcbi. 1006309.
Edward N Lorenz. Deterministic nonperiodic flow. Journal of atmospheric sciences, 2002):130-141, 1963.
Edward N Lorenz. Predictability: A problem partly solved. In Proc. Seminar on predictability, volume 1, 1996.
12

---

Niru Maheswaranathan, Alex H. Williams, Matthew D. Golub, Surya Ganguli, and David Sussillo. Reverse engineering recurrent networks for sentiment classification reveals line attractor dynamics. In Advances in neural information processing systems 32, 2019a. URL https: wwwww/wwww......o////. ncbi. nlm, nih, .ov.goov..ooo...oooooo.o..o.oo.l pmc/ a rt cetciiclesppccc.
6638/.
Niru Maheswaranathan, Alex H. Williams, Matthew D. Golub, Surya Ganguli, and David Sussillo. Universality and individuality in neural dynamics across large populations of recurrent networks. In Advances in Neural Information Processing Systems 32, 2019b. URL https:/ //wp::///wwww// ncbi.
nlm. nih. gov pmc/ aew.eetttttteeeeeeeereeeeeeeeett /MClees//MCCC (0MMCC44466399//
Bartlett W. Mel. Synaptic integration in an excitable dendritic tree. Journal of Neurophysiology,
1061 1993. ISSN 0022-3077, 1522-1598, doi: 10.1152/jn. 1993.70.3. 1086. URL
https: wwt.::wwwwp............ http://ppssooooooo org/ http:///.ooooy...ro/////........
(9922/3....9............... 3..3. 1086.
Bartlett W. Mel. Information Processing in Dendritic Trees. Neural Computation, 666):1031-1085, doi: 10.1 (tt.:/cccco1994.......0. URL bttts:://bi.orgg:
1994. ISSN 0899-76677
1162/ 1994. 6.6.033.
Bartlett W. Mel. Why Have Dendrites? A Computational Perspective. In Den-
drites. Oxford University Press, 1999. ISBN 978-0-19-17242--6. URL https:
//oxford. uni
Nima Mohajerin and Steven L. Waslander. Multi-step prediction of dynamic systems with recur- rent neural networks. arXiv: 1880600052 [cs), 2018. UtL. 00526.
Zahra Monfared and Daniel Durstewitz. Existence of n-cycles and border-collision bifurcations in piecewise-linear continuous maps with applications to recurrent neural networks. Nonlinear Dyn, 101(2):1037-152,, 2020a. ISSN 1573-269X. doi: 10.1007/110071-00000555111 URL
https: /ww.//oi.//////////// orgg/00 1007 /000//10000000000811-11--1
Zahra Monfared and Daniel Durstewitz. Transformation of ReL.Ubaaed recurrent neural networks from discrete-time to continuous-time, In Proceedings of the 37th International Conference on Ma- chine Learning, 2020b. URL sestp://proceedinngg mml.eeediiinnmm.........mmm. mlr. press/ /veess/vvvi vlvllll9999 9/wmnf/oooneeeeeeee html.
Alexander Norcliffe, Cristian Bodnar, Ben Day, Jacob Moss, and Pietro Liò. Neural ODE Processes. In Proceedings of the 9th International Conference on Learning Representations, 2021. URL
https: /tttps://pperreee //tp///eeeeeeeeeeee///eeeeeeeeee net wttp:rrrrweeeetttfooooooooo/o2/- f
Chethan Pandarinath, Daniel J. O' O'Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D. Stavisky, Jonathan C. Kao, Eric M. Trautmann, Matthew T. Kaufman, Stephen I. Ryu, Leigh R. Hochberg, Jaimie M. Henderson, Krishna V. Shenoy, L. F. Abbott, and David Sussillo. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(10):805-81 2018. ISSN 1 1548-7105 doi: 10. 101103/445992-00880 -11-01999.. URL https: http://ww...a.r..eeeo
Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In Proceedings of the 30th International Conference on Machine Learning. 2013, URL
.ml.......... www.mm.rre.ssm...mmmmmmmmmmmmmme press/ wt/pv/////////eee ppascanu13 3 -html.
Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott. Model-Free Prediction of Large Spatiotemporally Chaotic Systems from Data: A Reservoir Computing Approach. Phys. Rev. Lett. 12042) 2018. ISSN -901-9007, 1079-7114. doi: Ro.. 120.00004....1.0210.000. URL hwtph:ttt/tstttt aps. http://ooo100/// 1103 www.speessseeeeeeeeeeeee
Mahashweta Patra. Multiple Attractor Bifurcation in Piree-Dimensions Piecewise Linear Maps. Int. J. Bifurcation Chaos, 28(10):183003,, 2018. ISSN 0218-1274. doi: 10.
1142/802181274413300222
URL /
com/ doi /abs//
10. 10.11422/502111224188888
13

---

Barak Pearlmutter. Dynamic recurrent neural networks, 1990. URL https://kilthus _eerrreet__nura_____
cmu, anu.edu/articia journal hawaaaaaoooorbborooooooo/aaaaa ott.a///oo/o//aaaao/aa/aaaaaaaor networks/6605018/1
Panayiota Poirazi and Athanasia Papoutsi. Illuminating dendritic function with compu-
tational models. Nat Rev Neurosci, 21(6):303-321, 2020. ISSN 1471-003X, 1471-
0048. doi: 10. 11.1008//11183-020000000 641583-020-0301-
URL http: /wwww//www//////////////////.. www.naturee... www.uree..//riiiieiieeee
Panayiota Poirazi, Terrence Brannon, and Bartlett Mel. Pyramidal neuron as two-layer neural network. Neuron, 37(6))989-999) 2003.
Cunle Qian, Xuyun Sun, Shaomin Zhang, Dong Xing, Hongbao Li, Xiaoxiang Zheng, Gang Pan, and Yiwen Wang. Nonlinear modeling of neural interaction for spike prediction using the staged point-process model. (0(12):3189-3222 ISSN 0899-7667. doi: 10.1162/ncco.a01133 URL
https: //ooi. wwg///ggoo..ooo..ooooooooo waw.....11nnnoo... 01137.
Maziar Raissi. Deep Hidden Physics Models: Deep Learning of Nonlinear Partial Differential Equations. Journal of Machine Learning Research, 19(2))11224, 2018. URL http:// http:://mmmmmm
btggpaapes//v1//1 046 html.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Multistep neural networks for data- driven discovery of nonlinear dynamical systems. vrXiv:1801.0122 [nlin, physicx:physics, stat], 2018. URL http::////// /ttp//axvvvv.///// org/abs 110011111000 01236.
Danilo Rezende and Shakir Mohamed. Variational Inference with Normalizing Flows. Pro- ceedings of the 32nd International Conference on Machine Learning, 2015. URL http:
//proceedings..... ml www. r.e.ss.ppes/...... weteess/v3//reeeeeee
html.
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic Backpropagation and Approximate Inference in Deep Generative Models. In Proceedings of the 31st International Conference on Machine Learning. 2014. URL http: /ttp:///axvvvv/ee///e/eeeeeaevea org/ abs /1411 4082.
Samuel Rudy, Steven L. Brunton, Joshua L. Proctor, and Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances, 1602664,, 2017. ISSN 2375-2548. doi: 10.1 126/sciadv. 1602144. URL https: ptrttss:/advvaneeee /dwaaccssaaeaeeeeeeea sciencemag ozg/conteett///4 e 602614.
David E. Rumelhart, Geoffrey E. Hinton, and R. J. Williams. Learning Internal Representations by Error Propagation, volume 1, pp. 318-362. Bradford Books, Cambridge MA, 1986.
Tim Sauer, James A Yorke, and Martin Casdagli. Embedology. Journal of statistical Physics, 65(3): 579-616, 1991.
Andrew M. Saxe, James L. McClelland, and Surya Ganguli. Exact solutions to the nonlinear dynamics of learning deep linear neural networks. In Proceedings of the 2nd International Conference on Learning Representations, 2014. URL http: aawaaiiiivvvo....aa.. org/ alss /1112 61200.
Gerwin Schalk, Dennis McFarland, Thilo Hinterberger, Niels Birbaumer, and Jonathan R. Wolpaw. BC12000: puurpoe brain-computer interface (BCI) system. 51(6): 1034-1133, 2000. ISSN 0018-9294. doi: 10.1 109/TBME.2004.82777
Jackie Schiller, Guy Major, Helmut J. Koester, and Yitzhak Schiller. NMDA spikes basal dendrites of cortical pyramidal neurons. Nature, (04(6775) 289, 2000. ISSN 1476-4687. doi: 10.1038/ 35005094. URL https: www. nature. com carmmaaatiii arttiiiiiii... httm///artteeeeeeeeeee./
Dominik Schmidt, Georgia Koppe, Zahra Monfared, Max Beutelspacher, and Daniel Durstewitz. Identifying nonlinear dynamical systems with multiple time scales and long-ram rangrrage dependencies. In Proceedings of the 9th International Conference on Learning Representations, 2021. URL httt:///rxiie http::///ort.ooor//b//////// /11111111110.11..11111
14

---

Thomas Schreiber and Holger Kantz. Observing and Predicting Chaotic Signals: Is 2% Noise Too Much?, pp. 43-65. Springer Berlin Heidelberg. Berlin, Hcidelberg, 1996. ISBN 978-3-
642-80254-6. doi: 10. 10.1007/998-3-644-88222 978-3 001---366288054----
3. URL
1007/
A. Shalova and 1. Oseledets. Deep Representation Learning for Dynamical Systems Modeling. arXiv preprint arXiv:2002.05111, 2020a.
A. Shalova and 1. Oseledets. Tensorized Transformer for Dynamical Systems Modeling. arXiv preprint arXiv:2006. 00444,,00.03345, 2020b.
Martin Stemmler and Christof Koch. How voltage-dependent conductancees can adapt to maximize the information encoded by neuronal firing rate. Nature Neuroscience, 2(6):511-521 1999. ISSN 1546- 1726. doi: 10. 1038/917.. URL https::///// /www. nature com at wrw.iciiilss....... new.eee//nnnn999/
M. Storace and De Feo. Recewise-lineal approximation of nonlinear dynamical systems. IEEE Transactions Circuits and Systems I: Regular Papers, 51(4):830-84 842, April 2004. ISSN 1558- 0806. doi: 10.1 109/TCS1.2004.88 823444 Conference Name: IEEE Transactions on Circuits and Systems I: Regular Papers.
Robert Strauss.
Augmenting neural differential equations to model unknown dynamical systems with incomplete state information. arXiv:2008.08221 08226 [physics, q-biol, 2020. URL http: waaxii.orrgaae 2008 08226.
Floris Takens. Detecting strange attractors in turbulence. In Dynamical Systems and Turbulence, Warwick 1980, volume 898. pp. 366-381, Springer, 1981, ISBN 118-3-540-11171-99 978-3-540- 38945-3. URL http:
i ink. springer com 10. 1007 /B/////000099222
Sachin S. Talathi and Aniket Vartak. Improving performance of recurrent neural network with relu nonlinearity. In Proceedings of the 4th International Conference on Learning Representations,
2016. URLhtts aariiv.vvvv..././vv./.aaa///aava .rrg abs/ /oss/111... 03771.
Adam P. Trischler and Gabriele M.T. D'Eleuterio. Synthesis of recurrent neural networks for dynamical system simulation. Neural Networks, 80:67-78, 2016. ISSN 08936080. doi: 10.1016/ j 1.neumet.2016.04.000 URL
e elsevier. wwt./o/t////eeeeeeeeeeee rttm//etrieeeeii//////
50893608016300314
Pantelis R. Vlachas, Wonmin Byeon, Zhong Y. Wan, Themistoklis P. Sapsis, and Petros Koumout- sakos. Data-driven forecasting of high-ddmensional chaotic systems with long short-term mem- ory networks. Proc. R. Soc. A., 474(2213):201708444, 2018. ISSN 1364-5021, 1471-2946. doi: 10. 10.11088/ssa..2177.08 10. 1098 rapa. 2017.002200000000 0844.
URL https: htwp://oooooeeyyyulllllliiigggga org/ doi/
Grace Wahba. Spline models for observational data. SIAM, 1990.
Ronald J. Williams and David Zipser. A learning algorithm for continually running fully recurrent neural networks. Neural Computation, (2))277-280) June 1989. ISSN 0899-7667, 1530-888X. doi: 10.1 10.1162//eco.199 1989......99.1......0 URL https: : pentts::/direeee /ttptettttttee//ee miit mit 1/270-280/540
Enoch Yeung, Soumya Kundu, and Nathan Hodas. Learning Deep Neural Network Representations for Koopman Operators of Nonlinear Dynamical Systems. arXiv preprint /xXiv:7708.06850 2017. doi: 10.239 11A.1999//.....1881111
Yuan Yin, Vincent LE Guen, Jérémie Dona, Emmanuel de Bezenac, Ibrahim Ayed, Nicolas Thome, and Patrick Gallinari. Augmenting Physical Models with Deep Networks for Complex Dynamics Forecasting. In Proceedings of the 9th International Conference Learning Representations, 2021. URL https:t////o /pwprrvveeeeeeeeee/e eww. eww net/ / forum? htt.//ooommmmmmoooooooooom d=w.@@vrXTTBTT
Xun Zheng, Manzil Zaheer, Amr Ahmed. Yuan Wang, Eric P. Xing, and Alexander J. Smola. State Space LSTM Models with Particle MCMC Inference. arXiv preprint arXiv: 2017. URL deep.
15

---

Qunxi Zhu, Yao Guo, and Wei Lin. Neural Delay Differential Equations. In Proceedings of the 9th International Conference on Learning Representations, 2021. URL
2
16

---

6 APPENDIX
6.1 FURTHER METHODOLOGIICL, DETAILS
Manifold Attractor Regularization As proposed in Schmidt al. (2021), encourage the discovery of long-term dependencies and slow time scales in the data, a subset of Mreg S M states was regularized by adding the following term to the ELBO for the VI approach:
MaCK
Macx
Maxx
Creg A
(AA
Σ (A
1)2
W
(6)
This regularization pushes the regularized subset states toward a continuous set of marginally stable fixed points that tends form an attracting manifold in the full state space, which supports the learning of systems with widely differing time scales, such as the bursting neuron model (cf. Sec. 4). We found that for all datasets this regularization significantly helped to discover the underlying dynamics. To put on equal grounds with the regularization term, the ELBO was divided by the number of time steps T of given batch. Regularization settings used are summarized in Table $1 along other hyper-parameter settings.
To train deterministic version of the dendPLRNN, employ BPTT with scheduled version of TF (Williams & Zipser, 1989; Pearlmutter, 1990). To do so, we choose an "identity- mapping" for the observation model ã; Ize where I E X with Ix if k S N and zeroes everywhere else. This allows us to regularly replace latent states with observations to "recalibrate" the model and break trajectory divergence in case of chaotic dynamics. Consider a time series (21. T2. ( generated by DS we want to reconstruct. At times IT + E No. where 7 2 1 the forcing interval, we replace the first N latent states by observations & -t. k < N. The remaining latent states, & 2R, k: >> N, remain unaffected by the forcing. This means that we optimize the dendPLRNN such that a subspace of the latent space directly maps to the observed time series variables. The forcing interval т is a hyperparameter, with optimal settings varying depending on the dataset. The settings we chose are summarized in Table S1. With F (IT 1eN,, the dendPLRRNN updates then be written
BPTT-TF
(dendPLRNN f else
2++1
(dendPLRNN(z.)
The loss is calculated prior to the forcing, such that L, ||z, - I for every time step. To improve performance we employ mean-centred dendPLRNN (for details see next paragraph). In the evaluation phase, the trained dendPLRNN is simulated freely without any forcing. As the model deterministic, the initial condition 2 = [x | estimated from the first data point æ with a matrix L E NM)) which jointly learned with the other model parameters.
Ceenn-Cettered dendPLRNN Layer normalization has recently been developed as way of significantly improving RNN training (Ba et al., 2016). Here we adapt the idea of layer normalization to the piecewise-lineas nature of our dendPLRNN, Instead of fully standardizing the latent states at every time step before applying the activation function, we only mean-center them:
2 A (o(M()))))) + ha:
M Σ
where o6-) is given in Eq. 4 and M(2+-1) = 2t-1 M4-1 = 21-1
3 where
M ;=1
1 M vector of ones. Note that this mean-centering linear and can be rewritten as matrix-muliplication
M(2 = 21-1 M4-1
M
-1
M
2 = Mz-1-
(9)
M
-1
M
As Remark points out, all results about the tractability of the dendPLRNN also hold for the mean-centred dendPL.RNN.
17

---

State clipping Since the ReLU function used in the dendPLRNN is non-saturating, states may diverge to infinity. As Theorem guarantees, there is a simple and natural way to construct "clipped" dendPLRNN
2t = Az + W
on max(0, Ft-1 - -hs) (0, 2)= + ho.
(10)
Note that the results of Theorem also hold true when the manifold attractor regularization is applied. This is detailed in Proposition further below.
Approximate posterior and hyperparameter settings To estimate the true unknown pos- terior p(z(æ), we make Gaussian assumption for the approximate posterior go(z|æ) (ps(()) Ex))) where mean and covariance are functions of the observations. Without any simplifying assumptions, the number of parameters in Ex E R XMTxM7 would scale unacceptably with time series length T. We therefore made a mean field assumption and factorized qs (z x=)x) across time. Specifically, a time-dependent mean e4,0 and covariance Σt, were parameterized through stacked convolutional networks which take the observations (z.-v... } as inputs, with given by the largest kernel size. The mean is given by 4-layer CNN with decreasing kernel sizes (41, 31, 21 and 11, respectively), with the last layer the CNN feeding into a fully connected layer. For the diagonal covariance, the observations are mapped directly onto the logarithms of the covariance through a single convolutional layer (with kernel size of 41) feeding into a fully connected layer. The classical motivation behind using CNNs rests the assumption that the data contains translationally invariant patterns, and that this allows the recognition model to embed potentially meaningful temporal context into the latent representation (see e.g. Cui et al. (2016)). We note that while the mean-field approximation computationally highly efficient, makes potentially strongly simplifying assumptions (Blei et al., 2017; Bayer al.) that may limit the ability of the encoder model to approximate the true posterior.
To train the dendPLRNN in the VI framework, Adam (Kingma & Ba, 2015), with batch size of 1000 and learning rate of le - 3 was used as the optimizer. For the training with BPTT, we used the Adam optimizer with an initial learning rate of 1e that was iteratively reduced during training down to le 5. For each epoch we randomly sampled sequences of length T 500 (except for the Lorenz-63 runs, where seg 200 time steps were sufficient) from the total training data pool of each dataset, which are then fed into the reconstruction method in batches of size 16. Initial weights were drawn from a uniform distribution.
To find optimal haper-parameters we performed a grid search within Aveg E (0, 10) (VI), TIF {1, 5, 10, 20,50, 100) (BPTT-TF), M E (5, 10, 15, 20, 25,30, 35, 40, 45, 50, 75, 100]. and B E {0. 1, 2, 5, 10, 20, 35, 50}. Hyper-parameter reported in Table S1 below (note that these may not fully agree with the ranges initially scanned, as given above, since attempted to adjust them further in order to approximately match the number of parameters among models in Table S2).
chosen for the benchmarks in Sec. 4 are
Table S1: Hyperparameter settings for dendPLRNN VI/TF for the four different data sets from Sec. 4.
Dataset Lorenz-63
M
B
Miee/ /M
A
TTF
22
20
1.0/
1.0/
/25
Lorenz-96
42/50
50/30
1.0/-
1.0/-
-/10
Bursting Neuron
26
50/47
0.5/-
1.5/-
-/5
Neural Population Model
12/75
5/40
0.2/
1./-
-/5
EEG
117/128
50/50
0.8/0.5
1.0/1.0
-/10
6.2 PERFORMANCE MEASURES
Geometrical measure Dasp used for evaluating attractor geometries (Fig. 4) measures the match
between the ground truth distribution Prue (ar) and the generated distribution Pgen (2: 2) through the
18

---

discrete binning approximation (Koppe et al., 2019)
K
Datap (Prrue (a), Pgen (æ) z)) Σ
Pirue (ar) log
(k)
Pirue (z)
(11)
Pieen x 22)
where K is the total number of bins, and Pime (x) and Pgen (æ x) are estimated as relative frequen- cies through sampling trajectories from the benchmark DS and the trained reconstruction method, respectively. A range of 2 the data standard deviation on each dimension was partioned into m bins, yielding a total of K = m bins, where N the dimension the ground truth system. Initial transients are removed from sampled trajectories to ensure that the system has reached a limiting set. the bin size is chosen too large, important geometrical details may be lost, while it is chosen too small, noise and (low) sampling artifacts with many empty bins may misguide the approximation above. Here we chose a bin number of m = 30 per dimension as an optimal compromise that distinguished well between successful and poor reconstructions.
(k)
Since the number of bins needed to cover the relevant (populated) region of state space scales exponentially with the number of dimensions, for high-dimensional systems evaluating Dsup as outlined above is not feasible. We therefore resorted to an approximation of the densities based on Gaussian Mixture Models (GMMs), similar a strategy outlined in (Koppe et al., 2019). Specifically,
we approximate the true data distribution by a GMM Perue (z) R 1 Ei= W(x, Σ) with Gaussians hyper-parameter
T
centered on the observed data points {æt} and covariance Σ, which we treat that determines the granularity of the spatial resolution (similar to the bin size in Eq. 11). We can generate a likewise distribution by sampling trajectories from the trained models (or one very long trajectory) and placing Gaussians on the sampled data points, Pgen (a)) R Ei=i 21, E) (in the case of VI, rather than sampling, one could also use the model's distributional assumptions to build this posterior across the observations). For the Kullback-Leibler divergence between two GMMs efficient approximations are at hand (Hershey & Olsen, 2007). Here we employ a Monte Carlo approximation
xz. Σ) Detop (Ptrue (x). Pgen (xlz)) = Σ log 1/TE N(x (1) ;
(12)
V=1 1/I. W/ aZ Σ)
where n Monte- Carlo samples (0) are drawn from the GMM representing Ptrue- In practice, set the covariance Σ = I equal to a scaled identity matrix, with a single hyperparamter o². Scanning the range o² E {0.01, 0.02, 0.05, 0.1, 0.2, 0,,, 2,,,5), we found that values for o² 0.1 1.0 to differentiate best between good and bad reconstructions. We chose o² = 1.0 for numerical stability. For this setting, Daup as derived with the binning method and Dess computed through the GMMs also correlated highly on the low-dimensional benchmark systems (see Figure $1).
Power Spectrum Correlation The power spectrum correlations (PSC) were obtained by first sampling time series of 100,000 time steps, standardizing these, and computing dimension-wer wise Fast Fourier Transforms (using sci py . fft for both the ground truth systems and model-simulated time series. Individual power spectra were then slightly smoothed with a Gaussian kernel, normalized, and the long, high-frequency tails of the spectra, mainly representing noise, were cut off. Smoothing width and cutoff values scale linearily with the length of the time series used to compute the spectrum, and were chosen by visual inspection of the individual spectra. Dimension-wist correlations between smoothed power spectra were then averaged to obtain the reported PSC scores.
Mean Squared Prediction Error mean squared prediction error (PE) was computed across test sets of length T = 1000 by initializing the reconstruction model with the benchmark time series up to some time point f, from where it was then iterated forward by n time steps to yield a prediction at time step + n.. The n-step ahead prediction error (PE) is then defined the MSE between predicted and true observations:
T
PE(n) (T
| 2e
(13)
Note that for a chaotic system initially close trajectories will exponentially diverge, such that PEs for too large prediction steps not meaningful anymore (in chaotic system with noise, for large
19

---

-2
-4
log (tsp)))
-6
0.78
-8
-10
2.5
5.0
7.5
Dstsp
Figure S1: Correlation between the binning approximation (m 30) and the logarithm of the GMM approximation (o² = 1) to Datsp on the Lorenz-63 system for different noise realizations and variances. Similar as reported for the K approximation Koppe et (2019) we observed logarithmic relation between the GMM and binning based measures.
50
= 200
0.000
7.0
0.000
0:004
00089
11
30033
0.055
0.0000
0.000
2000
4000
4003
8009
10500
2000
4030
6000
6000
10600
020
Figure S2: Example power spectrum for different values of the smoothing factor o².
the PE may be high even when estimated from two different runs of the same ground truth model from the same initial condition; see (Koppe et al.. 2019)). How quickly this happens depends on the rate of exponential divergence as quantified through the system's 's maximal Lyapunov exponent (Kantz.& & Schreiber, 2004).
6.3 FURTHER RESULTS
(8)
Dstsp
20-step PE
(b)
PSC
20-step PE
# Parumeters
Parameters
1030
1.2
2500
1500
0.99
34
4000
0.00
1.0
12-
10f
10
0.0
0.06
30
25
10
40
0.93
10
15
number bases
number bases
number bases
number bases
Figure S3: Effect of basis expansion for dendPLRNN trained by VI. (a) Agreement in attractor geometries (left) and 20-step ahead prediction error (right) for the bursting neuron system as function of the number of bases (B) for fixed numbers of total parameters. (b) Agreement in power spectrum correlation (left) and 20-step ahead prediction error (right) for the Lorenz-96 system as function of the number of bases (B) for different numbers of latent states ( (M).
20

---

Epoch 10
Epoch 30
Epach 80
Epoch 150
-WWM0000
1.00
1,00 BA
1000000000000
wwwm
Time steps
Time steps
Time steps
Time steps
Figure S4: We observed that the dendPLRNY frequently underwent bifurcations hetween fixed point or various cyclic solutions until it reaches the chaotic behavior of the ground truth system.
Table S2: Comparison of dendPLRNN (Ours) trained by VI or BPTT+TF, and standard PLRNN (Schmidt et al., 2021), trained by VI or BPTT+TF on four DS benchmarks (top) and three challenging data situations (bottom). Values are mean + SEM.
Dataset/ / Setting
Method
PSC
Dep
20-slep PE
Dynamicul variables Parameters
ConPLLRNN' VI
0.997 + 0.001
0.80
+ 0.25
2.le-3 + 0.2c-3
1032
dendPLRNN TF
0.997
0.002
0.13
1 0.18
9.2r
2.8c
1032
Lorenz
PLRNN VI
0.91 0.004
16.6
1 0.4
1.8e-1
1: 0.1e-1
1032
PLRNN TF
0.994
0.001
1: 0.09
4.3e- 1 0.2e-3
1011
dendPLRNN VI
0.987 = 0.001
0.10
+ 0.01
3.1e-1 = 0.9c-1
42
138-
dendPLRNN TF
0.998
0.0001
0.04
0.01
4.1c-2 + 0.8c-2
4480
96
Lorenz-
PLRNN VI
0.93
0.002
1.68
1 0.03
2.le-3 1 0.2e
4260
PLRNN TF
0.996
0.0003
0.05
2.2c-1 0.2e-1
4700
dendPLRNN VI
0.35 = 0.13
7.5 0.4 0.4
6.1e-1 = 0.1c--1
26
2052
Bursting
dendPLRNN
0.76
2.20-2
2010
Neuron
PLRNN VI
0.54
0.01
17.5 + 0.5
1.17 + 0.14
2021
PLRNN TF
0.72 0.07
3.2 2.0
7.0e-2 2.Te-2
2021
Noural
dendPLRNN VI
0.45 1 0.05
0.56 1 0.05
0.82 1 0.09
12
821
Popula-
dendPLRNN
tion
PLRNN TF
PLRNN VI
0.48 = 0.01
0.38 0.15
11:65 + 1.32
7.3
0.68 = 0.09
19 23
13
832
12102
dendPLRNN
0.967
2.8e-2
2
Low
ReNPPLRNNT TF
0.97 + 0.04
6.9 + 5.3
1.5e- 2 + 0.9c 2
22
1032
data
PLRNN TF
0.96 = 0.04
9.0 + 5.4
1.8e-2 = 0.5e-2
90
1011
dendPI.RN
0.940
12.6
6.5e
1032
Purtially
dendPLRNN TF
0.993
0.003
0.54 + 0.16
6.3e-3 = 0.2c-3
22
1032
observed
0.944
0.2
2.7e-1
1020
PLRNN TF
0.994
0.003
0.56 1 0.34
5.0e-3 1 0.2c 3
30
1011
dendPLRNN
0.973
0.006
+ 0.75
3.5e- + 0e
1032
High
deodPLRNN TF
0.995
0.002
0.4 + 0.13
1.6e-3 = 0.1c-3
1032
PLRNN
0.004
+
6.4e-1 0.1e-1
1020
noise
PLRNN TF
0.994 0.002
0.5 4 0.08
4.3e
0.2c 3
22
1032
EEG Dataset Ellecroenccepaalogaan (EEG) data were taken from study by (Schalk et al., 2000) available at https: htpppppponeee oooet......oooo org/ content ceartt/eenntt /et////egmdiibee///ee/eee db/ 0.. 0.. These are 64-channel EEG data obtained from human subjects during different motor and imagery tasks. We trained the dendPLRNN using BPTT+TF on the "eyes open" baseline time series from subject O, which had total of 9760 time steps. The signal standardized and smoothed with Hann function, using numpy. hann ng and a window length of 15. Results for ground-truth and freely generated EEG signals from several brain regions shown in figure $5.
6.4 DETAILS ON DYNAMICAL SYSTEMS BENCHMARKS
Lorenz-63 system The famous 3d chaotic Lorenz attractor with the butterfly wing shape, originally proposed in (Lorenz, 1963), is defined by
21

---

EEG Delaavdd: Frantel P2 (24)
EEG: Detieeet Occipiaal 2822
generisted
4000
8500
10000
BEG Octaset: Temporal (41)
3000
20002
Figure S5: EEG recordings from frontal, occipital, parietal and temporal lobe vs. freely gener- ated trajectories, sampled from the dendPLRNN, trained with BPTT (M = 128, B 50, T 10, Mag/M 0.5, x 1.0).
dr
a x)) + del(()
dt
dt
dy
x((p 2) y + (1
(14)
df
dt
dæ ry 32 + de
Parameters used for producing ground truth data in the chaotic regime were a = 10, p = 28, and 8/3. Process noise was injected into the system by drawing from a Gaussian term de ~ (O, (0000.0ld 0.01dt x I).
Bursting neuron model The 3d biophysical bursting neuron model was introduced in (Durstewitz, 2009) and is defined by one voltage and two ion channel gating variables (one slow and one fast):
-CnV 9L (V - -EL) - Mo (V) (V
(V- - Ena)
+ gKn (V Fxx) + 9Mhh (V F.x)
(15)
9NNDA .33e -025VV
(V ENMDA)
h 110 (V)
(16)
à Tlo0 (V)
Tn
The limiting values of the ionic gates are given by
{moe, Dloo h2o ((Vaaa..ar.,aa
(17)
We borrowed parameter settings from Schmidt al. (2021) to place the system into the burst burst-firing regime:
C. 6FF, gr. 8mS, Er. -80m V, VNN. =00mS, ENa 60m V, VANa -20mV,
kNa = 15, gK = 10mS, F = 900mV, VAK = -25m kg = 5, T. = ms, Jas = 22mS
VAN 15mV kw 5n 200 ms. GNMDA
10.2mS
22

---

Lorenz-96 system The Lorenz- is a nigh-dimensional, spatially extended weather model, also introduced by Edward Lorenz (Lorenz, 1996):
dæ dl
dt, 1...N,
(ri41 Zi -i---tt2220 (4-2))8i
(18)
dl
with (constant) forcing term F. F 8 is common choice that leads to chaotic behavior. Process noise was added as for the Lorenz-63 system, de ~ (O0, 0.01dt x I). In our simulations we used N 10, but in principle the system allows for arbitrary dimensionality.
Neural population model A larger-scale neural population model was recently introduced in Landau & Sompolinsky (2018) to examine the effect of structured connectivity on top of randomly initialized network matrix. Specifically, an independently Gaussian distributed weight structure was combined with a rank-1 component with coupling strength J1. The dynamics of the single unit currents were defined as
dh = - Jo(hh))
Ev ph)).
(19)
VN
where ((h) tan)(h(t)). We produced 50-dimensional chaotic network model based the code provided in Landau & Sompolinsky (2018) using Ji = 0.09 and seeding the random number generator with 35.
The Lorenz-63 and Lorenz-96 systems simulated using sci scipy- py integrate, while for the bursting neuron and neural population model we used the code provided in Schmidt et al. (2021) and Landau & Sompolinsky (2018), respectively.
6.5 THEORETICAL ANALYSIS
Consider the PLRNN with linear spline basis expansion as defined by Eq. 1, Eq. 4, reproduced here for convenience:
B
2t = A2 W ar max(0, Zt-1 - h&) + ho + Cs + Ct2
(20)
where ~ N(0, Σ), E(Ex, -|ex,<]=00 0 for 1. 1,, E R are scalar weighting factors and hs E M different ReLU "activation thresholds", and all other parameters are as conventional PLRNNs (Koppe et al., 2019).
Defining
(b) R((--1) (2,-1 hs) = max(0, z, - -ha),
(21)
Eq. 20 can be rewritten
(()
21
A + W
12(-11)
#1-1
b=1
B
W
(b)
(-ha) + ho + Cs, + cc.
(22)
11(-1)
b=1
where D (b)
D((t-1) = diag
bb)
(b) 1 are diagonal binary indicator matrices with
1-
2
((6)
dm if 2m,t-1 Am,b and otherwise.
23

---

6.5.1
FIXED POINTS AND n-CYCLES OF SYSTEM EQ. 22
Defining
D() (b)
Di(--1)
=
Ot
(ttt1))
b=1
B
Otb D (b)
12(-1)
(- hb),
(23)
b=
WBB Rit-1)
A W DB
0(t-1))
and considering the autonomous system (i.e., without external inputs or noise terms), Eq. 22 can be rewritten as
Z1 W 1) 2 W (t) + ho-
(24)
Fixed points and cycles of Eq. 20, and their eigenvalue spectra, can now be computed in a way analogous standard PLRNNs. Specifically, solving the equation F(z11) 2"1, fixed points the dendPLRNN are given by
(I- (I- (T-WR WRix.))
W m + ho
(25)
where =11
Zt-1 21-1 and det Wo( Pwe
(1) # 0, i.c. WBB steet)
has no M(e)))
eigenvalue equal to 1.
For 1, n-cycle with periodic points {2nh, F(z*)), (2()"), Fx"1(>"))) of map F can be obtained by solving F" F="="") = nm. Therefore, in order to find the periodic points, we first compute F in the following way:
2: F(z--1) = W 1) 1 + W B (((--1 + ho,
2+11 (2( =((z) W W R( 21-1 (Wi()) W h h W h
+Wiil -Iho,
w 2142 = (z)-1) = F(z++1) = ( (()) 1(t-1) 2 + (W (( W W What-1
+ W 11)+1) Who(r) + Wh + (W 1(t+1) W + W (((+1) + I)ho.
n+1 2z4(n-1) = (z(-1) = II
n j+2
B
a)++n-1 21-1 + Σ
II
W
W B ((++n-) (0(+)-3))
1=2
3
n-j+2
+ Wh B
+
Σ 1=2
II
W
Rit+n-))
+ I1) ho:
+1)o.
(26)
1=2
where
W R(t+n-1) II 2=2
B Wit+n-2
BB W(+n-3))
B (tt-1)
Defining + =: nh the periodic point "h of the n-cycle of F can now be obtained as the fixed point of the n-times iterated map Fe as
n
2 I
W B I--
(
Σ
II
W
WhB Whe
Wh
((tmm-))
1=1
24

---

++
+
II
W
a(tm-)) + +1)
(27)
-=2
=
where nd = Zp = Zt ns if IT Wil (()=--)) is invertible, i.e.
det -- w (0
PTI W 9am-ng
(1) + 0,
which implies W == II. W ()) has no eigenvalue equal to
Remark 1. These results about fixed points and n-cycles also hold for the mean-centred dendPLRNN. This can easily be seen by defining W (it-1) == A + W DB Bit-1) and noting that the elements
of (it-1) are now determined by the mean-centred latent states. That is d m,t-1 =1 if 2m,1-1
M 2j=1 j=1 Pj >
1 Σ
bmmb and O otherwise. The rest of the calculations then proceeds as above.
6.5.2 SUB-REGIONS AND DISCONTINUITY BOUNDARIES CORRESPONDING TO SYSTEM EQ 22
Consider system Eq. 22 without external input and noise terms. Denoting hs h2,b ... MM T in Eq. 22, for = 1, 2. B. we can order the elements h2,1 h22, - hj for every 21,2, M.). Without loss of generality, let
1 h22 ... hg,B, j = 1,2 M.
(28)
Then, for every j. we define the intervals I as follows:
1ú11 = (-00, 11,11l
15, = (15 b= 2,3, ... B,
(29)
1,,41 (hj.B, +o0).
By definition of (+)
((t-1) in Eq. 22, the phase space separated into (B + 1)M sub-regions by MB 1) hyper-surfaces discontinuity boundaries. Every sub-region
be defined by the thresholds h2 as Cartesian product of suitable intervals in Eq. 29 for E f1,2, ... .M). (Note that in Eq. 28 we had
instead of strict inequalities' < " obviously the number of intervals, hence sub-regions, would decrease.) each sub-region the matrices
(( 11 = 1, ... B, have a different configuration. Therefore, in Eq. 24 there are (B 1) diffferent forms for D and so for W 1)) and ( as well. Hence, indexing D 1()- W (i(-1) and ((t-1) as D and h B for (1,2, ... (B + 1)M }. Eq. 22 can be written as (r)
(r) W (r)
21 = W 2 W h + hn-
(30)
To visualize the sub-regions and their borders, let for example M = and B = 2. In this case there 9 sub-regions divided by 12 borders. As illustrated in Fig. 6.5.2, there are different matrices (b) 1) 1,2, and D (r) D / 1,2, ... ,,,,, for each sub-region.
6.5.3 BOUNDED ORBITS ARE COMPATIBLE WITH THE MANIFOLD ATTRACTOR REGULARIZATION
Proposition 2. The results of Theorem 2 are also true when the manifold-attractoo regiaiaaii Eq. 6, is strictly enforced for the dendPL_RNN, Eq. 10.
Proof. Assume A, W. d(z11) (see proof of Theorem in Appx. 6.5.6 for the definition) and ho have the partitioned forms
Treg
TT
O Oreg
TO
A =
W =
Anreg
S
W W
vnrey
25

---

Pote-
(8 1
(11
Par-n
(b. 0
D2j = ( (00 )
D6.
my 4
-22a
0 (012)
1 0.)
0
11
1
0 1
52
11 4
11 1
1
Dia n
1
1
on
Dia
oo + )
-h211
(u )
0)
(u (00))
2 n
(1 8.) b)
22
(0 0
(b 11:0 0
(b 0
P
(00)
0
0000 D6 0
Di - (" 8)
Dig
+
))
0
11 0
D/
0 an)
1 1in
21 h11
21 h12
Figure S6: Example of different sub-regions and related matrices D (b)
R)( b 1, and D D(yo r
1, 2, ... ,9, for M = 2 and B 2. Here, it is assumed that the components of h1 = (2.11,2.11 and h22 (11,2, h² satisfy Eq. 28 with <.
ho heg
Treg (Zt-1)
ho-
H(ze-1)
(31)
hi
preg
Snreg (2t_1)
where Iddrep XMreg Ireg Mreg x MMes, Oldrey x Meeg =: Oreg E RrMryy O,S R(M BMM-Mr (MMeM M Meeg, the sub-matrices A A. MreM +1:M =: Anreg E (RMM-M (M- Mreg)
and W W:Mantt:MMMMro respectively, Furthermore, h0 @rre (z) E R and h (Mr Mre, Mrre+1:MM +iM. (2-1) =: @nreg (2t-1) E Mree,
Wnreg E RM (MMrr))x((MMMMMMMM)
are diagonal and off off-diagonal
ho
nreg
In this case |||| = mmax ( (A) = max (maa((Anne)) and
o
A W 0(zT-1-5))
A Aeg S @ureg (Ze-1) + ALee Wneg Pnreg (2+-1) 1)
= Ahrg S (nrre((1-1)) + Aling Wneg &nreg (2(11)
1(
O
| W ho
A neg
S hy
preg
+ A neg h nreg
A S tho
A W neeg ho
ureg
(32)
Thus, for max (Anrog) <
T-2
T-2 ||
1zr || I || ||>||| ||
+
Σ
A w (i
A hol
j=0
j=u
T-2
< ||||| + (e+ || )) (1|| + || |l)
Σ
Anegll
j=0
26

---

2 + [|holl) (1[))] Wneggl) (
<
(33)
- Anegl
6.5.4 PROOF OF PROPOSITION
Proof. For A = (aig) E RMXM W = (wi)) E RMXM = (E1,r, 62,,, ",E))) 8,
(81, 82,1, SM,)) and C = (ci)) E MxM writing Eq. 22 in scalar form yields
M
M
B
21,t
0fj2j,t11
Wij
Ob j 1z,,-1 hgp 11,0
(b)
Σ j=1
Ctyy
8j, + E1,t
j=1
f=1
M
M
Σ
41j2j-11 + Wij Σ
O j.6- 3,0-112j,- |=,,0-1 ) + 11 +
(b)
Σ
Ctj Sja + GtP
b=1
1=1
M
M
=: Σ f15(2))-1)
+ 01o +
Σ
aj + E1t =: Fi(2-11) 1 1-,,2,-.. ,MM (34)
Using this, we can write Eq. 22 in the vector form
2, (Fi(zi-1); Fe(2)-1),.. FAr(2)-1))
(35)
We show that every F; is continuous and Eq. 22 is continuous PWL map. For this purpose, by Eq. 34, it suffices to prove that every (1z,(1)))1) is continuous. According to the definition of the
intervals Jj, Eq. 29, for any j E [1,2, M) have
3j 1111
(b)
0
4b= 1,2, ..,,,,.,4.. B,
(0)
jt-1 = 1, b= b=1,2 ... ,8
2j,t-1 Ljs +
(b)
1--1 = 0, b= 1, ... B
8 2, 3, ... B,
2j Tj,B41
d Ab= = ,..
(36)
Hence, for 1,j 1.j=1,22 1,2, ... M,, each function (z can be stated as
113 = atj 2j,t-1i
2j E T5
[ = (ar + on wis) Rjæ-1 wiyla,15
2j,-1 E 11
fi5(2)001)= :
113 = (ar + Wij Σ B-1 ab) 21,1-1 - -Wis EB. ag hiM
2j,-1 E I
(B+1) (aij + wij E ab) 2j,1-1 Wij E ab hj.b
f13
2j,011 15,411
(37)
Since for every 1, 1... ...
B,
lim 1 (25,-11) 31--1 4,y56
(b)
lim (2,e 113 (h2,)), 2y--1-hha.
(38)
each function fi, (2 is continuous. Hence, Eq. 22 a continuous PWL map z (but has discontinuities in its Jacobian matrix across the borders). Because of these properties, all the results established for standard PLRNNs in Monfared & Durstewitz (2020a;b); Schmidt et al. (2021) apply the dendPI_RN as well, only that the sub-regions and discontinuity boundaries are different. [1
27

---

6.5.5 PROOF OF PROPOSITION
Proof. Defining &, as identical copies of 21.
1
,,,
21
:
(39)
EM+11
:
BMx1
and likewise
hu
ho,1
ha
ho.2
h1
ho
:
ha
ho
hM
ho
ho.
has+1
ho.M+1
hB/
ho
BMX1
BMx1
hAM
ho,BM
ABMXBM diag( AMXxs, AMXM. AMXM)
on W MXM
0.2 W MxM
0g W MMM
a1 WMXM
02 W M
ag W M
WBMXBM
on WMXM
022 WM
On W WMKM
Č$1,t
ē1,t
Č82,1
6,,1
Cst
Cs
Ex
Cs
EBMA
E
ēM.
(40)
ÊM+1,
Cst BMx1
BMx1
EBM.E
one can rewrite the dendPLRNNfrrom Eq. 20 as
& Aë W max(0, ž, - h) + tho + Cs. + Ec
(41)
Now performing the substitution
Vt
21 t 21 h,
(42)
Eq. 41 can be rewritten as the M x B-dimensional "conventional" PLRNN Eq. 5 with
ho (Ā 1)h- hg-
(43)
28

---

6.5.6 PROOF OF THEOREM 2
Proof. It can easily be shown that for every i (1,2, .. .
aashh 0) if sgn(as) = sgn(h);,)
ab max(maxx00
h4,b) 24, 21,1-1 max(0022,0-11)
00 aphi, else
(44)
By defining
0bb max 2c-1 ho) max(0,22111) 0(z-1)= (ōi(zt-);, (
(45)
and
if sgn(as) = sgn(hs,t)
aaphe if sgn(as) = sgn(hr,,)
4.6
i.b
(46)
aghib else
0 c
can conclude that
Low S (z)(2)1)) S '.
where C,
Jown /up
ΣE b= Cib
Jow/up For ca max (c/om 1,1 |c we have
(2)-1) S c7
and so letting c = maxfcr, .2. cM3 yields
M Σ (ōi(z)-))) s
Σ 2 c
(47)
Since
2, = A A + 8(z-11) + ho,
(48)
for T and 2.... ,T, computing 22, "3 2r recursively leads to
B2 Az + W 0(z1) + ho
23 A 21 + A AW )(z)) W x(22) + [A+ + I ho
T-2
T-2
2T = A T
21 +
Σ j=u
A w 0(z +
Σ
A ho-
(49)
1=0
Therefore, by Eq. 47, for every T
2, we have
T-2
T-2
| < ||A|]"
|| +
W|
Σ
||||
Σ
||
|
(50)
1=
Tmax (A) < 1, then lim
|A||T- and
T T+000
W
||
+ hool
lim
|zT|| < ē
W
|
$ + Σ
|||
|holl =
< 00.
(51)
T
1
A
==
1=
29


Reviews:
[{"review_id": "FS_psdtXdqH", "reviewer": "Reviewer_VV4w", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Post-rebuttal: Thanks for addressing my concerns. The authors have updated the manuscript and added experiments. I will not change my recommendation.\n\n#####\nPros: \n\n+ Capacity to capture nonlinear dynamics\n\n+ Competitive performance\n \n+ Simple form and better interpretability \n\nConcerns: \n\n* As a modification to PLRNN, PLRNN could be a baseline and comparing to it would be confirmative.\n\n* SINDy and LSTM require fully observed data. Do the authors consider to compare with other SSM methods like LFADS?\n\n* To demonstrate the expressive power, apart from chaotic attractors, there are many dynamical systems such as fixed attractors, continuous attractors and etc that are commonly used to describe e.g. neural computational (Wang 2006, Mante 2013). Such experiments could be useful for qualitative demonstration. \n\n* Lack for real-world data example.", "summary_of_the_paper": "The paper proposes to modify the nonlinearity of neuron couplings in PLRNN with a linear spline basis expansion. This modification is claimed to boost PLRNN's capacity of capturing arbitrary nonlinear dynamics in low-dimensions. The paper also introduces two training frameworks, variational inference and BPTT. The proposed method is evaluated on various nonlinear dynamical systems.", "summary_of_the_review": "Overall, I vote for accepting. The paper is well written. The method is evaluated on various data and compared to other methods. Hopefully the authors can address my concern in the rebuttal period. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "D9tUf1qjFmU", "reviewer": "Reviewer_B5hP", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The motivation of the model design in this paper is ambiguous. The proposed model is a simplified hybrid linear and nonlinear dynamic system. The reasoning of coupling this model with dendritic computation is not clear. The claimed interpretability of the proposed model lacks both theoretical and experimental justification. The tractability of the training and inference method seems to stem from the VI and BPTT, thereby making the paper lack technical contribution.\n\nEq. (2) and (4) refer to the same function denoted by \\phi, which incurs confusing definitions. Which function is exactly used in the model? \n\nThe proposed dynamical system is termed with piecewise-linear RNN. In time series analysis, \"piecewise\" often refers to the method that separately models the segments of time series [1]. In this paper, the implication of \"piecewise\" is not well explained.\n\nOne of the training methods applied in this paper is based on VI. Since VI on time series/sequence is heavily studied in recent research works, e.g. [2-5], the authors are expected to detail the design of their VI based training and the comparison to existing methods. \n\nIn the experiment section, the result is unconvincing in that the proposed method failed to exhibit consistent out-performance on the metrics/datasets, i.e. mostly only on half of the datasets, the proposed methods can outperform.  \n\n[1] 2007, SDM. A better alternative to piecewise linear time series segmentation\n\n[2] 2015, NeurIPS, A Recurrent Latent Variable Model for Sequential Data\n\n[3] 2016, NIPS, Sequential Neural Models with Stochastic Layers\n\n[4] 2017, AAAI, Structured Inference Networks for Nonlinear State Space Models\n\n[5] 2021, ICLR, Mind the Gap when Conditioning Amortised Inference in Sequential Latent-Variable Models\n", "summary_of_the_paper": "This paper focuses on the nonlinear dynamical system on time series data. It proposed a piecewise-linear (PL) recurrent neural network (RNN) and attempted to couple it with dendritic computation, a computing paradigm from the neural computation discipline. The training method of the proposed model is through variational inference, and BPTT with teacher forcing.\n", "summary_of_the_review": "The motivation of this paper needs more clarification. The proposed model is a simple hybrid dynamic system and the training method is based on existing VI and BPTT and thus the paper lacks novelty and technical contribution. The experiment is not solid enough to support the claim advantages of the model, i.e. interpretability and approximating ability. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "A8FKMDz3ttQ", "reviewer": "Reviewer_caaB", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I think the paper is well-written and easy to follow. The spline basis extension is straightforward and explained well by comparing it with PLRNN. However, I think the overall novelty is not that strong enough. The spline basis expansion is quite common in many problems, like in mixture models, instead of using one model, people extend to mixture models with a linear combination; in kernel study, instead of using a single kernel, people use a mixture of kernels. In early papers such as [Chan et al ACC 1998], people have studied spline bases to model the nonlinear relationship between input and outputs in RNN. Splines are a common choice of basis functions to capture nonlinear functional relationships in a variety of neural systems [Qian et al Neural Computation 2018, Huang et al Neural Computation 2009, Frank et al Journal of Neuro 2002]. Therefore, although spline basis performs well in the context of PLRNN, I don't think the technical contribution is much significant to the community of ICLR.\n\nSome detailed comments:\n\n1. It's claimed that \"we can essentially reformulate a high-dimensional PLRNN in terms of an equivalent lower-dimensional dendPLRNN\". I think the equivalence is only true when the high-dimensional PLRNN has a unique structure of the connection weight W as shown in eq. 40. Therefore, although a dendPLRNN can be formulated as a high-dimensional PLRNN, the opposite direction is not always true. \n\n2. Following point 1, at the end of page 6, \"it becomes clear that the basis expansion enables to reduce the model’s overall dimensionality without compromising performance.\" I think this is only true if the structured assumption of W is valid. If the authors want to compare PLRNN and dendPLRNN in the same space, e.g., M*B, it's needed to show that how much approximation error there will be. \n\n\n\n\n\n\n\n\n\n\n\n\n\n", "summary_of_the_paper": "The authors proposed an extension of PLRNN with a dendrite-like formulation of the nonlinear activation function. They showed that the proposed dendPLRNN maintains the mathematical tractability and dynamical systems interpretation as PLRNN. They compared with alternatives for DS reconstruction with an overall better performance.\n", "summary_of_the_review": "The spline basis expansion works empirically but the technical novelty is not strong enough. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "1: The contributions are neither significant nor novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: Inspired by dendritic nonlinearity, this paper extends previous work on PLRNN/PWL dynamical system modeling by Durstewitz's group. The extension replaces the ReLU nonlinearity with a linear combination of ReLUs. This preserves the theoretical properties of PLRNN, however, the dimensionality of the latent dynamics remains the same, increasing the expressive power of prior PLRNNs. I (area chair) actually read this paper since not all reviewers provided high-quality reviews and one key reviewer is having a personal emergency. Though I appreciate the premise, detailed numerical evaluations, and the inference approach, the novelty is marginal and I do not buy the theoretical advantage of this class of models as presented (see below). Therefore I cannot recommend this paper to appear at ICLR at this time.

Some additional weaknesses that reviewers did not point out:
1. Dendritic nonlinearity is summarized as a point nonlinearity; It lacks the interesting phenomena of dendrites such as nonlinear summation and calcium spikes with its own internal dynamics.
2. The many analytical properties of PLRNN may sound nice on paper, but very impractical. To search for the fixed points and cycles, the amount of required computation exponentially increases as the number of neurons and cycle length increases. In addition the boundary effects cannot always be ignored. In general detailed analysis can become quite non-trivial quickly, e.g., https://arxiv.org/abs/2109.03198
3. High-dimensional PLRNN that approximates a low-dimensional dynamical system due to model mismatch won't have the same topological stability structures. Theoretical analysis of higher-dimensional DS may be very misleading.
Decision: Reject