Multi-level Interaction Modeling
for Protein Mutational Effect Prediction

Traditional methods for $\Delta\Delta G$ prediction broadly fall into two categories: energy-based and statistics-based methods. Energy-based methods utilize physical energies, such as van der Waals, electrostatic potential and hydrogen bonding for estimation schymkowitz2005foldx ; alford2017rosetta ; barlow2018flex , while statistics-based methods use handcrafted features such as motif properties and evolutionary conservation to predict the mutational effects on binding pires2016mcsm ; li2016mutabind . Both of these approaches highly rely on human expertise, and fail to capture intricate interactions between proteins.

In contrast, deep learning-based methods have shown better promise, which can be categorized as supervised and pre-training-based models. The supervised methods extract features from wild-type and mutant complexes, and directly predict $\Delta\Delta G$ shan2022deep using labeled data in an end-to-end manner. Nevertheless, their development is hindered by the scarcity of labeled protein mutation data. Pre-training-based methods attempt to mitigate this challenge with various pre-training objectives. Some employ masked modeling on protein structures, training models to classify masked amino acids. The predicted probabilities of amino acid types before and after mutation are then used to estimate $\Delta\Delta G$, but the correlation between residue type probabilities and changes in binding free energy tends to be mild hsu2022learning ; bushuiev2023learning . Other methods leverage protein representations pre-trained with sidechain modeling objectives to capture mutation-related features, achieving state-of-the-art performance for $\Delta\Delta G$ prediction liu2021deep ; luo2023rotamer ; liu2023predicting . However, these methods only incorporate sidechain-level interaction modeling, so their effectiveness can be impaired when mutations lead to significant changes in protein backbone conformations and overall binding affinity. To this end, our approach fully models PPI across different levels to facilitate a comprehensive understanding of mutational effects on binding.

Protein-protein interaction modeling has been studied at each level for decades. Since mutations hierarchically impact interactions across different levels, comprehensive modeling of each level is required to fully characterize the mutational effects on PPIs.

Protein-level interaction modeling typically refers to the task of identifying potential protein binding pairs. Experimental methods are often costly and time-consuming fields1994two ; burckstummer2006efficient , so recent computational approaches use deep learning to predict protein interactions. These approaches extract informative features from different modalities of protein data to predict interactions zhang2019sequence ; zhao2023semignn ; bryant2022improved ; gao2023hierarchical ; wu2024mape . Some latest works also explore correlations within PPI networks to infer unknown interactions yang2020graph ; lv2021learning . In the context of mutational analysis, predicting protein interactions is challenging, since a small number of mutations on proteins can lead to significantly different binders hashimoto2010mechanisms ; li2014predicting . Therefore, our method adopts a contrastive manner, which has been successfully verified in multimodal learning li2021align ; radford2021learning , to model protein-level interactions in order to emphasize learning the relationship between inter-protein differences and resulting binding patterns.

The representative task for backbone-level interaction modeling is protein docking, which involves predicting the 3D structures of complexes from unbound states. Classical docking software typically suffers from intensive computational burden. schindler2017protein ; sunny2021fpdock ; torchala2013swarmdock ; vakser2014protein . Deep learning-based methods have emerged as better alternatives, capturing evolutionary constraints and geometric features from growing protein data ganea2021independent ; evans2021protein ; wang2024injecting . Recent approaches focus on rigid-body protein-protein docking, predicting rotation and translation to attain the bound structure ketata2023diffdock ; yu2024rigid ; wu2023neural . Nevertheless, predicting rigid transformations cannot fully capture the complicated spatial changes in inter-molecular distances caused by mutations. Consequently, for backbone-level interaction modeling, our method differs from previous docking approaches by predicting the relative distance between proteins, which implicitly models the flexibility of backbone conformations.

