Simulation of Mixed Self-Assembled Monolayers on Gold: Effect of

Article Cite This: J. Phys. Chem. B 2018, 122, 7699−7710

pubs.acs.org/JPCB

Simulation of Mixed Self-Assembled Monolayers on Gold: Effect of Terminal Alkyl Anchor Chain and Monolayer Composition Eric Schulze and Matthias Stein* Molecular Simulations and Design Group, Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstr. 1, 39106 Magdeburg, Germany

J. Phys. Chem. B 2018.122:7699-7710. Downloaded from pubs.acs.org by UNIV OF SOUTH DAKOTA on 08/29/18. For personal use only.

S Supporting Information *

ABSTRACT: Self-assembling monolayers provide a reproducible synthetic microenvironment for tethering lipid bilayers to incorporate proteins and lay the ground for numerous applications in nanotechnology and biomedical engineering. Although the structure of single-component monolayers is well investigated, there is far less insight into the molecular behavior at the interface of mixed monolayers at different mole fractions. Here, we present and apply a novel procedure to simulate and analyze multicomponent self-assemblies of alkanethiols over a wide range of mole concentrations of anchoring compounds. In particular, the structural features of monolayers consisting of a matrix compound and either a short (C8) or a long (C16) anchor compound on Au(111)-like surfaces were investigated first using coarse-grained and subsequently full-atomistic molecular dynamics simulations. Different scenarios of spatial distributions (random vs clustering) of anchoring molecules on flat surfaces were probed. The results of the simulations are in excellent agreement with the experimental data from ellipsometry and infrared reflection absorption spectroscopy. For short anchoring molecules, a random spatial distribution in the matrix is obtained. At low, experimentally relevant anchor compound mole fractions < 0.1, only for long-chain (C16)-terminal alkyls, phase segregation and self-association of the anchoring molecules can be observed, which are also seen in experiment.

■

INTRODUCTION The phenomenon of an adsorbing (self-assembling) monolayer (SAM) onto a metal surface was first observed more than 80 years ago by Bigelow and co-workers,1 and it later regained interest as a model system to investigate the fundamentals of intermolecular interaction and adsorption.2 The ability to modify both the head and tail groups of the layering molecules makes SAMs excellent systems to probe the numerous competing effects such as hydrophilicity versus hydrophobicity and order versus disorder. Apart from silanes or fatty acid derivatives on hydroxylated surfaces, a large number of publications are devoted to the studies of organothiol/disulfide assemblies on semiconductor or metal surfaces, in particular gold (for relevant reviews, see refs 2 or 3). The high degree of control and reproducibility of SAM formation make them an ideal model system to design a controllable microenvironment in nanotechnology and biology. For example, SAMs with properly designed terminal groups (herein referred to as “anchor”’ groups) serve as an excellent starting point for the development of tethered bilayer lipid membranes (tBLMs).4 Such tethered bilayers prepared by the fusion of liposomes5 can be regarded as a first-order “synthetic” model of a biological cell membrane6 as they can provide a close-to-native and protective environment for membrane-associated proteins.7 For a review article about © 2018 American Chemical Society

the structure and function of supported BLMs, see ref 8. Highly reproducible and standardizable protocols are available, which turn SAM−tBLMs9 into an attractive platform for the structural and functional characterization of membrane proteins and also for the screening of drug candidates targeting these proteins.10,11 For example, tBLMs have proven to be useful for monitoring the dynamics of a ternary cytokine receptor protein complex.12 In addition, SAMs and tBLMs on solid support are also well-suited for an extensive characterization using established surface analytical techniques such as surface plasmon resonance spectroscopy, X-ray and neutron reflectometry, contact angle goniometry, infrared (IR) spectroscopy, fluorescence microscopy, scanning probe microscopy, atomic force microscopy, and electrochemistry, to mention only a few.13,14 Organothiol SAMs are typically prepared by immersing a gold substrate overnight in a dilute ethanolic solution of the alkanethiol compound. The alkanethiol can be functionalized with soluble polymer/oligomer moieties, for example, oligo ethylene glycols (OEGs), and such SAMs have been extensively used for fundamental protein adsorption studReceived: May 28, 2018 Revised: July 18, 2018 Published: July 20, 2018 7699

DOI: 10.1021/acs.jpcb.8b05075 J. Phys. Chem. B 2018, 122, 7699−7710

Article

The Journal of Physical Chemistry B ies.15−17 OEG−alkanethiols (as a matrix compound) also have been proposed to function as a soft cushion for tBLMs by Lee et al.18 In this particular case, the matrix compound is mixed with a longer intercalating molecule that serves as an anchor for the lipid bilayer (see Figure 1).

molecular level cannot be reported. A detailed interpretation of the molecular conformations to rationalize the spectral differences, though, remains elusive. Here, molecular dynamics (MD) simulations can provide additional insight into the structure and dynamics of SAMs in atomic detail. In the past, it led to a clearer understanding of the architecture of the monolayers of alkanethiol chains with different terminal groups.21,22 Also, large-scale MD simulations of the alkanethiol self-assembled monolayers have been performed to study the effects of temperature and packing density on the structural parameters.23 Sufficient sampling and convergence to equilibrium configurations are difficult to achieve because of slow molecular reorientations, and for the reproduction of lateral diffusion, advanced conformational sampling algorithms are required.24,25 The accuracy of the different force fields varies to a great extent, especially in the reproduction of chain length-dependent tilt and twist angles.26 Nevertheless, MD is the appropriate tool to study the SAM structure and dynamics at high spatial and temporal resolution. To overcome the complexity and the necessity of long equilibration periods, especially in systems involving membranes or vesicles, the Martini coarse-grained (CG) force field has been of extensive use recently,27 and is thoroughly introduced in the Methods section. Besides the mentioned application on phospholipid bilayers28,29 and polymers,30 PEGylated lipids31 and alkanethiol-covered nanoparticles32 demonstrate the feasibility of this CG force field for molecular systems comparable to ours. Although most of the previous MD studies focused on single-component matrix molecule monolayers, we here present a new approach for the systematic simulation of a multitude of mixed monolayers at different mole fractions of different anchoring molecules using various representations and a concise transformation between them. Our simulations are most helpful to complement the experiments because they are able to probe mixed concentrations and mole fractions that are experimentally hard or impossible to access. In this article, we address for the first time the structure and dynamics of mole fraction-dependent orientational and conformational transitions of the alkyl anchor compounds in a matrix environment attached to a gold (111)-like surface. We here show that structural parameters such as lattice constant and tilt angle are almost independent of the type of anchoring molecule. The anchor alkyl-chain orientation, however, is critically controlled by the composition and its chemical nature. Although the short-chain C8 anchor molecules adopt a random, disordered conformation at low anchor densities and transition into a more ordered conformation with increasing anchor density, the C16 anchors already adopt highly ordered conformations at low anchor densities, which reoccur in the monolayers of higher anchor density. The results of our largescale, two-step, multiscale procedure are in excellent agreement with the recent experiments by Lee et al. and explain the observed spectral features. Finally, our newly composed modeling procedure forms the basis for future computational investigations on SAM−tLBMs and SAM-tethered lipid vesicles of various anchoring molecules and monolayer compositions.

Figure 1. Schematic representation and nomenclature of the modules of gold−alkanethiol SAM functionalized with OEG and covalently linked anchor compounds C8 and C16. The matrix compound is colored gray, short C8 anchoring compound green, and the long C16 anchoring compound blue.

In a recent experimental work of Lee et al., mixed SAMs and mixed SAM−tBLMs have been prepared from a solution of matrix and anchoring compounds. In their study, they synthesized an alkanethiol matrix compound (HS(CH2)15CONH(CH2)2OH) and two different anchoring compounds with the formulas HS(CH 2 ) 1 5 CONH((CH2)2O)6CH2CONH(CD2)7CD3 for the short C8 anchor and HS(CH2)15CONH((CH2)2O)6CH2CONH(CD2)15CD3, for the longer C16 anchor (see Figure 2). A detailed

Figure 2. Chemical structures and CG mapping schemes of the matrix and C8 and C16 anchoring compounds. The colors represent the Martini bead types. Yellow: C5, gray: C1, blue: P5, red: SN0, purple: SP2. The −(CH2−)15 matrix alkyl chain is reduced in the anchor compounds for clarity of display.

