Deep Signature: Characterization of Large-Scale Molecular Dynamics

## Abstract

Abstract Understanding protein dynamics are essential for deciphering protein functional mechanisms and developing molecular therapies. However, the complex high-dimensional dynamics and interatomic interactions of biological processes pose significant challenge for existing computational techniques. In this paper, we approach this problem for the first time by introducing Deep Signature, a novel computationally tractable framework that characterizes complex dynamics and interatomic interactions based on their evolving trajectories. Specifically, our approach incorporates soft spectral clustering that locally aggregates cooperative dynamics to reduce the size of the system, as well as signature transform that collects iterated integrals to provide a global characterization of the non-smooth interactive dynamics. Theoretical analysis demonstrates that Deep Signature exhibits several desirable properties, including invariance to translation, near invariance to rotation, equivariance to permutation of atomic coordinates, and invariance under time reparameterization. Furthermore, experimental results on three benchmarks of biological processes verify that our approach can achieve superior performance compared to baseline methods.

## 1 Introduction

Biological processes are fundamentally driven by the dynamical changes of macromolecules, particularly proteins and enzymes, within their respective functional conformation spaces. Typical examples of such processes include protein–ligand binding, molecule transport and enzymatic reactions, and modern computational biologists investigate their underlying functional mechanisms by molecular dynamics (MD) simulations (Dror et al., ; Lewandowski et al., ). Built upon density functional theory (Car & Parrinello, ), MD has demonstrated remarkable capability in providing accurate atomic trajectories in three-dimensional (3D) conformational space and consist agreement with experimental observations (Frenkel & Smit, ).

The computational analysis of MD data has been a subject of extensive research for decades, with the goal of characterizing systems from trajectory information. However, due to the main challenge posed by intricate interatomic interactions over large-scale systems across inconsistent timescales, many existing works resort to oversimplified setups that incorporate biophysical priors to analyze certain aspects of dynamics such as protein fluctuations, relaxation time, stability, and state transitions (Law et al., ; Qiu et al., ). More recently, empowered by the parallel processing ability of GPUs, machine learning especially deep learning sheds new light on this field as it can discretize macromolecules as particles distributed in a 3D voxel grid and automatically learn their relations in a data-driven fashion (Li et al., ; Rogers et al., ). Parallel to this, surface modeling-based approaches have emerged, firstly utilizing mathematical models to restore protein surfaces and then applying deep learning to analyze the chemical and geometrical features of surface regions around binding sites (Gainza et al., ; Zhu et al., ). Despite the great potential of these approaches in automatic drug discovery, their computational complexity would increase linearly with the number of time stamps when processing MD data, struggling with application to long-time simulations. Besides, these methods commonly build upon coarse grained dynamics for accelerating computation. Nevertheless, selecting an optimal coarse graining mapping strategy that effectively simplifies the representation of the system while preserving essential features remains an open research problem (Jin et al., ; Majewski et al., ).

Another limitation of current MD analysis methods is the deficient utilization of structural bioinformatics for the largely increased difficulty in handling high-order interatomic interactions during dynamic processes. However, such structural bioinformatics, manifested in various covalent and non-covalent bonds, plays a pivotal role in molecular design for its capability of propagating local perturbations to facilitate conformational dynamics and alter biological function (Tsai et al., ; Otten et al., ).
An illustrative example would be dihydrofolate reductase, which has been widely studied as important antitumor and antibacterial targets for treating tuberculosis and malaria. There exist four common mutations that confer drug resistance to antibiotics, proceeding in a stepwise fashion. Among them, the P21L mutation acts in a dynamical loop region associated with long range structural vibrations of the protein backbone, rather than directly on the active sites as other mutations (Toprak et al., ). Therefore, ignoring such interatomic interactive dynamics facilitated by molecular structure for a critical protein and counting the effects of active sites solely can result in biased assessments of designed drugs. Nonetheless, since the integration of structural bioinformatics into MD analysis would further introduce at least quadratic complexity with system size, existing works have not yet investigated this crucial aspect, highlighting a critical gap in our ability to comprehensively analyze and predict molecular behavior in drug design and resistance studies.

To this end, we aim to develop a computationally efficient framework that incorporates the structural bioinformatics with coarse graining mapping for automatically analyzing protein trajectory dynamics.
In particular, we first introduce a graph clustering module that learns to extract coarse grained dynamics by approximating soft spectral clustering. With the clustering assignment function implemented by a graph neural network and parameters learned automatically, we circumvent the need for manual selection of coarse graining mapping. Subsequently, we introduce a path signature transform module served as a feature extractor to characterize the interatomic interactive dynamics after coarse graining. Path signature is a mathematically principled concept that utilizes iterated integrals to describe geometric rough paths in a compact yet rich manner (Lyons, ), thus suitable for our tasks where molecular trajectories are highly sampled and non-smooth. After attaching with a task-specific differentiable classifier or regressor, we devise an end-to-end framework, named Deep Signature, for efficiently characterizing the complex protein dynamics. Notably, due to the existence of considerable random fluctuation in simulated trajectories, ideal features ought to maintain symmetry respecting certain geometric transformations. We provide theoretical analysis that our extracted features exhibit invariance to translation, near invariance to rotation, equivariance to permutation of atomic coordinates, and invariance under time reparametrization of paths. Finally, we target our task on predicting functional properties of proteins from MD data, a fundamental task for developing novel drug therapies. We consider three benchmarks including gene regulatory dynamics, epidermal growth factor receptor (EGFR) mutation dynamics, and G protein-coupled receptors (GPCR) dynamics for performance evaluation. In summary, we make the following three key contributions:

- •
We develop Deep signature, the first computationally efficient framework that characterizes the complex interatomic interactive dynamics of large-scale molecules.
- •
We theoretically demonstrate that our approach preserves symmetry under several geometrical transforms of atomic coordinates in 3D conformational space. Additionally, our method remains invariant under time reparameterization.
- •
We provide empirical results to show that our Deep Signature model achieves superior performance compared to other baseline methods on gene regulatory dynamics, EGFR mutation dynamics and GPCR dynamics.

## 2 Related Works

Molecular Representation. Encoding essential structural characteristics and biochemical properties into molecular representations is a long-standing research field in molecular biology, with wide applications in various drug discovery processes including virtual screening, similarity-based compound searches, target molecule ranking, drug toxicity prediction, and etc. (Li et al., ). One of the most widely used categories of molecular representation is two-dimensional fingerprints that extract the substructure, topological routes and circles solely from molecular connection tables. Owing to their ease of generation and usage, these 2D fingerprints are still extensively utilized as input for machine learning algorithms in modern drug discovery applications (Gao et al., ). In recent years, there has been a surge in the development of 3D structure-based molecular representations for their much more fine-grained characterization of interatomic interplay in 3D conformational space. Surface modeling-based approaches leverage mathematical models to restore protein surfaces and encode geometrical and chemical features of surface regions around binding sites into representations, exhibiting great potential in automatic drug discovery (Gainza et al., ; Zhu et al., ). In addition to surface-based representations, significant efforts have been devoted to learning accurate deep learning force fields for accelerating MD simulation (Schütt et al., ; Batatia et al., ; Batzner et al., ). While these methods can be adapted for molecular property prediction, they are limited to generating representation for a static frame, thus not suitable for our tasks where interframe interactions along timescales are crucial for understanding protein function.

The current investigation into characterizing molecular dynamics, especially interatomic interaction dynamics, remains limited, with only a few studies close to ours. Among them, Endo et al. (), Yasuda et al. () and Mustali et al. () are unsupervised methods that build local dynamics ensembles for manually selected particles and inspect each particle’s contribution to the dynamics independently. Li et al. () proposes to convert MD conformations into images and employ convolutional neural network to identify diverse active states from MD trajectories. Nevertheless, these approaches cannot recognize subtle interatomic interactions among atomic pathways for dissecting protein function, nor is the complexity of their representations independent of time stamps.

