Team:Aalto-Helsinki/Model Setup

Aalto-Helsinki




Simulation Setup

Molecular modeling methods ultimately aim for direct comparison with experimental measurements. As such, a good model of molecular interaction is essential. Quantum chemistry based ab initio molecular dynamics methods aim to reduce the amount of fitting and guesswork required for accurate modeling of molecular interactions. However, such approaches are generally limited to small systems and short timescales due to the added computational demand [17]. In classical molecular dynamics, molecules are described using stick-and-ball models: spherical atoms are connected by springs that represent bonds. As such, molecular forces can be described by simple mathematical models. For example, Hooke's law can be used to describe bonded interactions, while non-bonded interactions can be described by a Lennard-Jones potential. [14, 1]

Molecular dynamics (MD) is based on numerical, step-by-step, evaluation of Newton's equations of motion. Due to the many-body nature of the problem, Newton's equations of motion are discretized and solved numerically. MD trajectories, that describe the time evolution of the system, consist of both position and velocity vectors of the particles in the system. The position and velocity vectors are reevaluated according to a finite time interval by using numerical integrators. The position vectors define the geometric configuration of the system while the velocity vectors define the kinetic energy and temperature of the system. [14, 1, 17]

We investigated the behavior of DCD-1L in aqueous solvents at varying salt concentrations and temperatures using molecular dynamics. the simulations make use of the GROMOS forcefield parameter set 53A6 [11]. GROMOS53A6 is considered a united-atom forcefield, which maps a carbon and its associated apolar hydrogens as one interaction center. The parametrization of the forcefield is based on accurate reproduction of free enthalpies of hydration and apolar solvation of a wide range of compounds. The relative free enthalpy of solvation between apolar and polar environments is an important property in many biological phenomena, including protein folding and lipid membrane formation.

System Setup and Parameters

We probed the behavior of a single DCD-1L solvated in water with 2 added Na+ counter-ions to neutralize the net negative charge of the peptide, as well as added salt (NaCl), with the salt concentrations varying from 20 mM to 500 mM. The salt simulations were carried out at 310 K. Additionally we probed the temperature dependency of DCD-1L structure by simulating a single peptide in water with 2 Na+ counter-ions at temperatures of 290 K, 300 K, 310 K and 320 K. All of the simulation trajectories were calculated using the Gromacs v.4.6.7 open-source simulation package [13, 7].

As a starting structure for our simulations, chose the helical crystal structure of dermcidin derivative DCD-1 derived by Song et al. [15] available in the RCSB protein Data Bank (entry 2YMK). Due to the presence of missing atoms in the original model by Song et al., we constructed a homology model based on the DCD-1L sequence used in our laboratory constructs using the SWISS-MODEL automated protein structure homology-modeling server [3, 10, 2, 6]. The resulting structure was then mapped to the GROMOS53A6 forcefield using the Gromacs 4.6.7 tool pdb2gmx while capping the N-terminal and C-terminal charges of the peptide.

In all simulations, the dimensions of the simulation box were cubic with a side length of 11.2 nm. The peptide was inserted and centered into the simulation box and solvated in water. For all simulations, the simple point-charge (SPC) explicit water model was used. After solvation, two of the water molecules were replaced by 2 Na+ ions to act as counter-ions and effectively neutralize the charge on the simulated system. Additionally, for simulations involving added salt, an adequate number of water molecules were replaced by an equal number of Na+ and Cl- ions to mimic the desired salt concentration of either 20 mM, 50 mM, 70 mM, 100 mM, 120 mM, 150 mM, 200 mM, 300 mM or 500 mM. A visual overview of the simulation setup (simulation cubicle) can be seen in figure 1.

Figure 1. A typical simulation setup visualized using VMD of DCD-1L solvated in water with NaCl. For clarity, the water molecules have been left invisible.

Prior to the production run, each simulation system was energy minimized using the steepest descent algorithm. This energy minimization was followed by a 0.5 ns NVT equilibration with a time step of 2 fs. During NVT equilibration, the temperature was controlled by the stochastic velocity rescaling thermostat developed by Bussi et al. [4] with τT = 0.1 ps. The NVT equilibration was followed by NPT equilibration for 1 ns with a time step of 2 fs. During the NPT equilibration, the pressure of the system was set to 1.0 bar using the isotropic Parrinello-Rahman [12] pressure control with τp = 2.0 ps. The temperature was controlled by the stochastic velocity rescaling thermostat by Bussi et al. [4] with τT = 0.1 ps. During both NVT and NPT equilibration runs, position restraints were applied to the protein structure.

