Comparison of Biomolecular Force Fields for Alkanethiol Self

Subscriber access provided by READING UNIV

Article

Comparison of Biomolecular Force Fields for Alkanethiol Self-Assembled Monolayer Simulations Pratiti Bhadra, and Shirley Weng In Siu J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b08092 • Publication Date (Web): 07 Nov 2017 Downloaded from http://pubs.acs.org on November 10, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Comparison of Biomolecular Force Fields for Alkanethiol Self-Assembled Monolayer Simulations Pratiti Bhadra and Shirley W. I. Siu* Department of Computer and Information Science University of Macau Avenida da Universidade, Taipa, Macau

Corresponding author: [email protected] Phone number: +853-8822-4452

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 42

Abstract

Reasonable all-atom or united-atom biomolecular force fields have been developed to represent the properties of proteins and lipid membranes in molecular dynamics simulations. However, since they have not been parametrized for self-assembled monolayers (SAMs), their utility in simulating SAMs and protein-SAM systems has not been confirmed. Here, we compared six popular biomolecular force fields, Lipid14, GAFF, L-OPLS, CHARMM36, Slipids and GROMOS54a7, to simulate alkanethiol SAMs of short to long chains (C10-C18). Our results show that none of these force fields reproduce the chain length dependence of the tilt angle, tilt direction and twist angle. Although the droplet contact angles on SAMs are well represented by CHARMM36, L-OPLS, Slipids and GROMOS54a7, only GAFF and Lipid14 yield phase transition temperatures that are reasonably close to the experimental values. Overall, our comprehensive comparison suggests that GAFF and Lipid14 are better choices for SAM simulations; further improvements in the force field parameters for SAMs are required. (152 words)


2

Page 3 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Techniques to control the properties of solid surfaces, such as wetting, lubrication, corrosion and catalysis, have attracted attention in recent years for the development of nanotechnologies1. One effective way to achieve controlled and tunable properties of surfaces is by modifying them with self-assembled monolayers (SAMs). SAMs of organic molecules are well-ordered assemblies that form spontaneously on a substrate during the deposition process. Each constituent molecule, the ligand, is usually a linear structure that consists of three components: a head group (such as a thiol (-SH), a silane (-SiH4), or a phosphonate (C-PO(OH)2)), a spacer group (such as an alkyl chain) and a terminal functional group (such as -OH, -CH3, -NH2, -COOH, -CN, and -CF3). While the choice of the head group depends on the type of substrate material to which it needs to attach, the functional group depends on the desired wetting and interfacial properties of the modulated surface2. Without a thorough understanding of the SAM surface at the molecular level, the surface design can only be approached by experimental trial and error, which is costly and inefficient.

The molecular dynamics (MD) method has been proven to be a powerful tool for studying biomolecules at a great level of detail. At the heart of MD are the functional forms and parameters (referred to as force field) that describe interactions between atoms of different types in the simulation system. Empirical force fields have been actively developed in the research community, and parameters have been constantly revised to improve the accuracy. Popular MD force fields, such as AMBER, CHARMM, GROMOS, and OPLS, have been parameterized for proteins, nucleic acids and lipids. Additional force fields have been designed for small organic molecules (GAFF from AMBER, CGenFF from CHARMM) and non-biological macromolecules (such as DL_POLY and COMPASS).

Compared to proteins and lipids, the development of force fields for SAMs has been less active. The usual method of simulating SAMs is to modify and combine the parameters of existing MD force fields. ACS Paragon Plus Environment

3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 42

For example, Collier et al.3 and Szefczyk et al.4 simulated SAMs by adopting atom types and partial charges from selected amino acid side chains with matching atomic arrangements. Vellore et al. simulated SAMs using parameters from CGenFF and the rest of the system (protein and peptide) using standard CHARMM22 parameters5. Since the molecular structure of a SAM ligand highly resembles the structure of a lipid tail, a recent study by Peng et al.6 took the lipid parameters from CHARMM27 (lipid force field) and derived parameters for SAM ligands. A pioneering work for reparameterizing the force field specifically for SAM-protein interactions was performed by Abramyan et al., who tuned the CHARMM22 parameters for a high-density polyethylene (HDPE) SAM to reproduce the experimental protein adsorption free energies7.

Force field comparison studies provide an impartial opinion regarding the performance of a force field in representing molecular structures and behavior in MD simulations with respect to other force fields. These studies have been routinely conducted for protein and lipid force fields (such as in Ref. 8-10). For example, the latest work by Pluhackova et al. tested how four recent atomistic force fields, namely, CHARMM36, GROMOS54a7, Lipid14, and Slipids, represented four models of biomembranes10. They found that all of the force fields provided satisfactory results for the overall structural characteristics of phospholipid bilayers, although each had problems reproducing different aspects of the membrane. For SAMs, only one force field comparison study by Collier and co-workers has been reported to date. In their work, they compared the performances of three of the primarily used all-atomic force fields at that time (CHARMM22, AMBER94, and OPLS-AA) for experimentally reproducing the observed conformational behavior of LK peptides on SAM surface adsorption3.

