TRENDy: Temporal Regression of Effective Nonlinear Dynamics
Authors: Matt Ricci
Guy Pelc
Zoe Piran
Noa Moriel
Mor Nitzan
Venue: ICLR 2025 Poster (2025)
Keywords: dynamical systems; neural ODEs
representation learning
Primary area: learning on time series and dynamical systems
TL;DR: We learn reduced order models of PDEs for robust bifurcation prediction.

Abstract: Spatiotemporal dynamics pervade the natural sciences, from the morphogen dynamics underlying patterning in animal pigmentation to the protein waves controlling cell division. A central challenge lies in understanding how controllable parameters induce qualitative changes in system behavior called bifurcations. This endeavor is particularly difficult in realistic settings where governing partial differential equations (PDEs) are unknown and data is limited and noisy. To address this challenge, we propose TRENDy (Temporal Regression of Effective Nonlinear Dynamics), an equation-free approach to learning low-dimensional, predictive models of spatiotemporal dynamics. TRENDy first maps input data to a low-dimensional space of effective dynamics through a cascade of multiscale filtering operations. Our key insight is the recognition that these effective dynamics can be fit by a neural ordinary differential equation (NODE) having the same parameter space as the input PDE. The preceding filtering operations strongly regularize the phase space of the NODE, making TRENDy significantly more robust to noise compared to existing methods. We train TRENDy to predict the effective dynamics of synthetic and real data representing dynamics from across the physical and life sciences. We then demonstrate how we can automatically locate both Turing and Hopf bifurcations in unseen regions of parameter space. We finally apply our method to the analysis of spatial patterning of the ocellated lizard through development. We found that TRENDy's predicted effective state not only accurately predicts spatial changes over time but also identifies distinct pattern features unique to different anatomical regions, such as the tail, neck, and body--an insight that highlights the potential influence of surface geometry on reaction-diffusion mechanisms and their role in driving spatially varying pattern dynamics.

Published as a conference paper ICLR 2025
TRENDY: TEMPORAL REGRESSION OF EFFECTIVE NONLIN- EAR DYNAMICS
Matthew Ricci, Guy Pelc, Zoe Piran & Noa Moriel School of Computer Science & Engineering The Hebrew University of Jerusalem 9190401 Jerusalem, Israel ma teeeew ri icci @mai 1. hu jijii. aa..c........ ac. il
Mor Nitzan
School of Computer Science & Engineering
Racah Institute of Physics
Faculty of Medicine
The Hebrew University of Jerusalem
9190401 Jerusalem, Israel
mor. nit annmmaii........amm hujiiuiji.. ac
ABSTRACT
Spatiotemporal dynamies pervade the natural sciences, from the morphogen dynamics underlying patterning in animal pigmentation to the protein waves controlling cell division. A central challenge lies understanding how controllable parameters induce qualitative changes in system behavior called bifurcations. This endeavor is particularly difficult in realistic settings where governing partial differential equations (PDEs) are unknown and data is limited and noisy. To address this challenge, we propose TRENDY (Temporal Regression of Effective Nonlinear Dynamics), an equation-free approach to learning low- dimensional, predictive models of spatiotemporal dynamics. TRENDy first maps input data to a low-dimensional space of effective dynamics through a cascade of multiscale filtering operations. Our key insight is the recognition that these effective dynamics can be fit by neural ordinary differential equation (NODE) having the same parameter space as the input PDE. The preceding filtering operations strongly regularize the phase space of the NODE, making TRENDY significantly more robust to noise compared to existing methods. We train TRENDy predict the effective dynamics of synthetic and real data representing dynamics from across the physical and life sciences. We then demonstrate how we can automatically locate both Turing and Hopf bifurcations in unseen regions of parameter space. We finally apply our method to the analysis of spatial patterning of the ocellated lizard through development. We found that TRENDy's s predicted effective state not only accurately predicts spatial changes over time but also identifies distinct pattern features unique to different anatomical regions, such as the tail, neck, and body-an insight that highlights the potential influence of surface geometry and their role in driving spatially varying pattern dynamics.
caation-diifusion mechanisms
INTRODUCTION
Nonlinear partial differential equations (PDEs) govern fundamental processes in nature, from the protein waves underlying bacterial cell division (Meinhardt & de Boer) (001) to the morphogen dynamics giving rise to spatial patterns in animal skin (Kondo et al. 2009). Often, the important aspects these equations lie not in the fine-graines detail of their solutions but rather in their low-dinensional. qualitative behavior, sometimes referred to as their effective dynamics (Wilson 1965] Kupiainen 2015) Vlachas et al.) 2022); Will the cell divide not? Will the animal have stripes spots? And, crucially, how do these effective dynamics change or "bifurcate" as a function of PDE parameters? Understanding the relationship between effective dynamics and system parameters is major scientific challenge, and analytical approaches struggle in the

---