introduction to different molecules, their synthesis, and characteristics can be found in ref 19. The anchor alkyl chains are fully deuterated (d-alkyl) to make their assignment from IR spectroscopy unambiguous and distinguishable from the alkyl chains in the matrix region. A series of SAMs with increasing relative concentration of the anchoring compound were prepared from solutions with the following mole fractions: xA,Sol = 0.01, 0.05, 0.1, 0.15, 0.30, 0.5, and 1.0 and subsequently characterized using ellipsometry, contact angle goniometry, and IR reflection adsorption spectroscopy (IRRAS).18,20 With IRRAS, differences in phase behavior for the two model systems were found. The indications of phase segregation and domain formation were more prominent in the SAMs prepared from the longer compound, but structural insight at the

■

METHODS Coarse-Grained MD Setup. Martini is a CG force field with a mapping of an average of four heavy atoms into single interaction sites, also simply called beads. Here, the mapping of 7700


Article

The Journal of Physical Chemistry B atoms belonging to chemical groups into CG interaction sites was performed with the introduction of the original Martini force field,28,29 its extensions to proteins33 and polyethyleneglycols (PEGs).30 In fact, as the investigated molecular constituents are unbranched and their main components (long-chain alkanes and PEGs) are extensively studied, the bead mapping was straightforward. Our mapping schemes (Figure 2) employ C5 beads for the thioethyl head groups and C1 beads for each n-butyl group. The acetamide groups, linking the alkyl chains and OEG, are mapped into P5 beads. The oxyethylene groups are mapped into SN0 beads. In case of terminal ethanol groups, SP2 beads are used instead. The nonbonded parameters are determined by the bead types and were employed with only two minor changes: According to the authors of the Martini PEG extension, the SN0−SN0 and SN0−P4 Lennard Jones interactions were adjusted to the level of Nda particles.30 More specifically, the interactions were set to be more attractive than in the original Martini force-field parametrization. The bonded parameters for the alkyl and OEG chains were readily applied from previously published studies.28,30 Here, only the amide linkage between the PEG and alkyl chain was parametrized. This was achieved by fitting the results from CG MD (bond length and bondangle distributions) to united-atom simulations with the GROMOS 53A6 force field34 (see the Supporting Information). We note that the dihedral potential terms were omitted from the PEG model for the sake of stability when using larger integration time steps. Because of the short polymeric chain length of the PEG, this simplification allows integration time steps of 0.3 ps with only a marginal loss of accuracy. As a generally accepted parameter set for a polarizable gold surface is not available in Martini, we used a flat surface model, which was successfully applied earlier.22 Therefore, the thiol group beads are restrained to a plane via harmonic potentials with a force constant of 5000 kJ/(mol nm2). Initial Configurations of the Monolayers. The monolayers were built with the Python tool insane.py, initially developed by Wassennar et al.35 It was originally designed for an automated modeling of differently composed CG lipid bilayer membranes and, for example, used for the simulation of the plasma membrane36 and the early endosome membrane.36,37 The layering molecules are placed in a simplified, extended conformation onto a rectangular grid. The user has options to control the box size, the number of molecules and their relative abundance, choice of solvent, ionization states, and several more options. In case of multiple component systems, the spatial positioning is random following a uniform distribution. Our modifications of the script involve the introduction of the matrix and both C8- and C16-anchoring compounds as well as the molecular positioning onto a hexagonal grid with defined nearest neighbor distances (see Figure 3). An isotropic, hexagonal lattice (√3 × √3 R30°) is considered sufficient as the contribution of the c(4 + 2) superlattice at room temperature is negligible.38 Apart from the uniform anchor distribution, an additional optional “clustered” positioning of anchoring molecules was implemented. We here only address the spatial and orientational distribution of one type of matrix compound with one of the two types of anchoring molecules. An overview of all simulated systems is given in Table 1. A more detailed description of the simulated systems and the placing algorithm for the clustered monolayers can be found in the Supporting Information. In all instances, the molecules

Figure 3. Details of initial positioning of random uniform (A) and clustered (B) monolayer coverages (here at xA = 0.1) of C8 and C16 anchoring molecules. Top: side view into the simulation box; bottom: top view of the initial positioning in a hexagonal grid. The water beads are only displayed for half of the box.

Table 1. List of All MD Simulations for Mixed Matrix and Different Anchor Molecule SAMs Systemsa CG (20 × 20 nm2, 50 ns) C8 xA

R

0.00 0.05 0.10 0.15 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00

● ● ● ● ● ● ● ● ● ● ● ● ●

AA (10 × 10 nm2, 10 ns)

C16 C

R ● ● ● ● ● ● ● ● ● ● ● ● ●

C8 C

R

C16 C

R

C

● ● ●

● ●

● ● ●

● ●

●

●

●

●

●

●

●

●

●

●

a

R stands for random spatial distribution, C for clustered spatial distribution. CG = coarse-grained and AA = fully atomistic MD simulations.

were placed into a hexagonal arrangement with an initial nearest neighbor distance of 5.4 Å (see Supporting Information Figure S4). This value was chosen as a compromise between the optimum reproduction of tilt angles and packing and computational stability. Moreover, the simulation boxes had an initial normal dimension of 13 nm and were solvated with polarizable Martini water.39 The water beads were placed by insane.py into the standard rectangular grid. Overlapping beads were removed. The monolayer systems were not ionized. The two initial configuration scenarios are presented in Figure 3. The extended monolayer modeling script is available upon request. CG Simulation Protocol. For all CG simulations, the currently published Martini simulation parameters for polarizable water were used40 and are briefly described in the Supporting Information. Initially, the systems were minimized with the steepest descent algorithm for 5000 steps, followed by a short NVT equilibration for 10 ps with a 0.01 ps time step. Initial velocities were generated according to the Maxwell− Boltzmann distribution at 298.15 K. NPT equilibration was employed for 10 ns with a 0.02 ps time step and semi-isotropic 7701


Article

The Journal of Physical Chemistry B Berendsen pressure coupling.41 The coupling constant was set to 5 ps−1 and lateral and normal compressibility to 0 and 4.5 × 10−5 bar−1, respectively. The reference coordinates for position restraints were scaled accordingly. The production simulation was run for 60 ns with a time step of 0.03 ps. Here, a Parrinello−Rahman barostat42 with a coupling constant of 12 ps−1 was employed. The first 10 ns of the production were discarded as additional equilibration and thus not considered in the trajectory analysis. Full Atomistic MD Refinement. Following the CG system setup, initialization, equilibration, and production, more accurate fully atomistic MD simulations were performed to refine the results from CG simulations. The resolution transformation of CG topologies to full atomistic detail was performed with backward.py, as published by Wassennar and Pluhackova.43 The full atomistic CHARMM3644,45 topologies were generated via the CHARMM-GUI,46 which performs analogy searches for the parameters in the CHARMM General Force Field (CGenFF).46,47 The maximum penalties were 32 for the partial charges and 41 for the bonded parameters, which are found in the amide linkage between the alkyl and OEG chains. The generated CGenFF topologies are available upon request. The whole refinement procedure is similar to the one published in the backward.py article and is accurately described in the Supporting Information. Production was performed for 20 ns with a 0.002 ps time step, where the first 10 ns were excluded from the analysis. Trajectory Analysis. The major structural features that describe the packing and orientation of the SAMs are the lattice parameter ⟨a⟩, the alkane thiol tilt angle ⟨θ⟩, the monolayer thickness ⟨d⟩, and the surface hydrophobicity ⟨η⟩ (see Table 2).

local thickness and the lateral thickness distribution, the simulation box was divided into 50 × 50 cells, in which each the normal direction difference between the bead with the maximum z-coordinate and the average thiol bead (from all thiols beads) z-coordinate was calculated. This value was averaged over the full production trajectory runs. Exemplary density profiles can be found in the Supporting Information. The surface hydrophobicity ⟨η⟩ is estimated from the fraction of the solvent-accessible areas of the hydrophobic moieties relative to the total surface-accessible area of the monolayer. The lower surface, as formed by the thiol groups, was not taken into account. The solvent-accessible areas were computed using the GROMACS tool gmx sasa. The van der Waals radii of the Martini beads were set to 0.26 nm for the standard ones and 0.23 nm for the small ones. The probe radius was left default at 0.14 nm. Additionally, in the full atomistic simulations, the C−C bond order parameters SCC were calculated. As in ref 49, they can be defined as SCC =