Coarse Graining (CG). CG is a widely adopted technique with the objective of preserving the crucial characteristics and dynamics inherent to a molecular system. This is achieved by grouping sets of atoms into CG beads, thereby enabling high-throughput MD simulations over larger time and length scales. Existing CG methods can be broadly categorized into two types: chemical and physical intuition-based approaches and machine learning-based approaches. The methods of first type construct the CG mapping by incorporating various biochemical properties derived from expert knowledge, for example, mapping each elaborately constructed cluster of four heavy (nonhydrogen) atoms into a single CG bead (Marrink & Tieleman, ) or simply assigning one CG bead centered at the $\alpha$-carbon for each amino acid (Ingólfsson et al., ). Machine learning-based methods can rapidly learn accurate potential energy functions for reduced structures of MD by training on large databases. Recent advancements in this area include multiscale coarse graining that optimizes to maximize a variational force matching score (Wang & Gómez-Bombarelli, ), relative entropy minimization (Thaler et al., ), and spectral graph approaches that account for structural typologies of proteins (Webb et al., ; Li et al., ). However, their transferability to unseen molecules remains suspicious, and the representation capability for complex macromolecular systems without increasing dimension and complexity is still underexplored (Khot et al., ).

## 3 Methodology

Figure: Figure 1: An overview of our proposed Deep Signature method.
Refer to caption: x1.png

### 3.1 Problem Formulation

Consider various molecular systems $\smash{S^{(k)}}$ that are distinct in their molecular behaviors. An MD simulation trajectory on $\smash{S^{(k)}}$ provides the trajectories for all
$\smash{N_{k}}$
atoms constituting the molecules, and can be represented as a sequence of snapshots
$\mathbf{X}_{1:T_{k}}^{(k)}=\{\mathbf{X}^{(k)}_{1},\mathbf{X}^{(k)}_{2},\dots,
\mathbf{X}^{(k)}_{T_{k}}\}$
. Here,
$\mathbf{X}^{(k)}_{t}\in\mathbb{R}^{N_{k}\times 3}$
indicates the atomic positions in 3D conformational space at time step $t$ for $t\in\{1,\dots,T_{k}\}$ and $\smash{T_{k}}$ is the total number of frames. To describe the structure of molecules within $\smash{S^{(k)}}$, we define a molecular graph $\smash{\mathcal{G}^{(k)}=\{\mathcal{E}^{(k)},\mathcal{V}^{(k)}\}}$, $\smash{|\mathcal{V}^{(k)}|=N_{k}}$, where $\smash{\mathcal{V}^{(k)}}$ is the node set corresponding to the atoms and $\smash{\mathcal{E}^{(k)}}$ is the edge set corresponding to chemical bonds. The adjacency matrix of $\smash{\mathcal{G}^{(k)}}$ is represented by $\smash{\mathbf{A}^{(k)}\in\mathbb{R}^{N_{k}\times N_{k}}}$, with
$\smash{\mathbf{A}^{(k)}_{i,j}}=1$
if $\smash{v_{i},v_{j}\in\mathcal{V}^{(k)}}$ and $\smash{(v_{i},v_{j})\in\mathcal{E}^{(k)}}$. For the molecular property prediction task, we have access to MD trajectories from $K$ molecular systems, each endowed with property labels denoted as
$\smash{\{(\mathbf{X}_{1:T_{k}}^{(k)},y^{(k)})\}_{k=1}^{K}}$
. The objective is to train algorithms to accurately predict the property label when provided with an previously unseen MD trajectory
$\smash{\mathbf{X}_{1:T_{k}}}$
.

### 3.2 Deep Signature

Our proposed method, referred to as Deep Signature, consists of a deep spectral clustering module that uses GNNs to extract coarse grained dynamics from raw MD trajectories, a path signature transform module that collects iterated integrals to characterize interatomic interactions along pathways, and a classifier to enable property prediction. The overall architecture is illustrated in Fig. 1.

Deep spectral clustering with GNNs. Given the MD trajectory $\smash{\mathbf{X}_{1:T}}$ and molecular graph $\mathcal{G}$ for a molecular system $S$, we start with extracting the reduced trajectory
$\smash{\tilde{\mathbf{X}}_{1:T}}$
by coarse gaining $\smash{\mathbf{X}_{1:T}}$ using deep spectral clustering. Specifically, we first obtain node representations via GNN layers as

$$ $\mathbf{H}^{l}=\sigma(\tilde{\mathbf{D}}^{-1/2}\tilde{\mathbf{A}}\tilde{ \mathbf{D}}^{-1/2}\mathbf{H}^{l-1}\mathbf{W}^{l-1}_{\mathrm{GNN}}),$ (1) $$

where $\smash{\mathbf{H}^{l}}$ denotes the node feature matrix at the $l$-th layer, $\smash{\mathbf{H}^{0}=\mathbf{X}_{1:T}}$, $\smash{\tilde{\mathbf{A}}=\mathbf{A}+\mathbf{I}}$ is the adjacency matrix $\mathbf{A}$ plus the identity matrix $\mathbf{I}$, $\smash{\tilde{\mathbf{D}}}$ is the degree matrix of $\smash{\tilde{\mathbf{A}}}$,
$\smash{\mathbf{W}^{l-1}_{\mathrm{GNN}}}$
are the learnable parameters of GNNs, and $\sigma$ is a nonlinear activation function. After then, we compute the cluster assignment matrix $\mathbf{Q}$ for the nodes using a multi-layer perceptron (MLP) with softmax on the output layer

$$ $\mathbf{Q}=\textit{Softmax}(\mathbf{W}_{\mathrm{MLP}}\mathbf{H}^{l}+\mathbf{b}),$ (2) $$

where $\mathbf{W}_{\mathrm{MLP}}$ and $\mathbf{b}$ are trainable parameters of the MLP. For the assignment matrix $\mathbf{Q}\in\mathbb{R}^{N\times M}$, where $M\ll N$ specifies the number of clusters, each row of it represents the node’s probability of belonging to a particular cluster. We employ the normalized-cut relaxation (Bianchi et al., ) as our clustering objective for minimization

$$ $\mathcal{L}_{u}=-\frac{Tr(\mathbf{Q}^{T}\tilde{\mathbf{A}}\mathbf{Q})}{Tr( \mathbf{Q}^{T}\tilde{\mathbf{D}}\mathbf{Q})}+\left\|\frac{\mathbf{Q}^{T} \mathbf{Q}}{||\mathbf{Q}^{T}\mathbf{Q}||_{F}}-\frac{\mathbf{I}_{M}}{\sqrt{M}} \right\|_{F},$ (3) $$

where the first term promotes strongly connected components to be clustered together, while the second term encourages the cluster assignments to be orthogonal and have similar sizes.

Upon leveraging $\mathbf{Q}$ from Eq. (2) for clustering, the corresponding reduced feature embedding matrix $\smash{\mathbf{H}^{\prime}}$ and adjacency matrix $\smash{\mathbf{A}^{\prime}}$ can be derived as follows

$$ $\mathbf{H}^{\prime}=\mathbf{Q}^{T}\mathbf{H};\quad\mathbf{A}^{\prime}=\mathbf{ Q}^{T}\tilde{\mathbf{A}}\mathbf{Q}.$ (4) $$

Since our model takes the sequence $\smash{\mathbf{X}_{1:T}}$ as input, $\smash{\mathbf{H}^{\prime}}$ inherently maintains the temporal order in the form of $\smash{\mathbf{H}^{\prime}_{1:T}}$. We utilize another MLP with the parameters $\smash{\mathbf{W}^{\prime}_{\mathrm{MLP}}}$ and $\mathbf{b}$ to map $\smash{\mathbf{H}^{\prime}_{1:T}}$ back into a reduced conformational space with the resulting dynamics and ground truth dynamics expressed as

$$ $\tilde{\mathbf{X}}_{1:T}=\mathbf{W}^{\prime}_{\mathrm{MLP}}\mathbf{H}^{\prime} _{1:T}+\mathbf{b}^{\prime};\quad\tilde{\mathbf{X}}_{1:T}^{\mathrm{pool}}= \mathbf{Q}^{T}\mathbf{X}_{1:T}.$ (5) $$

To ensure the fidelity of the coarse grained dynamics towards the original high-dimensional system, we further introduce a temporal consistency constraint, defined through a mean absolute error loss function with the form

$$ $\mathcal{L}_{t}=\frac{1}{T}\sum_{i=1}^{T}\big{|}\tilde{\mathbf{X}}_{i}-\tilde{ \mathbf{X}}_{i}^{\mathrm{pool}}\big{|}.$ (6) $$

