Summary of responses to all reviewers

We thank all the reviewers for their constructive reviews and relevant feedback. The comments regarding the theoretical contribution were generally positive, but several of the reviewers raised concerns regarding the scope of our empirical evaluations. Although we believe that the primary contribution of our paper is the theoretical connection between inverse folding likelihoods and thermodynamic stability, we agree that expanding the empirical evaluations would further strengthen the paper. We have therefore initiated the following experiments:

1. Repeat experiments with ProteinMPNN, a leading inverse-folding model,
2. Replace simplistic site-independent model for $p(\mathbf{a})$ with a pLM,
3. Rerun our protocol using ensembles obtained from BioEmu, thereby avoiding the computational bottleneck of molecular dynamics – which would allow us to scale to larger datasets.

We hope to have preliminary results for 1) ready within the rebuttal period, and 2) and 3) will be added before camera-ready. Finally, we will change the paper title to "Zero-shot protein stability prediction by inverse folding models: a free energy interpretation" – to more specifically reflect its contents.

Response to Reviewer NjCR

We thank Reviewer NjCR for their detailed feedback to our paper. We give a point-by-point response to the identified weaknesses and questions below:

Assume the same beta (and thus temperature) across the data source (PDB)

We agree that the assumption about constant $\beta$ (inverse temperature) across all of PDB (eq. 13) is quite strong. However, we have since realized that this assumption can be substantially relaxed. It is sufficient to assume (cf. eq. 11) that each structure in the database is sampled from its own Boltzmann distribution with a potentially different (and unknown) $\beta$, that is

$$p_{\operatorname{D}}(\boldsymbol{\chi} | \boldsymbol{a}, \beta) \approx p(\boldsymbol{\chi} | \boldsymbol{a}, \beta).$$

Crucially, the downstream analysis does not require access to the individual $\beta$-values. In fact, we can also include the inverse temperature $\beta$ in our model from eq. 12, meaning that we assume a joint model $p_\theta(\mathbf{a}, \mathbf{x}, \beta)$ over sequence, structure and inverse temperature. This weaker assumption above (i.e., allowing for variation in $\beta$ across the database) then implies (cf. eq. 13) that

$$p_{\theta}(\mathbf{x} | \mathbf{a}, \beta) \approx p_{\operatorname{D}}(\mathbf{x} | \mathbf{a}, \beta) \approx p(\mathbf{x} | \mathbf{a}, \beta) ,$$

Using this assumption, we can now explicitly condition all factors in eq. 17 on the specific $\beta$ beta that the structure is observed and simulated at, and we get

$$\frac{p(\mathbf{x} | \mathbf{a}', \beta)}{p(\mathbf{x} | \mathbf{a}, \beta)} \approx \frac{p_\theta(\mathbf{x} | \mathbf{a}', \beta)}{p_\theta(\mathbf{x} | \mathbf{a}, \beta)} = \frac{\frac{p_\theta(\mathbf{a}' | \mathbf{x}, \beta) p_\theta(\mathbf{x}| \beta) } {p_\theta(\mathbf{a}'| \beta)}}{\frac{p_\theta(\mathbf{a} | \mathbf{x}, \beta) p_\theta(\mathbf{x}| \beta) } {p_\theta(\mathbf{a}|\beta)}} = \frac{p_\theta(\mathbf{a}' | \mathbf{x}, \beta)}{p_\theta(\mathbf{a} | \mathbf{x}, \beta )} \frac{p_\theta(\mathbf{a}|\beta)}{p_\theta(\mathbf{a}'|\beta)}$$

We now make the same weak assumption that is implicitly made when training an inverse folding model: that, conditioned on a given structure, the distribution over amino acid sequences is independent of the temperature at which the structure was recorded. In other words, we treat the structure as a complete description of the structural state, which means that $p_\theta(\mathbf{a} | \mathbf{x}, \beta) = p_\theta(\mathbf{a} | \mathbf{x})$. We make the similar weak assumption that the sequence distribution is independent of the temperature, i.e., $p_\theta(\mathbf{a} | \beta) = p_\theta(\mathbf{a})$. This gives us that

