Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries: Insights from the Poincaré–Hopf Theorem
Authors: In Huh
changwook jeong
Muhammad Alam
Venue: ICML 2025 poster (2025)
Keywords: neural ODE
dynamical system
bifurcation
symmetry breaking
out-of-domain
Poincaré index
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: Neural ODEs trained exclusively on pre-bifurcation data can forecast symmetry-breaking post-bifurcation behaviors, which can be interpreted through Poincaré index theory.

Abstract: Out-Of-Domain (OOD) generalization is a significant challenge in learning dynamical systems, especially when they exhibit bifurcation, a sudden topological transition triggered by a model parameter crossing a critical threshold. A prevailing belief is that machine learning models, unless equipped with strong priors, struggle to generalize across bifurcations due to the abrupt changes in data characteristics. Contrary to this belief, we demonstrate that context-dependent Neural Ordinary Differential Equations (NODEs), trained solely on localized, pre-bifurcation, symmetric data and without physics-based priors, can still identify post-bifurcation, symmetry-breaking behaviors, even in a zero-shot manner. We interpret this capability to the model's implicit utilization of topological invariants, particularly the Poincaré index, and offer a formal explanation based on the Poincaré–Hopf theorem. We derive the conditions under which NODEs can recover—or erroneously hallucinate—broken symmetries without explicit training. Building on this insight, we showcase a topological regularizer inspired by the Poincaré–Hopf theorem and validate it empirically on phase transitions of systems described by the Landau–Khalatnikov equation.

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries: Insights from the Poincaré-Hopf Theorem
In Huh Changwook Jeong 34 Muhammad Ashraful Alam
Abstract
structing accurate mathematical models that seamlessly in- tegrate observational data with physical principles. Properly modeled dynamical systems allow for the prediction rare events or even unobserved phenomena, providing scientific intuition and opportunities for groundbreaking discoveries. Recently, learning dynamical systems directly from data has emerged as promising alternative traditional model- ing approaches, offering partial automation of scientific discovery (Brunton et al., 2016; Fotiadis et al., 2023; Huh et al., 2020; Kirchmeyer et al.. 2022: Mouli et al., 2024; Nzoyem et al., 2025; Yin et al., 2021a:b). To serve as vi- able alternatives to physics-based models, these approaches must demonstrate robust forecasting capabilities beyond their training domain, making Out-Or-Domain (OOD) gen- eralization a critical challenge (Göring al., 2024).
Out-Of-Domain (OOD) generalization is sig- nificant challenge in learning dynamical systems, especially when they exhibit bifurcation,, a sud- den topological transition triggered by model parameter crossing a critical threshold. A prevail- ing belief is that machine learning models, unless equipped with strong priors, struggle to general- ize across bifurcations due the abrupt changes in data characteristics. Contrary to this belief, we demonstrate that context-dependent Neural Or- dinary Differential Equations (NODEs), trained solety on localized, pre- pre-bifurcation, symmetric data and without physics-base priors, can still identify post-bija mast-bifurcation, symmetry-breaking be- haviors, even in zero-shol manner. We interpret this capability to the model's 's implicit utilization of topological invariants, particularly the Poincaré index, and offer a formal explanation based on the Poincaré-Hopf theorem. We derive the condi- tions under which NODES can recover-or erro- neously hallucinate- broken symmetries without explicit training. Building on this insight, we showcase topological regularizer inspired by the Poincaré-Hopf theorem and validate it empir- ically on phase transitions of systems described by the Lanndu-Khalltnikor equation.
Formally, a continuous dynamical system is represented by phase space Ordinary Differential Equation (ODE):
x(t) = fx(()):)) x M, μ MER", f MXR" TM,
where x(t) a phase space state time l n- dimensional model parameter vector, f(: -) is vector field of the p-parameterized ODE, and M is phase manifold. Here, the model parameters refer to factors that charac- terize the physical environment of the system, such as the mass a pendulum'. Trajectories of ODEs are given by
x(T; = x(0) + L. f(x(t):
1. Introduction
where x(0) is called the initial condition. The goal of leam- ing dynamics is to use machine learning models like Recur- rent Neural Networks (RNNs) (Brenner et al., 2022; 2025; Vlachas et al., 2018) or Neural ODEs (NODEs) (Chen et al., 2018; Rubanova et al., 2019) to accurately approximate the unknown vector field ((::) based on the trajectory data, Clearly, for given trajectory over T,, there two degrees of freedom: the initial condition x(0) and the model parame- ter Therefore, OOD challenges in the learning dynamics problem can be categorized into two types: generalizations to the unseen x(0) and to the unseen (Mouli et al., 2024).
The laws of physical systems are frequently expressed by differential equations rooted in dynamical systems theory. Across all scientific disciplines, a key objective lies in con-
'Elmore Family School of Electrical and Computer Engi- neering, Purdue University 'CSE Team, Samsung Electronics "Graduate School of Semiconductor Materials and Devices Engi- neering, UNIST Artificial Intelligence Graduate School, UNIST. Correspondence to: In Huh <ihuh purdue.edu>> Changwook Jeong <changwook, jeong @ unist.ac.kr>>, Muhammad Ashraful Alam <alam @ purduc.ao purduc.edu>. edu>.
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
Therefooee we will use the terms environment and parame- interchangeably, though the latter specifically referring to numerical value that characterizes the environment.

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
The challenges of OOD learning dynamics become par- ticularly prominent in the presence of bifurcations, where the system undergoes abrupt qualitative transitions in be- havior as a certain parameter of crosses a critical thresh- old (Arnold et al., 2013; Strogatz, 2018). Among the var- ious types of bifurcations, spontaneous symmetry break- ing-where symmetric states lose stability and give way to asymmetric states-stand out for their profound impli- cations across diverse fields, including the Higgs mecha- nism in particle physics (Bernstein, 1974), phase transi- tions in condensed matter field theory (Aranson & Kramer, 2002), switching operations in nanoelectronic (Alam & Za- gni, 2024) and nanophotonic devices (Hamel et al., 2015), and even the dynamies of generative diffusion models (Raya & Ambrogioni, 2023). The wammetry--reakingg bifurea- tions induce significant topological changes in the dynami- cal system due to the creation and annihilation of equilib- rium states, thereby intertwining OOD challenges related to both initial conditions and parameters. Given these dramatic changes, purely data-driven models are widely believed to face challenges when extrapolating across bifurcation points, largely due to a lack of invariant features (Ye et al., 2021) be- tween pre- and nost-bifurcation trajectories. Hence, existing studies have emphasized integrating Mhyyscs-innfrmed pri- ors (Gareía Pérez et al., 2023; Ghadami & Epureanu, 2018; Kalia et al., 2021) or collecting cross cross-paraan parameeer datasets (Li & Yang, 2024) to improve the identifiability of bifurcations.
x()) (x)()))) 0 I....
This formulation describes a system operating within a sin- gle environment characterized by a fixed physical model parameter ft, which implicitly encoded in a single weight vector 0. To generalize NODES to represent multiple envi- romments, 8 should be expressed dependent on μ such that (t) f (x(;;; (())). In practice, however, the exact values of are often unknown. In such cases, the NODEs must also seek a latent vector & that captures the infor- mation of the specific physical model parameter μ = Pe corresponding to each given environment indicator e:
xt) = f(((t); 0g = O(E.)).
The environment-aware weight 0.- O(Ea) can be imple- mented using the Feature-wise Linear Modulation (FiLM) layer (Perez et al., 2018) or hypernetwork structure (Ha et al., 2017). For a specific example, the Context-informee Dynamics Adaptation (CoDA) model (Kirchmeyer et al., 2022) employs a hypernetwork with low-rank decomposi- tion, such that 0, ((E.) 6, + W Here, 0, is the centered weight that shared across all trajectories while & referred to as the context vector, is specifically optimized for each set of trajectories associated with a distinct environ- ment Typically, dim & < dim O = m and it is often as- sumed dim &. = dim μ = n. The matrix W, which is also shared across all environments, maps the low-dimensional context vector & to the m-dimensional weight space.
However, in this work, we make the surprising discovery that context deeeex--eependet free-form NODEs, trained exclu- sively on highly localized pre-bifurcation data, without any physics-base -baiics--aaddd-e priors, can unexpectedly detect the presence of commelry-breakinn bifurcation behavior of hidden symmetries, in a zero-shot manner. We demonstrate how purely data-driven NODE models can identify the bifurcation phenomena through the topological invariant known as the Poincaré index (Brasse- let et al., 2009; Strogatz, 2018), revealing that the models implicitly learn and leverage it as a key invariant feature of the data. To support this observation, we formalize our explanation using the Poincaré-Hopp Hopf theorem (Hopf, 1927: Milnor& & Weaver, 1997) and derive the conditions under which NODEs can either infer or, intriguingly, hallucinate broken symmetry without explicitly learning them. Based on this insight, we present proof-of-concept study intro- ducing novel Poincaré Hore regularization, demonstrated through its application to the theory (Landau & Khalatnikov, 1954) as an exemplar case.
Then, this context-informed NODE is trained with some standard loss functions such as the Mean Squared Error (MSE) between the predicted and ground truth trajectories:
bifurcations and recover the post-
T/At N -2>=>>>>>>
C(0 D&)= =
ΣΣ {x{(t)) - (0), 1
J=1 1=1
(( (W.)) = $ ['(8(0)0 +WE))
for an environment-speeific dataset D. consists of N differ- ent (T/ Lengghh discretized trajectories. Because there are multiple environments E,, sampled with different param- eter values in the training dataset, the final objective is given by a summation over Ext:
Σc(e, + WE D2) + R.W, (2)
min We,W EEE.
(L.K)
where R(.) is a regularizer defined as R(W Ee) AellEe 112 + do IlW.ll, which induces sparsity on WE. and encourages the influence of 9g to be sufficiently large across different environments (Kirchmeyer et al., 2022). de and X are hyperparameters. Note that both the architecture and loss function of the model do not include any specific physical priors, aside from standard assump- tions such as sparsity. After training, the model can simulate
2. Preliminary
Coneext-dependent free-form NODEs. Free-form NODES refer to vanilla NODEs (Chen et al., 2018: Rubanova et al., 2019), where the vector field is a neural network parameterized by trainable weights 0:
2

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
0.75 141
0.75
0.75
0.75
Canonical momenium (p)
0.50
0.50
0.50
0.50
025
0.25
0.25
0.25
0.00
0.00
-0.25
-0.25
-0.25
-0.25
oretts
-0.50
-0.50
-0.50
Saddla point
-0.50
Saddle point
Carrier of orbitt
Center ortall
Homoskiic orte!
Homoclinic
-0.75
-0.75
-0.75
-0.75
-10
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5 1.0
-1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
Canonical position (q)
Canonical position (9)
Canonical position (q)
Canonical position (q)
Figure 1. Phase portraits of the Hamiltonian system (3) for (a) - 0.5, (b) -0.1, (c) 0.1, and (d) lr 0.5. Background
contours represent the min-max normalized values of the Hamiltonian H(g, p: p1) (3). Red indicates higher relative values.
trajectories in new environment either by adapting & based on small number of new observations (few-shot adapta- tion) or by directly modulating &e without further training (zero-shot exploration, analogous to latent space traversal in generative models (Song al., 2023; Wei al., 2024)).
space. Consequently, the OOD condition for parameters arises when a model is trained on parameters ja < Perit but the support of the test data is ue > Perit
Symmetry-breakinng bifurcation, though purely arising from variations in parameters, inherently intertwines with OOD challenges in both model parameters and initial conditions. This complexity emerges because the broken symmetries introduce new asymmetric stable points with separatrices, complicating OOD generalization. Figure shows a rep- resentative example such a bifurcation, a Hamiltonian system with = p2/2 /2 ua² /2 + ///4, known as the double-well potential. The system's dynamics is given by
OOD in initial conditions. The OOD condition for initial states is defined scenario where the model is trained initial conditions x (0) p² (x))), but the initial con- ditions of the test data fall outside the training support supp(p"") (Mouli et al., 2024). As highlighted in recent publications (Göring et al., 2024), OOD challenges in initial conditions become particularly significant when the dynam- ical system contains separatrices. These separatrices parti- tion the phase space M into distinet regions, where trajecto- ries originating from different domains exhibit qualitatively distinct behaviors, such converging to different limit sets (e.g., in multistable systems) or transitioning between di- vergent motion patterns (c.g., from libration to rotation in pendulums). Consequently, each subdomain presents a chal- lenging OOD scenario relative to the others. Formally, one can decompose such system into K disjoint subdomains M; such that M U M.UMMM where θM is the separatrix such that its Lebesgue measure is zero. Notably, if the entire training dataset is collected from one subdomain supp(r"") C Ma but the test data come from a different subdomain supp(p festt C MMisk. this naturally creates the OOD condition, as the model will never encounter that par- ticular regime during training, regardless of the sample size N and time horizon T' of trajectories in the training dataset.
(3) where (q. represent the canonical coordinates, and μ is the model parameter. This system has critical value of Perit 0: shown in Figure I (a-b), when μ 0, the system exhibits single family of orbits centered at the sta- ble, symmetric center (0, 0) (black dots). However, when I > O, the bifurcation occurs as depicted in Figure (c-d); the center of orbits at (0, 0) becomes an unstable saddle point (white dots), and two new centers of orbits emerge at (", 0) (black dots), forming double-well structure. Note that after this bifurcation, any infinitesimal oscilla- tion around (0, 0) will drive the system toward one of these newly emerged wells, thereby spontaneously breaking the original symmetry; although the system remains symmetric under (g, p) - (-4, -p), this hidden symmetry cannot be observed in practice. Furthermore, the homoclinic orbits (dashed lines) originating from the saddle point (0, 0) delin- eate the separatrix, separating the left well, the right well, and the larger outer orbit structures. It naturally derives the OOD boundary under initial conditions.
x (q. -0,H,-0,,71 (p. -
OOD in model parameters. The separatrix and multi- stable structures of phase spaces naturally define the bound- aries of OOD problems terms of initial conditions. Simi- larly, in the context of model parameters, bifurcations serve as key markers for delineating these boundaries in the pa- rameter space. Formally, a parameter value perit is called bifurcation point of the family of parameterized ODEs * = [(x; if, for every neighborhood U of Perit in the parameter space, there exist parameters μ1, #2 E U such that x = f(x; 14)) and x = f(x: (2) are not topologically (or qualitatively) equivalent in some neighborhood of the corresponding fixed points (or invariant set) in the phase
3. Motivating Empirical Observations
3.1. Bifurcation Identification using NODEs
Training procedure. We trained the context-informed NODE (2) the Hamiltonian system (3) under the pre- bifurcation regime exclusively. Specifically, we randomly
For simplicity, this assumes codimension-1 bifurcation with n = 1. For higher-dimensional cases, refer Definition B.1.
3

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
0.75
0.75
(c)
0.76
Canonical momentum (p)
21
0.50
0.50
0.50
0.00
0.00
0.00
-0.25
-0.25
-0.25
22
-0.50
0.50
Treining dela
Adeptalon data
Vanila NODF
Contlcc-inlomeed NO0E
-0.55
0.75
-0.75
2
-1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
-1.
-0.5
0.0
0.5
1.0
Canonical position (q)
Canonical position (q)
Canonical position (q)
Canonical position (q)
Figure Examples of (a) training (μ = 0.5) and (b) adaptation (μ 0.5, the left well case) trajectories. Predicted phase portraits
of the (c) vanilla NODE and (d) context-informed NODE models, with background contours representing the numerically computed
Hamiltonian profiles. In (c), the Hamiltonian is physically meaningless as the constructed dynamics is not conservative system.
Table 1. MSEs (x10 between the ground truth and models.
0.1
0.75
Founs
Liner tring
0.50
SCENARIO
VANILLA
CONTEXT
Exrapo
0.0
0.25
LEFT WELL
2427 #00.
0,1517
1.0.668
Context
0.00
RIGHT WELL
1724 1365.1
0.2504
0.033
-0.1-
-0.25
OUTER ORBIT
3943 2004.8
0.1457 10.065
-0.50
Ground Rith Conicat infommss shat
-02 (8)
-0.75
sampled four initial conditions from the uniform distri-
-20 -1.5
-1.0
-0.5
0.0
0.5
-1.0
-0.5
0.0
0.5 1.0
Parameter (4)
Canonical position (e)
bution (q(0), p()) ~ 2([--2.0,20033 for each of the eight parameter values 1 (-2.0, -15, 1.25, 0.255- For each sampled initial con- dition, simulated the dynamics with time horizon = 2.0 and a time step Al = 0.1. This results in |D.| = 4 trajectories per value of u., yielding a total of [ x8 8 = 32 training trajectories. It is important to note that all train- ing data consists solely of single-orbit trajectories from the pre-bifuraation regime that exhibits the symmetric single- well structure. Figure (a) shows the phase portraits of the training trajectories corresponding μ = -0.5 (refer to Appendix D for other cases and experimental details).
Figure 3. (a) Relationship between f4g and Ee. (b) The phase por- trait of the context-informed model (jee = 0.5) extracted from (a).
-1.0, 00755 0.5,
yet remaining confined to the left well, does not improve the performance of this model. In contrast, the context-informed NODEs successfully reconstruct the entire phase topology across all scenarios. Remarkably, they successfully recover the hidden symmetry in the post-bifurcation regime, despite being trained exclusively on pre-bifurcation data and fine- tuned using only a single, chosen trajectory resulting from the broken symmetry, Table I shows the MSE between 32 test trajectories and NODE-predictee trajectories.
One-shot adaptation. After training the model with pre- bifurcation data, adapted it using a single trajectory for = 0.5, representing post-bifurcation data. We consid- ered the following three broken symmetry scenarios, namely adaptations using: (i) a trajectory confined to the left well, (ii) one confined to the right well, and (iii) one outside the separatrix, traversing the outer orbit. These scenarios limit the model's exposure to the global structure of the phase space during adaptation, reflecting the symmetry breaking observed in real-world situations. For comparison, we also trained vanilla neural ODE models (I) under each scenario.
Zero-shot exploration. Figure 3 (a) illustrates the rela- tionship between the ground truth parameters /e and the cor- responding context & constructed during the pre-bifurcation training. The plot reveals nearly linear correlation between the two sets of values. Building on this observation, we ex- trapolated the context value for the post-bite bost-bifurcation case where j4c = 0.5, obtaining & = -0.2058, and integrated it into the context-informet model. Note that this process was conducted without utilizing any post-bifurcatio data adaptation, thus we refer to zero-shot exploration. Figure 3 (b) presents the resultant phase portrait discovered in zero-shot manner, showing it reproduces the ground truth one accurately (refer to Appendix D for different values).
Figure 2 (b-d) compares the phase portraits of the ground truth (3), vanilla NODEs, and context-informed NODES for the left well adaptation scenario near (-vp, 0) (refer to Appendix D for other scenarios). The vanilla NODES struggle to accurately capture the double-well structure, replicating the topology only within the region covered by the adaptation data. This result is consistent with recent theoretical findings suggested in (Göring et al., 2024). Note that it is a fundamental limitation of vanilla NODEs in OOD: increasing the number of post-bifund bifurcation samples,
Automated generation of bifurcation diagrams. Lever- aging the zero-shot exploration capability and fully differen- tiable structure of NODEs, we can automatically construct bifurcation diagrams. To achieve this, we employ the neural- adjoint method (Ren et al., 2020) commonly used for inverse problems with neural networks. Specifically, for a given &, we identify fixed points x" such that f(x*) = by solving

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Compspaecing parameterr (e)
Comaaponding no parumeter
-0.5
0.0
0.5
1.0
2,0
-0.5
0.0
0.5
1.0
2,0
Location of fixed paints (q
Comter (Model) O-nter (ST)
Center (Modal)
Conter (ST)
Location fixed points (q)
Center (mode)
Center (OT)
Conter (modet)
Conttr (ST)
Socdie (Mode) -aado (BT)
Sacde (Mode)
Boode (GT)
Soddle (mlode))
Soddle (GT)
Sacdle (mode) Saddle (GT)
0-
0 0-0000
-2
(bb
-2 (8)
-2
(b)
0.1
0.2
0.3
0.4
0.06
0.10
0.15
0.20
0.25
-0.5
0.0
0.5
1.0
1.5
2.0
-0.5
0.0
0.5
1.0
1.5
2.0
Context (-8)
Context (-0
Parameter (v)
Parameter (u)
Figure 4. Bifurcation diagrams generated using context-informed models trained on (a) default and (b) localized domains.
Bifurcation diagrams generated by context-informed models trained under two realistic settings: (a) the noisy setting, where small Gaussian noise (a 0.02) is added to the trajectories to mimic measurement noise. (b) the limited data setting, where only two initial conditions are used per parameter, compared four in the default setting. Each setting is repeated five times with random initializations, indicated by distinct symbols.
Figure 5.
the inverse problem x* argmin xlgminx)lf((x;))) 2 using the neural-adjoint method. Once the fixed points x' are found, we compute the eigenvalues of the Jacobian matrix of ;x( :) at these points. The fixed points are then clas- sified based on the characteristics of their eigenvalues, such as sign and the presence of real or imaginary components. By modulating Ç and repeating this procedure, we construct bifurcation diagram that intuitively illustrates the emer- gence or disappearance of fixed points. Refer to Appendix C for detailed description of the proposed method. Figure 4 (a) presents the bifurcation diagram generated fully auto- matically, illustrating that a yymmetry-beakinn bifurcation occurs upon crossing the critical threshold, Eerit-
(xx,, we find ((x) = (0, 7g), with y = -8.3 which corre- sponds closely to the slope shown in Figure (a).
This analysis confirms that the model can approximate the parameter modulation as a feature map ((-). However, to fully model the symmetry-beakingg bifurcation, the model must accurately capture not only μi but also q². This raises pertinent question arises: Can the model accurately recon- struct the bifurcation of the dynamical system when trained on a highly constrained phase domain, such as the immedi- ate vicinity of (0, 0)? This presents a significant challenge, it genuinely involves simultaneously addressing both parameter OOD and initial condition OOD issues.
Robustness test. To evaluate whether context- context-informeed NODEs can identify the bifurcation structure under realistic conditions, including noisy observations and limited training data, we conducted additional robustness tests, as shown in Figure 5. These challenging settings revealed that the model remains capable of accurately identifying symmetry- breaking bifurcations, albeit with slightly increased variance compared to the ideal noise-free scenario.
3.2. Bifurcation Identification with Localized Domain
Experimental setting. We trained the context-informed NODE (2) on the Hamiltonian system (3) for pre-bifurcation scenarios, focusing on localized training domains. The ex- perimental setup is consistent with the previous configura- tion, except that initial conditions were sampled from nar- rower uniform distribution, (a(0), ~ U([-0.5, 0.5]²) (see Appendix E for details).
Context-informeet NODES identify broken symmetry. The above presented experiments reveal that the context- informed model can infer the symmetry breaking behavior in the on pre-bifurcation data in a purely data-driven manner. It showcases the model's capability to capture how parameter modulation influences the vector field, even without direct exposure to appropriate data. The most plausible expla- nation for this OOD learnability is that the context-aware model accurately captures the normal form of (3). Note that by using the Taylor expansion near & = Eerit,
regime, despite being trained solely
Remarks on the experimental setting. Within this re- stricted domain, the cubic term qa in (3) becomes less sig- nificant compared to the linear term sag, causing the training trajectories to resemble linear oscillators (ÿ. p) = (p. μq). This leads the model to predominantly observe linear be- havior, potentially causing it to underestimate the role of the cubic term. Consequently, this limitation is expected to make the model prone to overlooking the formation of wellle-well structure caused by the influence of the cubic term. Note that, for the linear system (q. P) (p, pq), bi- furcation still occurs at (0, 0) when Merit = 0, transitioning from a center of orbits to a saddle point (refer to Remark B.2 of Appendix B for details). However, unlike (3), this system does not produce additional double centers of orbits, homoclinic loops, and outer orbit structure.
fxx; ((g)) = f(x; ((Gert)) Vef(x; 0(E)) (emi))eeei (5 Ecritt
= (((:((()))) + (xx)- (g - Serit),
where ((5) 0 + WE and ((x)
Vef(x; 0g
W W is the feature map. If the model is properly de- composed such that (x) ~ (p. -²") and (xx) ~ (0, 19) where y is a constant, then the model can accurately repli- cate the bifurcation behavior. Indeed, when we compute
5

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
0.75
0.75
Canonical momenium (p)
1.0
0.50
0.50
0.50
052
0.25
0.25
0.00
-0.25
-0.25
-0.5
-0.5
-0.50
-0.50
-0.50
Contodieniomed, 0-00.51
Conmot-ifomeed 5--0 100
formeal-fforee 6- 134
-1.0
Conteck nlammed {- -0.218
-0.75
-0.75
0.75
-10
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
1.0
Canonical position (q)
Canonical position (q)
Canonical position (q)
Canonical position (4)
Figure 6. Phase portraits of the context-informed model trained with (q(0), ((0) ~ 0.5) x [-0.5, 0.5).. and constructed
zero-shot manner by linearly extrapolating &., correspond to (a) 0.5, (b) F.c 0.1, (c) pt. = 0.5, and (d) μse 1.5.
Context-informer NODEs identify broken symmetry, unexpectedly. We conducted the same zero-shot explo- ration procedure in Section 3.1. Figure 6 (b-d) presents the post-bifurcation phase portraits reconstructed in zero- shot manner by inputting the explored context values & for He = 0.1, 0.5, 1.5. We also plot the learned rrebbifurattion phase portrait at μ 0.5 as reference in Figure (a).
during bifurcations? We then interpret the observations in Section 3.2 through the lens of the Poincaré index.
Definition 4.1. (Poincaré index) Let M be an oriented smooth d-manifold and let f -M- + TM be a smooth vector field. Suppose x" E -M is isolated zero of J, i.e. f(x") O and there is neighborhood U of with [xx + for all E U {x'}. Choose an oriented coordi- nate chart around x" and a closed d-ball D C -M centered at x" such that x" is the only zero of f in and f(x) # O for all x OD on the boundary 0D. The Poincaré index of I at x", Ind(j: x"), defined as the topological degree of the map 4 OD 5d-1 D(x) =fxx)/ll(x)lII
Remarkably, even though the model was trained only localized linear-like data, it still reconstructs the broken symmetry, even for je = 0.5, where the double centers are located at v = 0.707 > 0.5 (Figure (c)). As noted earlier, it not surprising that the model identifies the center-to-saddle transition, because a similar bifurcation present in the linear system. If the model learns the local behavior near (0, 0) with respect to parameter modulation, it can naturally capture this transition. What is truly notable, however, is that the model accurately predicts the emergence of the double-wer structure along with the homoclinic loop.
Ind(f. deg EZ.
This integer is independent of the choice of chart and D. Remark 4.1. For = the Poincaré index of a zero x" coincides with the winding number of the vector field around x². Intuitively, this measures how many times and in which direction the vector field rotates around x" one traverses simple closed test contour dD that encircles x* counterclockwise. For f(x) = (fg(x), fx(x)) it is equal to
Another noteworthy observation is that, as jee increases (i.e., Ee increases), the model's double-well topology collapses, reverting single linear saddle point (see Figure 6 (d), and bifurcation diagram Figure 4 (b)). Note that the model from Section 3.1 retains the double-wel well structure throughout the range of > 0, as shown in Figure 4 (a).
1 $ do 1 1f -fdff + fydff
1
Ind xx)x 2T JOD
2x Jan
83 +
Remark 4.2. For = an isolated sink, source, center of orbits, and spiral each have a Poincaré index of +1, whereas a saddle point has a Poincaré index -1-
Additionally, we found that the context-informed model, even when trained on linear system (q.p) p) = (p. μq) with- out the q² term, surprisingly produces spurious double- well structure well. A detailed presentation and dis- cussion of this finding are provided in Section 4.3. These findings suggest the presence of an intrinsic mechanism that drives the model to spontaneously generate a double-well structure, rather than correctly learning the functional form.
Remark 4.3. For M R2, if test contour I contains mul- tiple isolated zeros Xx]. ..... xx), I' is Ind Ind(fx This follows because I can be continuously deformed into a new closed curve I' that consists k small closed loops (aD1, OD2...., aDk). each surrounding nected by two-way bridges (see Figure 7 (a)). Considering that the contributions from the bridges cancel each other out, Ind(f, E Ind(f,x)) holds.
the index of f along
of the zeros X2........ x1.), con-
4. Insights from the Poincaré-Hopp Theorem
I. Poincaré-Hopf Index Theory
Theorem 4.1. (Poincaré-Hopf Theorem) Let M be com- paci, oriented, smooth munifold without boundary. and let f M + TM be a smooth vector field on M with finitely many isolated zeros {x;, N5.x,....,,5 ...,xx,) Then, the sum of the Poincaré indices of f at these zeros is equal to the Euler characteristic x(M) of M: Σ Ind (fxx)) = x(M).
In this section and Appendix A, we briefly introduce the Poincaré index, an integer that can characterize topology of vector field not only locally, but also globally, depending on the choice of a test contour (Strogatz, 2018). 8). It helps answer questions such as: what types of fixed points merge
6

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
-0.6
comeszording (ammeter())) (2)
0.0
0.5
2.0
1 from the newly formed saddle. These additional fixed points serve centers of two orbits (see Figure 6 (b-c)).
(a)
1.0
6)
Local Global
Poincaré Index
&
f(x)
However, if the outer orbit containing fixed points collapses due to a strong perturbation, the constraints from Corollary 4.1 no longer apply, and there is no need to maintain these additional centers (Figure 6 (d)) We computed the Poincaré index as a function of & in Figure 7 (b), which predicts the lifetime of the double-well observed in Figure 4 (b) (highlighted in yellow for both plots).
0.0
-0.5
-1.0
0.10
0.15
0.20
0.25
Context (-9)
Figure 7. (a) Poincaré index of closed orbit r contains multiple fixed points. (b) Computed local (r = 0.1) and global (r = 10)
Poincaré indices for the model in Section 3.2 (see Figure 4 (b)), where r is the radius of a test contour centered at (0, 0).
Based on these observations, in Proposition 4.1, we for- mally derive the condition under which the context-informed NODE exhibits the -ymmetry-breaking bifurcation.
A formal proof can be found in textbooks on differential topology, such as (Milnor & Weaver, 1997). In addition, although Theorem 4.1 is presented in generalized form, the following 2-dimensional statement is more directly ap- plicable to interpret the situation observed in Section 3.2.
Proposition 4.1. Let f R x R - IR2 be a Hamiltonian vector field of the form (( P) (p. be + pP))), where P(g) is smooth function such that P(0) P0) = 0.
Consider another smooth vector field g: .R² x R - R2 for which there exists smooth bijective map Φ R + R such that I and g are 5-close in the C2 sense:
Corollary 4.1. (Poincaré-Hopf for Closed Orbits) Let I: R² + I²² be a smooth vector field and let I be simple (non-self-inters ielffiintersectingg closed orbit of the dynami- cal system x= = f(x). Suppose that all fixed points (all zeros of f) inside I are isolated. Then, the sum of the Poincaré indices of all fixed points inside I is equal to +1: Ex-emn(r) and(f.x'). =
suplig(x; p())) - -(xx 5, xEU
in some open interval μ E (-c,c), where U E R2 is some neighborhood of (0, 0). Suppose that the dynamical sys- tem x = g(x, E) admits at least one closed orbit T that encloses the isolated fixed point x* ~ (0,0) of g. for some neighborhood of Eecit = ((pra)) Then, 9 undergoes the (generalized) symmetry breaking, at least locally near Eerit-
Ex Eint(T)
Proof. A detailed proof can be found in standard text- books, such as (Strogatz, 2). Intuitively, because r is
an actual trajectory of the dynamical system x = fxx), fx) is tangential at every x E (see Figure 7 (a)). It ensures In) is equal to +1. Then, from Re-
See Appendix B for an informal proof. Proposition 4.1 points to intriguing perspective: it only requires that the ground truth f(= undergoes a center- center-to-sadee taddeee bifurca- tion, not necessarily the symmetry-breaking one. Such a transition could, for instance, result from a simple linear system as discussed. However, once the context-informed model f(: learns this local behavior and exhibits an orbit encompassing this local region, the Poincaré-Hopt theorem compels the model to produce the additional fixed points.
mark 4.3, Ind(f.1") = E (nd(f,x)) = +1, where xi, .x....,,;, are fixed points lying in the interior T. D
4.2. Poincaré-Hopf as an Implicit Regularization
Note that any NODE represents a smooth vector field when using smooth activations such as tanh or swish. Then, as- sume that the context-informa model effectively captures the transition where the center of single orbits at (0, 0) trans- forms into a saddle point. In this case, the emergence of the saddle point implies change in the local Poincaré index from +1 1 to- - 1, according to Remark 4.2.
4.3. Hallucinated Broken Symmetry
Inspired from Proposition 4.1, we present an experiment where the context-informed NODE identifies incorrect bi- furcation behavior. In this experiment, we consider learning linear system defined as (q.p) = (p, pq). As briefly explained in Section 3.2, this system has a simple orbit structure for μ < 0 and exhibits a center-to-saddle bifur- cation at (0, 0) when p > Perit = 0. but does not pro- duce the double-well structure. We simulated this linear system over T 2.0 with At = 0.1, using four sam- pled initial conditions from u1[ [/[[-0.3 0.3, )2 for each na e { {-0.35, -0.25, -0.15. - -0.05, 0.05, 0.15, 0.25, 0.35]. as training data (see Appendix F for details). Note that, in this setting, the model explicitly learns from data in the post- bifurcation regime, contrast to the previous experiment.
The context-informed model learns a foliation of closed orbits for { < Eerit from the single-well training data (Figure 6 (a)). Now, assume that the model preserves at least one orbit I encireling (0, 0). under small perturbations in Eerie 3 Then, the summation of the Poincaré indices within I' must be +1 by Corollary 4.1. Thus, to preserve the total index of +1 inside Т, the model must generate additional fixed points whose combined Poincaré index sums to +2, offsetting the
Empirically, we observed that this assumption often holds. It he partially justified hy noting that small perturbation does not completely collapse the foliated orbits of conservative systems.
7

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
nn parameter
Conreponndnng parameter (1)
-0.
0.0
0.5
1,0
1,5
2,0
-0.
0.
0.5
1.0
1.5
20
0.10
1.0
Canonical momentum (p)
2-
(mose))
(GT)
(mode)
Sudde (modee
-Satche (01)
Secdle (indae) Saddo (GT)
Location points
0.05
0.5
0.0
0.00
0-
2-00000
0.0
00000
-0.1
-0.05
-0.5
Fres form. C,0005
Hom Boniar ( 0,00007
(d)
-02
-066
-0.3
0.0
0.3
0.6
-0.3 -0.2 -0,1 0.0 0.1 0.2 0.3
-0.1
0.0
0.1
0.2
0.3
0.4
-0.1
0.0
0.1
0.2
0.3
Canonical position (q)
Canonical position (q)
Context (5)
Context (-E)
Figure 8. Phase portraits of (a) free-form and (b) Hamiltonian models trained with (q, p) = (p, pg),, and constructed in zero-shot manner
by linearly extrapolating &e, correspond to ee 0.01. Bifurcation diagrams of (c) free-form and (d) Hamiltonian models.
We evaluated two different architectures of context- informed NODEs: the default free-form and its Hamiltonian version (Greydanus et al., 2019). The latter is defined as
Localion of fixed paints (q "
---C
C 2)
C
1Su 3) 3:5
5 (Bun 1)
--
0
0 110000
C f)
x = (q. - f(x 0o) = (ap (x; 0z), - 0,,4(xxx,)))
8 (Sub
where 7((x; = 7((x; 0, + W. is a context-informet learnable Hamiltonian (see Appendix F for details). Since (g. p) = (p. pq) is also a Hamiltonian system, applying the Hamiltonian bias appears fundamentally appropriate.
-2
C (Su0 3)
(a)
(b)
8 (Sub 3)
-0.5 0.0
0.5
1.0
1.5
2.0
-0.5
0.0
0.5
1.0
1.5
2.0
Parameter (v)
Parameter (v)
Figure 9. Ensembled bifurcation diagrams trained on (a) the double-well (Section 3.1) and (b) the linear system (Section 4.3).
Figure 8 (a-b) illustrates the phase portraits explored by context-informed models trained on the linear system. As shown, despite being trained with both pre- and post- bifurcation data, the free-form and Hamiltonian NODES alike misinterpret the bifurcation, incorrectly identifying it as a gummetry-breakiing transition, resulting spurious double-well structure. Notably, the foulle-well structure in the free-form model rapidly collapses with respect to & (highlighted in purple), reverting to normal saddle point (Figure 8 (c)), whereas the Hamiltonian model retains the misidentified structure more persistently (Figure 8 (d)).
Each ensemble contains three sub-NODEs trained independently on bootstrapped datasets with different initializations.
scenarios: involving correctly identified bifurcation (from Section 3.1) and the other involving a hallucinated one (from earlier results). In Figure (a), the symmetry- breaking is correctly captured, and the estimated variance remains minimal, indicating stable detection of the bifurea- tion. In contrast, the ground true bifurcation in Figure (b) is center-to-saddle rectly generates two centers of orbits within 0 < μ < 0.5, representing hallucinated behavior. This results in low vari- ance (5.343 x 10 at μ 0.2) near the actual transition but significantly higher variance (0.4605 at μ = 0.2) in the spurious region, reflecting its structural instability.
q" O. However, the model incor-
This phenomenon can be explained using Proposition 4.1: both models successfully learn the local center-to-saddle bifurcation of the linear system. They also preserve outer orbit structure that encapsulates the local training domain at ~ Eerit- As result, they hallucinate extra fixed points, mimicking broken symmetry. The Hamiltonian model, how- ever, preserves the outer orbit more robustly due to its sym- plectic structure (Strogatz, 2018), leading to the persistent double-well configuration. Appendix G, we also provide an experiment conducted under the standard pre-bifurcation training, which yielded the similar conclusion.
5. Applications to Broader Systems
5.1. Identifying the Cusp Catastrophe
Catastrophe theory (Thom, 1977; Zeeman, 1979) describes how small, continuous changes in parameters can lead to sudden, discontinuous shifts in system behavior. Among its five fundamental types (Saunders, 1980), the cusp catas- trophe is one of the most well-known and extensively stud- ied. The cusp catastrophe is a codimension-2 - bifurcation system given by * (q. = (p, v + pa -² kp). with parameters (u.v) and fixed damping = 0.5. It is simple extension of (3), but the inclusion of the term explicitly breaks symmetry, making it a canonical model for studying hysteresis phenomena. We investigated whether context-informet NODEs can identify the cusp bifurcation and hysteresis loop, despite being trained only on the pre-
4.4. Diagnosing Hallucinated Bifurcations
Previous findings highlight the importance of determining whether model hallucinates bifurcations. We propose simple yet effective criterion: assessing the variance in bi- furcation diagrams generated from multiple independent training runs, similar to the deep ensemble method (Lak- shminarayanan et al., 2017). Figure 9 illustrates the esti- mated variance from 3-fold bootstrapped ensembles in two
8

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Valvoriio sator
Horosbble Tisking
laa change
0 pare
porto (Poene 10 15
(b) (ust 12.00120. 22.-05)
IAYI (2.....)
0.4
0.4
(b)
0.80
Nyzasmo kep
Double-well L Single-well
Double-well
Single-well
0.75
Locaon of frae pomns
Parameter
0.60
O
0.0
0.0
0.46
-02
Double
Trp
-0.2
Double-well
Triplewwll
0:30
0.15
-0.4
-0.4
Pranndr d) -vanneeee An
-04
-0.2
0.0
02
0.4
-3.4
-0.2
0.0
0.2
0.4
Parameter
Parameter
Figure II. Consistency scores (a = 0.5) for learning (4) with (a)
Figure 10. (a) Bifurcation surface generated by the context- informed model trained on cusp bifurcation system. (b) Compar- ison between the phase portraits of the ground truth (black lines) and the learned model (red lines) near the catastrophic transition.
vanilla and (b) regularized models, averaged over five runs.
for A + 0. In this case, the system must satisfy the con- ditions stated in Corollary 4.1. For nonzero K, there is no strict guarantee that the total Poincaré index will always be constrained to +1. However, the introduction of a lin- ear damping does not completely disrupt the overall phase topology. Consequently, the distribution of Poincaré indices is expected to be preserved. Inspired by this, we propose the following Poincaré Hopf regularization:
bifurcation regime. Specifically, training was conducted using parameters (HL v²") E {-2.0, -1.5, -1.0, -0.53" which lie entirely within the monostable region. For each training parameters, four training trajectories were gener- ated with ((O) 4((-2.0, 2.0)3). T = 2.0 and At = 0.1. After training, the bifurcation surface was gen- erated over mesh grid (1t ) E [-2.0, 2.0)², using the neural-adjoing method. Figures 10 (a) and (b) respec- tively show the identified bifurcation surface and the phase portraits of the learned vector fields near the catastrophic transition (see Appendix H for other cases). These results demonstrate that the model successfully captures the com- plicated bifurcation structure without explicitly learning it.
(Pn(Oe, W, Ee) [Ind(f(: 0. + WG.), [PH) - XPH 1
(5)
where IH is a proper test contour and XPH is the desired index (see Appendix for detailed description of (5)).
We trained the conventional model and its regularized coun- terpart with (5) on the system (4). For training data, we simu- lated (4) over T = 2.0 with At = 0.1, using four initial con- ditions sampled from (+) v' ~ ((00200 2.0]²), for each combination of (a (85) E (-0.4, -0.22 0.2, 0.4) To ensure the physical relevance of (4), = 0.05 and k = 0.3 were fixed. After training, we evaluated each model's ability to capture the long-term behavior and phase topology by assessing the consistency between its limit sets and those of the ground truth dynamics (Göring al., 2024). Specifically, 32 trajectories were sampled and simulated during T = 100.0, for each parameter vector on mesh grid over tasst , 3!eest E [-0.5, 0.5]2. For each (atest, ; gteest, ), the corresponding gtest was constructed using linear regres- sion fitted on (a (tr)), and then used as input to the model to simulate long-term trajectories also. the ground truth and model dynamics converge within a distance of a, we consider them consistent. Figure 11 confirms (5) enhances the long-term predictability across the parameter plane. De- tailed results are in Appendix J, with exhaustive ablation studies of the proposed regularization in Appendix K.
5.2. Identifying the Landau-Khalatnikov Theory
The LK theory (Landau & Khalatnikov, 1954) is widely used mathematical framework for describing phase tran- sition dynamics in various physical systems, particularly ferroelectric materials (Khan et al., 2015; Lo, 2003; Shi et al., 2016). The dynamics are governed by the following second-order ODE / + vab = F(bb) where & is an order parameter characterizes the phase of the system (e.g.. polarization of ferroelectric materials), x is fixed damping coefficient, and F(v) is the free energy expressed nomial expansion F(b)) = v2 + & v/ 2/06 +.... The coefficients a, Y,.. are phenomenological parameters typically expressed as functions of some physical quantities like temperature, e.g.. a = aT - ). The phase tran- sition of so occurs as these parameters change with T. The LK model can be formulated in Hamiltonian-lin form as
poly-
(1b, v/) (V', anb Bub² - 778 + +...), (4)
except for the damping term -KyL': if x - O, then the sys- tem approaches precise Hamiltonian equations of motion.
6. Conclusions and Limitations
We show that NODEs can identify symmetry- bi- furcations solely from symmetric data, interpreted through the Poincaré index. We also propose novel regularization inspired by the Poincaré Hopf theorem. While the current method is limited to 2-dimensional flows with nearly closed orbits, Theorem 4.1 reveals that the concept can extend to a broad class of flows, paving the way for further research.
The exact form of F, such how many orders should be included, is often unknown. However, one essential require- ment is that F(e) must not diverge to negative infinity to preserve physical consistency. This constraint ensures that the trajectories of (1,, ) in (4) remain bounded, forming a closed outer orbit that encapsulates the system, especially

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Acknowledgements
Brenner, M., Weber, E., Koppe, G., and Durstewitz, D. Learning interpretable hierarchical dynamical systems models from time series data. In International Conference on Learning Representations, 2025.
We thank Jabir Bin Jahangir, Jisu Ryu and Young-Gu Kim for their helpful comments on an earlier draft of this work. Changwook Jeong acknowledges support by the National Research Foundation of Korea (NRF), funded by the Min- istry of Science and ICT (MSIT) (RS-2023-00256666 RS- 2024-00458251, RS-2024-00408 180); by the Ministry of Trade, Industry and Energy (MOTIE) (RS-2023-0023 1996,,, P0023703); by the Institute for Information & Communica- tions Technology Planning& & Evaluation (ITT)), funded by MSIT (RS-2020-11201 136,, AIGS); by the research project fund of UNIST (1 1.255004.01; 1.220104.01); and by Sam- sung Electronics. Muhammad Ashraful Alam acknowledges support by the Jai N. Gupta Distinguished Chair Professor- ship fund.
Brunton, L., Proctor, J. L., and Kutz, J. N. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the Na- tional Academy of Sciences, 113(15):3932-3931 2016.
Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K. Neural ordinary differential equations. Advunces in Neural Information Processing Systems, 31:6572-6583, 2018.
Crawford, J. D. Introduction to bifurcation theory. Reviews of Modern Physics, 63(4):991, 1991.
Impact Statement
Fotiadis, S., Valencia, M. L., Hu, S.,, Garasto, S., Cantwell, C. D., and Bharath, A. A. Disentangled generative mod- els for robust prediction of system dynamics. In Inter- national Conference on Machine Learning, pp. 10222- 10248. PMLR, 2023.
This paper explores the potential of NODE models in ad- dressing OOD challenges in learning dynamics, focusing on bifurcation identification and generalization through the lens of topology and dynamical systems theory. The content is primarily theoretical, and no ethical concerns have been identified.
Gareía Pérez, J., Sanches, L., Ghadami, A., Michon, G., and Epureanu, B. I. A physics-informed data-driven ap- proach for forecasting bifurcations in dynamical systems. Nonlinear Dynamics, 111(1))11773-111
2023.
References
Ghadami, A. and Epureanu, B. Forecasting critical points and post-critical limit cycles in nonlinear oscillatory sys- tems using pre-critical transient responses. International Journal of Non-Linear Mechanics, 101: 156, 2018.
Alam, M. A. and Zagni, N. Physics of phase transition switches. IEEE Electron Devices Reviews, 2024.
Aranson, I. and Kramer, The world of the complex Ginzburg-Landau equation. Reviews of Modern Physics, 74(1):99, 2002.
Göring, N. A., Hess, F., Brenner, M., Monfared, Z., and Durstewitz, D. Out-of-domain generalization in dynami- cal systems reconstruction. In International Conference Machine Learning, pp. 1607 111114. PMLR, 2024.
Amold, V. I. Mathematical methods of classical mechanics, volume 60. Springer Science & Business Media, 2013.
Greydanus, S., Dzamba, M., and Yosinski, J. Hamiltonian neural networks. Advances in Neural Informatio Pro- cessing Systems, 32, 2019.
Arnold, V. I., Afrajmovich, V., yashenko,, Y. S., and nikov,, L. Dynamical systems V: bifurcation the- ory and catastrophe theory, volume 5. Springer Science & Business Media, 2013.
Guillemin, V. and Pollack. A. Differential topology, volume 370. American Mathematical Soc., 2010.
Bernstein, J. Spontaneous symmetry breaking, gauge theo- ries, the higgs mechanism and all that. Reviews of Modern Physicx, 46( (1)7,,111 1974.
Ha, D.,. Dai, A. M., and Le, Q. V. Hyperneworks.. In International Conference on Learning Representations, 2017.
Brasselet, J.PP., Seade, J., and Suwa, T. Vector fields on singular varieties, volume 1987. Springer Science & Business Media, 2009.
Hamel, P., Haddadi, S., Raineri, F,, Monnier, P., Beaudoin, G., Sagnes, I., Levenson, A., and Yacomotti, A. M. Spon- taneous mirror-symmetry breaking in coupled photonic- crystal nanolasers. Nature Photonics, (5):31133 2015.
Brenner, M., Hess, F. Mikhaeil, J. M., Bereska, L. F., Mon- fared, Z., Kuo, P.-C., and Durstewitz, D. Tractable den- dritic RNNs for reconstructing nonlinear dynamical sys- tems. In International Conference on Machine Learning, pp. 2292-2320. PMLR, 2022.
Hopf, H. Vektorfelder in n-dimensionalen mannig- faltigkeiten. Mathematische Annalen, 96:225 2249, 1927.
10

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Huh, I,, Yang, E., Hwang, S. J., and Shin, J. Time-reversal symmetric ODE network. Advances in Neural Informa- tion Processing Systems, 33: 19016-190227,, 2020.
Perez, E., Strub, E,, De Vries, H., Dumoulin, V., and Courville, A. FiLM: Visual reasoning with a general conditioning layer. In Proceedings of the AAAI Confer- ence on Artificial Intelligence, 2018.
Kalia, M., Brunton, S. L.,, Meijer, H. G., Brune, C., and Kutz, J. N. Learning normal form autoencoders for data-driven discovery of universal, parametee-deeeeddet governing equations. arXiv preprint arXiv:2106.05102. 2021.
Ramachandran, P., Zoph. B., and Le, Q. V. Searching for activation functions. arXiv preprint arXiv:1710.05941, 2017.
Khan, A. I., Chatterjee, K.. Wang, B., Drapcho, S., You, L., Serrao, C., Bakaul, S. R., Ramesh, R., and Salahuddin, S. Negative capacitance in a ferroelectric capacitor. Nature Materials, 14(2):182-188) 186, 2015.
Raya, G. and Ambrogioni, L. Spontaneous symmetry break- ing in generative diffusion models. Advances in Neural Information Processing Systems, 36:66377-66389, 2023.
Ren, S., Padilla, W., and Malof, J. Benchmarking deep inverse models over time, and the neural-adjoin, method. Advances in Neural Information Processing Systems, 33: 38 48, 2020.
Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015.
Kirchmeyer, M. Yin, Y.. Donå, J., Baskiotis, N., Rako- tomamonjy, A., and Gallinari, P. Generalizing to new physical systems via context-informed dynamics model. In International Conference on Machine Learning, pp. 11283-11301. PMLR, 2022.
Rubanova, Y., Chen, R. T., and Duvenaud, D. K. Latent ordinary differential equations for irreulaarl-ssapple time series. Advances in Neural Information Processing Systems, 32, 2019.
Saunders, P. T. An introduction catastrophe theory. Cam- bridge University Press, 1980.
Lakshminarayanan, B., Pritzel, A., and Blundell, C. Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in Neural Information Processing Systems, 30, 2017.
Shi, P.-P., Tang, Y.-Y, Li, P.-F., Liao, W.-Q.. Wang. Z.- X., Ye, Q., and Xiong, R.-G. Symmetry breaking in molecular ferroelectrics. Chemical Society Reviews, 45 (14):381 1-3827, 2016.
Lamb, A. M., Goyal, A., Zhang. Y., Zhang, S., Courville, A. C., and Bengio, Y. Professor forcing: A algorithm for training recurrent networks. Advances in Neural In- formation Processing Systems, 29, 2016.
Song, Y., Keller, T. A., Sebe, N., and Welling, M. La- tent traversals in generative models as potential flows. International Conference on Machine Learning, pp. 32288 32303. PMLR, 2023.
Landau, I.. and Khalatnikov, I. On the anomalous absorption of sound near second order phase transition point. Dokl, Akad. Nauk SSSR, 96(469):25, 1954.
Strogatz, S. H. Nonlinear dynamics and chaos: with appli- cations to physics, biology, chemistry, and engineering. CRC press, 2018.
Li, S. and Yang, Y. Data-driven modeling of bifurcation sys- tems by learning the bifurcation parameter generalization. Nonlinear Dynamics, pp. 1-12, 2024.
Thom, R. Structural stability, catastrophe theory, and ap- plied mathematics. SIAM review, 19(2):189-20) 1977.
Lo, V. C. Simulation of thickness effect in thin ferroelec- tric films using Banndau-Khalatike theory. Journal of Applied Physics, 94(5):3353-3333 2003.
Vlachas, P. R., Byeon, W., Wan, Z. Y., Sapsis, T. P., and Koumoutsakos, P. Data-driven forecasting of high- dimensional chaotic systems with long short-term mem- ory networks. Proceedings of the Royal Society A: Math- ematical, Physical and Engineering Sciences, 474(2213): 20170844, 2018.
Milnor, J. W. and Weaver, D. W. Topology from the differen- tiable viewpoint, volume 21. Princeton University Press, 1997.
Mouli, S. C., Alam, M., and Ribeiro, B. MetaPhysiCa: Im- proving OOD robustness in physics-informa machine learning. In International Conference on Learning Rep- resentations, 2024.
Wei, G., Huang, Y., Duan, C., Song, Y., and Du, Y. Navigat- ing chemical space with latent flows. Advances in Neural Information Processing Systems, 37, 2024.
Nzoyem, R. D., Barton, D. A., and Deakin, T. Neural context flows for meta-learning of dynamical systems. In International Conference on Learning Representations, 2025.
Ye, H., Xie, C., Cai, T., Li, R., Li, Z., and Wang. L. Towards a theoretical framework of distrriuitn generaliza- tion. Advances in Neural Information Processing Systems, 34:23519-299 23531,, 2021.

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Yin, Y., Ayed, I., de Bézenac, E., Baskiotis, N., and Galli- nari, P. LEADS: Learning dynamical systems that gener- alize across environments. Advances in Neural Informa- tion Processing Systems, 34:7561-7573, 2021a.
Yin, Y., Le Guen, V., Dona, J., de Bézenac, E., Ayed, I., Thome, N., and Gallinari, Augmenting physical models with deep networks for complex dynamics forecasting. Journal of Statistical Mechanics: Theory and Experiment, 2021(12): 124012, 2021b.
Zeeman, E. C. Catastrophe theory. In Structural Slability in Physics: Proceedings of Two International Symposia on Applications of Catastrophe Theory and Topological Concepts in Physics Tübingen, Fed. Rep. of Germany, May 2-6 and December 1 11-14, 1978, pp. 12-22. Springer, 1979.
12

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
A. A Brief Overview of the Poincaré-Hopf Index Theorem
This section provides a brief introduction to the Poincaré-Hopf index theory. For a more comprehensive treatment, we refer readers standard textbooks on dynamical systems theory (Strogatz, 2018) and differential topology (Guillemin Pollack, 2010; Milnor& & Weaver, 1997).
Definition A.1. (Poincaré index) Let M be an oriented smooth di d-manifold and let . + TM be a smooth vector field. Suppose x' M is an isolated zero of f. i.e. f(x") and there is neighborhood U of x" with f(x) + 0 for all EU {x*}. Choose an oriented coordinate chart around x" and closed d-bull D C -M centered at x" such that x" is the only zero of f in D and f(x) & foror all x E OD on the boundary OD. The Poincaré index of f at x', Ind(f,x") is defined as the topological degree of the map:
f(x)
(x)
:0 + = 11((x)]
that is,
Ind(f[x"))= deg (b) E Z.
This integer is independent of the choice of coordinate chart and the ball D, and thus defines a topological invariant associated to the isolated zero x". A visualization of the Poincaré index in TR. is shown in Figure 12.
Remark A.1. Intuitively, the map ((x) fxx) normalizes the vector field, capturing only the direction of eacc vector and discarding its magnitude. Thus, (xx) maps each point x on ÔD gá-1 to point on the unit sphere gd-1, representing the direction of [x).. The topological degree deg( ()) is an integer obtained by counting, with orientation signs, the preimages in 0D of any regular value in Sd-1; deg(4) Exce-1(s) sign(detD@) (x)) for a regular value y gd-1.
This integer can be interpreted as an oriented wrapping number of OD around the sphere Sd-1 under the map 4. In this sense, Definition A.I generalizes the classical winding number of a planar vector field in TR2 to ligheer-dmeessonall cases.
Remark A.2. For an isolated, non-degenerate zero x" of vector field f, the Poincaré index can be computed as
Ind(f, = (dig(det: D Dx())))
where Df(x*) denotes the Jacobian matrix of f evaluated x". This means that the index of an isolated and non-degenerate zero must be either +1 or 11 Intuitively, the sign of the Jacobian determinant indicates whether the vector field locally preserves orientation (Ind(J, [f, x*) ) = +1) or reverses it (Ind(J, = -1) near the zero.
Remark A.3. In the language of differential topology, the topological degree of the map tb : 0D - Sd-1 can be computed using the pullback of a volume form. Specifically,
/ w'w,
(e)
w
Jz
BD
where w is the volume form on the unit sphere gd-1 and w denotes its pullback via the map Q. Since the integral of the volume form over the unit (d 1)ssherr is simply given by Jga-1 w = Vol(Se-2) = 2/2 2 the Poincaré index can be computed as
T(d/2)
Ind(f,DD) = deg(() (b) Vol(Sd-1) JaD
F' JaD
(6)
2md/2
Remark A.4. When evaluating (6), the expression is well-defined and yields correct result as long as the test contour OD does not pass through zero of the vector field; that is, f(x) # O foror all x E OD. If multiple isolated zeros lie inside the region enclosed by OD, then (6) computes the total Poincaré index of all those zeros, defined as the sum of the indices of all enclosed zeros:
Ind(f.OD) = Σ Ind(f.x").
X Eint((D)
For = ?R? this can be understood intuitively as follows:
Suppose the contour 0D encloses several isolated zeros {xi, xxixxx...,,x?] The contour 0D can be continuously deformed into new closed curve 0D' consisting of k small loops DD,,,. ,,Dk,, each encircling exactly one of the zeros Xx.,, and connected by two-way narrow bridges. Since the contributions from these connecting bridges cancel out pairwise, the total index remains unchanged. This intuitive picture extends naturally to dimenssoal cases, where the topological invariance of the index ensures that the total index depends only on the enclosed zeros, not on the specific shape of the test surface.
13

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
A vector flldd flx) on 60 55 (a)
The itted 5 via exx: (iv/ - (b)
A vector teld f0o) on ab
The lited S via eix): fXxV/ffi ::
Ind +0
Ind +1
A vector field Axt on aD 5 (c)
The lited 5 va eix: Ax:::A/ARxx :5 St
(d)
A vector feld fix) on aD 5
The lited 5 via 6(x): AXV/[fIX :5
Ind +2
Ind +1
(e)
vecior flxt aD
The iieed va eix! (xx/Xx) 55
A veclor Tux) on dD 5²
5 via e(x) AXY]|AX
Ind
=
Figure Visualization six vector field configurations f(x) on the circular test contour 0D = S' (left column), and their corresponding images on g1 via the map ((x) (x) (right column). The test contour aD is traversed in the counterclockwise direction throughout. A color gradient from blue to red indicates the progression along the path. (a) uniform vector field along OD maps the entire circle to single point 8) the resulting Poincaré index is 0. (b) A radial vector field makes the contour wrap around S! in the counterclockwise direction, giving an index of +1. (c) Traversing 0D once counterclockwise, the vector field winds twice about the origin in the counterclockwise direction, the image covers S' twice and the index is +2. (d) A tangential vector field rotates counterclockwise along the contour, again producing an index of +1. (e) A saddle-like configuration winds once clockwise relative to the counter-clockwise traversal, producing index of 1. (f) A more complex field winds twice clockwise, resulting index of -2.
Remark A.5. From the definition of the Poincaré index (6), its two-dimensional version can be computed as
$ -Sndf + Sydfe
Ind(J,DD)
JaD f2+
While this expression can be derived directly from (6), it is more straightforward obtain it using the standard definition of the winding number. For real vector field f(x (q-r)) (Ją(x), Sp (x)), let its complex representation be F(x) fx(x) + +x)((x). Because d0 Im(dF/F), the index be written
1 f do
1 2π
f Jan
dF
1
(dfa + idf - iifp)
1 2
-f fadfo
Ind(f, (D)
2n Jop
Im
2 JD
1m
12+
JaD
+
Theorem A.1. (Poincaré-Hopf Theorem) Let M be compact, oriented, smooth manifold without boundary, and let M + TM be smooth vector field on M with finitely many isolated zeros X2,,x,,..,, Then, the sum of the Poincaré indices of ff at these zeros is equal to the Euler characteristic x(M) of M:
Ind(f.x)) x(M).
14

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Corollary A.I. (Poincaré-Hopf for Closed Orbits) Let f R + +R² be smooth vector field and let T be a simple iton-self-intersecting closed orbit of the dynamical system x f(x). Suppose that all fixed points (all zeros of I) inside T are isolated. Then, the sum of the Poincaré indices of all fixed points inside I' is equal to +1:
Σ
Ind(J,x")= =1.
Proof. Because M ?R and I is actual trajectory of the dynamical system x f(x), fxx) is tangential at every T. ensures Ind(f,T) T) equal to +1. Then, from Remark A.4, Ind(S.T) > Ind(f,x} x²) = +1, where xi, x2...... X fixed points lying in the interior of l'
B. Mathematical Details and Proofs
Assumption B.1. Consider a separable Hamiltonian H IR² x R2 + R:
p²
1a
f + +C,
H(g. P. = (()) + V2(
-
2
2
where P IR + R. is smooth function satisfying P(0) P()) assumed to be Hamiltonian vector field of H. and expressed as
0. The ground truth vector field f R2 x R - R2
f(g.p:)) = (.H(g. p: p)
07(q. 0H(q.P: = (p, un + P(4)),
(g. = =
Op
0q
Remark B.1. The Hamiltonian vector field f defined in Assumption B.1 has an isolated fixed point at (0,0) 0) because
(0 = (0, P(0)) (). Note that the Jacobian of f at (0, is
O9q 0,4 = (N (u O(P()) !) - >)(0,0:m) -00 (1 1) Oup 0pd)
0
Df(g.r:
The characteristic polynomial is
det()(())) -1) =
x=
Therefore, the fixed point (0, 0) is a center for μ < 01A = =ivllu) and saddle for > 0 = + +V(A)) which means the local bifurcation is the center-to- saddle type.
Remark B.2. Consider P(g) = 0. It satisfies Assumption B. In this case, (0.0) is the only fixed point. The system exhibits a center-to-suue saddle bifurcation at (0, 0) simply.
Remark B.3. Consider P(g) = -q It satisfies Assumption B.I. The fixed points of f are determined as follows:
[4 = (p, g(u -2)) 0),
(g.p) (0,0), (+v/,,
The additional fixed points (q.P) (+1,,, exist only for μ > O. Thus, the system undergoes 1 m 3 bifurcation O, generating two additional fixed points. For these fixed points, the Jacobian is given by
(n-
0 39²
1)
0 D)(+/T,0;,)
2p 0)
Di(a.P.))
The characteristic polynomial for these fixed points is
Xs(A, p) + x= -= p1pp
where μ > Thus, the emerging fixed points are center-like, and the system exhibits a commetry-breakiin bifurcation.
Lemma B.1. Any fixed point of a Hamiltonian system in a 2-dimensional phase space must either be a center of closed orbits or a saddle point.
15

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Proof. is direct application of the volume-preserving property of Hamiltonian dynamics (Arnold, 2013). The divergence of Hamiltonian system holds
0q + 0i.
0H
a 0H
V f(g.p)
0,
Oq Op
dq 0p
Op dq
for any ()) and je, regardless of the specific form of H. This implies that the trace of the Jacobian must be zero:
Ir(DS(u:P) p)) = v. f(.
This ensures that the characteristics equation at (q. p) is given by
xy(A, = X² + det (DI(g,)))-
Consequently, the fixed points (g.p) p) of Hamiltonian system must fall into one of the following two categories: a center when det (q. p)) det (Df(1 p)) (g((q.P) (q)
saddle when det (Df(g. p)) 0. Moreover, in the case of separable Hamiltonians,
Lemma B.2. Let g R22 x R IR2 is smooth vector field, Assume there exists smooth bijective map R - R such that f and are 8-close in the C. sense:
suplig o(H)) f(x; + sup Dg(x; 0(p)) DD(x S xeU
xeU
in some open interval E (-e,e), where f is as defined in Assumption B.1 and U E I2 is some neighborhood of (0, 0). Then, there exists a fixed point x" = (0,0) which a near-center when μ < Perit = and a saddle when p. > Perit 2
Proof. It follows directly from (Crawford, 1991). From the statement, we can rewrite g as
g(x; (())) [(x; p) (x:;), ||(x; < 0
Similarly, the Jacobians are given by
Dg(x; (())) = (xx + D)(x;P), |[D( S O
First, we want to show that for each small μ , there exists fixed point of g(: p(1)) closed to (0, Note that
g(0: ((1)) = f(: μ1 + A(0:)) A(0: μ))
By denoting G(x, μ) = g(x; 6(4)), we have [GG(O, - C8. Consider the Jacobians at x = (0,0):
D0 ((x) (0: DA
Because (ct( = 44, D0((: μ1) is invertible for all μ + 0. Then, for a sufficiently small 8,
det (g(0) = -4 + + O()) 40,
thus D9(x; remains invertible and μ + 0. Therefore, from the given conditions
||G[0, A S Cτ, det(DGG(O + O,
the continuity and implicit function theorem imply that for μ * there exists an unique solution G(x(8, (), g(x(6, p), ((H)) 0 satisfying x(6, μ)) - (0, 0) as 8 - 0. directly gives g has an isolated fixed point x(6, p)) close to (0, 0) whenever μ + O.
We now show that (x(6, p), μ)) of g has unique smooth curve (A, 1)) of eigenvalues near Xg for μ in a neighborhood of Ho # O. Consider the Jacobian at x(8, p):
Dg(x(6, p); ((n)) J6. = D)((;)) + =
The characteristic polynomial of J, (a) is given by
(g u1 x2 - -- (T. - 1T(2(600))
16

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
where ir(Ts) < O()). For (A, p1) (Ao, satisfying ( Mo) = 0 and (Ao: 1o) + (0,0), have Xo (Ao. no: 6) = -tr(Ja (6, Ho))Ao,
thus (oo 6)1 S C8. Then, consider the partial derivative of Xa with respect to A at (Ao,
0Xg (Ao: Po: 8) 2X tr trT (8, n)):
0x
which remains nonzero for sufficiently small 8. Observe the given conditions
Xa ax
Ixg Mo. 5)1 S C
(Ao, No. 22 C6 > 0.
Again, by the continuity and implicit function theorem, there exists unique curve (A(o), p(o)) such that Xg (())) 0. (()5), + (Au, as 8 + O for # 0 and a sufficiently small 5.
Now consider fa - 11| < 0 case. From the characteristic equation Xy discussed above,
Xo p. 0) X² + Il (r( n))a = == tr(Ja(o, = =ivlpl [p]
2
thus Re(A) ( Therefore, for μ O, the fixed point x" 0 of g behaves like a center, although it may exhibit infinitesimal growth damping depending the sign of tr (I[A). Next, for μ |a|| > O, we have
(A, 14, 6) = tr[]_t((J = A= tr tr(Ja (r( 1)) = +Vlll.
2
Given that a small perturbation cannot alter the signs of X, the fixed point x* 0 of g is classified as saddle point for μ 0, though the rate of convergence or divergence may vary slightly along each eigenvector direction.
Proposition B.1. Let f and g satisfy the assumptions stated in Lemma B.2. In addition, suppose that the dynamical system
x = ==g(x;;= = x(H)) admits al least one closed orbit I that encloses isolated fixed point XI = (0, 0) for a neighborhood of (. Then, 9 undergoes a (generalized) mommetry-breaking bifurcation, at least locally near Earit-
Proof. From Lemma B.2, the system x = g(x, (())) has a fixed point (p) near (0, 0), and undergoes a center-to-saddle bifurcation. It means that, as μ passes through Perit ~ O, the Poincaré index at (1), Ind(g, (1)) changes from Ind(g, < )(rit)= = +1 (a center) to x51 > Meert) = -1 (a saddle), resulting in net index change of -2.
By assumption, there exists closed orbit I encircling x" (p) for neighborhood of ((Perit). According to Corollary A. the sum of Poincaré indices of all fixed points enclosed by I' must be +1 near μ Merlt. Consequently, to preserve the total index of +1 within I, there must emerge additional fixed points, xi,...,,?; X within I, whose combined Poincaré indices sum to +2:
Ind(f. x;) = +2.
As established in Lemma B. a Hamiltonian system, thus its divergence V f(x; p)) for all x and μ. This implies that the trace of the Jacobian (r( wi)) = 0 for all x and μ. This constraint ensures that the eigenvalues and accordingly the types of fixed points must fall into one of the following two categories: (1) purely imaginary (a center) or (2) real with opposite signs (a saddle).
Then, as established in Lemma B.2, because g is 5-close to f in the C' sense, it ensures || (xx: ا1) = |tr(t +4) S (8) within U. For sufficiently small 8, this ensures that no fixed points with large real parts in their cigenvalues can spontaneously arise. Consequently, any non-saddle fixed points of g in must be near-centers, retaining a Poincaré index of +1 by continuity. Thus, under generic scenario, the newly created fixed points x xi,...,X, x are two centers, contributing total index of +2.
Definition B.1. (Parameter OOD) Let P R" be a parameter space and let B C P be bifurcation set. The complement P \ B admits a decomposition into I, disjoint connected components: P B = U P., where each P; is muximal connected subdomain. Intuitively, each P; is one qualitatively uniform parameter regime. Then, parameter OOD condition in learning dynamics arises when the support of training distribution is supp(p" (μ1) C Pt, but the support of the test distribution supp(p'.st test (μ)) S Pags. Equivalently, there exists no continuous path ? : [0, 1] + P B connecting p. and n test without crossing B.
17

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Ground (u= -2.0)
Ground -1.75)
Ground truth (u= -1.5)
Ground tuuth (u
Canonical momentum (P)
0-
0
2
-2
-3
-3
Canonical position (4)
Canonical position (9)
Canonical position (4)
Canonical position (4)
(u
Ground (y -0.75)
Ground (j = -0.5)
Ground (u = -0.25)
Canonical momenium (p)
2
1
0
D
-1
-1
-3
-3
-3
Canonical position (4)
Canonical position (q)
Canonical position (q)
Canonical position (4)
Figure 13. Examples of training trajectories from the experiment descrihed Section 3.1.
C. Automatic Generation of Bifurcation Diagrams
To automatically generate the bifurcation diagram, the fixed points of the learned context-informed NODE must be determined for each variable ç. This corresponds exactly to solving the inverse problem of f(x; {) given by
x" = f-=f=((0:))
for given ç. Because / is fully differentiable, the neural-adjoint method (Ren et al., 2020) provides one of the simplest solutions to this inverse problem. The neural-adjoint method employs gradient descent to iteratively update the candidate x! over iterations, aiming to satisfy fxx : 0:
x -
where is a learning rate and C(f(x).)) is a conventional loss function, such as MSE, designed for inverse problems of the form f (y) = x. In this case, since we aim to find fixed points, y = and the loss function simplifies to the squared L2 norm. This procedure performed on a batch of size with initial values (=) x" ~ p(x"), where is a prior distribution, typically defined uniform distribution (NN) over an initial domain N. In our experiment, we set K 1, 000, 7 10-2, N 100, and p(x") ([[22.0, 2.0 x [-2.0,
After completing the iteration procedure (i K), discard points where the squared L2 norm Ilf(x* E exceeds threshold c. For the remaining candidates x", we compute the Jacobians (x* E). Using the eigenvalues of Df (x**x)) we classify the fixed point candidates based on linearization theory (Strogatz, 2018). Specifically, if the eigenvalues A and A2 are both real with opposite signs, the fixed point classified as a saddle point. If X and X2 are purely imaginary and have opposite signs, the fixed point is classified as center of orbits and so on. For further details on the classification rule, refer to (Strogatz, 2018). After classifying all found fixed points, if any two points within distance have the same classification label, then they merged by averaging. We used 10-8 and 8 10-2 for the experiments.
D. Experiment Details: Section 3.1
Data preparation. We trained the context-informed NODEs (2) on the Hamiltonian system described Remark B.3 under the pre-bifurcation regime exclusively, Specifically, we randomly sampled four initial conditions from the uniform distribu- tion (g(0), ([_22.0, 2.0) 2-2.0 2.0)) for each of the eight parameter values n E {-2.0, -1.75, - -1.5, -1.25, -070, 0.5, -0.25). For each sampled initial condition, we simulated the dynamics with a time horizon T 2.0 and time step At 0.1. This results in De = 4 trajectories per value of p. yielding total of [D_| 8 32 training
II practice, a fixed point with complex eigenvalues dominated by imaginary parts of opposite signs is regarded as center-like.
18

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Canonical momentum (p)
[6]
Contesti
0.50
0.50
0.50
0.50
0.25
0.25
0.25
0.5
0.00
0.00
0.00
-0.50
0.50
-0.50
-0.50
-0.75
0075
-0.75
-0.75
-1.0
0.0
-0.5
-0.5
Canonical position (q)
Canonical position (q)
Canonical positin (4)
Canonical position (q)
Figure 14. Predicted phase portraits of the (a) vanilla NODE and (b) context-informed NODE models for the right well adaptation scenario.
Predicted phase portraits of the (c) vanilla NODE and (d) context-informed NODE models for the orbit adaptation scenario.
0.6
Canonical momentum (9)
1.3
0.4
1.0
0.2
0.5
0.5
0.0
0.0
-0.2
-05
-0.5
-0.4
Ground Tuth Corted-nfommdt zwroiho:
Ground tutt
Ground Irmm
Ground tuth
-1.0
Contertnfomes
Contest-- thit
Commat infomied
115
-0.5
0.0
0.5
-2
-2
2
Canonical position
Canonical position (4)
Canonical position (4)
Canonical position (9)
Figure Phase portraits constructed in zero-shot manner by linearly extrapolating context vector & based on pre- bifurcation training data, corresponding to different values of Ne: (a) μe 0.1, (b) μe = 1.0, (c) = 1.5, and (d) μr 2.0.
trajectories. It is important to note that all training data consists solely of single-orbit trajectories from the pre-bifurcation symmetric regime. The phase portraits of the training trajectories are visualized in Figure 13.
Architecture. We basically followed the settings outlined the original CoDA paper (Kirchmeyer et al., 2022). We employed 4-layer multi-layer perceptrons (MLPs) with hidden layers of width 64 and swish activation functions (Ra- machandran et al., 2017). We used the fourth-order Runge-Kutte Kutta (RK) method as an ODE solver.
Training details. We basically followed the settings outlined in the original CoDA paper (Kirchmeyer et al., 2022). We trained the context-informed NODEs using the Adam optimizer (Kingma & Ba, 2015) with default settings for 50,000 epochs, using learning rate of 10 and full-batch size. To improve training stability, employed exponential scheduled sampling for teacher forcing (Lamb et al., 2016), starting with an initial probability of 0.99 and applying decay rate of 0.99 every 30 epochs. Ag 10-1 and An 10 were used for the sparsity regularizer R(W, E.) in (2).
For the one-shot adaptation, We adapted the context-informed NODEs using the Adam optimizer with default settings for 1,000 epochs, using a learning rate of 10-3 and full-batch size. We employed exponential scheduled sampling for teacher forcing, starting with an initial probability of 0.95 and applying decay rate of 0.95 every 30 epochs. In comparison, the vanilla counterpart trained for 3,000 epochs under the same settings.
One-shot adaptation. After training the model with biffuuuuuation data, adapted it using single trajectory for 0.5, representing post-bifurcation data. We considered the following four broken symmetry scenarios, namely adaptations using (i) trajectory confined to the left well, (ii) a trajectory confined to the right well, and (iii) trajectory outside the separatrix, traversing the outer orbit. These scenarios limit the model's exposure to the global structure of the phase space during adaptation, reflecting the spontaneous symmetry breaking observed in real-world situations. For comparison, we also trained vanilla neural ODE models (1) under each scenario. Figure 14 compares the phase portraits of vanilla NODEs and context-informed NODEs: (a-b) for the right well adaptation scenario and (c-d) for the outer orbit adaptation scenario (see Figure 2 for the left well scenario). The context-informed NODEs successfully reconstruct the entire phase topology across all scenarios.
Zero-shot exploration. We extrapolated the context values for post-bifuree costbbifucattin cases with = 0.1, 1.0, p. = 1.5, and Me = 2.0 using a linear regression model fitted to the relationship between pre-bifurcation & and (4p (see Figure (a)). Figure 15 illustrates the resultant phase portraits obtained in zero-shot manner, showing that this method effectively reproduces the ground truth phase portraits with reasonable accuracy.
19

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Ground truth -2.0)
Ground -1.75)
Ground
Ground =11225)
1.0
1.0
1.0
1.0
Canonical momentum (p)
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
-0.5
-0.5
-0.5
-0.5
-1.0
-1.0
-1.0
-1.0
-0.5
Canonical position (q) fruth
Canonical posilion (4)
Canonical position (q)
Canonical position (q)
(u -0.75)
(μ -0.25)
1.0
1.0
1.0
1.0
Canonical momentum (p)
0.5
0.5
0.5
0.5
0.0
0.3
0.0
0.0
005
-0.5
-0.5
-0.5
-1.0
-1.0
-1.0
-1.0
-0.5
-0.5
Canonical position (q)
Canonical position (4)
Canonical position
Canonical position (q)
Figure 16. Examples of training trajectories from the experiment descrihed in Section 3.2.
Ground truth (u -0.35)
Ground truth -0)25)
Ground truth (u (0.10)
Ground truth -0.05)
1.0
1.0
Canonical momentum (8)
0.5
0.5
0.0
0.0
0.5
-0.5
-0.5
-05
-1.0
-1.0
1.0
-0.5
0.0
0.5
1.0
-1.0
-0.5
0.0
0.5
-1.0
-0.5
0.5
1.0
110
-0.5
0.0
0.5
Canonical position (q)
Canonical position (q)
Canonical position (g)
Canonical position (q)
Ground 0.05)
Ground truth
Ground truth 0.25)
Ground truth 0.35)
1.0
1.0
1.0
1.0
Canonical momentum (P)
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
-0.5
00.5
-0.5
-0.5
-1.0
-1.0
-1.0
-..0
Canonical posilion (q)
Canonical posilion (q)
Canonical position (q)
Canonical pasition (q)
Figure 17. Examples of training trajectories from the experiment described Section 4.3.
E. Experiment Details: Section 3.2
Data preparation. We trained the context- context-informed NODES (2) on the Hamiltonian system described Remark B.3 under the pre-bifureation regime with restricted training domain. Specifically, we randomly sampled four initial conditions from the uniform distribution (())),((0)
([,50.5, 0.5) [-0.5, 0.5)) for each of the eight parameter values μº= (-2.0, -1.75, -1.5- 11255, 1.0, -0775, -0.5, -0.25). For each sampled initial condition, simulated the dynamics with time horizon 7 = 2.0 and time step Al = 0.1. This results in Del = 4 trajectories per value of n. yielding total of |D_| x 8 32 training trajectories. It is important to note that all training data consists solely of single-orbit trajectories from the pre-bifurcation symmetric regime. The phase portraits of the training trajectories are visualized in Figure 16.
Architecture. We basically followed the settings outlined the original CoDA paper (Kirchmeyer et al., 2022). We employed 4-yayer MLPs with hidden layers width 64 and swish activation functions. We used the fourth-order RK method as an ODE solver.
20

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Ground truth -2.0)
Ground buth -1.75)
Ground
Ground =11.25)
1.0
1.0
1.0
1.0
Canonical momentum (p)
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
-0.5
-0.5
-0.5
-0.5
-1.0
-1.0
-1.0
-1.0
0.5 Canonical posilion (4)
-0.5
Canonical position (q) fruth
Canonical position (q)
Canonical position (q)
truth (u -0.75)
-0.5)
(j -0.25)
1.0
1.0
1.0
1.0
Canonical momenium (p)
0.5
0.5
0.5
0.5
0.0
0.3
0.0
0.0
-05
-0.5
-0.5
-0.5
-1.0
-1.0
-1.0
-1.0
-0.5
-0.5
Canonical position (4)
Canonical position (4)
Canonical position (4)
Canonical position (q)
Figure 18. Examples of training trajectories from the experiment described in Appendix G.
Training details. We basically followed the settings outlined in the original CoDA paper (Kirchmeyer et al., 2022). We trained the context-informed NODES using the Adam optimizer with default settings for 10,000 epochs, using a learning rate 10-3 3 and full-batch size. To improve training stability, we employed exponential scheduled sampling for teacher forcing, starting with an initial probability of 0.99 and applying decay rate of 0.99 every 30 epochs. Aς 10 and An 10-6 were used for the sparsity regularizer R(W, Ee) in (2).
F. Experiment Details: Section 4.3
Data preparation. We trained the context-informed NODES (2) and their Hamiltonian version (Greydanus et al., 2019) on the Hamiltonian system described in Remark B.2 considering both pre- and post-bifurcation regimes within restricted training domain. Specifically, we randomly sampled four initial conditions from the uni- form distribution (g(0),p(0))
u([-0.3,0. 0.3] x (-3) for each of the eight parameter values μe {-0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35]. For each sampled initial condition, with time horizon T 2.0 and time step Al = 0.1. This results in Dal = 4 trajectories per value of n yielding total of D.l x 8 = 32 training trajectories. The phase portraits of the training trajectories are visualized in Figure 17.
simulated the dynamics
Architecture. Note that for the Hamiltonian context-informed NODEs, a neural scalar Hamiltonian function H. is parameterized in the context-informed sense:
x (4.1) f(x, 0. WC.) (x,H(x;; WEe), (0,))--,,Hen 0x WW&))
Then, this model is trained using (2) in the same manner as described in Section 2. We employed 4-yayer MLPs with hidden layers of width 64 and tanh activation functions. We used the fourth-order RK method as an ODE solver.
Training details. We basically followed the settings outlined in the original CoDA paper (Kirchmeyer et al., 2022). We trained the context-informed NODES using Adam optimizer with default settings for 10,000 epochs, using a learning rate of 10 and full-batch size. To improve training stability, we employed exponential scheduled sampling for teacher forcing, starting with an initial probability of 0.99 and applying decay rate of 0.99 every 30 epochs. xg 10-4 and Ag 1 were used for the sparsity regularizer R(W, Ee) in (2).
G. Linear System Identification with the Standard Pre-Bifurcation Training Setting
Data preparation. We trained the context-informen NODEs (2) and their Hamiltonian version (Greydanus et al., 2019) on the Hamiltonian system described Remark B.2 under the pre-bifurcation regime with restricted training domain.
21

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Gomesponding parameter
Conssgonding paramater 0.0
-06
0.0
0.5
1.0
1.5
2.0
-0.5
0.5
1.0
2.0
0.4
0.2
Location of fiod paints (q (")
1.0
Canonical momentum (p)
(117)
(node) (8T)
Suddic (nssel)
-aadli (GT)
Sodall
Sodale (OT)
0.2
0.1
0.5
0.0
0.0
0.0
decee
-0.2
-0.5
2
-04
-0.2
0.00
0.05
0.15 0.20
0.25
Canonical pasition (9)
Canonical position (4)
Context (-5)
Context (5)
Figure 19. Phase portraits of (a) free-form and (b) Hamiltonian models trained with (á. = (p, jg)), and constructed in a zero-shot
manner by lincarly extrapolating &c. correspond to uee 0.1. Bifurcation diagrams of (c) free-form and (d) Hamiltonian models.
Specifically, randomly sampled four initial conditions from the uniform distribution (())), 4[[-033,03 x [-0.3, ). for each of the eight parameter values "²" {22.00- -1.75, -1.5, 11.25, -1.0, 00.75 -0.5, -0.25). For each sampled initial condition, we simulated the dynamics with a time horizon T = 2.0 and a time step A = 0.1. This results in [De] = trajectories per value of tr yielding a total of [De| x 8 = 32 training trajectories. It is important to note that all training data consists solely of single-orbit trajectories from the pre-bifurcation symmetric regime. The phase portraits of the training trajectories are visualized in Figure 18.
We employed 4-layer MLPs with hidden layers of width 64 and tanh activation functions, following Appendix F. We used the fourth-order RK method as an ODE solver.
Architecture.
Training details. We basically followed the settings outlined in the original CoDA paper (Kirchmeyer et al., 2022). We trained the context-informed NODEs using the Adam optimizer with default settings for 10,000 epochs, using learning rate of 10 -3 and full-batch size. To improve training stability, we employed exponential scheduled sampling for teacher forcing, starting with an initial probability of 0.99 and applying decay rate of 0.99 every 30 epochs. Xg 10-4 and An 10-6 were used for the sparsity regularizer R(W, &e) in (2).
Hallucinated broken symmetry. Figure 19 (a-b) illustrates the phase portraits explored by context-informer models trained on the linear system under the setting described in Appendix G. As shown, both the free-form and Hamiltonian NODEs misinterpret the bifurcation, incorrectly identifying sition, resulting in spurious double-well structure. The deulle-well structure in the free-form model collapses with moderately large &, reverting to a typical saddle point (Figure 19 (c)). The Hamiltonian model retains the misidentified double-well structure persistently (Figure 19 (d)). In addition, interestingly. the Hamiltonian model locally exhibits a triple-well structure at ę Eerit (that corresponds to He ~ O, see Figure 20), though collapses very rapidly. It can be understood in similar way that, if the model locally generates two saddle points incorrectly, with their Poincaré indices summing to 2, while (0, 0) remains center with index of +1 and the outer orbit structure is preserved, the model must generate two additional centers with total index of +2 to offset the discrepancy, ensuring that the overall index sums to - + 1 + 2 = +1 on the outer orbit. However, as the model correctly transitions the center at (0, 0) into saddle point, this structure quickly collapses, giving way to a (still spurious) double-well structure to maintain a total index of +1.
Canonical momentum (p)
0.10
0.05
symmetry-breakee
tran-
0.00
-0.05
Hamilonian, {= 0.0853
-0.10
-0.2
0.0
0.2
Canonical position (q)
Figure 20. The triple- well phase por- trait of the Hamiltonian NODE model. correspond to μie 0.01.
H. Experiment Details: Section 5.1
Data preparation. We trained the context-informet NODEs (2) the cusp bifurcating system under the pre- rre-bifurcation
regime exclusively. Specifically, we randomly sampled four initial conditions from the uniform distribution (g(0), p(0) u([-2.0, 2.0] x [-2.0, 2.0] for each of the 16 parameter vectors (ntt !² E {-2.0, -1.5, -1.0, -0.5)2. For each sampled
initial condition, we simulated the dynamics with a time horizon T = 2.0 and a time step At = 0.1. This results in D. = 4
trajectories per vector of (nt², 1/"), yielding a total of Del x 16 = 64 training trajectories. is important to note that all training data consists solely of monostable trajectories from the preubifurration regime.
22

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
(u.v)
(u.v)= (-2.0.00.))
[.V)=(2200088
(U.)= (.2.0.2.))
Momentum (p)
Ground buuh
Ground trrn Confosclifommed
Ground truth
Groune Futh
Conterii
informed.
Connent incrmedd
-2
-2
Position (q) (U.) (-0.5.-0.5)
(q)
Pasition (q) (-0.5,2.))
(-0.5.22.)
(L.) (-0.5.055)
Momentum (p)
Ground muth
Ground Inth Conteseentommedd zero
Ground Inch Commest-ntremmn shro
Groune Puth Commun-nttoeedd
Contersiromre
-2
-2
Position
Posilion (q)
(q)
Position (q) (0.5.2.0)
(u.) (0.5.-2.))
(U.) (0.5.-0.5)
(AA) (0.5.0.5)
Momentum (P)
-1-
Ground but
Ground Inis
Grouns buth Corment-eoormmd zama-sool
Drsune butt
Cintosderaomma, zorb-shol
Corrast-nfomedd cotommedddrrrtttte zerovation
-2
Position
Posttion (q) (u.v) (2.0.0.5)
Position (q)
Position (q) (μv) (2.0.2.0)
(U.) (2.0.-20)
(LV) (2.0.0.5)
Momentum (p)
D-
Ground suth
Cround tur?
Bround buth
Ground Futh
Contucinfomrss 28ro-thise
Collestiformed zcre-shit
Comiat- zero-shet
Contert- zero-shol
-2
-2
Position
Position
Position (q)
Position
Figure 21. Phase portraits constructed in zero-shot manner by linearly extrapolating the context vector & based pre-bifurcation training data, corresponding to different values of (pe.v.).
Architecture. We basically followed the settings outlined in the original CoDA paper (Kirchmeyer al., 2022). We employed 4-layer MLPs with hidden layers of width 64 and swish activation functions. We used the fourth-order RK method as an ODE solver.
Training details. We basically followed the settings outlined in the original CoDA paper (Kirchmeyer et al., 2022). We trained the context-informed NODEs using the Adam optimizer with default settings for 50,000 epochs, using learning rate 10-3 and a full-batch size. To improve training stability, we employed exponential scheduled sampling for teacher forcing, starting with an initial probability of 0.99 and applying a decay rate of 0.99 every 30 epochs. Ag = 10-4 and Ao 10-6 were used for sparsity regularizer R(W, {x) in (2).
Evaluation details. After training the model, we constructed mesh grid over (p test (etest) E [-2.0, 2.0 with intervals of Apttr 0.1 and Agtestt 0.1, resulting total 41 x 41 1681 parameter combinations, For each parameter vector (ptest 1yest), the vector field of the learned context-informed NODE model was obtained by inputting the corresponding context vector glest, which was estimated using linear regression model fitted on the relationship between the training parameters (jobr oft) and their associated context vectors ttr Then, the bifurcation surface shown in Figure 10 (a) was generated by identifying the fixed points of the modeled vector field at each gtest using the neural-adjoint method described in Appendix C. Figure 10 (b) directly compares the ground truth and the modeled vector fields for some selected parameter vectors near the catastrophic transition. In addition, an extended comparison over a broader range parameters (He, ve) is presented in Figure 21.
23

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
I. Poincaré-Hopf Regularization
The Poincars-Hopf regularization is formulated to minimize the following topological discrepancy by guiding the model's Poincaré indices toward the desired values:
RPH (0. W,E)
Σ Er (Txxp(( [P:Ind(f(:Oe: (D 1)),
where D--II-) represents distance function, X; is the i-th desired index value, & is the context vector, and pr is a distribution over the test contours r. Because the appropriate test contour depends on the desired index, pr is conditioned on X (e.g.. +1 for global test contour or -1 for a local test contour the case of Kouble-well potentials). The model's index [nd r) computed by numerically evaluating the two-dimensional contour integral described in Remark 4.1:
M
Sq. + Sq..
Ind(f,r)
2T
|| +
where M is the number of segments used to discretize I and S: (Sa.s fn.)² denotes the vector field at the i-th segment of l.. The terms A f = fi(p).(++ (1moMM - fu(p), represent the finite differences of the field components. The small constant c = 10-4 is added to avoid division by zero when [ - 0.
In the experiment described in Section 5.2 (and detailed Appendix J), we used a simplified version of the regularizer by fixing e E Eur and setting to predetermined global and local test contours as follows:
RPH W, Ex)= ΛPH (Rin(Oes W.) + (ppn(0e, W.
where (0p Ee) and Rfp W. respectively measure the global and local index mismatches, defined as R²H
Σ [Ind([(;:e++ W&)), (H)) 112.
|Eeel
eEEer
Rpm
[nd[f(: Wte), Fin) + 11/2 + Kela > 0H
1
[[nd([(; + WE.). 11|2. e|a>0
Σ
Hela < OH elaco
Here, PH represents a global contour centered at (0, 0), with a major axis length of r. = 4.0 and minor axis length of Tig = 0.5. Meanwhile, PHE defines a local contour as a circle with a radius of 0.5. Note that the desired index for the global regularization RpH is fixed at +1, following the Poincaré-Hopl theorem. For the local regularization R2H. the desired index is set to +1 for a > 0 (i.e., the single- and triple-well cases) and to -1 for a < 0 (i.e., the double-well dynamics), which requires more specific prior knowledge. In Appendix K, we compare the full regularizer defined in (7), i.e., RpH APn[R] + RpH)), with its global-only variant RPH RpnRiaH The hyperparameter λPH is set to 10
J. Experiment Details: Section 5.2
Data preparation. We trained the conventional context-informed NODEs and their topologically regularized counterpart (LK) system (4). For training, we randomly sampled 4 initial conditions from the uniform distribution (v²(0). o' (0)) 2([-2.0, 2.0] [-2.0, 2.0) for each combination of o (-0.4, 0.2, 0.2, 0.4] and f E {-0.4, -0.2, 0.2, 0.4), resulting in 16 different parameter vectors. To ensure the physical relevance of (4), y 0.05 and = 0.5 were fixed. For each sampled initial condition, simulated the dynamics over T = 2.0 with Al = 0.1. This results in Dre
with (7) on the
trajectories per value of n. yielding a total of [De] x 16 64 training trajectories.
Architecture. We basically followed the settings outlined in the original CoDA paper (Kirchmeyer et al., 2022). We employed 4-layer MLPs with hidden layers of width 64 and swish activation functions. We used the fourth-order RK method as an ODE solver.
Training details. We basically followed the settings outlined in the original CoDA paper (Kirchmeyer et al., 2022). We trained the context-informed NODES using the Adam optimizer with default settings for 50,000 epochs, using learning rate of 10-3 33 and a full-batch size. To improve training stability, we employed exponential scheduled sampling for teacher forcing, starting with an initial probability of 0.99 and applying a decay rate of 0.99 every 30 epochs. Ag = 10-4 and Aa = 10-6 were used for the sparsity regularizer R(W, {e) in (2). The regularizer (7) described in Appendix is added to (2) for the topologically regularized models. We repeated the experiment five times using different random initializations.
24

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Parameler B
02-
Deu
Bingle-well
0.2
Double-well
Single-well
0.75
Parameler p
0.2
Double weell
Single-well
0.2
weul
wenglo
0.7
0.0
0.0
0.0
0.0
0.45
0.45
Double-well
Trinfewee
Double-well
Trill-wwll
Doubleweel
Double-well
Triplewwell
-0.4
-0.4
-0.4
-04
0.00
0.00
0.4
02
-0.2
Parameter
Parameter
Parameter
Parameter a
0.2
0.5
Parameler B
02-
Double well
Single-wer well
0.2
Double well
Single
0.75
Paramelar p
0.2
Double well
Single-well
0.
Double welllewwlll
wingie-well
0.75
0.0
0.0
0.45
0.45
-0.2
Double-well
Triple
-0.2
Dauble-well
Triple-well
0.30
-0.2
Double-well
Triple
-02
Double-well
Triple-well
0,30
-04-
-0.4
-0.4
-0.4
0.00
0.00
0.4
Parameter
Parameter
Peremeter
Parameter
(ii)
0.7
(iv)
1.0
Figure 22. Contour plots of the consistency scores for (a) the conventional model and (b) the fully regularized model using (7).
Evaluation details. After training the models, we evaluated the models' ability to capture the long-term behavior and phase topology of the dynamical system by computing the L2 distance between the limit sets of the ground truth and those generated by each model. We constructed a mesh grid over (a,est, g [-0.5, 0.5]2 with intervals of An 0.05 and A, = 0.05, resulting in total of 21 x 21 441 parameter combinations. For each parameter vector (atest , 3test), we randomly sampled N = 32 initial conditions from the uniform distribution (v²)),, v (0)) ~ u([-2.0, 2.0] x 2-2.0, 2.0)) and simulated long-term trajectories with T' = 100.0 and At = 0.1. Context-informee NODE models and their regularized counterparts were simulated under the same and At. These simulations utilized 441 context vectors, constructed using linear regression fitted on the relationships between the training parameters (a ( and the corresponding context vectors Se The limit set for each dynamics model is defined as the converged attractors (v²", +(²)) xx ~ x(T"). We empirically found that T = 100.0 is sufficient to approximate the limit sets accurately. After obtaining x(T) for both the ground truth and the tested NODE models, the topological consistency score of e-th test parameter vector (atr, is calculated as 8c N/ ΣN [[]x{ (T) (T) o), where a threshold distance, and x" (T') and x (T') represent the state vectors T of the ground truth dynamics and the NODE models, respectively, starting from the i-th test initial condition for the e-th parameter vector.
Figure 22 presents contour plots of the topological consis- tency scores as the threshold a varies. As shown, the con- sistency scores of conventional context-informed NODE models show limited improvement despite increasing σ. particularly in the triple-well regions (a 0, 3 < 0). This suggests that the vanilla context-informed NODE struggles capture the correct topology, leading to di- verging behavior. On the other hand, the topologically regularized model demonstrates enhanced long-term pre- dictability. We also report the Mean Absolute Percentage Error (MAPE), which provides more reliable basis for performance comparison across different parameters than MSE (Kirchmeyer et al., 2022), shown Figure 23. This further validates the effectiveness of the proposed Poincaré- regularization.
Parameter B
Bouble-well
Single-well
Double-well
Single-well
46
36
0.0
0.0
27
-0.2
Boublewell
Teppee (eplewwell
-0.2
Double well Triple-wlll
18
-0.4
-02
-0.4
0.2
Parameter
Parameter a
Figure 23, Contour plots of the MAPEs for (a) the conventional model and (b) the fully regularized model using (7).
K. Ablation Study of Poincaré-Hopf Regularization
Topological regularization without a local prior. The proposed regularization (7) requires prior knowledge at both global (R&h)) and local (Reu) levels. Local regularization typically demands more detailed, system-specific information. For example, computing RpH in (7) necessitates prior knowledge about the transition behavior of the local index from
25

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Parameler B
02-
Deublowell
Bingle-well
0.2
wouble-well
Single-well
0.75
Parameler p
0.2
Double weell
Single-well
0.2
weullewelll
Single-well
0.75
0.0
0.0
0.0
0.0
0.45
0.45
Double-well
Trinfewee
Double-well
Tripl-well
-0.2
Doubleweell
Double-well
Triple-well
-0.4
-0.4
-0.4
-04
0.4
0.0 Parameter a
0.00
02
-0.2
0.00
Parameter
Parameter
Parameter a
0.2
0.5
Parameler B
02-
Double well
Single-wer well
0.2
Double-well
Single
0.75
Paramelar p
0.2
Double well
Single-well
0.
Double-well
welgie
0.75
0.0
0.0
0.0
0.45
0.45
-0.2
Double-well
Triplee
-0.2
Dauble-well
Triple-well
0.30
-0.2
Double-well
Triple
-02
Dauble-well
Triple-well
0.30
-04-
-0.4
-0.4
-0.4
0.00
0.00
0.4
Parameter
Parameter
Peremeter
Parameter
(ii)
0.7
(iv)
1.0
Figure 24. Contour plots of the mean consistency scores for (a) the conventional context-informed NODE and (b) the regularized version with only the global constraint.
as a function of a. In contrast, applying only global constraints requires significantly less prior information. As stated in Theorem 4.1, the global regularization term Rih relies only on coarse- -rannned knowledge about the phase manifold where the ODE is defined; for instance, in the L.K experiment, Riz only requires that the global index constraint equals + +1, regardless of parameters. Despite its simplicity, this type of global information can be especially valuable when detailed prior knowledge about the target system is limited.
To evaluate the effectiveness of topological regularization under such minimal assumptions, we revisited the LK experiment discussed in Appendix using only the global term, i.c. RPH Ee) λри RipH W, E.). All other experimental settings remained unchanged. Figure 24 compares the mean consistency scores between the vanilla model and the globally regularized model across different values of a, following the evaluation protocol described earlier. The results show that the globally reg- ularized model (Figure (b)) significantly outperforms the vanilla model (Figure 24 (a), identical to Figure 22 (a)), even when relying solely on global information. While its performance is slightly lower than that of the fully regularized model with both local and global constraints (Figure 22 (b)), this reflects trade-off between the cost of incorporating prior knowledge and achieving optimal performance. We also report the MAPE profiles in Figure 25, which lead to similar results. These findings suggest that imposing global topological constraints via the Poincaré-Hopf theorem provides a cost-effective and robust form of topology aware regularization, particularly in settings with limited prior knowledge.
(b)
0.4
Parameter B
0.2
Bouble-well
Single-wor well
0.2
Double well
Single-wol
46
36
0.0
0.0
27
-0.2
Toubledwell
Telpleswell
-022
Double-well
Triple-well
-0.2
0.4
-04
0.2
Parameter
Parameter
Figure 25. Contour plots of the MAPEs for (a) the conventional model and (b) the regularized version with only the global constraint.
Training Dynamics Analysis. The proposed Poincaré- regularization described in Appendix I based on a finite difference method, resulting in computational overhead proportional to the number of discretized points M. In all our experiments, we discretized the test contours using M = 128 points, which led moderate increase in runtime under identical computational environments (vanilla: 52.9 5.5 ms/epoch; global-only regularized: 65.7 + 5.1 ms/epoch; fully regularized: 77.6 + 4.7 ms/epoch). Despite this additional cost, the regularized model demonstrates faster convergence in practice and significantly improves test accuracy, thanks to the guidance provided by topological regularization.
Figure 26 compares the vanilla, globally regularized, and fully regularized models in terms of the following metrics, all plotted over normalized wall-clock time: (a) the global regularization value Ry (b) the local regularization value Rpgss (c) the test MSE on a linear time axis, and (d) the test MSE on a logarithmic time axis. As defined in (7), the regularization values directly reflect index mismatches between the model and the ground truth, making them useful for evaluating whether
26

---

Context-Informed Neural ODEs Unexpectedly Identify Broken Symmetries
Global regulortzttin
10"
Loosl regilaitzalio
10"
10
Varkii Regularions
10"
10
Test MSE
Test MSE
10"
-0"
10"
10
10"
Vhrdli Regulartzdd
Variis
arbodd
Regubaroodd (giatal ont)
(1)
(b)
(c)
(UD
0.8 1,0
0.0
0.2
1.0
0.0
0.2 0.4
1.0
10
10
10"
10"
10"
Normalized wall-clock 6me
Normalized wall-clock time
Normalized woll--ook time
Nommlized wall-dock time
Figure 26. Comparison between the vanilla and regularized models in terms of (a) the global regularization value RpH, (b) the local
regularization value Ryn, (c) the test MSE on a linear time axis, and (d) the test MSE on logarithmic time axis. In (a) and (b), the black dashed line indicates the discretization error limit corresponding to the ground truth dynamies.
the model adheres to the ground truth topology. Note that while all models compute both Ri and RipHs only the fully regularized model is explicitly trained minimize both terms. The vanilla model does not use any form of regularization, whereas the globally regularized model is constrained solely by the global term Ri The test MSE was computed in the same way as the training MSE loss, but using 32 trajectories per parameter that were different from those used during training. All of these loss profiles are from the previously discussed LK experiment.
As shown in Figure 26 (a), both the globally and fully regularized models maintain low and stable global index mismatch, converging to x 10-3 with minimal variance, as they jointly minimize data loss and the global topological loss. is worth noting that even the ground truth vector field incurs a nonzero global topological loss of 4 x 10 due to discretization error introduced by the finite approximation of the contour integral. In contrast, the vanilla model, which lacks any topological constraints, exhibits significantly higher global topological eπror exceeding 10 10-?, along with substantial variance. These results suggest that the vanilla model fails to reliably capture the correct global topological structure, highlighting the importance of incorporating topological regularization.
In Figure 26 (b), the fully regularized model, which explicitly minimizes the local regularization term, exhibits a sharp decline in this value and maintains it at a consistently low level of 2 x 10-4, which corresponds to the discretization limit of the ground truth. Meanwhile, the vanilla model shows gradual decrease, as it implicitly captures certain topological features from the data, but its learning efficiency significantly lower than that of the regularized model. Interestingly, the global-only model, despite lacking the local regularizer, also achieves noticcable and consistent reduction in the local index discrepancy, though less effectively than the fully regularized model. This phenomenon can be attributed to the discussion in Section 4: the global topological constraint- the total Poincaré index of +1- 1a1-naturallyy encourages the model to learn the correct local topological structure that governs bifurcation, in accordance with the Poincaré-Hopf theorem.
This difference in topological learning leads to faster convergence and better predictive accuracy in the regularized cases, as illustrated in Figure 26 (c-d). Between the globally and fully regularized models, the latter demonstrates superior performance, consistent with the earlier discussion of the trade-off between incorporating prior knowledge and achieving optimal performance. Overall, these findings underscore the effectiveness and efficiency of topology-awan regularization.
27


Reviews:
[{"review_id": "Ppdgz8Vpll", "reviewer": "Reviewer_RFdK", "summary": "This paper explores the use of context-informed NODEs to identify symmetry-breaking bifurcations in dynamical systems without relying on physics-based training data. The authors demonstrate that NODEs trained solely on symmetric (pre-bifurcation) data can predict post-bifurcation behaviors in a \"zero-shot\" manner, meaning they do so without prior experience with such data.\n\nThe key to this capability lies in the NODEs' ability to utilize topological invariants, specifically the Poincaré index, which helps them grasp the underlying structure of the data despite significant changes from bifurcations. The authors also introduce a novel topological regularizer inspired by the Poincaré-Hopf theorem to enhance NODE performance.\n\nThe paper challenges the assumption that data-driven models struggle with bifurcations and presents a promising approach that demonstrates how these models can extract meaningful dynamics from limited training data, contributing to their robustness and applicability in real-world complex systems.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "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": "We sincerely appreciate the reviewer’s thoughtful and constructive feedback. Below, we address each comment carefully. Additional experimental results are available in the **README** at https://github.com/anonymous-account123/icml2025-7637 and will be thoroughly incorporated into our revised paper.\n***\n**Q1. Do NODEs genuinely utilize topology? Is the hallucinated broken symmetry general?**\n\nYou highlight a critical point about distinguishing whether context-informed NODEs genuinely leverage topological constraints or merely approximate functions that happen to exhibit topological properties. Indeed, the double-well scenario poses a dilemma: better approximations naturally reflect correct symmetry-breaking behavior. Our hallucinated bifurcation scenario in Section 4.3, grounded explicitly in Proposition 4.1, is intentionally designed to address this challenge. If NODEs were merely approximating a functional form, a simple center-to-saddle bifurcation would emerge. Instead, the model exhibits an incorrect symmetry-breaking bifurcation, confirming the model genuinely utilizes (in this scenario, misuses) topological constraints, precisely as Proposition 4.1 predicts. Hence, our findings demonstrate that NODEs indeed exploit topological invariants. In addition, the hallucinated symmetry-breaking can arise whenever Proposition 4.1's conditions hold. We will highlight this more clearly in the revised manuscript.\n***\n**Q2. Robustness of bifurcation identification under noisy conditions.**\n\nWe acknowledge that approximation errors may increase under noisy perturbations or with limited training samples, as you pointed out. To address this, we conducted further experiments under these realistic conditions. Despite the increased variance, the model consistently identified the symmetry-breaking bifurcation, demonstrating its robustness. The results are illustrated in **Figure 1(b-c)** at the provided link, which will be included in the revised manuscript.\n***\n**Q3. Non-Hamiltonian and codimension-higher bifurcations.**\n\nFollowing your suggestion, we tested the context-informed NODEs on a non-Hamiltonian, codimension-2 cusp bifurcating system: $(\\dot{x}, \\dot{y}) = (y, b + a x - x^3 – d y)$, with variable parameters $(a,b)$ and fixed $d=0.5$. This system serves as a canonical model for capturing catastrophic transitions and hysteresis phenomena. Despite its dissipative, non-Hamiltonian nature, its bounded dynamics imply a preserved total Poincaré index of +1. Training was conducted exclusively on simple, monostable pre-bifurcation conditions within the parameter range $(a,b) \\in [-2.0, -1.5, -1.0, -0.5]^2$, yet the model successfully identified the cusp bifurcation surface, confirming broader applicability (**Figure 3** at the provided link). This result will be included in our revised paper.\n***\n**Q4. Generalization of Proposition 4.1.**\n\nOur analysis primarily targets closed-orbit systems (e.g., Hamiltonian flows), but the underlying principle of Proposition 4.1—the topological invariance of the total index from the Poincaré–Hopf theorem—has broader applicability. Indeed, Proposition 4.1 follows directly from Theorem 4.1 under closed orbit assumptions. Yet, the essence of the theorem lies in global topological invariants, specifically the Euler characteristic of the phase space $\\mathcal{M}$ (e.g., $\\chi(S^{2n}) = 2$, $\\chi(S^{2n+1}) = 0$, $\\chi(T^{n}) = 0$, ...), allowing its application even when the associated ODE defined on $\\mathcal{M}$ does not exhibit closed orbits or conserved quantities. We will explicitly address this generalization potential in our revision.\n***\n**Q5. Experimental designs.**\n\nWe have provided a detailed description and implementation of our proposed topological regularization approach in Appendix Section G. Additionally, Figure 17 in Appendix illustrates consistency scores across various threshold values ($\\sigma$) ranging from 0.2 to 1.0, for the results of Section 5. The regularized model consistently outperforms vanilla models, maintaining scores close to 1.0 regardless of $\\sigma$. Furthermore, Appendix Section H thoroughly summarizes our experimental setup (e.g., 16 training and 441 testing parameter combinations over an $(\\alpha,\\beta)$ grid).\n***\n**Q6. Why four training samples and one for adaptation?**\n\nOur experimental setup closely mirrors Kirchmeyer et al. (2022), who utilized four initial conditions per parameter for training on the Lotka–Volterra system—another 2D closed-orbit system analogous to ours. For adaptation, a single confined trajectory was selected intentionally to simulate an extreme symmetry-breaking scenario. This deliberately restricts the model’s exposure to the global phase space structure, effectively testing simultaneous parameter and initial condition OOD challenges.\n***\n**Q7. Minor comments.**\n\nWe will clarify the manuscript as recommended. For the clarified definition of OOD in parameters, please refer to our response to Question 3 from Reviewer rQs3.", "claims_and_evidence": "Overall, their claims are almost well-supported by the evidence presented. In particular, the authors provide a solid theoretical framework grounded in the Poincaré-Hopf theorem and demonstrate the efficacy of their approach through empirical validation with the Landau-Khalatnikov system. The results show that context-informed NODEs can identify symmetry-breaking behavior without direct training on post-bifurcation data, which testifies to the power of their model. However, there are some Unclear Claims:\n\n- The authors claim that NODEs implicitly learn topological invariants like the Poincaré index, but it is unclear if they truly \"learn\" them or if this just happens as a side effect of function approximation. The authors could make this claim stronger by showing exactly how the model learns these features.\n\n- The paper presents an example in which NODEs \"hallucinate\" broken symmetry when trained on a linear system (Sec. 4.3). But, it is unclear whether this phenomenon is a general characteristic of NODEs or if it is specific to the chosen experimental setup.\n\nThere are also some concerns that could impact the robustness of their claims more generally in more complex scenarios: \n\n- The study predominantly investigates **codimension-one** bifurcations, focusing solely on **2D** and **Hamiltonian** systems, which narrows the applicability of their findings. Although Theorem 4.1 hints at potential extensions to more complex systems or higher dimensions, it would strengthen the argument if the authors elaborated on potential strategies for tackling **codimension-higher bifurcations**, the challenges they might face in 3D or higher-dimensional systems, and the adaptations required in their approach for non-Hamiltonian systems. It would be helpful to provide insights into how these generalizations might occur and under what specific conditions it would hold true.", "methods_and_evaluation_criteria": "**Proposed Methods**: First off, they have nicely discussed the OOD generalization problem by categorizing it into two aspects: OOD in initial conditions and OOD in model parameters, to better clarify the two different challenges for learning dynamics. The use of context-informed NODEs and Poincaré-Hopf regularization makes sense for studying bifurcations and symmetry-breaking in dynamical systems. \n\n**Evaluation Approach**: The paper evaluates NODEs on well-defined 2D dynamical systems (e.g., Hamiltonian and Landau-Khalatnikov models), which are relevant for studying bifurcations. The experiments mostly involve Hamiltonian or nearly conservative systems and focus only on codimension-one bifurcations. Yet, real-world systems often include dissipation and forcing, which can alter bifurcation behavior and lead to codimension-higher bifurcations. So, evaluating on non-Hamiltonian, higher-dimensional, and multi-parameter bifurcations ($ \\mu \\in \\mathbb{R}^n, n>1 $) would provide a more comprehensive validation.", "theoretical_claims": "I tried to verify the correctness of the theoretical claims in the main text as well as the proofs in Appendix A.\n\n**Re the proof of Proposition 4.1**: \n- The proposition outlines a condition for local symmetry breaking in a Hamiltonian vector field via a center-to-saddle bifurcation. This condition requires a smooth bijective map that maintains vector fields close in the $C^1$ sense, potentially overlooking complexities in real-world systems where perturbations are not small/easily managed.\n\n- The proof logically follows from bifurcation theory, Implicit Function Theorem and $C^1$ closeness assumptions between the learned and true **Hamiltonian** vector fields. The proof assumes the learned vector field closely matches the true system, but in practice, Neural ODEs are approximations and may not perfectly capture the system's dynamics. If there are errors in approximation (due to limited training data, network capacity, or optimization issues), the conditions stated in the proof might not always hold, making the model’s predicted bifurcations less reliable in real-world applications.\n\n- The proof in Proposition 4.1 is specifically derived for Hamiltonian systems, where energy conservation and structured dynamics naturally enforce certain constraints. However, many real-world systems are non-Hamiltonian which can significantly alter bifurcation behavior. Since the proof relies on preserving closed orbits and the Poincaré-Hopf index, its conclusions may not directly apply to non-Hamiltonian systems where these properties don't hold. The authors should clarify how far their results, particularly Proposition 4.1, can be generalized to non-Hamiltonian systems. This includes discussing the limitations of their findings and what further theoretical or empirical validation is needed to adapt their conclusions for these systems.", "experimental_designs_or_analyses": "Re soundness/validity of experimental designs or analyses: \n\n**Topological Regularization**: They introduce a Poincaré-Hopf regularization aimed at enhancing the NODEs' learning capabilities. However, it lacks detailed information on how this regularization is specifically implemented and tuned, which is essential for assessing its effectiveness and reproducibility.\n\n**Model Evaluation**: The evaluation metric relies on checking if the model predictions converge within a distance (σ) from the true dynamics. If this distance threshold is not well-defined or justified, it could lead to misleading conclusions about the model's predictive capabilities and overall performance.\n\n**Parameter Exploration**:  A mesh grid approach is used for parameter sampling, but the methodology for how this grid is constructed and the resolution of the sampling are crucial.", "relation_to_broader_scientific_literature": "The key contributions of the paper relate to several areas, particularly regarding modeling dynamical systems and understanding bifurcations through data-driven approaches. It has implications for theoretical advancements in dynamical systems and OOD generalization challenges. Specifically, it discusses OOD in model parameters, which is another category of OOD generalization challenges, in addition to OOD in initial conditions as noted in  (G¨oring et al., 2024), which is very important.\n\n**Key Points**:\n\n- Data-Driven Dynamics: Builds on Brunton et al. (2016), demonstrating effective discovery of dynamical systems without physics-based priors.\n\n- Topological Invariants: Utilizes invariants (Brasselet et al., 2009) to enhance NODEs by linking topology and dynamics.\n\n- OOD Generalization: Challenges beliefs about model generalization across bifurcations (Ye et al., 2021), showing robust NODE performance without diverse training.\n\n- Symmetry Recovery: Contrasts with García Pérez et al. (2023), indicating NODEs' ability to model symmetry-breaking without physics-informed priors.", "essential_references_not_discussed": "-", "other_strengths_and_weaknesses": "**Other Strengths**:\n\n- Originality:  The paper introduces a fresh approach by showing how context-dependent NODEs, trained on localized data, can effectively identify symmetry-breaking bifurcations. This highlights a creative use of topological invariants in a data-driven setting.\n\n- Significance: By tackling how NODEs can extrapolate behaviors beyond their training domain, this research makes a meaningful contribution to both dynamical systems and machine learning. It could open up new possibilities for automating scientific discovery and deepens our understanding of complex systems.\n\n- Clarity: The paper is almost well-written and well-structured and explains complex ideas in a way that's easy to follow.", "other_comments_or_suggestions": "**Some minor comments**:\n\n- Page 1, second column, line 21: It is mentioned \"Formally, a dynamical system is represented by a phase space Ordinary Differential Equation (ODE)\". But, it is better to say \"Formally, many **continuous-time** dynamical system ... \" as discrete-time dynamical systems can be represented by recursive maps. Also, some continuous-time dynamical systems can be represented by PDEs, not necessarily ODEs.\n\n- Page 3, in Sect. \"OOD in model parameters\":  The definition of bifurcation parameter (starting from line 151) is for the general case $\\mu_{crit} \\in \\mathbb{R}^n$. But, the sentence \"Consequently, the OOD condition for parameters arises when a model is trained on parameters $\\mu_e^{tr} < \\mu_{crit}$ but the support of the test data is $\\mu_e^{test} > \\mu_{crit}$\"  is only valid for $\\mu_e^{tr} ,  \\mu_{crit} ,  \\mu_{crit} \\in \\mathbb{R}$ (codimension-one bifurcation). Please clarify it in the text.\n\n- Page 6, Definition 4.1: The topological degree of a map should be better clarified/explained.", "questions_for_authors": "1. In Section 3.1, could you elaborate on the rationale for randomly sampling **four** initial conditions? Furthermore, what is the justification for adapting the model using a **single** trajectory after training it with pre-bifurcation data? \n\n2. In Section 2, you introduce a Poincaré-Hopf regularization to enhance the NODEs' learning capabilities. Coulld you provide more detailed information on the specific implementation and tuning process of this regularization? How do you ensure that it effectively contributes to the learning without introducing unwanted biases?\n\n3. In Section 5, the evaluation of model performance is based on the convergence within a distance (σ) from the true dynamics. How is this distance threshold (σ) determined, and what justification or heuristics do you provide for selecting its value? Could different choices of σ significantly affect your conclusions regarding the model's predictive capabilities?\n\n4. Could you clarify how far your results, particularly Proposition 4.1, can be generalized to non-Hamiltonian systems? What specific characteristics or conditions of non-Hamiltonian systems would need to be considered for the applications of Proposition 4.1 to remain valid, and how might the presence of dissipation or external forcing impact the conclusions drawn from your study regarding symmetry-breaking behavior?\n\n5. The authors claim that NODEs implicitly learn topological invariants like the Poincaré index, but it is unclear if they truly \"learn\" them or if this just happens as a side effect of function approximation. Could elaborate on this further and make this claim stronger by showing exactly how the model learns these features?\n\nI am happy to increase my score if the authors can address my main concerns.", "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": "CMA45hZDWh", "reviewer": "Reviewer_rQs3", "summary": "The paper demonstrates that context-dependent Neural Ordinary Differential Equations can identify post-bifurcation behaviors, even when trained only on pre-bifurcation data. It then provides an interpretation for this phenomenon based on the Poincaré-Hopf theorem, and proposes a topological regularizer that mitigates the hallucination for post-bifurcation behaviors.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "We sincerely appreciate the reviewer’s thoughtful and constructive feedback. Below, we address each comment carefully. Additional experimental results are available in the **README** at https://github.com/anonymous-account123/icml2025-7637 and will be thoroughly incorporated into our revised paper.\n***\n**Q1. As acknowledged by the authors, the scope of the investigation is limited.** \n\nOur analysis primarily targets closed-orbit systems (e.g., Hamiltonian flows), but the underlying principle of Proposition 4.1—the topological invariance of the total index from the Poincaré–Hopf theorem—has broader applicability. Indeed, Proposition 4.1 follows directly from Theorem 4.1 under closed orbit assumptions. Yet, the essence of the theorem lies in global topological invariants, specifically the Euler characteristic of the phase space $\\mathcal{M}$ (e.g., $\\chi(S^{2n}) = 2$, $\\chi(S^{2n+1}) =0$,  $\\chi(T^{n}) = 0$, ...), allowing its application even when the associated ODE defined on $\\mathcal{M}$ does not exhibit closed orbits or conserved quantities. We will explicitly address this generalization potential in our revision.\n***\n**Q2. The regularization term seems to require detailed knowledge of the dynamical system beforehand, such as the global and local contour. This may limit the applicability of the regularization.**\n\nWe agree with your comment that topological regularization requires a certain level of prior knowledge about the system. In particular, the amount of required knowledge depends on whether we apply only a global constraint or both global and local ones. Local regularization typically requires more detailed, system-specific information. In contrast, using only global constraints demands relatively less prior knowledge. As stated in the Poincaré–Hopf theorem, it depends only on coarse information about the phase manifold on which the ODE is defined. For example, in a 2D closed-orbit system, the total Poincaré index is fixed at +1. Thus, as long as we know the system is non-divergent, this global constraint can often be applied without detailed knowledge of the vector field itself. In other setting—such as when the phase manifold is a sphere—the total index is determined by the Euler characteristic. For instance, when modeling vector fields on Earth (as in climate models), it is typically reasonable to assume a global index of $\\chi(S^2) = +2$.\n\nThis type of global information can be especially useful when the available training data is confined to a restricted domain of initial conditions. To assess the effectiveness of topological regularization under minimal prior assumptions, we revisited the experiment in Section 5 using only the global constraint. In this additional study, we limited the training domain to $[-1.0, 1.0]^2$ to mimic a restricted coverage scenario. **Figure 5** at the provided link illustrates that the topologically regularized model exhibits improved performance, even when relying solely on global information. We will include this result in the revised paper.\n***\n**Q3. On page 3 line 160, the left column, the definition of OOD condition for parameters in not clear for n-dimensional parameters.**\n\nFollowing your suggestion, we formally describe the definition of parameter OOD conditions with codimension-$k$ bifurcation in $n$-dimensional parameter space:\n\nLet $\\mathcal{P} \\subset \\mathbb{R}^n$ be a parameter space and let $\\mathcal{B} \\subset \\mathcal{P}$ be a bifurcation set. The complement $\\mathcal{P} \\setminus \\mathcal{B}$ admits a decomposition into $L$ disjoint connected components: $ \\mathcal{P} \\setminus \\mathcal{B} = \\bigcup_{i=1}^L \\mathcal{P}_i$, where each $\\mathcal{P}_i$ is a maximal connected subdomain. Intuitively, each $\\mathcal{P}_i$ is one qualitatively uniform parameter regime. Then, a parameter OOD condition in learning dynamics arises when the support of training distribution is $\\mathrm{supp}(p^\\mathrm{tr}_e(\\mu)) \\subseteq \\mathcal{P}_l$, but the support of the test distribution is $\\mathrm{supp}(p^\\mathrm{test}_e(\\mu)) \\subseteq \\mathcal{P}_i$ for some $i \\neq l$. Equivalently, there exists no continuous path $\\gamma: [0, 1] \\to \\mathcal{P} \\setminus \\mathcal{B}$ connecting $\\mu^\\mathrm{tr}_e$ and $\\mu^\\mathrm{test}_e$ without crossing $\\mathcal{B}$.\n\nThis generalized notion of parameter OOD parallels the concept of initial condition OOD, enabling a unified perspective on OOD in learning dynamics. Accordingly, we will revise our paper to correctly describe this more general setting beyond the one-dimensional case.", "claims_and_evidence": "The claims seem well supported by the experiments.", "methods_and_evaluation_criteria": "The evaluations make sense for the investigated problem.", "theoretical_claims": "I read the proof section in the appendix.", "experimental_designs_or_analyses": "The experiments are sound and well motivated.", "relation_to_broader_scientific_literature": "This paper studies context-dependent Neural Ordinary Differential Equations in depth on its behaviors when predicting OOD data.", "essential_references_not_discussed": "I am not aware of any essential references not being discussed.", "other_strengths_and_weaknesses": "Strengths:\n1. The paper is well-written and easy to follow.\n2. The problem is well-motivated and the arguments and experiments are coherent.\n\nWeaknesses:\n1. As acknowledged by the authors, the scope of the investigation is limited.\n2. The regularization term seems to require detailed knowledge of the dynamical system beforehand, such as the global and local contour. This may limit the applicability of the regularization.", "other_comments_or_suggestions": "1. On page 3 line 160, the left column, the definition of OOD condition for parameters in not clear for n-dimensional parameters.", "questions_for_authors": "I do not have additional questions.", "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": "J5v6JJIUPO", "reviewer": "Reviewer_XQYv", "summary": "The paper finds that context-informed Neural Ordinary Differential Equations (NODEs) can identify symmetry-breaking bifurcations in dynamical systems (DS). By leveraging topological invariants like the Poincaré index and the Poincaré-Hopf theorem,  the paper demonstrates conditions under which context-informed NODEs can out-of-domain (OOD) generalize without explicit physics-based priors. The manuscript further introduces a regularization based on these findings and empirically demonstrate its use on the Landau-Khalatnikov system.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "We sincerely appreciate the reviewer’s thoughtful and constructive feedback. Below, we address each comment carefully. Additional experimental results are available in the **README** at https://github.com/anonymous-account123/icml2025-7637 and will be thoroughly incorporated into our revised paper.\n***\n**Q1. Brenner et al.** \n\nThank you for pointing out the relevant reference, Brenner et al. (2024), which also enables both interpolation and extrapolation of model parameters via linear modulation of subject-specific features. However, as you noted, it does not address the bifurcation-driven OOD, specifically the case where training is conducted solely in the pre-bifurcation regime while the post-bifurcation regime is explored at test time. We will cite this work and include a discussion in the revised paper.\n***\n**Q2. How to detect and avoid hallucinated bifurcations.**\n\nThe hallucinated symmetry-breaking experiment in Section 4.3 is intentionally designed as a pathological case to illustrate how the model mistakenly classifies a simple center-to-saddle bifurcation as a symmetry-breaking one, as predicted by Proposition 4.1. In the linear system considered in Section 4.3, if the NODE model is merely approximating a known functional form, it should undergo a simple center-to-saddle bifurcation. However, the experimental results show that the NODE model instead undergoes a spurious symmetry-breaking bifurcation. This indicates that the model is leveraging (or in this case misusing) topological constraints. Therefore, its behavior should be interpreted using the topological arguments provided by Poincaré–Hopf theorem and Proposition 4.1.\n\nWe completely agree that identifying whether a model hallucinates bifurcations is vital. We propose a straightforward yet robust criterion: evaluating the variance in bifurcation diagrams generated by multiple independent training runs. Correctly identified bifurcations show minimal variance, indicating stable and consistent identification. Conversely, hallucinated bifurcations yield significantly higher variance due to their structural instability. **Figure 2** at the provided link illustrates this distinction clearly, highlighting how increased variance serves as an indicator of hallucinated behavior. We will include this result in the revised version, along with appropriate discussion.\n***\n**Q3. Non-Hamiltonian DS.**\n\nFollowing your suggestion, we tested the context-informed NODEs on a non-Hamiltonian, codimension-2 cusp bifurcating system: $(\\dot{x}, \\dot{y}) = (y, b + a x - x^3 – d y)$, with variable parameters $(a,b)$ and fixed $d=0.5$. This system serves as a canonical model for capturing catastrophic transitions and hysteresis phenomena, which are fundamental in science and engineering. Despite its dissipative, non-Hamiltonian nature, its bounded dynamics imply a preserved total Poincaré index of +1. Training was conducted exclusively on simple, monostable pre-bifurcation conditions within the parameter range $(a,b) \\in [-2.0, -1.5, -1.0, -0.5]^2$, yet the model successfully identified the cusp bifurcation surface, confirming broader applicability (**Figure 3** at the provided link). This result will be included in our revised paper.\n***\n**Q4. Regularization without physical priors.**\n\nWe agree with your comment that topological regularization requires a certain level of prior knowledge about the system. In particular, the amount of required knowledge depends on whether we apply only a global constraint or both global and local ones. Local regularization typically requires more detailed, system-specific information. In contrast, using only global constraints demands relatively less prior knowledge. As stated in the Poincaré–Hopf theorem, it depends only on coarse information about the phase manifold on which the ODE is defined. For example, in a 2D closed-orbit system, the total Poincaré index is fixed at +1. Thus, as long as we know the system is non-divergent, this global constraint can often be applied without detailed knowledge of the vector field itself. In other setting—such as when the phase manifold is a sphere—the total index is determined by the Euler characteristic. For instance, when modeling vector fields on Earth (as in climate models), it is typically reasonable to assume a global index of $\\chi(S^2) = +2$.\n\nThis type of global information can be especially useful when the available training data is confined to a restricted domain of initial conditions. To assess the effectiveness of topological regularization under minimal prior assumptions, we revisited the experiment in Section 5 using only the global constraint. In this additional study, we limited the training domain to $[-1.0, 1.0]^2$ to mimic a restricted coverage scenario. **Figure 5** at the provided link illustrates that the topologically regularized model exhibits improved performance, even when relying solely on global information. We will include this result in the revised paper.", "claims_and_evidence": "All claims are supported by evidence, even though some more evidence would be helpful (see **Weaknesses**).", "methods_and_evaluation_criteria": "Methods and evaluation criteria are appropriate.", "theoretical_claims": "The paper's theoretical findings are based on established DS theory, such as the Poincaré-Hopf theorem. Appropriate references thereof are provided (e.g. Strogatz 2018). I did not check the proof for Proposition 4.1.", "experimental_designs_or_analyses": "Experimental Analyses seem appropriate.", "relation_to_broader_scientific_literature": "Training on a limited window of ground-truth (GT) control parameter values using a hierarchical/meta-learning approach and extrapolating beyond this range has also been done in a previous study, albeit not across bifurcations [1]. However, similar linear relationships between context vector and GT control parameters have been observed (see Fig. 3a). Also, as mentioned in the current manuscript, [2] showed that *vanilla* DS reconstruction methods (such as vanilla NODEs) struggle to generalize across state space in a multistable (and hence) OOD scenario. In this context, this work is novel as it shows conditions under which OOD generalization is indeed possible.", "essential_references_not_discussed": "Not essential per se, but [1] backs the findings in that they also find that a similar hierarchical approach to CoDA, but using an recurrent neural network backbone, leads to linear relationships between GT control parameter and context vectors with capabilities of inter-and-extrapolation.", "other_strengths_and_weaknesses": "**Strengths**: I think the paper does a great job of motivating, explaining and presenting each and every experiment in a detailed fashion.\n\n**Weaknesses**: A current weakness of the manuscript is that the hypothesized mechanism for the possible OOD generalization, i.e. the Poincaré-Hopf theorem,  is only validated for NODE backbone based methods and tested on Hamiltonian systems. To remedy this weakness, the authors could try different DSR backbones, e.g. RNNs, and apply the regularizer to different, non-Hamiltonian DS.", "other_comments_or_suggestions": "**Typo**: “This raises a pertinent question arises:” (p. 5, ll. 273-274)\n\n**References**\n\n[1] Brenner, M. et al. (2024). Learning Interpretable Hierarchical Dynamical Systems Models from Time Series Data. (ICLR 2025)\n\n[2] Göring, N.A. et al. (2024). Out-of-Domain Generalization in Dynamical Systems Reconstruction. (ICML 2024)", "questions_for_authors": "1. The choice of $\\chi_{PH}$ (desired index) is a quantity that does depend on physical priors we have about the observed DS, correct? If so, how can the regularizer help in cases where this information is a priori *not* available?\n2. I’m struggling a bit with section 4.3; without knowledge of the GT system, is there a way to identify whether the model is hallucinating the wrong bifurcation?", "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": "xH9YXQPL8b", "reviewer": "Reviewer_xWSr", "summary": "The paper introduces Context-Informed Neural ODEs (CI-NODEs), a framework designed to learn dynamical systems exhibiting bifurcation behaviors, particularly symmetry-breaking bifurcations. The authors claim that CI-NODEs, despite being trained solely on pre-bifurcation, symmetric data, can predict post-bifurcation symmetry-breaking behaviors without explicit physics-based priors. They attribute this phenomenon to the implicit use of topological invariants, specifically the Poincaré index, and provide a formal explanation via the Poincaré-Hopf theorem. Additionally, a novel topological regularization term inspired by this theorem is proposed and tested on the Landau-Khalatnikov system to enhance generalization.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "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": "We sincerely appreciate the reviewer’s thoughtful and constructive feedback. Below, we address each comment carefully. Additional experimental results are available in the **README** at https://github.com/anonymous-account123/icml2025-7637 and will be thoroughly incorporated into our revised paper.\n***\n**Q1. Reliability of the zero-shot identification.**\n\nWe fully agree with your comment that evaluating the reliability of bifurcation identification is of critical importance. To statistically validate the robustness of zero-shot bifurcation identification, we conducted repeated experiments (five runs) using different random seeds and initial conditions. Consistently, the model reliably identified post-bifurcation behaviors, demonstrating statistical robustness and reproducibility (**Figure 1(a)** at the provided link). \n\nIn addition, we performed additional robustness tests under realistic scenarios, including noisy observations and limited training data. These challenging settings revealed that the model remains capable of accurately identifying symmetry-breaking bifurcations, albeit with slightly increased variance compared to the ideal noise-free scenario (**Figure 1(b-c)** at the provided link).\n***\n**Q2. Hallucinated bifurcations and empirical validation of the theorem.**\n\nThe hallucinated symmetry-breaking scenario in Section 4.3 is purposefully designed to illustrate the predictive capability of Proposition 4.1. In the analyzed linear system, a simple center-to-saddle bifurcation is expected if the model merely approximates the functional form. However, our empirical results demonstrate a spurious symmetry-breaking bifurcation, as predicted by Proposition 4.1, confirming that the model genuinely leverages topological constraints rather than simple functional approximation. We will clarify this critical interpretation in the revised paper as you suggested.\n***\n**Q3. How to detect and avoid hallucinated bifurcations.**\n\nWe completely agree that identifying whether a model hallucinates bifurcations is vital. We propose a straightforward yet robust criterion: evaluating the variance in bifurcation diagrams generated by multiple independent training runs. Correctly identified bifurcations show minimal variance, indicating stable and consistent identification. Conversely, hallucinated bifurcations yield significantly higher variance due to their structural instability. **Figure 2** at the provided link clearly demonstrates this differentiation approach.\n***\n**Q4. Other type of bifurcations or dynamical systems.**\n\nFollowing your suggestion, we investigated a non-Hamiltonian dynamical system exhibiting a codimension-2 cusp bifurcation: $(\\dot{x}, \\dot{y}) = (y, b + a x - x^3 – d y)$, with variable parameters $(a, b)$ and fixed $d = 0.5$. This system serves as a canonical model for capturing catastrophic transitions and hysteresis phenomena, which are fundamental in science and engineering. Remarkably, despite being trained only on the simple, monostable pre-bifurcation regime $(a,b) \\in \\[-2.0, -1.5, -1.0, -0.5 ]^2$, the model successfully identified the cusp bifurcation surface, validating the model's broader applicability. Detailed visualizations are available in **Figure 3** at the provided link.\n***\n**Q5. Comparison with PINNs.**\n\nWe acknowledge that PINNs effectively utilize explicit physical priors. However, they become inapplicable if physical laws are unknown or incorrectly assumed. Our proposed approach circumvents this limitation by leveraging general topological invariants. For instance, models representing vector fields on a sphere (e.g., climate models) naturally adopt a global index constraint of $\\chi(S^2) = +2$. Such universally applicable prior knowledge positions our approach as complementary and potentially advantageous in scenarios where explicit physical laws are unavailable. This point will be discussed in the revised manuscript.\n***\n**Q6.  Computational complexity.**\n\nModeling $N$ environments using vanilla NODEs with $M$ parameters requires $N \\times M$ parameters. Conversely, the CI-NODEs reduce this requirement to approximately $M(1 + K)$, where $K$ represents the context dimensionality, showing substantial parameter efficiency, particularly when $N > K$. \n\nRegarding computational overhead from topological regularization, our experiments confirm a moderate increase in computation time (from 53.5 ms/epoch for vanilla models to 80.9 ms/epoch for regularized models under the same environment). Nevertheless, the regularized model benefits from improved learning stability and accordingly demonstrates faster empirical convergence, as shown in **Figure 4** at the provided link. We will clearly address this trade-off in the revision.\n***\n**Q7. Clarity of the paper.**\n\nWe appreciate your suggestions for clarity improvement. We will refine the main text and Appendix to highlight essential details clearly, and will include the additional discussions and experiments made during the rebuttal.", "claims_and_evidence": "The claim that CI-NODEs can \"identify\" bifurcations and predict post-bifurcation behaviors is intriguing but lacks rigorous justification. The results indicate that the model sometimes hallucinates bifurcations, which contradicts the assertion that the approach reliably generalizes.\n\nThe use of the Poincaré-Hopf theorem to explain the model’s behavior is an interesting theoretical contribution, but the practical implications remain unclear. There is no empirical validation showing that the theorem correctly predicts when the model will succeed or fail.\n\nThe claim that the proposed regularization enhances generalization is weakly supported. The experiments show some improvements, but they lack robustness tests, such as evaluations across a broader range of bifurcation scenarios.", "methods_and_evaluation_criteria": "The experimental setup is well-structured but suffers from narrow validation. The datasets used for training and testing are limited to specific types of bifurcations, making it unclear how well the model generalizes to other dynamical systems.\n\nThe chosen evaluation metrics (e.g., Mean Squared Error, trajectory consistency) are appropriate but insufficient to fully assess model reliability. Additional evaluations such as robustness to noise, long-term stability, and sensitivity to initial conditions would be beneficial.\n\nThe ablation studies focus mainly on the impact of context-informed modeling but do not adequately assess individual architectural choices, such as the specific design of the fusion mechanism.", "theoretical_claims": "The theoretical contribution is a highlight of the paper but remains somewhat speculative. While the Poincaré-Hopf theorem provides an insightful perspective, the connection between the theorem and the model’s emergent properties is not formally established.\n\nThere is no rigorous proof explaining why CI-NODEs should be able to generalize beyond training data in the specific manner observed. The argument based on topological constraints is heuristic at best.\n\nThe discussion of the relationship between NODEs and dynamical systems theory is useful but lacks depth in explaining how these theoretical insights could be leveraged to improve practical performance.", "experimental_designs_or_analyses": "The experiments provide some compelling qualitative results, but the lack of quantitative robustness testing is a significant weakness.\n\nThe comparisons with baseline methods, such as traditional NODEs, are somewhat superficial. More rigorous benchmarking against established bifurcation detection techniques would strengthen the claims.\n\nThe zero-shot generalization experiments are interesting but suffer from a lack of statistical analysis. How often does the model correctly infer post-bifurcation behavior versus hallucinating incorrect structures?", "relation_to_broader_scientific_literature": "The claim that this work offers new insights into OOD generalization in dynamical systems is overstated. While the results are interesting, they do not establish a broadly applicable principle for OOD learning.", "essential_references_not_discussed": "None.", "other_strengths_and_weaknesses": "**Strengths**:\n\nThe paper presents a novel perspective on learning bifurcating dynamical systems using NODEs.\n\nThe use of topological invariants as an implicit learning signal is an interesting theoretical contribution.\n\nThe proposed topological regularization is conceptually novel and may inspire further research in constrained learning for dynamical systems.\n\n**Weaknesses**:\n\nThe empirical validation is not rigorous enough to support the paper’s ambitious claims.\n\nThe theoretical arguments, while insightful, are largely heuristic and lack formal proof.\n\nThe experimental design does not adequately explore failure cases, making it difficult to assess the model’s reliability.\n\nThe paper is dense and difficult to follow, with crucial details buried in the appendix.", "other_comments_or_suggestions": "(1) Provide statistical significance testing for key experimental results.\n\n(2) Clarify how computational complexity scales with increasing problem size and longer time horizons.\n\n(3) Conduct additional experiments to test the model’s robustness to noise and distribution shifts.\n\n(4) Improve clarity by restructuring the theoretical discussion to make the key insights more accessible.", "questions_for_authors": "(1) Can you provide more evidence that CI-NODEs do not hallucinate bifurcations in cases where no symmetry breaking occurs?\n\n(2) How does CI-NODE compare to physics-informed methods, such as PINNs, in terms of accuracy and interpretability?\n\n(3) What guarantees, if any, can be provided regarding the reliability of CI-NODEs in predicting post-bifurcation behaviors?\n\n(4) Have you tested CI-NODEs on higher-dimensional bifurcation problems, and if so, how does it perform?\n\n(5) How does the choice of context representation affect the model’s ability to generalize?", "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:

Meta-review: After the discussion, the reviews are all positive. The paper proposes a method for NODE network to learn symmetries in dynamical systems as supported by a complete theoretical analysis. It clearly reaches the level of quality and novelty for ICML.
Decision: Accept (poster)