where φ is the angle of a C−C bond vector relative to the surface normal. In particular, for the matrix ordering, the C−C bonds of the alkanethiol part, for OEG the C−C bonds in the PEG chain, and for the anchor compound the C−C bonds of the anchor alkyl chain were considered. Usually, the time and ensemble average of a particular bond order parameter is used, but as we investigate the order parameters of multiple similar bonds (which leads to multimodal distributions), we refrain from averaging and instead show the bond order parameter distributions for selected chemical bonds separately. Torsion angles were computed from the atomistic coordinates. The angles between 150° and 180° as well as −150° and −180° are considered as trans. The angles between 50° and 150° and −50° and −150° are considered gauche. The helicity is computed as described in ref 50. It is the number of consecutive, same-directed O−C−C−O gauche torsions in the OEG chain divided by the number of all O−C−C−O torsions (5, in our case) and adopts values between 0 and 1. The hydrogen bonds were calculated with GROMACS gmx hbond. If not stated otherwise, all analyses were performed with Python 2.7.1. The scripts are available upon request. All visualizations and renderings were created with VMD V. 2.9.3.51

Table 2. Frequently Used Symbols and Abbreviations quantity

unit

description

xA xC8

anchor mole fraction of either compound anchor mole fraction of C8 compound

xC16

anchor mole fraction of C16 compound

⟨a⟩

Å

⟨θ⟩ ⟨d⟩ ⟨η⟩ SCC ⟨rtrans⟩ ⟨rhelix⟩

deg nm

3 cos2 φ − 1 2

lattice constant, i.e., averaged nearest neighbor S−S (thiol−thiol bead) pinning distance averaged collective matrix alkyl chain tilt angle monolayer thickness averaged surface hydrophobicity C−C bond order parameter average ratio of trans torsion angles average helicity

■

RESULTS AND DISCUSSION CG Simulations. To assure consistency and reproducibility of the newly generated CG model setup, three replicates each of 50 ns production time for 14 different anchor mole fractions xA of the C8 and C16 anchor compounds in the matrix compound were performed (see Table 1). The anchor compounds were randomly positioned according to a uniform probability distribution. The nearest neighbor pinning distance of the thiol beads ⟨a⟩, the tilt angle of the alkane thiol chain ⟨θ⟩, the monolayer thickness ⟨d⟩, and the surface hydrophobicity ⟨η⟩ were constantly monitored during the simulations. The averaged results of each 30 ns production together with the standard deviation from the three independent simulation runs are presented in Figure 4. The results are consistent and well-reproducible in the three replica simulations. The standard deviation of the three replicates is

The value of the lattice constant ⟨a⟩ is computed as the time and ensemble average of the mean of the six shortest lateral sulfur−sulfur (thiol−thiol bead) distances, taking the mirror images into account. The tilt angle ⟨θ⟩ is calculated as the time and ensemble average of the angle between the surface normal and the first principal axis of the matrix alkyl chain carbon atoms, or C1 beads, respectively. The monolayer thickness ⟨d⟩ is determined from the time-averaged density profiles in normal direction, as the distance between the intercepts of the profiles of the solvent and bilayer densities, that is, between the effective phase boundaries. This approach is used to make the simulation results directly comparable to the experiments from ellipsometry.48 The density profiles were computed with the GROMACS tool gmx density by dividing the box in 50 (CG) or 100 (AA) equal slabs in normal direction. To calculate the 7702


Article

The Journal of Physical Chemistry B

matrix, 49.4 ± 0.3 Å for the pure C8 anchor, and 58.6 ± 0.4 Å for the pure C16 monolayers, respectively. In the mixed monolayers, ⟨d⟩ increases only slightly for xA less than or equal to 0.2, independent of the anchor chain length. Above xA of 0.3, ⟨d⟩ however increases monotonically with the anchor compound mole fraction and eventually levels out beyond xA of 0.8. The increase in thickness with xA is larger for the longer C16 anchor compound, which is expected because of the larger final layer thickness of 56 Å. The large increase of monolayer thickness ⟨d⟩ with xA between 0.2 and 0.8 is due to an orientational phase transition in which the anchor alkyl chains rearrange from a surface-parallel orientation at low mole fractions to an extended, densely packed perpendicular orientation at larger mole fractions (see Figure 5).

Figure 4. Lattice constant ⟨a⟩, tilt angle ⟨θ⟩, thickness ⟨d⟩, and surface hydrophobicity ⟨η⟩ of SAMS with C8 and C16 anchoring compounds as functions of molar fraction. Orange filled circles: C8 anchors; purple filled triangles: C16 anchors. The shaded area is the standard deviation from the mean of three replicate trajectories.

not visible to the eye. The lattice constant ⟨a⟩ as well as the tilt angle ⟨θ⟩ do not significantly depend on the anchor mole fraction xA, as they are the intrinsic properties of the matrix compound. The monolayer thickness ⟨d⟩ and the hydrophobicity ⟨η⟩ of the layer, however, do increase with xA, as the intermolecular interactions become more prominent. The dependency of the tilt angle ⟨θ⟩ on the lattice constant ⟨a⟩ is known from previous studies. It decreases from around 30° at a ∼5 Å spacing to 0° at 4.2 Å,52,53 as a result of the decrease in the intermolecular space. To reproduce the experimental tilt angles, ⟨a⟩ was initially modeled to 5.4 Å and restrained during the simulation by setting the lateral compressibility to 0 bar−1 (see the Methods section and the Supporting Information). Because of the structural rearrangements, it finally converges to 5.35 Å. Thus, a tilt angle of ∼24°, independent of the mole fraction and chain length of the anchor compounds, is obtained for the matrix alkyl chains. This is in good agreement with the experimental findings for the tilt angles of alkanethiols on gold (111) surfaces, ranging from 20 to 35°,2,54,55 depending on the experimental details. We note, however, that with the standard Martini bead definitions, it is not possible to exactly reproduce the experimental nearest neighbor sulfur−sulfur pinning distances of 4.97 Å for gold (111) surfaces,52,53 while still obtaining tilt angles very close to the experiment. Apparently, the van der Waals radii of the Martini beads are too large to allow a close alkyl-chain interlocking as observed in real alkanethiol monolayers. This has to be considered during the resolution transformation of simulated CG monolayers with tethered lipid bilayer membranes and/or vesicles. Here, a simple coordinate scaling by a factor of 0.93 during the backmapping protocol is used to resolve this issue. The simulations gave monolayer thicknesses of 23 Å for pure matrix molecules (xA = 0), 46 Å for the pure C8 anchor monolayer (xA(C8) = 1), and 56 Å for the pure C16 anchor monolayer (xA(C16) = 1) (see Figure 4). These results are in excellent agreement with the ellipsometric thickness measurements by Lee et al., who reported 24.6 ± 0.4 Å for the pure

Figure 5. Snapshots of CG monolayer MD simulations with (A) xC8 = 0.05 and (B) xC8 = 0.6. With higher xC8, both the anchor alkyl and the OEG chains adopt an extended conformation in vertical orientation. Matrix and anchor compounds are depicted as licorice, water beads in van der Waals representation. Only a small, representative excerpt of the entire box is shown.

This reorientation can be quantified and is thoroughly discussed in the section “Comparison to IR Reflection Absorption Spectra”. The monolayer thickness ⟨d⟩ of the mixed monolayers is more difficult to compare with the ellipsometric measurements because of the experimental inaccessibility of high mole fractions and uncertainties in the determination of the adsorbed anchor mole fraction. The simulations, however, give a very detailed picture of molecular assemblies that explains why the monolayer composition assessed by ellipsometry is most likely an underestimation, that is, not fully comparable between different anchor alkyl chain lengths. The surface hydrophobicity ⟨η⟩, on the other hand, increases almost linearly with xA and is only slightly smaller for the C8 than the C16 anchor SAMs. It takes xC8 of 1.0 to fully cover the surface, whereas it takes only 0.9 for xC16 because of its larger radius of gyration. Especially, this linearity is a strong indicator for the simulations to be in good agreement with experiments. Experimentally, the water contact angle is used as a metric for the anchor mole fraction. We note that the water contact angles have been successfully computed from MD simulations as, for example, in refs.56−58 The Martini polarizable water, however, does not form stable droplets in our simulation setup. Thus, we calculated the surface hydrophobicity directly from the solvent-accessible areas of hydrophilic and hydrophobic molecular parts. The fact that the surface hydrophobicity (as opposed to the monolayer thickness) is linear to xA makes the 7703