Published as a conference paper ICLR 2025
presence of strong nonlinearities (Guckenheimer Holmes 1983) and scale dependencies (inneyy et al. 1992) commonly found in real systems.
(0)
u(a, y. T:
('(
())
a(0)
a
Ramk
gn(a; 0)
a(t)
0
T
T
u(r,, yr 0;8)
y): 0). with parameters 9 and states u(x, y) taking values on a square domain S2 C R?. Spatial features are measured at each I using recursive, multiscale filtering via scattering, + (see Secs. 3 A.1.3). and stacked into single representation, (t) (t) (?),..... where superscripts represent the order of recursion. t thus maps u to reduced-order trajectories which can be further reduced by ranking and subsampling (red). This "effective" representation is controlled by unknown temporal dynamics, à. (Middle) Those unknown dynamics are modeled as a neural network, In. having learnable weights, æ, and which depend on the known, true parameters, 0. (Right) Simulated trajectories from the estimated effective dynamics, ā(t), are initialized on the true (0) (using the initial scattering coefficients St (0)) and regressed against true effective trajectories, a(t), with pointwise loss, L.
Figure 1: TRENDy. (Left) Observed dynamics solutions to PDEs, 14 -N(u<v;
Recently, data-driven approaches have emerged as promising supplement to analytical methods in the study of PDES (Rudy et al. 2017), and substantial progress has been made the fine-grained solution both forward and inverse PDE problems from data (Yang et al. 2202 Li al. (202T). However, comparatively little attention has been given to learning effective PDE dynamics and less still to the parametric case, For example, applied Koopman analysis seeks to explain nonlinear dynamics in terms of the behavior of a few coherent spatiotemporal modes, but it does so at the cost of increasing, rather than decreasing, the dimensionality of the state space (Folkestad et al. 2020). Neural operator methods have been previously applied to parametric models YYang et 2021 ), but focused on fine-gra -reinea reconstruction, rather than effective modeling. Work by Vlachas et. al. (Vlachas et al. 2022), for their part, modeled effective dynamics of PDEs with multiscale features and recurrent neural networks but the learned effective dynamics had no relation to system parameters. Lacking tool to extract the effective dynamics from PDEs a parametrically interpretable way, the ability to identify and control nonlinear systems is fundamentally limited.
A model of effective dynamics should ideally obey few desiderata. First, following the standard accuracy- complexity trade-off, the model should be able to represent the dynamics of the underlying PDE in compact manner; namely, the phase space of the effective system should be not only finite- but comparatively low- dimensional, while preserving sufficient information about the observed dynamics. Second, a model of effective dynamics should be predictive; it should be able to generalize its knowledge these low-dimensional dynamics to novel parameter and state settings not observed directly from data. Such capability is desirable for realistic experimental settings where direct observations of the spatiotemporal dynamics are limited and assessing the system's 's dynamics over the complete configuration space infeasible. Finally, the phase space of the effective model should be interpretable; the relatively few effective phase variables should represent meaningful spatiotemporal variables whose evolution is predicted by the effective model.
In this spirit, we propose TRENDy (Temporal Regression of Effective Nonlinear Dynamics, Fig. T]. a data-driven approach modeling effective dynamics in parameterized PDEs. TRENDy works by mapping
2

---

Published as a conference paper ICLR 2025
the dimeiseddimeniional state of the observed PDE low-dimensional space via recursive, multiscale measurements. These low-dimensional measurements are then fit with neural ODE sharing a parameter space with the observed data. Importantly, TRENDy's 's latent, effective dynamics allow for direct interpretable predictions about spatiotemporal dynamics, circumventing additional generative components, such (learned) encoder or a decoder. We demonstrate this capability on data sets taken from several synthetic and real test cases.
Our principal contributions are
A formal presentation of TRENDy, in particular its use of the scattering transform (Mallat 2012) for extracting multiscale measurements.
Experiments indicating TRENDy can predict qualitative, topological changes in the fitted system's 's behavior, i.e., bifurcations. We demonstrate this for two models of spatiotemporal dynamics, the Gray Scott model and Brusselator, benchmarking performance in several noise and feature conditions.
Demonstrations that TRENDy can predict and explain the emergence of complex spatial patterns like those hypothesized to appear through morphogenesis in biological tissue.
An analysis of spatiotemporal data taken from nigh-resolution videos of the development of the ocellated lizard. These results indicate a link between body geometry and pattern growth and showcase TRENDy's relevance real, noisy data in important biological test case
2 RELATED WORK
Learning parameterized dynamics From a data scientific perspective, parameterized models are especially useful as they provide natural means of generalization new regions of phase space. Classical methods like Galerkin projection (Holmes et al., (012) or finite difference schemes (Strikwerda) (004) have been used derive reduced-order models or solve for system behavior, although these assume prior knowledge of the governing equations, which is often not available. More recently, data-driven methods like SINDyCP (Sparse
Identification of Nonlinear Dynamical Systems with Control and Parameters (Nicolaou et al. 2023 ) have gained prominence. SINDyCP extends the otherwise non-parametric SINDy (Brunton et al. 2016) framework by learning control parameters which can be used to steer the system towards target functions. Deep learning approaches include (Tenachi et al., 2023) who showed how deep reinforcement learning could be used to build parametric models of astrophysical phenomena in an unsupervised manner. None of these methods, however, sought learn effective or reduced-order models, which limits their applicability to settings with strong noise or partial data.
Effective modeling of PDEs The study of effective dynamics seeks reduced descriptions of complex systems that retain essential features. Classical approaches from physics like coarse- coarse-graining (Ehrenfest 1907) and Wilsonian renormalization methods (Wilson 1965) have long been used in statistical mechanics to model large systems by averaging out microscopic fluctuations. Farly data-driven approaches like the equation- free framework (EFF) (Kevrekidis & Samacy| (2009), the heterogencous multiscale method (HMM) (Weinan &
Engquist 2003) and the Flow Averaging Integrators (FLAVORS) (Tao et al. 2010) approximate the interaction
between observed and effective spatial scales with projective integration. Building on this approach (Vlachas)
ct al. 2022) learned latent variable models that capture key dynamical behaviors, enabling efficient long- term predictions in high-dimensional systems like those resulting in turbulence. Contemporary data-driven approaches to learning effective dynamics, howFever, do not learn parametric models, which limits their applicability to the detection of bifurcations, or qualitative (sometimes catastrophic) changes in system's 's behavior, our main interest.
'The TRENDy codebase is available https//gitthh..commmmiaabbbbblll
3

---

Published as a conference paper ICLR 2025
3 METHODS
Our goal is to learn predictive models of data generated by spatially extended, autonomous PDEs of the form
(1) where N is a nonlinear operator, 0 E O C IR* is a vector of parameters, and w(r, y) takes values on the
Ut N [u(xz,))
square domain & E R with periodic boundaries (Fig. [ left). We assume that Eq. undergoes bifurcation (see Sec. A.I.I) at 0 = 0*, whereby an equilibrium point changes stability and the qualitative, long-term behavior of the system changes. We further assume that we only have access to a finite data set of solutions D (z, y, t; } indexed by different parameters, 0 and initial conditions, u(.,-, 0). We seek to build a low-dimensional model of this system which can predict behavior as function of 0 for 0 such that u o (for more details, see Sec. A.1.2..
We observe that, for any Frechét differentiable operator : U + R", the image a(t) 2 (()((t) i classically differentiable function of t. Hence, à exist and a(t) is the solution to
L==(a;;),
an ordinary differential equation sharing the same parameter space as Eq. I (Fig. I right). We refer to Eq.
the effective dynamics and A the effective phase space. The true effective dynamics of Eq. 2 is the reduced-order model we would like to learn but which we only know as an image of our data set 1)
We approximate the effective dynamics with neural ordinary differential equation (NODE) (Chen et al. 2018
gT 10)
defined on the effective phase space and having learnable weights π (Fig. small rounded rectangle). The NODE is instantiated as a multi-layer perceptron whose weights and biases comprise и. The dependence of the dynamics on and 0 is effected by augmenting the input to the NODE to be the vector [a 0J, where denotes concatenation. In all experiments, the NODE a multilayer perceptron with four layers, with 64 hidden units in each layer, and with zero-rectification nonlinearities. It is always initialized on the true a(0).
Each training iteration proceeds first by computing effective initial conditions ag = *[uo]] sampled from the data set, D. Then, Bq.Bis solved with those initial conditions and their associated parameters, 0. The resulting solutions, a, and their derivatives à(t), are compared to their ground truth counterparts, a(t) = 4 [u] (t), by integrating across time:
La a) = T - J. -
1
|la(t) a(t)ll2 + Bl(a(t) a(2 dt.
(4)
Here, J is burn-in in period and B is a regularizer on the derivative term. This loss is minimized across D to approximate the argmin weights π. The approximate effective dynamics, 9n. are then evaluated on A outside of D. particularly 0 for which a bifurcation is expected.
The choice of & is an important consideration in the TRENDY construction. In order to keep TRENDy both expressive and interpretable, we set d to be the scattering transform (Mallat) 2012). The scattering transform is a cascade of multiscale band-pass filters with interleaved nonlinearities and averaging. is provably stable to smooth deformations of the input signal and tends to capture most of the signal energy in its low-order coefficients. Scattering coefficients are computed by iteratively band-pas, filtering, taking a modulus, and averaging, For family of band-pas filters X with hyperparameters, À E A together with low-pass filter, o, examples of zeroth-order first- and second- coefficients are computed as
8d =uto, s (Iu SA,l ②) x)E
S (1lu 2xx 2x * b) dnddEas
(5)
We use dot notation for ODEs and subscript notation for PDEs.
4

---

Published as conference paper ICLR 2025
Clean
Boundaries
Patches
0.090
D0F
K-
Train Tast
0.085
8.00
0.075
906
0.070
0.2
0.063
031
0.055
110
Figure 2: Bifurcation prediction in the GS model with TRENDy. (Left): the bifurcation landscape. The GS model transitions to patterning near the orange curve in F k space. Under the curve, the system is spatially homogeneous (u 1, yellow inset). Near the curve, it produces wide variety of patterns (stripe inset), After the curve, the system is homogeneous again (u x O, teal inset). Training data was taken from a thin region of width .01 and centered on .054. A rectangle of height .01 was held out as test data centered on the bifurcation value of k' .662. (Righr): Detection performance. Numerical continuation was performed on TRENDY trained on five different measurements plus SINDyCP (SG = spatial gradient, Td= TRENDy with d coefficients, FV=Fourier vector, SCP=SINDYCP, details in main text). Estimated k- values (y axis) plotted for each of the three noise conditions (Clean as circles, Boundaries as squares and Patches as triangles) for all models. TRENDy with 10 scattering coefficients had the best performance and noise degradation, See main text for details.
We use Morlet wavelets as band-pass filters and a Gaussian for lowpass. Hyperparameters A, comprise Morlet orientations and scales, and except where otherwise noted
= 8 orientations and j 6 scales. Recursive filtering means that higher-order coefficients are more numerous. Coefficients across orders are stacked into a long vector (Fig. [ blue. yellow, green). Since many higher-order coefficients tend be near zero, we take the top d most actived coefficients across time as our representation of input signal (Fig. [ red bar). The scattering transform is lossy by design but can nevertheless be used for reconstruction given enough coefficients, circumventing the need for a direct decoder when d is large (see Figs. A.5
4 RESULTS
Gray Scott The Gray Scott model (GS) is simple PDE giving rise rather complex spatial patterning dynamics (Pearson 1993). Reaction-diffusion systems, like those of Gray Scott, were investigated in early work by Alan Turing (Turing) 1990) potential mechanisms for morphogenesis in biological systems, and they remain the subjects of intense theoretical study today (McGough & Riley 2004; Wang et al. 2016). Nevertheless, it remains challenging to detect craction-diffusee dynamics real-world systems, as measurement noise and the preponderance of nuisance variables work to obscure the patterning mechanism.
that spirit, we sought showcase TRENDy's 's ability to predict the emergence of spatial patterns the Gray Scott model in strong noise conditions. The Gray Scott model describes the evolution of two chemical species via the coupled PDEs:
D,, 2 un + F(1 u) 2 + (F + k)v
1/ V2
(6)
5

---

Published as a conference paper ICLR 2025
which depend on parameters F, k >> 0 collectively denoted 0 E θ. This system has spatially homogeneous fixed point 1, which loses its stability via pitchfork bifurcation at k -FF44- 4 F, producing two new fixed points which are spatially heterogeneous (Fig. 2 left panel, middle inset). Such a bifurcation resulting patterning called Turing bifureation.
Our goal was to learn reduced-order model of these dynamies that accurately predicted this Tur- ing bifurcation. We were particularly interested in the robustness of this prediction in three noise conditions: "boundaries", whereby randomly rotated polygonal boundary having between 4 and 6 sides was used to mask all but the interior of the spatial domain; "patches", whereby between and 6 square patches zeroed out random locations in the dynamics; and "clean", representing the noise- less data (examples in Fig. AA.4). Gray Scott parameters (1000 train, 250 test) were drawn from a narrow strip in F-k space comprising the rectangle [0.045, 0.055] 00,5) C θ (Fig. left, shaded regions). Test data was drawn from a small region of length .01 around the bifurcation.
TRENDY was fit to the top 2, 5 or 10 scatter- ing coefficients on aver- age across the data set (conditions T2, T5, TI0). We also compared our results to other features, namely, "spatial gradi- ent" (SG): the spatial average of (Vu, Vv); "Fourier vector" (FV): the Fourier power in 10 wave numbers spaced from the maximum value and the Nyquist limit; and "SINDyCP" (SCP): the aforementioned para- metric sparse regression method used for fitting parametric PDEs. The bifurcation for all meth- ods detected using numerical continuation (further training details Sec. A.1.5).
Class
1.0
0,070
Class
0.8233
0.8578
0.065
Test F1
0.6
0.5882
0.060
0.4
0.055
0.2
0.050
Theoretical
0.0400 0.0425 0,0450 0,0475 0,0500 0.0525 0,0550 0.0575 0,0600
0.0
50
12
15
110
Figure 3: Predicting effective patterning dynamics. (Left): true patterning land- scape. Samples were generated from the test region shown in Fig. aad classified with 4-way k-means using all scattering coefficients orientations, scales, up to order 2). K-means classes (example in insets) were (I) dense spots, (2) homo- geneity, (3) sparse spots, and (4) stripes. (Right): 4-way support vector machine was trained on the predicted state of TRENDy using five types of measurements (see main text). theoretical maximum Fl score was computed by classifying all features used to generate the original labels (dashed line). FI scores on test data in each condition are plotted.
After training, found that TRENDy performed better than competing benchmarks (est. K =0024) and degraded gracefully in noisy training conditions (Fig. 2 right). Performance was naturally strongest on the noiseless condition and was comparatively weaker in the boundaries and especially the patches condition. SINDyCP also performed well (est. k2 .0651) in the noiseless condition but degraded sharply with noise. TRENDy performance improved with more coefficients.
In order to demonstrate the interpretability of our approach, we performed a subsequent pattern classification experiment. Understanding the types of patterns that emerge with different parameters is an important area of study in biophysics, systems biology and applied PDEs (Ermentrout) 1991]. To that end, we generated a more focused data set within the testing regime of the former experiment and sought decode pattern classes from the effective dynamics learned by TRENDy in this patterning zone. Ground truth pattern labels were
6

---

Published as a conference paper at ICLR 2025
generated by performing 4-way k-means clustering on the full scattering spectrum of each image (which we confirmed was sufficient to reconstruct pattern quality; see Figs. A.5 A.6 Pattern classes in F-k space are shown in Fig. 3
TRENDY was trained in the same measurement settings (SG, T2, T5, T10 and FV and its effective steady state was used to decode the true label using a 4-way support vector machine. We found that TRENDy's T10 configuration achieved the highest F1 score on test data. Notably, SG suffices to detect spatial heterogeneity in general (Fig. right, "SG"", but fares poorly on distinguishing different patterning regimes. This confirms that TRENDy's effective coefficients can actually be used to categorically decode back to the observable state space. What's more, we confirmed that scattering coefficients can be directly interpreted to explain spectral distinctions between spots and stripes patterns (see Fig. A.7pp
The Brusselator The Gray Scott results were encouraging, but TRENDy was only trained and evaluated in a narrow strip in parameter space. It also sought model equilibrium dynamics, leaving open the question whether or not our framework need be dynamical at all. To resolve those issues, we used TRENDy to learn an entire bifurcation manifold for the Brusselator model, which exhibits a transition to nonstationary dynamics. The Brusselator is Frattionddiffssio the Belousov- Zhabotinsky reaction (Prigogine & Lefever 1968). Its dynamics are given by the evolution of the concentration of two chemical species:
equation used to model chemical oscillations like those observed in
Uc = De 2 u + A+ u²v B( +
1/2 = D.V* 0 - u2v + Bu
System behavior is governed by parameters (A, = e O. The Brusselator has a stable equilibrium e = (A, B/A) as long as B < 1 + A², after which Hopf bifurcation occurs, destabilizing the equilibrium and birthing oscillations. We denote the manifold of Hopf bifurcations by 7 = {AA, B) E B = 1 + A²). Under periodic boundary conditions, the oscillations in the Brusselator begin as plane waves which can grow into traveling waves which interfere in complex ways.
We sought to test how well TRENDy could learn this basic bifurcation behavior from data, especially under strong noise conditions. To that end, we generated sample solutions to the Brusselator with both the aforementioned noise conditions and a condition with a holdout factor, we designed the training set to only have solutions whose parameters 8 = (A, B) satisfied d, = ming- Ey | 11 = 0.15, 0.5, 1.0. The test set always had (9) < c 0.15. We used only the zeroth order coefficients, Si and S2, which simply measure the average concentration of the species u and v throughout the spatial domain. Since oscillations begin as spatially homogeneous plane waves, we reasoned that a global average would be a sufficient detector. We left coefficients in a linear scale as there was no need to co-register coefficients across several scattering orders. Full training details are in Sec. A.1.5
113 > e for
We found that TRENDY achieved low forecasting error on test data in the noiseless, low hold-out (E = 0.15) condition which gradually degraded as the noise and hold-out conditions intensified (Table (TableA..)). An example of the trained model exhibiting Hopf bifurcation is depicted in Fig. 4 Just as in the Gray Scott case, performance was worst in the patches condition and was intermediately worse in the boundaries condition.
Favorable results in solution forecasting do not necessarily imply that our framework accurately predicts the underlying bifurcation structure of the Brusselator. For example, small oscillations around the true steady state would result in relatively low regression error without accurately matching the intended equilibrium behavior. To investigate this deeper question, we used numerical continuation to locate the estimated bifurcation manifold, 7. learned by TRENDy in each condition. To that end, for each condition, we ran 100 numerical continuation experiments corresponding to evenly spaced A E [1.5, 3.5). All points, 8, (if any) determined to mark Hopf bifurcation as A was varied were collected and quadratic fit was applied to
SINDYCP does not have latent state, does not lend itself naturally to this classification problem.
7

---

Published as conference paper ICLR 2025
1.5
B 3.0
B 4.0
4
23
2-
1-
500
1000
1500
2000
500
1000
1500
2000
500
1000
1500
2000
[
No noise
emo.1
Boundaries
Patches
z-10
e 1.0
True: 9-1 4
2=1++A
0:15
4-
1.50 3.75
225
2.50
2.75 3.00
350
150 1.75
2.00
2:25
250
2.75
3500
32s 3.50
Figure 4: TRENDy learns reduced-onder models of the Brusselator. (Top row) TRENDy was trained on zeroth-order scattering coefficients, (S: in solid green, S2 in solid orange) measured from the evolving Brusselator model (Eq. 7 The trained model (here depicted with hold-out parameter .15 and without noise) closely predicted the true effective dynamics across bifurcation boundary (depicted: and 3.5 in 100 steps, numerical continuation was performed on the trained TRENDy model. Whenever a Hopf bifurcation was detected at B, that point (A, B) was tabulated and plotted. Then, a quadratic polynomial was fit to all points within one holdout (c) and noise condition. The true Hopf manifold is in orange at B + A Approximation of the manifold is qualitatively correct across conditions, but worsens in the boundaries and patches conditions.
L.5, 1.5, 3.0, (Bottom row). each A between 1.5
estimate ý. Note that, since the minimum e used for holdout was 0.15, TRENDy never observed an actual bifurcation during training.
Mirroring the forecasting results, TRENDY located the true bifurcation manifold most accurately in low-noise, low-holdout conditions and this capability slowly degraded in the more challenging conditions (Fig. 4. The lower panels of Fig. 4 show, for each c, the estimated bifurcation curves 7 along with the parameters to which they were fit. In the no noise condition, despite some spurious detections. TRENDy provides good fit to the true manifold. The fit weakens in the noise conditions which also tend to exacerbate the effect of e. The estimated manifold is almost always concave up and is always monotonically increasing, matching the qualitative behavior of the true manifold.
Spatial patterning in the ocellated lizard As a final proof of concept of our method, we trained TRENDy to predict the spatial patterning dynamics of the ocellated lizard (T. lepidus). This species of lizard has distinctive skin patterns made from a mosaic of black and green scales (Fig, [77s, upper inset taken from Fofonjka & Milinkovitch (202I)). The dynamics of these patterns through development has been the subject numerous computational and biological studies (Milinkovitch et al., 2023). Importantly, recently work by Fofonjka & Milinkovitch (201)); Manukyan et al. (2017) used combination of computational modeling and real data collection using optical high-resolution episcopic microscopy. to build extremely detailed reaction-dis diffiinniffusion models of scale patterning through development. The authors argue that growth and other
8

---

Published as conference paper ICLR 2025
Cranio-caudial
Medio lateraltrl
2.00
1.75
PCA
1.50-
log(S,)
1.25
1.00
0.75
PCA 2
0.50
PCA 1
0
50
100
t 200
150
250
Figure 5: TRENDY learns scale dynamics of the ocellated lizard. (Left, top inset). A cropped version of the adult ocellated lizard's torso in medio-lateral, cranio-cauda coordinates. Four regions are highlighted in colors corresponding to their quadrant label in this coordinate system. (Left, lower). Dynamics of one patch of scales taken from the origin from juvenile to adult state. Features represented spatially averaged scattering coefficients. Solid lines are the true dynamics and dashed are TRENDy estimates. The insets depict the patch at é = 0, 125, 250. (Right, top inset). An SVM was trained on the final, 10-d state of TRENDy with labels given by the patch's s quadrant. This inset depicts prediction accuracy. (Right, lower). Together with the classifier results, the moderate clustering of TRENDy states in PCA space indicates relation between scale dynamics and anatomical coordinates.
the test data, where TRENDy achieved 95%
mechanical processes influence patterning, much like how work by Murray (1981 explained why stripe emergence in animal fur occurs preferentially in tubular structures like limbs and tails.
that spirit, we used TRENDy learn a reduced-order model of the evolution of skin patches single ocellated lizard through development from its juvenile to adult stages, and then investigate the trained model for relation between pattern dynamics and body coordinates.
We first analyzed data generated by Fofonjka & Milinkovitch (2021); Manukyan et al. (2017) ). and acquired a dataset of patch evolutions by randomly sampling locations ocellated lizard. We fit TRENDy to d 10 scattering coefficients and used patch coordinates (upper left corner) as parameters, 0, of the NODE. More details are found in Sec. A.1.5
resouuuootttoon video of the developing
As before, the trained model closely tracked the true coefficients (Fig. 5A dashed vs solid lines, respectively). order to observe how TRENDY encoded the underlying animal geometry, we then labeled each sub-video by the quadrant from which it was sampled: in image coordinates, north-west, north-east, south-east sesth-west. We then investigated whether quadrant labels could be decoded from the final states, a(T), of the estimated dynamics. If this were to be the case, would demonstrate that TRENDy had learned relationship between scale dynamics and the underlying animal geometry. To that end, we did an 80-20 split of the a(T) into training and testing sets and fit an SVM to the training samples.
We found that this SVM had high testing performance (accuracy .95, weighted F1: .95), indicating that TRENDy had learned differentiate scale dynamics based on patch location and geometry. This is especially
9

---

Published as a conference paper ICLR 2025
clear if we color patch coordinates according to their labels (Fig. 5B, upper inset) and observe how the colors closely respect quadrant boundaries. PCA on the final states, a(T), also reveals significant clustering by quadrant label (Fig- 5B, lower panel, colors corresponding to quadrant labels in panel A). These results demonstrate that, not only can TRENDy fit dynamical data, but it can also highlight the potential influence of surface geometry on reaction-diffusion mechanisms and their role in driving spatially varying pattern dynamics. We note that scale dynamics have also been modeled with tools from statistical mechanics (Zakany
et al. 2022 and we include an application of TRENDy to this setting in Sec. A.1.5 Ising Model.
5 DISCUSSION
TRENDy, our approach to learning predictive models of spatially extended PDEs, synthesizes techniques from reduced-order modeling (Vlachas et al. 2022) with methods from parameterized modeling fitting (Nicolaou ct al. 2023). We have demonstrated this approach in benchmarked experiments on both synthetic and real data, and have used predict bifurcations and visualize the relation between dynamies and geometry in spatial patterning- Among its advantages, TRENDy's 's use of hardwired scattering features makes it robust to noise and avoids feature learning through seff-supervision. which can be costly, Importantly, we have also shown how existing methods for fine-graine uine-granedd spatial prediction are not necessarily performant on bifurcation and classification problems, especially in the presence of noise.
We envision several areas for improvement upon the current framework's 's limitations, Notably, though we have demonstrated that scattering coefficients provide an expressive, interpretable and efficient form of effective measurements, they are by no means the only possible choice. Future work could investigate learning these measurements, possibly with the inclusion of an explicit decoder and reconstruction loss. Here, challenges will be twofold: first, the reconstruction of observable dynamics is underdetermined since TRENDy's 's compression is lossy, though one could regularize decoding with symbolic regression as in Nicolaou et al. (2023); ; second, learned features may be less interpretable than scattering coefficients. Furthermore, while the current iteration of TRENDY performs well on real data obeying equilibrium dynamics, its application to real oscillatory data is limited. Indeed, high-frequency, noisy data simply lead to repeated, conflicting learning signals about the equation of motion at a particular state (Oh et al. (024). Additionally, subsequent studies could investigate frameworks which explicitly account for noise, for example (Frishman & Ronceray 2020).
While our investigation of spatial patterning in the ocellated lizard is largely a proof of concept, one of our principal interests the application to real data representing bifurcations in realistic noise settings and sample sizes. For example, various studies in synthetic biology have produced high- reghhreeolution videos of spatial dynamics of the Min protein system (Glocck et al., 2019; Rajasekaran et al. 2024), , an important regulatory mechanism for cell division and morphogenesis. Investigating how this and related systems undergo qualitative changes in behavior biology, and will require similar data-driven frameworks for robustly modeling the relation between system parameters and function which we have begun to study here.
necessary step for understanding fundamental processes in
ACKNOWLEDGMENTS
We thank the Nitzan lab for their thoughtful feedback. This work was funded by the European Union (ERC, DecodeSC, 101040660) (M.N.). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.
10

---

Published as conference paper ICLR 2025
REFERENCES
Mathieu Andreux, Tomas Angles, Georgios Exarchakis, Roberto Leonarduzzi, Gaspar Rochette, Louis Thiry, John Zarka, Stephane Mallat, Joakim Andaen, Eugene Belilovsky, Joan Bruna, Vincent Lostanlen, Muawiz Chaudhary, Matthew J. Hirn, Edouard Oyallon, Sixin Zhang, Carmine Cella, and Michael Eickenberg. Kymatio: Scattering transforms in python. Journal of Machine Learning Research, 21(60):1-6, 2020. URL htm_
http:// org/ wwtp://apers//vv11/
Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. Julia: A fresh approach to numerical computing. SIAM Review, (9661)65-999 98, 2017.
James J. Binney, Nigel J. Dowrick, Andrew Fisher, and Michael E. J. Newman. The Theory of Critical Phenomena: An Introduction to the Renormalization Group. Oxford University Press, Oxford, UK, 1992.
Steven I. Brunton, Joshua I. Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113(15): 3932-3937, 2016.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Pu Ehrenfest. Begriffliche Grundlagen der statistischen Auffassung in der Mechanik. Springer, 1907.
Bard Ermentrout. Stripes or spots? nonlinear effects in bifurcation of reaction-diffusion equations on the square. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 434: 413-417, 1991.
Anamarija Fofonjka and Michel C Milinkovitch. Deaction-diiffuiion in growing 3d domain of skin scales generates a discrete cellular automaton. Nature Communications, 12(1):2433, 2021.
Carl Folkestad, Daniel Pastor, Igor Mezic, Ryan Mohr, Maria Fonoberova, and Joel Burdick. Extended dynamic mode decomposition with learned koopman eigenfunctions for prediction and control, 2020.
Anna Frishman and Pierre Ronceray. Leamning foree fields from stochastic trajectories. Phys. Rev. 10: 021009, Apr 2020. doi: 10. wtw.///hysee..1102222.... 1009. URL htttt:: ops.org//01/100 1103/PhysRevX. 10. 021009
Roy J. Glauber. Time-dependent statistics of the ising model. Journal of Mathematical Physics, 4(2)29443077, 1963.
Philipp Glock, Beatrice Ramm, Tamara Heermann, Simon Kretschmer, Jakob Schweizer, Jonas Mücksch, Gökberk Alagöz, and Petra Schwille. Stationary patterns two-protein reaction-diffusion system. ACS Synthetic Biology, -157 2019. doi: 10. 10.10221aa htel///ssnbbo..bo/00 Publication Date: December 20, 2018.
URL httpps://oo....... org
10. 10.1021/acssynbioo. 8b00415
John Guckenheimer and Philip Holmes. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer- Verlaa, New York, NY, USA, 1983.
Philip Holmes, John I.. Lumley, Gahl Berkooz, and Clarence W. Rowley. Galerkin projection, pp. 106-129. Cambridge Monographs on Mechanics. Cambridge University Press, 2012.
Ioannis G Kevrekidis and Giovanni Samaey. Equation-free multiscale computation: Algorithms and applica- tions. Annual review of physical chemistry, 60:321-344, 2009.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint 1412v:1412.69888 2014.
11

---

Published as a conference paper ICLR 2025
Shigeru Kondo, Motoko Iwashita, and Motoomi Yamaguchi. How animals get their skin patterns: fish pigment pattern as a live turing wave. The International journal of developmental biology, 53:851-6, 02 2009. doi: 10. 1137//jdd..07250025.
Antti Kupiainen. Renormalization group and stochastic pde's, 2015.
Zongyi Li, Nikola Kovachki, Kamyar Arizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations, 2021. URL
httpss//aaxiv.org///
08895
Stéphane Mallat. Group invariant scattering, 2012.
Liana Manukyan, Sophie A Montandon, Anamarija Fofonjka, Stanislav Smirnov, and Michel C Milinkovitch. A living mesoscopic cellular automaton made of skin scales. Nature, 544(7649):133-77,, 2017.
Jeff' McGough and Kyle Riley. Pattern formation in the gray-scott model. Nonlinear analysis: real world applications,
2004.
H. Meinhardt and P. D. de Boer. Pattern formation escherichia coli: A model for the pole-to-pole oscillations of min proteins and the localization of the division site. Proceedings of the National Academy of Sciences the
the United States of America, 98: 14202 14207, 2001. doi: 10. (0.1077/pnaa.251116599
Michel C Milinkovitch, Ebrahim Jahanbakhsh, and Szaboles Zakany. The unreasonable effectiveness of reaction diffusion in vertebrate skin color patterning. Anmual Review of Cell and Developmental Biology, 39(1):145-111 2023.
Robert Munafo. The pearson classes
pearson classes. html Accessed: 2024-07-26.
James Dickson Murray, On pattern formation mechanisms for lepidopteran wing patterns and mammalian coat markings. Philosophical Transactions of the Royal Society of London. B, Biological Sciences, 295 (1078):473-496, 1981.
Zachary G Nicolaou, Guanyu Huo, Yihui Chen, Steven L Brunton, and Nathan Kutz. Data-driven discovery and extrapolation of parameterized fattern-forming dynamics. Physical Review Research, 544):104 2023.
YongKyung Oh, Dongyoung Lim, and Sungil Kim. Stable neural stochastic differential equations in analyzing irregular time series data. The Twelfth International Conference on Learning Representations, 2024. URL /ttpps:
net foee//oommmid==V NNu
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang. Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang. Junjic Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 8024 8035. Curran Associates, Inc., 2019.
John E. Pearson. Complex patterns in a simple system. Science, 261(5118) Jul 1993. doi: 10.1 1226
Ilya Prigogine and René Lefever. Symmetry breaking instabilities in dissipative systems. Journal of Chemical Physics, 48(4): (8(4))6695-1)) 1990-1700,-0 1968.
12

---

Published as conference paper ICLR 2025
Rohith Rajasekaran, Chih-Chia Chang, Elliott W. Z. Weix, Thomas M. Galateo, and Scott M. Coyle. programmable reaction-diffusion system for spatiotemporal cell signaling circuit design. Cell, 187(2): 145-359.e166 16, Jan 2024. doi: 10.1016fjeell..2233...0000 URL /ttpp:///ol..crg/1..0000E cell. 2211.2023..0 12. 007
Samuel H. Rudy, Steven L. Brumton, Joshua L. Proctor, and J. Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances, 164)x166 2017. doi: 10. 1t..11/6sciadv.6000 1602614. URL
https:// www. sci eecieenncrrrgg. doi /aaaa//s.00. 126/ sci adv. 1602614
Ashkan Shekaari and Mahmoud Jafari. Theory and simulation of the ising model. arXiv preprint
arXiv:2105.00841 2021.
John C Strikwerda. Finite difference schemes and partial differential equations. SIAM, 2004.
Molei Tao, Houman Owhadi, and Jerrold E Marsden. Nonintrusive and structure preserving multiscale
integration of stiff odes, sdes, and hamiltonian systems with hidden slow dynamics via flow averaging. Multiscale Modeling & Simulation, 8(4): (909:1226-13220 2010.
Wassim Tenachi, Rodrigo Ibata, and Foivos I Diakogiannis. Deep symbolic regression for physics guided by units constraints: toward the automated discovery of physical laws. The Astrophysical Journal, 959(2):99, 2023.
Alan Mathison Turing. The chemical basis of morphogenesis. Bulletin of mathematical biology, 52:153-197, 1990.
Romain Veltz. BifurcationKit.jl, July 2020. hal-02902346
URL Bttp:://bal.rrcclveeeevevessss
Pantelis R. Vlachas, Georgios Arampatzis, Caroline Uhler, et al. Multiscale simulations of complex systems by learning their effective dynamics. Nature Machine Intelligence, 4:359-366, 2022. doi: 10. 1038/
0222566022-0044--- URLL:tss///iii..r/110000003 342/5///4222255 0tt.://///222255222202000-00 -22
Jinfeng Wang. Junjie Wei, and Junping Shi. Global bifurcation analysis and pattern formation in homogeneous diffusive predator-prey systems. Journal of Differential Equations, 260(4):3495-3525 2016.
Weinan and Bjorn Engquist. The heterognous multiscale methods. Communications in Mathematical Sciences, 1(1):87-132, 2003.
Kenneth G. Wilson. Model hamiltonians for local quantum field theory. Kenneth G. Wilson Papers, #14- 22-4086. Division of Rare and Manuscript Collections, Cornell University Library. Box 1, Folder 7 1965.
Liu Yang, Xuhui Meng, and George Em Karniadakis. B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data. Journal of Computational Physics, 425:109913, January 2021. ISSN 0021-9991. doi: 10. jop. 109913. URL htthhhtp///te http:///xxxddd..ooo.....o /016/j.
109913
Szaboles Zakany, Stanislav Smirnov, and Michel C. Milinkovitch. Lizard skin patterns and the ising model. Phys. Rev. Lett. 128:048 102, Jan 2022. doi: 10. 10.110./PhysRevL.er. 048102 URLhttts 128. 048102)
7/link aps. o7Timk.aps..rrg/////// http:///....1....//o//.ooooo
13

---

Published as a conference paper ICLR 2025
A APPENDIX
A.1.1 PDES, TOPOLOGICAL CONJUGACY AND BIFURCATIONS
We will focus on PDEs of the form
Ut = DAu + fu(()) = Nolul,
(A.1)
for f a smooth (polynomial) operator on states, u(r; 1), taking values on square domain D E TR2 with periodic boundaries. The subscript 0 denotes set of parameters (D, 11) E θ C RR. Assume that the solution to Eq. A.1 with initial condition Mig u(0) admits solution in the Sobolev space (D) for each t.
We are interested in the long-term behavior of solutions, u, to Eq. A. Let Φ(t, un) be the flow associated to the system, which maps an initial condition uo to its state at time 7. Then, the w-limit set of uo is given by
w (wo) {u E (D) 3 [ + co such that @ft uo) + u
(A.2)
In other words, ((u)) is the set of all states visited by the system starting at ug after infinite amount of time. The set of all w-limit sets of Eq. [ for a given parameter @ is denoted P2g. It may comprise fixed points, limit cycles, strange attractors or yet more esoteric structures invariant under d. It is precisely the nature of Og as a function of 0 which we wish to understand and predict from data.
In many cases, changing 0 will have effect on Szg For example, consider two partial differential equations
Ut NNe, [u] U = [u].
(A.3)
representing two different parameter settings, 0, and 0z. The PDEs ut No, [u] and u, = No, [u] are said to be topologically conjugate if there exists homeomorphism h H12 + H²2 such that the following diagram commutes for all u and t:
o Φo, (t, = 6o, (1))(
where de, (t, u) and Φe, (t, u) are the flows for their respective PDEs. In this case, we write Φe, 60.- If two systems are topologically conjugate, then their orbits (as functions of f) can be mapped homeomorphically to one another so that their dynamics are topologically identical. Importantly, it can be shown that, if co, Q02 then Ωo S0g
However, if there exists a 6* such that, for all c > O there exists 0 such that ||0 1*112 < and 2g * Zo- then we say that the system exhibits bifurcation at 0". Informally, the parameter 0" is a transition point between two qualitatively different (topologically non-conjugate) behavioral regimes of the system. These transitions may represent the emergence of periodic behavior, spatial patterning or other dynamical regimes. Predicting the location of these bifurcations from data is a focus of this paper.
A.1.2 REDUCED-ORDER MODELING AND EFFECTIVE DYNAMICS
Reduced order modeling (ROM) aims to reduce the complexity of nigh-dimensional systems while preserving their effective dynamics. Analytical approaches, such as the method of weighted residuals and Galerkin projection, approximate the solution space by projecting the governing equations onto a lower-dimensional subspace defined by set of basis functions (Holmes al., 2012). Data-driven approaches, including proper orthogonal dccomposition (POD) and dynamic mode decomposition (DMD), leverage observed data to construct reduced models that capture the dominant features of the system (Brunton al., 2020). These methods are increasingly integrated with machine learning techniques to enhance model accuracy and generalization (Hesthaven & Ubbiali, 2018) (Vlachas et al. 2022).
14

---

Published as conference paper ICLR 2025
A
B
Figure A.1: Spatial structures in PDEs. All panels represent solutions the Gray Scott caction-ddiffsionn model. (A) A traveling wave of period equal to about 1/6 of the width of the spatial domain D. Fourier coefficients are good detector the Hopf bifurcation leading to these spatial structures. (B) However, if want to localize Turing bifurcation in terms of both its critical parameter, 8* and the location of its new equilibrium, we must be able to resolve spatial structures with features more sophisticated than Fourier modes. For example, stripes like those depicted here be aliased in Fourier space with spots other visually different structures. (C) Fourier features also struggle to model localized structures like pulses and solitons which occur in systems with strong nonlinearities and coupling across spatial scales. Figures publically available at Munafo
Our goal in this paper is to investigate which reduced order models lead generally to good predictive models of bifurcation structure in PDEs. We will conceive of such models operators (or measurements), 4, mapping from the observable state space of the original dynamics to finite- and indeed low-dimensional space of features. That is, wish know, given arbitrary PDE which bifurcates 0", which operators on PDE data lead reliably to an effective model which also bifurcates at 0*? This bifurcation is the "effect" we want the model to capture.
A few examples are instructive. For instance, imagine that = Nuu; 0 undergoes Hopf bifurcation at 0* whereby an initially stable spatially uniform fixed point, uo(r) c, loses stability preceding the emergence of uniform (plane) waves. The trivial averaging operator
(uu[(1) = 1/2 / u(r, 1) dr
(A.4)
gives a one-dimensional, spatially global signal a(t) = Φ (t))(t) which is differentiable in since + differentiable in the sense of Fréchet. Hence, is the solution the ODE
a₂ g 0)
(A.5)
for some giving the law of motion of the spatial average as function of A. Clearly, Eq. A.5 has fixed point at dg = T [n]] and undergoes a Hopf bifurcation rarely has a closed form, but, if we had a good model of g. we could reliably predict the onset of plane waves.
same (* the underlying PDE. Note that very
Now, imagine instead that the bifurcation 0* leads the onset of traveling waves with equilibrium position and having characteristic spatial wave number & (see Fig. AA. A). As the wave is centered symmetrically about e, we have D[u] = both before and after the Hopf bifurcation: i.e., simple average does not detect the relevant spatial structure associated to the bifurcation. In this case, more appropriate measurement
would be power associated to given spatial frequency, k <
k-2 + h- For instance, if
a(k,, ky, t) = F(u) (Ka ky:
(A.6)
15

---

Published as conference paper ICLR 2025
and define the spectral density (key kyy via
P(kg, ky: (2 Ry t)p²
(A.7)
then
F (1) = -1 2x
P(k cos ksin th, 1)k deb
(A.8)
is a reasonable detector for periodic spatial structures off frequency k. Indeed it is precisely these Fourier coefficients which have been used historically to solve systems like Eq. so we should expect a₁ g(a, 0) to be a good model of the effective dynamics. Even if we don't know the true k, using a composite measurement involving few different Fourier modes would help detect the relevant wave scale while remaining relatively low-dimensional.
4 [r]
di
Yet, being spatially global signal, spectral density cannot reliably detect changes in local structure, illustrated by the emergence of pulses, solitons or other non- non-periodic structures in bifurcating systems with strong nonlinearities or coupling across spatial scales [cite] (Fig. A.IC). From a pure signal processing perspective, we observe, for example, that the cosine grating
I (x v) = cos(Nig.t + kyy)
(A.9)
and the windowed grating
x² + y²
12 v) A exp
1 11 v)
(A.10)
20
To? Hence, any bifurcation heralded by the emergence of this sort of local structure is liable be missed by (global) signal like spectral density.
have the same spectral density at k V 1-2 + kg when A
Evidently, we need to rely on measurement operator which captures information across multiple spatial scales, preferably in an interpretable manner, so that particular spatial structures can be identified as being associated with bifurcation. One such measurement, investigated in subsequent experiments, is described next.
A.1.3 THE SCATTERING TRANSFORM
The scattering transform of Mallat and colleagues is mathematical framework designed extract stable and informative features from signals or images while preserving critical structures, such as texture and edges, that are typically lost under traditional signal processing methods like the Fourier wavelet transforms. The scattering transform extends the concept of wavelet transform by incorporating nonlinear operations, enabling it to capture higher-order interactions within signal.
Formally, the scattering transform is defined as a cascade of wavelet transforms followed by modulus nonlinearity and an averaging operation. Given a wwo-dimensional input signal u(z, y), set of wavelet filters (ıxx h hea is applied, where A denotes the set of scales and orientations of the wavelets. The wavelet transform of u then given by the convolution,
[u + v) / u(v, w)/((x - w) dodw.
(A.11)
The modulus operator |- is applied to the transformed signal to introduce nonlinearity, which crucial for capturing higher-order features:
16

---

Published as conference paper ICLR 2025
(A) = lu t/xl=
This operation is iterated to produce higher-order coefficients:
U (mm) (A1 (A....., Am) (y(m-1)()) ..,..., Am-1) 12mm
The final scattering coefficients obtained by applying low-pass filter, typically averaging operation, each ((m), yielding the scattering transform:
g(m) u = ==U(m)
where @ is a low-pass filter that ensures the invariance of the scattering coefficients to translations of the input signal.
The scattering transform can be discretely implemented as convolutional network where each layer performs a wavelet transform followed by a modulus operation and subsampling. Given a discretely sampled signal u jr, y/. the wavelet filters 1/A are applied, with indexing both orientation, 0, and scale, j < J. For the case of Gabor filter, this could be
1
x² + y²
13.0 Fr y)
exp (i(x cos 0 + ysin 0)) exp
(A.12)
/22
2
the discrete implementation, after applying the modulus operator | the result is subsampled by factor 21, where s J. This process captures multi-scale, multi-orientation features of the signal, constructing higher-order coefficients through iterative application. The final scattering coefficients are obtained by applying a low- low-pass filter, typically a convolution with scaling function ensuring translation invariance.
A.1.4 NUMERICAL CONTINUATION
Numerical continuation methods aim to trace the set of solutions x(u) of parameterized system of ordinary differential equations (ODEs):
dr = f(,
x E EER,
EER,
dt
the parameter μ varies. The goal is to compute the evolution of solutions x(p)) and detect bifurcations (see Sec. A.I.Ip. A solution Zo at a parameter value satisfies the steady-state condition
f po) 0.
To compute the continuation of this solution as i changes, one solves the augmented system:
F(r, f(x 0.
Starting from a known solution (zo, mo),, numerical continuation extends the solution curve by incrementally varying A typical method for this is the predictor-cone coorecor- approach. In the predictor step, a linear approximation the solution branch is made by using the tangent direction, which can be computed by differentiating F(x, 14) = 0 with respect to f.. In the corrector step, Newton's method (or another iterative solver) is used to refine the approximation that f(x, 11) 0 holds exactly.
A bifurcation occurs when the qualitative structure of the solution set changes. This can often be detected by examining the Jacobian matrix D (xr u1). At bifurcation point, eigenvalue of D o
μ) crosses
17

---

Published as a conference paper ICLR 2025
the imaginary axis, indicating a change in stability or the emergence of new solution branches. Numerical continuation algorithms can detect such points by tracking the eigenvalues of the Jacobian as is varied.
When a bifurcation condition, such as the zero crossing of an eigenvalue, occurs within an interval of parameter values, bisection algorithm can be used to refine the location of the bifurcation point. The bisection method works by iteratively halving the interval and evaluating the bifurcation condition at the midpoint, continuing this process until the bifurcation point is found to within a specified tolerance.
A.1.5 EXPERIMENTAL DETAILS
Simulated reaction-dee -ifftoooniffuiion data was created by solving a forward Euler scheme on a square domain of side length 64 and spatial discretization dx 1.0. Ising data was solved on the same domain but with a Glauber dynamics Glauber (963). All scattering coefficients were computed with L orientations and = 6 scales; i.c. up to the maximum scale of 2 64.
Scattering coefficients were averaged across the data set and time and subsequently ranked to produce the top d for use in fitting. For benchmarking, we chose d = 2,5, 10. We also compared these representations to
Spatial gradient (SG): For a solution with channels u, v, this is given by
1 Vu(z, y) dx dy, [S]]
Vo(z,y) y) dæ dy
(A.13)
[22]
where [2] is the area of the square domain. This is meant to detect spatial homogeneity in very general sense.
Fourier vector (FV): this is a vector of 10 wavenumbers. To compute the power spectrum at 10 evenly spaced wavenumbers between the Nyquist limit and the maximum spatial scale in a two-dimensional
square domain, excluding the DC component, we define Kxyquis
N for N = 64, the number of
2
grid points in one dimension of the square domain, and Kmin drx64
1 Then, ten evenly
64
spaced wavenumbers are kn = kmin + n Ak, == 1....... 10, where Ak- Aispiiakmaa Each dimension of the FV representation is given by power at the wavenumber, kn
P(kn) 2 ff |F(ke, kyy1 ds,
(A.14)
where C, is the circle given by V 1/2 + hg k, and F(key ky) is the discrete Fourier transform of the two-dimensional field. To exclude the DC component, the term corresponding to kg ki = 0 is explicitly omitted. We then converted P(k,,) log base 10.
SINDyCP (SCP): SINDyCP (Nicolaou et al. 2023) is sparse regression method for fitting differen- tial equations to data. We used it in its "weak" formulation, whereby equations are fit patches of data integrated against test functions. We used 500 subdomain patches and sequentially -terssolddd least squares (STLSQ) optimizer with threshold 1 x 10-
In all cases, TRENDy used four-layer MLP with rectified nonlinearities for its NODE module. This NODE was run with dt .01 for total duration of = 1.0. Note that since the number of time steps in TRENDy's s solution and the true solution were not always the same, we co-registered the corresponding time series and omitted intervening time steps during loss computation. For the Brusselator experiment. we set the derivative regularizer to be B 10 and it 0 otherwise. We used an Adam optimizer (Kingma & Ba 2014 with learning rate 10-4.
All numerical continuation experiments were performed using pseudo-arclength continuation with Newton-
Raphson correction having threshold x 10-5. In parallel, a bisection algorithm with 15 steps was used to detect bifurcations using eigenvalue threshold of 1.0.
18

---

Published as conference paper ICLR 2025
No Noise
Boundaries
Patches
TRENDy Train
4.23 10-1
7.61 x 10-1
9.44 x 10
TRENDy Test
4.64 x 10
7.59 x 10-1
1.01 x 100
Table A.1: Mean square error between true and predicted reduced-order dynamics in the reduced domain (TRENDy).
Noise
Holdout
No noise
Boundaries
Patches
(E)
6 .015
4.34 x 10-4
7.94
10
1.05 x
ç 0.5
8.91
10-4
1.11
10
7.85 10
= 1.0
1.01 x 10-3
4.25
10 3
1.05
10-2
Table A.2: TRENDy's Forecasting error across noise and holdout (e) conditions.
Gray Scott. Data were generated with parameters F, l drawn uniformly on [.045..055] .055) [0.0, For F -0., the true bifurcation occurs k- .062. We generated 1000 training samples with [ k- > .01 and 250 testing samples with |k- - k-| < .01.. Samples were generated by initializing (u, v) uniformly at (1,0) plus pointwise gaussian noise of mean 0 and standard deviation ,01. The initial condition was thresholded
0. Solutions were produced from forward Euler scheme with di 1, da 1.5 and which was solved for T 4500 steps. The size of the spatial domain was 64 x 64.
For the bifurcation experiment, TRENDy was fit over 2000 epochs, each of which took approximately 31.48 seconds. We used a burn-in period of 10 time steps. RMSE results are Table. A.I
For the classification experiment, we generated 5000 samples from within the testing regime and then performed a new 80/20 train/test split with class balancing. Ground truth labels were created by 4-way k-means clustering the full scattering spectrum. We confirmed these were visually sensible groupings and that they matched standard analyses of the Gray Scott model (Munafo). TRENDY was then fit over 5000 epochs to each feature type (SG, T2-10, FV). SINDyCP learns explicit equation for the PDE in the observable space, and, as we wished to evaluate reduced feature representations of the data, it was not appropriate benchmark. TRENDy's predicted steady state for each feature condition was used a representation of the data to be then classified by a support vector machine (SVM). We used sklearn's 's model out-of-the-boo with no modifications.
Brusselator We used 2000 training and 500 testing samples in the parameter regions described in the main text and controlled by the holdout factor, Solutions were produced with di = .01 and were run until T 20. RMSE on testing data for different holdout and noise values are found in Table A.2
Lizard patterning Data was acquired by randomly sampling 100 locations in the high-resolution, 300- frame video of the developing lizard. We took the d 10 most activate coefficients for training. Dynamics were fit in log scale in order to account for natural scale differences between zeroth, first and second-order coefficients.
Classification was performed on an 80/20 split of TRENDy's 's final state using an out-of-the-bon sklearn SVM which had an 2 regularizer of C 1.0.
19

---

Published as conference paper ICLR 2025
The Ising model. The Ising model simple, standard model of magnetism from statistical physics (Shekaari & Jafari 2021). It has been applied widely in the physical and life sciences, notably in recent work on spatial patterning living matter (Zakany et al. 2022). Unlike other systems investigated the current study, the Ising model is stochastic and discrete-time. Nevertheless, TRENDY can still be used to fit the effective behavior this model, which we demonstrate here.
A typical setting for the two-dimensional Ising model consists of large grid of spins S = s:/=1 where each & can assume a state in {-1,1}. The energy associated to this configuration is
((S) J Σ 8;8j.
(A.15)
(i.5)
where J/ is scalar and the sum is taken over nearest neighbors. The probability of finding the system in a given configuration is
Bis)
P(S) x
(A.16)
where kn is Boltzmann's 's constant, T is temperature scalar, and the constant of proportionality is given by the so-called partition function.
0,200
==
0.175
j-1
0.150
j-3
0.125
==4
0,100
0.075
0.050
0.025
2
6
8
10
t
Figure A.2: Example fit of Ising scales. Solid lines depict the evolution of six scattering scales. TRENDy fit in dashed lines with corresponding colors. Time arbitrarily scaled to 10 total units. The temperature for this sample T 6.46.
There exists critical temperature, Tori at which the system tends in the long term towards a magnetized state in which spins align on average. Above this temperature, magnetization is zero and correlations between spins decay exponentially temperature, these correlations tend to increase, resulting in the formation of so-called magnetic domains, iee. regions of aligned spins. We should expect these domains to grow larger in scale the closer one gets the critical temperature.
function of their distance in the grid. However, as one approaches the critical
We sought examine if TRENDY could learn the relation between temperature and predict domain size. We generated 1000 simulations of a 64 x 64 Ising model which was simulated for 100 timesteps using MCMC. Temperatures were chosen uniformly randomly [4.5, 7.5). range above the critical temperature. We
20

---

Published as conference paper ICLR 2025
True
Mean Scale
Est.
1.4
1.2
0.8
4.334
4.667
5.000
5.333
5.667
6.000
6.333
6.666
7.000
T
Figure A.3: Scale versus temperature. The average scale present after the evolution of the system was measured across range of temperatures. Mean values in ten bins from 4.5 to 7.5 were measured directly from scattering scales the Ising model and are plotted red with standard deviation the shaded region. The TRENDY estimate is in blue. Upper insets show samples Ising configurations three temperature values from left to right: 4.5, 6.0 and 7.5.
expect domain size to grow near the left end of this interval. As expect domain size to be reflected in first order spatial statistics, we restricted the scattering measurements to first order coefficients. Further, as there is complete rotational symmetry in the Ising model, we further averaged over all orientations. That is, TRENDy's effective state was given by
L
ag (1)
T. Σ
(u(,() * vee.i* x
(A.17)
for L = 8 orientations and j = ..... 6 scales, since 26 64 was the size the grid. Each the 1000 Ising samples thus converted into 100 time-step. 6 dimensional time series representing the evolution of these scales. These time series were quite noisy, so we smoothed them with a moving average having window size 10.
TRENDY was fit to 800 these samples with temperatures drawn uniformly from the aforementioned range (e.g. Fig. A.2 On the test data, measured the average scale of TRENDy's prediction the end of its evolution at time T. This average scale was given by E, ja,y (T) for e E, aj (T). As expected, we find an inverse relation between temperature and average scale (Fig. A.3 red), indicating the formation of magnetic domains nearer the critical temperature. This trend is fit closely by TRENDy (blue). Notably, see variance average seale also increase at lower temperatures. This is to be expected since, near the critical temperature, correlation length diverges and single scale dominates.
21

---

Published as conference paper ICLR 2025
(a) Clean
(b) Boundarics
(c) Patches
Figure A.4: Gray Scott samples with noise.
0
20
20
True
40
20
40
20
40
60
Recon
20
40
60
60
20
40
20
40
60
Figure A.5: Example reconstruction: spots. Reconstruction made from minimizing [[((u) S(uest l| where Vest was initialized white noise. Reconstruction was based on all scattering coefficients: 8 orientations and j = 6 scales up to order 2. We tried fitting these coefficients directly with TRENDy, but bifurcation prediction was substantially worse result. Only the T2, T5, T10 cases were both predictive of bifurcations and patterning class.
A.1.6 COMPUTE DETAILS
TRENDy was trained with an Adam optimizer (Kingma & Ba 2014) and with the Kymatio implementation of scattering (Andreux et al. 2020). Synthetic PDE data and the TRENDy NODE solutions were both acquired
by forward Euler schemes (details for each system appear in the main text). TRENDy was trained in PyTorch (Paszke et al., 2019) and numerical continuation was run on BifurcationKit (Veltz 2020) in Julia (Bezanson et al. 2017).
A.2
SUPPLEMENTARY FIGURES
22

---

Published as conference paper ICLR 2025
KS
True
60
60
Recon
40
40
60
60
Figure A.6: Example reconstruction: stripes. Reconstruction made from minimizing S(()) S((ust) where West was initialized as white noise. Reconstruction made from minimizing [[S(u) S((ues)ll where West was initialized as white noise.
23

---

Published as conference paper ICLR 2025
10
064
0.820
0.055
00
0.000
50
50-
0.803
00210
0.004
D.0225
0.03
0,0200
0,015
5.000
0.01
5.005
-0.02
12
Figure A.7: Scale comparison between two Gray Scott samples. (Top Row): Spots, and stripes images shaped 64 64, sampled from class and respectively (see main text). Scattering transform was performed with == orientations and J scales. Rightmost panel: average activity in each of the six scales in first order coefficients, where the average taken across orientations. Both patterns have same dominant spatial frequency (which is controlled by the diffusion constants, identical for these systems), but the dispersion of frequencies different: spots have more power higher scales; stripes, more power lower. (Bottom row): Second order scattering coefficients, again averaged over orientation. Arrays show activations in coefficients generated by filtering first at scale j1 and then at scale j2. First two panels correspond to the patterns above (spots/stripes). Rightmost panel is difference between these activations, which indicates that the 2, j2 = 5, a distinguishing feature. This makes sense since the spots spectrum overall flatter (so that small scales and large scales have similar activations).
24


Reviews:
[{"review_id": "q1AQxuqOZS", "reviewer": "Reviewer_ahWW", "summary": "The paper introduces TRENDy, a framework designed for learning low-dimensional surrogates of complex dynamical systems where underlying PDEs are unknown, and data is noisy or limited. The main contributions of TRENDy include:\n\n1. Modeling Effective Dynamics: TRENDy maps underlying PDE using multiscale filtering into a reduced space and models the reduced representation with a neural ODE. This NODE captures the system's behavior based on its governing parameters, which enables TRENDy to predict system dynamics in new parameter spaces.\n\n2. Predicting Bifurcation: moreover, for a parameter-dependent system, the framework is able to predict bifurcations, where a sudden qualitative change in its behaviors. TRENDy also shows robustness to noise in bifurcation localization.\n\n3. Application to Real-World Data: authors also used the patterning in the ocellated lizard as an example to illustrate how the framework's latent space captures meaningful biological features.", "questions": "Besides several concerns that mentioned in the weakness part, here are several questions regarding the paper details:\n\nFigure 1: what is $S_i(0)$ here? And why do they have different heights?\n\nEq 1: it is better to use $u(x, y)$ rather than $u(r)$, as you are talking about 2D space now.\n\nUsage of subscripts (Line 138 and other notations): the subscripts sometimes are very misleading. e.g. $u_{\\theta}$ and $u_0$.\n\nLine 139: for $u_0 \\notin D$, do you mean interpolation and/or exploration?\n\nLine 141: needs explanation of what $U$ is.\n\nLine 148: similar as what mentioned previously, why do you assume $\\Phi$ is hardwired and unlearned? Can the multiscale filtering parameters be learnable?\n\nFigure 2: this figure needs more details to explain. For example, you should say the inset squares means PDE solutions (otherwise it is misleading).\n\nLine 207: the approximately equal symbol here is incorrect. And moreover, what is SINDyCP? Formula of it? Does it have any assumptions? Have you cited it?\n\nLines 230-232: may need a figure to illustrate your conditions. For example, what is “patches”?\n\nLine 266: need to show why $S_{1, 2}$ almost equals to $<u>$ and $<v>$.\n\nLine 306: use def eq symbol “:=” here in $d_{\\gamma} (\\theta)$.", "limitations": "", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. Originality: \n\nThe approach that combines scattering transform and neural ODEs to model the effective dynamics is novel, especially given its application to bifurcation prediction, a challenging task where data is limited and governing equations are unknown.\n\n2. Quality and Clarity: \n\nThe paper shows rigorous methodology and fruitful details in various experiments. Explanations on filtering operations, the NODE structure, and training details led the model's design to be crafty and reliable. \n\n3. Significance: \n\nTRENDy addresses a crucial question in modeling systems governed by unknown or complex PDEs, where direct analytic solutions are impractical to get. The framework’s adaptability to new parameter spaces may also have numerous applications in real-time system control and scenario exploration. In a nutshell, the authors have shown that TRENDy has the potential to significantly advance research in fields like synthetic biology, physics, climate change and ecology, where such questions regarding complex dynamical systems are pretty common.", "weaknesses": "1. Multiscale Filtering:\n\nThe use of multiscale filtering (e.g., scattering transforms) is central to TRENDy, while the specific choice and design of the filtering process are not fully explored in the paper. Authors should provide more why they prefer this type of dimension reduction technique rather than others (for example, do ablation studies on other type of techs and show the one you mentioned is the best). Moreover, compared with too  many experimental details in the main text (better go to supporting materials), it is necessary to say more on multiscale filtering details, e.g. effects of the choice on scattering coefficients. Such explanations / experiments are essential to keep novelty of the paper, since they are numerous papers working on PDE + DL topics (and some of them should be acknowledged, e.g. PDE-net by Long et al. [1], PINNs by Karniadakis et al. [2], and other papers focusing on effective dynamics, see [3] and [4]). \n\n2. Reconstructing State Space:\n\nJust like lifting and restriction in the equation-free approach, TRENDy should have the module which maps the latent dimensions back to the full PDE state space. Without such an explicit decoder, the ability to verify the reduced dynamics against full state predictions will be limited. It will also become an obstacle for researchers in other fields to explore the explainability by utilizing your model. It seems adding a mechanism for decoding reduced dynamics back into full spatial states or maybe explaining why this is not feasible in your scope is essential.\n\n3. Miscellaneous:\n\nI suggested the reviewers consider the following issues, and if time allows, do some elaboration.\n\na) Extending the experimental scope (e.g. systems with chaotic attractors, or discrete-time systems).\n\nb) Discussing the model’s performance on large datasets and its computational demands in both training and inference.\n\nc) Implementing interpretability techniques (e.g., parameter sensitivity, feature importance) to provide insights on multiscale filtering.\n\n[1] https://arxiv.org/pdf/1710.09668\n\n[2] https://arxiv.org/pdf/1711.10561 \n\n[3] https://www.nature.com/articles/s41467-024-48024-7 \n\n[4] https://pubs.aip.org/aip/cha/article-abstract/34/6/063128/3298062/Tipping-points-of-evolving-epidemiological?redirectedFrom=fulltext", "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": "", "extra_scores": {}}, {"review_id": "A67W8WIwXr", "reviewer": "Reviewer_MJPy", "summary": "The framework introduced by the authors captures the dynamics of spatio-temporal parameterised PDE systems by relying on established reduced order modelling techniques and training a NODE solver. The author(s) showcase that the framework can capture the bifurcation points in various test and real datasets.", "questions": "- Did the authors trial any alternatives to the scattering transforms? To showcase the touted benefits of the chosen approach.  \n- Is there a reason the baseline is not implemented for The Brusselator and real-world examples?\n  - Have alternative types of bifurcations being considered (rate-induced bifurcations)?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The paper is well presented and clearly structured and introduces a framework which could prove to be impactful in the area of parameterised PDE modelling.\n- The authors have included a strong baseline model (SINDyCP) for the Gray Scott example.", "weaknesses": "- The lack of decoder reduces the interpretability of the approach. This point is highlighted by the authors themselves but should be made earlier in the work and ways in which it will be tackled should be outlined.\n  - All references should be checked, there are various instances of spelling mistakes.\n- The addition of legends to Figure 3 would improve the presentation of the results.", "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": "", "extra_scores": {}}, {"review_id": "ylY1qSt6QG", "reviewer": "Reviewer_Lune", "summary": "The paper introduces TRENDy as a data-driven, low-dimensional predictive modeling approach for spatiotemporal dynamics. TRENDy uses multiscale filtering to map data into effective dynamics, then fits a NODE to this reduced space. This process robustly regularizes the NODE phase space, making TRENDy resilient to noise. The so-called TRENDy is examined on synthetic and real data across physical and life sciences. It's posed as a versatile tool for studying and controlling spatiotemporal dynamics.", "questions": "The main question is comparing with the outstanding class of operator theoretic settings.", "limitations": "", "rating": 6, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "N/A", "weaknesses": "The significance of the work is missing by authors specifically compared to the state of the art approaches. In addition, the complexities of fitting by a neural ordinary differential equation need to be justified.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)