In view of the recent improvements in MD force fields, in particular lipid force fields, the objective of this work was to compare the suitability of these force fields in simulating SAMs. We selected five protein-compatible lipid force fields, i.e., the all-atom CHARMM3611, Lipid1412, Slipids13, L-OPLS14, ACS Paragon Plus Environment

4

Page 5 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the united-atom GROMOS54a715, and one organic molecular force field, GAFF16. Our comparative study focused on the molecular structure, wettability and thermal stability of hydrophobic alkanethiol SAMs (SH-(CH2)n-CH3) of various chain lengths on gold. The molecular structure of a SAM was characterized by the tilt angle, tilt direction and twist angle of the ligand molecules, and the wetting property of the monolayer was characterized by the droplet contact angle. To measure the thermal stability of a SAM, we simulated the solid-to-liquid phase transition by simulated annealing and computed the phase transition temperature. All simulated structural and dynamical properties were compared against experimental observations of SAMs with moderate to long chain ligands.


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 42

2. Methods 2.1 Molecular Models SAM systems: Five SAM systems with different thicknesses were constructed from alkanethiolates (SH-(CH2)n-1-CH3 where n = 10, 12, 14, 16, 18). Each monolayer consisted of 90 (9 × 10) thiol ligands with initial configurations that were all trans and with their hydrocarbon chain axis normal to the surface. The dimensions of the SAMs were 39.05 Å × 42.72 Å × Z Å, where Z ranged from 13.90 Å to 23.75 Å (see Supporting Information Table S1). The Au(111) substrate underneath the monolayer was modelled as a closely packed cubic face-centered lattice structure with a lattice constant of 4.08 Å and a nearest neighbor distance of 2.88 Å, according to experiments2, 17. In total, four layers of gold atoms (1080 atoms) were explicitly modelled with a vertical separation of 2.44 Å from the monolayer. Finally, the SAM surface was solvated with 1,850 water molecules, for an approximate 30 Å water layer above the surface (see Figure 1). The total number of atoms in the model systems ranged from 9,510 to 11,700 atoms. To mimic the experimental condition, a thick vacuum slab was included above the water phase. Thus, the total box length in the Z dimension is 2 times the height of the water-SAM system

Droplet-SAM systems: To investigate the wetting property of the SAMs, water droplet simulations were performed. The monolayer and gold substrate with 360 (18 × 20) ligands and 4320 Au atoms were arranged in the same way as described above. We used a droplet of 1800 water molecules which is consistent with previous wetting simulation studies18-19. The droplet was initially equilibrated in vacuum and was placed at 5 Å above the monolayer. The radius of the droplet was ~21 Å, and the XY dimensions of the monolayer were 82.27 Å × 89.50 Å (see Supporting Information Figure S1).

All simulation models were constructed by the Packmol tool20, Visual Molecular Dynamics program (VMD)21 and in-house developed scripts.


6

Page 7 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: Snapshot of the 1-Octadecanethiol SAM: initial structure (left) and equilibrated structure after 15 ns of simulation time (right). The system consisted of three main components (from top to bottom): the water layer, the alkanethiolate monolayer, and the gold substrate.

2.2 Simulation Details Force fields: In the context of a MD simulation of biomolecules, the term force fields refers to the combination of mathematical formulas and associated parameters that calculate the potential energy of the system as a function of atomic coordinates. In this study, we focus on the widely used biomolecular force fields CHARMM3611, Lipid1412, Slipids13, GROMOS54a715, and L-OPLS14 and the general AMBER GAFF16 to represent the properties of the SAMs. Lipid14 and Slipids were integrated into Amber ff14SB22, and L-OPLS was integrated into OPLS-AA23. Different water models were selected to use with different force fields based on their original parameterization or better performance reported in previous MD studies: TIP3P24 with Lipid14, modified TIP3P25 (also named TIPS3P, in which the vdW parameters of the water hydrogen atom are non-zero) with CHARMM36 and Slipids, SPC26 with GROMOS54a7, and TIP3P-MOD27 with L-OPLS. For the sulfur atom at the substrate-facing terminal of


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 42

a SAM ligand, the force field parameters were taken from the CYS residue of the respective protein force field. We employed the Au parameters (ε = 0.039 kcal mol – 1 and σ = 2.934 Å for the LJ 6-12 potential) of Zhao et al.28 because they showed better performance than did the other Au parameters (e.g., in Ref. 29-32

). For completeness, we listed the selected atom types for a SAM ligand in all six force fields in the

Supporting Information Table S2.

Simulation conditions: All systems were simulated using periodic boundary conditions in all three dimensions. The particle-based cutoff method (Verlet) was employed for neighboring searching. Longrange electrostatics were handled by the particle-mesh Ewald (PME) method. Non-bonded cutoff parameters for CHARMM, AMBER (Lipid14 and Slipids), and GROMOS were referenced from Pluhackova et al.10, and for L-OPLS, they were referenced from Siu et al.14; the same parameters used for Lipid14 were used for GAFF (see Supporting Information Table S3 for a list of the parameter values). To simulate the chemisorbed SAM on gold, ligands were placed at the three-fold hollow sites of the Au lattice, and S atoms were positional restrained with a force constant of 1000 kJ mol-1 nm-2 throughout the simulation. The substrate Au atoms were kept fixed. Simulations were performed at 300 K in a canonical ensemble (NVT) using GROMACS 5.0.733. A time step of 2 fs was used for all equilibration and production dynamics, in which data were collected at 1-ps intervals.

