Replica exchange molecular dynamics (REMD) is a kind of enhanced sampling algorithm that has been widely used in the study of protein conformation changes with their free energy landscape. Due to its rigorous principle and high sampling efficiency, researches on the optimization and development of the conventional REMD method have also sprung up, and have been in the ascendant recently, which greatly enhances its sampling efficiency and expands its application. However, those REMD variations usually have specific application, making it difficult to choose suitable REMD variations in practical application. Therefore, it is of great significance to summarize the REMD variations for understanding their advantages and disadvantages and for choosing suitable REMD variations. Here, we review recent development of the REMD variations from the perspective of their principle, and classify them into six types of methods: I) improving the approach and swapping rate of replica exchange, II) reducing the potential energy on exchange attempts to reduce the number of replica, III) Hamiltonians replica exchange molecular dynamics, IV) adjusting simulation process to improve sampling efficiency, V) changing the sampling methods, VI) adapting to heterogeneous and distributed computing environment. We hope that this review may be helpful in the understanding, application, and further improvement of the various REMD variations.
Contents
1 Introduction
2 Temperature replica exchange molecular dynamics
3 Improve the approach and swapping rate of replica exchange
4 reduce the potential energy on exchange attempts to Reduce the number of replica
Protein is one of the key biological macromolecules in life activities, and its function is often accompanied by conformational changes[1,2]. For example, protein kinases, epidermal growth factor receptor and tyrosine kinases, which are important drug targets, have distinct conformations of activation loops in activated and inactivated States[3]. In addition, the specific conformation of targeted proteins is also a highly selective source of some drug molecules[4⇓⇓⇓~8]. Therefore, the study of protein conformational changes is essential for drug discovery and design. However, the existing experimental techniques of structural biology, including X-ray crystallography and cryo-electron microscopy, are difficult to capture the dynamic behavior of proteins. Although NMR can provide dynamic information of protein in aqueous solution, it is difficult to be widely used due to the size of the research system. Among the computer simulation methods, Molecular Dynamics (MD) simulation, which is based on the principle of Newtonian mechanics, can observe the dynamic changes of proteins at the atomic level, thus becoming a general method to study the conformational changes of proteins[2,9]. In recent years, with the development of MD algorithm and computer hardware, the MD simulation scale can reach microseconds or even tens of microseconds, but it is still smaller than the actual time scale of important protein conformational transitions (generally milliseconds to seconds). Moreover, the conventional MD simulation method has low sampling efficiency and is easy to fall into local minima, which greatly limits its wide application in the study of large-scale protein conformational changes.
In order to further improve the sampling efficiency of MD, the development of enhanced sampling algorithms is in the ascendant, and various enhanced sampling algorithms are flourishing. These enhanced sampling algorithms can be classified into two categories according to whether they need to impose a limiting potential on the Collective Variables (CVs, also known as reaction coordinates): reaction-coordinate-dependent restricted enhanced sampling algorithms and reaction-coordinate-independent unrestricted enhanced sampling algorithms. The restricted enhanced Sampling algorithm is to apply additional potential energy to the simulation system on the predetermined CVs, such as atomic distance, dihedral angle, root mean square deviation, etc., so as to enhance the sampling in the CVs space. Typical methods include Umbrella Sampling (US), Steered MD, Targeted MD, Metadynamics (MTD), etc. This type of approach usually requires an in-depth knowledge of the study system to meet the pre-set CVs. Moreover, there are no suitable general CVs for complex protein conformational changes and the binding and dissociation of ligand molecules. In addition, the additional potential energy applied will limit the sampling to the preset CVs space, so that the obtained reaction path is limited. Moreover, for the simulation system with rugged potential energy surface, the reaction energy barrier may not exist in the selected CVs space, but in other CVs space, which is prone to the problem of hidden energy barrier. Therefore, whether the reaction coordinate is appropriate or not directly determines the sampling effect of this kind of method, which makes it difficult to apply to complex protein conformational changes. The unconstrained enhanced sampling algorithm, which does not depend on the reaction coordinate, uses some generalized physical quantities (such as temperature and system potential energy) as the enhanced sampling parameters, and can avoid the problem of CVs selection. Typical methods include simulated tempering and Replica Exchange molecular dynamics (Replica Exchange MD, REMD), among which REMD has been widely used in the conformational sampling of biological macromolecules, such as the conformational changes of proteins such as adenylate kinase and GPCR, since it was first applied to the folding simulation of methionine enkephalin chain in 1999[10⇓~12].
Temperature-based REMD is also known as Parallel Tempering, that is, N replicas are simulated in Parallel under the gradient temperature Tn(n=1,2,…,N), and the simulation temperature of each replica is exchanged regularly, so that each replica can be repeatedly traversed between the highest temperature and the lowest temperature to achieve the effect of enhanced sampling. In the dynamic simulation process, the system energy E is the sum of the system potential energy Epot depending on the atomic coordinate (X) and the system kinetic energy Ekin depending on the atomic velocity (V) (Formula 1).
According to the Metropolis criterion, the exchange probability (Pab) between neighboring temperature replicas (Ta and Tb) is calculated periodically.
Where β=1/kBT(kB is the Boltzmann constant) is the temperature factor, and . In order to ensure that the atomic velocity still obeys the Maxwell distribution after the replica exchange, the kinetic energy should be adjusted by using the temperature adjustment factor (Equation 3).
Where is the kinetic energy of the exchange of replica a to replica B and, correspondingly, is the kinetic energy of the exchange of replica B to replica a. After adjustment, the kinetic energy difference before and after the exchange of adjacent replicas is zero, and the exchange probability is only determined by the potential energy of the system.
In Equation 4, the exchange probability Pab is closely related to the temperature of the replica, and the smaller the temperature difference. DELTA.T between the replicas is, the greater the exchange success rate is.
The REMD method is often used to sample the conformation of small organic molecules and biomacromolecules, and even to study the structural transitions of Lennard-Jones clusters[13]. A large number of studies have shown that REMD method is an ideal tool for exploring the conformational distribution of organic molecules[13⇓~15]. Our REMD simulation on the macrocyclic compound loratinib also found that this method could accurately predict the conformational distribution of loratinib in different polar environments, and was consistent with the results of NMR experiments[14]. The application of REMD method in biological macromolecular systems such as proteins is also common[16⇓⇓~19]. For example, Pitera et al. successfully simulated the folding process of Trp-cage protein by REMD method[20]. In addition, the application of REMD in large protein systems, such as the prediction of the conformational transition energy barrier of the spike protein of COVID-19, also shows that REMD has the potential to study large-scale conformational changes[21]. REMD is also found in the simulation of complex chemical reactions. For example, Xin Liang et al. tested the simulation robustness of the REMD method on a series of reaction model systems with different activation energies based on the ReaxFF force field[22]. Due to the development of REMD and its integration with other methods in recent years, REMD can also be used to study specific problems such as biomineralization and binding free energy calculation[23][24].
In practice, the lowest temperature (Tmin) of the replica is generally set to the temperature of interest, such as room temperature 300 K, while the highest temperature Tmax is the lowest temperature that can help to cross the energy barrier of the system. In order to obtain a better sampling effect, it is necessary to ensure that the total average exchange rate between replicas is in the range of 0.2 to 0.5, and that each simulation replica can repeatedly traverse all temperatures. Therefore, it is necessary to ensure that the ratio of the potential energy difference ΔE to the energy fluctuation δE of all replica pairs is approximately the same. The energy fluctuation δE of the system at temperature T is proportional to the square root T of the total degree of freedom f of the system, while the potential energy difference ΔE is proportional to fΔT. Therefore, ΔE/δE is proportional to the /T, that is, the replica pair ΔT is proportional to the . When the total degree of freedom of the system increases, ΔT decreases, and for REMD simulations covering a certain temperature range (Tmin and Tmax are determined), the number of replicas increases accordingly. Therefore, for biological macromolecule systems with large degrees of freedom, especially in explicit solvent model systems, ordinary computing platforms are difficult to bear the huge computational cost caused by a large number of copies. Therefore, improving REMD, reducing its computational resources and further improving its sampling efficiency have become one of the important directions of REMD method development.
In recent years, the optimization of REMD methods has mushroomed. Review articles related to REMD are also frequently reported. For example, Ostermeir et al. 2013 focused on the application of REMD in protein and polypeptide flexibility simulation research, and introduced some REMD deformation methods to reduce the number of copies, such as hybride REMD, H-REMD (Hamiltonians REMD), etc[25]. Bernardi et al mentioned some mainstream REMD deformation methods, such as H-REMD and M-REMD (Multiplexed REMD), when summarizing the enhanced sampling algorithm in 2015[26]. Sugita et al introduced the Replica Exchange with Solute Tempering (REST) method in detail in Chapter 7 of Biomolecular Simulations, and listed other commonly used REMD deformation methods[27]. However, with the rapid development of REMD optimization, many new REMD deformation methods have emerged in the study of protein conformational changes in recent years, such as multi-dimensional REMD integrated with the latest sampling technology. In order to further understand the principle and application scope of various REMD deformation methods and make better use of these new methods in practical applications, it is necessary to summarize the various REMD deformation methods, especially the new trends of recent development, so as to help the application expansion and continuous improvement of REMD.
3 Improve replica exchange mode and exchange rate
The conventional REMD adopts the neighbor replica exchange mode, and only the even replica pair is considered in the even exchange, and only the odd replica pair is considered in the odd exchange. For example, for a REMD simulation of 2 N parallel copies, odd swaps compute whether the copy pairs 1 and 2, 3 and 4, …, 2 N − 1 and 2 N are swapped, and even swaps compute whether the copy pairs 2 and 3, 4 and 5, …, 2 N − 2 and 2 N − 1 are swapped. The exchange probability is obtained by the Metropolis algorithm such that the replica exchange satisfies the Detailed Balance Condition and the Markov Chain Monte Carlo state transition Condition. The successfully exchanged replica pairs facilitate the crossing of the potential barrier by exchanging atomic coordinates or simulation conditions. Therefore, improving the success rate of replica exchange can improve the sampling efficiency of each replica, which is also an important direction of the development of REMD methods. For example, Kouza et al. Used speed regulation to achieve zero rejection rate of replica exchange, but the study found that this forced exchange strategy did not significantly improve the efficiency of REMD sampling[28]. In addition, the Suwa-Todo algorithm proposed by Suwa and Todo introduces a graphical method called Weight Allocation to solve the detailed balance condition, which can minimize the exchange rejection rate under the global balance condition. However, due to the limited number of exchanges, the advantages of Suwa-Todo algorithm are not fully reflected in the neighborhood replica exchange mode of conventional REMD. In the Replica-Permutation Method (RPM), which is not limited to neighbor swapping, the advantages of Suwa-Todo algorithm are reflected[29,30]. In RPM, swapping can be performed between more than two replicas, that is, for a simulation with N replicas, there are N! -1 possibilities for the replica state of RPM after swapping. In recent years, the RPM strategy has also been gradually adopted by other REMD deformation methods[31]. However, when the number of replicas increases, the time cost of RPM method in the calculation of replica exchange probability will increase significantly. In order to reduce the number of combinations of Replica pairs without loss of sampling efficiency, Okumura et al. Further proposed the Replica Sub-Permutation Method (RSPM)[32]. In RSPM, the replica exchange only considers the exchange with two adjacent replicas at the same time, which significantly reduces the number of RPM replica combinations and saves the computational cost of RPM exchange probability.
On the other hand, only increasing the replica exchange frequency can also improve the sampling efficiency[33]. For example, in the Gibbs Sampling Replica Exchange method, all replicas can be exchanged in pairs with each other[34,35]; Brenner et al. Improved the Metropolis exchange probability algorithm based on Calvo's full replica pair exchange strategy (Calvo's all-Pairs exchange), thereby enhancing the sampling efficiency[36]; Multiscale Implementation of Infinite-Swap REMD (REMD-MSSA) combines with heterogeneous Multiscale methods to speed up the replica Swap frequency and reach the so-called Infinite swap limit, thus increasing the sampling efficiency[37]. Another method of Infinite swap REMD (Infinite Swapping REMD) is to increase the frequency of replica Swapping by using the limit equation obtained by infinite swap limit analysis[38,39]. However, the calculation of the coefficients involved in the limit equation needs to cover all possible permutations and combinations of the copies, which increases the computational cost accordingly, and sometimes the benefit is very small, and it may also cause a higher exchange rejection rate, so the calculation of the limit equation needs some approximation in practical application. The Designed-Walk replica exchange method (DEWREM), developed by Urano et al., implements a new replica exchange method by performing the exchange in even-odd order, that is, the replicas follow a pre-planned path in the temperature space instead of a random walk[40]. In this way, the walk between high and low temperatures is forced to be periodic, that is, the exchange of all odd replica pairs (1 and 2, …, 2N-1 and 2N) is carried out first. When one replica pair is successfully exchanged, the two replicas stop exchanging until all odd replica pairs complete one exchange, and then all even replica pairs are exchanged in the same way, and the cycle is repeated. It is worth noting that due to the large difference of the potential energy of each replica at different temperatures, although the replica exchange rate can be properly improved by various algorithms to improve the exchange mode, the exchange efficiency will not be significantly improved.
4 Reduce the total potential energy of the system to reduce the number of copies when performing the exchange
Since the number of copies is proportional to the square root of the total degrees of freedom of the system, the simulation of large proteins such as transmembrane GPCRs in explicit solvents requires a large number of copies. Therefore, reducing the number of replicas is the urgent need and improvement goal of REMD optimization.
4.1 Treatment of the degree of freedom of solvent in the system
Solvent atoms often account for more than 60% of the simulation system of biological macromolecules, which greatly increases the total degree of freedom of the simulation system and leads to a sharp increase in the number of copies. Therefore, how to deal with explicit solvents to reduce the number of replicas has become one of the important directions of REMD optimization. Compared with the explicit solvent, the implicit solvent models (including PB and GB models) are simple and fast, and can accurately reflect the solvation effect of the system, so they are widely used in the early molecular dynamics simulations, but the free energy surface obtained by the implicit solvent simulation is still different from that of the explicit solvent simulation. Therefore, how to use the implicit solvent in the REMD simulation while obtaining more accurate simulation results or how to use the explicit solvent while reducing the number of replicas is a problem that needs to be considered in the further optimization of REMD.
Because the explicit solvent leads to a significant increase in the total degree of freedom of the system, many REMD optimization methods adopt explicit-implicit hybrid simulation (hybrid REMD), that is, the explicit solvent is used to simulate the system, but the implicit solvent or mixed solvent model is used to calculate the energy during the exchange, thus reducing the number of simulation copies. Because the system is still simulated in the real explicit solvent, the simulation results are less affected by the implicit solvent model during exchange. For example, in the simulation of kernel REMD (kREMD), all conformations are sampled by explicit water molecules, while the energy obtained from the continuous solvation model (GB) is used as a scoring function to rank the conformations in the exchange process, and the kernel exchange probability is used instead of the Metropolis probability[41]. The present two different continuous solvation models (PB and GB models) can also be used in the hybrid REMD method. For example, Mu et al. Used PB implicit solvent model to calculate the solvation energy of solutes[42]. Chaudhury et al. Used the GBMV2 implicit solvent model to further reduce the calculation amount of the PB model[43]. In addition, in order to further reduce the number of replicas, the solvation energy can also be adjusted to further reduce the degree of freedom of the solvent, and the representative of such methods is REMshH (Replica-Exchange Molecular Dynamics Simulation with Scaled Hybrid Hamiltonians) method[44]. In REMshH, the solvation energy obtained by the GB calculation is further scaled down by an adjustment factor of 0.2 to reduce the number of replicas required for the simulation. This method can obtain accurate results in simple simulation systems, but the application of this method in more complex simulation systems is greatly limited because the implicit water model can not accurately consider the hydrogen bond interaction of water molecules and the implicit water model overestimates the electrostatic interaction. Therefore, a mixed solvent model has been developed based on this method, that is, when exchanging, a small number of solvent molecules closely related to the solute molecules are used in the explicit water solvent model, while the outer solvent molecules are used in the implicit solvent GB model, so that the hydrogen bonding and other special solvent effects of the explicit solvent are taken into account to improve the accuracy of the simulation[45].
The hybrid REMD method uses two different solvent models, that is, the explicit and implicit solvent models are used in the simulation and exchange respectively, but the system energies obtained by the two solvent models are not one-to-one correspondence before and after the replica exchange, so it is likely that the simulation system will obtain an incorrect ensemble. In order to overcome this defect, our group developed the REMD method of explicit-implicit solvent model based on velocity optimization, which correlates the difference between the energies obtained by the two solvent models by means of velocity adjustment (formulas 5 and 6, where Pis is the solvent-related energy calculated by the implicit solvent model), so as to obtain accurate simulation results[46,47]. In order to further remove the defects of the implicit solvent model, our group developed the vsREMD (Velocity-Scaling Optimized REMD) method by using the explicit water solvation energy to calculate the exchange probability on the basis of Velocity adjustment, so as to obtain more accurate simulation results[11].
4.2 Degree-of-freedom treatment of solute in the system
In the PREMD (Partial Replica Exchange Molecular Dynamics) method, the simulation system is divided into the Central Group and the Surrounding Group, so the Hamiltonian of the simulation system can be divided into three parts, namely, the internal energy of the Surrounding region Partial Replica Exchange Molecular Dynamics, the energy of the surrounding region-key region Esc, and the internal energy of the key region Ecc[48]. By coupling the temperature control system with two heat baths, PREMD only exchanges the simulation temperature for the key domain in the system, while the surrounding domain maintains the baseline temperature, so the exchange probability of PREMD can be approximately determined by (Esc+Ecc). LREMD (Local Replica Exchange Molecular Dynamics) further assumes that the surrounding domains share a set of atomic coordinates, so the Exchange probability of LREMD is completely determined by (Local Replica Exchange Molecular Dynamics)[80]. All in all, PREMD and LREMD are based on the assumption of weak coupling between the key domain and the surrounding domain, and the conformational change of the surrounding domain is small during the simulation. However, this assumption is usually not valid in many important biological processes, so the application of this kind of method in macromolecular biological simulation system is rare.
5 REMD Method Based on System Hamiltonian
The above REMD deformation methods all use temperature as a variable to exchange copies. Although the simulated temperature, as a physical quantity independent of the reaction coordinate of the system, is easy for users to operate, its enhanced sampling efficiency is not obvious for temperature-insensitive conformational changes. Therefore, REMD based on Hamiltonian energy (H-REMD) changes the temperature variable on the basis of REMD, and uses part of the Hamiltonian energy as a variable, such as changing the non-bonding interaction and pH of the system, which play a key role in protein folding and protein-small molecule binding[49].
5.1 Adjust the potential energy of the system to improve the sampling efficiency
The advantage of H-REMD is that the potential energy of the system is used to replace the temperature as the enhanced sampling variable[49]. Although the effect of adjusting the potential energy of the system can be regarded as changing the temperature, because the potential energy of the system can be decomposed into multiple energies, different energy components can be used as enhanced sampling variables to meet various simulation needs. For example, the potential energy of the system can be decomposed into three parts according to the molecular composition: the internal interaction energy of protein (Epp), the protein-solvent interaction energy (Epw) and the internal interaction energy of solvent (Eww).Adjusting the first two energies (with an s-adjustment factor for each replica corresponding to the effective temperature) is equivalent to controlling the simulated temperatures of Epp and Epw, while Eww maintains the baseline temperature, which effectively reduces the number of replicas required for REMD Equations 7 and 8.
The force field energy can be divided into a bonding term Ebonded and a non-bonding term Enon-bonded,In the general case, the bonding terms include the bond stretching term Ebonds, the bond angle bending term Eangles, the dihedral angle twisting term Edihedrals, and the impropers term Eimpropers.The non-bonding term includes the Coulomb interaction Eele and the van der Waals interaction EvdW (Equation 9).
van der Waals replica exchange (van der Waals REM, vWREM) controls the van der Waals action of the system by exchanging the van der Waals radius adjustment factor (Equation 10) and mimics the chaperone effect by reducing the hydrophobic interaction, thus achieving sufficient sampling[50]. The Coulomb Replica-Exchange Method (CREM) changes the proportion of the electrostatic term in the total energy of the system by exchanging the adjustment factor of atomic charge (Equation 11), so that the system can be fully sampled[51]. For example, Itoh et al. Successfully simulated the dimerization process of β-amyloid protein by Coulomb replica exchange method[52].
The H-REMD method can also utilize the bonding term energy as an enhanced sampling variable. For example, functional conformational changes in proteins primarily involve changes in torsional degrees of freedom, so MD in torsional space can sample the most relevant degrees of freedom. The Replica Exchange Torsion Angle MD (REX-TAMD) method is to enhance the sampling by exchanging the Torsion Angle term adjustment factor[53,54]. In addition, Okumura et al. Also introduced two parameters to control the skeleton dihedral angles of α helix and β sheet respectively, and changed the simulation parameters by replica exchange[55].
Q-REMD uses two parameter factors to nonlinearly adjust the empirical potential energy (Equation 12), where the value of the parameter Q starts from 1 and gradually increases with the increase of the replica temperature, so as to modify the potential energy function to reduce the potential barrier and achieve the effect of enhanced sampling[56].
The REST method proposed by Berne, based on the H-REMD principle, divides the total potential energy of the system into three parts: the internal interaction energy of solute (Epp), the interaction energy between solute and solvent (Epw) and the internal interaction energy of solvent (Eww)[57]. The commonly used REST is a modified version of REST2, which adjusts the solute-dependent potential energy according to Equation 13. It is equivalent to that the solute atoms of each copy are simulated at the set temperature gradient (Tm), while the remaining atoms are moving at the baseline temperature (T0). Therefore, the Eww is cancelled in the calculation of the exchange probability, thereby removing the influence of the explicit solvent on the number of copies (Equation 14), but in many cases, this treatment will cause the phenomenon of hot solute and cold solvent.
REFT (Replica Exchange with Flexible Tempering) divides proteins into rigid regions and Flexible regions on the basis of REST. The flexible region of proteins is set as the solute part, and the solvent and rigid part of proteins are set as the solvent part[58]. The gREST (Generalized REST) proposed by Kamiya et al. Further defines the solute as a part of the protein structure or a part of the potential energy of the system, such as the dihedral angle energy term or the van der Waals energy term[59]. If a certain part of the potential energy is taken as the solute, the total potential energy of the system can be expressed as Equation 15, where li is the maximum number of atoms involved in the i-th solvent-solute interaction. For example, I is 2 considering the bond energy, van der Waals or Coulomb interaction as the solute, I is 3 considering the bond angle, and I is the number of solute atoms participating in the i-th solvent-solute interaction 4;ki considering the dihedral angle, and I is 1 considering the van der Waals or Coulomb interaction as the solute. This method is suitable for the simulation of local motion in the critical region of the system.
FREM (Fragment Replica-Exchange Method) divides the protein into multiple fragments dispersed in different copies[60]. Therefore, for each copy, the total energy of the system can be divided into fragment energy and non-fragment energy (Equation 16).
The segment energy is scaled by an adjustment factor s, and the s-parameter of each replica is not uniform, so the swap probability calculation formula can be written as Equation 17.
MFREM (Multiple Fragment-size Replica-Exchange) is based on the FREM method, which considers the s-factor of temperature and fragment energy together as an enhanced sampling variable[60]. ResEx (Resolution Exchange) adds some coarse-grained model copies with significantly fewer atomic degrees of freedom to the regular REMD copies, and crosses the energy barrier through these coarse-grained models[61].
5.2 Replica exchange introduces other enhanced sampling methods
Because the replica-swapping strategy in REMD can enlarge the sampling space while maintaining a regular simulation ensemble, many other types of enhanced sampling algorithms and free energy calculation methods have also introduced the replica-swapping strategy. Based on the original idea of enhanced sampling, these new methods change the parameters or CVs through replica exchange. On the one hand, they further increase the sampling efficiency, on the other hand, they also reduce the difficulty of parameter setting and CVs selection.
MTD is a widely used restricted enhanced sampling algorithm, which imposes a time-dependent bias potential (V) on CVs to enhance sampling in the space of CVs. Bias-exchange MTD (BE-MTD) introduces MTD into replica exchange, which further improves the sampling efficiency of REMD (Equation 18)[62].
Essential Dynamics is an enhanced sampling in the Essential set mode space of key substructures. Kubitzki et al. introduced it into replica exchange, that is, in each replica, only the essential set mode of the key structure is coupled to the high temperature, and the rest of the system is simulated at the baseline temperature, thus effectively reducing the number of replicas, and the increase in the temperature difference between adjacent replicas also improves the sampling efficiency[63]. Sugita et al. Also significantly reduced the number of REMD copies by integrating REMD with the multicanonical algorithm[64,65][66,67].
Accelerated Molecular Dynamics (AMD) adds a non-negative boost potential on the potential energy surface of the system to help the system cross the energy barrier (Equation 19), thereby expanding the sampling space, where Accelerated Molecular Dynamics is a preset energy threshold, and α is an acceleration factor. Fajer et al. Also combined AMD with the replica exchange method, and introduced the replica exchange by changing its energy threshold or α parameter, thereby improving the sampling efficiency[68]. Compared with AMD, Gaussian Accelerated Molecular Dynamics (GaMD) adds a biased potential that follows a Gaussian distribution to improve the energy weight distribution (Equation 20), where K is the harmonic force constant, and Gaussian Accelerated Molecular Dynamics and Emin are the maximum and minimum values of the potential energy of the system, respectively. rex-GaMD is a combination of replica exchange method and GaMD, which can realize two different versions of enhanced sampling methods: force constant rex-GaMD and energy threshold rex-GaMD.Rex-GaMD has the following advantages: 1[69]. The copies of the force constant rex-GaMD have different harmonic force constants and the same energy threshold parameters, while the copies of the energy threshold rex-GaMD have different energy thresholds and the same harmonic force constants.
Oshima et al. further combined the replica exchange strategy, GaMD and US into the GaREUS enhanced sampling algorithm. When performing the exchange, the system potential also needs to add the bias potential (Equation 21) added by umbrella sampling and Gaussian acceleration, where I refers to replica I[70].
Umbrella sampling achieves sufficient sampling within a window by applying a harmonic potential within a predetermined window (Equation 22), where refers to a predetermined reaction coordinate. Twenty years ago, Sugita et al. Introduced the replica exchange strategy into the US and developed the REUS enhanced sampling algorithm[24].
In recent years, the replica exchange strategy has also been introduced into the free energy calculation method. Enveloping Distribution Sampling (EDS) is able to calculate the free energy barriers of multiple endpoint States in a single simulation, and its reference state R contains N (≥ 2) Hamiltonian endpoint States. The accuracy of energy barrier calculation can be further improved by introducing replica exchange (Equation 23), where s is the potential energy smoothing parameter and ER is the energy shift[71].
MSλD allows the simultaneous evaluation of free energy differences between multiple compounds in a single strip simulation, introducing the Biasing Potential Replica Exchange Multisite λ-Dynamics (BP-REXMSλD) method[72]. BP-REXMS λD perturbs the flexible group to calculate the relative binding free energy of the compound, and the system potential energy is as shown in Equation 24, where X and xi represent the surrounding atomic coordinates and the atomic coordinates of ligand I, respectively, Msites refers to the total number of substituent sites S, Ns refers to the total number of substituent I at a substitution site, V is the bias potential, and λ is the perturbation parameter.
5.3 Multidimensional REMD
Because of its simple principle and considerable sampling efficiency, REMD is often combined with other enhanced sampling methods to form multi-dimensional REMD to further improve the sampling efficiency of molecular simulation (Table 1). For example, RSE-MTD, which introduces the replica exchange strategy in the BE-MTD method, can further enhance the sampling efficiency of BE-MTD and REMD[73⇓⇓~76][77,78]. In order to further reduce the number of replicas and improve the sampling efficiency of the integrated method of REMD and MTD, our group developed the ossPTMetaD method by introducing REMD into a new MTD method (Well-tempered MTD) and eliminating the influence of the internal interaction energy of the solvent when the replicas are exchanged. Tested on protein model systems with different motion types, it is found that ossPTMetaD requires only 30% of the computational resources of T-REMD and is effective in simulating the challenging kinase DFG turnover motion[79]. In addition, the combination of REMD with AMD (w-REXAMD), REUS, H-REMD and other methods is also often used in the study of protein conformational changes[80⇓~82][83,84][85]. ABMD (Adaptively Biased Umbrella Sampling) allows the system to be fully sampled within a preset window by imposing an adaptive bias potential. The combination of ABMD and REMD makes it possible to simultaneously enhance sampling in the two-dimensional space of CVs and temperature, so its Metropolis algorithm for calculating the exchange probability also takes into account both the bias potential and the temperature variable (Equation 25)[86].
In addition, the replica exchange strategy can also include other variables such as temperature when introducing the free energy calculation method. For example, in the envelope distribution sampling mentioned above, after introducing multiple copies, setting different external pH for each copy can accelerate its conformational change[87,88]. Free Energy Perturbation (FEP) is one of the most rigorous and accurate Free Energy calculation methods, which is widely used to predict the binding affinity of protein-ligand. In the principle of protein-ligand binding free energy calculation, the FEP method is to set a series of windows λ (from 0 to 1) to complete the conversion between two ligands and calculate their relative binding free energy. The combination of FEP and REMD deformation method further improves the accuracy of free energy calculation of FEP, which has attracted much attention in recent years[89⇓~91].
It is worth mentioning that the Tq-REM method, which is a combination of REMD and Q-REMD, is different from other multi-dimensional REMDs in the setting of replica parameters[92]. Each copy of the traditional multi-dimensional REMD has a different multi-dimensional variable setting, and the temperature is also arranged in a gradient from low to high. In the Tq-REM method, Lee et al. Replace the high-temperature replicas with Q-REMD, that is, the temperature gradient range only covers 1 to N-M (M ≥ 1) replicas, and the Q of these replicas is the same, while the replicas N-M to N are set to the same temperature but different Q. In addition, the multidimensional REMD method needs to calculate the multidimensional variables at the same time, which increases the complexity and computation of the swap probability formula. The Multidimension al Replica Exchange Simulations proposed by Cruzeiro is different from others in the calculation of Exchange probability, because it is not the integration of multiple variables, but only considers one dimensional variable in each exchange, although it uses temperature, pH and redox potential (Multidimension al Replica Exchange Simulations) to enhance the sampling of Multidimensional variable space[93]. The exchange variable is calculated in turn, for example, using Equation 26 when exchanging the pH, where the is the number of protons in the copy; The redox potential is exchanged using Equation 27, where F refers to Faraday's constant, refers to the redox potential of replica a, and refers to the number of electrons of the replica; When the temperature is exchanged, the exchange probability of the conventional REMD is used to calculate Equation 4.
5.4 Introducing a physical quantity of a simulation process as an enhanced sampling variable
In addition to replacing or adding variables directly related to the system potential for enhanced sampling, many other physical quantities involved in the simulation process are also related to the sampling efficiency, which can be combined with the system temperature to form a multi-dimensional REMD enhanced sampling method, thus reducing the number of copies required. Moreover, these variables often do not affect the potential energy formula of the system, so there is no need to modify the conventional REMD source code extensively. For example, V-REMD adds solvent viscosity as a replica exchange variable in conventional REMD simulation methods[94]. The lower the Viscosity of the solvent, the faster the protein chains diffuse, thereby enhancing the sampling rate of the conformational space. Surface tension is a force that exists at the interface of vacuum/water, oil/water, film/water and other liquids. Positive surface tension increases the surface area of the lipid bilayer, accelerates the lateral diffusion of lipid molecules, and decreases the thickness of the lipid bilayer due to volume incompressibility. Surface-Tension REMD can effectively simulate the dynamics of biomembrane systems by introducing Surface tension as a replica exchange variable through the NPγT ensemble[95]. Because the explicit lipid bilayer can not maintain the structural uniformity at high temperature, it affects the accuracy of the membrane system treated by conventional REMD. MSREMD (Mass Scaling REMD) method adjusts the particle Mass according to the reference temperature of each replica. The higher the temperature is, the heavier the particle is, which improves the stability of the simulation. It also avoids the problem of speed adjustment after exchange and saves part of the calculation time of REMD[96]. Volume-Temperature REMD (VTREMD), which exchanges Temperature and density between replicas, is suitable for conformational change simulations where the free energy landscape has many local minima[97]. VTREMD still performs the exchange according to the Metropolis algorithm, but takes into account that the potential energy after the exchange is affected by the volume V of the simulation box (Equation 28). Where P is the pressure and A is the number of molecules in the box.
6 Adjust the simulation process to improve the sampling efficiency
6.1 Temperature control process
A series of copies of REMD are simulated at a temperature gradient, and the success rate of exchange depends on the energy distribution of the simulated system at different temperatures. However, all replicas of TIGER (Temperature Intervals with Global Energy Reassignment) are exchanged at the same baseline Temperature Temperature Intervals with Global Energy Reassignment, which can significantly reduce the number of replicas[98]. The way of TIGER realizing TES lies in its three cyclic steps of replica exchange: (1) each replica is heated from the baseline temperature (T0) to the designated temperature (Tm), the atomic momentum in the replica is adjusted again with the and the system is equilibrated.(2) molecular dynamics sampling at a constant temperature Tm, (3) quenching of the system by energy optimization to a baseline temperature T0, and a short time of equilibration before exchange is performed (fig. 1). This scheme is based on certain assumptions to satisfy the detailed equilibrium condition, that is, the hydrazine conformation of the system is unchanged during the quenching process (such as dihedral angle) and the heat energy consumed during the quenching process is independent of the state of the system, but these assumptions are not valid in complex systems, such as peptide folding systems simulated by implicit solvents. TIGER2 improves step 3 by adjusting the momentum of the atoms in the replica to reach the T0 temperature by , thus satisfying the weak equilibrium condition[99]. In addition, in recent years, TIGER2 has successively solved the problem that the sampling method is not sensitive to the energy of the solute due to the large fluctuation of the solvent energy (TIGER2A), and further used an explicit-implicit mixed solvent to reduce the number of replicas (TIGER2h) and added an explicit solvent shell to consider the steric effect, solvent polarization and reaction field energy (TAGER2hs)[100][101][102]. TIGER3 is to reduce the van der Waals radius of atoms in high temperature simulation, so that protein chains can pass through each other at higher temperature, and to restore the atomic radius to normal value and rebalance before quenching, so as to eliminate the chain entanglement problem[103].
The effect of enhanced sampling can also be achieved by readjusting the replica temperature according to the replica exchange. For example, the Adaptive Temperature-Based Replica Exchange (T-REMD) enables the simulation to achieve the maximum average Exchange rate (EE-REMD) and the maximum Temperature range (CM-REMD) by resetting the temperature during the simulation[104][105]. In addition, Statistical Temperature Molecular Dynamics (STMD) uses an adaptive weighting factor to expand the simulation Temperature range for sampling, and combines with REMD to make each replica undergo STMD, which can significantly reduce the number of replicas[106].
Modified REMD 1 divides the replica into an active replica and an inactive replica (fig. 1), where the active replica is identical to the conventional REMD replica, but the simulation of the inactive replica is frozen to preserve the behavior of the simulation system at that temperature[107]. The addition of inactive copies can enlarge the temperature setting range and reduce the temperature setting interval, but it will also increase the computational cost. Modified REMD 2 records the simulated potential energy of each active temperature replica into a potential energy database, and randomly selects a potential energy from the corresponding temperature database to give to the inactive replica when performing the exchange, thereby avoiding the simulation of the inactive replica and saving the computational cost[108].
6.2 Parallel copy condition settings
In the conventional REMD simulation, the parallel copies have the same settings such as temperature and pressure control methods, except for the different simulation temperatures. To further improve the sampling efficiency of REMD and help it cross the energy barrier, Reservoir REMD sets the highest temperature replica as a pre-generated conformation library, while the other replicas operate normally[109⇓~111]. Therefore, the advantage of Reservoir REMD is that the exchange between the highest temperature replica and the adjacent replica can introduce a conformation with a longer spatial distribution distance to the adjacent replica, thus better helping it to cross the energy barrier and expand the sampling range. REDS (Replica Exchange with Dynamical Scaling) is to set some copies in all copies, which can dynamically adjust the potential energy in the system when executing the Exchange, so as to improve the success rate of the exchange[112]. Convective REMD increases the sampling range by setting one of all replicas to ramp up or down in the simulation. In addition, in REMD simulations, constraints imposed on the backbone hold a protein conformation of interest in a basin of space, and the constraints along the temperature gradient shrink from full strength to zero strength, so that the highest temperature becomes a conformationally non-specific reference state[113].
6.3 Parallel replica multilayer simulation
Multiscale MD (MuSTAR MD) sets each parallel replica to at least three layers, including an all-atomic fine-grained model, at least one coarse-grained model, and an ensemble variable model that interacts with other model layers in the same replica through a coupling energy term[114]. The coarse-grained model is used to facilitate the efficient sampling of large-scale conformational changes, and the fine-grained model is used to accurately sample the conformational changes. The collective variable model can not only be used as an intermediary between the coarse-grained model and the fine-grained model, but also can improve the sampling efficiency through high temperature. In Multiplexed REMD, each temperature replica layer also runs multiple independent parallel simulations, and exchanges are performed within the replica layer and between layers[115]. Based on the idea of Multiplexed REMD, mQ-REMD (Multiplexed Q-REMD) improves the conventional MD simulation of the replica layer to Q-REMD simulation mode, which further improves the sampling efficiency and reduces the number of replicas[116].
6.4 The change path expands the sampling range for the initial conformation
Chen et al. Selected the conformation from the HIV-1 protease change path obtained by NUMD as the initial conformation of the REMD copy, on the one hand, to help REMD cross the energy barrier and improve its sampling efficiency, on the other hand, to avoid the "path trap" problem of guessing the initial conformation[117][118]. Fujisaki et al. Also used the Onsager-Machlup equation of motion for path sampling as the initial conformation of REMD to fully sample the changing paths[119].
7 Change the sampling method
Different sampling methods have their own merits. For example, the Gibbs sampling replica exchange method mentioned above inherits the advantages of Gibbs sampling and can realize full replica exchange. Tsallis-REMD (Combination of Tsallis Weight Sampling and REMD) Sampling follows Tsallis Weight, which is different from the conventional MD sampling following Gibbs-Boltzmann weight, and can accelerate the convergence of sampling[120]. RXSGLD (Replica Exchanging Self-guided Langevin Dynamics) uses the advantages of Self-guided Langevin Dynamics, does not need high temperature to accelerate conformational movement, and its conformational distribution perturbation is not as large as that of high temperature, so the exchange rate is higher[121,122]. QM/MM-REMD (Quantum Mechanical/Molecular Mechanics-REMD) can be used to study the complex chemical processes of biological macromolecular systems in solution by taking advantage of QM/MM and using QM precision in key areas[123]. Sampling methods such as BOMD (Born-Oppenheimer MD), SE (Semiempirical) -BOREMD and CPMD (Car-Parrinello MD) introduce replica exchange to make them suitable for chemical reaction simulation[124][125][126]. For large-scale conformational changes, Monte Carlo replica exchange and ICMD (Internal Coordinate MD) help the system to cross higher barriers[127,128][129]. Kamberaj modifies the particle momentum expression in Newtonian mechanics by introducing the information of the global optimal conformation of all replicas, so that the replicas can communicate with each other individually or globally, and further optimizes the Newtonian equation of motion to obey the Maxwell-Boltzmann distribution[130][131].
8 Adapt to Heterogeneous and Distributed Computing Environment and Solve the Problem of Resource Demand
Grid computing network composed of heterogeneous and distributed computing environment and homogeneous high-performance cluster can solve the needs of large-scale REMD simulation. greedy REMD (gREMD) uses a pre-computed inter-replica exchange table to decentralize replica management, making it suitable for heterogeneous computing grids[132]. Platania et al. Proposed a Hadoop-based Replica Exchange (HaRE). Hadoop is a distributed system infrastructure developed by the Apache Foundation, which can run efficiently on homogeneous or heterogeneous distributed network infrastructures (DCIs) of different sizes and configurations[133]. The basic idea of asynchronous Replica swapping (Replica Exchange Molecular Dynamics (REMD) Simulations Asynchronously, ASyncRE) is to assign all replicas to a running or waiting list and allow a subset of replicas on the waiting list to perform the swap independently of other replicas on the running list[134]. Because replica swapping in REMD does not rely on a centralized synchronization step, the ASyncRE algorithm can be scaled to any number of processors and avoids the maintenance requirement of a static processor pool.
9 Conclusion and prospect
The enhanced sampling algorithm improves the efficiency of MD sampling through various principles, and has become one of the mainstream technologies to simulate the functional conformational changes of proteins. REMD is a kind of MD sampling method in which information is exchanged through multiple simulated replicas, which conforms to strict physical principles and can efficiently cross energy barriers to improve sampling efficiency, so it is widely used in the study of important protein conformational changes, and its morphology is constantly improved. Starting from the basic principles of REMD, this review summarizes the various new deformation methods of REMD that have emerged in recent years (Table 2). Since the number of copies of REMD is proportional to the square root of the total degrees of freedom of the system, the main improvement direction of REMD is to reduce the number of copies. Dealing with the explicit solvent and adjusting the potential energy of the system are the main improvements to reduce the number of REMD replicas in recent years. In addition, the advantages of REMD also attract the attention of other enhanced sampling algorithms. On the one hand, the combination of REMD and other enhanced sampling algorithms further improves the sampling efficiency of REMD and expands the application scope of REMD. On the other hand, it reduces the difficulty of parameter setting or set variable selection of other enhanced sampling algorithms. In addition, for the further improvement of REMD method, it is necessary to correctly handle the system state changes occurring in the large-scale conformational changes of proteins, including protonation state, disulfide bond cleavage and formation, while considering the reduction of the number of copies and the improvement of sampling efficiency. In a word, although a large number of new REMD methods have enriched the choice of conformational change simulation, there is no perfect method that is "universally applicable" to the actual simulation, and it is still necessary to take into account the computational cost and sampling efficiency to select the method suitable for the research system.
PDF
