Simple approach for ranking structure determining residues

Statistics

Principal component analysis was performed using the R software (R Core Development Team, 2013). The thermal unfolding experiments of the 6aJL2 mutants were repeated in triplicate.

Crystallographic structures

The PDB codes, size, and protein class of the selected structures are as follows: chymotrypsin inhibitor (2CI2; 65 residues; 16% α + 21%β); 6aJL2 (2W0K; chain A, 111 residues; 5% α + 46%β); apoflavodoxin (1FTG; 168 residues; 35% α + 19%β); dimeric Arc repressor (1ARR; 108 residues; 62% α + 9%β); DNA-binding domain of human estrogen receptor α in complex with the estrogen response element DNA duplex (1HCQ; chains A, B; 128 residues; 26% α + 9%β, and chains C and D, 56 nucleotides); cold shock protein from Bacillus subtilis (1CSP; 67 residues; 4% α + 55%β); cold shock protein from B. caldolyticus (1C9O, 66 residues; 4% α + 62%β); and the largely disordered N-terminus of suppressor of cytokine signaling 5 mammalian suppressor of Janus Kinase interaction region (2N34, 70 residues; 12% α). Minor structural modifications were performed on PDB files such as fulfilling N-terminal domains and incomplete lateral chains by using Swiss PDB software (Guex & Peitsch, 1997). Furthermore, all nodes were renumbered in continuous order avoiding repetitions; i.e., in PDB 2CI2 the first residue begins as residue 19. Consequently, all residues were systematically renumbered so the first residue was identified as residue 1 and so forth. Solvent-exposed area was calculated with the NOC software (Chen, Cang & Nymeyer, 2007).

Grade

Residues can be regarded as nodes and contacts as edges. Hence, an edge was defined when any two non-hydrogen atoms from a pair of residues are within distance of 5 Å. To study the topology of the residue contact network, we measured the degree of node-i, Ki, as the number of neighbors of node-i. Chain A of crystallographic structures was selected. In the case of quaternary structures, all the structure was considered to measure the grade of each node but only chain A was selected for further comparisons.

Molecular dynamics

The query protein was prepared using the Protein Preparation Wizard in Maestro 9.2 package (Schrödinger LCC, NY, USA) and included in a 10-Å water box that contained 14513 SPC-type water molecules for 2CI2, 16667 for 6aJL2, 20051 for 1FTG, 20990 for 1ARR, 28142 for 1HCQ, 12219 for 1CSP, and 12033 for 1C9O. Neutralizing ions were added, and other metal ions already present in the protein structure were left at the same place. The simulated annealing calculations and data analysis were conducted using the Desmond and Maestro programs, respectively (Maestro-Desmond Interoperability Tools, version 3.0; Schrödinger, NY, USA). The OPLS˙2005 force field was used for every molecular dynamics simulation. Cubic periodic boundary conditions were used for most of the proteins (53.9 × 53.9 × 53.9 ų for 2CI2, 56.8 × 56.8 × 56.8 ų for 6aJL2, 60.9 × 60.9 × 60.9 ų for 1FTG, 61.2 × 61.2.9 × 61.2 ų for 1ARR, 51.1 × 51.1 × 51.1 ų for 1CSP, and 50.8 × 50.8 × 50.8 ų for 1C9O); due to the size of the complex 1HCQ, rectangular cuboid boundary conditions were used (58.5 × 59 × 91 ų for 1HCQ). Each simulation was adjusted with an NPT ensemble by weak coupling to an external bath temperature at constant pressure of 1 atm and relaxation time of 2 ps, regulated by Berendsen barostat (Berendsen et al., 1984). All short-range interactions were computed using a 9 Å cutoff, and for long-range interactions (electrostatic and Van der Waals), a smooth particle mesh Ewald method with a tolerance of 1 × 10⁻⁹ was applied (Essmann et al., 1995 ). To ensure that our simulations started from local minima, a simulated annealing algorithm was performed. This method started the simulation at high temperature (400 K) to overcome thermodynamic and conformational barriers, followed by gradual cooling (annealing) to reach low energy regimes. It is widely used for the optimization of structures from experimental methods, comparative protein modeling, or studying the conformational dynamics of protein or peptide folding and unfolding (Mori & Okamoto, 2009). The full system was heated at 10 K for 30 ps, 100 K for 100 ps, 300 K for 200 ps, 400 K for 300 ps, and 400 K for 500 ps and then cooled to 298 K for 1,000 ps. To ensure that heating at 400 K did not affect protein or protein/DNA structure, the RMSD of heavy atoms (C, N, O, S, P) derived from the annealed structure was compared against corresponding crystallographic structure. If RMSD standard deviation was ≤1.0 Å, then it was assumed that the annealing algorithm did not changed the whole structure or denature it. In fact, since the displacement of heavy atoms above 1 Å is considered as a conformational change, the algorithm can be trusted in finding the local minimum. A lineal interpolation step between two adjacent time points was employed. After the sixth step, a production of 25 ns was achieved with an integration time of 1 fs. The entire analysis was performed using trajectory coordinates, and the energies were written to a disk every 1.2 ps. A frame was extracted every 0.25 ns throughout the simulation, and the overall frames were saved as a PDB file.

Dynamical entropy