Sidechain-level interaction modeling focuses on predicting rotamers given the protein backbone structures. Traditional methods for rotamer prediction operate by minimizing the energy function across a pre-defined rotamer library xu2006fast ; krivov2009improved ; huang2020faspr ; leman2020macromolecular . Recent methods employ various deep learning techniques such as 3D convolution network misiura2022dlpacker , SE(3)-Transformer mcpartlon2022end and diffusion models zhang2024diffpack ; liu2023predicting , to achieve more accurate rotamer prediction. Among them, several approaches utilize sidechain modeling to predict mutational effects and exhibit superior performance, including RDE luo2023rotamer employing a flow-based generative model to estimate the probability distribution of rotamers, as well as SidechainDiff liu2023predicting using a Riemannian diffusion model to learn the generative process of rotamers. These methods inject mutation-related knowledge into pre-trained protein representations via sidechain-level interaction modeling.

Notations In protein-protein interaction, we denote a complex as $\mathbf{c}\in\mathbb{A}^{n}$, where $n$ represents the number of residues, and $\mathbb{A}$ is the alphabet of amino acids $\{1,\dots,20\}$. Complexes with two binding proteins consist of a receptor protein and a ligand protein, represented as disjoint sets of indices $\mathcal{P}^{r}$ and $\mathcal{P}^{l}$, respectively. Each residue in the complex is characterized using its type $a_{i}\in\mathbb{A}$, position $\boldsymbol{p}_{i}\in\mathbb{R}^{3}$, orientation $\boldsymbol{O}_{i}\in SO(3)$, and the sidechain torsion angles $\boldsymbol{\chi}_{i}=(\chi_{i}^{(k)})_{k=1}^{t}$ where $\chi_{i}^{(k)}\in[0,2\pi)$ and $t$ denotes the number of torsion angles $(i\in\{1\dots n\})$.

Model Overview The overall framework of ProMIM is illustrated in Figure 2. In the early stage of attaining residue-wise representations, we follow RDE luo2023rotamer to capture information about each residue and its structural context. Specifically, given a complex, two protein binders are randomly assigned as the receptor protein and ligand protein. To simulate an unbound state, the receptor protein is kept at its initial location, while random transformations, including rotation and translation, are performed to the ligand protein to change its location. Next, single and pair features are extracted from the complex. Single features include the residue type, backbone dihedral angles, and local atom coordinates, while the pair features include the relative position and residue type pairs. We denote the single embeddings and pair embeddings as $\boldsymbol{e}_{i}$ and $\boldsymbol{z}_{ij}$ respectively, where $i,j\in\{1\dots n\}$. These two features are then transformed into hidden representations $\boldsymbol{h}_{i}$ with the Invariant Point Attention Module (IPA), introduced by AlphaFold2 jumper2021highly . The resulting residue-wise representations $\boldsymbol{h}_{i}$ are then employed to conduct the subsequent interaction-aware pre-training process across three levels.

In this section, we introduce how ProMIM designs pre-training objectives for modeling each level of interaction, to fully capture the protein mutational effects.

Protein-level Interaction Modeling Mutation can change the binding relation between protein pairs, i.e. binding affinity. In this paper, we formulate PIM as a matching problem to capture this binding relation in a coarse grain, namely whether a pair of proteins should bind together or not, which reduces the need for data containing affinity labels. Matching proteins is challenging since proteins with highly similar sequences can participate in different interactions, and even a small number of amino acid mutations can lead to different interaction partners hashimoto2010mechanisms ; li2014predicting . This is similar to the problem of matching documents and images li2021align ; radford2021learning in multimodal learning, where the superficial differences between words and pixels are not the matching key, but the underlying semantics are. Advanced works in this field, such as CLIP radford2021learning , inspire us to adopt the contrastive learning objective to model the protein-level interactions.

Specifically, for two proteins $\mathcal{P}^{l}$ and $\mathcal{P}^{r}$ belonging to a complex $\mathbf{c}$, we apply max-pooling to the residue-wise hidden representations $\boldsymbol{h}_{i}(i\in\{1\dots n\})$ obtained from the encoder to get a global structure representation for each protein, denoted as $\boldsymbol{H}^{l}$ and $\boldsymbol{H}^{r}$. Given a batch of protein complexes $\mathbf{c}_{k}(k\in\{1\dots N\})$ with batch size $N$, we attain global representations for each pair of proteins in the batch using our encoder. Subsequently, we combine them to form $N^{2}$ representation pairs $(\boldsymbol{H}^{l}_{i},\boldsymbol{H}^{r}_{j})$ where $i,j\in\{1\dots N\}$. When $i=j$ it is a positive pair, and when $i\neq j$ it is a negative pair. Two contrastive losses are introduced to facilitate the training process:

$$\mathcal{L}_{k}^{l}(\boldsymbol{H}_{k}^{l},\{\boldsymbol{H}_{i}^{r}\}_{i=1}^{N})=-\frac{1}{N}\log\frac{\exp(s(\boldsymbol{H}_{k}^{l},\boldsymbol{H}_{k}^{r})/\tau)}{\sum_{i}\exp(s(\boldsymbol{H}_{k}^{l},\boldsymbol{H}_{i}^{r})/\tau)},$$

(1)

$$\mathcal{L}_{k}^{r}(\boldsymbol{H}_{k}^{r},\{\boldsymbol{H}_{k}^{l}\}_{i=1}^{N})=-\frac{1}{N}\log\frac{\exp(s(\boldsymbol{H}_{k}^{l},\boldsymbol{H}_{k}^{r})/\tau)}{\sum_{i}\exp(s(\boldsymbol{H}_{i}^{l},\boldsymbol{H}_{k}^{r})/\tau)},$$

(2)

where $\tau$ is the temperature parameter that controls softmax distribution, and $s(\cdot)$ is a function to calculate the similarity score of a given representation pair. Here we adopt cosine similarity. Both two contrastive losses are utilized to identify the binding protein for a given one and enforce symmetry.

The final contrastive loss for a batch is formulated below:

$$\mathcal{L}_{PIM}=\frac{1}{2}\sum_{k=1}^{N}(\mathcal{L}_{k}^{l}+\mathcal{L}_{k}^{r}),$$

(3)

Backbone-level Interaction Modeling At a finer level, mutation can change the backbone conformation of a complex, and predicting this conformation pertains to the protein docking problem. General docking methods typically predict a rotation and a translation on the input unbound structures to obtain their bound state ganea2021independent ; wang2024injecting . However, mutation can lead to complicated variations in inter-molecular distances, which cannot be adequately represented solely by rotation and translation. Therefore, we use the relative distance between protein pairs as the supervision signal for BIM, which can implicitly model backbone conformation such as its flexibility.

To be specific, given the residue-wise representations $\boldsymbol{h}_{i}(i\in\{1\dots n\})$ of a protein complex $\mathbf{c}$ formed by $\mathcal{P}^{l}$ and $\mathcal{P}^{r}$, we firstly transform them into a pairwise representation $\boldsymbol{h}_{lr}^{(0)}$ using a MLP, which serves as the attention bias added to the attention weight. Subsequently, the M-layer Transformer takes both residue-wise and pairwise representations as inputs and generates the final pairwise representation $\boldsymbol{h}_{lr}^{(M)}$. We denote the ground-truth distance map of the original complex structure in the bound state as $\mathcal{D}$, where each element $d_{ij}\in\mathcal{D}$ represents the distance between the $i$-th $C_{\alpha}$ atom in the ligand protein $\mathcal{P}^{l}$ and the $j$-th $C_{\alpha}$ atom in the receptor protein. The predicted distance map $\hat{\mathcal{D}}$ is obtained using a one-layer MLP to transform the final pairwise representation $\boldsymbol{h}_{lr}^{(M)}$. We regard the distance map prediction as a regression task and employ the Mean Square Error (MSE) loss in the training:

$$\mathcal{L}_{BIM}=\mathcal{L}_{MSE}(\mathcal{D},\hat{\mathcal{D}}),$$

(4)

Sidechain-level Interaction Modeling When mutations occur, changes in protein-protein affinity and backbone structure actually stem from alterations in the nearby interaction environment xu2023opus ; eyal2003protein , caused by the changes in sidechain conformations. Previous methods have well studied SIM by modeling rotamer distribution as the pre-training task and obtained great performance on mutational effect prediction. We adopt the flow-based modeling proposed by RDE luo2023rotamer , where a conditional normalizing flow serves as a flexible method to model the complex probability density of rotamers by transforming it into a simple base density in an invertible manner durkan2019neural .

