Protein Fluctuations in Response to Random External Forces

Abstract

Elastic network models (ENMs) have been widely used in the last decades to investigate protein motions and dynamics. There the intrinsic fluctuations based on the isolated structures are obtained from the normal modes of these elastic networks, and they generally show good agreement with the B-factors extracted from X-ray crystallographic experiments, which are commonly considered to be indicators of protein flexibility. In this paper, we propose a new approach to analyze protein fluctuations and flexibility, which has a more appropriate physical basis. It is based on the application of random forces to the protein ENM to simulate the effects of collisions of solvent on a protein structure. For this purpose, we consider both the C^α-atom coarse-grained anisotropic network model (ANM) and an elastic network augmented with points included for the crystallized waters. We apply random forces to these protein networks everywhere, as well as only on the protein surface alone. Despite the randomness of the directions of the applied perturbations, the computed average displacements of the protein network show a remarkably good agreement with the experimental B-factors. In particular, for our set of 919 protein structures, we find that the highest correlation with the B-factors is obtained when applying forces to the external surface of the water-augmented ANM (an overall gain of 3% in the Pearson’s coefficient for the entire dataset, with improvements up to 30% for individual proteins), rather than when evaluating the fluctuations obtained from the normal modes of a standard C^α-atom coarse-grained ANM. It follows that protein fluctuations should be considered not just as the intrinsic fluctuations of the internal dynamics, but also equally well as responses to external solvent forces, or as a combination of both.

1. Introduction

The B-factors of a protein, the Debye-Waller factors or temperature factors, are measures of the atomic displacements about their equilibrium position, i.e., atomic fluctuations [1,2,3], but also the effects of multiple conformations as well as errors in the structures. They are generally accepted to be mostly the result of internal protein dynamics and any static disorder [4]. They have also been shown to be associated with protein flexibility and to correspond closely to protein mechanisms [5,6,7,8,9]. B-factors have been associated with protein flexibility, which is strictly related to protein action and function [10,11,12]. The experimental B-factors obtained from X-ray crystallography have been reproduced fairly accurately by various computational models.

One of the most widely used computational methods for investigating protein dynamics and fluctuations has been molecular dynamics (MD). MD simulations have proven their usefulness for investigations of protein folding, enzyme catalysis, and protein mechanisms in general [13,14,15]. Also, it has been shown that the MD-derived atomic fluctuations due to the internal protein motions show some degree of agreement with the experimental B-factors [16,17]. However, due to the high computational burden of MD simulations, these can sometimes be expensive for investigating the large molecular complexes, especially regarding the slowest protein motions, accessible only at long simulation times. These slow motions are in fact usually the ones most closely related to the functional mechanisms of the protein and can take place on longer time scales than may be accessible in standard MD simulations. The harmonicity assumption has been exploited for the extraction of the low-frequency protein dynamics [18,19,20,21]. Normal mode analysis (NMA) came into play as a simplified yet powerful tool to investigate the slower protein motions and for evaluating protein fluctuations and mechanisms [22,23], even in torsional space [24].

The seminal work of Tirion [25] showed that even a single-parameter harmonic potential, only based on the elastic properties of a network of Hookean springs connecting the protein atoms, was sufficient to reproduce the slow dynamics in good detail. All of the elastic network models are essentially entropic models since there is not usually any distinction of atom or amino acid types, i.e., all springs are taken to be similar in character. A further step towards simplification came with the coarse-graining development for these elastic network models (ENMs). Among the ENMs, the gaussian network model (GNM) was developed to obtain insights into protein dynamics and fluctuations simply by diagonalizing the Kirchhoff matrix, built by using the network contacts between close neighboring C^α atoms [26,27,28,29,30,31]. Despite the remarkable correlations obtained between the GNM-based fluctuations and the experimental B-factors, the GNM lacks the information about the directions of motions, since it assumes that the motions are fundamentally isotropic in all directions [28]. The anisotropic network model (ANM) was then developed to include the three-dimensional directionality in the calculation of protein motions [32]. The ANM was then improved by various research groups to achieve higher correlations between the computed fluctuations and the experimental B-factors [33,34,35,36,37]. These elastic models were subsequently used to study the conformational changes of proteins arising from sets of low-frequency modes [38,39,40,41,42,43,44,45,46,47,48] as well as to generate feasible pathways between two known conformations [48,49,50,51,52,53,54].