Article


a fully atomistic representation was performed after CG modeling and equilibration (see Figure 7). The full atomistic

monolayer anchor mole fraction, as calculated from the contact angles, more reliable. For the two-dimensional, lateral distribution of the monolayer thickness, we report a remarkable nonuniform distribution and the formation of self-aggregates (see Figure 6). From a top perspective, there is a clear lateral separation,

Figure 7. Snapshots of an xC16 = 0.1 monolayer before (A) and after (B) resolution transformation in van der Waals representation. The anchor compound carbon atoms (beads) are colored green.

systems were subsequently carefully minimized and equilibrated before the production simulations, as described in detail in the Supporting Information section, were performed. We ensured that the systems converge to the same orientations and conformations, independent of the number of CG steps preceding the resolution transformation. The properties obtained from the CG simulation are wellreproduced in the full atomistic simulation (Figure 8), which

Figure 6. Spatial distribution of the monolayer thickness (in Å) for increasing anchor mole fractions (from left to right) and the two anchor compounds. Top: C8, bottom: C16. The average thickness here appears smaller than that in Figure 3 by means of different definitions (interface area distance from density profile vs bead to bead center of mass distance).

which is attributed to the surface segregation of the anchor chains. We have to stress that this segregation is not the result of a phase separation of the individual molecules on the gold surface, but rather a consequence of self-aggregation of the dalkyl chains above the monolayer surface. Comparing the identical mole fractions of C8 and C16 anchors, we find larger aggregates for the C16 monolayers. This can easily be explained by the larger radius of gyration of the C16 anchor and thus its higher probability to interact with others. We note that the formed aggregates are stable during the entire 50 ns of production in all three replicates. Only single anchor molecules at the boundary regions shuttle between different aggregates. It can be speculated that the interactions of the superficial alkyl anchor aggregates are the driving forces for the slow phase segregation on the gold surface in a timescale of hours to days. As these interactions will clearly be stronger for longer alkyl chains (in identical anchor mole fraction), the hypothesis of phase segregation exclusively for the C16 compound is supported by the simulations. As the timescales of the SAM phase segregation are too long, even for CG simulation, initially clustered monolayers were modeled and analyzed. The corresponding results can be found in the section “Chemical Bond Order Parameters” and the following sections. To summarize, the CG model is able to yield accurate structural parameters of the C8 and C16 anchoring compounds in a matrix environment over a large range of mole fractions in a fast, automatable, and reproducible manner. We consider the CG model to be convenient for high-throughput modeling and equilibration of differently composed SAM and SAM−tBLM systems in future applications. Here, we want to point to a recent study by Hoiles et al.,59 who also applied CG MD to investigate the structural and dynamical parameters of tethered lipid membranes depending on the bilayer composition and tether density. However, their studies do not include the matrix of alkanethiols in which the anchoring/tethering molecules are embedded. Full Atomistic Simulations. To further refine the structural insight into anchoring molecules in a monolayer matrix, a resolution transformation from the four-atom bead to

Figure 8. Structural monolayer properties computed from full atomistic simulations after resolution transformation. Orange circles represent the C8 anchors and purple triangles the C16 anchors.

underlines the accuracy of the CG model, retrospectively. The lattice constant ⟨a⟩ is consistently 4.97 Å for the entire range of anchor mole fractions and in perfect agreement with the experimentally reported value (4.97 Å for the alkanethiol monolayers on gold (111)52,53). This confirms the accuracy of the applied f lat surface model. The matrix alkyl-chain tilt angle of all simulated monolayers is between 25° and 28°, which decreases monotonically with an increasing anchor mole fraction, with one exception: the tilt angle is the lowest for the xC16 = 0.5 monolayer. The most plausible explanation for such a decrease is that a smaller tilt angle allows more and stronger hydrophobic interactions between neighboring C16 anchor chains. This exceptional 7704


Article

The Journal of Physical Chemistry B behavior only at the xC16 = 0.5 monolayer and its absence in the C8 anchor compound further support this hypothesis. Because of the large surface aggregates of the anchor alkyl chains, as shown for the CG model in Figure 6, the crystal-like packing of the matrix alkyl chains is perturbed. The shorter C8 anchor chains do not exhibit such strong interactions, and the xC16 = 1.0 monolayer is too tight to introduce such perturbations. Such an effect is not seen in the CG model, in which the tilt angles are independent of the monolayer composition. In the AA simulation, the monolayer thickness ⟨d⟩ increases with xA as in the CG model. The experimental thickness values of pure monolayers are also well-reproduced. We find ⟨d⟩ = 25 Å for the pure matrix compound monolayer and ⟨d⟩ values of 51 and 62 Å for the pure C8 and C16 anchor compound monolayers, respectively. In comparison to the thickness values of 24.6 ± 0.4, 49.4 ± 0.3, and 58.6 ± 0.4 Å, as determined via ellipsometry by Lee et al., our full atomistic simulations only slightly overestimate the values, especially for the pure anchor compound monolayers. This is most likely because of a slight underestimation of helicity in the OEG chains (PEG in a helical conformation is about 20% shorter than that in a zigzag conformation17), as more thoroughly discussed in the sections below. Interestingly, in the full atomistic model, a significant onset of difference in ⟨d⟩ between the anchor compounds has already been detected at the mole fractions of around 0.2, which is slightly lower than that in the CG model. At even lower anchor densities, ⟨d⟩ remains constant and is independent of the anchor compound type. This observation further supports that the interaction and formation of small, bulk-like aggregates of the anchor alkyl chains occur at a lower mole fraction for the C16 than for the C8 anchoring compound. The surface hydrophobicity ⟨η⟩, as computed from the full atomistic simulations, is close to the CG model. Its value linearly correlates with the anchor mole fraction. For xA ≤ 0.2, ⟨η⟩ is identical for the C8 and C16 compounds. For larger anchor mole fractions of xA > 0.2, ⟨η⟩ of the C16 anchor composite monolayers is larger by an amount of 0.1. The pure anchor compound monolayers display a fully hydrophobic surface. The full atomistic simulations further increase the accuracy of our simulation procedure and reveals small additional differences and details. In particular, because of the full atomistic refinement, the lattice parameter decreased from 5.35 to 4.97 Å, the tilt angles increased from ∼24° to ∼26°, and the thickness systematically increased by 2−5 Å depending on the anchor mole fraction, whereas the hydrophobicity of the monolayer surface is not affected. Chemical Bond Order Parameters. The chemical bond order parameters are used to describe the orientation of chemical bonds relative to the surface normal. In Figure 9, the average order parameters of all backbone C−X bonds are shown for two xC16 = 0.05 monolayers, one in a random and one in a clustered spatial distribution. The order parameter can take values between −0.5, corresponding to a bond orientation perfectly perpendicular to the surface normal, and 1, corresponding to a full parallel orientation. We find that the matrix alkyl chains are highly ordered and in zigzag conformations in both monolayers. The odd-numbered bonds have an order parameter of 0.1 and the even ones of 0.6, corresponding to the angles of 50° and 25° from the surface normal. The sulfur−carbon bond close to the surface displays a tilt angle slightly larger than the carbon−carbon

Figure 9. Backbone order parameter profile of a monolayer with 5% C16 anchors in random (dashed line) and clustered distribution (solid line). The background colors indicate the molecular regions. From left to right: gold−−surface; yellow−−thiol; light gray I−−matrix alkyl chains; blue−−amide linkage; red−−OEG; light gray−−anchor alkyl chains.

