Article Information

• PDF (501 kb)
• Full Text with Large Figures

PubMed

• Articles by Kei Moritsugu
• Articles by Jeremy C. Smith

Related Articles

• Plastocyanin conformation. An analysis of its near ultraviol...

  Plastocyanin conformation. An analysis of its near ultraviolet absorption and circular dichroic spectra Biophysical Journal, Volume 49, Issue 4, 1 April 1986, Pages 891-900 J.E. Draheim, G.P. Anderson, J.W. Duane and E.L. Gross Abstract The near-ultraviolet absorption and circular dichroic spectra of plastocyanin are dependent upon the redox state, solution pH, and ammonium sulfate concentration. This dependency was observed in plastocyanin isolated from spinach, poplar, and lettuce. Removal of the copper atom also perturbed the near-ultraviolet spectra. Upon reduction there are increases in both extinction and ellipticity at 252 nm. Further increases at 252 nm were observed upon formation of apo plastocyanin eliminating charge transfer transitions as the cause. The spectral changes in the near-ultraviolet imply a flexible tertiary conformation for plastocyanin. There are at least two charge transfer transitions at approximately 295–340 nm. One of these transitions is sensitive to low pH's and is attributed to the His 87 copper ligand. The redox state dependent changes observed in the near-ultraviolet spectra of plastocyanin are attenuated either by decreasing the pH to 5 or by increasing the ammonium sulfate concentration to 2.7 M. This attenuation cannot be easily explained by simple charge screening. Hydrophobic interactions probably play an important role in this phenomenon. The pH and redox state dependent conformational changes may play an important role in regulating electron transport. Abstract | PDF (1371 kb)

• Higher plants contain a modified cytochrome c6

  Higher plants contain a modified cytochrome c6 Trends in Plant Science, Volume 7, Issue 6, 1 June 2002, Pages 244-245 Jürgen Wastl, Derek S Bendall and Christopher J Howe Full Text | PDF (36 kb)

• Tracking the function of the cytochrome c6-like protein in h...

  Tracking the function of the cytochrome c6-like protein in higher plants Trends in Plant Science, Volume 8, Issue 11, 1 November 2003, Pages 513-517 Martin Weigel, Paolo Pesaresi and Dario Leister Abstract The contention that plastocyanin is the only mobile electron donor to photosystem I in higher plants was recently shaken by the discovery of a cytochrome -like protein in and other flowering plants. However, the genetic and biochemical data presented in support of the idea that the cytochrome homologue can replace plastocyanin have now been challenged by two complementary studies. This re-opens the debate on the real function(s) of cytochrome in the chloroplasts of higher plants. Abstract | Full Text | PDF (349 kb)

• …more

Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 95, Issue 4, 1639-1648, 15 August 2008

doi:10.1529/biophysj.108.131714

Biophysical Theory and Modeling

REACH Coarse-Grained Biomolecular Simulation: Transferability between Different Protein Structural Classes

Kei Moritsugu^{*, †, ‡} and Jeremy C. Smith^*, ,  

^* Center for Molecular Biophysics, University of Tennessee/Oak Ridge National Laboratory, Oak Ridge, Tennessee
^† Computational Molecular Biophysics, Interdisciplinary Center for Scientific Computing (IWR), University of Heidelberg, Heidelberg, Germany
^‡ Computational Science Research Program, RIKEN, Wako, Japan

Address reprint requests to Jeremy C. Smith, Center for Molecular Biophysics, University of Tennessee/Oak Ridge National Laboratory, 1 Bethel Valley Road, Oak Ridge, TN 37831. Tel.: 865-574-9635; Fax: 865-576-7651.

Molecular dynamics (MD) computer simulation is a useful tool for characterizing protein dynamics. The accuracy of an MD simulation depends on that of the protein three-dimensional structure providing the initial coordinates, on the force field defined by the potential energy function, and on the simulation methodology.

Atomistic force fields, such as CHARMM 1, AMBER 2, and GROMOS 3, have been widely used for protein MD studies. However, simulating biological systems on long timescales (e.g., microseconds) is difficult because of limitations of computational power. To reduce computational costs for large systems, coarse-graining methods, in which the number of degrees of freedom is reduced, have been derived. For example, a protein molecule can be represented by point residues centered at the C_α atoms.