Structural elastic models, particularly the ANM, were applied widely for the investigation of protein dynamics, fluctuations, and mechanism. However, they are also well-suited for the analysis of the protein structural responses from the application of external perturbations. Based on the work from Ikeguchi et al. [55], who showed that protein conformational changes upon ligand-binding could be analyzed based on linear response theory, the perturbation-response scanning (PRS) method was proposed by the Atilgan group [56,57]. Randomly oriented forces were applied at selected residues, and the corresponding response of the ferric binding protein [56] and another 24 proteins [57] were found to agree fairly accurately with the experimentally detected conformational change. A similar study was conducted by Gerek and Ozkan [58] to study the allosteric network in PDZ domains. A PRS-based technique, coupled with energy-based Metropolis Monte Carlo (MMC) simulations, was carried out by Liu et al. [59] to simulate the closed-to-open conformational change of a GroEL subunit due to directional forces presumed to originate from exothermic ATP hydrolysis. Interestingly, some of the apparent conformational changes being attributed to the binding of ATP or ADP may originate from the exothermic forces generated by hydrolysis. Recently, it was also shown that the application of forces in a dynamic fashion is able to drive the conformational change with a strong directionality correlation [60]. Eyal and Bahar [61] investigated the mechanical response of protein structure to external pulling forces in order to detect the anisotropic mechanical resistance to explain the outcomes of single-molecule manipulation techniques. More recently, we made use of a similar pairwise force application methodology in order to measure the overall protein flexibility by using the engineering concepts of structural compliance and stiffness [62].

Most of the works based on the coarse-grained ENMs include only one or a few representative atoms of the amino acids in the protein network, e.g., the C^α atoms. Remarkably it has been seen that this geometric coarse-graining at the level of one point per amino acid yields almost exactly the same motions as from a full atomic elastic model. This result is believed to be the result of the dense packing leading to the strong stability of protein domains [63]. However, most of these models do not explicitly account for the protein surface, which is the part most exposed to the surrounding environment. Water and small molecules can often be tightly bound to the protein surface, thereby affecting what is actually considered to be the surface of a protein structure. The role of such tightly bound crystalized waters in protein dynamics has been studied in the last few decades [64,65,66,67]. There have also been investigations of the solvent network surrounding the protein and its effect on the dynamics [68,69,70]. The inclusion of water molecules in the structure yields some increases in the quality of calculated enthalpies and of the residue interaction network [71,72]. This is one of the important reasons why all-atom MD simulations usually include these explicit waters.

This paper presents a novel method based on random perturbations applied to protein ENMs to assess protein fluctuations and flexibilities. Random forces are applied both throughout the complete protein elastic network and also separately to only the protein surface, which is exposed to the surrounding environment. In addition, a water-enriched ENM is considered, where the water molecules whose coordinates are given in the Protein Data Bank (PDB) files [73] are used as additional nodes in the elastic network. These latter force application simulations aim to mimic the random collisions occurring on the protein structure due to the interaction with the solvent and other solutes. From the calculation of the displacements within the protein network, i.e., the protein responses, we show that a good correlation is found with the experimental B-factors, thus leading to a good prediction of the protein flexibilities. The correlations with the usual mode-based fluctuations are also reported for comparison. It is also found that, in most cases, applying random forces on the surface of the water-enriched protein ENM leads to the highest correlation between the resulting displacements and the experimental B-factors. This demonstrates that the protein fluctuations may reflect more than the internal dynamics alone, and also include some effects from the continuous random bombardments or restraints by the surrounding solvent on the protein structure.

2. Methods

In this section, we briefly recount the fundamentals of the Anisotropic Network Model (ANM) [32], that is commonly used for generating the fluctuations in terms of the normal modes, and then we describe the computational framework related to the presently adopted force applications on the elastic networks.

2.1. Anisotropic Network Model (ANM) and the Calculation of Normal Mode-Based Fluctuations

ANM relies on the assumption that proteins can be modeled as simple elastic networks, made up of point nodes connected by linear elastic springs, allowing insights regarding fluctuations and global mechanisms [32,34,74]. For a system of N points, e.g., N residues in the one-bead-per-residue coarse-grained representation, the 3N × 3N Hessian matrix of the system takes the following form:

$$\mathit{H}=\mathbf{\left[}\begin{array}{ccccccc}{\mathit{H}}_{\mathbf{1},\mathbf{1}}& \dots & {\mathit{H}}_{\mathbf{1},\mathit{i}}& \dots & {\mathit{H}}_{\mathbf{1},\mathit{j}}& \dots & {\mathit{H}}_{\mathbf{1},\mathit{N}}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ {\mathit{H}}_{\mathit{i},\mathbf{1}}& \dots & {\mathit{H}}_{\mathit{i},\mathit{i}}& \dots & {\mathit{H}}_{\mathit{i},\mathit{j}}& \dots & {\mathit{H}}_{\mathit{i},\mathit{N}}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ {\mathit{H}}_{\mathit{j},\mathbf{1}}& \dots & {\mathit{H}}_{\mathit{j},\mathit{i}}& \dots & {\mathit{H}}_{\mathit{j},\mathit{j}}& \dots & {\mathit{H}}_{\mathit{j},\mathit{N}}\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ {\mathit{H}}_{\mathit{N},\mathbf{1}}& \dots & {\mathit{H}}_{\mathit{N},\mathit{i}}& \dots & {\mathit{H}}_{\mathit{N},\mathit{j}}& \dots & {\mathit{H}}_{\mathit{N},\mathit{N}}\end{array}\mathbf{\right]},$$

(1)

where each 3 × 3 submatrix H_i,j contains the stiffness coefficients of the springs connecting nodes i and j. The off-diagonal submatrix H_i,j is computed based on the harmonic potential of the elastic spring with force constant γ_i,j. The diagonal submatrices H_i,i are calculated as the summation involving all the nodes linked to ith node with a negative sign [32]. The model, and consequently the Hessian matrix, depends on some numerical parameters: the usual model uses a cut-off limit in distance to define the network topology. The original ANM was developed by considering equal spring constants for all connections, i.e., γ_i,j = γ, and a geometrical cut-off r_c was applied in order to consider springs placed only between close nodes. Typical values of r_c in the ANM are around 15 Å. Later on, distance-dependent force constants were introduced [33,34], as:

$${\gamma }_{i,j}\propto \frac{1}{{(r_{i,j}{}^0)}^p},$$

(2)

where p represents an inverse number for the decay parameter that allows connecting all points in a structure, with springs with variable strength—lower spring constants for longer inter-node distances r_i,j⁰. This distance-dependent spring network was shown to provide an improved agreement between the results and experimental data [33,34].

Once the Hessian matrix is computed based on the protein coordinates from the PDB file [73] and the spring connectivity, the 3N eigenvalues λ_n and 3N eigenvectors U_n are obtained by solving the eigenvalue-eigenvectors decomposition. Due to the fact that the protein structure is usually not externally constrained, the first six eigenvalues are found to be zero, with corresponding mode shapes associated with the six rigid-body motions of the whole molecule. Therefore, these six motions are factored out and singular value decomposition is used to obtain the normal modes. Hence, the fluctuations based on the normal modes can be easily calculated as [22,75,76]:

$$B_i=\frac{8}{3}{\pi }^2{k}_{B}T{{\displaystyle \sum }}_{n=7}^{3N}\frac{U_{i,n}{}^2}{{\lambda }_n},$$

(3)

where B_i represents the computed B-factor for residue i, k_B is the Boltzmann constant, T is the absolute temperature, U_i,n stands for the displacement of node i in the nth mode, and λ_n is the eigenvalue for the nth mode.

2.2. Force Application on Elastic Networks

Here, we propose a new approach for the calculation of protein flexibility and fluctuations. This approach is based on the application of random forces on protein elastic networks. These forces are intended to simulate the external perturbations that arise from the external environment, i.e., protein-solvent interactions, Brownian motions, collisions of molecules, etc. The reality, however, is that the environment is not usually known, but cryoEM has the promise of providing some information about this.

In this work, we use two different ENMs for modelling the protein structure. The first one is the parameter-free anisotropic network model (pfANM) [33], where the C^α atoms are the only nodes used to build the protein network and all the C^α-C^α connections are considered to be linked with distance-dependent springs. In the second model, the water molecules contained in the PDB file are also added to the network as additional nodes. Additional springs are correspondingly created that connect the water molecules to all the other nodes of the network. We refer to this second model as the water-pfANM (wpfANM). Both models are built by considering a decay exponent p for the spring constant equal to 3 (see Equation (2)), based on results shown to yield the best results analyzed in our previous work [62].

The response of the protein structure to external perturbations is evaluated by applying forces to the nodes of the network and consequently computing the corresponding elastic displacements. Various force application patterns are considered here. For the pfANM, the perturbations are applied both to the complete structure, i.e., on all C^α atoms, and separately only to the nodes lying on the external protein surface. For the wpfANM, three different force patterns are considered: (1) forces acting on the entire network, i.e., all the C^α atoms and water molecules, (2) only on the nodes lying on the protein-water network surface, and (3) only on the water molecules.