Path signature transform. We now adopt the path signature method to characterize the interatomic temporal interactions among the coarse grained dynamics $\tilde{\mathbf{X}}_{1:T}\in\mathbb{R}^{T\times 3M}$. The basic idea of path signature is that, for a multidimensional continuous path, we can construct an ordered set consisting of all possible path integrals and combinations involving the path integrals among various individual dimensions as a comprehensive representation for this path (Lyons, ). Striving for a more precise definition, consider our coarse grained trajectory $\tilde{\mathbf{X}}_{1:T}$ with $(\tilde{\mathbf{X}}^{1}_{t},\tilde{\mathbf{X}}^{2}_{t},\dots,\tilde{\mathbf{X}
}^{3M}_{t})$ for $t\in\{1,\dots,T\}$, let us define $\widehat{\mathbf{X}}:[1,T]\rightarrow\mathbb{R}^{3M}$ as a piecewise linear interpolation of $\tilde{\mathbf{X}}_{1:T}$ such that $\widehat{\mathbf{X}}_{t}=\tilde{\mathbf{X}}_{t}$ for any $t\in\{1,\dots,T\}$, and a sub-time interval $[r_{i},r_{i+1}]$ corresponding to a time partition of $[1,T]$ with $\smash{1=r_{1}<r_{2}<\cdots<r_{\tau}=T}$. The depth-$D$ signature transform of $\smash{\tilde{\mathbf{X}}}$ over the interval $\smash{[r_{i},r_{i+1}]}$ is the vector defined as

$$ $\mathrm{Sig}^{D}_{r_{i},r_{i+1}}(\tilde{\mathbf{X}})=\left(1,\left\{S(\tilde{ \mathbf{X}})_{r_{i},r_{i+1}}^{j}\right\}_{j=1}^{3M},\ldots,\left\{S(\tilde{ \mathbf{X}})_{r_{i},r_{i+1}}^{j_{1},\ldots,j_{d}}\right\}_{j_{1},\ldots,j_{d}= 1}^{3M}\right),$ (7) $$

where for any $\smash{(j_{1},\ldots,j_{d})\in\{1,\dots,3M\}^{D}}$,

$$ $S(\tilde{\mathbf{X}})_{r_{i},r_{i+1}}^{j_{1},\ldots,j_{d}}=\underset{1<t_{1}< \ldots<t_{d}<T}{\int\dots\int}d\widehat{\mathbf{X}}_{t_{1}}^{j_{1}}\ldots d \widehat{\mathbf{X}}_{t_{d}}^{j_{d}}.$ (8) $$

After that, we take the logarithm of the signature transform features presented in Eq. (7) to eliminate redundant elements according to the shuffle product identity (Lyons et al., ) and acquire some minimal collection of the stream over a time interval. Specifically, given the vector space described by depth-$D$ signature in formal power series form as

$$ $\mathrm{Sig}^{D}_{r_{i},r_{i+1}}(\tilde{\mathbf{X}})=\sum^{D}_{d=0}\sum_{j_{1} ,\ldots,j_{d}\in\{1,\dots,3M\}}S(\tilde{\mathbf{X}})_{1,T}^{j_{1},\ldots,j_{d} }e_{j_{1}}\dots e_{j_{d}},$ (9) $$

The depth-$D$ log-signature transform corresponds to taking the formal logarithm of Eq. (9), which can be expressed as

$$ $\mathrm{LogSig}^{D}_{r_{i},r_{i+1}}(\tilde{\mathbf{X}})=\sum_{n=1}^{3M}S( \tilde{\mathbf{X}})^{n}_{r_{i},r_{i+1}}e_{n}+\sum_{1\leq n<m\leq 3M}\frac{1}{2 }(S(\tilde{\mathbf{X}})^{n,m}_{r_{i},r_{i+1}}-S(\tilde{\mathbf{X}}))^{m,n}_{r_ {i},r_{i+1}}[e_{n},e_{m}]+\dots$ (10) $$

where $[e_{n},e_{m}]$ is the Lie bracket (Reizenstein, ) defined by $[e_{n},e_{m}]=e_{n}\times e_{m}-e_{m}\times e_{n}$. Generally, log-signature possesses a lower dimension compared to the original signature feature while carrying exactly the same information. Through this transform, the interatomic dynamic correlations in 3D space are embedded as geometric areas into features. The resulting signature sequence, denoted as $(\mathrm{LogSig}^{D}_{r_{1},r_{2}}(\tilde{\mathbf{X}}),\mathrm{LogSig}^{D}_{r_
{2},r_{3}}(\tilde{\mathbf{X}}),\dots)$, can be regarded as a discretization of dynamic data. In this work, we further utilize LSTM to tackle the nonlinear interactions between them.

Classifier. Finally, we introduce a classifier $f$ implemented as a two-layer MLP which outputs the predicted label $\hat{y}=f(\mathrm{LogSig}^{d}_{1,T}(\tilde{\mathbf{X}}))$ for molecular property prediction. To establish the classification objective function, we leverage cross-entropy loss defined as follows

$$ $\mathcal{L}_{c}=-y\log(\hat{y})-(1-y)\log(1-\hat{y}),$ (11) $$

where $y$ is the ground truth label for the system and $\hat{y}$ is the predicted label generated by our model.

By combining the loss terms from Eqs. (3), (6) and (11), we arrive at the overall loss function for training our Deep Signature model as

$$ $\mathcal{L}=\lambda_{1}\mathcal{L}_{u}+\lambda_{2}\mathcal{L}_{t}+\lambda_{3} \mathcal{L}_{c}.$ (12) $$

Here, $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are scaling parameters to balance the contributions of these three loss terms.

It is worth noting that, recent research has delved into the integration of path signatures with GNNs for spatial-temporal modeling, and marked significant strides in traffic forecasting scenarios (Choi & Park, ; Wang & Zang, ). These approaches utilize path signatures to process temporal information for each node individually at a low level, while employing GNNs empowered neural differential equation at a higher level to ensure the preservation of spatial topological consistency during dynamic forecasting. Nonetheless, they fall short in characterizing the temporal interactions among different nodes, making them unsuitable for molecular property prediction tasks where the interactions among dominant kinetic pathways play an essential role in protein function (Araki et al., ; Mustali et al., ).

### 3.3 Equivariance & invariance of log-signature

Physical properties often exhibit well-defined symmetry characteristics, and integrating such equivariance into our learned feature space can enhance interpretability and mitigate learning difficulty. In the subsequent analysis, we rigorously examine the equivariance of our deep signature in relation to various geometric transformations including translation, rotation, and permutation on atomic coordinates, as well as reparameterization over time. The complete proofs are provided in Appendix A.

Translation invariance.
The coarse-grained dynamics acquired by a linear mapping in Eq. (5) demonstrate that $\smash{\tilde{\mathbf{X}}_{1:T}}$ maintain equivariance with respect to translations on input trajectory $\smash{\mathbf{X}_{1:T}}$. Besides, since the path signature is composed of iterated integrals, it inherits the property of translation invariance. As a result, log-signature features are translation-invariant for the input trajectory.

Rotation invariance. When a rotation is applied to the 3D conformational space, an equivalent rotation exists in the coarse grained conformational space,
indicating the equivariance of coarse grained dynamics $\smash{\tilde{\mathbf{X}}_{1:T}}$ to rotations of the input trajectory $\smash{\mathbf{X}_{1:T}}$.
Furthermore, while research on rotation invariance of path signatures for higher depth and multi-dimensional paths is still limited (Diehl, ), we demonstrate that the majority of elements constituting our Deep Signature features are rotation-invariant, particularly in large-scale systems like the molecular systems studied in this work.

Permutation equivariance. GCNs that aggregate contributions from neighboring atoms are invariant to permutations of those atoms, which implies that our coarse-grained dynamics are equivariant to permutations of atom indices. Additionally, the multi-indices used to index the iterated integrals in a signature are sorted in ascending order during implementation, ensuring that our signature features remain equivariant to permutations. This property is essential for facilitating the interpretability of our method, as it allows us to trace the dominant kinetic pathways that contribute to protein function.

Time-reparametrization invariance. Path signature possesses a powerful property that it remains invariant under time-reparameterization of the underlying stream (Lyons, ). This property substantially decreases the complexity of certain challenges by eliminating dependence on the sampling rate while preserving all other relevant information within the stream. Furthermore, it enhances the robustness of our method against deviations in molecular dynamics occurring over time scales.

### 3.4 Implementation

