Molecular dynamics simulations of the thermal stability of tteRBP and ecRBP

Molecular dynamics simulations were performed for investigating the thermal stability of the extremely thermophilic Thermoanaerobacter tengcongensis ribose binding protein (tteRBP) and the mesophilic homologous Escherichia coli ribose binding protein (ecRBP). The simulations for the two proteins were carried out under the room temperature (300 K) and the optimal activity temperature (tteRBP 375 K and ecRBP 329 K), respectively. The comparative analyses of the trajectories show that the two proteins have stable overall structures at the two temperatures; further analyses indicate that they both have strong side-chain interactions and different backbone flexibilities at the different temperatures. The tteRBP 375 K and ecRBP 329 K have stronger internal motion and higher flexibility than tteRBP 300 K and ecRBP 300 K, respectively, it is noted that the flexibility of tteRBP is much higher than that of ecRBP at the two temperatures. Therefore, tteRBP 375 K can adapt to high temperature due to its higher flexibility of backbone. Combining with the researches by Cuneo et al., it is concluded that the side-chain interactions and flexibility of backbone are both the key factors to maintain thermal stability of the two proteins. An animated Interactive 3D Complement (I3DC) is available in Proteopedia at http://proteopedia.org/w/Journal:JBSD:22

The thermophilic Thermoanaerobacter tengcongensis ribose binding protein (tteRBP) was isolated from an extremely thermophilic bacterium Thermoanaerobacter tengcongensis which was isolated from a Chinese hot spring (Yanfen, Yi, Liu, Ma, & Zhou, 2001). The mesophilic homologous protein Escherichia coli ribose binding protein (ecRBP) was derived from the E. coli periplasm. The backbone structure of tteRBP (optimal activity temperature is 375 K) and ecRBP (optimal activity temperature is 329 K) (Cuneo, Tian, Allert, & Hellinga, 2008) is very similar; the root mean square deviations (RMSD) of the fit of the two all-atom's structures is 0.58 Å and have 76% similarity of the amino acid sequences (as Figure 1) (Cuneo et al., 2008), only some loops and turns that connect alternating β-strands and α-helices have modest differences, but significantly different in thermal stability.
The tteRBP and ecRBP are both members of the periplasmic binding protein (PBP) superfamily. PBP plays many roles in prokaryotic ABC transport (Boos & Shuman, 1998), chemotaxis (Davidson, Shuman, & Nikaido, 1992), and intercellular communication systems (Neiditch et al., 2006). The PBP fold consists of two three-layered α/β/α sandwich motif domains which are linked by a flexible hinge. The hinge which contains two or three β-strands divides the protein structure into split shaped and controls the bending of the two domains.
In order to investigate the thermal stability of tteRBP and ecRBP, Cuneo et al., (2008) used the method of cloning overexpression and purification to experimentally determine the X-ray crystal structure of tteRBP and the optimal activity temperatures of tteRBP and ecRBP by linear extrapolation, and further compared the hydrogen bonds, the salt bridges, and the divergence in overall structures; they speculated that the amino acid side-chain interactions can contribute to the different thermal stability of the two proteins. Now, we study the dynamics changes of the two proteins to further understand the mechanisms of thermal stability by means of molecular dynamics (MD) simulation. The results show that the overall structures of the ecRBP and tteRBP are well maintained at the two temperatures, respectively. The radius of gyration, solvent accessible surface area, salt bridges, side-chain interactions, and hydrogen bonds of the two proteins are almost unchanged at the different temperatures, their conformations also have no obvious changes; moreover, the root-mean-square inner product (RMSIP) reflects that the overall motion of tteRBP or ecRBP move along the same direction of the essential fluctuations at the different temperatures, respectively. The further analysis shows that the tteRBP and ecRBP both have strong side-chain interactions to maintain structural stability at the high temperatures, and they have different patterns of the intramolecular motion and flexibility at the two temperatures. At the room temperature and the optimal activity temperature, the tteRBP has higher flexibility and more high flexible regions than the ecRBP. These are the most differences between tteRBP and ecRBP. So the higher flexibility of tteRBP 375 K is more advantageous to maintain overall structure stability and adapt to high temperature by the fine-tune structure with higher flexibility. According to literature of Matthew J. Cuneo et al., we can conclude that the strong side-chain interactions and flexibility of backbone are the key factors to maintain thermal stability of the two proteins. In addition, the simulations of the tteRBP 387 K and ecRBP 345 K (387 and 345 K are the optimal activity temperatures of the ribose complex) (Cuneo et al., 2008) also further verify the above results.