bonds. The order parameters derived from our simulations are in full agreement with the earlier MD studies of SAMs.21 In the OEG region, the random and clustered assemblies differ in their bond ordering. The order parameter decreases from 0.25 to −0.1 in the randomly distributed monolayer and to 0.1 in the clustered one. In the randomly distributed monolayer, there is a long-range ordering effect, and the preferred zigzag conformation vanishes only beyond the 25th chemical bond, that is, the second ethylene glycol (EG) unit, in favor of a more disordered structure. On average, in the OEG region, the chemical bonds are tilted from the surface normal by around 60°. Similar ordering parameters can be reported for the terminal anchor alkyl chains, which are less ordered than the matrix alkyl chain. In the clustered monolayer, the order increases again after the 30th bond. The anchor alkyl chains display an ordering similar to that of the matrix alkyl chains. In the next section, we present an in-detail analysis of order parameter distributions for the three modules (matrix, OEG, and anchor) in differently composed monolayers. Additionally, torsional angles are used to distinguish between the trans versus cis conformations and helical orientations of the anchoring molecules. Comparison to IR Reflection Absorption Spectra. When it comes to IRRAS measurements, our data from atomistic simulations can rationalize observations, confirm and dismiss certain hypotheses, and help designing further experiments. On the basis of MD simulations and quantum chemistry calculations, for example,18 prominent spectral features can be assigned to certain orientational and conformational attributes of the molecules in the monolayers. We here refer the discussion of results from our simulations to experimental findings from IRRAS by Lee et al. To have a numerical means to quantify the chemical bond orientations, bond order parameter distributions were computed for the different molecular regions: matrix, OEG, and anchors. The molecular conformation is expressed as the ratio of specific torsion angles: for the matrix and anchoring alkyl chains, the fraction of trans torsional angles is considered, whereas in the case of OEG, we compute the helicity (see the Methods section). The ordering parameters and molecular conformations were calculated for the three individual single-component monolayers (matrix, C8, and C16), as well as for the mixtures at 0.05, 0.1, 0.2, and 0.5 anchor mole fractions of C8 and C16 in the matrix. For each degree of coverage, the two anchor chain lengths (C8 and C16) and the two spatial distributions (random and clustered) were investigated individually. The calculated molecular orientations, expressed as bond ordering parameter 7705


Article


Figure 10. Violin plot of CC bond order parameter distributions in the different modules of the SAMs depending on the anchor molecule chain length (C8: orange, C16: purple), spatial distribution (random: darker colored left-hand side; clustered: bright colored right-hand side) at various molar fractions.

distributions, are given in Figure 10 and the molecular conformations in Figure 11. The MD snapshots of selected monolayer compositions are shown in Figure 12 to visually represent our findings.

Figure 11. Molecular conformations of the different molecular regions in SAMs (top: matrix region, middle: OEG, bottom: anchor molecules) as a function of anchor mole fraction xA, chain length (C8 orange, C16 purple), and spatial distribution (left column: random distribution, right column: clustered distribution). Figure 12. Representative snapshots of C8 and C16 monolayers at xA = 0.1 after 20 ns full atomistic simulations. (A) C8 random, (B) C8 clustered, (C) C16 random, and (D) C16 clustered. Only a 1 nm slab is shown in this representation. Molecules are displayed as licorice. Hydrogen atoms are omitted for clarity. The matrix compound is displayed with gray and the anchor compound with green carbon atoms. As the direction of view is in the y-axis, the tilt direction appears different from snapshot to snapshot.

In the matrix region, two populations of bond order orientations can be seen independent of the chain length of the anchor compound, its mole fraction, and spatial distribution (Figure 10, top row). The two populations can be assigned to the odd-numbered chemical bonds with an average order parameter of −0.1 and the even numbered bonds with an average order parameter of 0.8. The peaks are of Gaussian shape, similar in magnitude, well-separated, and exhibit standard deviations of 0.1. This indicates highly ordered anchor molecule chains, which are mutually tilted by 30° from the surface normal. Additionally, the amount of gauche defects in the matrix region is less than 2% in all

simulated monolayers, corresponding to a trans fraction rtrans of 0.98, as shown in Figure 11 (top row). To summarize, the following results can be obtained: 7706


Article

The Journal of Physical Chemistry B • the matrix region is in a highly ordered, densely packed, all-trans conformation and • the structure of the matrix region is not significantly affected by the OEG or anchor phases. These are in excellent agreement with the reported IR spectra of Lee et al. Especially, the vibrational bands at 2918 and 2850 cm−1 are prominent for all mixtures of monolayers and originate from the alkanethiol chains in the matrix region of the monolayer. These bands are characteristic for the densely packed, all-trans alkyl chains. The order parameter distributions of the C−C bonds in the OEG region (Figure 10, central row) clearly show the phase transition of the anchor chains from a close-to-surface to an upright conformation, as already described in the previous sections. For low anchor mole fractions, there is a large population of bonds with an order parameter close to −0.5, that is, oriented perpendicular to the surface normal (parallel to the surface). The magnitude of this population is not affected by the anchor chain length and is present for all mixed monolayers in a random spatial distribution. The order parameter distributions of the pure anchor compound monolayers, however, exhibit a broad peak at 0.5, corresponding to a 30° bond tilt relative to the surface normal. Apparently, the ordering in the OEG chains occurs only at mole fractions greater than 0.5. In the clustered distribution monolayers, the order parameters look completely different. The population of order parameters close to −0.5 is strongly reduced. Yet, a single broad peak appears, which becomes even more prominent with increasing anchor mole fractions. The order parameter distribution finally converges to the same distribution as the pure anchor compound. Interestingly, there is a consistent shift of the peak toward higher ordering parameters with increasing anchor mole fractions, which shows a slight spatial extension of the OEG chains. Overall, the OEG portion displays a higher degree in disorder than the anchor alkyl chains, independent of the anchor molecule chain length and its spatial distribution, which can be seen in the broad distributions. This corresponds well with the IRRAS results, which show very similar IR spectra for the mixed monolayers in the fingerprint region of wavenumbers between 800 and 1800 cm−1, which are indicative of similar carbon bond orientations in the mixed monolayers. It has been shown experimentally60,61 and by quantum chemical calculations62 that the PEG chains partially adopt helical conformations, which can also be seen in the spectra obtained by Lee et al., and in our simulations (Figure 11, central row). We note, however, that the level of helicity in the simulation is generally lower as compared to that in the IRRAS experiment. This might arise from the fact that the simulations are performed in the liquid phase, whereas the experiments are done in the gas phase. In the simulation, the water molecules disrupt the intramolecular, stabilizing OEG hydrogen bonds and lead to a more disordered OEG phase.63 In our simulations in solution, a slight decrease in helicity with increasing anchor mole fraction for both anchor molecules and both spatial distributions from 0.5 to 0.3 can be reported. The difference between the alkyl anchor chain lengths is negligibly small. As for the different spatial distributions, we note that the largest decrease in OEG chain helicity can be detected between xA of 0.5 and 1 in randomly distributed anchoring molecules, whereas in the clustered distributions, it occurs between 0.2 and 0.5. This decrease in helicity, expressed as the amount of

