was further developed by Miyazawa and Jernigan (MJ) [2,3] in a mean-field spirit of an effective medium. Over the years, a number of knowledge-based contact potentials [1–11] have been re-examined and redeveloped to understand the folding dynamics of a protein. For example, Betancourt and Thirumalai [7] have examined the MJ contact matrix and the potential matrix by Skolnick et al. [11] and proposed their own contact potential matrix (BT). By selecting an appropriate reference solvent (Thr) within the Miyazawa and Jernigan scheme [2,3], they find [7] that the BT interaction matrix gives ‘hydrophobicities that are in very good agreement with experi- ment.’ Bastolla et al. (BFKV) [8] have examined some of these knowledge-based interaction potentials and presented a scheme to guarantee optimal stability for most representative structures. They have pointed [8] out that ‘the optimized energy function guarantees high stability and a well-correlated energy landscape to most representative proteins in the PDB database.’ We have recently implemented [33] the classical MJ interaction matrix to examine the thermal response of two histones (H3.1, H2AX) [36,37]. These proteins are of comparable size (with 136 and 143 residues, respectively) but respond very differently to temperature. Whereas H3.1 exhibits a systematic transition from a random coil to globular conformation (see below), H2AX shows non-mono- tonic dependence (expanded conformations followed by compact structures) with a maximum as a function of temperature [33]. Because of the phenomenological nature of the knowledge-based interactions, it is worth re-analyzing the thermal response again with the tested and improved potentials such as BT [7] and BFKV [8] potentials in addition to classical MJ potential to understand similarity and differences in the conformational response to temperature. In context to extracting the optimal weight associated with the knowledge-based residue-residue contact matrix elements, Pokar- owski et al [35] have analyzed 29 contact potentials and concluded that ‘one-body approximations of the contact potentials could be useful in some applications.’ They have pointed out the ‘opportunities’ to develop different further types of potentials (perhaps multibody)’. Such an extensive analysis clearly under- scores the fact that the knowledge-based contact potentials are phenomenological measures and are somewhat adhoc and that the reliability of results about the general features and specific response should be carefully examined. In absence of compre- hensive analysis based on fundamental hierarchical interactions, knowledge-based interaction do provide a feasible mechanism to incorporate some specificity of residues in understanding the structure of a protein. We focus here on the conformation of histone H3.1 [36,37] as a function of temperature. In order to identify common results (reproducible by different potentials) and differences, we use three contact matrices as an input to a phenomenological potential (see below) to investigate the effect of temperature on the conformation of histone H3.1: 1M 2A 3R 4T 5K 6Q 7T 8A 9R 10 K 11S 12T 13G 14G 15K 16A 17P 18R 19K 20Q 21L 22A 23T 24K 25A 26A 27R 28K 29S 30A 31P 32A 33T 34G 35G 36V 37K 38K 39P 40H 41R 42Y 43R 44P 45G 46T 47V 48A 49L 50R 51E 52I 53R 54R 55Y 56Q 57K 58S 59T 60E 61L 62L 63I 64R 65K 66L 67P 68F 69Q 70R 71L 72V 73R 74E 75I 76A 77Q 78D 79F 80K 81T 82D 83L 84R 85F 86Q 87S 88S 89A 90V 91M 92A 93L 94Q 95E 96A 97C 98E 99A 100Y 101L 102V 103G 104L 105F 106E 107D 108T 109N 110L 111C 112A 113I 114H 115A 116K 117R 118V 119T 120I 121M 122P 123K 124D 125I 126Q 127L 128A 129R 130R 131I 132R 133G 134E 135R 136A. Random Coil to Globular Thermal Response (H3.1) Model and Methods The histone H3.1 consists of 136 residues in a unique sequence [36,37]. It is represented [31–33] by 136 nodes tethered together by the fluctuating bonds [38] on a cubic lattice in our coarse- grained description. A node (residue) is represented by a unit cubic cell (occupying its eight lattice sites) with excluded volume constraints [39]. The bond length between consecutive nodes can vary between 2 and !10. Such a bond fluctuation description is known to capture the computational efficiency while incorpo- rating ample degrees of freedom and extensively used in complex polymer systems [39], multi-component nano-composites [39,40] and protein chains [31–33,41]. Note that the lattice model with fluctuating (i.e., expanding and contracting) covalent bonds between consecutive residues has more degrees of freedom that the minimalist HP model used for sensitivity test by Betancourt and Thirumalai [7]. The large-scale simulations become feasible with such simplified coarse-grained representation of a residue without the all-atom details while capturing the specificity of each residue via residue-residue interactions. Each residue interacts with the neighboring residues within a range (rc) with a generalized Lennard-Jones potential, " Uij ~ Deij D s rij 12 s zeij rij 6# , rij vrc where rij is the distance between the residues at site i and j; rc = !8 and s = 1 in units of lattice constant. The potential strength eij is unique for each interaction pair with appropriate positive (repulsive) and negative (attractive) values (see below). The Metropolis algorithm is used to move each tethered residue randomly. For example, a residue, for instance, at a site i, is selected randomly to move to a neighboring lattice site, j. As long as the excluded volume constraints and the limitations on changes in the covalent bond length are satisfied, the residue is moved from site i to site j with the Boltzmann probability exp(2DEij/T), where DEij = Ej – Ei is the change in energy between its new (Ej) and old (Ei) configuration and T is the temperature in reduced units of the Boltzmann constant and the energy (eij). A unit Monte Carlo step (MCS) is defined as attempts to move each residue once. During the course of simulation, we keep track of a number of local and global physical quantities including energy of each residue, its mobility, mean square displacement of the center of mass of the protein, radius of gyration and its structure factor. Simulations at each temperature are performed for a sufficiently long time (typically for ten million time steps) with many independent samples (typically 150 samples for long runs and 1000 samples for short runs) to estimate these quantities. The data presented here are generated on a 643 lattice although different lattice sizes are used to assure that there is no finite size effect on the qualitative variations of the physical quantities and our conclusions. As mentioned above, we use three knowledge-based contact matrices (figure 1), i.e., the classical MJ [2], BT [7], and BFKV [8] for the residue-residue pair interaction (eij). On a first look, these matrices appear somewhat similar in general apart from the difference in magnitudes. There are however differences that are easier to spot with a closer look, e.g. elements 1–10, 190–200, etc., which may show in the final results on the thermal response (see below). Results Some insight into the structural relaxations and equilibration of the conformation at local and global scales of protein can be PLOS ONE | www.plosone.org 2 November 2012 | Volume 7 | Issue 11 | e49352