An analytical algorithm encoded in the Perl language was developed to generate three files using the 100-frames PDB-file as the template. In the first step, a single file is generated for each model, and this single file contains the atom coordinates, the corresponding molecular weight, and the name of the node (amino acid or nucleotide base). The second step calculates the coordinates of the center of mass of the node. The third step calculates the distance between the center of mass of each pair of nodes and their normalized distance described by Eq. (1), (1)${\displaystyle \overline{d_{ij}}=\frac{d_{ij}}{r_{vd{W}_i}+r_{vd{W}_j}}}$ where $\overline{d_{ij}}$ is the normalized distance, d_ij is the distance between the center of mass of node i and that of node j, and r_vdW is the Van der Waals radius of the respective node considering all their atoms. For amino acid residues, values of their Van der Waals radiuses comprehending the whole residue were obtained from Darby & Creighton (1993). For nucleotide bases, values were obtained from Voss & Gerstein (2005). The mean distance between each pair of nodes of the 100 structures sampled was calculated. Since the next steps involve eigenvector and eigenvalue measurement properties, the mean value of the inverse normalized distance was calculated, so the largest weight represents the closest distance between a pair of nodes and the smallest weight represents the longest distance.

A weighted adjacency matrix A = (a_ij) ≥ 0 of size N × N, where N is the number of nodes, was constructed. In this case, matrix A is symmetric (a_ij = a_ji) and undirected. Following the mathematical strategy described in Demetrius & Manke (2004) we now assume that the stochastic process is given by a Markov Matrix P = p_ij where p_ij ≥ 0 and ${\mathrm{\Sigma }}_j{p}_{ij}=1$. The stationary distribution of matrix P is described by Eq. (2), (2)${\displaystyle \pi P=\pi }$ where π is defined as the left-hand eigenvector associated with the largest eigenvalue solved with Mathcad 15™ software. The dynamical entropy of this process of each node, H_i, is described by Eq. (3), (3)${\displaystyle H_i=-\sum _j{\pi }_i{p}_{ij}log{p}_{ij}}$ where the term p_ijlogp_ij is the standard Shannon entropy.

6aJL2 mutants

The synthesis of single-point mutants of 6aJL2 (Arg25His, Ile30Gly, Tyr36Phe, and Gln6Asn) was performed using recursive PCR (Prodromou & Pearl, 1992). The obtained DNAs were cloned into the pSyn1 expression vector (Schier et al., 1995). All of the constructions were verified by nucleotide sequencing (Sanger, Nicklen & Coulson, 1977). The variants were expressed in Escherichia coli BL21 (DE3) and purified as described previously (Del Pozo-Yauner et al., 2008). The protein purity was verified using SDS-PAGE electrophoresis, and the protein concentration was determined spectrophotometrically at 280 nm in 6.5 M GdnHCl and 20 mM sodium phosphate buffer, pH 7.5, using molar extinction coefficients calculated from the amino acid sequence using the ProtParam software, which is available at the ExPASy website (Artimo et al. al., 2012).

Unfolding

Samples containing 50 µg/ml of protein in phosphate-buffered saline (PBS), pH 7.5, were placed into a 3-ml quartz cuvette. Changes in the tryptophan fluorescence were measured using a LS50B Perkin Elmer spectrofluorometer with an excitation wavelength of 295 nm (2.5-mm bandwidth) and an emission wavelength of 355 nm (5-mm bandwidth). The temperature was increased from 298 to 350 K at a rate of 1 K/min and then samples were cooled to 298 K at the same rate. The data were analyzed using the thermal unfolding Eq. (4), which was obtained from Eftink (1995). (4)${\displaystyle F_{Trp}=\frac{\left(y_n+m_{n}T\right)+\left(y_d+m_{d}T\right)e^{\left(\frac{\mathrm{\Delta }{H}_m}{R{T}_m}-\frac{\mathrm{\Delta }{H}_m}{RT}\right)}}{1+e^{\left(\frac{\mathrm{\Delta }{H}_m}{R{T}_m}-\frac{\mathrm{\Delta }{H}_m}{RT}\right)}}}$ where F_Trp is the tryptophan fluorescence, T is the temperature, T_m is the temperature of the midpoint, $\mathrm{\Delta }{H}_m$ is the enthalpy at T_m, y_n and m_n describe the pre-transition phase, and y_d and m_d describe the post-transition phase. Non-linear regression was performed using the OriginPro8™ software. The change in the Gibbs energy of the wild-type versus the mutant $\left(\mathrm{\Delta }\mathrm{\Delta }G\right)$ was calculated for temperature- and denaturant-induced unfolding processes using the following equations.

Thermal unfolding (Becktel & Schellman, 1987), Eq. (5): (5)${\displaystyle \mathrm{\Delta }\mathrm{\Delta }G=\frac{{\mathrm{\Delta }H}_{m\mathrm{WT}}}{T_{m\mathrm{WT}}}\left(T_{m\mathrm{MUT}}-T_{m\mathrm{WT}}\right)}$ where WT refers to the wild-type values and MUT refers to the mutant values of the melting temperature T_m, and $\mathrm{\Delta }{H}_m$ is the enthalpy value.

Chemical unfolding (Creighton, 1990), Eq. (6): (6)${\displaystyle \mathrm{\Delta }\mathrm{\Delta }G=m_{\mathrm{WT}}\left(C_{m\mathrm{MUT}}-C_{m\mathrm{WT}}\right)}$ where m_WT is the transition slope of the wild-type, and C_m is the denaturant concentration at which $\mathrm{\Delta }G=0$.