An early simplified residue-based coarse-grained (CG) force field, the elastic network model (ENM), uses a nearest-neighbor harmonic approximation and enables collective vibrational normal modes to be rapidly calculated 4,5,6,7,8. Although extremely simple, ENM has been usefully applied in various studies of protein function 9,10,11,12,13. More sophisticated force fields, including virtual local (bond, angle, torsion, etc.) and nonlocal interactions, have been proposed for CGMD simulations 14,15,16,17,18.

Parameters for CG models have been calculated using the Boltzmann inversion of the MD-derived probability distribution 14,15 or by iteratively matching the internal coordinate fluctuations 16 or the forces on the CG sites 17 to all-atom MD. In another approach, the implicit degrees of freedom are separated dynamically from the CG degrees of freedom using a united-residue energy function 18. Although the above approaches all have merits and hold promise for determining dynamic aspects of macromolecular function using CGMD 14,15,16,17,18,19,20,21,22, the parameters of the CG model are determined a posteriori so as to reproduce, e.g., vibrational amplitudes. As a result, the derived parameters do not have clearly attributable physical origins, in contrast to all-atom force fields.

In previous work 23, we introduced a methodology, REACH (Realistic Extension Algorithm via Covariance Hessian), for deriving residue-scale ENM force constants from the variance-covariance matrix calculated from atomic-detail MD simulation. The REACH method is a direct mapping, requiring no iterative fitting and no input of experimental data. The REACH normal modes were shown to reproduce the amplitudes and frequencies of MD-derived fluctuations in myoglobin 23.

The aim of this study is to examine the transferability of the REACH CG force constants among protein molecules. To do this, the REACH method was applied to three model proteins of different structural classes 24: myoglobin (an α-fold of eight α-helices), plastocyanin (a β-fold of eight β-strands), and dihydrofolate reductase (DHFR) (an α/β-fold, comprising 4 α-helices and 10 β-strands). The structures of the three proteins are shown in Fig. 1.

Display large version of this figure

Figure 1

Three-dimensional structures of myoglobin, DHFR, and plastocyanin. Drawn with the program MOLSCRIPT 38.

The REACH force constants for the three proteins are calculated from 10-ns atomistic MD trajectories. The interactions within α-helices and β-strands, and the interactions between β-strands (forming β-sheets) are focused on. Analytical model functions are derived representing the distance dependence of the force constants. The REACH mean-square fluctuations are calculated and compared with those from the atomistic MD. In a further development, the method is extended from normal mode analysis to MD. CGMD simulations are performed with the REACH force field. The mean-square fluctuations, the vibrational densities of states, and the dynamic cross-correlation matrices are calculated from the atomistic and coarse-grained MD and compared. A generic, transferable analytical REACH force field is derived that reproduces well the fluctuations from proteins of different structural classes.

MD simulations of myoglobin, DHFR, and plastocyanin were performed. The starting structures were obtained from the Protein Data Bank 25: 1A6G 26, 1RX1 27, and 1PLC 28, respectively.

The model systems were constructed as follows: A rectangular primary simulation box was built of dimensions, 60 Å×60 Å×60 Å for myoglobin, 62 Å×52 Å×48 Å for DHFR, and 62 Å×54 Å ×52 Å for plastocyanin. Then 7181 TIP3P water molecules 29 were placed in the box for myoglobin, 3173 for DHFR, and 4413 for plastocyanin, together with one chloride ion for myoglobin, six potassium ions for DHFR and seven potassium ions for plastocyanin. This procedure constructed electrically neutral systems of 21,110 atoms for myoglobin, 12,094 for DHFR, and 14,712 for plastocyanin. Periodic boundary conditions were imposed on the boxes.

The simulations were performed using the program NAMD2 30. The CHARMM all-atom parameter set 22 31 was used for the potential function. Electrostatic interactions were calculated with a dielectric constant of 1 using the particle mesh Ewald method 32, i.e., without truncation. The systems were energy minimized with 1000 steps of the conjugate gradient method. Then, each simulation system was uniformly heated to the target temperature of 120K over 30ps and equilibrated for 100ps with velocity scaling in the NVE ensemble. Subsequent simulations were carried out at constant temperature (120K) and pressure (1atm) conditions (the NPT ensemble) for a 500-ps equilibrium run and 10-ns production runs. The temperature of 120K was selected so as to reduce the contribution of anharmonic protein motion: the previous study showed a clear decrease of the force constants with temperature corresponding to softening of protein dynamics 23. Subsequent work will focus on construction of a force field capable of accurately simulating anharmonic dynamics at physiological temperatures. Langevin dynamics was used to control the temperature and pressure 33. All production runs were performed for 10ns. The atomic coordinates were saved every 50fs for analysis.

