A new approach for extracting information from protein dynamics

We retrieved crystal structure for FimH_L wild type and mutant, Siglec-8 apo and holo, and SARS-CoV-2 spike protein from the Protein Data Bank, as detailed in Table 1. For FimH_L, we used crystal structures of the lectin domain without ligand for both the wild type and the mutant. To compare dynamics with and without the disulfide bond as a local perturbation, we used Visual Molecular Dynamics (VMD) to define the bond or two cystines for FimH_L. For Siglec-8, we simulated the ligand 6’S sLe^x without the 3-amino-propyl linker, which is not thought to interact with the binding pocket ²⁹.

For the SARS-CoV-2 spike protein, we started from the refined structure on the CHARMM-GUI archive ³¹, ³⁰. To focus on one hinge system that is thought to be different between the down and up states, we isolated the receptor binding domain (RBD) and subdomain 1 (SD1) protein subunits without the glycans. For the up state, we used chain A where the RBD is accessible for binding the ACE2 receptor on human cells ³⁰, and for the down state, we used chain B. We used the structure

We prepared all systems using VMD version 1.9.3 ³². We solvated each protein with at least 16 Å of TIP3P water molecules on each side to prevent interactions with itself through the periodic boundary conditions. We added sodium and chloride ions to neutralize the system and achieve the desired salt concentration in Table 1.

We performed all-atomistic MD simulations using NAMD ³³, with the CHARMM force field³⁴. Our NAMD simulation parameters and system details are listed in Table 2.

After observing differences in correlated protein motions between replicates, we performed three replicates of over 200 ns each for wild type FimH_L. Due to the tradeoff between the number replicates and simulation length, we performed six replicates of 20 ns of FimH_L to make comparisons of wild type and mutant FimH, as well as wild type FimH with and without the Cys3-Cys44 disulfide bond.

Since twenty lowest-energy structures are reported for Siglec-8, we performed a single replicate of 50 ns for each structure to compare apo and holo Siglec-8. For the spike RBD-SD1 domains, we performed six replicates of 60 ns.

We used dihedral angles to capture protein dynamics because dihedral angles identify localized regions responsible for the collective displacement of regions distal from the angular rotation, such as in hinged motion ³⁵. Dihedral angles are also an internal coordinate system that avoids the structure alignment step when using Cartesian coordinates, which can introduce artifacts ²⁷. We use both backbone ($\phi,\psi$) and sidechain ($\chi_{1}-\chi_{5}$) dihedral angles. We extracted protein features with MDTraj 1.9.5.

In the literature, the covariance matrix is one approach used to identify protein regions with motions that are related to the motions of many other regions; in particular it is used to identify correlated motions between distant regions in allostery ²², ³⁶. However, constructing networks from the covariance matrix, even with a threshold to remove weak correlations, is susceptible to induced correlations when two nodes (e.g. A and C) are not directly connected but share a connection with a third node (B) ¹⁴. Borrowing from the field of network reconstruction, we use the inverse of the covariance matrix to identify the connections and weights, or edges, between nodes ¹⁴. This approach is consistent with finding the inverse of a covariance matrix based on C$\alpha$ positions, which fits the Hessian matrix describing an elastic spring network with anisotropy ²⁵, ²⁶. We have found that the anisotropic elastic network model (ANM) has large errors when used to describe the motion of FimH, which is consistent with errors for hinge-motion described in literature ³⁷. As a result, we use dihedral angles. This approach is similar to the torsional network model (TNM) which uses equal spring constants to describe dihedral angles across the protein ³⁸. In contrast, the inverse of the covariance matrix uses the variances and the covariances of angles to calculate spring constants for a network of torsional springs. The nodes in our network are dihedral angles, the edges are like linearly coupled torsional springs, and the inverse of the covariance matrix is the Hessian matrix for a TNM.

To construct our network, we calculate the Moore-Penrose pseudo-inverse of the covariance matrix using both the backbone and sidechain dihedral angles. We use this approach to understand the relative contributions of backbone and sidechain dynamics to collective motion. Since the sign describes whether the angles turn in the same direction, we take the absolute value to get the interaction strength. We do not apply distance filters. While we use the 97^th percentile as a value threshold for selecting strong interactions or visualizing the network on the protein, we do not use any thresholds for network comparisons. More generally, we recommend caution for applying thresholds to these networks for analysis.