Figure: Figure 2: The architecture of deep spectral clustering module.
Refer to caption: x2.png

Architecture. Our introduced framework consists of three fundamental modules: a deep spectral clustering module, a path signature transform module and a classifier. We adopt a hierarchical pooling architecture to implement the deep spectral clustering module, which is schematically depicted in Fig. 2. As shown, it contains a stack of $L$ graph pooling layers with each layer consisting of a GCN layer for obtaining node embeddings and an MLP layer for cluster assignment, such that it would gradually coarsen the dynamics from the atomic level. The last GCN layer and MLP take the node feature matrix and adjacency matrix processed by the $L$-th graph pooling layer as input and output the final coarse grained dynamics $\tilde{\mathbf{X}}_{1:T}$.

Figure: Figure 3: The architecture of path signature transform module.
Refer to caption: x3.png

The implementation of the path signature transform module is illustrated in Fig. 3. It takes $\tilde{\mathbf{X}}_{1:T}$ as input, which is first partitioned into equally spaced intervals. The resulting segments are subsequently processed through the $\mathrm{LogSig}(\cdot)$ function to extract log-signature features.
This is implemented based on Signatory(^1^11[https://github.com/patrick-kidger/signatory](https://github.com/patrick-kidger/signatory)), a Python package that facilitates differentiable computations of the signature and log-signature transforms on both CPU and GPU. The resulting signature sequence is then fed into an LSTM layer to capture their interactions and followed by layer normalization (Ba et al., ) to standardize the feature values. For the classifier, we utilize a two-layer MLP together with a sigmoid activation function to enable property prediction. It is noteworthy that the overall architecture is fully differentiable, thus we can leverage off-the-shelf optimization techniques to train our model in an end-to-end fashion.

Cross validation and independent test. We employ a five-fold cross-validation strategy for training our model. In contrast to conventional random division, we take the temporal nature of the trajectory data into consideration for data partition. Specifically, each trajectory is divided into 5 groups with the same time interval according to its temporal order. Subsequently, we further partition the data within each group into five folds. The validation set is constructed by selecting one fold from each group and is employed for model selection, while the remaining four folds within each group are gathered to form the training set. This process is repeated five times sequentially, resulting in the creation of the five-fold cross-validation dataset. Moreover, for each running, we evaluate the prediction accuracy of our method on an independent unseen test set and report the averaged results.

## 4 Experiments

In this section, we conduct experiments on representative molecular dynamic systems to evaluate the effectiveness of our proposed method. We begin with a synthetic dataset that reports gene regulatory dynamics (Gao et al., ). We then assess the performance on two large-scale MD simulation datasets, including epidermal growth factor receptor (EGFR) mutant dynamics (Zhu et al., ) and G protein-coupled receptors (GPCR) dynamics (Rodríguez-Espigares et al., ). More details on dataset construction can be found in Appendix B.

### 4.1 Experiments on Gene regulatory dynamics

For gene regulatory dynamics, we generate 100 trajectories that describe the interactive dynamics between genes and transcription factors. These dynamics are categorized into degradation type or dimerization type, with an equal number of trajectories for each type. The systems are simulated over a period of $2$s and with a time interval of $0.004$s, resulting in 500 frames per trajectory. Besides, each system encompasses 100 nodes with randomly generated Power-law network structure to describe their relationships. We implement the deep spectral clustering module in our approach using two graph pooling layers that first coarsen the dynamics into 60 nodes and then into 30 nodes. The coefficients of loss terms are set as $\lambda_{1}=1$, $\lambda_{2}=0.01$, and $\lambda_{3}=10$. We optimize our model using Adam with an initial learning rate of 5e-4 and a weight decay of 1e-4. For comparison, we consider several baseline approaches which aggregate nodes’ dynamics by averaging without taking into account their interactions. These methods incorporate a deep spectral clustering model with an MLP to process the first frame only (Head), the last frame only (Tail), the first and last frame (Head & Tail), and all frames with an additional LSTM layer to aggregate nodes’ temporal dynamics (GraphLSTM). The overall comparison results with baseline methods are presented in Table 4.

Figure: Table 1: Comparisons of classification performance on gene regulatory dynamics. Results are averaged over 5 runs.
Refer to caption: x4.png

From the table, we observe that our Deep Signature method achieves the best performance with a clear margin over baseline approaches. Besides, compared to GraphLSTM which aggregates atomic dynamics by averaging them directly, our Deep Signature exhibits a substantial performance boost by incorporating the log-signature transform into our method design, demonstrating the importance of capturing interatomic interaction dynamics for property prediction. Moreover, we validate the scalability of our approach to dynamical system across various time scales and different system sizes. As shown in Fig. 4(a), when extending the time duration of dynamics from $1$s to $5$s, our Deep Signature consistently obtains the best performance across these settings. Furthermore, when we change the system size via increasing the number of nodes that would enhance the complexity of dynamics, the results in Fig. 4(b) show that our method can still perform well under these conditions. These results showcase the good scalability of our method for dealing with dynamical systems of different timescales and sizes. In contrast, while GraphLSTM yields competitive results for systems with fewer than 250 nodes, a notable performance degradation appears as the number of nodes keeps increasing. We conjecture that this is because the averaged dynamics become less representative of the global dynamics as the system scale grows.

### 4.2 Experiments on epidermal growth factor receptor mutation dynamics

**Table 2: Comparison of different methods for the classification results on the EGFR dynamics.**
| Method | Dihedral angle | C$\bm{\alpha}$-dihedral angle | Head | Tail | Head & Tail | GraphLSTM | Deep Signature |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Accuracy | 67.467 | 61.333 | 61.467 | 59.733 | 59.467 | 65.600 | 69.333 |
| Recall | 0.024 | 0.040 | 0.024 | 0.028 | 0.060 | 0.084 | 0.220 |

Figure: Table 3: Ablation study on different loss items for Deep Signature. We report accuracy and recall for performance evaluation.
Refer to caption: x6.png

In this study, we investigate the binding process between EGFR mutations and EGFR-receptor tyrosine kinases (RTK) during which their interactions can contribute to drug resistance mechanisms. The combinations between four RTK partners and five mutation types together with wild type as a reference are considered. To acquire the dynamic data, we follow the pipeline presented in (Zhu et al., ) that each system first undergoes energy minimization prior to simulation, then it is heated for 100 ps, followed by density equilibration for 100 ps and constant pressure equilibration for 5 ns, and finally runs simulation on the equilibrated structures for 50 ns. This results in 24 trajectories, each comprising 1,000 frames describing the temporal interactions among approximately 5,000 atoms. Each trajectory is labeled according to its sensitivity towards the drug. Here, our deep spectral clustering model consists of three layers that progressively coarsen the dynamics into 400, 200, 50 nodes. We set the scaling parameters as $\lambda_{1}=1$, $\lambda_{2}=0.01$, and $\lambda_{3}=10$. The model is trained for 200 epochs with an initial learning rate 5e-5 and a weight decay of 1e-5. In addition to the deep learning-based baselines introduced above, we further consider two baselines that extract dihedral angles as geometrical descriptors, followed by PCA to reduce the dimension of features and SVM for classification. The overall comparison results are summarized in Table 2.

As can be seen, our Deep Signature method consistently obtains the best classification results on EGFR dynamics, achieving an accuracy of 69.333% and a recall ratio of 0.220. This validates the practical applicability of our method for its capability to tackle complex dynamics inherent in large-scale molecular systems. In addition, in contrast to baseline methods that achieve high accuracy but low recall which often due to their reliance on bypasses caused by the varying sample sizes for different classes, our method demonstrates a substantial improvement in recall, indicating its robustness under class imbalance issues. Furthermore, when compared to conventional geometrical descriptor methods, deep learning-based approaches, with the exception of ours, tend to underperform on this dataset. This discrepancy can be attributed to the intricate non-linearity and atomic fluctuations characteristic of EGFR dynamics, which pose significant challenges for deep learning techniques to extract meaningful patterns from a limited amount of data.

Ablation study. We conduct an ablation study to assess the necessity of each loss term used in our model, with the results for EGFR dynamics presented in Table 5. We begin by training our Deep Signature model using only the classification loss $\mathcal{L}_{c}$, which yields an accuracy of $67.6\%$. Upon incorporating the temporal consistency loss $\mathcal{L}_{t}$, the accuracy increases to $67.967\%$. However, when we replace $\mathcal{L}_{t}$ with the spectral clustering loss $\mathcal{L}_{u}$, performance slightly degrades to $66\%$. This decrease may be attributed to the tendency of $\mathcal{L}_{u}$ to guide the spectral clustering module towards hard assignments, thereby impairing the representation ability of the coarse grained dynamics. Training with all three loss terms produces the best results, demonstrating the validity of our method design.

Figure: Figure 6: The effect of the deep spectral clustering module for C797S and G719S mutated EGFR dimers. (a) are origianl 3D conformations, (b) are coarse grained graphs, (c) are the RMSD curves of raw trajectories and coarse grained trajectories.
Refer to caption: x7.png

Analysis on deep spectral clustering module. To understand the impact of our spectral clustering module on the graph structures and trajectories for EGFR dynamics, we investigate two cases: a drug-resistant C797S mutant and a drug-sensitive G719S mutant respectively dimerized with EGFR partners, with their results visualized in Fig. 6. By comparing (a) and (b), we observe that coarse grained graphs produced by our module can maintain the overall structure of original conformations, and the nodes mainly distribute in the middle region, which is valid since interatomic interactions that facilitate protein function commonly occur in this area. Besides, as shown in (c), the root mean square deviation (RMSD) curves for both the raw and coarse grained trajectories generally follow the same trend, indicating a high fidelity of the coarse grained dynamics to the original dynamics.

Analysis on path signature transform module. We visualize the feature space via t-SNE to explore the discriminability of learned features after path signature transform. The results of all data samples with colors indicating their classes are presented in Fig. 5. As seen, both the non-learnable features, which compute dihedral angles for all C$\alpha$ atoms, and the learned features derived solely from the first frame tend to gather tightly in feature space, making it difficult to establish a decision boundary to distinguish samples belonging to different classes. In contrast, features extracted via our Deep Signature method are more uniformly distributed, and we can find such a decision boundary easily, verifying the good discriminability of these features.
We further use green circles to indicate training samples and orange dashed circles to denote test samples in Fig. 5. These two circles typically exhibit limited overlap due to distribution shift between the training and test samples, which arises from heterogeneity in the number of atoms and molecular topological structures. Nevertheless, Deep Signature exhibits an ability to learn generalizable features for unseen samples. This is because learning features that respect symmetry enhances generalization capability (Schütt et al., ), and incorporating layer normalization helps mitigate the distribution discrepancy (Ba et al., ).

Figure: Figure 7: The visualization of critical pathways and interatomic interactions that contribute to drug sensitiveness on the EGFR dynamics.
Refer to caption: x8.png

Analysis on interpretability. Our approach also exhibits notable interpretability for its ability to identify every possible type of interaction among atoms that are essential for protein functioning. This capability stems from the deep spectral clustering module which only involves linear mapping for cluster assignment, and log-signature features that maintain permutation equivariance as discussed in Section 3.3. In practice, we apply the Gradient $\odot$ Input method (Shrikumar et al., ) to quantify the contribution of each element in log-signature feature to the final prediction output and then identify three key atoms whose interactive dynamics play a pivotal role in the progression of EGFR. The identified atoms, along with their interatomic distances, are illustrated in Fig. 7. Notably, these identified atoms fall within the hinge region that comprises all ATP binding sites, demonstrating strong consistency with the experimental observation (Kaufman et al., ).

### 4.3 Experiments on G protein-coupled receptors dynamics

**Table 4: Comparison of different methods for the classification results on the GPCR dynamics.**
| Method | Dihedral angle | C$\bm{\alpha}$-dihedral angle | Head | Tail | Head & Tail | GraphLSTM | Deep Signature |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Accuracy | 44.000 | 50.000 | 46.667 | 45.467 | 48.533 | 53.333 | 64.200 |
| Recall | 0.355 | 0.667 | 0.366 | 0.366 | 0.385 | 0.360 | 0.413 |

The GPCR superfamily is a major therapeutic target as its functioning regulates nearly every physiological process in the human body. To create a dataset for its analysis, we select 26 structures with their MD simulations from the Molecular Dynamics Database for GPCRs (GPCRmd) (Rodríguez-Espigares et al., ), including 13 active state structures and 13 inactive state structures. The task is to identify these distinct active states. Each simulation is conducted for 500 ns with a time interval of 200 ps, therefore each trajectory consists of 2,500 frames that describe the conformational dynamics. Here, we only consider C$\alpha$ atoms. Experimental results are provided in Table 4.

It is evident that Deep Signature achieves the highest performance on the GPCR dynamics, with an accuracy of 64.200% and a recall ratio of 0.413. The results are significantly better than the compared baselines. Traditional geometrical descriptor methods, such as the dihedral angles and C$\alpha$-dihedral angles, exhibit relatively lower classification accuracy and recall compared to our method. This again verifies that geometrical descriptors are inadequate for capturing dynamic patterns when the dynamics are highly nonlinear due to the existence of atomic fluctuations in GPCR dynamics.
Moreover, relying solely on the first or last frame fails to account for the dynamic nature of interactions, thereby limiting the ability to discern useful patterns for predicting structural states. Although the GraphLSTM method outperforms the geometric descriptor methods, achieving a moderate accuracy of 53.333% and a recall ratio of 0.360, it still falls short of the superior classification capabilities exhibited by the Deep Signature method, demonstrating the effectiveness of our approach in characterizing the complex dynamics of GPCR systems.

## 5 Conclusion

In this paper, we introduce a novel deep learning framework, Deep Signature, to deal with the long-standing challenge of understanding protein dynamics in large-scale biological systems. It comprises soft spectral clustering to aggregate cooperative dynamics and log-signature transformation to characterize global interactive dynamics. Theoretically, our method exhibits desirable properties such as invariance to translation, near-invariance to rotation, and equivariance to atomic coordinate permutation. Experimental results on three biological process benchmarks verify the effectiveness of Deep Signature in capturing the complex interactive dynamics of large-scale molecular systems. We hope our work can offer a promising new direction for the analysis of protein dynamics.

## Appendix A Proof for Properties

In the beginning of our proof, we would first decompose our feature extracting process into two operations $g_{\mathrm{GNN}}(\cdot)$ and $\mathrm{LogSig}(\cdot)$ as presented in Section 3.2 such that the overall feature transform is their composition $g_{\mathrm{GNN}}\circ\mathrm{LogSig}$. This decomposition offers us the opportunity to analyze the properties of component separately. Besides, since the logarithm map is bijective, which implies there is one-to-one correspondence between the signature and the log-signature (Lyons et al., ), thus we can turn to analyze signature transform for intuitive demonstration.

### A.1 Translation Invariance

Given trajectory data $\mathbf{X}_{1:T}\in\mathbb{R}^{T\times N\times 3}$, let $\mathcal{T}_{B}$ represent a translation matrix $B\in\mathbb{R}^{N\times 3}$ on the trajectory data such that at a certain time stamp $t$ we have $\mathcal{T}_{B}(\mathbf{X}_{t})=\mathbf{X}_{t}+B$. For the coarse grained dynamics acquired by $\tilde{\mathbf{X}}_{t}^{\mathrm{pool}}=\mathbf{Q}^{T}\mathbf{X}_{t}$ as presented by Eq. (5), as $\mathbf{Q}^{T}$ is a linear matrix, we can get that $\mathbf{Q}^{T}\mathcal{T}_{B}(\mathbf{X}_{t})=\mathbf{Q}^{T}\mathbf{X}_{t}+
\mathbf{Q}^{T}B=\mathcal{T}_{\mathbf{Q}^{T}B}(\mathbf{Q}^{T}\mathbf{X}_{t})$, thus $\tilde{\mathbf{X}}_{1:T}$ maintain equivariance with respect to translation $\mathcal{T}_{B}$ on input trajectory $\mathbf{X}_{1:T}$. Besides, as presented in Eq. (8) that the path signature is composed of iterated path integrals, it inherits the properties of translation invariance that $\mathrm{Sig}^{D}_{a,b}(\tilde{\mathbf{X}}+c)=\mathrm{Sig}^{D}_{a,b}(\tilde{
\mathbf{X}})$. In the whole, our path signature features $\mathrm{Sig}^{D}_{a,b}(\tilde{\mathbf{X}})$ are translation-invariant with respect to input trajectory data $\mathbf{X}_{a,b}$.

### A.2 Rotation Invariance

Given the coarse grained dynamics acquired by $\tilde{\mathbf{X}}_{t}=\mathbf{Q}^{T}\mathbf{X}_{t}$, where $\mathbf{X}_{t}\in\mathbb{R}^{N\times 3}$ and $\tilde{\mathbf{X}}_{t}\in\mathbb{R}^{M\times 3}$. The cluster assignment matrix $\mathbf{Q}\in\mathbb{R}^{N\times M}$ holds that $\mathbf{Q}\mathbf{1}_{M}=\mathbf{1}_{N}$ and $\mathbf{Q}^{T}\mathbf{Q}=\mathbf{I}_{M}$. For any orthogonal matrix $\mathbf{R}_{\theta}\in\mathbb{R}^{3\times 3}$ performed on $\mathbf{X}_{t}$ termed as $\mathbf{R}_{\theta}\mathbf{X}_{t}=(\mathbf{R}_{\theta}\mathbf{X}_{t}^{1},\dots
,\mathbf{R}_{\theta}\mathbf{X}_{t}^{N})$, we proof below that rotating the input results in an equivalent rotation of the output $\mathbf{R}_{\phi}\tilde{\mathbf{X}}_{t}=\mathbf{Q}^{T}\mathbf{R}_{\theta}
\mathbf{X}_{t}$, where $\mathbf{R}_{\phi}\in\mathbb{R}^{M\times 3}$ is an orthogonal matrix performed on $\tilde{\mathbf{X}}_{t}$. Let we start with the right-hand side $\mathbf{Q}^{T}\mathbf{R}_{\theta}\mathbf{X}_{t}$. Here, $\mathbf{Q}^{T}$ has the form $\mathbf{Q}^{T}=(\mathbf{Q}^{T}_{1},\dots,\mathbf{Q}^{T}_{M})$, where the assignment function for $m$-th cluster is $\mathbf{Q}^{T}_{m}\in\mathbb{R}^{N}$, $m\in\{1,\dots,M\}$, and the corresponding aggregated coordinate is $\mathbf{Q}^{T}_{m}\mathbf{X}_{t}$ which is a linear combination in each dimension. Due to the fact that $\mathbf{Q}^{T}_{m}\mathbf{R}_{\theta}\mathbf{X}_{t}$ can always be represented by $\mathbf{R}_{\phi}\mathbf{Q}^{T}_{m}\mathbf{X}_{t}$, which implies that for a rotation transform $\mathbf{R}_{\theta}$ on the original conformational space, there always exists an equivalent rotation transformation $\mathbf{R}_{\phi}$ in coarse grained conformational space. We then accomplish our proof that $\mathbf{R}_{\phi}\tilde{\mathbf{X}}_{t}=\mathbf{Q}^{T}\mathbf{R}_{\theta}
\mathbf{X}_{t}$.

After that, we turn to analyze the effect of rotation on features acquired by signature transform. Unfortunately, not all elements of a path signature feature exhibit rotation invariance. As demonstrated by Diehl (), rotation invariants only exist on levels of even order. For a 2-dimensional continuous path $X:[a,b]\rightarrow\mathbb{R}^{2}$ as presented in Fig. 8, the depth-1 terms corresponding to the variations for each channel over the interval are represented by $\Delta X_{1}$ and $\Delta X_{2}$, while the depth-2 term corresponding to the signed area between the chord connecting the endpoints and the real path is denoted as $A$. The calculation of depth-2 term corresponds to the coefficient of the polynomial $[e_{n},e_{m}]$ in Eq. (10)

$$ $A=A_{+}-A_{-}=\frac{1}{2}\left(S^{1,2}_{a,b}(X)-S^{2,1}_{a,b}(X)\right).$ $$

Here, the sign of area exactly indicates the orientation of acceleration, while an increase in the absolute value of surrounding accelerations would increase the proportion of the represented area accordingly. Obliviously, the area $A$ is rotation invariant, while $\Delta X_{1}$, $\Delta X_{2}$ do not hold this property. As suggested by Diehl (), a depth-2 term in Eq. (10) is rotation invariant only if it involves the iterated integrals over different channels. Although further investigation into rotation invariants for path signature with higher depth over beyond two dimensional paths is still lacking, current conclusion on rotation invariance acquired for depth-2 log-signauture can readily extend to higher dimensional paths.

In our experimental setup, given the coarse grained dynamics $\tilde{\mathbf{X}}_{a,b}\in\mathbb{R}^{T\times 3M}$, the dimension of depth-$D$ log-signature is

$$ $\mathrm{dim}(\mathrm{LogSig}^{D}_{a,b}(\tilde{\mathbf{X}}))=\sum_{d=1}^{D} \frac{1}{d}\sum_{i|d}\mu\left(\frac{d}{i}\right)(3M)^{i},$ $$

where $\mu$ is the M$\ddot{\mathrm{o}}$bius function defined as

$$ $\mu(n)=\left\{\begin{array}[]{rcl}&0&\mbox{if $n$ has one or more repeated prime factors}\\ &1&\mbox{if $n=1$}\\ &(-1)^{k}&\mbox{if $n$ is the product of $k$ distinct prime numbers}\end{array }\right.$ $$

When we specify $d=2$, the ratio of rotation-invariant elements over the whole log-signature feature can then be calculated as

$$ $\gamma=\frac{3M-1}{3M+1},$ $$

which indicates the majority of features are rotation-invariant for a large-scale molecular system as we focus on in this paper.

Figure: Figure 8: Geometric intuition for the first two levels of the log-signature on a 2-dimensional path. We can observe that the depth-1 terms represent the change in each of the coordinates over the interval, and the depth-2 term corresponds to the $L\acute{e}vy~{}area$ of the path, shown as the signed area enclosed by the curve and the chord connecting its start and endpoints.
Refer to caption: x9.png

### A.3 Permutation Equivariance

For the deep clustering module implemented by GNNs $g_{\mathrm{GNN}}(\cdot)$ with the formula presented in Eq. (1) as follows

$$ $\mathbf{H}^{l}=\sigma(\tilde{\mathbf{D}}^{-1/2}\tilde{\mathbf{A}}\tilde{ \mathbf{D}}^{-1/2}\mathbf{H}^{l-1}\mathbf{W}^{l-1}_{\mathrm{GNN}}).$ $$

For one layer, the calculation of node embeddings per node involves sum over contributions from different atoms. As a result, node embedding matrix $\mathbf{H}^{l}$ naturally exhibit equivariance with respect to the permutation symmetries of graphs. In addition, the linear mapping shown in Eq. (5) only works for dimension alignment, making no influence on the permutation equivariance. So the coarse grained dynamics $\tilde{\mathbf{X}}_{1:T}$ is permutation equivariant with respect to an input set of atoms.

We then delve into log-signature transform by examining Eq. (10). Notably, the monomials such as $e_{n}$ and $[e_{n},e_{m}]$ constitute the basis vectors of the vector space, and their coefficients are calculated using iterate integrals over paths indexed by these monomials. Ideally speaking, log-signature features hold permutation invariance since we can arbitrarily select and order the monomials. Such phenomenon becomes even clearer when we revisit the signature transform in Eq. (7). In the signature transform, the superscripts indicating the paths traverse the set of all multi-indexes denoted as

$$ $W=\{(j_{1},\dots,j_{d})|d\geq 1,j_{1},\dots,j_{d}\in\{1,\dots,3M\}\}.$ $$

Hence, the order of the multi-indexes has no impact on the captured dynamic information. In practical implementations, we organize the multi-indexes in ascending order to collect these iterated integral terms (Kidger & Lyons, ). Consequently, a permutation on the indices of input atoms would yield a predictable permutation transform on the indices of log-signature terms. As a result, we drive the conclusion that our Deep signature method exhibits permutation equivariance with respect to input atom indices.

### A.4 Time-reparametrization invariance

Since coarse graining will maintain temporal consistency, we only need to analyze the invariance under time-reparametrization for coarse grained dynamics. A reparametrization of a coarse grained pathway $\mathbf{X}:[a,b]\rightarrow\mathbb{R}^{3M}$ is a path $\breve{\mathbf{X}}:[a,b]\rightarrow\mathbb{R}^{3M}$ where $\breve{\mathbf{X}}_{t}=\mathbf{X}_{\psi_{t}}$ where $\psi$ is a subjective, continuous, non-decreasing function $\psi:[a,b]\rightarrow[a,b]$. We have the following theorem

###### Theorem 1

Let $\breve{\mathbf{X}}_{t}$ be a reparametrization of $\mathbf{X}$, then we have $\mathrm{Sig}(\breve{\mathbf{X}})=\mathrm{Sig}(\mathbf{X})$.

Consider two real-valued paths $X,Y:[a,b]\rightarrow\mathbb{R}$ and a surjective, continuous, non-decreasing reparametrization function over time $\psi:[a,b]\rightarrow[a,b]$. Then we have reparametrized paths $\breve{X},\breve{Y}:[a,b]\rightarrow\mathbb{R}$ by $\breve{X}_{t}=X_{\psi_{t}}$ and $\breve{Y}_{t}=Y_{\psi_{t}}$

$$ $\int_{1}^{T}\breve{Y}_{t}d\breve{X}_{t}=\int_{1}^{T}Y_{\psi_{t}}dX_{\psi_{t}}= \int_{1}^{T}Y_{\psi_{t}}\dot{X}_{\psi_{t}}\dot{\psi_{t}}dt=\int_{1}^{T}Y_{\psi _{t}}\dot{X}_{\psi_{t}}d\psi_{t}$ $$

After replacing $u=\psi_{t}$, we have $\int_{1}^{T}\breve{Y}_{t}d\breve{X}_{t}=\int_{1}^{T}Y_{u}dX_{u}$, which means path integrals are invariant under a time reparametrization of both paths.

Since every term of the signature $S(\mathbf{X})^{j_{1},\dots,j_{d}}_{a,b}$ is defined as an iterated path integral of $\mathbf{X}$, it follows from the above that

$$ $S(\breve{\mathbf{X}})^{j_{1},\dots,j_{d}}_{a,b}=S(\mathbf{X})^{j_{1},\dots,j_{ d}}_{a,b},\quad\forall k\geq 0,~{}j_{1},\dots,j_{d}\in\{1,\dots,3M\}$ $$

This complements the proof.

###### Theorem 1

## Appendix B Experimental Settings

### B.1 Gene regulatory dynamics

The dynamics for gene regulatory networks are governed by Michaelis-Menten equation as follows,

$$ $\frac{\mathrm{d}\bm{x}_{t}(v_{i})}{\mathrm{d}t}=-b_{i}\bm{x}(v_{i})f+\sum_{j=1 }^{n}\mathbf{A}^{(i,j)}\frac{\bm{x}^{h}(v_{j})}{\bm{x}^{h}(v_{j})+1},$ (13) $$

where the first term models the degradation when $f=1$ or dimerization when $f=2$, and the second term represents genetic activation, with the Hill coefficient $h$ determining the level of cooperation in the regulation of the gene.

Data generation. For the gene regulatory dynamics, we curate a dataset consisting of 100 trajectories that delineate the intricate interactive dynamics between genes and transcription factors. These trajectories are divided into two distinct categories: degradation type ($f=1$) and dimerization type ($f=2$), each encompassing an equal number of trajectories. To commence a simulation, we first initialize a graph network featuring 100 nodes using a Power-law network generator to elucidate the structural interconnections between these nodes. Subsequently, we employ the Dormand-Prince method to numerically solve the gene regulatory system described by Eq. (13), with a simulation duration of 2 seconds and a time interval of 0.004 seconds. This computation simulation yields trajectories comprising 500 frames each, laying the foundation for our subsequent experiments.

Model architecture. We implement the deep spectral clustering module in our approach using two graph pooling layers, which first coarsen the dynamics into 60 nodes and subsequently into 30 nodes. The dimension of hidden states for the stacked GCN layers is kept as 10. For a tiny implementation, we simply extract the log-signature features for the coarse grained dynamics and then employ a two-layer MLP for property prediction.

Training details. We train our model with mini-batches of size 48 for 100 epochs using Adam with the initial learning rate of 5e-4 and a weight decay 1e-4. The coefficients of loss terms are set as $\lambda_{1}=1$, $\lambda_{2}=0.01$, and $\lambda_{3}=10$.

### B.2 Epidermal growth factor receptor mutation dynamics

Data generation. In this investigation, we delve into the intricate binding dynamics between EGFR mutations and RTK, exploring how their interactions can trigger mechanisms of drug resistance. The study encompasses the amalgamation of four RTK partners and five distinct mutation types, alongside the wild type serving as a reference point. To capture the dynamic essence of these interactions, we adhere to the methodological framework outlined in (Zhu et al., ), wherein each system undergoes a meticulous pre-simulation process. This involves an initial energy minimization phase, followed by a 100 picosecond heating stage, subsequent density equilibration spanning 100 picoseconds, and a further 5 nanoseconds of constant pressure equilibration. The equilibrated structures then undergo a simulation period of 50 nanoseconds, resulting in the generation of 24 trajectories. Each trajectory comprises 1000 frames, detailing the temporal interplay among approximately 5,000 atoms. Notably, these trajectories are categorized based on their sensitivity or resistance to the administered drug. For the MD simulations, we leverage the explicit-solvent model integrated within the Amber software suite, utilizing the *Ff99SB* and *gaff* force fields to drive the simulations.

Model architecture. For the deep spectral clustering module, we utilize three graph pooling layers that progressively coarsen the dynamics into 400, 200, 50 nodes. The dimension of hidden states for the stacked GCN layers is kept as 20. For the path signature transform module, we partition the input coarse grained dynamics into four segments for the subsequent application of the log-signature transform.

Training details. The model is trained with mini-batches of size 16 for a total of 200 epochs. We employ the Adam optimizer with an initial learning rate of 5e-5 and a weight decay 1e-5 to optimize the model parameters. Additionally, the scaling parameters are specified as $\lambda_{1}=1$, $\lambda_{2}=0.01$, and $\lambda_{3}=10$.

### B.3 G protein-coupled receptors dynamics

Data generation. We download our data from the GPCRmd (http://gpcrmd.org/) (Rodríguez-Espigares et al., ) database, which is an online platform that incorporates web-based visualization capabilities and shares data. This database includes at least one representative structure from each of the four structurally characterized GPCR classes, and holds more than 600 GPCR MD simulations from GPCRmd community and individual contributions. To create our dataset, we select 26 trajectories of the $\beta$2AR-rh1 GPCR inactive (2RH1) and active (3P0G) receptor state with a full agonist. The receptor consists of 282 and 285 amino acids for inactive and active state respectively. Each simulation is conducted for 500 ns with a time interval of 200 ps, therefore every trajectory consists of 2,500 frames that describe the atomic 3D positions over time.

Model architecture. We implement the deep spectral clustering module in our approach using two graph pooling layers, which first coarsen the dynamics into 100 nodes and subsequently into 50 nodes. The dimension of hidden states for the stacked GCN layers is kept as 20. In the path signature transform module, the input coarse grained dynamics are partitioned into four segments to facilitate the subsequent application of the log-signature transform.

Training details. The model is trained with mini-batches of size 16 for a total of 200 epochs. We employ the Adam optimizer with an initial learning rate of 5e-5 and a weight decay 1e-5 to optimize the model parameters. Furthermore, the scaling parameters are specified as $\lambda_{1}=1$, $\lambda_{2}=0.01$, and $\lambda_{3}=10$.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. Existing equivariant neural networks with geometric inductive biases have been outperformed, both in terms of performance and efficiency, by Transformer-style architectures (ie, something as simple as torch.nn.TransformerEncoder). Was any ablation done that compares Deep Signature with such architectures that "tokenize" the trajectory intervals?   

2. I'm curious what the learned coarse-grained beads/maps look like for GPCR proteins – is it creating CG beads only based on atoms within a locality? How does this compare to naively taking atoms within some radial ball and considering them a bead? More often than not, these naive CG choices work well in practice, without the need to overcomplicate it with a learnable method. I'd like to see whether this "learned GCN-based coarse-graining" is really necessary.

3. How sensitive is the path signature transform to very dynamic conformational changes in a short duration of time? For instance, as I mentioned above, fast-folding proteins undergo significant changes in 3D shape in just short simulations. Can this method capture this dynamics well enough?
soundness: 3
strengths: 1. Incorporates appropriate symmetries for (temporal) 3D point cloud learning. In particular, the invariance to time reparameterization is key because the underlying MD simulations might be prone to random restarts and miscellaneous artefacts preventing them from being a smooth "video" (ie, the MD itself might be erratic with the same state sampled repeatedly).

2. Deep Signature _learns_ the ideal CG beads relevant to the task, bypassing the need to manually remove degrees of freedom (eg: CA-level coarse-graining, as done in existing methods). The use of the signature transform allows for local and global temporal interactions to be captured well, as opposed to learning representations on separate frames. 

3. I'm confident this method can be used for good representation learning of trajectory information for tasks beyond the ones mentioned in the paper (eg: interpolation in the latent space of protein conformations, perhaps for ensemble generation).

4. Demonstrates relatively good performance compared to weak and strong temporal graph learning baselines (like GraphLSTM).
summary: The paper introduces Deep Signature, a framework for characterizing dynamics on graphs. Deep Signature involves three key components: (1) graph coarsening using a deep spectral clustering module, (2) computation of path signatures to capture global inter-node interactions using iterated integrals over time, and (3) a two-layer MLP for property prediction.
The authors tested Deep Signature on three datasets: (1) gene expression dynamics, classifying into degradation or dimerization types, (2) GPCR dynamics, distinguishing between active and inactive states, and (3) EGFR mutant dynamics, predicting drug sensitivity. They conducted an ablation study validating the contribution of each loss component in the model and performed limited comparisons with methods based on static structures, dihedral angles, and GraphLSTM.
weaknesses: The paper does not clarify the level of coarsening necessary to make path signature computations feasible. For instance, in gene regulatory dynamics, the graph was reduced from 100 nodes to 30, while the EGFR dataset was reduced from approximately 5000 atoms to 50 nodes, but the level of coarsening for the GPCR dataset is unspecified. This raises questions about whether 50 nodes is a computational limit for the method. Additionally, while protein structure can be intuitively coarsened into backbones, sidechains, and motifs, the interpretation for coarsened graphs in gene expression dynamics is unclear. It is uncertain if the coarsened nodes correspond to gene hubs or other biologically meaningful groupings.

The equation used to model GRN dynamics (Eqn 13 in Section B.1) appears incorrect. Specifically, for dimerization (when f=2), the concentration should should be squared under Michaelis-Menten kinetics, rather than simply doubling the decay rate. This could lead to inaccurate modeling of gene expression dynamics and affect the results in this section.

The comparisons used in the experiments are limited and static, primarily involving the first and last frames (head, tail, and head & tail) which do not capture temporal dynamics. The authors do not benchmark against dynamic approaches that consider time-varying information, such as MDTraj [1], Timewarp [2], and DSR [3]. Including these comparisons would provide a more rigorous evaluation of Deep Signature's effectiveness relative to established tools for molecular trajectory analysis.

The authors claim that the coarsened dynamics in Fig 6c follow the same trend as the original dynamics, yet this similarity is not quantified. Providing a quantitative measure, such as correlation coefficients between the original and coarsened dynamics at various coarsening levels, would better support this claim. Additionally, the paper would benefit from comparisons of the authors' coarsening strategy against other deterministic and learnable methods for protein graph coarsening to demonstrate its effectiveness and fidelity in preserving dynamics.

The description of the cross-validation process and test set creation is confusing. The authors mention “for each running, we evaluate the prediction accuracy of our method on an independent unseen test set and report the averaged results,” but it is not clear how this set is constructed. If results are averaged, please report the standard deviation in all the tables.
[1] MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. Biophysical Journal, 2015.
[2] Timewarp: Transferable Acceleration of Molecular Dynamics by Learning Time-Coarsened Dynamics. NeurIPS 2023
[3] DSR: Dynamical Surface Representation as Implicit Neural Networks for Protein. NeurIPS 2023.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 3
questions: 1. I think the classification task using evolution trajectories is not a common setting. Could you explain why this problem setting is reasonable comparing to the one frame classification setting?  Do we really need the complete trajectory to do classification?

2. I do not understand your setting in the EGFR dynamics experiment since I'm not an expert on this domain. Could you please explain why the trajectory can be labeled according to its sensitivity towards the drug? I think the sensitivity should be a number rather than a 0/1 label.
soundness: 2
strengths: 1. The presentation of the method is very clear and easy to understand.
2. The motivation of applying signature transform is reasonable.
summary: In the paper, the authors introduce the Deep Signature framework to capture complex dynamics using the evolution trajectories. The Deep Signature framework includes spectral clustering, signature transform and a classifier. Additionally, the authors show that the framework satisfies the desired symmetry constraints. Experiments on  gene regulatory dynamics, EGFR mutation dynamics and GPCR dynamics exhibit the empirical performance of the framework.
weaknesses: 1. Lack of experiments on large dataset. As the paper claims, the Deep Signature framework can capture large-scale complex dynamics. So I think experiments on datasets with large amount of data and system size are necessary. But the paper only includes the experiments on datasets with large system size.

2. I think the baseline in this paper is too weak. For example, the author should compare the strong baseline with graph transformer architecture[1] based on the first/last frame of the trajectory. Comparations between these strong baseline may strengthen the empirical performance of the framework.

[1].Ying, Chengxuan, et al. "Do transformers really perform badly for graph representation?." Advances in neural information processing systems 34 (2021): 28877-28888.

---

confidence: 2
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. How were the hyperparameters chosen, such as the loss coefficient parameters $\lambda_i$,  the number of nodes in deep spectral clustering model, and time interval $[r_i, r_{i+1}]$ in path signature transform?
2. When visualizing critical pathways and interatomic interactions on the EGFR dynamics, why were three atoms identified specifically? Could more than three be selected?
3. Have the authors considered comparing their approach with advanced time series classification algorithms？
4. How does the computational efficiency of the proposed method compare to baseline methods?
soundness: 2
strengths: 1. Authors develop an end-to-end framework to characterize interatomic interactions and dynamics of large-scale molecules, it shows improvement on three benchmarks and provides interpretability. 
2. The size of the system is reduced by deep spectral clustering module without any expert knowledge. 
3. The framework's desirable properties are supported by theoretical analysis.
summary: The paper presents a framework, Deep Signature, designed to analyze the dynamics and interatomic interactions of large-scale molecular system. The method uses a deep spectral clustering model to capture coarse grained dynamics, a path signature transform module to characterize interatomic interactions through iterated integrals, and a classifier for property prediction. Theoretically, Deep Signature is shown to maintain desirable symmetry properties. The method demonstrates improved accuracy across three benchmarks and demonstrates interpretability.
weaknesses: 1. The authors compare their approach to baseline methods, but a more comprehensive comparison with some SOTA baselines would provide a more robust evaluation. 
2. The manuscript would benefit from a comparison of deep spectral clustering module with existing coarse graining methods. 
3. An analysis of the model’s sensitivity to hyperparameters would provide insights into its robustness and reproducibility.

Meta-review: This paper introduces Deep Signature, a novel and computationally efficient framework for characterizing complex dynamics and interatomic interactions based on their evolving trajectories. The authors propose an approach that combines soft spectral clustering for local aggregation of cooperative dynamics and a signature transform that collects iterated integrals to provide a global representation of non-smooth interactive dynamics. Theoretical analysis establishes that Deep Signature possesses several desirable properties. The experimental evaluation on three benchmark biological processes demonstrates that the proposed framework outperforms baseline methods, highlighting its potential for applications in dynamic systems analysis.

Strengths:

1. The method is presented in a clear and easily understandable manner.

2. The proposed approach is well-motivated and thoughtfully designed, incorporating appropriate symmetries, learning relevant coarse-grained (CG) beads for the task, and requiring no expert knowledge.

3. The framework’s desirable properties are backed by solid theoretical analysis.

4. The method demonstrates strong potential for effective representation learning of trajectory information, with competitive performance relative to both weak and strong temporal graph learning baselines.

Weaknesses:

1. There is a lack of experiments conducted on large datasets.

2. The experimental evaluation and comparison to baseline methods remain insufficient.

3. Several typographical errors and unclear presentations need attention.

Overall, the proposed method is novel and offers significant value to the molecular dynamics (MD) field, with desirable performance that suggests its potential for impactful downstream applications. After reviewing the rebuttal, all reviewers provided positive assessments. Therefore, I recommend acceptance.

Decision: Accept (Poster)