Simulation system
Our MD simulations are based on the X-ray structure from the Protein Data Bank (PDB code: 2IOY (Alexander et al., GROMACS Tutorial for Drug-Enzyme Complex) (tteRBP) and 2DRI (Bjorkman et al., 1994) (ecRBP)). The geometry of ribose was optimized to default tolerances at the B3LYP/6-31G ⁄⁄ level of theory using Gaussian03 (Frisch et al., 2004) (Table S1). The restrained electrostatic potential (RESP) method was used for charge fitting for each atom. The tteRBP and ecRBP have a total charge between À7 and À2. The remaining force field parameters including bond, angle, torsional angle, and VDW radii were taken from the generalized Amber force field (Cornell et al., 1995). The tteRBP and ecRBP systems were solvated in a bath of TIP3P water molecules (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983) and placed in a 4.828 Â 4.358 Â 6.346 nm 3 and 5.053 Â 6.133 Â 3.886 nm 3 cubic boxes with 20,930 and 17,070 water molecules, respectively. Also to keep charge at neutrality, seven or two sodium ions were subsequently added to the tteRBP and ecRBP system.

MD simulations
The simulations were carried out in the NPT ensemble using GROMACS 3.3.1 (Berendsen, Van der Spoel, & van Drunen, 1995;Lindahl, Hess, & van der Spoel, 2001;Spoel et al., 2005) with AMBER-03 all-atoms force field (Duan et al., 2003). We do six MD simulations, namely tteRBP 300 K, tteRBP 375 K, tteRBP 387 K, ecRBP 300 K, ecRBP 329 K, and ecRBP 345 K. The temperatures were kept constant by coupling to a Berendsen thermostat with a coupling time of T = 0.1 ps (Berendsen, Postma, van Gunsteren, Dinola, & Haak, 1984). The pressure was kept constant at 1 bar by a compressibility of 4.5 Â 10 À5 bar À1 and coupling to a Berendsen barostat with p = 1.0 ps. All bonds were constrained using the LINCS algorithm (Hess, Bekker, Berendsen, & Fraaije., 1997), whereas the water molecules were constrained using SETTLE algorithm (Miyamoto & Kollman, 1992). An integration step of 2 fs was used. Nonbonded interactions were calculated using a cutoff of 1 nm. Long-range electrostatic interactions were calculated by Particle mesh Ewald summation (Darden, York, & Pedersen, 1998;Essmann et al., 1995) with grid spacing of 0.12 nm and cubic interpolation. Initially, all the heavy atoms of the complex of protein-ribose and the crystallographic water oxygen atoms were restrained to their crystallographic positions using a harmonic potential with a force constant of k = 1000 kJ mol À1 nm À2 , whereas surrounding TIP3P water molecules were first minimized. Then, the resulting system's energy was minimized without any restraints, the counterions was introduced and after another minimization, the goal temperatures were achieved by 100 ps uniform gradually heated from 300 K. Then, an additional 1.0 ns simulation was run to equilibrate the entire system at this temperature. The MD simulations were then performed for 40 ns for all six systems.