To identify interactions that are stronger in one protein state than another, we compare each edge. We select for large differences between groups, relative to the variability within each group. To do this, we filter for differences larger than twice the standard deviation for each group. To compare an edge $e$ between states $a$ and $b$, each with an ensemble of $m$ and $n$ networks, this is $|\langle e_{a}\rangle_{m}-\langle e_{b}\rangle_{n}|>2\sigma_{e_{a}}$ and $>\sigma_{e_{b}}$. We apply this rule without determining statistical significance with corrections for multiple comparisons, in order to see the full effects of comparing all interactions on the matrix. In Supplementary Figure 11, we also show an example of comparisons without filtering for large differences, in order to illustrate the persistence of the contact-map pattern. We perform the network comparisons in two ways: 1) for every edge on the network (see Supplementary Figure 11), 2) accounting for the multi-layer structure of the network by collapsing the backbone-backbone interactions into residue-residue interactions (Figure 4). We analyzed data in python, using SciPy (RRID:SCR_008058) and custom packages.

Our approach for constructing a network representation of the dynamics of a given protein is comprised of three steps. In the first step, we obtain temporal dynamics for the nodes, which are the backbone ($\phi$, $\psi$) and sidechain ($\chi_{1-5}$) dihedral angles for each residue ²⁴, ³⁵. In the second step, we calculate the circular covariance for dihedral angles ³⁵, which can be thought of as the linearization of the interactions captured by mutual information. In the third step, we invert the covariance matrix using the Moore-Penrose pseudo-inverse to calculate the best fit for a linear coupling system that can give rise to the observed covariance matrix ¹⁴.

Similar to mutual information calculations conducted on other proteins ²⁴, ¹⁷, we find that backbone-backbone interactions computed from the covariance are weak compared to sidechain-sidechain interactions (compare red and blue boxes in Figure 1B and the mutual information in Figure S4). In these networks, some dihedral angles have a banding pattern, suggesting long-range interactions with many other dihedral angles (Figure 1B and Figure S1). In contrast, the inverse covariance matrix has localized and specific interactions. In addition, the stronger backbone-backbone interactions have a repeating pattern that resembles the contact map of the protein, and this pattern appears to repeat more weakly in backbone-sidechain interactions.

The banding pattern in the covariance matrix is also widespread in mutual information matrices, where they have been interpreted as long-range interactions important in protein allostery ²⁴, ¹⁷. The large number of long-distance edges produce ”hairball” networks, which led to the use of pruning algorithms ⁴⁵, ⁴⁶ or distance filters ²², ³⁶, ⁴⁷ in prior studies, in order to make network analysis tractable. Thus, we wondered whether the long-range interactions are capturing a physical feature of the dynamics. To answer this question, we investigate the reproducibility of the covariance, correlation, and inverse covariance matrices extracted from different replicates of MD simulations.

In Figure 2A, we contrast the matrix for one replicate in the upper-diagonal with the second replicate in the lower-diagonal and quantify the similarity in Figure 2B. For replicate MD simulations, we used the same initial protein structure with randomized solvation and initial velocities. Both covariance matrices have banding patterns suggesting hotspots that interact with many residues across the protein. However, each replicate has its own banding pattern, with interaction strengths that are over ten times greater than those found in the other replicate, indicating high variability in the networks that one would construct from replicate simulations.

Since the banding pattern is associated with dihedral angles with high variance, we also consider the correlation matrix, which normalizes the covariance matrix by the variance of each dihedral angle. It is visually apparent that the banding in the covariance matrix is not simply due to high variance because there are still bands in the correlation matrix. While normalizing to the correlation matrix uncovers some interactions in a contact map pattern, they are weak compared to the banding pattern (see Figure S5 for a lower maximum value on the color map scale).

In contrast to the irreproducible results obtained with the covariance matrix, for the inverse of the covariance matrix, we find a pattern that is visually similar to the 12Å contact map (Figure 2C). The diagonally symmetric contact map pattern in blue indicates similarly strong interactions for two replicates (Figure 2A). After quantifying the robustness across three replicates using the Jaccard similarity index, we find that the inverse covariance has higher similarity (59-72% shared edges) than the covariance (8-10%) or the correlation (13-16%) for $\psi-\psi$ backbone interactions (Figure 2B). The $\chi_{1}-\chi_{1}$ similarity values for inverse covariance are lower, but still higher than for the other two methods. These data clearly demonstrate that networks inferred from inverse covariance analysis are more robust than networks constructed from correlational measures.

Prompted by the strong visual resemblance between the inverse covariance network and the 12Å contact map and the complete absence of this pattern in the covariance network, we wondered if the inverse covariance matrix could be used to identify specific physical interactions. To answer this question, we overlaid the strongest edges ($\geq$ 97^th percentile) on the 12Å contact map, highlighting the backbone-backbone edges as blue dots and the sidechain-sidechain edges as red crosses (Figure 3A and Figure S10).