For each of the models considered (pfANM and wpfANM) and their force application patterns, the calculation is based on the generation of a random 3 × 1 force vector F_i^s for each node i to be perturbed for each simulation s. The three scalar components of this force vector, i.e., F_i,x^s, F_i,y^s and F_i,z^s, are sampled from a uniform distribution U in the interval [–1,1]. The complete force vector F^s is then generated by assembling all the 3 × 1 nodal vectors. Note that F^s is a 3N × 1 vector, with N the number of C^α atoms in the pfANM case or the total number of C^α atoms plus water molecules in the wpfANM. Once the force vector F^s is defined, it is straightforward to compute the 3N × 1 displacement vector δ^s for each simulation s that contains the elastic displacements of the nodes, i.e., the protein response, as follows:

$${\mathit{\delta }}^{\mathit{s}}={\mathit{H}}^{-\mathbf{1}}{\mathit{F}}^{\mathit{s}},$$

(4)

where H⁻¹ is the pseudo-inverse of the elastic network. H⁻¹ can be computed from the eigenvalues and eigenvectors of the Hessian matrix as:

$${\tilde{\mathit{H}}}^{-\mathbf{1}}={{\displaystyle \sum }}_{n=7}^{3N}\frac{{\mathit{U}}_{\mathit{n}}{\mathit{U}}_{\mathit{n}}{}^{\mathit{T}}}{{\lambda }_n}.$$

(5)

From the displacement vector δ^s, the total displacement of each ith node can be computed as:

$${\delta }_i{}^s=\sqrt{{\left({\delta }_{i,x}{}^s\right)}^2+{\left({\delta }_{i,y}{}^s\right)}^2+{\left({\delta }_{i,z}{}^s\right)}^2},$$

(6)

With δ_i,x^s, δ_i,y^s and δ_i,z^s being the three Cartesian components of the node displacements.

This procedure is repeated multiple times in order to generate a sample with different random force vectors F^s and evaluating the corresponding node displacements δ^s each time. The average displacement of each node i is then evaluated as the average of all the displacements δ_i^s over the total number of simulations S:

$${\delta }_i=\frac{1}{S}{{\displaystyle \sum }}_{s=1}^S{\delta }_i{}^s.$$

(7)

In this analysis, we have generated a sample of 10,000, i.e., S = 10,000. Then, the average displacement δ_i of the ith residue for the sample can be compared to the experimental B-factors available in the PDB file. Pearson’s correlation coefficient can finally be used to estimate the similarity between the two distributions, i.e., between the experimental B-factors and the simulated average displacements of the protein network due to the random perturbations. As a result, high Pearson’s coefficients would indicate a high degree of similarity between the computed protein fluctuations and the experimental B-factors.

We mention above that different force application patterns are considered in this study. Specifically, besides applying forces to the entire protein network in the pfANM and wpfANM, and to the water molecules alone in the wpfANM, we also apply forces only to the nodes lying on the external protein surface (pfANM) and on the protein-water network surface (wpfANM). The reason for this is due to the fact that the effect of random collisions is more likely to occur on the exterior protein surface, rather than in the interior. For this purpose, the surface residues were calculated by computing the boundary geometry of the set of 3D coordinates of the network points. The external nodes were defined as those lying on the boundary surface. In this analysis, the generation of this surface was dependent on a parameter, known as the shrink factor. The shrink factor characterizes the amount of shrinkage of the boundary geometry, with values ranging from 0 to 1: zero corresponds to the convex hull, one corresponds to the maximum shrunk boundary. Note that the shrink factors used here for the generation of the external surface correspond to the normalized alpha shape, recently used by us [77] to extract hinge-domain information from protein structures. By using the approach based on the shrink factor, different external surfaces are generated by varying this numerical parameter from 0 to 1, in steps of 0.1. The surfaces obtained are finally used to select the external nodes on which the external forces will be applied.

2.3. Protein Dataset and Summary of the Models and Analyses

The analysis was performed on the same dataset used in our previous work [62] that includes 921 high-resolution X-ray single-chain protein structures from the PDB. The resolution of the selected crystal structures is below 1.3 Å, with a maximum sequence identity of 30%. Two structures were removed from the dataset, i.e., 1IXH and 2BK9, because of errors found in the PDB file regarding the waters. The size of the 919 final proteins range from 101 to 1174 residues.

As mentioned in Section 2.2, the pfANM and wpfANM were built for all protein structures by neglecting and considering water molecules, respectively. For both models, after evaluating the Hessian matrix, the mode-based fluctuations were evaluated according to Equation (3) and compared to the experimental B-factors. Then, for the pfANM, random forces were applied: (1) to all the nodes of the network and (2) only to the external nodes. Similarly, for the wpfANM, forces are considered to be acting on: (1) the entire network, (2) only on the external nodes and (3) only on the waters. From each of these simulations, the nodal displacements are computed from Equation (7) and Pearson correlation coefficients with the experimental B-factors are also calculated. In