Simulated annealing for melting simulations: To study the phase behavior of SAMs from different force fields, we started from an equilibrated structure at 300 K and gradually increased the temperature to 420 K in 6 ns, with a rate constant of +1 K ps-1. Because the starting temperature of monolayer desorption from a gold substrate is below 400 K34-35, heating the system to 420 K should be sufficient to capture the melting process of SAMs. To obtain statistics, each simulation was repeated three times with different starting velocities.


8

Page 9 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.3 Analysis/properties Structural parameters of the SAMs: As shown in Figure 2, the orientation of a well-ordered SAM ligand can be described by three angles: the tilt angle, the tilt direction, and the twist angle. The tilt angle (α) is the tilt of the chain, which is defined as the angle between the hydrocarbon axis and the surface normal (Z-axis). The tilt direction (δ) of the chain, which is defined with respect to the next nearest neighbor (NNN) direction of the sulfur atom (SNNN) (see Supporting Information Figure S2), is computed as the angle between the projection of the hydrocarbon axis on the surface plane and the S-SNNN vector. Finally, the twist angle (χ), which is the rotation of the chain about its long axis, is computed as the angle between the carbon plane (composed of the first three atoms of the ligand: S, Cα, and Cβ) and the plane determined by the hydrocarbon axis and the surface normal.

Figure 2: Schematic representation of a SAM ligand and three geometric parameters describing its orientation on the substrate surface: tilt angle (α), tilt direction (δ), and twist angle (χ). Three atoms are important in computing the geometric parameters: Sulfur (S in yellow), the first carbon (Cα in red), and the second carbon (Cβ in purple). The hydrocarbon axis (also known as the backbone) is given by the vector S-Cβ, and the hydrocarbon plane (in deep red) is defined by S, Cα, and Cβ. The plane through the surface normal and hydrocarbon axis is indicated by the light red rectangular plane. ACS Paragon Plus Environment

9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 42

Gauche defect: The orderness or disorderness of a SAM can be characterized by its percentage of gauche defects along the hydrocarbon chains. While the alkanethiolate chains of SAMs at room temperature are predominantly all trans in arrangement, at higher temperatures, different amounts of gauche defects can be observed. In this study, the conformation is considered gauche when the torsion angle of four consecutive carbon atoms is less than 120° in the positive or negative direction36.

Melting temperature: The phase transition temperature, Tp, is determined as the onset temperature when the melting process begins. It is obtained by fitting the hyperbolic tangent function in Eq (1) to the percentage of gauche defects as a function of temperature (T):

2

2

1

where Tm is the middle temperature of the transition period, C is one-half of the temperature range in the transition period, gu and gl are upper and lower values of the percentage of gauche defects. After Tm and C were determined, we computed the phase transition temperature as Tp = (Tm – C). The curve fitting tool in the MATLAB software was used for the fitting.

Wettability of the SAMs: The wettability of a surface is its propensity to form a common interface with the contacting liquid. It can be measured as the contact angle of the liquid (droplet) at the threephase contact point of a sessile drop profile. To compute the contact angle of a droplet from the MD trajectory, we employed the geometric method proposed in Refs. 37-38. Assuming that the droplet adopts a hemispherical shape, the contact angle can be estimated from the height of the hemisphere and the droplet contact surface area. Details of the calculation are provided in Supporting Information Section II.


10

Page 11 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3. Results 3.1 Structural properties The ability to reproduce the experimental structural properties of SAMs was the first test on the six biomolecular force fields. In equilibrium, a SAM is a highly ordered film where chains are packed to minimize the surface free energy. The orientation of the SAM ligands under different conditions (such as the chain length and substrate material) can be characterized by three angles: the tilt angle (α), the tilt direction (δ), and the twist angle (χ). In this section, we compare the structural parameters for short- to long-chain SAMs from different force fields to the experimental values. Our simulation systems equilibrated rather quickly within 1 ns. Nevertheless, we extended the runs to 5 ns equilibration and 10 ns data production. A test of 100-ns long simulation on the GAFF C12 system was carried out and we found no significant differences in the results. Property means were calculated by ensemble averaging and the standard error was estimated by block averaging.

3.1.1 Tilt Angle According to the grazing incidence X-ray diffraction (XRD) experiment by Fenter and co-workers39, both α and δ of the SAM’s hydrocarbon chains on Au(111) displayed systematic dependence on chain length. Two distinct regimes were identified: the short regime, where n ≤ 14, 18°, and the long regime, where n ≥ 16, = ~30.5° and 4°
15)

GAFF

Short-Chain Regime (n

Comparison of Biomolecular Force Fields for Alkanethiol Self

Recommend Documents