To correct for the high variability across replicates we previously found for the covariance networks, we averaged networks across the three replicates. Despite the averaging, the covariance network shows strong interactions across distant protein regions and are dominated by $\chi_{1}-\chi_{1}$ interactions.

To understand how these contrasting patterns affect interpretation, we next visualize strong interactions as edges drawn on the protein structure (Figure 3B). Since drawing all edges would make the covariance network indecipherable, we only show the edges originating from Lys4. In contrast, the inverse covariance network uncovers edges that mostly connect physically close residues. Specifically, we find a $\chi_{1}-\chi_{1}$ edge between Cys3-Cys44 for the sole disulfide bond in FimH_L, and that this edge is missing from the covariance network (red vs green arrow in Figure 3A).

Examining the contact map pattern of the inverse covariance network in more detail, we compare edge weight with the distance between C$\alpha$ atoms (Figure S9). We find that for backbone-backbone interactions, the strongest interactions are between residues connected by a peptide bond, followed by hydrogen bonds within $\beta$ sheets, and then non-bonding interactions.

After examining backbone-backbone interactions, we next looked at the progressively weaker interactions involving sidechains distal from the backbone (Figure 1B and Figure S10). We find the contact map pattern is still apparent for $\phi-\chi_{1}$ or $\psi-\chi_{1}$ interactions, but becomes very weak for $\phi-\chi_{2}$ or $\psi-\chi_{2}$ interactions, and becomes indistinguishable from noise for interactions between proximal and distal sidechain dihedrals. The inverse covariance analysis thus suggests that backbone dihedral motion is most strongly coupled to nearby backbone dihedrals and has more dissipated effects on sidechain dihedrals. This relationship is consistent with how backbone motions can sterically trap or free sidechains, whereas sidechain motions are more limited in their impact on backbone motion ⁴⁸.

The different strengths of interactions for backbone-backbone and backbone-sidechain edges suggests that qualitatively different types of interactions have different properties. For two residues i and j, the backbone-backbone edge $\phi-\psi[i,j]$ are larger than the backbone-sidechain edge $\chi_{1}-\chi_{1}[i,j]$, which is consistent with the physical differences between these two edge types. Moreover, most backbone-backbone edges within the same residue, $\phi-\psi[i,i]$, are stronger than backbone edges connecting to other residues, $\phi-\phi$[i, j] and $\psi-\psi$[i, j] (Figure S9).

Like FimH_L, the human immune cell adhesion protein, Siglec-8, is also a lectin with an immunoglobulin-like fold with a single disulfide bond. For Siglec-8, we compare the apo (no ligand) and holo (bound to the native 6’S-sLe_x ligand) states. Due to the rigid binding pocket loop that only differs by a few sidechain rearrangements, apo and holo Siglec-8 have extremely similar structures²⁹. Intriguingly, the rigidity of the CC’ binding pocket loop in apo Siglec-8 occurs in the absence of stabilizing secondary structure motifs ²⁹ and plays a major role in recognizing specific ligands in the airway to avoid autoimmunity ⁴³.

One hypothesized mechanism for stabilizing the CC’ loop in apo Siglec-8 is that the Arg70 sidechain forms hydrogen bonds with the loop backbone at Pro57 and Asp60 (Figure 5A)²⁹. We identified Arg70-Pro57 or Arg70-Asp60 hydrogen bonds with occupancy between 10-33% in seven out of 20 simulations from the ensemble of NMR structures, and higher occupancy (80 and 83%) in only two simulations. This suggests the hydrogen bond is rarely present. Nonetheless, the unstructured CC’ loop has surprisingly small fluctuations across the ensemble of MD simulations (Figure 5B-E).

Using our inferred interaction networks, we identified strong interactions within the CC’ loop, as well as interactions from outside the loop to its hinges at Ala53 and Pro62 (Figure S12). Consistent with the hydrogen bonding analysis, we only identified strong interactions with Arg70 in a few MD simulations, which disappeared after averaging the networks across the ensemble. Counter-intuitively, in the two structures where Arg70 does form stable hydrogen bonds with Pro57 and Asp60, the CC’ loop has paradoxically larger fluctuations, especially at Tyr58 and Gln59 at the loop tip (Figure S12). Taken together, we suspect that steric hindrance plays a role in stabilizing the CC’ loop internally and externally at the hinge edges.