Essential dynamics analysis
Essential dynamics (ED) is a method of conformation analysis, which is based on the diagonalization of the covariance matrix (Amadei, Linssen, & Berendsen, 1993;Van Aalten, Findlay, Amadei, & Berendsen, 1995;Van Aalten et al., 1996), built from atomic fluctuations relative to their average positions of the MD simulation trajectories, where X represents the x-, y-, and z-coordinates of atoms, and < X > the average positions of the coordinated overall structure after they have been superimposed onto a reference structure to remove overall translational and rotational motion. Any subset of atoms is suitable for this analysis. Normally, however, only C α or backbone atoms are used since these subsets capture most of the conformational changes in the protein. The covariance matrix is diagonalized, yielding a set of eigenvalues and eigenvectors. The eigenvectors are directions in a 3Ndimensional space (where N is the number of atoms), and motion along single eigenvectors that correspond to concerted fluctuations of atoms. The eigenvalues repre-sent the total mean square fluctuation of the system along the corresponding eigenvectors. The eigenvectors are sorted by the size of their corresponding eigenvalues, the 'first' eigenvector being the one with the largest eigenvalue. In protein, there are always only a few ('essential') eigenvectors with large eigenvalues. Hence, the overall internal motion of the protein can be adequately described using only a few degrees of freedom.
We assume that the more biologically significant motion would correspond to the eigenvectors with larger eigenvalues. In this paper, we are interested in the intermolecular concerted motion. The similar degree of correlated motion in the two proteins at the different temperatures was evaluated by comparing the principal subspace (first 10 eigenvectors) of each structural trajectory by using the RMSIP (Amadei et al., 1999) between the first 10 eigenvectors. The RMSIP is defined as where g a i and v b j ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 10 are the ith and jth eigenvectors of the two different sets, respectively. The value of RMSIP equal to 0 means that the correlated motion of the two proteins are completely different, on the contrary equal to 1 means that the correlated motion are exactly the same.

Results and discussion
Simulation equilibrium and overall structural stability After four 40 ns MD simulations, the root mean square deviations (RMSD) of the two protein (tteRBP and ecRBP) backbones, calculated after a mass weighted superposition on the starting structure, are shown in Figure 2. It is obvious that the four systems can easily reach equilibrium phase after 15.0 ns and converged to values about 0.0956 nm (tteRBP 300 K), 0.1128 nm (tteRBP 375 K), 0.0877 nm (ecRBP 300 K), and 0.1294 nm (ecRBP 329 K). The values show a close similarity to the starting structure. The values of the cosine content (Berk Hess, 2002) (tteRBP 300 K, tteRBP 375 K, ecRBP 300 K, and ecRBP 329 K is 0.0358253, 0.0526691, 0.0450956, and 0.0299864, respectively) reflect that these trajectories are well convergence. These trajectories can be used for further analysis. In addition, the RMSD of tteRBP 387 K and ecRBP 345 K are shown in Figure S1. As shown in Table 1, at different temperatures, the radius of gyration (Rg) does not show significant changes for tteRBP and ecRBP, this reveals that the respective protein can keep overall structural stability; solvent accessible surface area (SASA) for tteRBP and ecRBP also do not obviously change, it means that almost no changes occur at the protein surface. Moreover, the comparative analyses of the hydrogen bonds of protein-protein and protein-ribose for tteRBP and ecRBP show that the overall structure of two proteins is very similar at different temperatures. Besides the hydrogen bonds of protein-solvent have apparent changes, this may be related to that at the high temperature, the water molecules cannot form stable hydrogen bonds with the protein, and this cannot affect the structure of the two proteins. In summary, Table 1 shows that at two temperatures, the two proteins well maintained the overall structural stability, respectively.
Then, we analyze the interactions between protein and substrate in the four MD simulations. The contacts were found by LPC program which bases on the surface complementarily approach (Sobolev, Wade, Vriend, & Edelman, 1996). The tteRBP and ecRBP has 41 and 43 contacts, respectively. The average distances of these contacts were calculated from 15 to 40 ns of the MD simulation trajectories. As shown in Table 2 (Due to lack of space, all data are not all listed), comparative analyses of tteRBP 300 K and tteRBP 375 K, the contacts of the residue ALA136 (has three contacts) with the substrate obviously enhanced, the contacts of the residue ASP215 (has six contacts) with the substrate slightly enhanced, but the contacts of the residue PHE164 (has seven contacts) with the substrate slightly reduced. Generally speaking, the interactions of residues-substrate have no significant changes. The comparative analysis of ecRBP 300 K and ecRBP 329 K also have the similar result.