$$\frac{p(\mathbf{x} | \mathbf{a}', \beta)}{p(\mathbf{x} | \mathbf{a}, \beta)} \approx \frac{p_\theta(\mathbf{a}' | \mathbf{x})}{p_\theta(\mathbf{a} | \mathbf{x})} \frac{p_\theta(\mathbf{a})}{p_\theta(\mathbf{a}')}.$$

This can now be substituted into eq. 16, yielding

$$\beta \Delta \tilde{G}\_{\mathbf{a}' \to \mathbf{a}}^{S} = \ln \frac{p(S | \mathbf{a}', \beta)}{p(S | \mathbf{a}, \beta)} \approx \ln \mathbb{E}\_{\mathbf{x} \sim p\_\theta(\mathbf{x} | S, \mathbf{a}, \beta)} \left[ \frac{p\_\theta(\mathbf{a}' | \mathbf{x})}{p\_\theta(\mathbf{a} | \mathbf{x} )} \right] \frac{p\_\theta(\mathbf{a})}{p\_\theta(\mathbf{a}')} .$$

This provides a temperature-dependent estimate of $\beta \Delta \tilde{G}\_{\mathbf{a}' \to \mathbf{a}}^{S}$, where the temperature is determined by that of the MD simulation, i.e., it is specified in the distribution $p_\theta(\mathbf{x} | S, \mathbf{a}, \beta)$. We have updated the paper to reflect this relaxed assumption.

The inverse folding model estimates $p(a|x)$, however is only trained on folded CATH structures

It is indeed not obvious that this is reasonable. Empirically, we found that inverse folding models do seem to provide somewhat reasonable amino acid preferences for unfolded states. A possible explanation is that in the autoregressive factorization of an inverse-folding model, the distribution of the $i$'th amino acid is typically conditioned only on the local structural environment surrounding it. Some of these local structural environments (e.g., coil or surface regions) will be similar to those of unfolded proteins, providing some degree of generalization to the unfolded state. We rely on similar generalizations even within the folded state: if $x$ actually represented the full structural state, CATH would provide us only with a single sequence for $\mathbf{a}$ given $\mathbf{x}$, and our model would predict a delta function for each amino acid, providing no information about alternative amino acids (or $\Delta\Delta G$).

The related work is short and narrowly scoped. It does not discuss other avenues of estimating protein stability beyond inverse folding methods.

Due to space constraints, we indeed ended up with a quite narrowly scoped related work section. With the extra page provided for camera-ready, we will broaden it to include a broader overview of the developments over the last decades, including methods such as FoldX and Rosetta.

How does equation 10 connect to the rest of the technical section?

Equation 10 is used in 3.3.2. We acknowledge that this was not very clear, since 3.3.2 refers to appendix A.2, which in turn refers to eq 10. We will add a direct reference to eq 10 from 3.3.2 to make the link clearer.

How does section 3.4 Change in Stability with sequence models tie into estimating protein stability using inverse folding models?

Section 3.4 was introduced to motivate the way we deal with the unfolded state (using amino acid frequencies from the unfolded state). However, the section indeed suggests a protocol relying entirely on sequence-only models, which we did not find time to test empirically. Since we will be evaluating pLMs for $p(a)$ anyway (see response to reviewer sDYD), we will add the results of this protocol as well. This should integrate section 3.4 more closely with the remaining story.

Why does the inverse folding model compute a meaningful $p_\theta(a|x)$ for unfolded structures, as it was never trained on unfolded structures?

See discussion above.

Could the authors elaborate on the position independent sequence model estimated from described in L248?

This is a simple amino acid frequency model estimated from the human proteome (UP000005640_9606.fasta). As was also pointed out by reviewer sDYD, it would be meaningful to estimate $p(a)$ using a pLM, and we will include this result in the new version of the manuscript.

Summary of responses to all reviewers

We thank all the reviewers for their constructive reviews and relevant feedback. The comments regarding the theoretical contribution were generally positive, but several of the reviewers raised concerns regarding the scope of our empirical evaluations. Although we believe that the primary contribution of our paper is the theoretical connection between inverse folding likelihoods and thermodynamic stability, we agree that expanding the empirical evaluations would further strengthen the paper. We have therefore initiated the following experiments:

1. Repeat experiments with ProteinMPNN, a leading inverse-folding model,
2. Replace simplistic site-independent model for $p(\mathbf{a})$ with a pLM,
3. Rerun our protocol using ensembles obtained from BioEmu, thereby avoiding the computational bottleneck of molecular dynamics – which would allow us to scale to larger datasets.

We hope to have preliminary results for 1) ready within the rebuttal period, and 2) and 3) will be added before camera-ready. Finally, we will change the paper title to "Zero-shot protein stability prediction by inverse folding models: a free energy interpretation" – to more specifically reflect its contents.