Of note, this example demonstrates an instance where the inferred network approach provides more clarity than contact or interaction energy networks, because the inverse of the covariance matrix identifies the degrees of freedom that most strongly affect the fluctuations at the CC’ loop.

After examining why the CC’ loop is similar in apo and holo Siglec-8, we next compared inferred networks for the apo and holo states. The holo state has stronger $\chi_{1}$-$\chi_{1}$ interactions corresponding to the Cys31-Cys91 disulfide bond (Figure 5C and Figure S13). The disulfide bond is conserved in the siglec family and is located on the sheet of the $\beta$-sandwich opposite the binding pocket ²⁹. Nearby, we also observe other changes in interaction strength involving Asp90. Although distant from the binding site, the Asp90-Cys91-Ser92 motif in Siglec-8 is a variant of the Asn-Cys-Ser or -Thr motif that is conserved in the rest of the siglec family ⁵¹. Differences in interaction strength between apo and holo states identify changes in the dynamics of this evolutionarily conserved region during ligand-binding.

In contrast, reducing the Cys31-Cys91 disulfide bond in the holo state in silico has a different pattern (Figure S13). Both ligand-binding and the presence of the disulfide bond correspond to stronger a $\chi_{1}$-$\chi_{1}$ interaction at Cys31-Cys91. Both conditions also stabilize Cys31, indicated by decreased backbone dihedral fluctuations and increased duration within an extended secondary structure as assigned by the DSSP algorithm, and shorter Cys31-Cys91 C$\alpha$ distance. Taken together, we find that the network rearrangement that occurs with ligand-binding increases Siglec-8 conformational stability near an evolutionarily conserved disulfide bond, even though the region is not near the binding site.

Next, we investigated whether there are differences in the networks inferred from the ‘up’, ‘down’, and ‘off’ states of the RBD-SD1 domains of the SARS-CoV-2 spike protein (Figure 1A). The RBD connects to SD1 (Fig 5E) via two hinge-like loops that are more flexible in the up state than in the down or off states ⁵². While the up state has a different orientation around the hinge than the down and off states, they all have similar SD1 structures (C$\alpha$-RMSD $\leq$ 0.64Å) and somewhat similar RBD structures (C$\alpha$-RMSD $\leq$ 1.55Å).

Opening the hinge angle in the down-to-up transition is thought to make the binding site on the RBD available to attach to human ACE2 ³⁰, ⁴⁴. To focus on the RBD-SD1 hinge, we isolated these domains from the rest of the spike protein and ignored glycosylated sugars. While these simplifications limit the strength of our conclusions into the function of the spike protein, the RBD-SD1 structures nonetheless provide a useful system for comparing dynamics in a system that initially resembles rigid-body motion around a hinge.

Since the trimeric spike protein with one exposed RBD has two hidden RBDs, we isolated one RBD-SD1 fragment in the ‘up’ state, and two fragments in the ‘down’ state. We obtained another RBD-SD1 domain in the ‘off’ state from a structure with all RBDs hidden and 3-fold rotational symmetry. We first compared the RBD-SD1 fragments in the down and off conformations. Despite the structural similarity, we nevertheless detect differences in the inferred networks (Fig 5E) based on dynamics. We find the networks for the two down conformations are more similar to each other than to the networks for the off conformation (Figure S15).

We next compared RBD-SD1 fragments in the down and off conformations to the fragment in the up conformation (Figure S15). Surprisingly, we find that the differences in inferred networks among the down and up conformations are more localized, while the networks for the off conformation are distinct from the others. The localized differences are in the hook of the RBD that is exposed in the up conformation, and at an $\alpha$-helix near the RBD-SD1 hinge. Intriguingly, we also identified differences in the structural orientation and the inferred network interaction at the Cys336-Cys361 disulfide bond. In this region, one down conformation is more similar to the up conformation, while the other down conformation resembles the off conformation (Fig 5E).

As a proof-of-concept, we investigated whether it is possible to infer an interaction network for the trimeric SARS-CoV-2 spike protein, even though it is much larger than the RBD-SD1 domains. However, the amount of data required for network inference using the inverse covariance method scales faster than the number of residues. As a result, our approach required more data than was available from the two sets of publicly-accessible simulations for which the up and down states are labeled ⁵³, ⁵⁴. Instead, we chose the longer simulations that explore the transition among the down, up, and open states published by the Bowman lab ⁵⁵. Using 100,000 snapshots, we were able to infer an interaction network for the 3,363 residues of the spike protein (Figure S14). The number of snapshots is an order of magnitude larger than those currently available for labeled states. Thus, it may be feasible to compare states for the entire S1/S2 complex of the spike protein if there is sufficient data, or by choosing a less data-intensive method for solving the inverse problem.

