provides atomic-detailed evidence of the ability of the selected mutations to modulate the dynamic properties of AHA, and offers unprecedented insight in the way this modulation takes place in a large protein. Interestingly, a correlation (with a Pearson corre- lation coefficient higher than 0.78) exists between the combination of mesophilic-like salt bridge interactions acquired by the AHA mutants (highlighted in Table S1 in Supporting Information S1) and the experimentally determined kcat/Km and Tm values of each mutant [37]. These interactions include, in particular, the restored PPA-like D150A-K190A salt-bridge, the loss of K177A- D203A salt-bridge and the modification of interactions involving K334A. In a broader context, the approach applied here, integrating different properties related to protein dynamics (in particular the parallel and integrated evaluation of correlated motions and electrostatic interaction network), could be useful to predict and rationalize the effects induced by a mutation on protein activity and stability of an extremophilic enzymes, if comparative studies with other temperature-adapted counterparts can be carried out. Materials and Methods System setup and molecular dynamics simulations The X-ray structure of AHA (PDB entry 1AQH [33]) was in silico modified to model 7 mutant forms: 5 single mutants AHANR (N12R), AHAQI (Q164I), AHAVF (V196F) and AHATV (T232V) [35]; one double mutant AHASS (Q58C A99C) and 2 multiple mutants AHA5 (N150D, Q164I, V196F, T232V, K300R) and AHA5SS (Q58C, A99C, N150D, Q164I, V196F, T232V, K300R) [37]. MD simulations were performed using GROMOS96 force- field with the GROMACS software (www.gromacs.org). The mutated forms of AHA were soaked in a dodecahedral box including counterions (Na+) and SPC [65] water molecules using periodic boundary conditions, with a minimum distance between the solute and the box of 0.5 nm (average box size of 527.66 nm3 including 49893 atoms) and the systems were equilibrated in several steps [50]. Productive MD simulations were performed in the NPT ensemble (300K, 1 bar and 2 fs time-step with LINCS algorithm [66]). Electrostatic interactions were calculated using PME [67] summation scheme. Van der Waals and Coulomb interactions were truncated at 0.8 nm and conformations stored every 2 ps. To improve conformational sampling, six independent 6 ns simulations were carried out for each system. MD simulations showing the slowest convergence of root mean square deviation (rmsd) values have been elongated up to 10 ns (run 2 of AHAQI and run 3 of AHA5SS) and the first 0.2–2 ns depending on the system have been discarded to assure stability of the trajectories. In fact, MD multiple- replica approach allows a wider conformational sampling than only few longer molecular dynamics simulations, especially in the case of large multi-domain proteins as AHA is, if the stability of the trajectory and convergence of the analyzed properties have been carefully verified [68–71]. The equilibrated portions of the simulations of the same system were joined into a macro-trajectory. The conformational ensemble collected for each mutant is compa- rable, in terms of simulations protocol, degree of essential motions captured by the first principal components and number of frames in the collected ensemble to that achieved for the wild-type cold- (AHA) and warm-adapted (PPA) adapted amylases previously published [50], in order to carry out a prompt comparison of dynamical properties. Molecular and structural properties Secondary structure was determined for all stored conforma- tions using DSSP [72] and the most frequently attained secondary Molecular Dynamics of Psychrophilic Enzyme Mutants structure for each residue was evaluated to obtain a residue- dependent persistence degree of secondary structure profile (PDSSP). Intramolecular hydrogen bonds (H-bonds) were evalu- ated using both a donor-acceptor distance cut-off of 0.35 nm and an acceptor-donor-hydrogen angular cut-off (#30u); ensemble averages were normalized using the total number of donors and acceptors. Surrounding hydrophobicity (Hp) was computed for each residue according to [73]. In particular Hp was computed for each residue as the sum of hydrophobicity indexes, obtained from thermodynamic transfer experiments [74], of the surrounding residues within a 0.8 nm distance cutoff: N Hp(i)~ X j=i H(rd {rij )hj where the summation runs over the whole protein, rij is the distance between the Ca atoms of residues i and j, H(x) is the Heaviside step function [H(x) = 1 if x$0, and zero otherwise], rd is the distance cutoff and hj is the hydrophobicity index of residues j in kcal/mol. The adopted cutoff has been demonstrated to be suitable to characterize the hydrophobic behaviour of amino acids and to accommodate both the local and non local interactions [73,75]. Hp is therefore a measure of the tendency of the environment of each residue to exclude water, computed on the basis of the 3D structure using experimental information. Computed on equilib- rium ensembles obtained from MD simulations, it gives an indication of the relative strength and stability of the hydrophobic clusters in a protein. It should be in principle able of identifying the effects of a mutation taking into account both the short-range effects, due to changes in the physico-chemical properties of the mutated residues, and the long-range effects of mutation on structure and dynamics. However, it is worth mentioning that the accuracy of the computed values is hindered, independently of the hydrophobicity scale adopted, by two strong assumptions. In fact, by summing up the contributions of different residues, the Hp calculation assumes that these contributions are independent and therefore additive, which may lead to misleading conclusions. On the other hand, the use of hydrophobicity indexes assumes that these values are transferable from the experimental study to the analysis of folded protein structures [76]. Residue-based comparison of the proteins was carried out according to a manually-checked structural alignment with DALI (Figure S8 in Supporting Information S1). Domains A, B and C of AHA are composed of residues 1–86/147–356, 87–146 and 357– 448, respectively, whereas domains A, B and C of PPA are composed of residues 1–98/170–404, 99–169 and 405–496. Amino acid numbers in the text are labelled with a superscript A or P, referring to the AHA or PPA sequences, respectively. The 8 b-strands and the 8 a-helices, composing the (b-a)-barrel of domain A, are referred to as b1 to b8 and a1 to a8 according to their order in the sequence. Protein dynamics The essential dynamics (ED) analysis reveals high-amplitude concerted motions in the equilibrated portion of the trajectories, based on the diagonalization of the Ca covariance matrix of the atomic positional fluctuations [77]. About 10 eigenvectors are required to obtain the description of 70% of the variance in the simulated proteins. The collection of the selected eigenvectors describing the collective motions is termed essential subspace and can describe protein motions at a reasonable level of accuracy. PLoS ONE | www.plosone.org 10 September 2011 | Volume 6 | Issue 9 | e24214