From Biased to Unbiased Dynamics: An Infinitesimal Generator Approach

We appreciate the reviewer's insightful evaluation and valuable comments. In what follows, we aim to address the highlighted weaknesses and respond to the reviewer's questions.

Summary & strengths:

• Thank you for your summary, while it captures the general idea of the paper, maybe it lacks one important aspect that, to the best of our knowledge, we are the first to show that covariance reweighting in transfer operator (TO) models built from biased simulations fails, while infinitesimal generator (IG) models do not suffer from this effect. Indeed we prove that one can learn eigenvalues and eigenfunctions of the IG via its resolvent in a scalable and statistically consistent way from a single biased trajectory.

• Related to the first strength, please refer to our general reply, part Clarity of presentation. Besides proposing Alg. 1 (methodological contribution), we also theoretically prove its statistical consistency. Let us provide the context: data-driven IG methods emerged relatively recently in the literature [1,15,20,21,46], and up to our knowledge, the only method that possesses finite-sample statistical guarantees for the proper spectral estimation is [21]. However, all of the above approaches analyzed learning from unbiased simulations, which is unfeasible in many realistic scenarios. Moreover, the method in [21] suffers from scalability issues inherent to kernel methods. On the other hand, a deep learning method for the IG [46] was applied to biased simulations , but its statistical consistency was an open question.

Weaknesses:

1. Noting that learning the spectral decomposition of IG is a problem of interest besides protein dynamics, e.g. crystallization, we agree with the suggested best practice. We did our best to meet such a requirement, as reported in Experiments of our general reply and the attached pdf.

• We believe that the main reason why IG and TO methods are not initially tested on larger scale systems is the complexity of the task. Namely, recovering proper time scales and the transition state ensemble is a complex task that highly depends on the descriptors of the system and their effective dynamics. Indeed for large proteins one of the main challenges is to design appropriate descriptors so that effective Langevin dynamics can reveal the true dynamics of a protein [D]. These are highly dependent on the protein, and their design is a problem per se. Once this has been resolved, the next challenging task is to explore the rare transitions by biasing dynamics and recover proper time scales and transition state ensembles. While our method is designed to solve the latter, the former is an active field of research that currently limits the development of the full pipeline for large scale complex molecular systems.
• In order to fully address the question on the capability of our method, we have performed an additional experiment. To the best of our knowledge, the most complex benchmark for all TO and IG methods is the (very long) unbiased simulation dynamics of a small protein, typically Chignolin [7,19]. So, we designed a biased simulation experiment based on our method, which is, to the best of our knowledge, the first experiment of its kind for generator learning.
• Concerning baselines for ALDP, since we don’t have ground-truth it would be just a qualitative comparison.

2. Most non DL CV are heuristics based on chemical/physical intuition, they are system dependent and often suboptimal. This has stimulated intense research towards ML based alternatives. In this context, TO/IG methods are well-motivated, since eigenfunctions of these operators are “optimal” CVs for biasing simulations, see e.g. [A]. While some classical ML methods are available for this [21,23], they are intractable at large scale since at each time step one needs to compute the kernel function with respect to the current ensemble of the simulation.
3. Our method is as expensive as standard CV based methods, but it reliably leads to more accurate results, which is important when one wants to reconstruct the free energy profile and use them iteratively to learn the generator. See also our additional experiment on Chignolin, see the global reply and the attached pdf where we highlight the computational speedup relative to unbiased bruteforce simulations.
4. We incorporated your suggestion in the general reply on the general pipeline. Here we address your related questions. First, since the input of the DNN encoder in the representation learning part are the descriptors, the input dimension depends on the problem. On the other hand, the output (latent features) dimension is typically very small, reflecting the number of slow time-scales to learn (typically around 3-5). Second, the choice of using separated DNNs per latent feature is made to aid learning of orthonormal features, since, depending on the architecture, shared weights may implicitly introduce prohibitive constraints. Third, eigenfunctions, defined in line 13 of Alg. 1, map descriptors of the system to a real number that, when properly scaled, has mean zero and variance 1 w.r.t. equilibrium distribution.
5. Eigenfunctions can be plotted in the plane of relevant features to analyze the data. For example, for alanine dipeptide we have for each configuration the values of the angles $\phi$ and $\psi$ and the values of the eigenfunction, therefore, we may plot the points in the $\phi-\psi$ plane colored according to the eigenfunction value.

