Explicit versus implicit consideration of binding partners in protein–protein complex to elucidate intrinsic dynamics

The binding of many proteins to their protein partners is tightly regulated via control of their relative intrinsic dynamics during the binding process, a phenomenon which can in turn be modulated. Therefore, investigating the intrinsic dynamics of proteins is necessary to understand function in a comprehensive way. By intrinsic dynamics herein, we principally refer to the vibrational signature of a protein molecule popularly obtained from normal modes or essential modes. For normal modes, one often considers that the molecule under investigation is a collection of springs in a solvent-free or implicit-solvent medium. In the context of a protein-binding partner, the analysis of vibration of the target protein is often complicated due to molecular interaction within the complex. Generally, it is assumed that the isolated bound conformation of the target protein captures the implicit effect of the binding partner on the intrinsic dynamics, therefore suggesting that any influence of the partner molecule is also already integrated. Such an assumption allows large-scale studies of the conservation of protein flexibility. However, in cases where a partner protein directly influences the vibration of the target via critical contacts at the protein-protein interface, the above assumption falls short of providing a detailed view. In this review article, we discuss the implications of considering the dynamics of a protein in a protein-protein complex, as modelled implicitly and explicitly with methods dependent on elastic network models. We further propose how such an explicit consideration can be applied to understand critical protein-protein contacts that can be targeted in future studies.