Response to Reviewer r63Z

We thank the Reviewer r63Z for their feedback on our paper. It is a valid point that we cannot rule out that a similar outcome could have been obtained through a different line of reasoning. However, the assumptions we make do not arise out of the blue. They are based on approximations of the underlying biophysical systems. In our view, the value of our work lies in the fact that our derivations make it possible to have very specific discussions about the types of assumptions we make implicitly when using inverse folding models – and potentially develop better models in the future by revisiting these assumptions.

As also discussed in our response to PxEM, we considered our empirical studies primarily as a sanity check on the underlying theory, and were comforted to see that the correction terms generally contribute positively to performance. It is true that many of the correction terms contribute only marginally to the $\Delta\Delta G$ prediction performance. However, apart from the noisy Guerois dataset, we consider the total performance gains as quite substantial (0.66 vs 0.72 for ProteinG and 0.51 vs 0.54 for VAMP-seq), in particular when compared to the small variations typically reported between the top performing methods on ProteinGym.

We acknowledge that stronger conclusions about the relative merits of the individual correction terms would require a more elaborate set of experiments. Scaling up to more systems has been hindered by the fact that we need to conduct full molecular dynamics simulations for each considered system. In an attempt to solve this issue, we are currently investigating the use of BioEmu as a replacement for full molecular dynamics simulations (for details, see the general post above).

Summary of responses to all reviewers

We thank all the reviewers for their constructive reviews and relevant feedback. The comments regarding the theoretical contribution were generally positive, but several of the reviewers raised concerns regarding the scope of our empirical evaluations. Although we believe that the primary contribution of our paper is the theoretical connection between inverse folding likelihoods and thermodynamic stability, we agree that expanding the empirical evaluations would further strengthen the paper. We have therefore initiated the following experiments:

1. Repeat experiments with ProteinMPNN, a leading inverse-folding model,
2. Replace simplistic site-independent model for $p(\mathbf{a})$ with a pLM,
3. Rerun our protocol using ensembles obtained from BioEmu, thereby avoiding the computational bottleneck of molecular dynamics – which would allow us to scale to larger datasets.

We hope to have preliminary results for 1) ready within the rebuttal period, and 2) and 3) will be added before camera-ready. Finally, we will change the paper title to "Zero-shot protein stability prediction by inverse folding models: a free energy interpretation" – to more specifically reflect its contents.

Response to Reviewer sDYD

We thank the Reviewer sDYD for their positive remarks and detailed feedback. We agree that it would be natural to use a language model for calculating $p(\mathbf{a})$, and will add this result to the paper. As stated above in response to reviewers u4zq and PxEM, we also agree that the paper would benefit from a more elaborate empirical evaluation, and we are currently setting up experiments using ProteinMPNN and BioEmu (see details above). Thanks for pointing out the difficulties reading sections 4.2 and 4.3. We will make an effort to rephrase these. And thanks for pointing out the list of typos, which we will incorporate immediately.

Please see the point-by-point answers to the questions below:

1. Equation 22: Should the second term be calculated with the same model as the first term or can different models be used? In principle, it should be calculated using the same model. So ideally we would want a joint model that supported marginalizing out the structural context, and we could use the same model for both terms. In practice, inverse-folding models do not support this, and we therefore typically employ different models for these two terms: an inverse-folding model for the first and a sequence-only model for the 2nd.

2. Can the authors clarify how the first term in equation 25 can be calculated or how this conditioning can be trained? The first term in equation 25 is the probability of observing an amino acid sequence, conditioned on being in an unfolded state. It is indeed not trivial how to evaluate this in practice – and we therefore dedicate section 4.3 to discussing various alternatives. The most simple solution that we found is to use the amino acid frequencies obtained from intrinsically disordered proteins for this quantity. Given your comment above about section 4.3, we acknowledge that this was not sufficiently clear in the paper, and will elaborate on this explanation in the new version of the paper.

3. (line 220) The authors mention that they use predicted structures for UR50. Does this mean that some of the structures used for the datasets analyzed are predicted structures? Apologies for the confusing sentence structure. The intention was to communicate that the original authors of ESM-IF trained their model on both on CATH AND on predicted structures for UR50. When evaluating the models, we only used crystal structures and structures from molecular dynamics simulations. We have rephrased the sentence to avoid the confusion.

4. What is the meaning of the sentence in lines 232-235: "While data generated … for protein stability." We apologize for this poorly formulated sentence. We meant to say that the abundance signal measured by VAMP-seq is a convoluted signal of thermodynamic stability, effects on local unfolding stability and interactions with the cell's quality control system, and we therefore do not expect to get as clear a signal on protein stability from this type of experiment. We will rephrase to make this clearer.

5. (lines 246-249) Why a position-independent sequence model? Why did the authors not test a masked protein LM in this case? The choice of the position-independent sequence model was made as a simple proxy which had similar sequence statistics as the training data of the ESM-IF model (which was trained on UniRef50 sequences). Protein language models like ESM2 are trained on similar sequence datasets, and should therefore indeed be directly applicable as well. We will do so in the new version of the paper.

6. The writing and explanation of section 4.3 are confusing. It needs additional explanation, intuition, and, if possible, a visualization. As there is no clear way to do that, readers should be clear about the reasoning of the authors for generating or modeling these unfolded ensembles. We acknowledge that this was not sufficiently clear. One of the main problems seems to be that we do not explicitly relate the discussion in 4.3 back to the terms we wish to approximate (e.g. the first term in eq 25). We will rework these sections to make this clearer, possibly adding a figure if necessary.

7. Can the authors clarify the differences between this work and the work done by Jiao et al, 2024? (currently being discussed in lines 301-304 and lines 319-321) Our work considers the more general case of the full structural ensembles, while Jiao et al. represent structural states by single conformations. In addition, Jiao et al. consider protein binding while we consider protein stability, but as discussed in our response to reviewer u4zq above, these tasks are closely related, and it would be possible to use our approach in the setting introduced by Jiao et al., by considering the full structural ensembles.

8. How can the foundations presented in the work be applied to generate better objectives for self-supervised losses when training models using datasets like the PDB? Do the authors think these can also help model more efficient losses for using datasets with augmented predicted structures? Indeed, it would be possible to construct an self-supervised objective, in a very similar way as Jiao et al. (2024), but in our case the loss would be evaluated over a full structural ensemble. It would be interesting to explore whether this leads to improved binding energy prediction. However, due to the considerable experimentation required, we consider it outside the scope of the current manuscript.

Summary of responses to all reviewers

We thank all the reviewers for their constructive reviews and relevant feedback. The comments regarding the theoretical contribution were generally positive, but several of the reviewers raised concerns regarding the scope of our empirical evaluations. Although we believe that the primary contribution of our paper is the theoretical connection between inverse folding likelihoods and thermodynamic stability, we agree that expanding the empirical evaluations would further strengthen the paper. We have therefore initiated the following experiments:

1. Repeat experiments with ProteinMPNN, a leading inverse-folding model,
2. Replace simplistic site-independent model for $p(\mathbf{a})$ with a pLM,
3. Rerun our protocol using ensembles obtained from BioEmu, thereby avoiding the computational bottleneck of molecular dynamics – which would allow us to scale to larger datasets.

