Performance of the General Amber Force Field in Modeling Aqueous POPC Membrane Bilayers BALAZS �� JOJART,1,2 �� �� TAMAS �� A. MARTINEK3 1Department of Chemistry and Chemical Informatics, Faculty of Education, University of Szeged, Boldogasszony sgt. 6, H-6725 Szeged, Hungary 2Institute of Pharmacodynamics and Biopharmacy, University of Szeged, Eotvos �� �� u. 6, H-6720 Szeged, Hungary 3Institute of Pharmaceutical Chemistry, University of Szeged, Eotvos �� �� u. 6, H-6720 Szeged, Hungary Received 7 January 2007 Accepted 7 March 2007 DOI 10.1002/jcc.20748 Published online 12 April 2007 in Wiley InterScience (www.interscience.wiley.com). Abstract: The aim of this work was to answer the question of whether the general amber force field (GAFF) is good enough to simulate fully hydrated POPC membrane bilayers. The test system contained 128 POPC and 2985 TIP3P water molecules. The equilibration was carried out in a nonarbitrary manner to reach the stable liquid-crystal- line phase. The simulations were performed by using particle mesh Ewald electrostatics implemented in molecular dynamics packages Amber8 (NPT ensembles) and NAMD2 (NP T ensembles). The computational results were assessed against the following experimental membrane properties: (i) area per lipid, (ii) area compressibility modu- lus, (iii) order parameter, (iv) gauche conformations per acyl chain, (v) lateral diffusion coefficients, (vi) electron density profile, and (vii) bound water at the lipid/water interface. The analyses revealed that the tested force field combination approximates the experimental values at an unexpectedly good level when the NP T ensemble is applied with a surface tension of 60 mN m 1 per bilayer. It is concluded that the GAFF/TIP3P combination can be utilized for aqueous membrane bilayer simulations, as it provides acceptable accuracy for biomolecular modeling. q 2007 Wiley Periodicals, Inc. J Comput Chem 28: 2051���2058, 2007 Key words: AMBER GAFF POPC lipid bilayer membrane simulation Introduction An understanding of the dynamic behavior of phospholipid bilayers at an atomistic level is of fundamental importance.1 The complexity involved in the fine balance between the order gener- ated by the self-assembly and the inherent disorder of fluidity renders molecular modeling an indispensable tool for the study of these systems.2,3 The pharmaceutical relevance of membrane proteins and the exponentially increasing number of experimen- tally resolved structures compels the medicinal chemistry com- munity to strive to model complex ligand���protein���lipid systems. Fortunately, the well-established methodological improvements (e.g., the application of periodic boundary conditions together with particle mesh Ewald electrostatics4) implemented in popu- lar modeling software packages and the increasing availability of high-performance computing resources facilitate the molecu- lar dynamic studies of fully solvated membrane protein systems even outside the leading laboratories. It is essential in the mod- eling of ligand���protein���lipid complexes that the force field pa- rameters for all the constituent chemical entities should be con- sistent. Internal consistency for the protein, lipid, and solvent has been achieved for only a few biomolecular force fields by fine-tuning the potential energy parameters: CHARMM,5 lipid- optimized GROMOS,6,7 and OPLS.8 The transferablility of these force fields to small organic molecules is limited, and explicit parametrization of a drug-like ligand requires special methodol- ogy and effort from experts in the field.9 Unfortunately, the pro- tocols are sometimes not fully described in the publications.10,11 Another situation arises when a general force field is available for drug-like molecules, such as the general amber force field (GAFF),12 which is consistent with the protein force field, but the lipid-optimized parameters are missing. As the literature Contract/grant sponsor: Gedeon Richter Pharmaceuticals, Hungary Contract/grant sponsor: HPC U-Szeged contract/grant number: ALAP4- 00092/2005 Contract/grant sponsor: Janos �� Bolyai Scholarship, Hungarian Academy of Sciences Correspondence to: T. A. Martinek e-mail: martinek@pharm. u-szeged.hu q 2007 Wiley Periodicals, Inc.