Introduction
Proteins are the basic building blocks of life, and by associating form large macromolecular machinery to perform a wide array of complex processes such as DNA replication, transcription, translation, and signal transduction etc. (Alberts et al. 2002). Protein-protein interactions (PPIs) control biological activities in a wide variety of settings, emerging as a complicated network of interactions as part of the "interactome" (Ardini et al. 2022). Aberrant PPIs are associated with many diseases and thus act as targets for therapeutics while advances in protein engineering to develop supramolecular assemblies in the fields of biomedicine and nanotechnology have also shown a need to understand multi-subunit dynamics for more effective rational design (Lu et al. 2020).
Characterizing protein-protein interfaces based on amino acid profiles has provided great insight into the identification of such potential hotspots. Since the early description of obligate and non-obligate interfaces in dimers (Jones and Thornton 1996) and the comprehensive structural analysis of a diverse set of protein-protein complexes (Conte et al. 1999), we have had an explosion in structural data, particularly from cryo-EM that have provided a more sophisticated picture of the amino acid profiles and structural characteristics of protein-protein binding (for a detailed historical perspective, refer to Kastritis and Bonvin (2013)). Molecular modeling and simulations have been instrumental in uncovering the dynamic properties of protein-protein interactions, laying the foundation for the CAPRI community-wide experiment to assess protein-docking methods that predict protein-protein interactions for the past 20 years (Janin 2002).
Characterizing dynamics has become an integral part of understanding protein function. Out of many established methods, elastic network models (ENMs) have been instrumental in probing the dynamics of protein-protein interactions for a number of studies, through structural perturbation methods (Zheng et al. 2007;Thirumalai and Hyeon 2018), and the normal mode analysis (NMA) of oligomers and large complexes (Marcos et al. 2014;Mishra et al. 2017). In 2013 and 2014, Dr. Bhaskar Dasgupta, Dr. Akira R. Kinjo and Prof. Haruki Nakamura observed that a counterbalance of different types of dynamics, and more specifically, the relative rigid-body motions of partner molecules, dominates protein-protein binding interactions in nine different Ubiquitin-based complexes through all-atom ENMs (Dasgupta et al. 2013(Dasgupta et al. , 2014b(Dasgupta et al. , 2020. More importantly, they found that the ligand-bound dynamics is encoded in the apo-form of Ubiquitin, exposing the power of capturing dynamics specific to its conformational state. In this review, we look at the developments made within ENMs to address protein-protein interactions, how such analysis has been performed implicitly, and what we can gain from explicit modeling. Specifically, we discuss the all-atom ENM implementation by Dr. Bhaskar Dasgupta, Dr. Akira R Kinjo andProf. Haruki Nakamura (2013, 2014b) to understand the intrinsic dynamic coupling between protein binding partners without considering specific chemical interactions.

Intrinsic dynamics of protein include functionally relevant motions
Since Tirion's original formulation (Tirion 1996), ENMs have been shown to capture functionally relevant motions in many systems (Dykeman and Sankey 2010). The lowest frequency (or energy) normal modes capture mechanical components of conformational change, whether by large amplitude motions or the propagation of allosteric signals (Tama et al. 2002;Tama 2003;Reuter et al. 2003;Gerek and Ozkan 2011;Rodgers et al. 2013;Mahajan and Sanejouand 2015;Fuglebakk et al. 2015;McLeish et al. 2015;Atilgan 2018). The success of the ENM can be elegantly explained by Halle's observation that much of the thermal fluctuations as captured by X-ray crystallography can be explained as a feature of local protein packing densities (Halle 2002). When well-parameterized against experimental data or more detailed simulations, the ENM reproduces the influence of both protein packing density and anisotropic effects (Poon et al. 2007;Hinsen 2008;Yang and Chng 2008;Fuglebakk et al. 2013, Zhou andLiu 2014;Koehl et al. 2021;Dubanevics and McLeish 2022).
Many studies have strengthened the idea that dynamics is more highly conserved than structural topology or amino acid sequence. This has given rise to the idea that much of protein function can be tweaked by mutations in vital regions of deformation or by interactions with other proteins or effectors. For instance, other studies have suggested that such patterns of intrinsic dynamics have influenced the evolution of protein structure, whereas structural topology is also an important determinant of intrinsic dynamics conservation regardless of function (Hollup et al. 2011;Fenwick et al. 2014;Reuter 2016, 2018). Similarly, the intrinsic dynamics of proteins have an influence on their interactions with partner molecules, giving rise to large complex machinery and structural components of the cell, to gain functionally relevant mobility.
As ENMs developed, efforts were made to reduce the size of the Hessian matrix to gain computing efficiency. The most often used ENM based on connecting C-alpha atoms was first introduced by Hinsen (1998), giving rise to C-alpha ENMs with a variety of cut-offs and benchmarking methods (Atilgan et al. 2001;Moritsugu and Smith 2007). The performance of many of these ENMs has been evaluated based on their agreement with MD simulations (Fuglebakk et al. 2013). Other more sophisticated ENMs have been developed: rigid backbone-enhanced (Srivastava et al. 2012;Dubanevics and McLeish 2022), as well as an ENM that uses machine learning to differentiate breaking from maintained contacts (Putz and Brock 2017). As an alternative to transforming 3D Cartesian coordinates, rotationally invariant internal coordinates such as the dihedral angle have also been explored (Kitao and Gō 1991;Wako et al. 2004;Mendez and Bastolla 2010;Lopéz-Blanco et al. 2011).

Normal mode analysis captures dynamics of a protein structure around its native conformation
Following various established formalisms (Wilson et al. 1955;Field 2007;Fuglebakk et al. 2013;Dasgupta et al. 2014a), we describe how NMA is performed using ENMs as a general method in this section. NMA captures vibrational signature of a given protein structure, assuming the harmonic approximation. In this method, the potential energy of the system is expressed as a Taylor series of the displacement of atomic coordinates truncated to the second differential term.
For a system of N atoms, the ith atom has coordinate (x i , y i , z i ), and corresponding displacements are given as x i , Δy i , and Δz i . The coordinate system can be defined as q k = Δx i , q k+1 = Δy i , q k+2 = Δz i where k = 1, 2, 3, …, 3 N-2, 3 N-1, 3 N describing a total of 3 N degrees of freedom.
The potential for the above-described coordinate system can take the following form, where the equilibrium condition can be imposed by limiting q, around the minimum potential energy dU∕dq i q=0 = 0 . We can further simplify the above equation by setting first term in the left-hand side to zero, giving, The potential energy curve is a parabolic function centered around the equilibrium position of the system. The partial derivatives in (2) constitute a Hessian matrix (H) of the size 3 N × 3 N, where each term is given by H ij = d 2 U∕dq i dq j q=0 . The time-evolution of the system given in (2) using Newton's equation of motion is given by, for each j, and j = 1, 2, 3, …, 3 N and m j is the mass of atom j. Thus, the Hessian matrix is constituted by 3 N simultaneous second-order differential equations, where a coupled harmonic oscillator is constructed for each degree of freedom. To solve the second-order differential equations, the following general form can be used, where A j is the amplitude of the jth oscillator with frequency ω j with phase ϕ j . By combining (3) and (4), the secular equation is as follows, where H is mass-weighted hessian matrix where each term ij is Kronecker delta-function, which is 1 only when i = j, and otherwise 0.
where each eigenvector has 3 N dimensions. By diagonalizing the Hessian matrix (Eq. (5)), the displacement of atom i from equilibrium position is obtained as a sum of contribution from individual oscillators, where K B is Boltzmann constant; T is temperature. The superscript "(0)" indicates the native conformation. The effective eigenvectors run from 1 to (3 N-6), where the first six are related to six rigid-body motions of the system (three translational and three rotational). Using Eq. (6), the covariance matrix (3 N × 3 N) can be calculated via the following, where v k,i indicates ith component of kth eigenvector, M is diagonal mass-matrix with k, k + 1 and k + 2 diagonal elements from atom i as 1∕ √ m i , and Λ g is diagonal matrix with 3 N-6 non-zero eigenvalues. The superscript "T" indicates the transpose operation.
For standard NMA, the system of atoms should be minimized so that the whole system reaches a minimum energy conformation prior to solving Eq. (5). However, for a reasonably large system, it is difficult to ensure such minimum energy conformation by energy minimization techniques, due to the complexity of the potential energy function in molecular mechanics force fields. A reasonable approximation is derived by assuming that experimentally determined coordinates capture a low-energy conformation (a local minima, indicating native conformation), and thus, the Taylor series expresses the potential energy at the energetic minimum, avoiding costly energy minimization.
Still, obtaining each Hessian matrix term remains arduous. This can again be bypassed effectively by using ENMs. In ENMs, the harmonic oscillator is envisioned as a Hookean spring (Tirion 1996;Fuglebakk et al. 2013), fluctuating around a native conformation. To define a spring, the distance between a pair of atoms i and j, d ij is used, for which a spring between atom i and j has energy, where c ij is force constant of the spring, and d (0) ij is distance between i and j in native conformation. The total potential energy is sum of individual Hookean oscillator, where K is a phenomenological constant. First note that, to calculate each Hessian term, d ij should be written in terms o f i n d i v i d u a l c o o r d i n a t e d i f fe r e n c e s , l i ke , d 2 ij = Δx ij 2 + Δy ij 2 + Δz ij 2 . Then, double differential of U is calculated with respect to change in coordinate in two different directions (say α and β where , ∈ {x, y, z} ) for two atoms i and j. This constitutes mass-unweighted Hessian terms, where A popular way to define the pairwise force constant is, where d cutoff is a cutoff distance specified by user. This also referred as distant-dependent force constant. Note that, for a system with N atoms, the Hessian matrix is 3 N × 3 N dimensional and thus, for large systems, diagonalization of Hessian is a memory-intensive procedure. There are several ways to increase computational efficiency without sacrificing significantly on accuracy: 1) A popular coarse-graining scheme is to define the elastic model based on only Cα-atoms. This reduces the computational burden by considerable extent. Moreover, in this case, all constituent atoms are identical so that massweighting of Hessian is redundant (Bahar and Rader 2005). 2) Another approach related to decreasing burden of diagonalization is rotational-translational-block (RTB) approach (Durand et al. 1994;Tama and (8)  . Here, the Hessian matrix is projected to a subspace, which is then diagonalized. The projection matrix is designed to have a dimension of 3 N × 6 M, where M is the number of blocks of the system. After projection, the Hessian matrix is 6 M × 6 M in dimension, which is much simpler as M is much less than N. If M is equal to number of residues, then each residue is included in a block, and RTB and Cα-based coarse-graining are of similar size. In such a condition, one could argue that Cα-based method treats the effect of other atoms in each residue implicitly.
3) The normal mode of protein structure can be obtained in dihedral angle space (Noguti and Gō 1983;Ishida et al. 1998;Lopéz-Blanco et al. 2011). For this, the above potential energy expression (1) needs to be changed to internal coordinate space, i.e., degrees of freedom of the system are defined in terms of bond lengths, bond angles, and dihedral angles. Then, one can assume that bond lengths and bond angles do not change in functionally relevant motions, and those degrees of freedoms are fixed. Thus, only dihedral angles are allowed to fluctuate. One limitation of this method is that the normal mode displacements (as in Eq. (6)) is defined in terms of dihedral angle, requiring a back-transformation to the usual Cartesian coordinate system (Wako and Endo 2011).
An important aspect of the covariance matrix (Eq. (7)) is that when diagonalized, it returns the normal modes of the system with the reciprocal of the eigenvalues as the approximate normal mode frequencies. This acts as an equivalent method to characterize intrinsic dynamics of a system via molecular dynamics (MD) simulation (Amadei et al. 1993). An MD simulation obtains a set of conformations of a protein, from which the covariance matrix ( C ij,MD ) can deduced from ensemble average, given by where N c is number of conformations in the ensemble, C ij,MD is i,jth element of the MD-based covariance matrix, x ik is ith coordinate out of 3 N degrees of freedom from kth conformation, and < x i > represents average of ith coordinate over all the conformation. For MD, diagonalization of the above covariance matrix gives us principal components (eigenvalue) as variance, where each eigenvalue represents magnitude of variation of the data along the principal eigenvector. The eigenvector includes coefficients for the linear combination of observed structures (often referred as loadings for the principal components).

Normal mode analysis for protein-protein complex
For a single protein chain, the above method of obtaining normal mode analysis is straightforward. However, to consider the interplay of intrinsic dynamics in a multi-chain assembly requires careful treatment of the normal mode calculation. In this section, we use the CBL-b and Ubiquitin (Ub) complex as an example system, as originally used in by Dasgupta et al. (Dasgupta et al. 2013, 2014b. (Fig. 1a).
Generally, for a dimer with chains A and B, there are rigid-body roto-translational motions between A and B.
Here, A and B include N A and N B atoms, respectively, giving rise to a system of N A + N B = N atoms. The number of degrees of freedom for the A-B dimer is thus 3 N-6, where 6 rigid-body motions of the dimer are subtracted. If treated individually, A and B chains include 3N A -6 and 3N B -6 intrinsic vibrational motions, which sum to 3(N A + N B ) -12 = 3 N-12 degrees of freedom. Therefore, when applying NMA on AB dimer, we retrieve an extra six degrees of freedom in the system, which represent the relative rotation and translation between A and B. We refer to these extra six degrees of freedom as external motions between A and B, which obscure the vibrational signatures present in the A-B dimer. Fortunately, there is a way to exclude these external motions by applying a suitable projection matrix to A and B, referred to as P A and P B , respectively (Field 2007;Dasgupta et al. 2014b).
Considering the covariance matrix obtained from usual NMA analysis of A-B, C complex , which have the following blocks (Fig. 1c, for the example case of CBL-b and Ub dimer, where chain A is CBL-b and chain B is Ub), the C AA part (or C BB part) is a symmetric square matrix that corresponds to the coordinates of chain A (or B) along its rows and columns, while the C AB part is a non-symmetric rectangular block that corresponds to the coordinates of A along its rows and coordinates of B along its columns. Thus, when projecting C AA to a subspace without the six external degrees of freedom between A and B, we obtain, and similarly, when projecting C AB to a subspace without the six external degrees of freedom between A and B, we obtain, The matrix C AA,int includes motions of A within the dimer which are due to the intrinsic flexibility of A itself, and those Fig. 1 a The dimeric system of CBL-b (Cb) and Ubiquitin (Ub) from the PDB ID 2OOB as an illustrative example. The alpha-helical and beta-strand secondary structures are numbered and annotated as "H" and "S," respectively. b Intrinsic dynamic characterization of a AB heterodimer can be approached with explicit and implicit consideration of partner. From the explicit treatment, we obtain self-and directly coupled motion of A and B. From the implicit treatment of the partner, we also obtain motions of A and B. c Correlation matrix obtained from the normal modes of AB heterodimer can be split into Ub and Cb, from which self-coupled motion can be obtained. The remaining rectangular part of the matrix, "CBL-b and Ub direct coupled," includes directly coupled motions. To obtain vibrational signatures from self-and directly coupled submatrices, projection matrices, annotated with "P" and subscript corresponding to the component, are used motions are referred to as "self-coupled motion". The matrix C AB,int is a rectangular matrix, from which one can obtain motion of A (and B) by performing singular value decomposition. The corresponding singular vectors of A include dynamical features of A due to the direct influence of B. We refer to this as "directly coupled motion" between A and B or as "ligand-coupled motion"-envisioning B as the ligand or partner of A. Here, note that self-coupled motion of A includes dynamics of A by itself, where B is present explicitly. The directly coupled motion of A includes cross-talk between A and B. Therefore, for a dimeric system, we obtain self-coupled motions from each component and directly coupled motion between two components. This classification is shown in Fig. 1b.
For a complex, perhaps the simplest process to the extract intrinsic dynamics of a subunit or chain (for example, A or B) from a dimer is to perform NMA on it after isolating its structural coordinates from its binding partner. By this process, we use the bound or holo-conformation of the component, where the influence of the partner molecule is present implicitly. Such an implicit consideration sometimes useful to discuss about the functional dynamics of the protein without great computational cost due to system size or the additional computation of projecting multiple subspaces, as performed in the reference (Perica et al. 2014). This implicitly considered dynamics is also included in the classification shown in Fig. 1b. How the different types of motions in a complex, viz. dynamics of A in implicit or explicit presence of B, and dynamics of A due to direct cross-talk with B, are related is complicated. Nevertheless, it is useful to understand the interaction between A and B. Below such salient dynamical features of a complex is demonstrated by the intrinsic dynamics of the ubiquitin and CBL-b protein complex (Fig. 1a).

Recovering self-coupled motion of components within a dimer
To demonstrate how different types of motion of components within a protein-protein complex shapes its dynamics, we use a dimer of ubiquitin associated domain (UBA) of CBL-b protein and ubiquitin (Peschard et al. 2007). CBL-b is a small ubiquitin ligase within CBL family and ubiquitinates the activated form of receptor tyrosine kinases to negatively regulate them. The CBL family of proteins includes a highly conserved N-terminal domain by which it binds to tyrosine. In the C-terminal region, CBL-b and c-CBL include the UBA domain, which is a tightly packed three-helix bundle with a hydrophobic surface mediating the interaction with ubiquitin. Here, we specifically analyze CBL-b C-terminal UBA domain (which we refer henceforth as Cb) complexed with ubiquitin (or simply Ub) (PDB ID: 2OOB, Fig. 1a) in three ways-(i) we recover dynamics of CBL-b UBA domain with the explicit and (ii) implicit presence of ubiquitin, as well as (iii) analyze dynamical cross-talk between CBL-b UBA domain and ubiquitin.
To obtain the intrinsic dynamics of Cb with the explicit presence of Ub, all heavy-atoms in the complex (N Cb = 350, N Ub = 574, N = 924) are included in the initial NMA to retrieve the full covariance matrix (C complex ) of the whole complex. From the full covariance matrix, the symmetric block matrix related to Cb, C Cb,Cb (3N Cb × 3N Cb ), is projected by the matrix P Cb , to obtain C Cb,Cb,int (see Eq. 13), which includes the vibrational signature of Cb with the explicit presence of Ub. By diagonalizing C Cb,Cb,int , we obtain the normal modes of vibration of Cb which are referred to as self-coupled motion of Cb (in total 1044 modes for 350 atomic system). For the example illustrated here, the 100 lowest frequency modes were used to reconstruct the selfcoupled covariance matrix by Eq. (7), where the sum runs from 1 to 100 modes.
The reconstructed covariance matrix ( C 100 Cb,Cb,int ) represents a reduced space representation of self-coupled motion of Cb and subsequently a correlation matrix is obtained from covariance matrix by normalizing off-diagonal elements by diagonal elements (Fig. 2a). In a similar manner, the Ub-specific part of the covariance matrix, C Ub,Ub (3N Ub × 3N Ub ) is projected by P Ub to obtain C Ub,Ub,int , and the following diagonalization provides 1716 vibrational self-coupled modes of Ub. The correlation matrix corresponding to Cα-atoms in Ub is obtained from the covariance matrix based on the 100 lowest energy normal modes, C 100 Ub,Ub,int , representing a reduced space representation of self-coupled motion in Ub (Fig. 2b).
The high positive correlation between Cb Cα-atoms in Fig. 2a demonstrates that vibrational motion relay along the alpha-helices. The alpha-helical pairs H1-H2 and H2-H3 move in sync in the same direction. Such collective motion of Cb in the explicit presence of Ub involve residues of Cb both near and far from the interface. In Fig. 2b, high positive correlations are observed between beta-strand pairs S1-S2, S3-S4, and S1-S4 in Ubiquitin which face Cb. In other words, Ub self-coupled motions are concentrated in the secondary structure element close to Cb, but the alpha-helices in Ubiquitin positioned at the distal surface remain relatively rigid. When comparing between Cb and Ub, the latter is more rigid overall in terms of its intrinsic flexibility. Thus, when examining self-coupled motion, we can find collective motions related to both the binding interface and distal regions as encoded by structural topology.

Recovering cross-talks between two components in a dimer
The dynamical features of cross-talk between two interacting protein subunits in a complex are recovered by analyzing the rectangular part, C AB, of the full covariance matrix (Fig. 1c). An element in this part of the covariance matrix includes the joint variability of two atoms from two different subunits. Such a covariance exists because in the preceding calculation of Hessian matrix, the interaction between atoms from the two subunits was included explicitly. However, the intrinsic dynamics within C AB includes not only the vibrational signature of the components but also the relative rotations and translations between them. Previously, we showed that such relative rotations and translations depend on the shape of the interacting molecules (Dasgupta et al. 2014b). To filter out these relative rigid-body motions, C AB needs to be projected by projection matrices involving both subunits (see Eq. 13b), following which singular value decomposition can be performed to obtain vibrational modes for both subunits (Fig. 3a). It is common to identify a few low-frequency normal modes (e.g.100 normal modes considered here) to reconstruct a covariance matrix that excludes highfrequency noisy motions.
In the example of Cb and Ub, we find some specific vibrational signatures that are directly linked to the crosstalk between them, which may be lost without explicit consideration. Here, when comparing the individual correlation matrices of Cb and Ub constructed from directly coupled normal modes (Fig. 3b, c), we observe that in Cb, both positive and negative correlations are easily identified, whereas for Ub, mostly positive correlations exist. Therefore, for directly coupled motion in Cb, a few clusters of atoms move in the same direction; a few residues within alpha-helices H1 and H3 (and in alpha-helices H2 and H3) move in the same direction (white circular regions in Fig. 3b), while such strong positive correlations do not exist between alpha-helices H1 and H2. The terminal residues are strongly anti-correlated with alpha-helices H2 and H3 (black rectangular regions in Fig. 3b).
In Fig. 3c, the correlation matrix from Ub directly coupled motion shows a contrasting set of positive and negative correlations as compared to the correlation matrix constructed from its self-coupled motion (Fig. 2b). In the Ub directly coupled motion, negatively correlated atoms mostly involve loop regions and helix termini residues, while positively correlated atoms correspond to betastrand residues, identical to the self-coupled profile. The apparent rigidity observed in self-coupled motion (Fig. 2b) is not present here as several groups of atoms move in sync, either in the same or opposite direction. These observations suggest that Ub intrinsic dynamics, due to interaction with CBL-b and vice versa, connect C-alpha atoms that are far apart in sequence. The contrast in negatively and positively correlated intrinsic dynamics in self-coupled and directly coupled motion of CBL-b and Ub is further demonstrated in Fig. 3d.

Specific contact signatures from explicit consideration of binding partners
The intrinsic dynamics captured in the directly coupled normal modes are related to molecular recognition in Cb to Ub. Focusing on the specific set of residues that contribute to the directly coupled correlation matrix of Cb to Ub and Ub to Cb (using first 100 normal modes), we can separate the regions based on positively and negatively correlated atoms (Fig. 2). Through more sophisticated methods, such as by setting a threshold on the positive correlation and performing clustering, we can identify potentially key-interacting regions. For example, in Fig. 4a and c, we used popular clustering method DBSCAN (or density-based spatial clustering of applications with noise) on correlation values to identify important interacting regions in the complex. The cluster identifiers and the residues comprising the clusters are shown in the structure (Fig. 4b, d).
Due to the collective nature of the intrinsic dynamics, the interacting regions identified are not confined to the binding interface between Cb and Ub. Therefore, even though the fluctuations at the binding interface are high due to intermolecular cross-talk, the motions propagate to other parts of the structure, to allow different parts of the whole polypeptide chain to move in sync. Therefore, a dynamic view of intermolecular cross-talk obtained by explicitly considering partner molecule enables us to probe conformational change far distant from the interface, which is important for targeting binding interactions allosterically.

Difference between explicit and implicit consideration of binding partners
We contrast the analysis in the previous section with the implicit consideration of presence of Ub in Cb dynamics (and similarly presence of Cb in Ub dynamics) using correlation matrices obtained from 100 lowest energy normal modes (Fig. 5a, b). In general, we find that the correlation values are more positive and contain less contrast between regions compared to the correlation matrices from the directly coupled normal modes (Fig. 6b). For both Cb (Fig. 5a) and Ub (Fig. 5b), only a few Cα atoms are positively correlated, while all other pairwise correlations are moderately positive. Therefore, these results show fewer regions moving in sync in each subunit under bound conditions, and that only highly significant correlations can be retrieved. Furthermore, the global similarity between implicitly coupled and self-coupled normal modes is evident from Table 1, which show the Bhattacharyya coefficients (BCs) between all-atom covariance matrices (Fuglebakk et al. 2012). The BCs indicate that directly coupled motions are distinct from self-coupled motions derived from explicit or implicit consideration of the binding partner. Between Cb and Ub, the BCs are higher in the case of Ub, which may be related to more rigid nature of Ubiquitin.
A very common measure, the root-mean-square fluctuations (RMSF) of Cα-atoms, gives use a per-residue view of the intrinsic dynamics of Cb and Ub (Fig. 6). The expected dampening of the fluctuations in the self-coupled and implicitly coupled dynamics of Cb relative to directly coupled dynamics indicates that the implicit confirmation encodes the presence of the binding partner. This explains the success of studies that rely on the initial conformation of individual subunits when studying oligomers (e.g., Perica et al. (2014)). In the Cb and Ub example, the self-coupled dynamical fluctuation pattern embedded in RMSF is strongly correlated to implicitly coupled motion. The fluctuations of directly coupled normal modes of Cb is strongly influenced by interfacial residues, as those residues show higher RMSF. Fig. 6 a Root-mean-square fluctuations (RMSFs) averaged over all residues for Cb are compared for directly coupled, self-coupled, and implicitly coupled motion. In the Pearson correlation (inset) between RMSF values from different types of motion are also compared. b The ratio of interface to non-interface average atomic fluctuations (RMSFs) are compared for directly coupled, self-coupled, and implicitly coupled motion. c For Cb (in magenta) interacting with Ub (in cyan), the residues that show high fluctuations in self-coupled motion are shown. Glu947 of Cb has high RMSF values in all three types of motion and is positioned close to the binding interface. d Cb-Ub dimeric interface shown in the context of apo Cb (2OOA) oligomeric complex (as predicted by PISA) (top). At the bottom of the structure, Arg951 and Gln957 are shown oriented towards a hypothetical interface between two Cb chains in a homo-tetramer There are a few residues showing higher than average fluctuation in implicitly coupled motion and also in self-coupled motion, which are exposed on the surface upon unbinding, e.g., Arg951, Glu954, Gln957, Asn959, and Glu961 (marked in red in Fig. 6c). The residue Glu947 shows high fluctuation in all types of motion and is positioned closer to the interface between Cb and Ub, near the "rim" region of the interface (Chakrabarti and Janin 2002).
Note that the residues that show higher than average fluctuation in the implicitly coupled motion of Cb are either interfacial residues (relatively damped but still higher) or residues positioned at the opposite end, distal to the interface (much higher). The packing of interfacial residues is expected to be directly impacted by the presence of a binding partner, which may account for the dampening of the fluctuations in an otherwise exposed region of the surface. In this way, plausible interface residues can be predicted directly from NMA without explicitly accounting for a binding partner. In fact, the apo-form of Cb (PDB ID 2OOA) is present in a homo-tetramer (as predicted by PISA (Krissinel and Henrick 2007)), where the tetrameric interface between Cb chains would be positioned at the opposite end of the heterodimeric interface between Cb and Ub (Fig. 6d,  top). Again, the implicitly derived self-coupled motion of Cb capture the profile of the residues at the binding interfaces within homo-tetramers (Fig. 6d, bottom).

Future perspectives
Elastic network models (ENMs) have been successful at uncovering the intrinsic flexibility of proteins as a single chain or cohesive complex assembly despite its simplicity. When considering all atoms, computational cost rise quickly due to the complexity of treating the many degrees of freedom conferred by individual chains in a complex. Furthermore, extracting intrinsic dynamics is complicated due to the complexity of the structure, which is heterogeneously packed in space. Here, we have discussed a more sophisticated method to uncover salient dynamical characteristics related to protein-protein interactions in multisubunit complexes using ENMs. Notably, Dr. Bhaskar Dasgupta, Dr. Akira R. Kinjo and Prof. Haruki Nakamura (Dasgupta et al. (2013(Dasgupta et al. ( , 2014b) previously implemented an ENM-based projection method to interpret specificity in intrinsic dynamics for Ub binding to an array of partner molecules, and showed that the self-coupled and directly coupled motion of Ub balance out each other on average. In a related study, Dasgupta et al. showed that the relative rigid-body motion, which depend on the shape of the interacting molecules, dictate directly coupled motion (Dasgupta et al. 2014b). However, the vibrational signatures present in directly coupled motion for a complex were not analyzed in detail. As an illustrative example to expose the power of interrogating protein-protein intrinsic dynamics features, we discussed the vibrational signatures present in directly coupled motion for the Ub and Cb complex.
The self-coupled motion is an important concept discussed in this review. The self-coupled motion of Ub (or Cb) is the motion intrinsic to the protein in the explicit presence of the partner. Different from approaches that aim to integrate out the dynamical characteristics of the partner by applying a projection based on directly coupled part of the Hessian matrix, as performed in some studies (Zheng and Brooks 2005;Li et al. 2020), we compare self-coupled motion with implicitly coupled motion to understand the effect of considering the partner molecule explicitly. Additionally, it was previously shown that the apo-Ubiquitin motion include some directly coupled motion, indicating that the recognition of partners by Ub is embedded in implicit motion of Ub or that Ub is intrinsically less flexible. Explicit consideration of a binding partner is useful for dispelling such dilemma.
The ENM is an interconnected network of springs, which treats the protein molecule as a solid substance, while in physiological conditions, it exists in solution, affecting its packing density and overall behavior (Liang and Dill 2001). An additional application of the ENM-based method was previously outlined by Dasgupta et al. (2014a) where a protein molecule is treated with more liquid-like properties and the effect of solvation was treated implicitly. The potential for implicitly treating non-protein partner molecules, such as water, to the intrinsic dynamics of proteins increases the potential for deeper insight into the mechanical influence of solvent on proteins. Comparing such methods with others that treat the effect of solvent implicitly could inform us further of how phenomena such as allosteric communication between key residues occur and should be investigated in the future.
Another fundamental area of further study would be to analyze directly coupled motion for a protein-protein interface related to a rough binding surface. A rough binding surface may include some atoms which are not directly contacting the partner molecule, i.e., having no change in accessible surface area due to binding, but are positioned within the close proximity of the binding partner (Dasgupta et al. 2011). Such interface roughness potentially affects protein-protein recognition, allowing local structural changes to affect protein-protein recognition.
For the immediate future, a practical extension of the current algorithm would be to study effect of direct-coupling and self-coupling in oligomeric assembly. For example, tetrameric systems such as the ones from the PyrR protein family (Perica et al. 2014) which already showed conservation of intrinsic dynamics pertaining to oligomeric state (Tiwari and Reuter 2016) can be studied in greater detail when explicitly treating partnering subunits and allosteric effectors (which shift the oligomeric state). In the Ub-Cb example, we observed that the fluctuations profiles from self-coupled and implicitly coupled motion can identify residues which form the binding interface and this observation requires meticulous investigation in higher-order complexes.
A striking observation in this review is that directly coupled motion for Cb are not confined to the binding interface, but also to residues positioned distally. This could indicate the presence of allostery, caused by changes to surface residues. In a method where the partner molecule is explicitly considered, the allosteric effect observed is directly related to perturbation at the interface. By calculating directly coupled motion via the intrinsic flexibility of the protein, we can design molecules that target the regulation of protein function through remote-control.
ENMs lend themselves well to coarse-graining and parameterization at the user's discretion. This feature complements well to fine-tuning methodology that investigates the dynamics of larger biological assemblies. With the rapid increase in data, especially from bioimaging (cryogenic electron microscopy and tomography), the need for hybrid approaches for investigative complex biological systems is ever-increasing (Ognjenović et al. 2019;Srivastava et al. 2020;Liao et al. 2022;He et al. 2022). At the same time, we envision that increasing computational capacities will allow for constructing larger Hessian matrices and largerscale analysis. Since the method outlined here depends on proper calculations of covariances and projections, it is also transferrable to analyzing molecular dynamics simulations. ENM-based analysis could then provide the basis for such an extension to MD simulations that will allow protein dynamics to be simulated in a time-dependent manner, and the effects of additional physicochemical parameters on directly coupled motion to be probed.
Funding SPT was supported by the HIRAKU consortium start-up funding.

Data availability
The raw data supporting the conclusions of this article will be made available by the authors upon request.