Questions:

1. If we understand well, two statements are equivalent. Indeed the noise of the feature is $z(X_{t+s}) - \mathbb{E}[z(X\_{t+s})\vert X_t] = z(X\_{t+s}) -[ {\cal T}\_s z ] (X_t)$, and from one trajectory we only observe an instance of the r. v. $z(X_{t+s})$, and not the output of TO.
2. Please see our general reply.

Ref: [D] Zhang et al., Effective dynamics along given reaction coordinates and reaction rate theory. Faraday Discuss. (2016)

Thank you very much for confirming some points, and for addressing the reference. I would like to take the freedom to raise my score from 6 to 8 and to increase my confidence.

We appreciate the reviewer's insightful evaluation and valuable comments. In what follows, we aim to address the highlighted weaknesses and respond to the reviewer's questions.

Weaknesses:

We hope that our general reply brings more clarity. We commit to incorporate all the suggestions and additional discussions in the revised version of the paper. If the reviewer would like us to clarify any other specific aspects, we are happy to address them during the discussion period.

Questions:

Thank you for raising this question, in the following we try to briefly review the usefulness of generator eigenfunctions (“modes”) for interpretability and control of Langevin dynamics. We plan to add it in the revision to further clarify the rationale of our method for the general reader. From the spectral decomposition of the generator in Equation (4) and its link to the transfer operator in Equation (5) , it can be seen that each mode $f_i(x)\langle f_i,f \rangle$ (essentially the eigenfunction $f_i$) is related to a timescale of the process ($-\lambda_i>0$). This means that eigenfunctions corresponding to the slow time-scales (eigenvalues close to zero ) give information on metastable states. Recalling that normalized eigenfunctions (in $\mathcal{L}_\pi^2$) have mean zero and unit variance, they tend to have constant values for ensembles $x$ in such meta-stable states (meaning at time-scale $-\lambda_i>0$ there is no dynamics), while it changes sign when ensembles move from one meta-stable state to another. This is what the reviewer called “order parameter”. But it is more than this, because contrary to an order parameter which only discriminates two states, the eigenfunction will also give valuable information on how the transition takes place, see illustrative example in appendix E.1 and corresponding Figure 5. This feature of the eigenfunctions is of particular relevance, since the knowledge of rare transitions can be used by experimentalists to accelerate or slow down the process, see e.g. Wei Zhang, Christof Schütte, Understanding recent deep‐learning techniques for identifying collective variables of molecular dynamics. PAMM 2023, 23 (4).

We appreciate the reviewer's insightful evaluation and valuable comments. In what follows, we aim to address the highlighted weaknesses and respond to the reviewer's questions.

Weaknesses:

Thank you for motivating us to present the overarching summary of our approach. While we discuss this in the general reply, we would like to provide a few more relevant details here.

• First, data-driven generator methods emerged relatively recently in the literature [1,15,20,21,46], and up to our knowledge, the only method that possesses finite-sample statistical guarantees for the proper spectral estimation is [21]. However, all of the above approaches analyzed learning from unbiased simulations, which is unfeasible in many realistic scenarios. Moreover, the method in [21] suffers from scalability issues inherent to kernel methods. On the other hand, deep learning methods for the generator were applied to biased simulations [46], but their statistical consistency was not proven. In sharp contrast to all of the above, our work is the first to show that one can learn eigenvalues and eigenfunctions of the infinitesimal generator via its resolvent in a scalable and statistically consistent way from biased simulations.
• To summarize, we extend the method in [21] by incorporating biased dynamics, developing appropriate representation learning, and designing the overall scalable and statistically consistent approach. Our method is twofold: first, we learn a basis set using Theorem 2 (representation), and then we learn the generator’s resolvent on this basis set using Theorem 1 (regression). This approach is common in many fields, notably in transfer operator approaches [24,29], but is applied for the first time ever, to the best of our knowledge, to the case of the infinitesimal generator. In doing so, we followed the core idea of [21] to formulate the problem of learning the resolvent so that the geometry of the process is efficiently exploited via the energy norm, see also reply to Q1. Practically, this means that one can avoid inverting the shifted generator by rescaling the norm in which one learns the resolvent. This allows one to work in the latent space where the state becomes represented as $(\eta I - L)^{1/2} z$ and all the relevant inner products and norms can be computed using the features $z$ and their gradients.