Analysis of structural changes
In order to find the structural changes, first, we comparatively analyze the respective conformational changes of tteRBP and ecRBP at different temperatures by the free energy landscapes of four trajectories. We can see that the protein tteRBP and ecRBP have different free energy changes at two temperatures ( Figure 3). The tteRBP 300 K has the two conformational regions which can be transformed into each other, but at 375 K has the only one conformational region. The conformational changes of ecRBP are not obvious at 300 and 329 K which have the only one conformational region. We compare the representative structure of tteRBP and ecRBP in the conformational regions with their corresponding crystal structure, respectively. We found that there are only a few differences at the amino acid side-chain locations, such as the surface residues 75, 94, 122, 153, and 225 of the ecRBP and the surface residues 56, 205, and 256 of the tteRBP (these residues are marked in blue) at the different temperatures, but the most of the residues and the overall structure have not obvious diversity. As an example, the RMSD of the fitting of all-atom tteRBP and ecRBP representative structures 1-5 with their corresponding crystal structure are 0.86, 0.76, 0.86, 0.84, 0.98 Å, respectively. These values also illustrate that these structures are very similar. Thus at the different temperatures, the two proteins did not undergo obvious conformational changes and the overall structure is stable. Second, the side-chain interactions of the two proteins have not evident changes at two temperatures ( Figure S2), so the global locations of the two protein side-chain residues are stable. Moreover, the tteRBP 387 K and ecRBP 345 K have similar results compared to low temperatures ( Figure S3). At the same time, we find that the salt bridge interactions are consistent with the changes of the side-chain interactions. Therefore, the tteRBP and ecRBP keep stable side-chain and salt bridge interactions at the different temperatures, respectively.
On the above analyses, we found that the free energy landscapes, the side-chain interactions and the salt bridge interactions all prove that at the different temperatures, the proteins ecRBP and tteRBP keep the overall structure and locations of side-chain residues stable. It means that from low to high temperature, they have strong side-chain interactions to maintain the overall structure. The free energy landscapes suggest that at the different temperatures, the two protein structure has some little conformational adjustments. In order to gain insights into the structural and conformational changes, we further analyze the flexibility and essential dynamics of the two proteins.