For the production run, the pressure of the system was set to 1.0 bar using the isotropic Parrinello- Rahman [12] pressure control with τp = 12 ps. The temperature was controlled by the stochastic velocity rescaling thermostat by Bussi et al. [4] with τT = 0.1 ps. A time-step of 2 fs was used. For all equilibration and production runs, electrostatic interactions were calculated using the Particle-Mesh Ewald method [5] and periodic boundary conditions were applied. The simulation time for all systems during the production run was 50 ns.

Analysis of the simulation trajectories was carried out using Gromacs 4.6.7 built-in tools and visualization of trajectories was accomplished using VMD [8]. Assignment of secondary structure elements of the peptide was done using a Gromacs interface for DSSP [16, 9].

Figure 2. Sisu supercomputer (Cray XC30) Ⓒ Samuli Saarinen /CSC. The simulation trajectories were calculated using CSC resources, including the Taito (eng. "Skill") computational cluster and the Sisu (eng. "Brute Force") supercomputer. The CSC data center is located in Kajaani, northern Finland.



References

[1] Adcock, S. A., & McCammon, J. A. (2006, May). Molecular dynamics: survey of methods for simulating the activity of proteins. Chem. Rev., 106(5), 1589–1615.
[2] Arnold, K., Bordoli, L., Kopp, J., & Schwede, T. (2006, Jan). The SWISS-MODEL workspace: a webbased environment for protein structure homology modelling. Bioinformatics, 22(2), 195–201.
[3] Biasini, M., Bienert, S., Waterhouse, A., Arnold, K., Studer, G., Schmidt, T., . . . Schwede, T. (2014, Jul). SWISS-MODEL: modelling protein tertiary and quaternary structure using evolutionary information. Nucleic Acids Res., 42(Web Server issue), W252–258.
[4] Bussi, G., Donadio, D., & Parrinello, M. (2007, January). Canonical sampling through velocity rescaling. The Journal of Chemical Physics, 126(1), 014101.
[5] Darden, T., York, D., & Pedersen, L. (1993). Particle mesh ewald: An nlog(n) method for ewald sums in large systems. J. Chem. Phys., 98(12), 10089-10092.
[6] Guex, N., Peitsch, M. C., & Schwede, T. (2009, Jun). Automated comparative protein structure modeling with SWISS-MODEL and Swiss-PdbViewer: a historical perspective. Electrophoresis, 30 Suppl 1 , S162–173.
[7] Hess, B., Kutzner, C., van der Spoel, D., & Lindahl, E. (2008, Mar). GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput., 4(3), 435–447.
[8] Humphrey, W., Dalke, A., & Schulten, K. (1996). VMD: visual molecular dynamics. J Mol Graph, 14(1), 33–38.
[9] Kabsch, W., & Sander, C. (1983, Dec). Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features. Biopolymers, 22(12), 2577–2637.
[10] Kiefer, F., Arnold, K., Kunzli, M., Bordoli, L., & Schwede, T. (2009, Jan). The SWISS-MODEL Repository and associated resources. Nucleic Acids Res., 37(Database issue), D387–392.
[11] Oostenbrink, C., Villa, A., Mark, A. E., & van Gunsteren, W. F. (2004, Oct). A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem, 25(13), 1656–1676.
[12] Parrinello, M., & Rahman, A. (1981, December). Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys., 52(12), 7182–7190.
[13] Pronk, S., Pall, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., . . . Lindahl, E. (2013, Apr). GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics, 29(7), 845–854.
[14] Schlick, T. (2010). Molecular modeling and simulation: an interdisciplinary guide: an interdisciplinary guide (Vol. 21). Springer Science & Business Media.
[15] Song, C., Weichbrodt, C., Salnikov, E. S., Dynowski, M., Forsberg, B. O., Bechinger, B., . . . Zeth, K. (2013, Mar). Crystal structure and functional mechanism of a human antimicrobial membrane channel. Proc. Natl. Acad. Sci. U.S.A., 110(12), 4586–4591.
[16] Touw, W. G., Baakman, C., Black, J., te Beek, T. A., Krieger, E., Joosten, R. P., & Vriend, G. (2015, Jan). A series of PDB-related databanks for everyday needs. Nucleic Acids Res., 43(Database issue), D364–368.
[17] van Gunsteren, W. F., Bakowies, D., Baron, R., Chandrasekhar, I., Christen, M., Daura, X., . . . Yu, H. B. (2006, Jun). Biomolecular modeling: Goals, problems, perspectives. Angew. Chem. Int. Ed. Engl., 45(25), 4064–4092.