O−C−C−O gauche torsion angles or their molecular areas, agrees well with other, earlier MD studies.64,65 In IRRAS experiments, bands at 2890, 1464, 1346, 1244, 115, and 963 cm−1 are found, which are interpreted as OEG chains adopting a helical conformation along the surface normal. They are not present at lower anchor mole fractions, which can be explained by either a disordering or only a low degree of surface coverage. We report a larger number of helical elements at lower anchor mole fraction monolayers. As these helices are oriented perpendicular to the surface normal, they may thus not be detectable by IRRAS or appear less prominent in the spectra because of the surface dipole selection rule. From our results, we suggest that the decrease in band intensity is rather due to an orientational transition than to an intramolecular disordering. For the anchor alkyl chains (Figure 10, bottom row), we observe a comparable behavior as for the OEG region. The pure C8 and C16 anchor monolayers form highly ordered alkyl chains in zigzag conformations, with the bond order parameter distributions similar to the matrix alkyl chains. Here, the distribution of the C16 alkyl chains is significantly more distinct than that of the C8 chains, indicating a slight disordering even in the pure C8 anchor monolayers. These observations are in full agreement with the current reflection absorption spectra obtained by Lee et al. For low anchor coverages of both the anchor chain lengths at xA < 0.5 in random spatial distribution, the C−C bonds are mainly disordered with a strong tendency for an orientation perpendicular to the surface normal. For the random distribution at xC8 = 0.5, two small peaks can be detected, indicating the onset of the orientational transition from parallel to perpendicular to the surface. For the randomly scattered C16 anchors, the distributions look similar, but these small peaks at xC16 = 0.5 are absent. In the clustered C8 monolayers, the order parameter distributions of the terminal alkyl chains are very similar in shape as that for the pure C8 anchor monolayer. Yet, the two peaks are less prominent, and a significant number of bonds are still oriented parallel to the surface (with the order parameters around −0.5). For xC8 = 0.05, only a single, broad, and centrally positioned peak is visible. For the clustered distribution C8 monolayers, the order parameter distributions indicate mutually tilted anchor alkyl chains with increasing order parameters for all anchor compositions. In the mixed monolayers, with clustered C16 anchor compound molecules, we find considerably less bonds with the order parameters around −0.5. Already at xC16 = 0.05, two distributions can be found. Unexpectedly, for xC16 = 0.1 and xC16 = 0.2, only a single peak at SCC = 0.5 is found. This, in fact, shows an even higher ordered conformation of the alkyl chains, in which the mutual tilt angle of the anchors is close to 0°, that is, an orientation which is perfectly perpendicular to the surface. This orientation is most prominent in the C16 monolayer at xC16 = 0.1, corresponding to a domain size of 40 molecules. With an increasing domain size, the C16 anchor alkyl chains begin to simultaneously tilt away from the surface normal. All of the above is supported by the ratio of the trans C−C−C−C torsion angles, as presented in Figure 11 (bottom row). In the randomly scattered anchor monolayers, the trans ratio increases linearly from 0.7, as also reported for bulk alkanes, to 0.95, similar to the matrix alkyl chains. In the clustered monolayers, the trans ratio is as high as 0.83 (C8) and 0.93 (C18) already for xA of 0.05. The trans ratio in the 7707


Article


necessary to form these highly organized anchor phases might even be smaller compared to the simulations.

clustered C8 anchor chains exceeds 0.9 for the monolayers of xC8 ≥ 0.5. In IRRAS experiments by Lee et al., asymmetric and symmetric CD2 stretching modes at 2193 and 2090 cm−1, and the corresponding asymmetric and symmetric CD3 stretching modes at 2216 and 2074 cm−1, are reported for the C16 anchor monolayers. These bands are interpreted as the d-alkyl chains to adopt a highly ordered, all-trans conformation. Even though these bands become less prominent and red-shifted with a decreasing anchor mole fraction, the overall shape of the spectrum is conserved. At the same time, the spectrum of the C8 anchor monolayers undergoes a change in spectral shape with decreasing anchor mole fractions. The most prominent bands are the CD2 asymmetric and symmetric stretching modes at 2189 and 2101 cm−1, as well as the corresponding asymmetric and symmetric CD3 stretching modes at 2221 and 2074 cm−1. With regard of our MD simulations, it appears not feasible for the C8 anchor molecules to form domains on the surface. The orientation and conformation of the randomly distributed C8 anchor molecules are in qualitative agreement with the experiment. Thus, the simulations clearly favor a random spatial distribution of the C8 anchoring compound. On the other hand, the C16 anchor molecules in a clustered spatial distribution deploy nearly all-trans conformation already at low mole fractions. The number of bonds that are parallel to the surface is significantly lower than that in the corresponding mixed C8 monolayers. For the C16 anchoring compound, the calculated features of only the clustered distribution monolayers can explain the IR reflection absorption spectra. Thus, our simulations strongly support the hypothesis of Lee et al. that phase segregation occurs exclusively for the longer C16 anchoring compound. As an explanation, we suggest that the C16 anchor compound initially adsorbs in a random spatial distribution and then self-associates. In the low anchor mole fraction regime, with xC8 smaller than 0.3, and xC16 smaller than 0.2, it is likely that the C8 alkyl chains do not interact, whereas the C16 alkyl chains already do (see Figure 6). This interaction at the surface and the aggregation of alkyl chains could be the driving forces for a slow phase separation on the surface. Place exchange diffusion has recently been observed on gold nanoparticles66,67 and might be enhanced by the anchor chain interactions. If we look at the molecular conformations of C8 and C16 at the same mole fraction of xA = 0.1 with different chain lengths and spatial distributions (see Figure 12), the differences become most apparent. For C8, the anchor alkyl chains aggregate only to a small extent at the surface. In the C16 monolayer, however, there are more and stronger interactions. The clustered C16 monolayer is already wellordered at the low anchor mole fraction. In the clustered monolayers, we observe a significant degree of disorder, especially in the bordering regions. Our simulations are conducted in an aqueous environment, whereas the IRRAS studies are performed with a dry sample. In air, the intra- and intermolecular interactions (hydrogen bonding, hydrophobic effects) in the OEG and anchor alkyl chains are more pronounced, resulting from a lack of competing interactions with the solvent molecules. Thus, in direct comparison to the experiments by Lee et al., our simulations slightly underestimate the alkyl chain ordering and OEG helix formation. We conclude that in a dry sample, the domain size which is

■

CONCLUSIONS SAMs provide an ideally suitable and well-controllable microenvironment for the incorporation and immobilization

Figure 13. Future application for the developed workflow: tBLMs (left) and tethered lipid vesicles (right). The CG simulation will be valuable to equilibrate bilayers and vesicles for long timescales before the model transformation to further simulations with an improved allatom accuracy.

of large biomolecules. Matrix molecules and long-chain anchoring components consisting of OEG (here 6 units) and C8 and C16 long alkyl chains provide the possibility for tethering phospholipid bilayers and vesicles. The design and setup of a consistent simulation protocol for full atomistic simulations of mixed self-assembled monolayers on flat model surfaces, based on previous CG initialization, modeling, and equilibration, were developed. The results of both the CG and the full atomistic models agree well with the experimentally reported structural parameters of the self-assembled monolayers. We demonstrate the feasibility of a CG model to bring any mixed SAM system including the membrane-anchoring compounds into thermodynamic equilibrium to start subsequent full atomistic simulations. By simulating different types of C8 and C16 long anchor molecules over the complete range of molar fractions, significant differences in packing and interactions can be obtained. The structural parameters such as lattice constant and tilt angle are almost independent of the type of the anchoring molecule. The anchor alkyl-chain orientation, however, is critically controlled by the composition and its chemical nature. Although the short chain C8 anchor molecules adopt a random, disordered conformation at low anchor densities and undergo a transition into a more ordered conformation with increasing anchor density, the C16 anchors already adopt highly ordered conformations at low anchor densities, which reoccur in the monolayers of higher anchor density. The results from our MD simulations are in excellent agreement with current experiments by Lee et al. and explain the observed spectral features. In addition, the simulations are able to probe mixed concentrations and mole fractions that are experimentally hard to access. Complementing the experimental IRRAS measurements, our simulations results provide additional insights into the arrangement and degree of ordering of the anchoring compounds in the monolayer. Depending on the initial spatial anchor compound distribution, we find distinct orientational and conformational phase transitions. Our simulations support the hypothesis that the C16 anchor compound exclusively forms self-aggregates, whereas the C8 anchor compound scatters randomly in the matrix. The accuracy of the results might even be optimized by using a more elaborate force field or surface model. Anyway, we consider the model accurate enough to be employed as a 7708


Article