Specifically, a rational quadratic spline flow durkan2019neural ; rezende2020normalizing is firstly utilized to model the distribution of a rotamer with one torsion angle, which is a piece-wise bijective function that contains $K$ pieces delimited by $K+1$ knots, with each piece formulated as follows:

	$\displaystyle f_{k}(x|x_{k,k+1},y_{k,k+1},\delta_{k,k+1})$	$\displaystyle=y_{k}+\frac{(y_{k+1}-y_{k})\left[s_{k}\xi_{k}^{2}(x)+\delta_{k}\left(1-\xi_{k}(x)\right)\xi_{k}(x)\right]}{s_{k}+[\delta_{k+1}+\delta_{k}-2s_{k}]\left(1-\xi_{k}(x)\right)\xi_{k}(x)},$		(5)
	$\displaystyle\mathrm{where~{}}s_{k}$	$\displaystyle=\frac{y_{k+1}-y_{k}}{x_{k+1}-x_{k}},\mathrm{and~{}}\xi(x)=\frac{x-x_{k}}{x_{k+1}-x_{k}}\quad(x\in[x_{k},x_{k+1}]),$

The bijective is parameterized by the coordinates and derivatives of the $K+1$ knots, denoted as $x_{k}$, $y_{k}$, and $\delta_{k}(k\in{1\dots K+1})$, and forms a strictly monotonically increasing function on $[0,2\pi]$, denoted as $f:[0,2\pi]\to[0,2\pi]$. These parameters are conditioned on the hidden representation $\boldsymbol{h}_{i}$ using neural network transformation, so the bijective can also be denoted as $f(x|\boldsymbol{h}_{i})$. The uniform distribution $p_{z}(z)=\frac{1}{2\pi}(z\in[0,2\pi])$ is employed as the base distribution, and $f$ maps from the target rotamer distribution to the base distribution. According to the change-of-variable formula in the probability density function, we attain the target rotamer density as follows:

$$\log p(x|\boldsymbol{h}_{i})=\log p_{z}\left(f(x)\right)+\log|f^{\prime}(x|\boldsymbol{h}_{i})|=-\log 2\pi+\log|f^{\prime}(x|\boldsymbol{h}_{i})|$$

(6)

Subsequently, coupling layers stack multiple bijectives to enable modeling distributions of rotamers with more torsion angles. The training process aims to minimize the negative log-likelihood of ground truth rotamers, where $n$ denotes the number of residues in the complex:

$$\mathcal{L}_{SIM}=-\frac{1}{n}\sum_{i=1}^{n}\log p(\boldsymbol{\chi}_{i}|\boldsymbol{h}_{i}),$$

(7)

ProMIM is pretrained on the PPIRef50K bushuiev2023learning and PDB-REDO joosten2014pdb_redo datasets. Please refer to the Appendix A for training details. Given a wild-type complex $\mathbf{c}_{wt}$ and its mutant complex $\mathbf{c}_{mt}$, we first obtain their residue-wise hidden representations $\boldsymbol{h}_{i}^{wt}$ and $\boldsymbol{h}_{i}^{mt}$ from the pre-trained ProMIM, which provides rich interaction knowledge across three levels. Next, we adopt the same encoder architecture in the pre-training stage for $\Delta\Delta G$ prediction. $\boldsymbol{h}_{i}^{wt}$ and $\boldsymbol{h}_{i}^{mt}$ are fused with the single embeddings using a one-layer MLP, following the same practice in RDE. A max-pooling layer is applied to the output of IPA module to attain global structure representations for the wild-type and mutant complexes, denoted as $\boldsymbol{H}^{wt}$ and $\boldsymbol{H}^{mt}$ respectively. Finally, these global representations are fed into another MLP to predict $\Delta\Delta G$. The MSE loss is used for training.