## 1 Introduction
[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 [PERSON] 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 [PERSON] and [PERSON] [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 [PERSON] 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]. [PERSON] and [PERSON] [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\({}^{\alpha}\) 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. Thisresult 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.