Table 1. The Designators for the Pearson Collection Coefficients for the Various Elastic Network Models.

Model	pfANM			wpfANM
Nodes in the network	Cα atoms			Cα atoms + water molecules *
Analysis	Mode-based fluctuations	Force application		Mode-based fluctuations	Force application
Nodes perturbed	-	All nodes	External nodes	-	All nodes	External nodes **	Water nodes ***
Correlation coefficient	ρFL	ρFR,ALL	ρFR,EXT	ρW,FL	ρW,FR,ALL	ρW,FR,EXT	ρW,FR,WAT

* All waters in the pdb are included, without checking whether they are exterior. ** All nodes lying on the external surface (boundary geometry) of the protein-water network are perturbed by random forces; *** Only water molecules in the protein-water network are perturbed by random forces.

, the summary of the models and simulations and their designators for their Pearson correlation coefficients are given.

3. Results and Discussion

3.1. Fluctuations and Force Application on the pfANM

In this section, the flexibility of the protein structure will be investigated by using traditional pfANM mode-based fluctuations as well as from the outcome of the random force applications to the protein. Figure 1 shows the correlation coefficients obtained from a comparison of the experimental B-factors with the mode-based fluctuations as well as the average displacements due to the applied random forces. Figure 1a,b report the Pearson coefficients obtained for the 919 single-chain protein structures, with the values ordered by ascending values of ρ_FL. Figure 1c shows the statistical distributions of the correlation coefficients, whose median values and standard deviations are reported in the keys.

Figure 1a,b shows the distributions of correlation coefficients ρ_FR,ALL and ρ_FR,EXT, that are due to the random force applications, and these are observed to oscillate near the population of ρ_FL. This means that the average displacements of the protein elastic network due to the force application are indeed well correlated with the experimental B-factors, with a similar agreement as for the traditionally used mode-based fluctuations. The same conclusions can also be drawn by looking at Figure 1c, where the three statistical distributions exhibit the same pattern and similar median values, i.e., 0.63, 0.63 and 0.64, for ρ_FL, ρ_FR,ALL and ρ_FR,EXT, respectively. Therefore, it cannot be concluded that applying random forces on the protein structures always enhances the correlation with the experimental B-factors, while it can be concluded that perturbing the protein structure by applying random forces leads to a good estimate of the experimentally determined fluctuations, at least as good as those found with the normal modes.

As mentioned in the previous section, the application of random forces on the external protein surface requires the selection of the nodes that lie on the exterior. For this purpose, various external boundaries were generated by changing the shrink factor of the surface, varying from 0 to 1 with steps of 0.1. Figure 2 shows an example of different external surfaces generated with shrink factors of 0, 0.2, 0.4, 0.6, 0.8 and 1 for the infrared fluorescent protein (PDB: 5AJG). As can be observed, increasing the shrink factor leads to considering a higher number of nodes lying on the surface, which in turn has a more detailed structure. Since the primary determinant of a structure’s dynamics is its shape, clearly the most detailed structure would be expected to be the best [78].

For each of these generated surfaces, the random forces were applied only to the nodes lying on the external boundary. The resulting network displacements were then evaluated from Equations (4)–(7). It follows that the correlation of the average displacements (from 10,000 samples) due to the application of the external forces, i.e., ρ_FR,EXT, also depends on the adopted shrink factor. The shrink factor that leads to the maximum value of ρ_FR,EXT for each protein is then selected as the optimal one. Figure 3 reports the statistical distribution of the optimal shrink factors obtained for the 919 single-chain proteins. As can be observed, the optimal shrink factor assumes almost all values, meaning that it is strongly protein-specific. Nevertheless, a slight preference towards shrink factors equal to 1.0 is observed.

In Figure 4 we show results for the example of the infrared fluorescent protein (PDB: 5AJG), where the correlation coefficients ρ_FL, ρ_FR,ALL and ρ_FR,EXT, are shown depending on the shrink factor for the external surface representation. From these calculations, we obtain a correlation between the B-factors and the traditional mode-based fluctuations ρ_FL equal to 0.56, a correlation with the displacements resulting from the application of forces to the entire structure ρ_FR,ALL equal to 0.62, and a correlation derived from perturbations only on the external surface ρ_FR,EXT which varies with the shrink factor and reaches a maximum value of 0.62 for an optimal shrink factor of 1.0. As can be seen from the results, in this case, applying random forces on the protein network enhances the correlation with the experimental B-factors of 6% (0.62 vs. 0.56) compared with the usual mode-based fluctuations. This result points out the cohesive nature of the protein structure, showing that the point of application of forces does not matter much, with the result of applying forces in all possible directions on the surface yielding nearly the same result as applying them in all directions throughout the structure when the surface representation is detailed enough.