Fluorescence-Interference Detection. Biophys. J. 2005, 88, 4289− 4302. (12) Lata, S.; Gavutis, M.; Piehler, J. Monitoring the Dynamics of Ligand−Receptor Complexes on Model Membranes. J. Am. Chem. Soc. 2006, 128, 6−7. (13) Zharnikov, M.; Grunze, M. Spectroscopic Characterization of Thiol-Derived Self-Assembling Monolayers. J. Phys.: Condens. Matter 2001, 13, 11333−11365. (14) Gavutis, M.; Lata, S.; Piehler, J. Probing 2-Dimensional Protein-Protein Interactions on Model Membranes. Nat. Protoc. 2006, 1, 2091−2103. (15) Benesch, J.; Svedhem, S.; Svensson, S. C. T.; Valiokas, R.; Liedberg, B.; Tengvall, P. Protein Adsorption to Oligo(Ethylene Glycol) Self-Assembled Monolayers: Experiments with Fibrinogen, Heparinized Plasma, and Serum. J. Biomater. Sci., Polym. Ed. 2001, 12, 581−597. (16) Pale-Grosdemange, C.; Simon, E. S.; Prime, K. L.; Whitesides, G. M. Formation of Self-Assembled Monolayers by Chemisorption of Derivatives of Oligo(Ethylene Glycol) of Structure Hs(Ch2)11(Och2ch2)Moh on Gold. J. Am. Chem. Soc. 1991, 113, 12−20. (17) Harder, P.; Grunze, M.; Dahint, R.; Whitesides, G. M.; Laibinis, P. E. Molecular Conformation in Oligo(Ethylene Glycol)-Terminated Self-Assembled Monolayers on Gold and Silver Surfaces Determines Their Ability to Resist Protein Adsorption. J. Phys. Chem. B 1998, 102, 426−436. (18) Lee, H.-H.; Ruželė, Z.; Malysheva, L.; Onipko, A.; Gutés, A.; Björefors, F.; Valiokas, R.; Liedberg, B. Long-Chain Alkylthiol Assemblies Containing Buried In-Plane Stabilizing Architectures†. Langmuir 2009, 25, 13959−13971. (19) Valiokas, R.; Malysheva, L.; Onipko, A.; Lee, H.-H.; Ruželė, Ž .; Svedhem, S.; Svensson, S. C. T.; Gelius, U.; Liedberg, B. On the Quality and Structural Characteristics of Oligo(Ethylene Glycol) Assemblies on Gold: An Experimental and Theoretical Study. J. Electron Spectrosc. Relat. Phenom. 2009, 172, 9−20. (20) Laibinis, P. E.; Fox, M. A.; Folkers, J. P.; Whitesides, G. M. Comparisons of Self-Assembled Monolayers on Silver and Gold: Mixed Monolayers Derived from Hs(Ch2)21x and Hs(Ch2)10y (X, Y = Ch3, Ch2oh) Have Similar Properties. Langmuir 1991, 7, 3167− 3173. (21) Hautman, J.; Klein, M. L. Simulation of a Monolayer of Alkyl Thiol Chains. J. Chem. Phys. 1989, 91, 4994−5001. (22) Hautman, J.; Bareman, J. P.; Mar, W.; Klein, M. L. Molecular Dynamics Investigations of Self-Assembled Monolayers. Faraday J. Chem. Soc., Faraday Trans. 1991, 87, 2031. (23) Vemparala, S.; Karki, B. B.; Kalia, R. K.; Nakano, A.; Vashishta, P. Large-Scale Molecular Dynamics Simulations of Alkanethiol SelfAssembled Monolayers. J. Chem. Phys. 2004, 121, 4323−4330. (24) Esteban-Martín, S.; Salgado, J. Self-Assembling of Peptide/ Membrane Complexes by Atomistic Molecular Dynamics Simulations. Biophys. J. 2007, 92, 903−912. (25) Mori, T.; Miyashita, N.; Im, W.; Feig, M.; Sugita, Y. Molecular Dynamics Simulations of Biological Membranes and Membrane Proteins Using Enhanced Conformational Sampling Algorithms. Biochim. Biophys. Acta 2016, 1858, 1635−1651. (26) Bhadra, P.; Siu, S. W. I. Comparison of Biomolecular Force Fields for Alkanethiol Self-Assembled Monolayer Simulations. J. Phys. Chem. C 2017, 121, 26340−26349. (27) Marrink, S. J.; Tieleman, D. P. Perspective on the Martini Model. Chem. Soc. Rev. 2013, 42, 6801−6822. (28) Marrink, S. J.; de Vries, A. H.; Mark, A. E. Coarse Grained Model for Semiquantitative Lipid Simulations. J. Phys. Chem. B 2004, 108, 750−760. (29) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The Martini Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (30) Lee, H.; de Vries, A. H.; Marrink, S.-J.; Pastor, R. W. A CoarseGrained Model for Polyethylene Oxide and Polyethylene Glycol: Conformation and Hydrodynamics. J. Phys. Chem. B 2009, 113, 13186−13194.

basis for future applications in the simulation of tethered bilayer membranes and vesicles, as presented in Figure 13. Finally, we consider the developed simulation protocol and atomistic model generation for mixed composition SAMs accurate, convenient, and extendable to be used in further simulations of tBLMs and vesicles.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.8b05075.

■

