Introduction
Proteins are the main biological macromolecules that carry out biological functions and are also the major components for building synthetic biological systems. The function of proteins depended on their specific three-dimensional spatial structure and specific intermolecular interactions (Anfinsen et al., Reference Anfinsen, Haber, Sela and White1961). The amino acid sequence determined the three-dimensional structure and interactions of the mass, thus determining the function of the protein (Zhang and Wang, Reference Zhang and Wang2023). The amino acid sequences of natural proteins have been selected over time by evolution and adapted to the functional needs of the corresponding organism (Marsh and Teichmann, Reference Marsh and Teichmann2015). The function of natural proteins might not be adequate for other required purposes, and even sometimes the available natural proteins cannot be found (Zhu et al., Reference Zhu, Avakyan, Kakkis, Hoffnagle, Han, Li, Zhang, Choi, Na, Yu and Tezcan2021). Therefore, design proteins with the target properties and functions are a significant breakthrough in the history of biology. Protein folding mechanism studies focus on the thermodynamic state and dynamic process from sequence to structure. The aim of de novo design of protein is to obtain proteins with a target structure. The process from structure to sequence can be regarded as the reverse process of protein folding.
The concept of de novo protein design had a long history. Most of the early protein design successes were based on the design of subordinate secondary structures. In the 1980s, William DeGrado and his colleagues designed a tetrameric structure with a helical sequence connected by three loops. The overall design strategy was relatively simple, and the obtained protein was the α4 protein (Regan and DeGrado, Reference Regan and DeGrado1988). Following this work, Hetch MH and his team reported their successful design of the four-helix Felix, in which more factors than the α4 protein were considered in the design process, such as hydrophobic interactions inside the helix, helical hydrogen bonding, etc. (Hecht et al., Reference Hecht, Richardson, Richardson and Ogden1990; Patel et al., Reference Patel, Bradley, Jinadasa and Hecht2009). These works had great implications for the de novo design of protein and accumulated rich experience. Based on these works, scientists then challenged the more difficult de novo design proteins. In 1994, Blanco et al. reported a successful design of β-hairpin, a motif containing 9 amino acids, which was the first reported successful design of β-hairpin (Blanco et al., Reference Blanco, Rivas and Serrano1994). In 1997, Mayo et al. reported their first design of a small protein full sequence design 1 (FSD-1), using computer-aided technology methods (Dahiyat and Mayo, Reference Dahiyat and Mayo1997). Another important development was that David Baker and his team designed a ββαβαββ structure, named Top7, whose thermodynamic stability was high. The denaturation temperature reached 69 °C in a guanidine hydrochloride solution, and it still maintained structural stability at 90 °C without denaturant, and its stability was much higher than that of natural proteins of the same scale (Huang et al., Reference Huang, Feldmeier, Parmeggiani, Velasco, Höcker and Baker2016).
Significant progresses have been made in de novo protein design, and protein scientists were studying the properties of these particular proteins to take advantage of them and provide functions that nature proteins do not have. For example, Correia et al. demonstrated that computational protein design could produce small, thermally and conformationally stable protein scaffolds with a neutralization epitope from respiratory syncytial virus, accurately mimicking the viral epitope structure and induce potent neutralizing antibodies (Correia et al., Reference Correia, Bates, Loomis, Baneyx, Carrico, Jardine, Rupert, Correnti, Kalyuzhniy, Vittal, Connell, Stevens, Schroeter, Chen, Macpherson, Serra, Adachi, Holmes, Li, Klevit, Graham, Wyatt, Baker, Strong, Crowe, Johnson and Schief2014). Marcandalli et al. reported that custom-designed vaccines were achieved by designing protein self-assembled nanoparticles as a skeleton to fix and present viral glycoprotein antigen complexes (Marcandalli et al., Reference Marcandalli, Fiala, Ols, Perotti, de van der Schueren, Snijder, Hodge, Benhaim, Ravichandran, Carter, Sheffler, Brunner, Lawrenz, Dubois, Lanzavecchia, Sallusto, Lee, Veesler, Correnti, Stewart, Baker, Loré, Perez and King2019). Sesterhenn et al. established the TopoBuilder system to design proteins that could stabilize complex structural motifs from scratch. Through this system, they obtained proteins that could present three antigens simultaneously (Sesterhenn et al., Reference Sesterhenn, Yang, Bonet, Cramer, Wen, Wang, Chiang, Abriata, Kucharska, Castoro, Vollers, Galloux, Dheilly, Rosset, Corthésy, Georgeon, Villard, Richard, Descamps, Delgado, Oricchio, Rameix-Welti, Más, Ervin, Eléouët, Riffault, Bates, Julien, Li, Jardetzky, Krey and Correia2020). Wannier et al. obtained red fluorescent proteins that maintained light intensity and were not easy to aggregate through surface mutation design and main-chain dependent side-chain rotamer sampling optimization (Wannier et al., Reference Wannier, Moore, Mou and Mayo2015). As we all know, protein functions are closely related to the conformational space and folding process. The commonly used experimental techniques for studying the mechanism of protein folding had restrictions because protein folding typically ranged from microsecond to milliseconds (Dill and MacCallum, Reference Dill and MacCallum2012). Studying the dynamics of protein folding and conformational exploration on such a fast time scale was extremely challenging, especially for de novo-designed proteins whose properties might be entirely different from the natural proteins (Korendovych and DeGrado, Reference Korendovych and DeGrado2020). Therefore, molecular modeling and simulation were important complements to experimental characterization for understanding the structural and dynamic properties of de novo design proteins at an atomic level (Wang, Reference Wang2021).
Molecular dynamics (MD) simulation provided the most detailed description of the dynamics and structures of de novo-designed proteins (Gonzalez et al., Reference Gonzalez, Li and McCully2022; Leong et al., Reference Leong, Lim, Ismail and Choong2018; Wang et al., Reference Wang, Hu and Zhang2016). However, big challenges also existed in the MD simulation method and were similar to the experimental biophysical approaches. Due to the special nature of de novo design proteins, the folding time scale was longer than that of natural proteins. In this case, the sampling efficiency of classical all-atom simulation methods would be reduced, resulting in difficulties in the study of folding dynamics (Zhu et al., Reference Zhu, Dai, Liang, Zhang, Gai and Lai2013). The stability of designed proteins in the simulations is also an important factor. The choice and efficiency of the sampling method and force field are the prominent elements in the MD simulations. In this review, we will focus on the sampling method and force field, as well as the stability and dynamics of MD simulations of de novo-designed proteins.
MD simulations of de novo-designed proteins
Sampling methods
The classical methods of traditional MD simulations have been wildly used for sampling conformation methods of peptides and large proteins (Bernardi et al., Reference Bernardi, Melo and Schulten2015; Kamberaj, Reference Kamberaj2018; Lazim et al., Reference Lazim, Suh and Choi2020; Wang, Reference Wang2021). To enhance the sampling efficiency and observe the complicated protein activities, one might need enhanced sampling methods rather than the traditional method due to the rugged free energy surface, especially for de novo-designed proteins whose protein folding funnels were more challenging than those of natural protein (Dill and Chan, Reference Dill and Chan1997; Leopold et al., Reference Leopold, Montal and Onuchic1992). One of the most effective enhanced sampling methods used in the designed proteins was replica exchange molecular dynamics (REMD). REMD was a simulation method that combined the replica exchange method with the MD method. The basic idea was to construct a series of interacting replica systems, one for each temperature, covering a wide range of temperatures from low to high, and then performed a separate MD simulation for each replica at each temperature (Bernardi et al., Reference Bernardi, Melo and Schulten2015; Li et al., Reference Li, Yang, Chen, Zhang, Wei and Zhang2023). The Monte Carlo Metropolis rule controlled the exchange of replicas adjacent in temperature and ensured the detailed balance, and that the structure ensemble at temperature of interested could be generated (Wang, Reference Wang2021; Zhou, Reference Zhou2007). At adjacent temperatures, every two copies can be exchanged according to the Metropolis standard, and based on this, the REMD method can make the low-temperature conformation space escape the local potential energy lowest point (Chen et al., Reference Chen, Xiao and Huang2015; Sugita et al., Reference Sugita, Kamiya, Oshima and Re2019). In practical applications, the minimum temperature of the replica (Tmin) was generally set as the temperature of most interest, such as room temperature (300 K), and the maximum temperature (Tmax) was the lowest temperature that can help to cross the energy barrier of the system. In order to obtain a better sampling effect, it was necessary to ensure that the total average exchange rate between the replicas was in the range of 0.2–0.5 and ensure that each analog copy can traverse all temperatures (Rosta and Hummer, Reference Rosta and Hummer2009). The wide application of REMD has also played a significant role in the field of protein folding. Lazim et al. explored the mechanism directing the structure variation from α/4β-fold protein to 3α-fold protein after mutation by conducting REMD simulation on 42 replicas and provided insights into the folding and unfolding pathway of 3α-fold protein and α/4β-fold protein, respectively, paving the way toward the understanding of the ongoing conformational variation (Lazim et al., Reference Lazim, Mei and Zhang2012).
In recent years, there have been significant developments in the optimization of REMD methods, with many related articles being reported. For example, Bernardi et al. mentioned mainstream REMD deformation methods such as H-REMD and M-REMD in their 2015 summary of enhanced sampling algorithms (Ostermeir and Zacharias, Reference Ostermeir and Zacharias2013). Mateusz Marcisz et al. also docked 24- and 48-meric glycosaminoglycans using an all-atomic repulsive-scaling Hamiltonian replica exchange molecular dynamics (RS-REMD) method. This advanced methodology was based on replicas with scaled van der Waals radii of interacting molecules and has performed well for proteins complexed with oligomeric GAGs (Marcisz et al., Reference Marcisz, Maszota-Zieleniak, Huard and Samsonov2021). Danial Pourjafar-Dehkordi studied the flexibility of native and pS111-Rab1b in complex with GTP or GDP using extensive MD simulations and an advanced sampling method called dihedral angle-biasing potential replica exchange molecular dynamics (DIA-REMD), which promoted backbone and side chain dihedral transitions along a series of replica simulations in selected protein segments and through exchanges also improves sampling in an unbiased reference simulation (Pourjafar-Dehkordi and Zacharias, Reference Pourjafar-Dehkordi and Zacharias2021). Thus, advanced REMD methods have made significant contributions to the simulation sampling of protein foldings, especially for proteins with rough potential energy surfaces.
In order to improve the sampling efficiency in protein simulations, there are many other sampling methods in addition to REMD, such as metadynamics, enhanced sampling method like ossPT-MetaD, population particle MD and other efficient sampling methods.
Metadynamics was a class of methods that improved sampling efficiency by introducing additional biased potential energy or force acting on certain degrees of freedom. These degrees of freedom, often called Collective Variable (CV), were generally functions of points in state space that distinguish between two or more thermodynamic states that needed to be studied. The free energy landscape along the collective variable was a low-dimensional projection of the complete free energy landscape in the direction of the CV, containing all the thermodynamic information we cared about, and the solution of it could be accelerated by the metadynamics (Giberti et al., Reference Giberti, Salvalaglio and Parrinello2015).
The ossPT-MetaD, as an enhanced sampling method, performed enhanced sampling in both energy and geometric space. By optimizing the sampling efficiency under REMD framework and combining the normal vibration mode as a set variable, the number of replicas could be reduced while maintaining a high exchange rate and simulation convergence. This significantly reduced the degree of freedom of the system and the number of copies required, improved the simulation convergence, and demonstrated high efficiency and rationality in a variety of protein systems (Peng et al., Reference Peng, Wang, Shi, Xu and Zhu2021).
The swarm population particle molecular dynamics (SPMD) was an enhanced conformation sampling method for computer simulation of protein folding, particularly suitable for use in combination with replica exchange methods (REM). In the framework of REM, the sampling efficiency was enhanced by introducing the concept of swarm intelligence. In SPMD, individual copies were distinguished not only by temperature but also by particle swarm optimization for efficient search and exchange in the conformation space. In the framework of REM, the sampling efficiency is enhanced by introducing the concept of swarm intelligence. In summary, SPMD was an advanced computational method that significantly improved the efficiency of MD simulations through enhanced sampling techniques, especially in protein folding and conformation space exploration. This method provided a powerful tool for efficient exploration and prediction of protein structure and contributed to a deep understanding of protein dynamic behavior and function (Kamberaj, Reference Kamberaj2018).
Force field
The molecular mechanics method, also known as molecular force field, is the application of classical physics to explain the microstructure and macroscopic properties of atoms and molecules, and it has become an irreplaceable research tool in scientific exploration. The accuracy of the molecular force field was essential for describing both bonded and nonbond interactions in proteins (Ding et al., Reference Ding, Yu and Huang2023; Latour, Reference Latour2014; Lopes et al., Reference Lopes, Guvench and MacKerell2015). The energy of the whole molecule was the sum of the bond stretching energy, bending energy, twisting energy, and nonbond interaction (Hwang et al., Reference Hwang, Lee, Lee, Ma, Kang, Cho, Kim, Kwon, Yoon, Kang, Yoon, Nam, Kim, In, Chai, Acree, Grant, Gibson, Jhon, Scheraga and No2020). The relevant parameters and expressions constituted the force field, that was, the parameters representing the interactions between molecules and atoms (Collier et al., Reference Collier, Piggot and Allison2020). A complete force field contained a set of functions and related parameters. The parameters in its potential function were derived from experimental results and the methods of quantum mechanics. Different force fields used different parameters and function forms, resulting in different systems and development directions of force fields. For example, the AMBER force field was first proposed by Kollman et al. in 1984 (Weiner et al., Reference Weiner, Kollman, Case, Singh, Ghio, Alagona, Profeta and Weiner1984), which initially only provided corresponding atomic types for protein and nucleic acid systems. In 1990, the force field parameters suitable for polysaccharide simulation were developed (Homans, Reference Homans1990), and around 2000 atomic types and parameters suitable for small organic molecules were added (Wang et al., Reference Wang, Cieplak and Kollman1999). Amber force field model was a relatively classical force field model in MD. The CHARMM force field was proposed by Karplus group in 1983 and is suitable for the calculation and simulation of various molecular properties, from isolated small molecules to solvated large molecules, but it was not suitable for metal complexes (Brooks et al., Reference Brooks, Bruccoleri, Olafson, States, Swaminathan and Karplus1983). In addition, the CVFF force field proposed by the Dauber-Osguthorpe group in 1988 applied to small molecules and protein systems and could be extended to simulate some inorganic systems, such as silicates, aluminosilicates, and phosphate-aluminum compounds (Brooks et al., Reference Brooks, Bruccoleri, Olafson, States, Swaminathan and Karplus1983), which were mainly used to predict the structure and binding free energy of molecules. The selection of force field should refer to the function form of force field, combine the calculation experience and the characteristics of the system, and select the most suitable force field. The AMBER force field and CHARMM force field are the most used in protein MD, especially in the simulation of de novo design proteins.
The overestimation of secondary structures and excessive conformations were the main problems of the early version force field (Best et al., Reference Best, Buchete and Hummer2008; Piana et al., Reference Piana, Klepeis and Shaw2014; Uversky, Reference Uversky2013). Many efforts were made in the development of the traditional force fields to improve protein simulation. Wu et al. developed a residue-specific force field (RSFF2) based on conformational free-energy distributions of the 20 amino acid residues from a protein coil library (Zhou et al., Reference Zhou, Jiang and Wu2015). They made modifications to the Amber 99SB force field parameters, and it successfully folded α-helical structures better than RSFF1 and Amber 99SB force field. It also provided folding enthalpies and entropies in reasonably good agreement with available experimental results. We believed that such progress could be very helpful in the simulating designed proteins, particularly those with specific entropies. Shaw et al. described a force field that substantially improved on the state-of-the-art accuracy for simulations of disordered proteins without sacrificing accuracy for folded proteins, thereby broadening the range of biological systems amenable to MD simulations (Robustelli et al., Reference Robustelli, Piana and Shaw2018). Regarding the CHARMM force field, many computational biologists have made efforts to modify the parameters to perfection. Lin et al. optimized the backbone parameters of the CHARMM36 force field, aiming to solve the gas properties of the alanine dipeptide, resulting in the Drude-2019 protein FF (Lin et al., Reference Lin, Huang, Pandey, Rupakheti, Li, Roux and MacKerell2020). Their results showed that the updated Drude-2019 protein FF yields smaller overall root-mean-square differences of proteins compared to the additive CHARMM36m and Drude-2013 FFs, as well as similar or improved agreement with experimental NMR properties, allowing for long time scale simulation studies of proteins and more complex biomolecular systems in conjunction with the remainder of the Drude polarizable FF. Apart from force fields mentioned above, examples with improved dihedral angle parameters include Amber ff14SB (Maier et al., Reference Maier, Martinez, Kasavajhala, Wickstrom, Hauser and Simmerling2015), Amber ff99SB-ILDN (Lindorff-Larsen et al., Reference Lindorff-Larsen, Piana, Palmo, Maragakis, Klepeis, Dror and Shaw2010), CHARMM36m (Huang et al., Reference Huang, Rauscher, Nawrocki, Ran, Feig, de Groot, Grubmüller and MacKerell2017), CHARMM22 (Piana et al., Reference Piana, Lindorff-Larsen and Shaw2011), OPLS3 (Harder et al., Reference Harder, Damm, Maple, Wu, Reboul, Xiang, Wang, Lupyan, Dahlgren, Knight, Kaus, Cerutti, Krilov, Jorgensen, Abel and Friesner2016), and OPLS-AA/M (Robertson et al., Reference Robertson, Tirado-Rives and Jorgensen2015).
Designed proteins had not been selected by evolution, the folding energy barrier was generally lower than that of natural proteins, and the composition of hydrophobic cores might be more complex than that of natural proteins (Cao et al., Reference Cao, Goreshnik, Coventry, Case, Miller, Kozodoy, Chen, Carter, Walls, Park, Strauch, Stewart, Diamond, Veesler and Baker2020; Silva et al., Reference Silva, Yu, Ulge, Spangler, Jude, Labão-Almeida, Ali, Quijano-Rubio, Ruterbusch, Leung, Biary, Crowley, Marcos, Walkey, Weitzner, Pardo-Avila, Castellanos, Carter, Stewart, Riddell, Pepper, Bernardes, Dougan, Garcia and Baker2019). All these features require high accuracy of protein folding simulation parameters. Therefore, the optimization and development of force field parameters are important aspects of de novo-designed proteins.
Water models
In the field of MD, the classification of force fields of water molecules is of unparalleled importance compared to other types of molecules. At the same time, the number and variety of classical force fields for water molecules are unmatched, with more than 100 different types. Similar to force fields for other molecules, different levels of approximation can be used to obtain different types of force fields for water molecules. The classical force fields for water molecules are roughly categorized into the following four types: rigid force fields, flexible force fields, polarizable force fields, and dissociative force fields.
The TIP3P water model has been widely used in MD simulations of proteins, including de novo protein design (Park et al., Reference Park, Robinson and Kim2022). However, the prevalence of TIP3P may vary depending on the research community and the specific application. In recent years, there has been a trend toward using newer water models, such as TIP4P-2005 and TIP5P, which more accurately reproduce the properties of liquid water. These models may provide better results for specific properties or phenomena that are sensitive to the water model. The TIP4P-2005 and TIP5P water models are more complex than TIP3P, with additional interaction sites that allow for a more accurate description of water’s tetrahedral structure and dynamics. These models have been shown to better reproduce the properties of liquid water, such as its density, heat capacity, and diffusion coefficient. Max F Döpke et al. demonstrated that the TIP4P-2005 model was among the best rigid water force fields, whereas many of the most successful ion parameters were optimized in combination with SPC-E, TIP3P, or TIP4P-Ew water (Döpke et al., Reference Döpke, Moultos and Hartkamp2020).
The TIP5P water model is a five-site model that was introduced to more accurately capture the tetrahedral structure of liquid water and the properties of water in various phases. One of the most significant features of the TIP5P model is its ability to better reproduce the tetrahedral arrangement of water molecules in the liquid phase. This is achieved by using four charge sites (two hydrogen atoms and two lone pairs) arranged tetrahedrally around the central oxygen atom.
Overall, the choice of water model depends on the specific research question and the desired level of accuracy. While TIP4P-2005 and TIP5P may provide improved accuracy for certain properties or phenomena, they also require more computational resources and may not be necessary for all applications. James Dix et al. concluded that, among the MD models tested, the TIP4P-2005 and TIP5P force fields were for now the most reliable when simulating water under confinement (Dix et al., Reference Dix, Lue and Carbone2018). Researchers should carefully consider the tradeoffs between accuracy and computational efficiency when choosing a water model for protein design and simulation.
MD Simulation to study protein thermal stability
Simulating proteins at different temperatures often results in different average properties and conformational systems. The general approach to studying the protein thermal stability using MD simulations involves performing simulated sampling at temperatures of interest, then comparing the sampled ensemble and analyzing the interactions. This can include calculating structural fluctuations with root mean square displacement (RMSD), root mean square fluctuation (RMSF), and cyclometric radius (Rg) of atomic positions (Ahamad et al., Reference Ahamad, Gupta and Kumar2022; Su et al., Reference Su, Sun, Wang and Shen2022). Additionally, the solvent accessibility area (SASA) can used to obtain more information about protein folding and hydrophobicity, while hydrogen bond statistics can be used to reflect local structural changes. Furthermore, free energy calculations can be utilized to determine the most stable conformation in the ensemble. The comprehensive application of these analytical methods is helpful for exploring the possible molecular mechanism of thermal stability.
There have been numerous MD simulation studies analyzing the heat resistance mechanisms of natural proteins in detail. For example, based on the simulation results of firefly luciferase and its mutants, Rahban et al. concluded that the increase in thermal stability of the enzyme was due to the increase of salt Bridges and hydrogen bonds in the molecule (Rahban et al., Reference Rahban, Salehi, Saboury, Hosseinkhani, Karimi-Jafari, Firouzi, Rezaei-Ghaleh and Moosavi-Movahedi2017). In another study, Kumar et al. utilized MD simulation to investigate the molecular basis of the thermal stability of the natural protein SazCA. They calculated hydrogen bond, SASA, Ramchandran plot to trace the folding and unfolding paths of the protein and tried to explain the highest thermal stability at 353 K using free energy calculation (Kumar and Deshpande, Reference Kumar and Deshpande2021). MD simulations can also be used directly to explore potential causes of the macroscopic properties of de novo-designed proteins. For example, the design of highly stable proteins often favored the encapsulation of hydrophobic residues (Kuhlman and Baker, Reference Kuhlman and Baker2004). Protein Top7 serves as a classical example (Kuhlman et al., Reference Kuhlman, Dantas, Ireton, Varani, Stoddard and Baker2003). Bochek et al. combined experimental and computational methods to investigate the effects of low PH, high temperature, and the structural stability and dynamics of Top7 and its mutants in order to obtain atomic details of their conformational changes under extreme conditions (Boschek et al., Reference Boschek, Apiyo, Soares, Engelmann, Pefaur, Straatsma and Baird2009; Soares et al., Reference Soares, Boschek, Apiyo, Baird and Straatsma2010). The results showed that the mutant of Top7 demonstrated obvious folding under different PH and temperature conditions. Although the helix content was lower at high temperature, the antiparallel β-strand could still fold and form extensive interaction with loop. It was concluded that the molecular basis of tolerance of the Top7 mutant was a very stable hydrophobic β-fold core.
In addition to studying hydrophobic interactions, MD simulations have also been used to investigate another class of mechanisms by which de novo-designed proteins maintain thermal stability: electrostatic interactions. Missimer et al. studied the static mechanism of de novo-designed protein (17 peptide) to maintain stability at high temperature through MD simulations (Missimer et al., Reference Missimer, Steinmetz, Baron, Winkler, Kammerer, Daura and van Gunsteren2007). The solvent model was used to simulate the protein at three temperatures, and the stability difference between the monomer and the trimer was compared. It was found that the RMSF of the trimer was essentially the same at 278 K and 330 K, while the RMSF of the monomer increased with temperature. It was also observed that the average occupancy of the salt bridge in the trimer helix was higher than that of the monomer, indicating that the salt bridge, as a key electrostatic interaction, significantly contributed to the thermal stability of the trimer. Additionally, the interaction between the trimer and water weakened with increasing temperature, while the monomer showed the opposite trend. This led to the suggestion that desolvation at higher temperatures allowed electrostatic interactions to occur in the form of salt Bridges, acting as a major factor in the thermal stability of proteins.
While many of the aforementioned studies aimed to use MD simulations to explore the thermal stability of de novo-designed proteins, the majority of research focused on proteins automatically designed by sequences based on physical energy functions. There is a lack of analysis regarding stability factors for proteins automatically designed based on statistical energy functions. Unlike physical energy models, statistical energy used for sequence design may capture implicit information in the training data and reflect it on the macroscopic properties of the designed protein. By conducting in-depth studies of these designed proteins, it is possible to gain insights into this information, and the new insights obtained can be iteratively incorporated into new design methods to improve and optimize the design process.
Free energy calculations have now become rapid enough to significantly impact protein design. The calculation of free energy represents an important research direction in the field of protein design, particularly for understanding and predicting protein-ligand interactions. In recent years, significant progress has been made in the field of protein design due to the development of computational technology, leading to advancements in free energy calculation. Binding free energy is the most critical physical quantity for describing the process of recognition and assembly. For example, in the process of drug design, the binding free energy of candidate drugs and receptor proteins in the database is often calculated through molecular docking or MD, and potential drugs can be identified based on the results of these calculations.
MD Simulation to study designed protein dynamics
Protein folding describes the process through which one-dimensional structural sequences attain their three-dimensional structure with biological functions, whereas protein design commences from the structure and delves into which sequence can achieve a specific structure. Consequently, research on designing protein folding provides deeper insights into these two aspects and offers more valuable information to the field of protein research.
Membrane proteins play an important role in nature, and the design and simulation of membrane proteins represent essential tools for studying this class of proteins. In recent years, scientists have made significant progress in studying the folding mechanism of newly designed membrane proteins. For instance, in terms of the de novo design, DeGrado et al. identified a stable packing motif in MD simulations of pentameric α-helical barrel phospholamban and designed a new artificial transmembrane parallel pentamer using the sequence pattern (Mravic et al., Reference Mravic, Thomaston, Tucker, Solomon, Liu and DeGrado2019). The MD study in this research indicated that van der Waal interaction contributed to the packing of apolar residues between α-helices and demonstrated that the geometry of helices in such packing was strictly determined, showing that MD simulations were useful for evaluating the dynamics and stability of membrane-embedded protein structures (Niitsu and Sugita, Reference Niitsu and Sugita2023). Ulmschnerder et al. illustrated that membrane binding, insertion, association, and dissociation, all of which were included in the folding dynamics events, can be sampled through MD simulations (Chen et al., Reference Chen, Melo, Berglund, Khan, de la Fuente-Nunez, Ulmschneider and Ulmschneider2020). The designed a 14-residue pore-forming peptide with sequence pattering and conducted high-temperature MD simulations for this de novo-designed protein, making a comparison of several sequences in the long-time MD simulations. The long-time MD simulations revealed split folding events, aiding scientists in understanding the folding mechanisms of the de novo-designed proteins (Chen et al., Reference Chen, Starr, Troendle, Wiedman, Wimley, Ulmschneider and Ulmschneider2019). The pioneering work of Voth and coworkers involved testing a hypothesis that protons were conducted through dry spots by forming transient water wires, often highly correlated with the presence of the excess proton itself in the water wire. To test this hypothesis, they used MD simulations to design transmembrane channels with stable water pockets interspersed by apolar segments capable of forming flickering water wires, providing a compelling example of a de novo proton channel that allowed exploration of the roles transient water wires through apolar regions play in proton selectivity and conduction using the MD tool (Kratochvil et al., Reference Kratochvil, Watkins, Mravic, Thomaston, Nicoludis, Somberg, Liu, Hong, Voth and DeGrado2023).
In addition to the study of membrane proteins, other de novo protein simulation studies have also made significant progress. For example, Schiffer et al. studied the dynamics of explicitly solvated designed proteins through all-atom MD simulations to gain insight into the causes leading to the low affinity or instability of most of these designs. Their simulations ranged from 500 to 1000 ns per replicate on 37 designed protein variants encompassing two distinct folds and a range of ligand affinities, resulting in more than 180 μs of combined sampling. The simulations provided retrospective insights into the properties affecting ligand affinity, proving useful in guiding further steps of design optimization and demonstrating that MD could act as a screening step in protein design, resulting in a more efficient process (Barros et al., Reference Barros, Schiffer, Vorobieva, Dou, Baker and Amaro2019).
Additionally, DeGrado, Handel, and colleagues designed a three-helix bundle, α3D, whose interior sidechains consisted of a diverse set of apolar residues that packed in a geometrically complementary manner (Korendovych and DeGrado, Reference Korendovych and DeGrado2020). Due to its relatively simple but cooperatively folded globular structure, α3D quickly became a widely studied protein for computational and experimental studies of protein folding. Labs that contributed to this work included Shaw, Daggett, and Brooks, and using MD simulations, they assessed the consistency of folding time and protein stability. α3D exhibited distinct characteristics, and the folding enthalpy obtained by MD simulation was in good agreement with the experimental estimation, highlighting the important role of MD simulation in the study of de novo design protein folding (Piana et al., Reference Piana, Klepeis and Shaw2014). In a recent paper, Ulas et al. used MD and free energy calculations to design low to sub-nM drug-binding proteins with a very high success rate. This contrasts with the very low success rates and weak binders obtained in recent studies in which researchers have relied exclusively on design calculations (Lu et al., Reference Lu, Gou, Tan, Mann, Yang, Zhong, Gazgalis, Valdiviezo, Jo, Wu, Diolaiti, Ashworth, Polizzi and DeGrado2023).
In another study, Zhang et al. utilized all-atom REMD simulations, through which several comparably stable intermediate states were observed along the folding free-energy landscape of the de novo design protein (DS119). MD simulations revealed that when two unfolded DS119 proteins bound together, most of the binding sites of the dimeric aggregates were located in the N-terminal segment, particularly residues 5–10, which were predicted to form β-sheet with its own C-terminal segment. The complex folding behavior of DS119 suggested the significant influence of natural selection on protein-folding kinetics, indicating the need for further improvement in rational protein design (Wang et al., Reference Wang, Hu and Zhang2016).
From the above examples, it was evident that in the field of de novo protein design, MD simulation technology had significantly enhanced people’s understanding of the relationship between protein structure and function, offering dynamic information at the atomic level. MD simulation has unveiled the molecular basis of protein function by investigating protein structure and dynamic behavior, capturing the fundamental features of atomic interactions that control movement within proteins to comprehend how proteins move and interact at the atomic level. Furthermore, it has provided dynamic information about conformational transitions, interactions, and control stability and function, which can play a pivotal role in improving the accuracy and efficiency of protein design.
The protein–ligand binding process involved complex structural changes that were challenging to fully capture within the time scales available with MD. To enhance efficiency and reduce the conformational space of the sample, many researchers have utilized geometric constraints, such as funnel-shaped constraints and spherical constraints, to improve the convergence rate of the simulation (Negron et al., Reference Negron, Fufezan and Koder2009). The binding free energy was then accurately calculated by subtracting the influence of constraints through post-processing. Free energy perturbation was a computational method used to estimate the free energy change resulting from small changes in the system (Jespers et al., Reference Jespers, Åqvist and Gutiérrez-de-Terán2021). The importance sampling method improved computational efficiency by altering the sampling distribution, especially when dealing with rare events. The combination of geometric constraints with free energy perturbation and importance sampling methods can enhance the accuracy of protein–ligand binding free energy calculation. The various methods of free energy computation not only aid in comprehending the fundamental properties of proteins but also guide drug design and protein engineering, providing powerful tools in the field of protein design to understand and manipulate the structure and function of proteins.
Conclusions and outlook
MD simulation is a computational technique used to study the dynamic behavior of proteins at the atomic level, serving as a powerful tool in the simulation field of de novo-designed proteins. It has played a crucial role in advancing our understanding of the dynamic behavior, stability, and function of designed proteins, leading to significant advantages in force field modeling and sampling methods.
In the current field of protein molecular simulation, the reliability of simulation results primarily depends on the accuracy of the molecular mechanical force field employed. However, the currently developed force fields are still not without imperfections, and the simulation results of some systems may not align well with experimental data. One method to address this issue is to utilize different force fields or optimize the force field parameters themselves. Additionally, due to limitations in simulation algorithms and computing power, many folding simulation runs are too short to accurately capture the various possible conformational states of artificially de novo-designed proteins. This limitation restricts our understanding of the conformational changes of artificially de novo-designed proteins at different temperatures. To tackle this problem, it may be valuable to consider using more sophisticated simulation algorithms or more powerful computers.
Numerous research examples have demonstrated that MD simulation can, to a certain extent, replica experimental observations regarding the stability and dynamic properties of de novo-designed proteins. It can provide insights into the stability and folding dynamics characteristics of proteins at a microscopic scale, offering a partial reflection on the folding properties of de novo-designed proteins. Taking into account the substantial workload and lengthy cycle involved in experimental verification of de novo-designed proteins, MD simulations can be conducted on the de novo-designed proteins at appropriate temperatures initially. The structural fluctuations can be observed as a preliminary assessment of their structural stability and folding characteristics, allowing for the screening of sequences that may fold into target structures. Subsequently, these sequences can undergo experimental verification.
However, MD simulation methods have limitations when it comes to exploring protein dynamics. One such limitation is sampling constraints: MD simulations are inherently stochastic, and sampling rare events can be a challenging task. This can result in an incomplete or biased view of the system’s dynamics, particularly for processes that occur over long timescales or involve multiple pathways. Another limitation is related to system size: MD simulations are typically restricted to relatively small systems, typically containing a few hundred atoms or less. This can pose a challenge when studying larger proteins or protein complexes, where important dynamic processes may occur on a larger scale. Additionally, computational resource limitations and force-field accuracy present further constraints. However, we believe that with the continuous development of computer software and hardware, along with the widespread adoption of artificial intelligence algorithms, these limitations will gradually diminish in the field of protein simulation. In conclusion, MD simulation has emerged as a powerful tool in the field of de novo protein design, playing a crucial role in advancing our understanding of the dynamic behavior, stability, and function of designed proteins. In particular, in terms of structure refinement, MD simulations have been found to be effective in enhancing the structures of de novo-designed proteins. They can assist in tasks such as relaxing initial models, resolving steric clashes, and optimizing side-chain conformations, contributing to the generation of more accurate and stable protein structures. MD simulations have provided valuable insights for the design of de novo proteins by helping researchers understand the dynamics and stability of protein structures. This technique can be utilized to predict the stability of de novo protein designs. By simulating the protein’s dynamics over time, researchers can identify regions of the protein that may be unstable or prone to unfolding. Additionally, by simulating the dynamics of protein–protein complexes, researchers can identify key residues involved in the interaction and optimize them for stronger binding. This provides valuable insights into the dynamics and stability of de novo proteins and can aid in the design of novel proteins with specific functions and improved stability.
MD simulations provide valuable insights into the stability of designed proteins by monitoring fluctuations in the protein’s backbone and side chains. This information is critical for identifying potential weaknesses and improving the design. Continued advancements in enhanced sampling techniques, such as accelerated MD and metadynamics, will enable researchers to explore rare events and long-timescale dynamics more efficiently.
In summary, MD simulation is poised to remain a cornerstone in the field of de novo protein design, playing a pivotal role in advancing our ability to create novel proteins with tailored functions and applications in diverse fields, from biotechnology to medicine. As computational methods and experimental techniques continue to evolve, the synergy between simulation and experimentation will undoubtedly lead to groundbreaking discoveries and innovations.
Funding
This research was funded by National Natural Science Foundation of China, grant number 82101922.
Competing interest
The authors declare no conflict of interest.