Figures S1–S5 in the Supplementary Material report similar results, obtained by adopting different values of p for the decay exponent of the ENM spring constants (p = 1, 2, 4, 6 and 12). As can be seen there, similar conclusions can be drawn for these cases. Note that, for this protein, higher values of p, e.g., p = 4 and 6, lead to a greater enhancement in the correlation with experimental B-factors when the ENM forces are applied, rather than just looking at the intrinsic dynamical fluctuations.

3.2. Incorporating Waters into the Computations

In this section, we show results obtained by also including the localized waters given in the crystal structure as part of the structure for defining the elastic network. There is some ongoing debate about whether or not these bound waters should be considered as an actual part of the structure. Each high-quality protein structure available in the PDB contains a substantial number of waters that were present in the crystal formed at low temperatures. The open question is whether these remain bound at higher temperatures. These molecules often typically form a network of hydrogen bonds with the side chains of polar amino acids on the protein surface and thus can appear to be quite stable. It follows that these waters may possibly cause some changes to the overall flexibility and dynamics of any given protein. Moreover, since we are interested in looking at the responses of each protein structure due to external perturbations, the inclusion of these external waters would be expected to affect the motions to some extent. Figure 5 shows a surface representation of the infrared fluorescent protein (PDB: 5AJG), with and without the addition of water molecules in Figure 5a,b, respectively. The protein structure is shown in light-gray, with the surface depiction highlighting the external surface and cavities. Water molecules available in the PDB are shown in Figure 5a as blue spheres. As can be seen from the comparison between the two figures, most of the crystallized water molecules are bound in concave parts of the protein surface, and thus act to smooth the structure [79]. This smoothing might restrict the flexibility of certain parts of the proteins that might cause problems for the mechanisms otherwise; this can be looked at as flowing liquid that represents protein motion, if a certain direction of the flow is restricted, it may change the overall flow of the water, therefore, to have an optimal flow path (specific functional protein motion), the restrictions (waters) are as important as the structure itself.

As mentioned in the previous section, the wpfANM is built as the usual pfANM with p equal to 3 (see Equation (2)), except that both the coordinates of the C^α atoms and the water molecules are now considered as a part of the whole structure. Based on the resulting wpfANM, the mode-based fluctuations are evaluated from Equation (3), whereas the average displacements resulting from the random perturbations are computed according to Equations (4)–(7). In the latter cases, as explained above, three types of force application are considered, as shown in Table 1. From the calculations, we then obtain four correlation coefficients to compare with the experimental B-factors, namely ρ_W,FL, ρ_W,FR,ALL, ρ_W,FR,EXT and ρ_W,FR,WAT, with these being described in Table 1.

Figure 6 shows the correlation coefficients for the 919 proteins investigated. Figure 6a–c report the correlations ordered by increasing values of ρ_W,FL, whereas Figure 6d displays the statistical distribution for all four Pearson coefficients, whose median values and standard deviations are shown in the key. From the results reported in Figure 6, it follows that applying forces on the protein network slightly (<10%) enhances the prediction of the B-factors with respect to the traditional mode-based fluctuations. As a matter of fact, the median value of ρ_W,FL was found to be 0.57 for the selected dataset, whereas the median values of ρ_W,FR,ALL, ρ_W,FR,EXT and ρ_W,FR,WAT were 0.60, 0.66, and 0.59, respectively. It is evident that applying random forces on the surface of the network (which now considers also the layer of water molecules) yields a significant 10% gain in the correlation with the experimental fluctuations, compared to the mode-based fluctuations.

Also in the case of the wpfANM, the surface of the network is not unique but depends on the adopted shrink factor. The water molecules contained in the network play a major role in the definition of the surface since they are mostly placed on the exterior of the structure (see Figure 5). As an example, Figure 7 shows the different surfaces generated for the infrared fluorescent protein (PDB: 5AJG, with water molecules included) with shrink factors equal to 0, 0.2, 0.4, 0.6, 0.8 and 1. As can be seen by comparing Figure 7 to Figure 2, other than the selected value of the shrink factor, the shape of the external surface also depends on the presence of the water molecules within the network. Similarly to what was shown in the previous section, Figure 8 shows the optimal shrink factors (with the best correlation with B-factors) obtained for the 919 water-enriched protein structures. Again, the distribution of the optimal shrink factor is rather uniform, although in this case a slight bias towards the convex hulls surfaces appears, i.e., with a shrink factor equal to 0. This probably reflects a preference for smoother structures when water molecules are added (see Figure 5).

