Interpretable Relational Inference with LLM-Guided Symbolic Dynamics Modeling
Xiaoxiao Liang * 1 Juyuan Zhang * 1 Liming Pan 1 Linyuan Lü 1
arXiv:2604.12806v1 [cs.LG] 14 Apr 2026
Abstract
Relational Inference
Inferring latent interaction structures from observed dynamics is a fundamental inverse problem in many-body interacting systems. Most neural approaches rely on black-box surrogates over trainable graphs, achieving accuracy at the expense of mechanistic interpretability. Symbolic regression offers explicit dynamical equations and stronger inductive biases, but typically assumes known topology and a fixed function library. We propose COSINE (Co-Optimization of Symbolic Interactions and Network Edges), a differentiable framework that jointly discovers interaction graphs and sparse symbolic dynamics. To overcome the limitations of fixed symbolic libraries, COSINE further incorporates an outerloop large language model that adaptively prunes and expands the hypothesis space using feedback from the inner optimization loop. Experiments on synthetic systems and large-scale real-world epidemic data demonstrate robust structural recovery and compact, mechanism-aligned dynamical expressions. Code: https://anonymous. 4open.science/r/COSINE-6D43.
Symbolic Regression 풅 = �( ) 풅
COSINE 풅 = �� �� �� ( 풅
Figure 1. Conceptual comparison of relational inference, symbolic regression, and COSINE. Relational inference recovers latent interaction graphs from dynamical data; symbolic regression discovers governing equations of dynamics given known structures; COSINE integrates sparse symbolic modeling with graph learning to jointly infer both structure and dynamics.
As illustrated in Figure 1, relational inference estimates pairwise dependencies from dynamical data. Classical methods, including Granger causality (Granger, 1969), mutual information (Wu et al., 2020), and transfer entropy (Schreiber, 2000), focus on statistical dependencies but neglect explicit dynamical modeling. Neural approaches, such as Neural Relational Inference (NRI) and its extensions (Kipf et al., 2018; Graber & Schwing, 2020; Pan et al., 2024; Löwe et al., 2022; Chen et al., 2021), jointly infer latent structures and dynamics. Attention-based models (Lu et al., 2023; Wu et al., 2025) improve efficiency via attention, they remain limited by a lack of mechanistic transparency. Critically, these black-box neural methods tend to exploit over-parameterization to fit spurious edges, which accommodates the observed data at the expense of discovering the true governing mechanisms.
1. Introduction Interacting structures among constituent entities are fundamental to understanding dynamical systems across physics, biology, neuroscience, epidemiology, finance, and sociology (Arenas et al., 2008; Zhu et al., 2022; Izhikevich, 2007; Pastor-Satorras et al., 2015; Moraffah et al., 2021; Castellano et al., 2009). In many settings, graphs are unobserved due to experimental, ethical, or economic constraints, while time-resolved dynamics are measurable, motivating inference of latent interactions from trajectories.
Symbolic regression (SR) (Schmidt & Lipson, 2009) instead targets explicit governing equations from dynamics (Udrescu & Tegmark, 2020; Cranmer et al., 2019; 2020; Wong & Cranmer, 2022; Kamienny et al., 2022). Concise expressions impose strong inductive bias and improve mechanistic readability; expression-level structure is further supported by recent compositional models, including KAN and related studies (Liu et al., 2024; Wetzel et al.,
*
1 Equal contribution School of Cyber Science and Technology, University of Science and Technology of China, Hefei, Anhui, China. Correspondence to: Liming Pan <pan [email protected]>, Linyuan Lü <[email protected]>.
Preprint. April 15, 2026.
1
Co-Optimization of Symbolic Interactions and Network Edges
2025; Makke & Chawla, 2024). This compactness acts as a functional constraint: incorrect edges are harder to absorb through over-parameterization.
2. Related Work
Despite these advantages, symbolic regression on network dynamics remains challenging when networks are unknown. As shown in Figure 1, most existing SR methods assume known network (Grayeli et al., 2024; Yu et al., 2025; Ruan et al., 2025; Tsoi et al., 2025). With latent structure and dynamics, discovery becomes a coupled search; SINDy reduces equation discovery to sparse regression over fixed libraries (Brunton et al., 2016a; Basiri & Khanmohammadi, 2024; Hoffmann et al., 2019; Stanković et al., 2020; Zhao et al., 2025; Goyal & Benner, 2022; Lai, 2021), but closedworld libraries require domain knowledge and can induce spurious edges and degrade robustness.
Inferring latent interaction structures from multivariate dynamical data has been widely studied. Classical techniques characterize interaction structures in observational data using statistical measures, including Granger causality (Granger, 1969), correlation and partial correlation (Peng et al., 2009), mutual information (Wu et al., 2020), and transfer entropy (Schreiber, 2000). These approaches focus on statistical dependencies, neglect explicit dynamical modeling, and fail under complex nonlinear dynamics.
2.1. Relational Inference
Neural Relational Inference (NRI) (Kipf et al., 2018) introduced a unified framework that jointly infers latent graphs and system dynamics via end-to-end learning. An encoder encodes graph from observed trajectories, while a decoder predicts future states conditioned on the inferred graph. Subsequent extensions address time-varying graphs (Graber & Schwing, 2020), heterogeneous interactions (Webb et al., 2019; Ha & Jeong, 2023), directed graph inference (Wang & Pang, 2022), and sparsity/degree priors (Li et al., 2019; Chen et al., 2021). Related work also integrates attention-based architectures to capture fine-grained, potentially asymmetric dependencies (Ha & Jeong, 2020; 2022; Lu et al., 2023; Wu et al., 2025). However, despite their expressive power, these methods largely treat dynamics as a black box, limiting their interpretability and making them susceptible to spurious associations with over-parameterization.
To address these challenges, we propose COSINE (CoOptimization of Symbolic Interactions and Network Edges), a differentiable framework for joint structure and mechanism discovery (Figure 1). COSINE decomposes dynamics into reusable message and update components as sparse combinations over shared symbolic libraries, enabling expression reuse and a search space independent of network size. To relax fixed-library bias, COSINE uses an outerloop large language model (LLM) as a symbolic hypothesis proposal mechanism to prune and augment terms from inner-loop feedback, without system-specific templates. Our contributions are summarized as follows: • COSINE Framework: We propose an end-to-end differentiable pipeline that jointly optimizes latent graph structures and sparse dynamical equations, enabling the discovery of interpretable dynamical models under unknown networks.
2.2. Symbolic and Interpretable Dynamical Modeling Sparse Symbolic Regression on Dynamics. Sparse Identification of Nonlinear Dynamics (SINDy) (Brunton et al., 2016a) identifies governing equations from dynamical systems by selecting a parsimonious subset of basis functions from a predefined library via sparse optimization. Numerous extensions improve its robustness under noise and higherorder integration schemes (Mangan et al., 2017; Brunton et al., 2016b; Goyal & Benner, 2022; Zhao et al., 2025). While network-aware variants like SINDyG (Basiri & Khanmohammadi, 2024; Gao & Yan, 2022) exist, they typically assume a fixed or known graph structure, making them difficult to apply when the network is hidden and no domain knowledge is available.
• Sparse Symbolic Message Passing: We cast the joint discovery problem as sparse regression over message/update libraries. This constraint prevents overparameterized dynamics from masking spurious edges, thereby enhancing structural identifiability. • LLM-Guided Library Evolution: We utilize a feedback-driven pruning/augmentation strategy to adaptively adjust the basis library without hand-crafted, system-specific templates. This approach decouples LLM-driven hypothesis generation from optimizationbased hypothesis selection, ensuring the discovered mechanisms remain numerically grounded.
Symbolic Discovery on Network Dynamics. To explore broader symbolic spaces, neural-guided symbolic regression has emerged, utilizing evolutionary search (Schmidt & Lipson, 2009; Cranmer et al., 2020), neural approximations (Udrescu & Tegmark, 2020), or Transformer-based generation (Kamienny et al., 2022; Biggio et al., 2021). For networked systems, recent methods like ND2 (Yu et al., 2025) and PSE (Ruan et al., 2025) decompose dynamics into node-, interaction-, and aggregation-level operators, while
• Extensive Experiments: Evaluations across diverse synthetic dynamical systems and real-world datasets demonstrate that COSINE consistently achieves stateof-the-art performance in relational inference while successfully recovering parsimonious and physically consistent symbolic governing laws. 2
Co-Optimization of Symbolic Interactions and Network Edges
LaSR (Grayeli et al., 2024) leverages large-model-assisted reasoning to refine hypotheses, introducing an LLM into its genetic algorithm to improve performance.
the temperature-controlled Gumbel–Softmax to map latent variables to edge-type probabilities: 0 1 Aaij = softmax τ1 [Ψ0ij + gij , Ψ1ij + gij ] , (1)
Neural relational inference methods often lack interpretability and can overfit spurious interactions. Symbolic regression offers stronger inductive biases for dynamical systems, but becomes difficult when the graph is unknown, since the graph structure and equation form are tightly coupled and the search space explodes. SINDy alleviates symbolic search via sparsity, yet relies on a pre-defined operator library and domain priors. COSINE addresses this gap by jointly optimizing a latent interaction graph and sparse symbolic dynamics, enabling accurate structure inference together with interpretable mechanism discovery.
a
where a ∈ {0, 1} denotes the edge types, τ > 0 is the temperature parameter controlling the sharpness of the disa tribution, and gij denotes Gumbel noise to improve graph diversity. As τ decreases, the edge distribution becomes increasingly peaked and converges to a deterministic graph in the limit. This construction naturally generalizes to multi-relation settings by assigning multi-dimensional latent scores to each edge and applying a Softmax to obtain a categorical distribution over relation types. In subsequent dynamical modeling, we use the soft adjacency Asoft ij as the weight for message aggregation, enabling differentiable structure learning during training.
3. Method 3.1. Overall Framework
3.3. Sparse Regression Modeling
COSINE (Co-Optimization of Symbolic Interactions and Network Edges) aims to jointly discover the latent interaction graph and governing mechanisms from observed trajectories X = (x0 , . . . , xT −1 ). We consider N interacting variables, where the state of node i at discretized time t is xti ∈ RD . Let xt = {xt1 , . . . , xtN } be the global state and xi = (x0i , . . . , xTi −1 ) be the trajectory of node i. We write t,(d) (d) xi for the d-th component (omitting t in tables as xi ).
Complex networked dynamics can often be written as a combination of self-dynamics and interaction dynamics: N X d xi = fi (x, A) = W (xi ) + Aij Q(xi , xj ), dt j=1
(2)
where W (·) describes the intrinsic state evolution of each node, Q(·, ·) captures interactions with neighboring nodes, and Aij ∈ {0, 1} indicates the network. In classical neural approaches, a surrogate message-passing GNN is leveraged to learn the unknown dynamics from data. A v → e transformation generates edge features, which are then aggregated back to nodes to update node states (e → v):
As shown in Figure 2, COSINE operates as a two-level closed-loop system. The outer loop employs an LLM as a symbolic supervisor to iteratively refine the basis function library Θ(·) based on training feedback. Within the inner loop, we adopt a non-amortized setting where the graph A (parameterized by Ψ) and dynamical coefficients W are co-optimized via gradient descent to minimize the prediction error. This design enables the discovery of parsimonious, interpretable governing mechanisms without requiring system-specific prior knowledge.
v→e : e→v :
Unless otherwise specified, we train COSINE using one-step supervision (predicting xt+1 from xt ) under teacher forcing, which provides a stable and comparable objective across continuous- and discrete-time systems. At evaluation time, the learned mechanisms can also be rolled out for multiple steps to inspect qualitative trajectory behavior, while the inferred graph is assessed directly by edge-level metrics.
zl(i,j) = fel [zli , zlj , x(i,j) ] , X zl+1 = fvl [ zl(i,j) , xj ] , j
(3) (4)
i∈Nj
where zli denotes the embedding of node vi at layer l, zl(i,j) denotes the embedding of edge e(i,j) , and xi and x(i,j) represent the initial node and edge features, respectively. Nj denotes the set of neighbors with incoming edges to node j, and [·, ·] indicates vector concatenation. The functions fv and fe are node- and edge-specific neural networks, respectively. The final node embeddings zL are then used to predict the next-step state of the dynamical system.
3.2. Graph Generation Module The graph generator encodes unknown interaction structures into a set of trainable latent variables and maps them to an adjacency matrix. For binary relations, we let a ∈ {0, 1} denote the edge type. For each ordered pair (i, j), we maintain two real-valued latent variables Ψ0ij , Ψ1ij ∈ R. To obtain differentiable edge sampling and soft weights, we employ
Inspired by the message-passing decomposition in Equation (3), we design a sparse regression module that employs two basis-function libraries to approximate the two stages of GNN modeling. Specifically, we decompose the fitting process into a message flow and an update flow. The message 3
Co-Optimization of Symbolic Interactions and Network Edges System Prompt
Basis Library
Workflow System Trajectories
You are an expert in discovering dynamical systems and their interaction structures using sparse symbolic regression. Your Task: Act as an expert in dynamical systems, sparse regression, and PyTorch expressions. Only design / refine the basis function library. State dimension D = {feature_dim} Budget: at most {max_terms} terms in each stream
Training Objective
model interactions between nodes: �(�� , �� ) Update terms: Functions are used to model
Optimize the graph and basis term weights to fit data： �풕풐풕 � 풕풐풕
= �� � ++ � �� =� �풏 ++��풌 �� � 풏
풌
�
�
�
Graph Guided Dynamic Modeling �=
Softmax
Context Feedback Training Feedback Basis Library and term weights Example prompt:
Message terms: Functions are used to
Graph Generator
�
�풊� �(�풊 , �� ) �(�풊 , �� )
��
� Message terms
LLM Reasoning
��+�
Update terms
Equation Discovery �(�)
System Prompt
�
×
Context Feedback
...
You are optimizing basis function selection based on training feedback. Inputs provided with each call include: Recent 3 rounds (each round has loss + Message terms + Update terms; each term includes name, weight, expr)
��+� = �� + �(�, �)
풊
= �
each node’s own state update: �(�, �) {
"message_terms": [{ "name": "prod_self", "expr": "xi * xj", "type": "vector" }， { "name": "lag_interact", "expr": "xj[..., 1:2] * (xj[..., 0:1] - xi[..., 0:1])", "type": "vector" }, { "name": "diff_norm_scaled", "expr": "diff / (torch.norm(diff) + 1e-6)", "type": "vector" }], "update_terms": [{ "name": "saturation", "expr": "torch.tanh(x)", "type": "vector" }, { "name": "self_decay", "expr": "x * h", "type": "vector" }, { "name": "history", "expr": "h[..., 1:2]", "type": "vector" }, { "name": "h_sigmoid", "expr": "torch.sigmoid(h)", "type": "vector" }] }
Figure 2. The COSINE architecture. (Left) LLM-based Reasoning refines the basis library Θ(·) via performance feedback. (Middle) Graph-Guided Modeling co-optimizes latent edges Aij and coefficients W through differentiable symbolic message-passing. (Right) Symbolic Basis Library bridges high-level reasoning with numerical discovery of governing mechanisms.
Update term. Given an update basis set Θupd = upd {ϕupd 1 , . . . , ϕU }, we construct update features for node t i using the current state aggregated message hti , and Pxi , the soft the node degree ki = j Aij :
flow generates interaction messages from neighbors, while the update flow updates node states based on aggregated messages. By modeling each component using a sparse linear combination of basis functions, we can explicitly capture key mechanisms and interpret the equation structure through coefficient sparsity: eij = MESSAGE(xi , xj ), xt+1 = UPDATE xti , AGGj∈Nin (i) etij . i
t t D ϕupd n (xi , hi , ki ) ∈ R ,
(6)
∆xti =
M X
t t D Wmsg [m] ⊙ ϕmsg m (xi , xj ) ∈ R ,
t t Wupd [n] ⊙ ϕupd n (xi , hi , ki ).
(11)
Thus, the one-step prediction is x̂t+1 = xti + ∆xti . i
(12)
For continuous-time systems, this can be viewed as an explicit Euler discretization where ∆xti approximates ∆t ẋti (with ∆t absorbed into the learned coefficients). For discrete-time systems, ∆xti acts as a residual update that directly models the transition xti 7→ xt+1 . i
(7)
where D is the node feature dimension. A linear combination with learnable weights yields the edge message vector etij =
U X n=1
Message term. Given a message basis function set Θmsg = {ϕ1 , . . . , ϕM }, for each directed pair (i, j) we compute m = 1, . . . , M,
(10)
and obtain the final increment
(5)
Here Nin (i) denotes the incoming neighbors of node i under the adjacency convention above (i.e., indices j such that Aij = 1).
t t D ϕmsg m (xi , xj ) ∈ R ,
n = 1, . . . , U,
Joint optimization manifests at two levels: (i) the structural layer, where Asoft ij represents the probability or strength of each edge; and (ii) the equation layer, where the sparsity of Wmsg and Wupd determines which basis functions truly participate in generating the dynamics.
(8)
m=1
where ⊙ denotes element-wise multiplication and Wmsg [m] ∈ RD is the coefficient vector for the m-th basis. Given the soft structure Asoft , messages are aggregated for node i as N X t hti = Asoft (9) ij eij .
3.4. Training Objective and Loss Given the one-step prediction x̂t+1 and the ground-truth observation xt+1 , we employ the Gaussian negative loglikelihood (NLL) as the primary objective, which is equivalent to minimizing the scaled squared error under a fixed
j=1
4
Co-Optimization of Symbolic Interactions and Network Edges
variance σ 2 . The prediction loss is formulated as: Lnll =
N 1 X t+1 2 x̂ − xt+1 . i 2 2σ 2 i=1 i
4. Experiments 4.1. Results on Synthetic Data
(13)
Networks. We conduct experiments on three standard graph families: Erdős–Rényi (ER), Barabási–Albert (BA), and Watts–Strogatz (WS). For ER graphs, the edge probability is set to p = 0.1; for BA graphs, each new node attaches to m = 2 existing nodes; for WS graphs, we start from a ring lattice with k = 2 nearest neighbors per node and rewire edges with probability p = 0.1.
To stabilize structure learning and impose a prior over edgetype distributions, we include a KL regularization term (computed over all ordered node pairs): LA = KL q(A) ∥ r , (14) where q(A) denotes the edge distribution induced by the graph generator, and r is a simple factorized prior over edge types (e.g., a Bernoulli prior favoring sparsity when appropriate; if no prior knowledge is available, we default to a uniform prior). To encourage sparsity in symbolic regression, we apply ℓ1 regularization to both sets of coefficients: LW = ∥Wmsg ∥1 + ∥Wupd ∥1 .
Dynamical systems. We evaluate COSINE on six representative systems covering biochemical reactions, diffusion, mechanical coupling, synchronization, opinion dynamics, and chaotic maps: (1) Michaelis-Menten kinetics (MM) (Karlebach & Shamir, 2008); (2) Diffusion (Diff); (3) Network of springs (Spr); (4) Kuramoto model (Kura) (Kuramoto, 1975); (5) Friedkin-Johnsen dynamics (FJ) (Friedkin & Johnsen, 1990; Abebe et al., 2018); and (6) Coupled map network (CMN) (Garcia et al., 2002). These systems, generated on ER, BA, and WS networks, span continuous and discrete time with diverse nonlinear interactions.
(15)
The total loss is Ltotal = Lnll + βKL LA + λW LW ,
(16)
where βKL and λW are hyperparameters. The full training process is shown in Section A.1. After each inner-loop optimization, we summarize the current library Θ performance together with the learned coefficient weights and feed this information to the LLM to guide an evolution of the basis library for the next inner-loop.
Evaluation metrics. We report AUC as our primary metric for relational inference, computed by comparing the predicted edge scores Asoft ij against the binary ground-truth adjacency. In addition, for mechanism discovery, we compute term accuracy as primitive coverage: for each module (Message/Update), we take the top K = 3 terms by absolute coefficient magnitude and measure the fraction of required ground-truth primitives recovered among them (using simple equivalence rules, e.g., matching sinusoidal coupling on xj − xi up to sign/constant scaling).
3.5. LLM-Guided Library Evolution The outer loop performs a lightweight search over symbolic space by editing the basis library. Each round uses inner-loop feedback (loss, sparsity, and residual patterns) to prompt the LLM to prune or augment terms, producing a candidate library. To avoid expensive multi-branch search while keeping evolution controlled, we use a best-so-far strategy that maintains a single optimal library and only accepts candidates that strictly improve prediction loss.
Baselines. We compare COSINE with various baselines from three categories: classical statistical methods (GC (Granger, 1969), MI (Butte & Kohane, 1999), and TE (Schreiber, 2000)), neural relational inference models (NRI (Kipf et al., 2018) and GDP (Pan et al., 2024)), and attention-based models (RIVA (Wu et al., 2025)). We adapt RIVA to the unweighted setting and tune hyperparameters for all baselines using their public implementations to ensure fair comparison in experiments.
Formally, we track the best library across rounds: Θ⋆ = arg
min Θ∈{Θ(0) ,...,Θ(r) }
Lval (Θ),
(17)
e only if Lval (Θ) e < Lval (Θ⋆ ). This and update Θ⋆ ← Θ strategy keeps computation predictable and reduces the risk of degradation from noisy edits, while leveraging strong numerical feedback to steer the symbolic space toward parsimonious, physically consistent mechanisms. In practice, the validation gate decouples symbolic exploration from numerical selection: the LLM only proposes candidates, while the inner-loop loss determines acceptance. If no candidate improves Lval over successive rounds, the search naturally terminates without accumulating regressions. In Section A.2 Section A.3, details of the library evolution and LLM (prompts, model) are presented.
Implementation details. For COSINE, the initial basis library is generated by the LLM without system-specific templates. We set the Gumbel–Softmax temperature to τ = 0.3, βKL = 0.1, and λW = 0.1, and optimize with Adam (lr 0.1 for structure and 0.005 for coefficients). Each outer round runs 1000 inner epochs for up to 10 rounds, using validation loss and coefficient statistics to refine the library under a controlled edit budget and a maximum library size. We use GPT-OSS (20B) as the default LLM. More details of synthetic experiments are presented in Section B. 5
Co-Optimization of Symbolic Interactions and Network Edges Table 1. Relational inference performance (AUC %) on different graphs with N = 50 nodes across various dynamical systems. In the VOLUME column, a×b corresponds to trajectories × sampled steps. Boldface marks the highest accuracy. DYN .
G RAPH
VOLUME
GC
MI
TE
NRI
GDP
RIVA
COSINE
MM
ER-50 BA-50 WS-50
50 × 10 50 × 10 50 × 10
60.24 62.79 80.33
77.84 88.04 97.77
53.66 63.96 57.56
96.25±2.22 82.67±1.81 99.83±0.15
98.31±1.41 93.02±3.94 56.77±2.55
52.17±1.15 52.83±1.22 53.65±1.40
99.63±0.12 98.07±0.45 99.10±0.25
D IFF
ER-50 BA-50 WS-50
50 × 10 50 × 10 50 × 10
63.84 75.54 73.63
56.00 72.06 76.39
57.63 61.71 90.04
91.87±7.80 94.16±11.62 99.31±0.55
93.44±4.87 94.41±3.23 56.68±3.12
53.02±1.35 52.89±1.18 53.72±1.25
99.36±0.10 96.55±1.05 100.00±0.00
S PR
ER-50 BA-50 WS-50
15 × 10 15 × 10 15 × 10
50.61 54.53 55.44
72.24 91.16 50.40
76.05 84.67 51.83
99.84±0.47 98.17±5.40 96.77±1.25
99.99±0.02 99.88±0.36 99.94±0.05
52.82±1.50 54.16±1.62 55.20±1.45
100.00±0.00 100.00±0.00 100.00±0.00
K URA
ER-50 BA-50 WS-50
30 × 30 30 × 30 30 × 30
54.94 51.50 59.27
64.69 55.46 100.00
64.76 61.87 96.63
82.09±19.14 69.70±18.16 50.17±1.05
94.93±12.94 90.13±12.38 100.00±0.00
75.98±2.34 76.06±2.15 98.19±0.55
99.99±0.01 99.85±0.14 100.00±0.00
FJ
ER-50 BA-50 WS-50
20 × 10 20 × 10 20 × 10
93.00 92.87 99.81
53.66 52.32 66.14
83.64 86.88 95.37
97.67±1.06 91.62±4.67 68.76±5.25
99.82±0.47 92.63±13.46 89.51±3.44
52.55±1.10 52.64±1.32 52.89±1.45
100.00±0.00 99.98±0.02 100.00±0.00
CMN
ER-50 BA-50 WS-50
20 × 10 20 × 10 20 × 10
52.89 53.44 61.08
87.39 87.84 98.94
64.35 71.51 86.70
89.76±2.59 83.35±2.30 69.21±4.12
97.58±3.38 88.83±6.19 51.51±1.10
50.84±0.95 53.87±1.67 55.96±1.88
99.75±0.15 93.23±2.45 99.05±0.55
Table 2. Dominant terms discovered by COSINE on ER graphs (Table 1), ranked by absolute coefficient magnitude. Term 1–3 denote the top three terms for Message and Update modules. Term Acc. reports row-wise primitive coverage scores. x: node state; hi : aggregated message; ki : degree; ϵ = 10−6 . DYN . MM D IFF S PR
G ROUND - TRUTH DYNAMICS P A xj ẋi = −xi + j∈N (i) kij i 1+xj
C OMPONENT
T ERM 1
T ERM 2
T ERM 3
T ERM ACC .
M ESSAGE U PDATE
(xj − xi )2 (0.31) hi (0.70) ki +ϵ
|xj − xi | (0.28) x (0.42)
xi xj (0.20) sin(x) (0.25)
0.00 1.00
M ESSAGE U PDATE
xj − xi (0.17) x+hi (0.31) ki +1+ϵ
xi xj (0.12) hi (0.30) ki +1+ϵ
cos(xj − xi ) (0.11) cos(x) (0.08)
1.00 1.00
M ESSAGE U PDATE
xj − xi (0.02) hi (0.04) ki +1
(0) (1)
j∈N (i) Aij (pi − pj )
xi xj
sin(xi xj ) (0.01) x (0.02)
j∈N (i) Aij κ sin(ϕj − ϕi )
M ESSAGE U PDATE
sin(xj − xi ) (0.53) hi (0.48) ki +1
|xj − xi | (0.16) xhi (0.19) ki +1
M ESSAGE U PDATE
xj − xi (0.25) hi (0.97) ki +1
M ESSAGE U PDATE
xi xj (0.33) hi (2.66) ki +1
ẋ = −βLx, p̈i = −k
P
P
K URA
ϕ̇i = ωi +
FJ
xt+1 = ki1+1 i
CMN
L = DIN − A
P
t j∈N (i) Aij xj + si
θit+1 = (1 − s)f (θit ) + ksi
t j∈N (i) Aij f (θj )
P
Relational Inference. Table 1 shows that COSINE achieves state-of-the-art or highly competitive performance, significantly outperforming statistical methods (GC, MI, TE) and matching or surpassing deep learning models (NRI, GDP, RIVA). While statistical approaches struggle with high-order nonlinearities, COSINE’s joint optimization of structure and symbolic dynamics ensures robust reconstruction across diverse regimes. Notably, in complex scenarios like Kuramoto system on BA graphs, COSINE maintains a near-perfect AUC (99.85), substantially exceeding NRI and GDP. Although performance margins can be influenced by specific graphs, in practice, across settings, COSINE consistently demonstrates superior effectiveness in disentangling interactions within heterogeneous networks.
(0.01)
x (0.03) ki +1
(0) (1)
1.00 1.00
xi xj (0.14) hi (0.08)
(1) (2)
1.00 1.00
xi xj (0.03) x (0.06)
xi xj (0.03) hi (0.06)
1.00 1.00
(xj − xi )2 (0.29) tanh(hi ) (1.18)
xj − xi (0.12) x (1.12)
0.50 1.00
(0) (1)
Mechanism Discovery. As shown in Table 2, COSINE consistently recovers mechanism-aligned primitives, with a clear strength in identifying aggregation and normalization patterns in the Update module. COSINE captures intended interaction primitives such as xj − xi and sinusoidal coupling, while the Message module occasionally includes predictive proxy nonlinearities when multiple bases can approximate similar interaction effects under finite data and noise. Overall, the alignment between inferred terms and groundtruth mechanisms indicates that COSINE goes beyond pure statistical association toward interpretable mechanism discovery; its goal is interpretable joint discovery of both latent structure and dynamics, providing sparse, human-readable interaction primitives and aggregation forms. 6
Co-Optimization of Symbolic Interactions and Network Edges Mobility Graph Inferred Graph
Message