As part of our ongoing efforts to support diverse force fields and simulation programs in CHARMM-GUI, this work presents the development of FF-Converter to prepare Amber simulation inputs with various Amber force fields within the current CHARMM-GUI workflow. The currently supported Amber force fields are ff14SB/ff19SB (protein), Bsc1 (DNA), OL3 (RNA), GLYCAM06 (carbohydrate), Lipid17 (lipid), GAFF/GAFF2 (small molecule), TIP3P/TIP4P-EW/OPC (water), and 12-6-4 ions, and more will be added if necessary. The robustness and usefulness of this new CHARMM-GUI extension are demonstrated by two exemplary systems: a protein/N-glycan/ligand/membrane system and a protein/DNA/RNA system. Currently, CHARMM-GUI supports the Amber force fields only for the Amber program, but we will expand the FF-Converter functionality to support other simulation programs that support the Amber force fields.
INTRODUCTION
Molecular modeling and molecular dynamics (MD) simulation methods have become popular and powerful tools in modern biomolecular science, as current state-of-the-art MD simulations and force fields (FFs) have been validated and mature enough to interpret experiments and guide new experiments with testable hypotheses. During MD simulations, one can follow the motion of a molecule or an ensemble of molecules comprising up to millions of atoms at the atomic resolution and at the femtosecond time resolution. Since such a resolution is difficult or sometimes even impossible to achieve in experiments, MD simulations give a complementary view on biomolecular systems, which is why these methods were even called a “computational microscope.”1 Molecular simulations, often in combination with experiment, can help to model and predict the structure of biomolecules and allow one to extract and analyze kinetics and thermodynamics of biomolecular systems.
However, these computational tasks often require the setup of a complex series of computational subtasks and scripts ranging from the preparation of a simulation starting structure, system equilibration, production simulations, and appropriate data gathering. Especially as larger spatial scales, longer time scales, and higher levels of realism become possible and necessary, generation of realistic biomolecular systems can be a time-consuming and error-prone process for beginners but also for expert users. CHARMM-GUI2,3 (http://www.charmm-gui.org) has at least partially, if not entirely, addressed these challenges by developing a web-based graphical user interface (GUI) workflow to interactively build complex multi-component biomolecular systems and prepare input files for state-of-the-art MD simulations using commonly used simulation packages (e.g., CHARMM,4 Amber,5 NAMD,6 GROMACS,7 GENESIS,8 LAMMPS,9 Desmond,10 and OpenMM11).
In recent years, several CHARMM-GUI modules have been developed to support the setup of complex membrane protein systems including bacterial outer membranes12–15 and the fitting of protein subunits into Cryo-EM or x-ray density maps,16 complex glycan systems,17–19 and coarse-grained simulation systems.20–22 Currently, CHARMM-GUI relies exclusively on the additive CHARMM FF23–30 for all-atom system building, as it covers most molecular building blocks of life and various small molecules. Therefore, so far, CHARMM-GUI has provided a comprehensive set of simulation system and input files using only the CHARMM FF for the aforementioned simulation packages.
In this work, we present an extension of CHARMM-GUI that allows users to use the Amber package in combination with a variety of Amber FFs31–40 (e.g., ff14SB/ff19SB, Bsc1/OL3, GLYCAM06, Lipid17, GAFF/GAFF2, TIP3P/TIP4P-EW/OPC, and 12-6-4 ions). The Amber package5 is one of the most widely used biomolecular simulation packages with program modules that very effectively employ graphical processing units (GPUs) to propagate MD trajectories, e.g., pmemd.cuda.41,42 CHARMM-GUI also supports hydrogen mass repartitioning (HMR)43 and appropriately handles input structures containing the Amber FF residue names. In addition to the generated multi-component biomolecular systems that can consist of various combinations of subunits, including multiple proteins, nucleic acids, carbohydrates, small molecules, and lipids, CHARMM-GUI provides all necessary input files for minimization, equilibration, and production simulations using the Amber suite of simulation programs.
In the following, the new CHARMM-GUI functionality for the Amber FFs and how it can be used for input and script generation are described. To demonstrate the robustness and usefulness of the new CHARMM-GUI extension, two exemplary systems are considered, and their simulation results are described: a protein/N-glycan/ligand/membrane system and a protein/DNA/RNA system.
METHODS
Conversion of CHARMM FF compatible system into Amber FF compatible system
The current CHARMM-GUI workflow prepares simulation systems with the CHARMM FF, i.e., all residue and atom names of the generated system follow the CHARMM FF. In order to support the Amber FFs (and later other FFs) within the CHARMM-GUI workflow, a best practical procedure was chosen to convert a CHARMM FF compatible system to an Amber FF compatible system by mapping the residue and atom names in the CHARMM FF to those in the Amber FF. A fundamental difference between the CHARMM and Amber FFs is that the CHARMM FF uses a patching scheme to modify a certain residue, whereas the Amber FF requires different residue types for every modification. For example, protonation of aspartate can be done by applying an ASPP patch in the CHARMM FF, but the Amber FF requires a different residue (ASH) for a protonated aspartic acid in addition to charged ASP. To make the conversion procedure robust, we have developed an in-house program (FF-Converter) that takes CHARMM PSF (protein structure file) and CRD (CHARMM coordinate file) as an input and provides Amber FF compatible structure files for the Amber package. Briefly, CHARMM PSF and CRD files are used to generate an internal data structure containing topology (atoms, residues, chains, bonds, angles, dihedral angles, improper angles, and CMAP) and coordinate information of the system, and the topology and coordinate information is then changed following atom and residue names in the Amber FFs.
More specifically, for protein molecules, the CHARMM FF and Amber FF share the same residue and atom names in most cases, so only mismatched residues are renamed. For example, a preferred tautomeric state of histidine in the CHARMM FF is HSD (proton on ND1), while the corresponding residue is HID in the Amber FF. In this case, the conversion can be done by simply renaming HSD to HID in the data structure: similarly, HSE (proton on NE2) → HIE and HSP (protonated) → HIP. Unlike histidine, protonation of other amino acids requires more considerations as their residue names are the same, regardless of protonation states in the CHARMM FF. For example, aspartic acid is called ASP, regardless of protonation states, but protonated ASP has an additional hydrogen atom (HD2) in the CHARMM FF. Therefore, once ASP HD2 is detected by FF-Converter, it is renamed to ASH for the Amber FF. In the Amber FF, there are three different residue types for cysteine: CYS for a standard cysteine, CYM for a deprotonated cysteine, and CYX for a disulfide bond. To properly identify these residue types, FF-Converter first checks if a hydrogen atom is attached to sulfur of cysteine in the CHARMM PSF. If not found, the cysteine residue is renamed to CYM. If a S–S bond between two cysteine residues is detected, the residue names are changed to CYX.
Similarly, the terminal capping of protein requires extra considerations as the CHARMM and Amber FFs use very different approaches. In the Amber FF, the protein termini are capped by individual residues especially designed for termini. For example, to cap the N-terminal alanine residue with , the Amber FF uses NALA residue that is a modified alanine residue for terminus. In addition, to cap the N-terminus of alanine with an acetyl group (–CHOCH3), the Amber FF uses ACE, a separated residue for acetyl group, and makes a link between the acetyl group and the alanine residue. However, to apply the same terminal capping, the CHARMM FF uses NTER or ACE patche to alanine. To identify which terminal group is used, FF-Converter uses an internal library that contains a list of atoms for each terminal group. Once a terminal group in the CHARMM FF is identified, each terminal residue is renamed or split into two residues depending on the terminal group type in the Amber FF.
For RNA/DNA molecules, the CHARMM FF uses the standard RNA residue names and a patch is applied to RNA residues to make DNA residues, whereas the Amber FF has separated residues for both RNA and DNA. Therefore, the system conversion is done by first checking the existence of a hydroxyl group at C2′ of ribose. If exists, the residue name is changed to a corresponding RNA residue in the Amber FF; otherwise, it is changed to a corresponding DNA residue.
For the conversion of lipid molecules, charmmlipid2amber.csv in AmberTools is utilized. The Amber lipid FF uses discrete residues for lipid head and tail groups and combines three residues to construct a single lipid molecule, whereas the CHARMM FF uses a single residue for each lipid molecule. Therefore, the conversion is accomplished by splitting a CHARMM lipid residue into three pieces using the residue/atom mapping information in charmmlipid2amber.csv file.
Handling carbohydrate residues and their connectivity properly is the most demanding process in converting a CHARMM FF compatible system to an Amber FF compatible one. The Amber FF uses the GLYCAM FF37 for carbohydrates. The GLYCAM FF has a distinctive three-letter naming scheme for sugar residues: one letter for linkage sites (see Table 2 in http://glycam.org/docs/forcefield/glycam-naming-2/), one letter for the glycan type (Table S1), and one letter for the anomeric state (A/B for pyranoses and D/U for furanoses). For example, VGA is alpha (A) d-glucose (G) having glycosidic linkages with other sugars at C3 and C6 positions (V). The CHARMM FF uses a unique residue name for each standard sugar type and applies patches for glycosidic linkages and chemical modifications. Note that Glycan Reader and Modeler in CHARMM-GUI takes care of reading carbohydrates with proper linkages and chemical modifications from a PDB file or modeling them in silico. The conversion of the CHARMM FF residue name to the GLYCAM FF residue name is achieved by (i) detecting linkage sites using bond connectivity information, (ii) identifying a glycan type corresponding to the GLYCAM FF, and (iii) determining an anomeric state of the sugar. Based on the determined residue name, heavy atoms are renamed if necessary. For terminal capping of reducing-end sugar, the GLYCAM FF uses a separated residue “ROH.” Therefore, the anomeric hydroxyl group of reducing-end sugar is split into a new residue (ROH) and stored in the data structure. In the case of N-/O-linked glycans, the glycosylated protein residue is renamed to the corresponding one in the GLYCAM FF: NLN, OLS, and OLT.
In most cases, residue and atom names are matched between the CHARMM and Amber FFs. However, there are a few cases that the names do not match between the FFs. To map the residue and atom names, we constructed an internal mapping library by comparing the FFs. Note that all hydrogen atoms of protein, DNA, RNA, and carbohydrates in a CHARMM FF compatible system are removed to simplify the internal mapping library, and then, hydrogen atoms are regenerated using internal coordinate information in the Amber FF after conversion is completed. If any non-supported residue is detected during the conversion, FF-Converter issues an error and terminates the conversion process. The converted data structure is then used to generate structure files (parm7 and rst7 files) for Amber simulation.
Web interface implementation
As mentioned above, to support the Amber FFs within the CHARMM-GUI workflow, a best practical procedure was chosen to convert a CHARMM FF compatible system to an Amber FF compatible system. Therefore, the FF selection option (Fig. 1) was added to the input generation page (STEP3 for the solution system and STEP5 for the membrane system). Upon user’s selection of the Amber FF, available FFs for individual system components are displayed, and users can select different FFs according to their preference (Table I).
Components . | Available FFs . |
---|---|
Protein | FF14SB,31 FF19SB |
DNA | BSC132 |
RNA | OL333 |
Lipid | LIPID17 |
Carbohydrate | GLYCAM 06j37 |
Ligand | GAFF, GAFF238 |
Water | TIP3P, TIP4P-EW, OPC35 |
Ions | 12-6-4 (on39,40/off36) |
CHARMM-GUI supports two different methods for ligand FF generation during the system building process with the CHARMM FF: CGenFF and Antechamber.44,45 Regardless of which method is chosen, a ligand FF is re-generated with the selected ligand FF method (GAFF or GAFF2) if users choose to use the Amber FF. The re-parameterization is performed using Antechamber with AM1-BCC charges. The final ligand structure from the system building process is used for the charge calculation. Therefore, it is the user’s responsibility to define the ligand protonation state during the system building process. The re-parameterized ligand FF is then used for the CHARMM to Amber system conversion with FF-Converter.
Hydrogen mass repartitioning (HMR)43 is a widely adopted simulation technique to accelerate the MD simulation by repartitioning atomic masses, allowing the increase in simulation time step up to 4 fs. CHARMM-GUI Input Generator now supports HMR for the Amber FFs, for which hydrogen atom mass increases by a factor of 3, and the increased mass is subtracted from the linked heavy atom mass.
Interactions involving multi-valent ions are often underestimated due to the neglect of the ion-induced dipole interaction. To overcome this issue, Li et al. introduced the 12-6-4 ion model,39,40 which has an additional r−4 term for the Lennard–Johns potential to describe the ion-induced dipole interaction. FF-Converter supports various 12-6-4 ion models in the Amber FFs when the option is activated.
RESULTS AND DISCUSSION
Energy comparison between the inputs prepared by LEaP and FF-Converter
FF-Converter produces parm7 and rst7 files without using LEaP in AmberTools.5 To ensure and validate the accuracy of the parm7 file generated by FF-Converter, the energies calculated using the inputs prepared by LEaP and FF-Converter were compared. The energies of 5O8F and 6O0Z systems were examined with the final simulation coordinates. The PDB file provided by FF-Converter was used to prepare the parm7 file with LEaP. As shown in Table S2, there is no difference in the energies except for a slight variation in the improper angle energy. Such a variation results from the difference in the order of atoms that comprises the improper angle. There are two different topology file types available in the Amber FFs; one is a LIB file and the other is a PREP file. While the latter contains an explicit improper angle list for each residue, the former does not contain improper angle information. Therefore, when a LIB file is used in LEaP, the improper angle is constructed by searching for potential improper angle candidates and verifying that the FF includes a parameter for the detected improper angle. While FF-Converter uses the same strategy, LEaP rearranges the non-central improper angle atom order after applying the FF parameters. For example, for the improper angle composed of C–H–*Ng–Cg atom types (* for the central atom), LEaP applies the corresponding FF parameter and rearranges the atoms to be C–Cg–*Ng–H. FF-Converter follows the atomic order in the FF parameter, resulting in a slight variation in the improper angle energy compared to LEaP.
Testing of Amber force fields ff14SB, Lipid17, GLYCAM, and GAFF2
PDB ID 5O8F46 is a crystal structure of a GABAA receptor, a glycosylated transmembrane protein in complex with several pregnanolone molecules. Running an MD simulation of such a complex target requires a combination of several different Amber FFs, such as ff14SB,31 LIPID17 (unfortunately, only LIPID1434 has been published so far), TIP3P,35 GLYCAM_06j,37 and GAFF2.38 For this reason, we chose PDB:5O8F to test the automatically generated system and input files for simulations with the Amber FF family.
GABAA receptors are large, pentameric, ligand gated ion channels, enabling permeation of Cl− ions through the cell membrane if activated by neurotransmitter GABA (γ-aminobutyric acid).47 A top view (extracellular side) of the protein complex is shown in Fig. 2(a), with individually colored protein chains. Upon switching from an apo to a ligand bound form, the protein complex opens the pore and allows ion permeation. This opening-process can be modulated by neurosteroids, such as pregnanolone.48 The binding sites of modulators and agonists are situated on different domains of the receptor, allowing simultaneous binding.49 This enables natural or artificial modulators to allosterically influence neurological processes, such as stress response (anxiety) or sedation. For this reason, GABA receptors are prime targets for psychoactive drugs, such as benzodiazepines.50
Membrane Builder was used to embed the GABAA receptor complex in a membrane bilayer consisting of 1016 POPC and 127 cholesterol molecules in a 150 mM KCl solution [Fig. 2(b)], a popular and simple mimic for eukaryotic plasma membranes in MD simulations.51–53 An online video tutorial explaining the building process in detail is available (http://www.charmm-gui.org/demo/amber_ff/1). The entire system was modeled using ff14SB31 (proteins), LIPID17 (lipids), TIP3P35 (water), GLYCAM_06j37 (carbohydrates), and GAFF238 (pregnanolone). The cutoff distance for non-bonded interactions was set to 9 Å, after which Coulombic interactions were treated with the particle mesh Ewald method,54 and effects of long-range van der Waals interactions were estimated by using a dispersion correction model. Periodic boundary conditions were employed. Using the equilibration and production inputs generated by CHARMM-GUI, we first ran the established seven-step minimization/equilibration process,3,12 and then, two independent 1-µs production runs were performed using the GPU-accelerated CUDA55 version of Amber18, pmemd.5 Specifics on each of the equilibration steps as well as all used input files can be found in the supplementary material (SI-5O8F.zip).
Since we only wanted to assess if the CHARMM-GUI outputted files can be used for simulations without further processing, we did not assign special protonation states to titratable residues but used the standard states instead (i.e., charged Asp, Glu, Lys, Arg, and His protonated on Nδ). Note that CHARMM-GUI is also compatible with PDB files utilizing the Amber residue naming conventions, allowing users to use the output of the PDB2PQR server (http://nbcr-222.ucsd.edu/pdb2pqr_2.0.0/) as an input PDB file in CHARMM-GUI. In other words, CHARMM-GUI properly handles different Amber residue names for different protonation states (e.g., ASP and ASH for charged and protonated Asp, respectively). An online video tutorial explaining the usage of the output from the PDB2PQR server in detail is available (http://www.charmm-gui.org/demo/amber_ff/3).
In order to enable a simulation time step of 4 fs, the HMR approach43 was used in combination with the SHAKE algorithm.56 For the production simulations, we used the current standard input file for NPT simulations, generated by CHARMM-GUI, i.e., temperature control by using the Langevin thermostat57 with a friction coefficient of 1.0 ps−1 and semi-isotropic pressure control by using the Berendsen barostat58 with a relaxation time of 1.0 ps. The target temperature and pressure were set to 303.15 K and 1.0 bar, respectively. As an alternative pressure control for membrane simulations, the semi-isotropic Monte Carlo (MC) barostat59 could be used and is more efficient for simulations on GPUs. However, we strongly recommend using the Berendsen barostat for membrane simulations with LIPID14 and LIPID17, since these lipid FFs have been parametrized and tested only with the Berendsen barostat in mind.34 According to the authors of the LIPID14 paper, LIPID14 produces a realistic membrane surface tension if combined with the Berendsen barostat. We believe that this behavior cannot be guaranteed if the MC barostat is used, so the appearance of artifacts concerning the shape of the membrane surface cannot be ruled out, especially if larger membranes are simulated with the MC barostat (see also Fig. S1). Since CHARMM-GUI will be updated as soon as new FFs become available, the standard setting for pressure coupling may be changed if the next generation of the Amber lipid FF can support the MC barostat.
For trajectory analysis, root mean square deviations (RMSD) were calculated using VMD,60 while root mean square fluctuations (RMSF) were evaluated using CPPTRAJ.5 For RMSD calculations, each frame in the trajectories was first aligned with the initial structure for the selected atoms (heavy atoms for ligands/glycans and backbone atoms for proteins). For RMSF calculations, the trajectories have been aligned along the protein Cα atoms. Snapshots were rendered using VMD,60 and the plots were prepared with xmgrace (http://plasma-gate.weizmann.ac.il/Grace/).
To assess the stability of the simulated system on longer time scales, after the initial equilibration phase, we performed two independent 1-µs long simulations. Both simulations featured stable densities of 1.02 g/cm3 and box sizes of about 144 × 144 × 161 Å3 [a time-series of the x (=y) dimension is depicted in Fig. S2]. Figure 3 shows RMSD time-series of each simulation with respect to the initial structure, suggesting that all parts of the simulated protein/glycan/pregnanolone complex were very stable throughout the sampling phase of both simulations. Mean heavy atom deviations of 2.1 Å and 2.4 Å (black lines) indicate the high structural stability of the system. The glycans (green lines) show the highest deviations with averages of 3.0 Å and 3.2 Å due to their flexible nature. All pregnanolone ligands remained docked to their respective binding sites throughout the entirety of both simulations, resulting in very low RMSD with averages of 1.2 Å and 1.1 Å (red lines). The low positional deviations of the ligands indicate that the automatically generated GAFF2 FF with empirical point charges was sufficient to enable stable binding of the ligands. The high stability of the simulated structure is further illustrated by the RMSF plots (Fig. 4). Only the loop regions exhibit higher fluctuations (especially on the N-terminus), with transmembrane domain helices remaining very stable. We also calculated the radius of gyration for both simulations and the PDB structure. With 41.0 ± 0.1 Å and 41.1 ± 0.1 Å, the radius of gyration is very similar to the PDB structure (40.4 Å) in both simulations, indicating a low mobility overall. With this illustrative simulation of PDB:5O8F in a bilayer, we demonstrate that the membrane simulation system and Amber equilibration and production input files with FF14SB, Lipid17, GLYCAM_06j, TIP3P, and GAFF2 can be reliably generated by CHARMM-GUI and used without further modification to perform μs time scale simulations of chemically very complex biomolecular systems.
Testing of Amber force fields ff14SB, BSC1, OL3, and TIP3P
PDB ID 6O0Z61 is a cryo-EM structure of a protein/DNA/RNA complex with sgRNA-guided cas9 endonuclease from Streptococcus pyogenes in a precatalytic form. The system is illustrated in Fig. 5. We simulated the complex using ff14SB31 (protein), BSC132 (DNA), OL333 (RNA), TIP3P35 (water), and K+/Cl− ions.36 We chose PDB:6O0Z to illustrate and validate the automated translation of the CHARMM FF into the respective Amber FFs in CHARMM-GUI. In particular, we note that PDB:6O0Z has various forms of RNA and DNA. Besides dsDNA and self-pairing RNA, it also contains ssDNA, ssRNA, and even a DNA–RNA hybrid, mostly in a protein environment.
The system was prepared using CHARMM-GUI Solution Builder. A video tutorial and a description of the setup procedure are available online (http://www.charmm-gui.org/demo/amber_ff/2). First, all structurally unresolved sequences were modeled using the functionality of modeling missing residues in PDB Reader and Manipulator. There are no more than 15 consecutive missing residues in the cryo-EM structure, which is an upper bound for building reasonable structures using GalaxyFill62 in CHARMM-GUI. The system was then solvated with at least 10-Å explicit water around the complex in a 150 mM KCl solution with additional K+ ions to neutralize the nucleic acids (340 K+ and 230 Cl− in a box size of 140 × 140 × 140 Å3). Using the Langevin thermostat57 with a friction coefficient of 1.0 ps−1, the equilibration and production runs were conducted at a target temperature of 303.15 K. While the equilibration was performed in the NVT (constant particle number, volume, and temperature) ensemble, the atmospheric pressure was maintained by using a Monte Carlo (MC) barostat59 with isotropic pressure scaling in the production run. The cutoff distance for non-bonded interactions was set to 9 Å in combination with the particle mesh Ewald method54 and periodic boundary conditions. We did not alter the default protonation states since our aim was to assess the stability of the nucleic acids during the simulation, not the specific dynamics of the complex structure.
The initial structure was minimized (5000 steps), equilibration was performed for 50 ns, and the production run was performed for 1 µs using the input, topology, and restart files provided by Solution Builder. During the minimization and equilibration, positional restraints were applied. SHAKE56 and HMR43 were utilized to increase the time step to 4 fs in the production run. Specifics on the related input files can be found in the supplementary material (SI-6O0Z.zip). To assess the stability of the system, the RMSD of the entire complex, as well as selected subunits, was calculated with respect to the initial structure. Also, the RMSF of the production run was determined after alignment using the backbone atoms of the protein and phosphorous atoms of the nucleic acids. All calculations were carried out using CPPTRAJ.5
The overall structure shows an RMSD of around 7–8 Å and was very stable after a fast transition at the beginning of the production run (Fig. 6). The magnitude of the structural deviation was mostly due to a movement of the recognition lobe that exhibited an RMSD of around 8 Å. This movement of the recognition lobe, together with that of the HNH domain, encloses the target DNA/sgRNA hybrid, as shown in Fig. 5. The divergence from the cryo-EM structure could be due to artifacts of vitrification, the buffer solution, or the FFs. Due to its natural high flexibility, the sgRNA deviated by around 4 Å. Other regions deviated significantly less (Fig. 6).
More importantly, the DNA strands show a stable RMSD of below 3 Å on average. The RMSF of the nucleic acids is well below 2 Å on average, and the termini and pairing regions of the sgRNA show higher fluctuations (Fig. 7). All evaluations suggest that the nucleic acids are quite stable in the complex protein environment. The results, therefore, demonstrate that the automatically prepared system and input files using Amber FFs ff14SB, BSC1, OL3, and TIP3P allow the production of microsecond-time scale trajectories without further modification by the user.
CONCLUSIONS
Preparing a realistic biomolecular simulation system often requires considerable knowledge of a complex series of computational subtasks and scripts, which can be a time-consuming and error-prone process even for an expert. To facilitate the simulation system preparation, many web-based and standalone tools have been developed, such as CHARMM-GUI,2 LEaP (AmberTools),5 PACKMOL-Memgen,63 Modeller (OpenMM),11 QwikMD (VMD plugin),60,64 HTMD (Acellera),65 and Discovery Studio (BIOVIA).66 However, there is no unified method or tool to setup a simulation system with different FFs, making it difficult to prepare a simulation with a FF that the user is not familiar with. As part of our ongoing efforts to overcome this issue, we have presented the development of FF-Converter to support the Amber FFs in CHARMM-GUI.
To support the Amber FFs (and later other FFs) within the CHARMM-GUI workflow, FF-Converter was developed to convert a CHARMM FF compatible system to an Amber FF compatible system by mapping the residue and atom names in the CHARMM FF to those in the Amber FF. Note that if any non-supported residue is detected during the conversion, FF-Converter issues an error and terminates the conversion process. For example, glycolipids and lipopolysaccharides are not officially supported in the Amber FF, and thus, conversion of any system containing glycolipid and lipopolysaccharide molecules will fail. In addition, any post-translational modification (e.g., palmitoylation), which is not officially supported in the Amber FFs, will be ignored. The system and input files with the Amber FFs in CHARMM-GUI were tested and validated using two illustrative systems: PDB ID 5O8F for a protein/N-glycan/ligand/membrane system and PDB ID 6O0Z for a protein/DNA/RNA system. Each system showed its stability for the microsecond-time scale of simulation. The energies calculated by the input prepared by CHARMM-GUI were identical to those prepared by LEaP, with a slight variation in the improper angle energy resulting from the difference in the non-central improper angle atom order.
Currently, CHARMM-GUI supports the Amber FFs only for the Amber program, but we will expand the FF-Converter functionality to support other simulation programs that support the Amber FFs. In addition, we will continue to support other Amber FFs that are not included in this work. CHARMM-GUI Amber FF support is only for standard MD simulations, but we are under development to expand this to be utilized for the free energy calculation with the thermodynamic integration method (Amber-TI). Ultimately, we are aiming to develop a unified protocol to prepare complex biomolecular simulation systems with any FFs for any MD simulation engines, lowering the boundary in preparing simulation systems and inputs.
AUTHORS’ CONTRIBUTIONS
J.L., M.H., and M.R. contributed equally to this work.
SUPPLEMENTARY MATERIAL
See the supplementary material for glycan residue names in the CHARMM force filed and the corresponding one letter glycan type in the GLYCAM force field (Table S1), a comparison of energies calculated using the inputs prepared with LEaP and FF-Converter (Table S2), and two archive files SI-5O8F.zip for the equilibration and production inputs of the 5O8F system and Si-6O0Z.zip for the equilibration and production inputs of the 6O0Z system.
DATA AVAILABILITY
The data that support the findings of this study are available within the article (and its supplementary material).
ACKNOWLEDGMENTS
This work was supported, in part, by grants from the NSF (Grant Nos. DBI-1145987 and DBI-1660380), a Friedrich Wilhelm Bessel Research Award from the Humboldt Foundation (W.I.), the Deutsche Forschungsgemeinschaft research unit FOR2290 (Grant No. Za153/27-1), and the Leibniz Supercomputing Centre (Grant No. pr27za; M.Z.).
The authors declare no competing financial interest.