The REACH CG method uses only C_α atom coordinates. Heteroatoms such as ions and ligands are not included. The mass of each CG residue is defined as the sum of the atomic masses comprising the corresponding residue.

The REACH method is described in detail in the literature 23. Here we give a brief summary. The potential energy in the ENM is written as follows 4:

(1)

(2)

(3)

(4)

In this study, the 10-ns MD trajectories were each separated into ten 1-ns trajectories, allowing the calculation of C from each 1-ns trajectory. Subsequently, the ten C matrices were averaged, and the associated force constants derived using Eqs. (3). To calculate the force constants for local interactions, i.e., the virtual 1−2 (between residues i and i+1), 1−3 (between residues i and i+2) and 1−4 (between residues i and i+3) interactions, segments of 20 residues were fitted individually to calculate submatrices of C. This avoids problems associated with the best fit to the overall structure arising from the incorporation of external motions of the segments, resulting in errors in the pairwise covariances 23.

An analysis was performed of the dependence of the REACH force constants on the structural fold, i.e., the secondary-structural elements. For this, the DSSP program 34 was used to define the residues comprising secondary structures. The interactions within α-helices and β-strands, and between β-strands, were calculated using the corresponding secondary structural segments: each secondary structure was fitted individually to calculate a submatrix of C, and then the force constants were derived. For α-helices, a virtual 1−5 (between residues i and i+4) interaction was included in the force field arising from the helical pitch of ∼3.6 residues.

Mathematical model functions describing the distance dependence of the force constants are useful for obtaining a simplified understanding of protein dynamics and for convenient application in CG simulations. Terms were constructed for each interaction type and include a dependence on the secondary-structural elements. The local (1,2, 1–3, 1–4, and 1–5) interactions and the inter-β-strand interaction were modeled with a single force constant. In contrast, the nonbonded interaction force constants were modeled as distance dependent with a double-exponential function. The resulting force constant parameters are as follows:

(5)

Using the potential energy (Eq. (1)) together with the equilibrium structure and the force constant models, residue-scale CG normal mode calculations and MD simulations were performed.

Normal-mode eigenvalues, and eigenvectors, were calculated by diagonalizing the Hessian matrices derived from the potential energy. The mean-square fluctuation of residue n, at temperature T is given by

(6)

CGMD simulations were performed using the program DL_POLY 35. For this, a Morse potential, was applied for the nonlocal interactions, replacing the harmonic potential in Eq. (1), i.e.,

(7)

Initial coordinates were taken as the C_α coordinates of the average structure over the 10-ns atomistic MD simulation. Both 1-ns equilibration and 1-ns production runs were performed at constant volume and temperature conditions. The time step applied was 50fs. Langevin dynamics 33 was used to maintain the constant temperature with a friction constant of γ=3ps⁻¹, a value chosen to be similar to the average friction exerted on the low-frequency motions < 100cm⁻¹ in atomistic solvated myoglobin at 120K 36. Nonlocal interactions were truncated at r_c<40 Å. The mean-square fluctuation of each residue was calculated from the trajectories derived and compared with the atomistic MD and CG normal-mode results.

The vibrational density of states, g(ω), was calculated from the MD trajectory as the Fourier transform of the velocity autocorrelation function; i.e.,

(8)

The dynamic cross-correlation matrix of atomic fluctuations consists of elements c_i,j defined as

(9)

The REACH force constants for the three proteins, myoglobin, DHFR, and plastocyanin, were calculated from the corresponding atomistic MD trajectories and are compared. The interactions within and between secondary-structural elements are also examined.

In Fig. 2, the REACH force constants, k, for interactions between the local residue pairs are plotted as a function of the pairwise distance, r. The distributions, k(r), of the virtual 1−2, 1−3, 1−4, and 1−5 interactions are very similar among the three proteins. The distributions of the 1−2 interactions are quite different from the others, with a larger k, peaked at∼3.83 Å, arising from virtual bonds between neighboring amino acids. The 1−2 interaction via cis peptide bonds has a smaller k and r than that via trans peptide bonds: cis peptide bonds occur at residues 95−96 in DHFR (k=491kJ/mol·Å² and r=2.97 Å), and residues 15−16 (k=612kJ/mol·Å² and r=3.05 Å) and 35−36 (k=587kJ/mol·Å² and r=3.11 Å) in plastocyanin.