reveals, there is not a straightforward way to set up a force field for nonpeptidic ligand���protein���lipid systems within a single functional framework. On the other hand, the impetus to perform the pharmaceutically relevant modeling of receptor-nonpeptidic ligand complexes in aqueous membrane environment has resulted in publications where GAFF was utilized to model the lipid environment without tests of its accuracy.13���15 In the present work, our goal was to elucidate the real per- formance of the GAFF combined with RESP charges16 and with the popular TIP3P water force field17 in the modeling of an aqueous lipid bilayer. The studied lipid bilayer contained POPC molecules as frequently used in protein-membrane modeling. Moreover, the unsaturated carbons in the oleyl chain facilitate testing of the corresponding force field parameters and sufficient experimental data have been accumulated, with which the mod- eled observables can be compared. Methods Force Field Parameters and RESP Charge Calculations The GAFF atom types were determined by using the module Antechamber in Amber8 package with default setup. For the charge calculations, 12 diverse conformations were selected from the liquid crystal structure of the POPC bilayer published by Heller et al. (Fig. 1).18 The conformations were minimized in a consecutive manner: (i) molecular mechanics minimization with CHARMM2719 force field implemented in Chemical Computing Group���s MOE20 until the gradient was 0.0001 (ii) semiempirical minimization with the PM3 by using the Gaussian0321 program package and (iii) ab initio optimization at the HF/3-21G level. The molecular electrostatic potential was generated at the HF/6-31G* level. The atomic charges were fitted to the molecular electrostatic potential with the aid of the module Resp in Amber8 program package.22 During the charge fitting, the following restraints were used: (i) charge of ��1 for the N(CH3)3 group (ii) charge of 1 for the PO4 moiety (iii) in the alkyl chains, both the car- bon atoms and the hydrogen atoms were treated as equivalent and (iv) charge of 0.2 was ascribed to the hydrogens in the N(CH3)3 moiety. Initial Structure and Equilibration The system consisting of 128 POPC molecules was prepared by using the fully ordered POPC bilayer structure published by Heller et al.18 The system was rotated in order to set the z axis normal to the bilayer, and was soaked with TIP3P water mole- cules. The number of water molecules was 2985 (23.32 water molecules per lipid), and the system therefore consisted of 26,107 atoms. The equilibration was carried out by utilizing the nonarbitrary protocol in the literature.23 The method for the membrane equilibration was as follows: (i) the initial structure was minimized in 5000 steps with fixed POPC (the force con- stant was set to 500 kcal 1 mol 1 A �� 2) (ii) full minimization in 10,000 steps (iii) heating from 0 to 310 K in 50 ps with restrained POPC (the force constant was set to 10 kcal 1 mol 1 A �� 2) (iv) equilibration at 310 K for 20 ps (v) heating from 310 to 510 K in 10 ps (vi) equilibration at 510 K for 80 ps (vii) cooling from 510 to 310 K in 10 ps and (viii) simulation at 310 K for 20 ns. The molecular dynamics simulations were performed with two program packages: Amber821 and NAMD2.24 Periodic boundary conditions were applied in three dimensions, and the electrostatic interactions were calculated via the particle mesh Ewald method, and a 10-A �� cut-off was used for the van der Waals interactions. The bonds between the heavy atoms and hydrogens were constrained with the SHAKE algorithm, and a time step of 1 fs was applied in the heating and crushing steps, and a value of 2 fs in the equilibration phase. The pressure was controlled anisotropically in each case. In the NPT Langevin molecular dynamics simulations, the volume compressibility of the system was set to 44.6 10 6 bar 1 and 150 10 6 bar 1, which correspond to the volume compressibilities of water and proteins, respectively,25,26 The collision frequency was 1.0 ps 1. The pressure was maintained at 1 bar, and the pressure coupling constant was set to 1 ps. In the NP T Langevin molecular dynamics simulations (NAMD2), surface tension ( ) values of 0, 45, and 60 mN m 1 per bilayer was used. The pressure was again maintained at 1 bar. A sum- mary of the simulations performed is presented in Table1. The calculations were performed either on eight dual AMD Opteron processors, running the Debian GNU/Linux 3.1 (Sarge) amd64 version, or on an SGI Altrix 3000 allocating 16 CPU/ job 1-ns molecular dynamic simulation took *24 h of proces- sor time. Figure 1. Overlay of the 12 conformations utilized for the multiple conformer RESP charge calculation. 2052 Jojart �� �� and Martinek ��� Vol. 28, No. 12 ��� Journal of Computational Chemistry Journal of Computational Chemistry DOI 10.1002/jcc