Figure 9 shows the correlation coefficients of the wpfANM of the infrared fluorescent protein (PDB: 5AJG), and how it depends on the adopted shrink factor for the surface. From the calculations, we obtain a correlation with the mode-based fluctuations ρ_W,FL equal to 0.84, a correlation with the displacements from the application of forces to the entire structure ρ_W,FR,ALL equal to 0.85, a correlation with the displacements from the application of forces only to the water molecules ρ_W,FR,WAT equal to 0.63, and a correlation with the displacements due to the perturbations only on the surface ρ_W,FR,EXT that varies with the shrink factor and reaches a maximum of 0.83 for the optimal shrink factor of 0.3. In this case, it is remarkable that by applying random perturbations on only 31% of the nodes (corresponding to a shrink factor of 0.3), we obtain a high correlation with the experimental data (0.83 Pearson coefficient). It should be noted that this correlation is found by comparing the experimental B-factors of all protein residues (both on the surface and within the core) to the computed average displacements due to the application of perturbations only on the external part of the elastic network. Thus, it follows that even perturbing a small portion of the protein surface (31%) allows us to predict fairly accurately the fluctuations of the entire protein. This result has its origin in the strong coupling throughout the elastic network model: since the ENM is a highly cooperative model, the perturbation of only a small part of the structure can indeed generate fluctuations over the entire protein. This arises from the specific features of the three-dimensional protein structure and all of the internal connections, which are both built into the ENM.

Figures S6–S10 in the Supplementary Material show similar results as those reported in Figure 9, but for different exponents p in the water-augmented ENM, namely, p = 1, 2, 4, 6 and 12. Despite some differences in the numerical values of the correlation coefficients, similar conclusions are reached as for those in Figure 9.

3.3. Comparison between pfANM and wpfANM Results

In the previous sections, it has been shown that perturbing the protein structure with random forces, either on the entire structure or on the surface, generally leads to a fairly accurate prediction of the protein fluctuations and flexibility. In a large number of cases, it has also been found that the agreement with the experimental B-factors was improved with respect to considering the traditional fluctuations of the usual elastic network. As an example, the results shown in the previous sections for the infrared fluorescent protein (PDB: 5AJG) are reported together in Figure 10 in terms of correlation coefficients with the experimental B-factors. Several observations can be made regarding this figure. First, it is evident that the traditionally employed internal protein fluctuations provide the lowest correlation with the experimental data, with a correlation of about 55% (ρ_FL). Conversely, applying random perturbations on the same elastic network leads to a 5% gain in the correlations with the B-factors (ρ_FR,ALL and ρ_FR,EXT). Furthermore, adding the PDB water molecules to the elastic network further improves the correlation with the experimental data. In this case, considering the internal protein fluctuations of the water-enriched elastic network or applying random forces to it yields correlation coefficients of about 85% (ρ_W,FL, ρ_W,FR,ALL and ρ_W,FR,EXT). Also, applying random forces only on these water molecules, i.e., not perturbing the protein molecule directly but only the water molecules attached to the network (see Figure 5a), leads to a correlation of about 60% (ρ_W,FR,WAT), which is still higher than the correlation obtained with the traditional mode-based internal protein fluctuations of the ENM (ρ_FL = 55%).

Figures S11–S15 in the Supplementary Material show similar outcomes, but are obtained by changing the exponent p of the ENM. As can be seen, similar conclusions are drawn. In all cases, it is found that perturbing the protein with random forces provides an enhancement of the correlation with experimental B-factors, rather than considering the traditional ENM with only intra-protein interactions. Moreover, the inclusion of waters in the ENM leads to further improvements in the correlation, with some gains in the Pearson coefficient being as high as 30–35%.

The example shown in Figure 10 obviously refers to one single case, but these results were found for a decent amount of protein structures. For other protein structures, the addition of water molecules to the protein network led to results which were quite similar to the ones obtained in the classical way, i.e., calculating the mode-based fluctuations of the traditional no-water elastic network.

Figure 11a reports the median values and standard deviations of the seven correlation coefficients obtained for the dataset of 919 single-chain protein structures, as reported in the keys of Figure 1c and Figure 6d. As can be observed, the median values lie in the range 0.60–0.65, and present a similar standard deviation (15–20%). However, the distribution with the highest median value and the lowest standard deviation was found for the analysis involving the application of random perturbation on the external surface of the wpfANM, i.e., ρ_W,FR,EXT. A direct comparison between the statistical distribution of ρ_W,FR,EXT and the one related to the traditional mode-based internal protein fluctuations of the no-water pfANM, i.e., ρ_W,FL, is reported in Figure 11b. From the direct comparison of the two distributions, it is clear that applying random forces on the surface of the wpfANM leads to an overall, yet slight, improvement of the correlation with the experimental data.