Display large version of this figure

Figure 2

REACH force constants, k, of virtual 1−2, 1−3, 1−4, and 1−5 interactions in (1) myoglobin, (2) DHFR, and (3) plastocyanin are plotted as a function of pairwise distance, r. Interactions within α-helices and β-strands are plotted with circles and triangles, respectively. The virtual 1−2 interaction via cis peptide bonds (residues 95–96 in DHFR and residues 15–16 and 35–36 in plastocyanin) is excluded from the figure.

The local interactions within secondary structures are also shown in Fig. 2. Small differences in k(r) among the three proteins are caused by differences in interactions within the secondary-structural elements. The distribution of interactions within α-helices in DHFR is narrower than that in myoglobin because there is more variety in the structure and dynamics of the eight helices (comprising 113 residues) in myoglobin than in the four helices (comprising 32 residues) in DHFR. The 1−2 interactions within the α-helices are slightly larger than those within the β-strands. The average equilibrium distance, 〈r〉, is also slightly larger for α-helices (3.84 Å for both myoglobin and DHFR) than for β-strands (3.82 Å for DHFR and 3.81 Å for plastocyanin). For the 1−3, 1−4, and 1−5 interactions, the distributions of force constants within α-helices are narrower than those within β-strands, with 〈r〉 being 5.49, 5.16, 6.22 Å for myoglobin and 5.51, 5.13, 6.10 Å for DHFR, respectively. For α-helices, the 1−4 interaction has a smaller 〈r〉 than the 1−3 interaction because of the helical pitch of ∼3.6 residues. In contrast, the 1−3 and 1−4 interactions within β-strands have larger 〈r〉 and also more scattered distributions because of the extended structure of β-strands.

Fig. 3 shows k for the nonlocal interaction residue pairs as a function of r. The distributions of the nonlocal interactions, including the inter-β-interactions, are again very similar among the three proteins. Because of the hydrogen bonds between β-strands forming β-sheets, the distributions of the inter-β-interactions are found, on average, at shorter pairwise distances than those within β-strands, with 〈r〉 for the former being 4.83 Å for DHFR and 4.81 Å for plastocyanin.

Display large version of this figure

Figure 3

REACH force constants, k, of nonbonded interactions in (1) myoglobin, (2) DHFR, and (3) plastocyanin are plotted with dots as a function of pairwise distance, r. Interactions between β-strands are plotted with triangles. In the insets the averages of k within bins of 1 Å width are shown together with the corresponding fitted model function (dashed curve).

Model functions of the distance dependence of the force constants were constructed. The model functions of the local interactions and the inter-β-interactions, which have narrow k(r) distributions and no significant distance dependence, were assumed to be constant with distance, as in the previous REACH implementation 23 and another CGMD study 14. In contrast, the model functions for the nonbonded interactions were assumed to be a double-exponential in form. The slower decaying exponential term, which is the smaller in magnitude, was found to be essential for reproducing the atomistic MD fluctuations (results not shown). The resulting parameters and standard-deviation errors are listed in Table 1,Table 2. Again, the force constant parameters are similar among the three proteins. The force constant values obtained for interactions within secondary-structural elements for the three proteins are very close to each other, and indeed, the values obtained from the three proteins by averaging the interactions within each α- or β-element, i.e., k_12α or k_12β for, e.g., the 1−2 interaction, are closer to each other than the values averaged over all the 1−2 interactions (i.e., k₁₂).

With the REACH force-constant model functions listed in Table 1,Table 2, residue-scale CG normal-mode analyses and MD simulations were performed for the three proteins, myoglobin, DHFR, and plastocyanin, and the mean-square fluctuations, x², were calculated.

Fig. 4 shows x² from the normal modes (CGNM,noss) calculated using Eq. (5). As summarized in Table 3, the magnitude of x² from the REACH normal modes is slightly larger than from the atomistic MD but is in satisfactory agreement given the simplicity of the REACH method. The value of x² for residue 95 in DHFR is much larger from the atomistic MD than from the REACH calculation because of the presence of the cis peptide bond between residues 95 and 96.

Display large version of this figure

Figure 4