Using dihedral angles enables us to easily account for sidechain dynamics avoid the artifacts introduced by structural alignment ³⁵. While protein dihedral angles have steep energy wells, we make the simplifying assumption that each dihedral angle has a harmonic potential energy function centered around one preferred orientation. This is only reasonable for short, picosecond-timescales, which is usually faster than the timescale for backbone dihedral flipping ⁵⁹. Moreover, let’s assume that dihedral angle pairs only have linear coupling, which neglects nonlinear interactions and many-body interactions^{27, 37, 38}.

We can describe these linearly coupled dihedral angles as the Hessian matrix, H, where $H_{ij}$ is the coupling constant between a pair of dihedrals, $i,j$. A positive coupling indicates rotation in the same direction. A negative coupling indicates rotation in opposite directions.

The harmonic coupling assumption means that $H_{ij}$ goes as $\frac{energy}{radians^{2}}$. Radians are a unitless measure, so $H_{ij}$ has units $J=Nm$, but we find it easier to use $J/rad^{2}$. Some combination of angular displacements from equilibrium, $\theta$, then produces the torque vector $f=-\textbf{H}\theta$. The potential energy difference from that at equilibrium is $\Delta U=\int-fd\theta=\frac{1}{2}\theta^{T}\textbf{H}\theta$.

We can also the use the Boltzmann relationship to describe the probability of a particular configuration in terms of angular displacements, $\theta$, as $p[\theta]\sim exp[-\Delta U/k_{b}T]$. The exponent is $-\theta^{T}\textbf{H}\theta/(2k_{b}T)$ and can be rearranged to $-\frac{1}{2}\theta^{T}(\frac{1}{k_{b}T}\textbf{H})\theta$. This yields the probability distribution function in terms of H (equation 1).

By recognizing the form of a multivariate Gaussian distribution with covariance matrix C (equation 2), we can see the inverse relationship between the Hessian and covariance matrices (3). This derivation mirrors the one for anisotropic elastic network models ^{25, 26}, which use displacements in Cartesian space, instead of an internal coordinate system of dihedral angles.

To demonstrate that the different patterns found between the two matrices apply to multiple MD trajectories, we show the covariance matrix diagonally across from its inverse for FimH (Supplementary Figure 1), Siglec-8 (Supplementary Figure 2), and the SARS-CoV-2 spike protein RBD-SD1 domains (Supplementary Figure 3). For each, we set the color scale maximum to the 97^th percentile. We show these in the same format as in Figure 1.

We highlight the backbone-backbone $\psi$-$\psi$ interaction in blue, and the sidechain-sidechain $\chi_{1}$-$\chi_{1}$ interaction in red. The covariance and mutual information (Supplementary Figure 4) matrices have stronger $\chi_{1}$-$\chi_{1}$ interactions than $\psi$-$\psi$. This is reversed for inverse covariance matrices.

While we do not use thresholds for comparing inferred interactions, we do use thresholds to visualize networks as the adjacency matrix and on the protein. In the matrix visualizations, we set the color scale maximum to the 97^th percentile. However, to better illustrate the pattern at lower values in the correlation matrix in Figure 2A, we also show a lower color scale maximum set to the 95^th in Figure S5. In contrast, for the inverse covariance matrix, showing lower values does not have a profound effect, since the contact map pattern is formed by inferred interactions with high values.

We also show how thresholds affect the Jaccard similarity. In Figure 2B, we used a threshold of 97^th percentile to create a mask where every edge above the threshold is set to 1, and edges below the threshold are set to 0. Here we show masks for thresholds at the 50, 95, 97, and 99^th percentiles, for the covariance, mutual information, and inverse covariance matrices (Figure S8). In Figure S7, we show that for strong interactions, the inverse covariance matrix has high similarity across replicates, while the other matrices show decreasing similarities. While other matrices have higher similarity around the 50^th percentile, the masks in Figure S8 show the extremely large number of edges included at these low thresholds.

The banded pattern indicates edges connecting one dihedral angle to many others, and an excessive number of these edges produces hairball networks. Thus, for the covariance, correlation, and mutual information matrices, there is a tradeoff where low thresholds with higher similarity produce hairball networks. In contrast, for the inverse covariance network, the high thresholds that produce higher similarity select for stronger interactions, producing physically interpretable networks that resemble the contact map.