We hope to have preliminary results for 1) ready within the rebuttal period, and 2) and 3) will be added before camera-ready. Finally, we will change the paper title to "Zero-shot protein stability prediction by inverse folding models: a free energy interpretation" – to more specifically reflect its contents.

Response to Reviewer PxEM

We thank the Reviewer PxEM for their positive evaluation of our theoretical contribution and their critical assessment of our paper. Our main motivation with this paper was to attempt to understand why log-odds score from inverse-folding models correlate so well with stability changes. From this perspective, we did not consider it a problem that not all correction terms provide substantial improvements, and considered the fact that these correction terms generally had a positive effect mostly as a sanity-check of the underlying derivations.

Of course, we agree with the reviewer that the paper would benefit from an expanded set of experiments – both to get some more statistically significant results and because the assumptions we make might hold to a different degree for different protein systems. As we have discussed above with reviewer u4zq, the main bottleneck was conducting properly converged molecular dynamics simulations. We agree with the reviewer that this computational complexity makes our approach less interesting for practical application. In an effort to overcome this, we are now repeating our experiments using ensembles from the BioEmu model, with the hope that the structural ensembles sampled from this model are of sufficient quality to make it an alternative to molecular dynamics. If this proves successful, it allows us to scale our analysis to a larger stability subset of ProteinGym (e.g., the Tsuboyama subset), and also suggests a path towards making our approach more useful for practical use.

We agree that it makes sense to compare to MPNN as well. We chose ESM-IF because it generally performs slightly better on zero-shot stability prediction than ProteinMPNN (according to ProteinGym), but the reviewer is right that the ubiquitous use of ProteinMPNN makes it a relevant comparison. We are currently running these experiments and hope to have them done before the rebuttal period ends.

We thank the reviewer for their comments regarding Figure 1 and the abstract/introduction, which we will incorporate in the new version of the manuscript.

It is a good question whether subsets of ProteinGym – or differences in experimental settings – would display different tendencies. We could clearly expect the most meaningful results for subsets categorized under the "Stability" label. In particular, the Tsuboyama megascale experiments are direct probes for protein stability and should thus behave well. The VAMP-seq datasets – like the one we use in our paper – measure abundance, but have previously been shown to correlate well with stability, and should thus be another meaningful subset of ProteinGym to consider (although likely displaying somewhat higher levels of noise).

We do see differences in the ranking of mutants between different ablations of terms, as reflected by the Spearman rank correlation scores in Figure 1 and Table 1. It would indeed be interesting to investigate whether there are any patterns to these effects – e.g., whether there is a clear difference between surface and core residues.

Summary of responses to all reviewers

We thank all the reviewers for their constructive reviews and relevant feedback. The comments regarding the theoretical contribution were generally positive, but several of the reviewers raised concerns regarding the scope of our empirical evaluations. Although we believe that the primary contribution of our paper is the theoretical connection between inverse folding likelihoods and thermodynamic stability, we agree that expanding the empirical evaluations would further strengthen the paper. We have therefore initiated the following experiments:

1. Repeat experiments with ProteinMPNN, a leading inverse-folding model,
2. Replace simplistic site-independent model for $p(\mathbf{a})$ with a pLM,
3. Rerun our protocol using ensembles obtained from BioEmu, thereby avoiding the computational bottleneck of molecular dynamics – which would allow us to scale to larger datasets.

We hope to have preliminary results for 1) ready within the rebuttal period, and 2) and 3) will be added before camera-ready. Finally, we will change the paper title to "Zero-shot protein stability prediction by inverse folding models: a free energy interpretation" – to more specifically reflect its contents.

Response to Reviewer u4zq