Flexibility analysis
Protein flexibility is directly related to the protein conformation and function. So we analyze the root mean square fluctuations (RMSF) which is an important parameter to show the protein flexibility information. The RMSF of four MD simulation trajectories is shown in Figure 4. We can see that the ecRBP or tteRBP has similar distributions of flexible regions and flexibility is increased from low to high temperature, but the extent and the increasing pattern are different. Combining with Figure 4 and the structure of tteRBP 300 K, tteRBP 375 K, ecRBP 300 K, and ecRBP 329 K shows that the two proteins locally adjust their structure by the individual residues enter and out of the α-helix or β-sheet.
At the same time, when ecRBP 329 K is compared with ecRBP 300 K (Figure 4(b)), the overall C α flexibility also increased. But the range of increase is far less than tteRBP. Only the flexibility of the residues N terminal-13, 28-45, 55-58, 92-95, and 235-250 has significantly increased and these residues are located in the up domain and a helix of the down domain ( Figure 5(b)). The number of the residues which have obvious flexible changes is less than the half of the number of tteRBP flexible increased residues. The configurations have minute adjustments and stable overall structure at the different temperatures also. Such as, the residues 226 into α-helix and residues 14, 51-53, 152, 184-185, and 259-260 dissociation from the α-helix, residues 28-29 dissociation from β-sheet and no residues enter to the β-sheet.
Then, we compare the flexibility of the tteRBP 300 K vs. ecRBP 300 K and tteRBP 375 K vs. ecRBP 329 K (Figure 4(c) and Figure 4(d)) find that the tteRBP have higher overall flexibility, especially the down domain residues 110-152 of Figure 4(c) and 105-180 of Figure 4 (d). When comparing the amino acid sequences of tteRBP and ecRBP in the residues 105-180, we find that there is a concentration of 41 different residues. So we extract the residues 105-180 from the representative structures 1-5 of Figure 3 and fit them each other (the RMSD are shown in Table 3. At the same time, the hydrogen bonds among the residues 105-180 of tteRBP 300 K, tteRBP 375 K, ecRBP 300 K, and ecRBP 329 K are 53. 88, 52.36, 53.66, and 53.97. These two data both reflect that the residues 105-180 of the two proteins which have the most obvious differences in flexibility   have stable overall structure. The result may be the reason that tteRBP adapt to extreme high temperature.
In addition, the tteRBP 375 K and ecRBP 329 K have higher overall flexibility than the tteRBP 300 K and ecRBP 300 K, respectively. At room temperature or at the optimal activity temperatures, tteRBP has higher overall flexibility than ecRBP. These are the most obvious structural differences between tteRBP and ecRBP. The tteRBP and homologous ecRBP have significant different thermal stability, so the flexibility of backbone is directly related to the thermal stability of the two proteins. In addition, at the two temperatures, the two proteins have the stable interactions of residues-substrate, these show that the flexibility play role in the complexes stability, such as the high flexible residues of the tteRBP 136, 164, and ecRBP 13 still maintain strong interactions with substrate. The higher flexibility of tteRBP 375 K and ecRBP 329 K have more advantages in maintaining the overall structural stability by fine-tuning their structure, especially the former.
We analyze the RMSF of the tteRBP 375 K vs. tteRBP 387 K ( Figure S4(a)) and ecRBP 329 K vs. ecRBP 345 K ( Figure S4(b)) to further confirm the above result. From Figure S4 (a), We can see that the overall C α flexibility of tteRBP has obvious increase. The higher flexible regions consist of the residues N terminal-42, 62-78, 89-110, 120, 133-143, 155-172, 195-208, and 251-258. But, Figure S4(b) shows that the overall C α flexibility of ecRBP shows no changes, only few residues flexibility increase or decrease. These verify that the two proteins have different patterns of C α flexibility and different flexible distribution.

Essential dynamics analysis
We implement ED analysis to gain further information about C α flexibility and conformation. Table 4 shows the first five eigenvalues of the four diagonal molecular dynamics simulation trajectory covariance matrixes. It can be seen that the first three eigenvectors represent the molecular motion that dominates the overall molecular internal motion. So we focus on the first three eigenvectors.
We firstly look at the motion of tteRBP described by the first three eigenvectors, as shown in Figure 6(a, c, e) and (b, d, f) are the first, second and third eigenvectors of the tteRBP 300 and 375 K, respectively). The overall intensity of these movements are not particularly strong, especially the motion described by the first tteRBP eigenvector. The differences of the both overall molecular internal motion described by the first tteRBP 300 K and 375 K eigenvector are not obvious (Figure 6(a) and (b)). Comparing the tteRBP 375 K and the tteRBP 300 K, the motion intensity of residue 121 is slightly weaker, the residues 27-32, 38-45, and 149-153 are obviously stronger. The residues 131-134, which are the part of the protein binding pocket in an opening state have almost the most strong movement, but the same region no obvious movement at 300 K.
From Figure 6(c) and (d), we know that the internal motion described by the second tteRBP eigenvector is not significantly lower than the first eigenvector, but they have more obvious differences. The tteRBP 375 K has more obvious strong motion than tteRBP 300 K. At 300 K, the C terminal residues 270-274 have strong movements and the other parts of protein have obvious movements, especially the N terminal domain has no visible movement. At 375 K, not only C terminal residues 268-274, but also residues 53-58, 131-136 (the part of the protein binding pocket in an opening state), 149-154, and 198-208 have strong movements, most of the other regions and N terminal have obvious movements. It can be seen that the internal motion enhanced obviously, especially in the opening of protein binding pocket at the 375 K.
Both of the internal motions described by the third tteRBP eigenvector have more differences (Figure 6(e) and (f)). At 300 K, only some residues near C terminal have obvious movements. At 375 K, the overall motion compared that the second eigenvector has not significantly reduced. Most of residues have obvious movements and some residues near C terminal and residue 178 have strong movements, especially the residues 132-139 have visible movements. The above comparison shows that the overall internal molecular motion has obviously increased at 375 K, especially in the part of protein binding pocket which has high flexibility. Table 3. The RMSD of the fit of the residues 105-180 which is extracted from the representative structure ( Figure 3   Then we look at the motion of ecRBP described by the first three eigenvectors, Figure 7(a, c, e) and (b, d, f) are the first, second, and third eigenvectors of the ecRBP 300 K and 329 K, respectively). The motion of the ecRBP 300 and 329 K described by the first, second, and third ecRBP shows obvious gradual decrease. We can see that the overall molecular internal motion described by the first ecRBP eigenvector is obviously stronger than the motions described by the first tteRBP eigenvector (Figure 7(a) and (b)). Comparing the ecRBP 329 K with ecRBP 300 K, the N terminal and the residues 53-58, 72-80, 94-98, 122-124, 147-153, 174-177, 221-228, and 240-271 have more strong movements, especially the residues 241-244. Figure 7(c) and (d) show the overall molecular internal motions described by the second ecRBP eigenvector is significantly decreased. The motion intensity of N terminal and the residue 53-58, 76-82, and 223-227 and C terminal at 329 K is much more than that at 300 K.
Both the overall motions described by the third ecRBP eigenvector are continually reduced (Figure 7(e) and (f)). At 300 K, only the residues 53-55 have strong motion and the other part of protein has slight motion. At 329 K, the C terminal and the residues 54-57 have strong movement. The residues 179-207 have slightly weaker motion and a lot of residues have obvious motion. The above comparison shows that the overall motion has distinct increase at 329 K.
The tteRBP 375 K and ecRBP 329 K have higher overall motion than tteRBP 300 K and ecRBP 300 K, respectively. At the same time, the tteRBP and ecRBP have different patterns of intramolecular motion. This is consistent with the flexibilities of tteRBP and ecRBP, which have different C α flexibility increase and distributed regions at the different temperatures.
At the same time, in order to verify these results, we calculate the values of B-factor and RMSIP. As shown in Figure 8, the tteRBP and ecRBP both have higher Bfactor values at the high temperatures. The tteRBP 375 K (tteRBP 300 K) has higher B-factor value than ecRBP 329 K (ecRBP 300 K), especially the atoms 120-180. This result is consistent with the flexibility analysis and ED analysis. Then the RMSIP of the tteRBP 375 K vs. tteRBP 300 K and ecRBP 329 K vs. ecRBP 300 K is 0.770266 or 0.824369, respectively, showing that the overall motion of tteRBP or ecRBP move along the same direction of the essential fluctuations at the different temperatures. It can be concluded that the two proteins keep the motion modes from low to high temperature; these also imply that the conformational changes of the two proteins highly overlap at the different temperatures, respectively. Moreover, the values of the RMSIP which are shown in Table 5 illustrate the similar results.
Combination with the above flexibility analysis, we can infer that the two proteins must depend on their own changes of backbone flexibility to maintain and even enhance their function at high temperature.