Mean-square fluctuation, x², of each residue, derived from atomistic MD (dashed curve) and REACH normal modes with (CGNM,ss, solid curve) and without (CGNM,noss, dotted curve) secondary-structure force constants: (1) myoglobin, (2) DHFR, and (3) plastocyanin.

Normal-mode analyses were also performed using a modified potential function in which the force constant parameters were derived separately for the secondary-structural elements (CGNM,ss). The residue dependences of x² with and without the secondary-structural force constants are slightly different (see Fig. 4), as are the magnitudes of x² and the correlation coefficients with the atomistic MD results (see Table 3). The agreement in x² magnitude with the atomistic MD is slightly better when the force constant model parameters defined separately for secondary structural elements. However, this may only be because the additional parameters enable better fitting to the k(r) distributions.

In Fig. 5, x² from the CGMD is shown for the three proteins. The averages and correlation coefficients with the atomistic MD are listed in Table 3. The average x² from the CGMD is smaller and in better agreement with the atomistic MD than that obtained from the normal modes. The finite timescale of 1ns may decrease protein fluctuations compared with the normal-mode results, which cover an infinite timescale, given the harmonic approximation to the potential energy. Again, the large x² from the atomistic MD at residue 95 in DHFR, probably because of the cis peptide bond between residues 95 and 96, is not reproduced by the CGMD even though in this case a smaller k=491 [kJ/mol·Å²] was used for the residue pair instead of the 1−2 interaction force constant. Additional parameters for the 1−3 and 1−4 interactions may be necessary for future modeling of cis peptide bonds.

Display large version of this figure

Figure 5

Mean-square fluctuation, x², of each residue derived from atomistic MD (dashed curve) and REACH CGMD with secondary-structure force constants (CGMD, solid curve). x² from REACH CGMD using REACH force constants averaged over the three proteins (see text in detail) is plotted as a dot profile (CGMD,ave): (1) myoglobin, (2) DHFR, and (3) plastocyanin.

Finally, we investigated whether a single set of REACH force constants can be used in CGMD of the three proteins, in an analogous way that atomistic MD simulations of the three proteins are performed with a single atomistic force field (in our study, CHARMM 22 31). To do this, a new set of force constant parameters was constructed by averaging over the three proteins. For the local interactions, the derived “constant” model parameters were averaged. For the nonlocal interactions, the three distributions of k_nonlocal(r) were superimposed and fitted with a new model function: the resulting function derived is, k_nonlocal(r)=4810 exp(−0.872r)+1.7 exp(−0.068r) [kJ/mol·Å²]. The new, common force-constant parameters were then applied to perform CGMD (CGMDave) of the three proteins, and x² was calculated.

The results, shown in Fig. 5 and Table 3, indicate that the adoption of the average set of force constant parameters increases the average fluctuations for the three proteins. However, the magnitude of the average x² remains similar to the atomistic MD, and the residual x² is still highly correlated with the atomistic MD.

Protein internal dynamics simulated by the residue-scale REACH CGMD is now compared with that obtained from the atomistic MD. To do this, the vibrational density of states, g(ω), and the dynamic cross-correlation matrix elements, c_ij, were calculated from the MD trajectories of the three proteins.

Fig. 6 shows g(ω) for the three proteins. In the low-frequency region of interest (ω<200cm⁻¹) g(ω) calculated from the atomistic MD by averaging over all the atoms is found to be similar to that averaging over only the C_α atoms (results not shown), a result following from the global, collective nature of low-frequency protein motion. The g(ω) from the CGMD lacks intensity at ω=100−150cm⁻¹ relative to the all-atom model; i.e., there is an absence of vibrational modes at high frequencies, as was also found in previous work demonstrating that the residue-scale CG force field cannot reproduce well >100cm⁻¹ all-atom protein dynamics 23. In contrast, the g(ω) intensity at ω<100cm⁻¹ is larger in the CGMD than in the atomistic MD spectrum, indicating that additional lower-frequency motions are added in the CGMD to satisfy the amplitudes present in the atomistic MD.

Display large version of this figure

Figure 6

Vibrational density of states, g(ω), derived from atomistic MD (dashed curve) and REACH CGMD (solid curve): (1) myoglobin, (2) DHFR, and (3) plastocyanin.