United atom simulation protocol, detailed description of the simulated systems, CG force field and input parameters, and all-atom refinement and simulation protocol (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +493916110436. ORCID

Matthias Stein: 0000-0001-7793-0052 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research is supported by the Max Planck Society for the Advancement of Science and the International Max Planck Research School (IMPRS) for Advanced Methods in Process and Systems Engineering, Magdeburg, Germany.

■

REFERENCES

(1) Bigelow, W. C.; Pickett, D. L.; Zisman, W. A. Oleophobic Monolayers. J. Colloid Sci. 1946, 1, 513−538. (2) Ulman, A. Formation and Structure of Self-Assembled Monolayers. Chem. Rev. 1996, 96, 1533−1554. (3) Schreiber, F. Structure and Growth of Self-Assembling Monolayers. Prog. Surf. Sci. 2000, 65, 151−257. (4) Knoll, W.; Frank, C. W.; Heibel, C.; Naumann, R.; Offenhäusser, A.; Rühe, J.; Schmidt, E. K.; Shen, W. W.; Sinner, A. Functional Tethered Lipid Bilayers. Rev. Mol. Biotechnol. 2000, 74, 137−158. (5) Naumann, R.; Schiller, S. M.; Giess, F.; Grohe, B.; Hartman, K. B.; Kärcher, I.; Köper, I.; Lübben, J.; Vasilev, K.; Knoll, W. Tethered Lipid Bilayers on Ultraflat Gold Surfaces. Langmuir 2003, 19, 5435− 5443. (6) Giess, F.; Friedrich, M. G.; Heberle, J.; Naumann, R. L.; Knoll, W. The Protein-Tethered Lipid Bilayer: A Novel Mimic of the Biological Membrane. Biophys. J. 2004, 87, 3213−3220. (7) Castellana, E. T.; Cremer, P. S. Solid Supported Lipid Bilayers: From Biophysical Studies to Sensor Design. Surf. Sci. Rep. 2006, 61, 429−444. (8) Andersson, J.; Köper, I. Tethered and Polymer Supported Bilayer Lipid Membranes: Structure and Function. Membranes 2016, 6, 30. (9) Oelmeier, S. A.; Dismer, F.; Hubbuch, J. Molecular Dynamics Simulations on Aqueous Two-Phase Systems - Single Peg-Molecules in Solution. BMC Biophys. 2012, 5, 14. (10) Lamken, P.; Gavutis, M.; Peters, I.; Van der Heyden, J.; Uzé, G.; Piehler, J. Functional Cartography of the Ectodomain of the Type I Interferon Receptor Subunit Ifnar1. J. Mol. Biol. 2005, 350, 476− 488. (11) Gavutis, M.; Lata, S.; Lamken, P.; Müller, P.; Piehler, J. Lateral Ligand-Receptor Interactions on Membranes Probed by Simultaneous 7709


Article

The Journal of Physical Chemistry B (31) Lee, H.; Pastor, R. W. Coarse-Grained Model for Pegylated Lipids: Effect of Pegylation on the Size and Shape of Self-Assembled Structures. J. Phys. Chem. B 2011, 115, 7830−7837. (32) Kyrychenko, A.; Karpushina, G. V.; Bogatyrenko, S. I.; Kryshtal, A. P.; Doroshenko, A. O. Preparation, Structure, and a CoarseGrained Molecular Dynamics Model for Dodecanethiol-Stabilized Gold Nanoparticles. Comput. Theor. Chem. 2011, 977, 34−39. (33) Monticelli, L.; Kandasamy, S. K.; Periole, X.; Larson, R. G.; Tieleman, D. P.; Marrink, S.-J. The Martini Coarse-Grained Force Field: Extension to Proteins. J. Chem. Theory Comput. 2008, 4, 819− 834. (34) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. 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. 2004, 25, 1656−1676. (35) Wassenaar, T. A.; Ingólfsson, H. I.; Böckmann, R. A.; Tieleman, D. P.; Marrink, S. J. Computational Lipidomics with Insane: A Versatile Tool for Generating Custom Membranes for Molecular Simulations. J. Chem. Theory Comput. 2015, 11, 2144−2155. (36) Ingólfsson, H. I.; Melo, M. N.; van Eerden, F. J.; Arnarez, C.; Lopez, C. A.; Wassenaar, T. A.; Periole, X.; de Vries, A. H.; Tieleman, D. P.; Marrink, S. J. Lipid Organization of the Plasma Membrane. J. Am. Chem. Soc. 2014, 136, 14554−14559. (37) Edler, E.; Schulze, E.; Stein, M. Membrane Localization and Dynamics of Geranylgeranylated Rab5 Hypervariable Region. Biochim. Biophys. Acta 2017, 1859, 1335−1349. (38) Vericat, C.; Vela, M. E.; Salvarezza, R. C. Self-Assembled Monolayers of Alkanethiols on Au(111): Surface Structures, Defects and Dynamics. Phys. Chem. Chem. Phys. 2005, 7, 3258−3268. (39) Yesylevskyy, S. O.; Schäfer, L. V.; Sengupta, D.; Marrink, S. J. Polarizable Water Model for the Coarse-Grained Martini Force Field. PLoS Comput. Biol. 2010, 6, No. e1000810. (40) de Jong, D. H.; Baoukina, S.; Ingólfsson, H. I.; Marrink, S. J. Martini Straight: Boosting Performance Using a Shorter Cutoff and Gpus. Comput. Phys. Commun. 2016, 199, 1−7. (41) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (42) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (43) Wassenaar, T. A.; Pluhackova, K.; Böckmann, R. A.; Marrink, S. J.; Tieleman, D. P. Going Backward: A Flexible Geometric Approach to Reverse Transformation from Coarse Grained to Atomistic Models. J. Chem. Theory Comput. 2014, 10, 676−690. (44) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of the Charmm All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. J. Phys. Chem. B 2010, 114, 7830−7843. (45) Huang, J.; MacKerell, A. D., Jr. Charmm36 All-Atom Additive Protein Force Field: Validation Based on Comparison to Nmr Data. J. Comput. Chem. 2013, 34, 2135−2145. (46) Lee, J.; Cheng, X.; Swails, J. M.; Yeom, M. S.; Eastman, P. K.; Lemkul, J. A.; Wei, S.; Buckner, J.; Jeong, J. C.; Qi, Y.; et al. CharmmGui Input Generator for Namd, Gromacs, Amber, Openmm, and Charmm/Openmm Simulations Using the Charmm36 Additive Force Field. J. Chem. Theory Comput. 2016, 12, 405−413. (47) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. Charmm General Force Field: A Force Field for Drug-Like Molecules Compatible with the Charmm All-Atom Additive Biological Force Fields. J. Comput. Chem. 2009, 31, 671. (48) Azzam, R.; Bashara, N.; Burns, D. T. Ellipsometry and Polarized Light; Elsevier: North Holland, Amsterdam, 1987; Vol. Xvii+, p 539. (49) van der Ploeg, P.; Berendsen, H. J. C. Molecular Dynamics Simulation of a Bilayer Membrane. J. Chem. Phys. 1982, 76, 3271− 3276.

(50) Oelmeier, S. A.; Dismer, F.; Hubbuch, J. Molecular Dynamics Simulations on Aqueous Two-Phase Systems - Single Peg-Molecules in Solution. BMC Biophys. 2012, 5, 14. (51) Humphrey, W.; Dalke, A.; Schulten, K. Vmd: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (52) Ulman, A.; Eilers, J. E.; Tillman, N. Packing and Molecular Orientation of Alkanethiol Monolayers on Gold Surfaces. Langmuir 1989, 5, 1147−1152. (53) Vericat, C.; Vela, M. E.; Benitez, G.; Carro, P.; Salvarezza, R. C. Self-Assembled Monolayers of Thiols and Dithiols on Gold: New Challenges for a Well-Known System. Chem. Soc. Rev. 2010, 39, 1805−1834. (54) Porter, M. D.; Bright, T. B.; Allara, D. L.; Chidsey, C. E. D. Spontaneously Organized Molecular Assemblies. 4. Structural Characterization of N-Alkyl Thiol Monolayers on Gold by Optical Ellipsometry, Infrared Spectroscopy, and Electrochemistry. J. Am. Chem. Soc. 1987, 109, 3559−3568. (55) Nuzzo, R. G.; Zegarski, B. R.; Dubois, L. H. Fundamental Studies of the Chemisorption of Organosulfur Compounds on Gold(111). Implications for Molecular Self-Assembly on Gold Surfaces. J. Am. Chem. Soc. 1987, 109, 733−740. (56) Shi, B.; Dhir, V. K. Molecular Dynamics Simulation of the Contact Angle of Liquids on Solid Surfaces. J. Chem. Phys. 2009, 130, 034705. (57) Koishi, T.; Yasuoka, K.; Fujikawa, S.; Zeng, X. C. Measurement of Contact-Angle Hysteresis for Droplets on Nanopillared Surface and in the Cassie and Wenzel States: A Molecular Dynamics Simulation Study. ACS Nano 2011, 5, 6834−6842. (58) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Noca, F.; Koumoutsakos, P. Molecular Dynamics Simulation of Contact Angles of Water Droplets in Carbon Nanotubes. Nano Lett. 2001, 1, 697− 702. (59) Hoiles, W.; Gupta, R.; Cornell, B.; Cranfield, C.; Krishnamurthy, V. The Effect of Tethers on Artificial Cell Membranes: A Coarse-Grained Molecular Dynamics Study. PLoS One 2016, 11, No. e0162790. (60) Alessi, M. L.; Norman, A. I.; Knowlton, S. E.; Ho, D. L.; Greer, S. C. Helical and Coil Conformations of Poly(Ethylene Glycol) in Isobutyric Acid and Water. Macromolecules 2005, 38, 9333−9340. (61) Azri, A.; Giamarchi, P.; Grohens, Y.; Olier, R.; Privat, M. Polyethylene Glycol Aggregates in Water Formed through Hydrophobic Helical Structures. J. Colloid Interface Sci. 2012, 379, 14−19. (62) Wang, R. L. C.; Kreuzer, H. J.; Grunze, M. The Interaction of Oligo(Ethylene Oxide) with Water: A Quantum Mechanical Study. Phys. Chem. Chem. Phys. 2000, 2, 3613−3622. (63) Zolk, M.; Eisert, F.; Pipper, J.; Herrwerth, S.; Eck, W.; Buck, M.; Grunze, M. Solvation of Oligo(Ethylene Glycol)-Terminated SelfAssembled Monolayers Studied by Vibrational Sum Frequency Spectroscopy. Langmuir 2000, 16, 5849−5852. (64) Zheng, J.; Li, L.; Chen, S.; Jiang, S. Molecular Simulation Study of Water Interactions with Oligo (Ethylene Glycol)-Terminated Alkanethiol Self-Assembled Monolayers. Langmuir 2004, 20, 8931− 8938. (65) Ismail, A. E.; Grest, G. S.; Stevens, M. J. Structure and Dynamics of Water near the Interface with Oligo(Ethylene Oxide) Self-Assembled Monolayers. Langmuir 2007, 23, 8508−8514. (66) Kassam, A.; Bremner, G.; Clark, B.; Ulibarri, G.; Lennox, R. B. Place Exchange Reactions of Alkyl Thiols on Gold Nanoparticles. J. Am. Chem. Soc. 2006, 128, 3476−3477. (67) Montalti, M.; Prodi, L.; Zaccheroni, N.; Baxter, R.; Teobaldi, G.; Zerbetto, F. Kinetics of Place-Exchange Reactions of Thiols on Gold Nanoparticles. Langmuir 2003, 19, 5172−5174.

7710


Simulation of Mixed Self-Assembled Monolayers on Gold: Effect of

Recommend Documents