Conclusion
From low to high temperature, the two proteins well maintained their overall structure; especially the evidences of the RMSIP support this. Further analyses indicated that the interactions of salt bridges and the side-chain residues of the tteRBP or ecRBP have not obviously changed. Even the comparisons of the sidechain residues interactions of the tteRBP 387 K or ecRBP 345 K with their corresponding to low-temperature proteins are the same. These data show that the two proteins have stable locations of the side-chain residues. At the two temperatures, the tteRBP or ecRBP has stable overall structure and locations of the side-chain residues; it means that the two proteins have the strong side-chain interactions to maintain structure stability at the high temperatures. This is consistent with the view of the research of Cuneo et al. (2008).  At the same time, the tteRBP 375 K and ecRBP 329 K have higher C α flexibility than their 300 K, respectively. At the two temperatures, the flexible regions of tteRBP or ecRBP are similar. Simultaneously, although some minute adjustments, the secondary structures and overall structures of the two proteins are still stable. From low to high temperature, tteRBP have higher flexibilities than ecRBP, and the RMSF of tteRBP 387 K or ecRBP 345 K also prove this. The different C α flexible distributions are the most obvious structural differences between tteRBP and ecRBP. The further information of ED shows that the tteRBP 375 K or ecRBP 329 K has higher intensity motion than tteRBP 300 K or ecRBP 300 K; this is consistent with the higher flexibility of tteRBP 375 K or ecRBP 329 K. At the two temperatures, the tteRBP and ecRBP have different patterns of the increased intramolecular motion; this should be the underlying reason that the two proteins have different flexibility. Therefore the ecRBP 329 K and tteRBP 375 K which have increased overall flexibility can better adjust favorable conformation to maintain the overall structural stability at the high temperature, especially the tteRBP 375 K.
At the different temperatures, the two proteins have strong side-chain interactions to assist the stable global structure. Furthermore, the comparison of the two homologous protein flexibilities illustrates that the higher flexibility maintains the thermal stability. So the sidechain interactions and flexibility of backbone are the key factors that contribute to thermal stability of the proteins tteRBP and ecRBP.

Supplementary material
The supplementary material for this paper is available online at http://10.1080/07391102.2012.721497.