We thank Reviewer u4zq for their positive remarks. We agree that the test sets used in our paper were not as extensive as one might have wished for. The primary reason for focusing on a limited set of systems was the molecular dynamics simulations: not only the computational cost, but also making sure that the molecular systems were set up correctly, relaxed properly in the forcefield, etc. We believe the current set of experiments is sufficient as a sanity check of the presented theoretical contribution. However, we do agree that it would be valuable to see whether the trends hold across a larger dataset, and to develop a protocol that would make the approach more robustly applicable in practice. We are therefore now in the process of running our experiments using ensembles sampled from BioEmu, and will attempt to scale this to the Tsuboyama et al dataset. We believe we can have this full evaluation completed in time for camera-ready. If the BioEmu ensemble-approach fails, we would be happy to at least extend our current analysis to several more systems using the full molecular dynamics protocol.

We also agree that it would be interesting to see if MPNN shows similar trends. We are setting up these experiments now. This does not involve new molecular dynamics simulations, we expect that we should have preliminary results ready before the discussion period ends, and hope to post the results when they are ready. A priori, it is unclear to us what to expect from MPNN. The ProteinGym zero-shot benchmark for the stability category reports ESM-IF likelihoods to have somewhat higher correlation to experimental values (0.624) than ProteinMPNN (0.565), suggesting that it might perform slightly worse, but the relative trends of the ablation might show very similar behavior.

It would indeed be interesting to investigate our approach in the context of binding energy. Theoretically, it is a simple extension, since the standard way to calculate mutational effects on binding is to look at the difference in mutational effects on the stability of the complex ($A:B$) and the monomers ($A$ & $B$) is $\Delta\Delta G = \Delta G_{A:B}^{\text{wt} \to \text{mut}} - \Delta G_{A}^{\text{wt} \to \text{mut}} - \Delta G_{B}^{\text{wt} \to \text{mut}}$ . In practice, it would require some work to do the analysis, since structural ensembles would need to be acquired for these different structural states, over a range of systems for which we have experimental results (e.g. from SKEMPI as you suggest). Because of the work involved in these practical matters, we defined this to be out of scope of the current paper, and focused our attention on stability. In the new version of the manuscript, we will, however, add a comment to the discussion to make it clearer how our approach applies to the binding case.

The encouraging results from the BioEmu experiments make it possible to scale our experiments across a larger range of proteins. We are now setting up experiments on a 95-protein subset of the Tsuboyama dataset, which was withheld as a validation set by the authors of the BioEmu paper [1], and should thus be a fair way to test the effectiveness of predicted ensembles in our $\Delta\Delta G$-protocol. It is unlikely that all these runs will be finished before the discussion period ends, but we'll post any updates we have prior to the discussion closing.

While still preliminary, we hope the results above provide the reviewers with additional evidence of the practical utility of our contribution and the robustness of our approach to 1) using alternative inverse-folding methods and 2) employing different procedures for generating structural ensembles.

Ultimately, the main contribution of our paper is theoretical, establishing a formal connection between likelihood ratios of inverse folding models and protein relative stabilities. The theoretical part seems to have been well received by all reviewers, but if anyone has any remaining concerns they would like to discuss, please let us know.

[1] Lewis, S., Hempel, T., Jiménez-Luna, J., Gastegger, M., Xie, Y., Foong, A.Y., Satorras, V.G., Abdin, O., Veeling, B.S., Zaporozhets, I. and Chen, Y., 2025. Scalable emulation of protein equilibrium ensembles with generative deep learning. Science, p.eadv9817.

Tsuboyama dataset

Regarding the scaling of our analysis to more systems, we also managed to complete ten out of the ~40 runs on the Tsuboyama dataset (from ProteinGym), using BioEmu to generate the structural ensembles.

The results look similar to the performances we reported for Protein G in the submitted version of the paper, with a consistent and significant increase in performance both when taking into account the structural ensemble (now through BioEmu) and when modelling the unfolded state (using IDP statistics). We will add a comprehensive overview and discussion of these new results in the paper (once the remaining experiments have been completed).

Thanks to the reviewers for their constructive feedback – we believe the additional experiments have significantly improved the empirical part of the paper, clearly showing the theoretical considerations to have significant practical effects as well.