Fig. 7 shows the cross correlations c_ij for the three proteins. The agreement in c_ij between the atomistic and CGMD is good for the three proteins. Particularly good agreement is seen for the correlations of C_α atom pairs with large positive values. Thus, the correlations of C_α atom pairs in atomistic MD simulations are well reproduced by the REACH CGMD. However, the magnitude of cross correlation from the CGMD is in general smaller than that from the atomistic MD. One reason for this may be that the cross-correlation matrix is related to the contact map, defined by the pairwise distance or the interaction strength. Atom pairs with larger magnitudes of c_ij are likely to be close in distance and to interact strongly. Coarse-graining eliminates explicit covalent interactions between C_α atom pairs and thus decreases the dynamic correlations.

Display large version of this figure

Figure 7

Dynamic cross-correlation matrix, c_ij, derived from atomistic MD (upper row) and REACH CGMD (lower row). In lower right (upper left), c_ij is plotted when c_ij is positive (negative). Pairs with c_ij>0.3 are plotted with dots using the same gray shade as c_ij=0.3: (1) myoglobin, (2) DHFR, and (3) plastocyanin.

This work examines the transferability of the REACH CG simulation potential among proteins of different structural classes: myoglobin (α-fold), DHFR (α/β-fold), and plastocyanin (β-fold).

The REACH force constants for the three proteins were obtained from the variance-covariance matrices calculated from the atomistic MD trajectories. The force constants, when plotted as a function of the residue-residue pairwise distance, exhibit distinct distributions for the local (the virtual 1−2, 1−3, 1−4, and 1−5) and nonlocal interactions. These distributions were found to be similar among the three proteins. Model functions of distance dependence of the force constants, constructed for each type of interaction, were again found to be similar for the three proteins. Both the REACH CG normal modes and MD simulations reproduce reasonably well the mean-square fluctuations from the atomistic MD. The introduction of separate model function parameters for describing interaction within α- and β-secondary-structural elements improves slightly the agreement in x² values. A further introduction of interaction model function parameters, e.g., including residue-type dependence, may decrease the width of the distributions in the distance-dependent force constants, increase similarity of the distributions among the three proteins, and also improve the x² agreement.

Force-constant model functions were constructed by averaging the results over the three proteins, with the goal of examining whether a single analytical residue-scale force field, transferable between structural classes, could be generally applied, analogously to atomistic force fields such as CHARMM 1, AMBER 2, and GROMOS 3. This force field averaged over the three proteins is also found to reproduce reasonably well the mean-square fluctuations for all the three proteins, indicating that a single force field for the residue-scale CG model can be applied to different protein classes. Therefore, the REACH force constants from the three proteins can be, for practical purposes, considered identical. One implication of this finding is that the REACH method indeed maps directly the atomistic MD force field onto a generic residue-scale CG potential energy function.

To examine whether not only fluctuations but also timescales and residue-pair correlations found in atomistic MD protein internal dynamics are also reproduced by the CGMD, the vibrational densities of states and the dynamic cross-correlation matrices were calculated and compared. Coarse graining decreases the spectral amplitude of higher-frequency motions (> ∼ 100cm⁻¹) as a result of the reduction in the number of degrees of freedom, and the magnitude of the C_α atom correlations is also reduced somewhat. Nevertheless, the correspondence between the atomistic and CGMD spectral and correlation results is satisfactory given the simplicity of the REACH method.

The REACH method is self-consistent, involving a direct mapping of atomistic MD results onto the CG model. This self-consistency is clearly evident in the finding of this study that the distributions of distance-dependent force constants and the associated model functions are similar for the three proteins and that the CG force field constants averaged over those of the three proteins lead to a generalized set of force fields that reproduce the mean-square fluctuations and vibrational frequencies from the atomistic MD for all three proteins.

In future work, the REACH force constants will be calculated from atomistic MD simulation of a larger number of proteins so as to obtain better statistics for generalized residue-scale force-constant functions. Further validation of the generic force field will be performed using an extended subset of proteins that are not used for derivation of the REACH model parameters. Furthermore, a more flexible CG potential energy function, including more extensive anharmonic effects, will be implemented in the REACH methodology, allowing the dynamics present in MD at physiological temperatures to be reproduced. As a future extension, an anharmonic potential for C_α atom pairs between secondary structural elements may be introduced while keeping the intrasecondary-structure interactions harmonic. Subsequent development may include the incorporation of anharmonicity into intrasecondary-structure interactions, allowing formation and deformation of secondary-structural elements to be simulated. Dynamic solvent effects will be taken into account, and a methodology for calculating the Langevin friction on the C_α atoms from atomistic MD trajectories will be implemented 36,37.