Questions:

• Q1: In fact, the estimator is not a Galerkin projection of the resolvent onto $\mathcal{H}$ as a subspace of $\mathcal{L}^2_\pi$. This would be the case if in (11) we replace the energy with the expectation w.r.t. invariant distribution. However, note that this is not feasible since we cannot compute/estimate the action of the resolvent. To solve this issue we change the domain of the problem from $\cal{L}^2_\pi$ (expectation) to $\cal{H}\^\eta\_\pi$ (energy). So, the reviewer can think of our estimator as a (regularized) Galerkin projection of the resolvent acting on $\cal{H}\^\eta\_\pi$ onto a subspace $\cal{H}$. The reason why we do so, as shown in (12), lies in the fact that Hilbert Schmidt norm of an operator $A$ on \cal{H}^\\eta_\\pi becomes $\rm{tr}(A^*(\eta I-L)A)$. That is, the features $z\in\mathcal{L}^2_\pi$ are transformed into features $\tilde z:=(\eta I - L)^{1/2}) z$, so inner products $\langle \widetilde z, \widetilde z’ \rangle = \langle z,(\eta I - L) z’ \rangle$ and $\langle \widetilde z, (\eta I - L)^{-1} \widetilde z’ \rangle = \langle z, z’ \rangle$ can be now estimated from data. Using this trick we neatly avoid estimating the action of the resolvent, and no differentiation in (11) is needed since integral in $\chi_\eta$ is essentially canceled out by differential operator $\eta I -L$. This is what we meant by contrasting the resolvent in line 182. To further clarify our approach to the reader, we will include this discussion in the revision.

• Q2: Do we correctly understand that the reviewer is referring to Figure 1 instead of 2, since there we compare to [46] (Zhang et al. 2022)? If so, let us note that, as pointed out, we can use both methods for representation learning (loss (19) or loss 39 from [46]) to identify a subspace and then apply our regression method. Based on the theoretical background given in [46], we believe that their method indeed reduces representation error in Theorem 1, making this approach sound. However, when we performed this additional test, the first eigenvalue went from 0.043 to 0.013, the second one went from 1.13 to 4.19. This indicates that even with regression the result does not significantly improve. To strengthen our point, we also did the opposite operation: we computed the Rayleigh quotient from our features without doing any regression, just as is done in the work [46] we obtain for the first eigenvalue 1.79e-4 and 0.36 for the second one, which are close to the ones found by regression and very close to the ground truth values. We are speculating that the main reason is in the efficiency of our loss and its robustness on the hyperparameter choices.

• Q3 and Q4: Thank you for pointing out typos, we will correct them. Also, sorry for the confusion, we intended to say that if $t$ is chosen too small w.r.t. the mixing time, the noise in the corresponding transfer operator regression dominates the signal, meaning that the learning becomes harder. We will correct this discussion in the revision.

Thank you for the detailed response. I have some follow-up points.

• Q1. I now agree that the estimator is a regularized Galerkin projection of the resolvent onto the space of basis functions, using the energy inner product. This interpretation seems more convenient from an expository standpoint. It would be helpful if the authors can incorporate this perspective into the paper and clarify the writing. In my experience, the current mathematical exposition, with its various operators, norms and formulas, did not provide a very pleasant reading experience.

• Q2. Yes, it is Figure 1. Thank you. The explanation and results are interesting. Do the authors have any insights regarding the superiority of the loss (19) for representation learning? It seems to me the loss is not that simple and easy to optimize or compute. There is a penalization term to make the $z_i$ normalized and one has the hyper parameter $\alpha$ to tune. Furthermore, based on the authors' response, do I understand correctly that the primary improvement stems from the loss function (19) for representation learning, which is more significant than the regression formulation proposed in the paper? If so I think these insights should be added to the paper.