Based on the correlation coefficients obtained for our entire dataset of protein structures, the analysis among the seven performed ones (Table 1) that gave the highest correlations with the experimental data was noted. Figure 12a shows the relative number of such occurrences for each type of analysis. It was found that, out of the 919 investigated protein structures, the mode-based fluctuations of the pfANM provided the highest correlation coefficient in 107 cases (11.6%), the application of random forces on the entire pfANM in 74 cases (8.1%), the application of forces on the external surface of the pfANM for 172 cases (18.7%), the mode-based fluctuations of the wpfANM in 60 cases (6.5%), the application of forces on the entire wpfANM in 70 cases (7.62%), the application of forces only on the water molecules of the wpfANM in 93 cases (10.1%), while the application of forces on the external surface of the wpfANM were provided in 343 cases (37.3%). As can be seen and has already been discussed concerning Figure 11, the application of random forces on the surface of the water-enriched protein network is statistically the best performing, although not the only one, with regards to better predicting protein fluctuations and flexibility in terms of experimental B-factors. Moreover, looking cumulatively at the analyses FR,EXT and W,FR,EXT, in 515 cases (56.0%), the application of random forces on the external surface of the protein network yields the highest correlation coefficients, thus confirming that perturbing the external protein surface can induce a response in good agreement with the experimental B-factors.

Figure 12b shows the occurrence of the highest correlation coefficients for the two models, i.e., the pfANM vs. wpfANM. From the outcomes, it was obtained that in 353 cases (38.41%) the highest correlation with the B-factors was obtained with the pfANM, whereas in the remaining 566 cases (61.59%) the wpfANM was allowed to reach the highest correlation with the experimental data. This confirms that considering the PDB water molecules might actually enhance the prediction of the protein fluctuations and therefore the numerical evaluation of the experimental B-factors.

4. Conclusions

Research carried out in the last decades has shown that protein fluctuations and flexibility, as measured by the experimental B-factors, mainly arise from the internal protein motions and inherent dynamics. The dynamics are known to originate from the tertiary structure, as recognized within the fundamental sequence-structure-dynamics-function paradigm of protein action. Therefore, it should be more appropriate to say that protein fluctuations and flexibility arise from the protein structure and can be mediated by its dynamics. As a matter of fact, in a recent work [62], we have shown that the overall protein flexibility, as measured by the experimental B-factors, can also be retrieved by applying pairwise static forces to the protein ENM and measuring the amount of compliance against these external pulling forces.

In this paper we have proposed an additional viewpoint as regards protein fluctuations and flexibility. We applied random static forces throughout the protein elastic network and evaluated the response of the network via the computation of average nodal displacements. From the comparison of these average displacements against the experimental B-factors, we found that the protein flexibility, and therefore its fluctuations, can indeed be elucidated with such a procedure. Also, we found that if these perturbations are applied on the protein surface, and if crystallized water molecules are also inserted into the model, higher correlations with experimental data can often be found. This suggests that protein fluctuations can also be seen as (fully or partly) the response of the protein structure to external forces, which might be induced by the continuous collisions of the solvent and other solute molecules around the protein structure.

It is important to mention that the goal of the analysis presented here was to propose a new methodology, and therefore a new perspective, to understand protein fluctuations. However, no additional work has been carried out yet as regards the optimization of the elastic model upon which the random perturbations are applied. This might eventually improve the correlation with the experimental data and is our next goal. As described in the text, the pfANM has been used for the standard ENM, whereas a wpfANM has been generated in order to account for the presence of crystallized water, where water molecules have simply been added as additional nodes to the network. We plan to optimize the elastic model by including different spring constants for the C^α-C^α connections, C^α-water connections and water-water connections, which should simulate more realistically the different atomic interactions (residue-residue, residue-water, water-water).

Attention has been paid to the external surface of the network. Thus, we are also planning to use a different version of ENM, where the contact topology is not generated by using the traditional cut-off limit, but using alpha-shapes associated with Delaunay tessellations. A recent work from Koehl et al. [36] showed that such a procedure is able to generate elastic models with enhanced agreement with experimental data. Applying external forces on such optimized models, and adding water molecules as well, might enhance the correlation with experimental B-factors, allowing for the better explanation of fluctuations, and therefore the way a protein moves and functions.