Peptide Recognition Capabilities of Cellulose in Molecular Dynamics

Sep 30, 2015 - Peptide Recognition Capabilities of Cellulose in Molecular Dynamics ... We notice that the rough and dynamic surface of cellulose favor...
4 downloads 0 Views 8MB Size
Article pubs.acs.org/JPCC

Peptide Recognition Capabilities of Cellulose in Molecular Dynamics Simulations Grzegorz Nawrocki,†,‡,∥ Pierre-André Cazade,‡,§ Damien Thompson,‡,§ and Marek Cieplak*,∥ †

Department of Biochemistry and Molecular Biology, Michigan State University, East Lansing, Michigan 48824, United States Department of Physics and Energy and §Materials and Surface Science Institute, University of Limerick, Limerick, Ireland ∥ Institute of Physics, Polish Academy of Sciences, Aleja Lotników 32/46, 02-668 Warsaw, Poland ‡

S Supporting Information *

ABSTRACT: Knowing the amino acid discriminating capacity of the surface of cellulose can be valuable information toward the rational design and re-engineering of the cellulosome in order to improve its catalytic properties. This aim can be achieved by the determination of the binding free energies of amino acids in molecular dynamics simulations. Conventional simulations do not always allow for sufficient sampling of the configuration space, especially in a system where large energy barriers occur. A better sampling can be obtained by replica exchange molecular dynamics (REMD). Here, we use REMD combined with umbrella sampling to determine the potential of the mean force for amino acids, analogues of their side chains, dipeptides, and tripeptides for the interactions with the (100) face of the crystalline cellulose Iβ. We also use REMD to characterize the adsorption dynamics of a small protein, tryptophan cage. Our results show that all 20 standard amino acids adsorb on the surface of cellulose with binding energies ranging from 3 to 10 kJ/mol with the specificity of approximately 30%. The largest affinity to cellulose is shown to belong to the aromatic residues. On the basis of simulations of selected dipeptides and tripeptides, we determine that cellulose-adsorption energies should be scaled down to perhaps about 20% to correctly describe adsorption of protein residues as opposed to free amino acid molecules. We notice that the rough and dynamic surface of cellulose favors adsorption of flexible parts of the protein, in contrast to the specificity expressed by flat and static surfaces of nonorganic solids. We show that the contact angle for a droplet of water on the (100) face of the cellulose is about 24°, reflecting the heterogeneous hydrophilicity of the surface.



choice of the force fieldwhether it yields hydrophobic or hydrophilic behavior and whether it incorporates the electric polarization of the metal. Another aspect of interest in the comparative studies of adsorption is the behavior of water near the surfaces studied. It can be characterized through the density and polarization profiles. These profiles have also been found to depend on the surface.14−16 In particular, the first water layer that forms near the strongly hydrophilic ZnO14 corresponds to packing at such a density that single amino acids are unable to penetrate it and reach the solid. In this paper, we extend our comparative studies to organic surfaces and consider the one formed by the crystalline cellulose Iβ19−21 that is stabilized by hydrogen bonds22,23 and exposes its (100) face. Cellulose,24 in its various forms, is the main component of plant cell walls and thus also of biomass.25 The enzymes known as cellulases and their complexes with scaffoldin, known as cellulosomes,26 can break cellulose down into simple sugars. Enhancing the efficiency of this process requires studies on how they do it and on the role of various supporting proteins. Such studies should then suggest

INTRODUCTION Understanding the mechanisms of protein adsorption on surfaces is important for many areas of fundamental and applied research.1−4 Examples include biomineralization,5−8 design of biosensors,9−11 behavior of biofunctionalized solid nanoparticles,12 and catalysis.13 A systematic comparison of the adsorptive properties of various surfaces can be accomplished by making assessments involving a fixed set of the reference biomolecules instead of considering specific applications at hand. We have initiated such a program by selecting the 20 natural amino acids and a small protein, tryptophan cage, as the reference systems, together with some selected dipeptides. Another reason to include amino acids is that proteins interact with surfaces primarily through their exposed amino acids. If the interactions with the residues are sufficiently specific, they may steer the adsorption process even if their strengths are modified by the presence of nearby residues. Our assessments were done through all-atom molecular dynamics simulations. The surfaces compared so far were four cleavage faces of ZnO,14 the (110) face of ZnS,15 the (111) surface of gold,16 andonly for selected proteinsmica.17,18 The adsorbing powers and the specificities have been found to be sensitively dependent on the surface, including the selection of the particular face. In the case of gold, they also depend on the © 2015 American Chemical Society

Received: July 22, 2015 Revised: September 28, 2015 Published: September 30, 2015 24404

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

The space above the slab was filled by the TIP3P molecules of water.35 The overall geometry at the water−crystal interface is illustrated by the snapshots shown in Figure 1, which also define the system of the coordinates used.

mutagenetic ways to design improved protein complexes that degrade biomass. Determination of the adsorptive characteristics of cellulose for biomolecules and water should provide important guiding principles for such studies. It should be noted that ordered water layers at the cellulose−water interface may hinder diffusion of the enzymes toward the substrate,27 so it is important to understand which residues favor movements toward the surface and which do not. We demonstrate that cellulose Iβ is mildly hydrophilica droplet of water makes a contact angle of about 24° with the surface. We also show that the planar distribution of water density is affected by the heterogeneous properties of the cellulose. The primary goal of this paper is to determine the potential of the mean force (PMF)28 for single capped amino acids interacting with the cellulose. We do this through replica exchange molecular dynamics (REMD), in which the cellulose is allowed to evolve in time. In our previous studies, we used standard molecular dynamics and the surfaces were considered to be frozen. We find that the residue averaged binding energy is 6.3 kJ/mol, which is comparable to that of the stronger hydrogen bonds, and that the specificity is on the order of 28%. The highest binding affinity is found to be associated with the aromatic amino acids Tyr, Phe, and Trp (10.3, 9.4, and 8.9 kJ/ mol, respectively). This is consistent with previous theoretical29 and experimental30,31 studies of several proteins, including the CBM (cellulose binding module). The weakest affinity is found for the charged Asp (3.2 kJ/mol). We probe these values by considering a range of amino acid analogues. It is interesting to note that the distinguished role of Tyr and Phe is also seen in the experimental studies of adsorption of amino acids to βcyclodextrins.32,33 Another goal of this paper is to consider dipeptides and tripeptides and to demonstrate that the binding energies are scaled systematically compared to the sum of the single amino acid contributions to binding. We then consider the tryptophan cage and study the involvement of the individual amino acids in the process of binding.

Figure 1. Snapshots of the molecules at the cellulose−water interface as viewed from three directions. The oxygen, hydrogen, and carbon atoms are colored white, red, and cyan, respectively. The rightmost panel shows the top layer of the cellulose chains and its first layer of water molecules.

The concentration of the background ions corresponded to the physiological 150 mM, i.e., there were six Na+ ions and six Cl− ions. If a biomolecule had a net charge on its side chains, extra ions were inserted to neutralize the system. The amino acids and dipeptides were capped by an acetyl group at the Nterminus and N-methylamide group at the C-terminus. These caps eliminate the terminal charges and mimic the presence of a peptide chain in which the amino acids adopt the un-ionized forms (see, for example, ref 36). Histidine is considered in its three possible protonation states: HIE (H on the ϵ N atom), HID (H on the δ N atom), and the positively charged HIP (H on both ϵ and δ N atoms). Histidine’s pKa value is 6, which means that at the neutral pH of 7, all three protonation forms coexist. In order to determine the relevance of a given chemical group in the process of adsorption, several analogues of amino acid side chains were considered. These are methane, which is an analogue of Ala; butylammonium and ammonium, both analogues of Lys; methanol, an analogue of Ser; ethanoate and methanoate, analogues of Asp; guanidinium, an analogue of Arg; and benzene, an analogue of Phe. When considering dipeptides, six representative amino acids were used: positively charged Lys, negatively charged Asp, polar Ser, hydrophobic Leu, aromatic Phe, and the smallest amino acid, Gly. We combine them into dipeptides in all possible ways. We also consider tripeptides, but only with identical repeats, those made exclusively of Asp, Ser, Leu, Phe, and Gly. As the reference protein, we take the tryptophan cage, with the PDB code of 1L2Y. This protein consists of 20 amino acids and has a net total charge of +1e. Simulation Settings. The simulations were performed using the GROMACS 4.6.5 package37 with the CHARMM36 force field38 that contained parameters for proteins,39−41 monosaccharides,42 and the bonds between them.43 The time integration of the equations of motion was performed using the leapfrog algorithm with a time step of 1 fs. The energy of the system was minimized using the steepest descent algorithm. The cutoff radius for both the electrostatic and van der Waals



METHODS Geometry of the System and the Biomolecules Studied. The crystalline structure of cellulose was generated in its Iβ form by using the Cellulose-Builder toolkit.34 The crystal was modeled as a slab of unit cells, each containing three layers of five nearly parallel chains that were eight-monomers long. The terminal molecules of glucose on each chain were bonded together across the periodic walls of the box. The layers extended primarily in the x−y-plane, and the x-direction runs along the cellulose chain from the C1 atom to the C4 atom. The direction normal to the surface is denoted as z. The exposed surface is (100). The slab was placed at the bottom of the simulation box, leaving a space of approximately 4 × 4 × 4 nm above it. For the simulations involving a water droplet (see below) an enlarged cell size of 20 × 20 × 20 nm was used. Unlike the systems with the solid surfaces considered by us previously,14−16 we do not freeze the atoms of the crystal. During the simulations involving umbrella sampling, the deepest layer of the chains was restrained harmonically to keep it in its bulk form, and the spring constant used was 1000 kJ/mol/nm2. In simulations of the protein, the periodic boundary conditions allow the protein to adsorb from both sides of the slab and then we restrain the middle layer. 24405

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

The weighted histogram analysis method (WHAM) algorithm47 was used to derive the PMF profile through umbrella sampling combined with the REMD approach. The statistical errors were estimated by bootstrap analysis.50 The first 1 ns was considered to be an equilibration part and was thus excluded from the analysis. In the MD approach, the force was time-averaged and then integrated over the z-coordinate. We have found that both MD and REMD produce almost identical functional forms of the PMF profiles. The depth of the lowest negative minimum in the PMF profile defines the energy parameter (ϵ) and its location indicates the optimum separation (σ). The parameter ϵ is an approximate measure of the Gibbs free energy of binding. Free energy profiles obtained by umbrella sampling and the WHAM method do not directly provide binding free energies. The latter correspond to the free energy difference between bound and unbound states,51 according to bonding constant

interactions was set to 1.2 nm. The truncated van der Waals potential was shifted so that it vanished at the cutoff distance. The particle mesh Ewald procedure44 with a grid spacing of 0.16 nm was used to treat the long-range electrostatics. The periodic boundary conditions were applied in all directions. Temperature was constrained, in most cases, to 300 K, by using the Berendsen thermostat with a stochastic term that ensures that a proper canonical ensemble is generated.45 The coupling to the isotropic pressure reservoir was implemented every 2 ps. The reference pressure was equal to 1 bar and compressibility to 4.5 × 10−5 bar−1. The LINCS (linear constraint solver) algorithm was used to constrain the H atoms. The VMD package was used for visualization of the molecules in the system.46 Energy Minimization and Equilibration of the System. We first considered a slab of cellulose in a vacuum. The energy was minimized and then the system was equilibrated for 1 ns. In the second stage, the slab was immersed in water and the energy was minimized again. This was followed by 10 ns of equilibration. Properties of the water, i.e., density and polarization, were determined from the last 1 ns of the equilibrated simulation. The final solvated surface structure was used as the starting point for the subsequent tests of cellulose binding to biomolecules. The biomolecules were placed in the box and the background salt concentration was modified whenever necessary. The energy was minimized and the system was equilibrated for 10 ps with restraints on the center of mass of the biomolecules. An additional equilibration was effectively performed “on the fly” when the biomolecules were approaching the surface, as in the umbrella sampling simulations. In simulations of the water droplet, the first 1 ns of the simulation was used as an equilibration to enable sufficient spreading on the surface. The subsequent 4 ns of molecular dynamics at room temperature were used in the analysis. Umbrella Sampling and Replica Exchange Methods. We determined the PMF by implementing the umbrella sampling method.47,48 The procedure involved two stages. In the first stage, a set of initial conformations was generated for representative values of z by pulling the center of mass of the biomolecule (an amino acid, an analogue, a dipeptide, or a tripepdtide) by a “dummy particle” moving along the z-axis toward the surface. The speed of pulling was 1 nm/ns and the starting vertical distance corresponded to 2 nm. The spring constant involved was 5000 kJ/(mol nm2). Conformations were saved every 0.1 ps from the resulting trajectory and those corresponding to the discretized values of z, separated by 0.05 nm, were selected for further studies. In the second stage, the selected conformations were used for additional 5 ns separate runs, where the z-location of the pulling particle was fixed, whereas the center of mass was observed to move within a sampling window of width Δz. The spring constant was now chosen so that Δz values from neighboring bins allowed for an overlap. We have used the umbrella sampling method in two modes: regular molecular dynamics (MD), as in our previous papers,16 and with REMD. The efficiency of umbrella sampling has been found to increase when using the sampling windows as independent replicas and by exchanging locations of the biasing potentials of the adjacent replicas every 50 ps. However, unlike the procedure used in, for example, ref 49, we do not introduce any criteria of acceptance for the exchange to take place.

ΔG bind = −kBT log(Kb) Kb =

(1)

∫V P(x ,y ,z) dx dy dz

Vbox

Vbox

∫V P(x ,y ,z) dx dy dz

b

u

(2)

Vbox is the volume of the simulation box, Vb the volume of the bound domain, and Vu the volume of the unbound domain. P(x,y,z) is the probability of finding the binding molecule at (x,y,z). Neglecting the surface roughness, P(x,y,z) is constant within a slice of thickness dz (parallel to the surface). Since the derived ΔGWHAM is averaged in the x−y-plane and z plays the role of the reaction coordinate, we get ΔG bind = −kBT log log

∫Z

⎛ ΔG WHAM(z) ⎞ exp⎜ − ⎟ dz + kBT kBT ⎝ ⎠ b

⎛ ΔG WHAM(z) ⎞ ⎟ dz kBT ⎠

∫Z exp⎜⎝− u

(3)

The bound domain, Zb, along the z axis stretches from the smallest value of z out to the first maximum, and the unbound domain, Zu, ranges from the distance at which the profile plateaus out to the largest value of z. In the present work, due to sometimes complex profiles, the best results are given by dividing the profile into two domains: from the smallest distance to 0.9 nm and from 0.9 nm to the largest distance. In our studies of adsorption of the tryptophan cage, the protein was allowed to evolve freely for 20 ns, during which time the adsorption events were observed. We then performed REMD for an additional 10 ns to enhance the efficiency of sampling of the conformational space. The system was simulated using 48 replicas simultaneously, each at a different temperature in the range between 300 and 407 K. The trajectory generated at 300 K was the only one analyzed. The whole procedure was repeated twice with different starting orientations of the protein with respect to the surface, so the results shown are averaged altogether over 2 × 10 ns.



RESULTS Water. Density Profile. The top panel of Figure 2 shows the number density profile of neat water above the surface of the cellulose. It is seen that the profile undergoes undulations similar to what is observed already in generic Lennard-Jones fluids for various strengths of interactions with atomic walls.52,53 Two density layers can be identified, and we denote 24406

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

Electric Polarization. The approximately uniform distribution of the water electric polarization in Figure 3 signifies that

Figure 3. Distributions of water polarization for the first two water layers; the darker color corresponds to layer I and the lighter to layer II, as defined in Figure 2. P is the magnitude of the polarization vector, and Pα denotes its α-component. nP is the number of water molecules with a given polarization divided by the number of all water molecules in the two layers and expressed as a percentage. The contents of the bins in the two layers sum to 100%. The typical errors for both layers are shown in selected bins.

Figure 2. Time-averaged number density, n, of water molecules above the surface of the cellulose. The averaging bins had a linear size of 0.075 nm in each Cartesian direction. This distance is approximately equal to half of the diameter of the water molecule. n ≈ 100 corresponds to bulk water at standard temperature and pressure. The top panel: the number density profile of water molecules as a function of z. This coordinate is measured from the center of mass of the topmost layer of the cellulose chains. The regions denoted by I and II correspond to the first and second layers, respectively. The bottom panel: n as a function of x and y as averaged over the first two layers. The bar on the right explains the shades of blue used for n. The topmost five chains of the cellulose illustrate the topography of the surface.

the directions of the polarization vectors are much more random than near ZnO, ZnS, and gold,14−16 despite the mildly hydrophilic character of the cellulose surface. Nevertheless, the maximum near Pz/P = −0.5 means that most of the water molecules in the first layer tend to turn their hydrogen atoms toward the surface in order to form hydrogen bonds with the oxygen atoms in the hydroxyl groups of the cellulose. We also observe a tendency for the water dipoles to align antiparallel to the axis of the chain. Wettability. An assessment of the wettability can be obtained by determining the contact angle that a droplet of water makes with the surface. An example of an equilibrated state of the droplet on the surface of the cellulose is shown in Figure 4. Like in ref 16, we have used the procedure developed by Ruijter et al.55 and adopted by Werder et al.56 First, the Cartesian space coordinates (x,y,z) of the droplet atoms are transformed to the planar coordinates (r,z), where r is the distance from the z axis passing through the center of mass of the droplet. Second, the r−z-plane is covered by a fine rectangular grid of points at which the time-averaged number density of atoms is evaluated. Third, a sigmoidal function is fitted to the resulting density profiles in each horizontal slice of the grid points. The average location of the boundary of the droplet in the slice is defined as the position of the center of the sigmoidal function. The top and bottom portions of the droplet are discarded in these fits. Finally, a circle is fitted to these points and the contact angle between the surface and a tangent to the circle is measured. We obtain 24° ± 2°, which is consistent with the experimental finding that a pure cellulose is hydrophilic and produces a contact angle around 20−30°.57 This value of the angle signifies that the cellulose is less hydrophilic than clean

them as I and II. The density in layer I is about 25% higher than the bulk density, which suggests existence of some hydrophilic properties. However, the degree of hydrophilicity is much weaker than in the presence of ZnO,14 for which the first layer nearly doubles the density of the bulk. At the same time, it is stronger than for ZnS, for which only the second layer exhibits some excess density relative to the bulk.15 The small error bars on the values in the density profile indicate that the general structure of water does not change much during the 10 ns of MD. Neither does the structure of cellulose surface, indicating the existence of a fairly ordered water−cellulose interface, which reflects the hydrophilic character of the cellulose surface. The bottom panel of Figure 2 shows the average density of water in the first two layers as viewed across the x−y-plane of the crystal. We notice that the molecules of water favor adsorption to the spots where hydroxyl groups of glucose are exposed, i.e., between the chains. As revealed by Bellesia et al.,54 alteration of dihedral angles involving these groups triggers a transition from cellulose Iβ to cellulose III, a polymorph with ≈5-fold higher enzymatic degradation rates. 24407

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

gold (0°)16 but much more hydrophilic than contaminated gold, for which the calculated contact angle is 150°.16 It should be noted that other surfaces of cellulose are expected to have different wettability properties and, in particular, to be hydrophobic. Amino Acids. Figure 5 shows the PMFs for the 20 standard amino acids as derived both by the MD and REMD methods. The potentials are plotted against the vertical position of the center of mass of the amino acid. The two methods are seen to yield mutually consistent results, but REMD provides better averaging and smoother distributions. The largest discrepancies occur for the three forms of His as well as for the positively charged residues, Lys and Arg. The REMD method is more reliable, as evidenced by the fact that it leads to nearly the same PMF profiles for the two uncharged forms of His (i.e., HIE and HID), unlike the simple MD method. A similar beneficial effect of using the REMD sampling is also seen for Leu and Ile, which are isomers and should behave alike. Many of the PMF distributions have a simple one-well profile, like that obtained for Val. However, there are some with

Figure 4. Final state of a water droplet placed at the surface of cellulose. The snapshot shows only a part of the whole system. The droplet initially has a spherical shape with a diameter of 5 nm. It gradually spreads over the crystal until it reaches equilibrium with a nonzero contact angle. The droplet contains about 2000 molecules of water. The contact angle, α, is determined by fitting an arc (blue line) to the average boundary (the green dots) of the droplet. The contact angle is measured as the slope of a tangent (black line) at the point where the arc intersects the surface.

Figure 5. Potential of mean force for capped amino acids in water as a function of the distance of the center of mass of the residue away from the surface, i.e., away from the center of mass of the topmost layer of the cellulose chains. The MD results are shown in green and those by REMD in blue. The vertical dashed lines delineate the boundaries of the first and second water layers, as defined in Figure 2. 24408

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

description is provided by the free binding energy, as described in the Methods section and listed in Table 1. The way that ϵ depends on the residue is very similar to the behavior of ΔGbind. In our previous papers, we have used the term “binding energy” in reference to ϵ instead of to ΔGbind, which would be more correct. However, to enable making comparisons with our results on other surfaces, we continue to focus on ϵ, which characterizes the PMF near its deepest minimum. The average ϵ is found to be similar to the one determined for contaminated gold,16 but the specificity is almost halved. The specificity parameter is defined here as the ratio of the dispersion to the mean in the set of the values of ϵ for the 20 amino acids. This is because this strongly hydrophobic surface of gold favors adsorption of nonpolar amino acids, whereas the surface of cellulose is more amphiphilic. Since cellulose comes with both hydrophobic (on the side of a chain) and hydrophilic (on the top of a layer) sites, the nonpolar Ala and Leu and the polar Asn and Cys are seen to bind with similar energies of about 7 kJ/mol. We note that the three highest binding energies were determined for the aromatic amino acids, i.e. Tyr, Phe, and Trp. On the other hand, the three weakest binders are the negatively charged Asp and the positively charged HIP and Lys. However, charged Arg binds well, which illustrates the amphiphilicity of cellulose. Figure 6 shows examples of the bound conformations of selected amino acids as obtained through the conventional MD umbrella sampling method. In these simulations, the amino acids are constrained to be at the distance z ≅ σ. (The figure shows only a single conformation for each of the amino acids, whereas the binding energy is an average obtained over 40 000 states.) The conformations shown represent a variety of different types of interactions with the cellulose. In the case of Asp, adsorption is screened by the first layer of water, which leads to a small value of ϵ. The conformation shown for Lys is at the second (higher) minimum in PMF (see Figure 5). We observe that the binding here indeed affects the cap and not the side chain. Leu (a single minimum) and Asn bind both through the cap and through the side chain (nonpolar for Leu and polar for Asn). Gln is an example of an amino acid that can bind to the cellulose both with hydrophobic and hydrophilic groups at the same time, in a striking illustration of the heterogeneous properties of the cellulose. Despite the roughness of the surface, aromatic amino acids can efficiently adhere their rings to flat regions, as illustrated by the two projections of the conformation shown for Phe. This explains the relatively strong adhesion of the aromatic amino acids. Examples of conformations of all amino acids in three different projections are shown in Figures 1−4 in the Supporting Information (SI). Dipeptides and Tripeptides. Figure 7 shows the REMDderived PMF for 22 selected dipeptides. The binding energies (ϵ) and optimum separations (σ) are summarized in Table 2. Examples of near optimal conformations are shown in Figures 6−8 of the SI. It is seen that the ϵ for any dipeptide is always smaller than the sum of the ϵ values corresponding to the constituting amino acids. This is because two side chains cannot be simultaneously optimal with respect to the surface because of their mutual interactions. When one plots the dipeptide binding energy vs the sum of the single amino acid energies, one gets a reasonable correlation with the effective slope of 54%, as shown in the left panel of Figure 8. This is close to the 57% slope obtained for hydrophobic gold,16 suggesting some near-universality.

more pronounced undulations, like Glu or Tyr. These more complex energy landscapes reflect changes in the typical conformations that seem to be coupled with the presence of the adsorbed layers of water. For instance, the local maxima observed, e.g., for Glu, Thr, Tyr, Ala, and Arg, are at z = 0.6 nm, which corresponds to the dip in the water density profile, between layers I and II. The values of ϵ and σ obtained by REMD are listed in Table 1. Almost all values of σ are around 0.4−0.5 nm, which Table 1. Values of the Binding Energy (ϵ) and the Optimum Protein−Cellulose Separation (σ) between the Center of Mass of an Amino Acid and the Surface of Cellulose, As Determined by the REMD Simulations, and the the Free Binding Energy (ΔGbind)a amino acid Asp Glu Cys Asn Phe Thr Tyr Gln Ser Met Trp Val Leu Ile Gly Ala Pro HIE HID HIP Lys Arg av ϵ dispersion in ϵ % specificity

ϵ (kJ/mol)

σ (nm)

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.7 0.5 0.4 0.4 0.5 0.5 0.4 0.5 0.5 0.4 0.4 0.5 0.5 0.5 0.4 0.4 0.4 0.4 0.4 0.6 0.7 0.5 0.5 0.1 17

3.2 4.9 6.5 7.1 9.4 4.7 10.3 5.4 5.2 8.3 8.9 6.4 6.8 6.4 4.6 6.7 7.4 6.4 5.9 3.9 3.7 6.0 6.4 1.8 28

0.5 1.0 0.9 0.6 1.1 0.5 0.8 1.0 0.4 1.0 1.3 0.6 0.7 0.6 0.5 0.4 0.9 0.8 1.0 0.4 0.4 1.3

ΔGbind (kJ/mol) −2.1 −3.7 −4.7 −4.9 −5.5 −2.9 −6.1 −3.9 −3.8 −5.3 −4.8 −4.1 −3.8 −3.7 −4.2 −4.7 −4.4 −4.4 −4.1 −3.1 −3.1 −3.6 −4.2 0.9 22

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.1 0.1 0.1 0.1 0.2 0.1 0.2 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

a The average and dispersion, calculated as the standard deviation, of ϵ, σ, and ΔGbind are listed at the bottom of the table. Results corresponding to the various forms of histidine are first averaged to form one entry when determining the overall averages.

indicates that the amino acids adsorb to the cellulose directly. However, Asp and HIP stay further away. Lys has two comparable minima, close to the surface and further away from the surface, but the deeper minimum corresponds to the latter and the higher minimum is artificial, as it is due to binding of the cap instead of the amino acid. For more hydrophilic surfaces like, say, ZnO,14 the first dense water layer intervenes and prevents direct attachment of all amino acids. Similar “nonfouling” behavior has been measured experimentally and confirmed using molecular dynamics simulations for organic monolayers containing polar OH groups, grown on silica58 and graphene59 surfaces. Our results suggest that cellulose is more likely to be “fouled” in a biological environment. The listed values of ϵ and σ pertain to the deepest minima in PMF. Whenever several minima are present, a more accurate 24409

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

The Amino Acid Analogues. We now consider molecules that are the analogues of the side chains in amino acids. They are listed in Table 3 together with their corresponding amino acid. Four of these amino acids have two analogues each, as explained in more details in ref 16. The plots of the PMF distributions for the analogues are shown in Figure 9 in the SI. They display more undulations than the amino acids because, being smaller molecules, they are more sensitive to the layering of density in water. The REMD-based results for the analogues are summarized in Table 3 together with the data for the corresponding amino acids. They are also plotted in Figure 10 in the SI. Examples of the near optimal conformations are shown in Figures 13 and 14 in the SI. Note that the differences between values of ϵ for amino acids and their analogues vary between 0.9 and 4.2 kJ/ mol. Despite the differences, Phe and its aromatic analogue, benzene, are recognized as the strongest binders. The second largest ϵ found for guanidinium is consistent with the relatively high value found for Arg. This can be explained by the high rigidity and flatness of the guanidinium group. The adhesion of the group does not decrease the entropy much and the size of the group is sufficiently small to adjust to the roughness of the cellulose surface. The adsorption of both benzene and guanidinium appears to be driven by the hydrophobic effect, but the net charge on guanidinium reduces the effectiveness of the process. Comparable energies to Arg were found for Ala and Ser, but the values of ϵ for their analogues are smaller by 50% because the side chains of Ala and Ser are smaller than those of Phe and Arg, which enhances the role of the caps. The other two analogues, butylammonium and ethanoate (both flexible and charged), can be considered as behaving consistently with the ranking for the corresponding amino acids, since their ϵ values agree within the error bars. Ammonium (an analogue of Lys) and methanoate (Asp) are found to have the weakest binding. Even water binds better. This observation confirms that the charged and nonplanar groups do not generate an affinity to the surface of the cellulose. Participation of the caps in the adsorption process of single amino acids does not necessarily mean that the results for the analogues are more reliable. One reason for this is that the amino acid backbone (−NH2−Cα−CO−) can itself participate in binding and, therefore, should be included in the simulations. Another is that the backbone, together with the caps (CH3− CO− and −NH2−CH3), enforce such orientations of side chains that are appropriate for a protein. Protein. We now consider a tryptophan cage in the vicinity of the cellulose surface. Its structure is that of a hairpin with one branch forming the α-helix (sites 2−8). There is also a turn at site 10 and the 3/10 helix at sites 11−14. As explained in the Methods section, we consider two MD trajectories that are followed by the REMD sampling. In the MD part, the protein gets to the surface through an unbiased diffusion. Even though tryptophan cage has 20 sites, it is made of only 12 different amino acids: Pro (four sites), Ser (three), Gly (three), Leu (two), and single copies of Asn, Tyr, Ile, Gln, Trp, Lys, Asp, and Arg. The question we ask is, which of these amino acids are most likely to drive the adsorption events? An example of such an event is shown in Figure 10 (obtained through the standard MD simulation). The top panel plots the time-varying vertical distances of, top to bottom, the highest atom, center of mass of the highest amino acid, center of mass

Figure 6. Snapshots of selected amino acids adsorbed at the surface of cellulose in the conventional MD simulation, as projected to the y−zplane. The inset in the panel for Phe corresponds to a projection to the x−z-plane. The color convention is like in Figure 1 (the oxygen, hydrogen, and carbon atoms correspond to spheres in red, white, and cyan) and the additional blue is for the nitrogen atoms. Amino acids are shown with the N- and C-terminus caps. The atoms of the amino acids are circled and those of the caps are not. The disconnected spheres indicate locations of the ions; Cl− is in green and Na+ is in orange. The amino acids are shown against the background of the water molecules. The bottom chains of cellulose are restrained during the simulations.

It is clear then that the net binding energy of a polypeptide is not a simple addition of the individual amino acid binding energies. Indeed, when the values of ϵ corresponding for six tripeptides are plotted against the sum of the individual well depths, one gets a statistical 43% reduction, as evidenced in Figure 9. Longer chains should lower the effective slope still further, perhaps in a weaker and weaker way. On the other hand, if the peptide sequence is a part of a larger protein, then its orientation is influenced by additional contacts with other segments of the protein. These surface-unrelated interactions are expected to reduce the opportunities for an optimal adjustment to the surface significantly. It is also important to note that the surface of the cellulose is highly heterogeneous and that one residue may be attracted in one region but be repelled in another, but an optimal adjustment may be again ruled out because of the chain connectivity. Perhaps, in a protein, the effective single amino acid binding energy could be as low as, say, 20% of the value determined here. Further results pertaining to the tripeptides can be found in the SI, where Figure 5 provides results on PMF and Figure 8 shows examples of the conformations. We just mention here that poly-Asp behaves differently than other chains in that the longer it is, the poorer the binding. This is because a single Asp binds merely by one of its caps, whereas its side chain is repelled. On increasing the chain length, the role of the caps is waning. 24410

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

Figure 7. Data points in black give the PMF for the selected dipeptides as indicated. The PMFs for their constituent amino acids are shown in green (the first component) and red (the second component if different from the first). They are taken from Figure 5. The vertical dashed lines represent the boundaries of the water layers.

The protein is slightly deformed relative to its native state. Even though the radius of gyration (Rg) for the protein at the surface, 0.62 ± 0.04 nm, is close to Rg of 0.60 ± 0.03 in bulk water, the protein becomes more elongated and thus less globular. This distortion can be characterized by a dimensionless parameter (w) defined in terms of the three main radii of gyration derived from the tensor of inertia,17,18,60 so that w = 0 corresponds to a perfectly spherical symmetry. We find w in bulk water to be 0.035 ± 0.085 and that at the surface to be 0.133 ± 0.081 . Further insight into the protein−cellulose interactions can be obtained through a statistical analysis of many adsorption events. This can be accomplished by switching to the REMDbased simulations. Figure 12 shows the combined outcome of the last 10 ns in two MD-switched-to-REMD simulations. We observe that the average information about the heights of the lowest atoms in the amino acids is quite close to what has been found for the single adsorption event (Figure 10), with the enhanced sampling of REMD more quickly finding additional binding modes than conventional MD. Another way to characterize adsorption is provided by Figure 13, which shows the percentage of the total adsorption time

of the protein, center of mass of the lowest amino acid, and the lowest atom. The identity of the lowest/highest atoms/amino acids typically varies in time. The bottom panel shows the vertical positions of the lowest atoms in the individual amino acids. It identifies residues 12−18 as being responsible for the adsorption, together with 3-Tyr and 6-Trp in the α-helix, which were found to have the largest value of ϵ in Table 1 (among the amino acids that are present in the sequence). Interestingly, Gly and Ser in segment 12−18 have an ϵ that is smaller than the average (6.3 kJ/mol), Arg is close to the average, and the three prolines have an ϵ that is larger. The fourth proline, 19-Pro, is involved in this event much more weakly. Notice that the helix got attached only by two strongly binding amino acids, whereas the elastic fragment was absorbed through many amino acids with the values of ϵ that are in the medium range. In the example shown, adsorption sets in already around 1 ns (in the other trajectory it takes place around 7 ns). One of the corresponding adsorbed conformations is shown in Figure 11 (obtained in REMD simulations). The C-terminal part is seen to be the closest to the surface, whereas the structured part, i.e., the α-helix, is exposed to the bulk water. 24411

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C Table 2. ϵ, σ, and ΔGbind, as in Table 1, but for Dipeptidesa dipeptides Asp-Asp Asp-Gly Asp-Leu Asp-Phe Asp-Ser Gly-Gly Gly-Ser Leu-Gly Leu-Leu Leu-Phe Lys-Asp Lys-Gly Lys-Leu Lys-Lys Lys-Phe Lys-Ser Phe-Gly Phe-Phe Ser-Gly Ser-Leu Ser-Phe Ser-Ser av ϵ dispersion in ϵ a

ϵ (kJ/mol)

σ (nm)

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.7 0.6 0.7 0.6 0.6 0.4 0.5 0.5 0.7 0.6 0.8 0.5 0.6 0.6 0.6 0.6 0.5 0.5 0.5 0.5 0.5 0.5 1.9 0.1

2.6 3.9 4.7 7.2 5.6 7.0 6.7 6.2 6.3 8.6 2.9 6.7 5.2 3.6 7.8 3.1 8.1 9.9 6.8 7.3 7.7 7.2 6.1 0.6

0.5 0.7 0.6 0.9 0.6 0.7 0.9 1.4 0.6 1.0 0.6 1.2 0.9 1.0 0.9 0.8 1.1 1.7 0.9 1.3 1.4 1.0

ΔGbind (kJ/mol) −0.9 −2.6 −2.8 −4.6 −3.6 −4.7 −4.1 −3.8 −3.7 −5.3 −0.8 −3.8 −2.7 −1.7 −4.7 −1.3 −5.2 −6.5 −4.3 −4.5 −5.2 −4.7 −3.7 1.5

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.2 0.1 0.1 0.1 0.2 0.1 0.1 0.2 0.1 0.2 0.2 0.1 0.1 0.2 0.1

Figure 9. Comparison of the binding energy of the tripeptides (ϵijk) and the sums of their constituent single amino acids (ϵi, ϵj, and ϵk) to the surface of cellulose. The scaling factor of the slope of the straightline fits to the data (least-squares method) is 43%.

by the water layer. 16-Arg goes through layer I and is in close contact with the surface. However, the center of mass of 16-Arg stays in layer II. 2-Leu and 3-Tyr exhibit only transient adsorption, whereas 6-Trp, 13-Ser, and 14-Ser stay at the surface most of the time. The average binding energy of amino acids for cellulose is very similar to that calculated for hydrophobic gold, 6.4 and 6.9 kJ/mol respectively.16 Among the 12 distinct amino acids of the tryptophan cage, only four differ by more than 3 kJ/mol between both surfaces: Asn (7.1 and 3.6 kJ/mol for cellulose and gold respectively), Ile (6.4 and 10.9), Leu (6.8 and 11.7), and Trp (8.9 and 12.9). Interestingly, binding of tryptophan cage to the hydrophobic gold occurs quite differently than to the cellulose, as the gold displays a stronger affinity for the αhelix. This affinity can be characterized by determining the

The data were obtained by the REMD method.

when the individual amino acids of the protein are in contact with the surface of cellulose. The contact can be defined either by the lowest atom or by the center of mass of the amino acidthe uniform or hatched bins in the figurestaying withing a given distance from the surface that can be discretized through boundaries of water layers. Here, we distinguish staying either in layer I or layer II, the darker and lighter bins, respectively. Irrespective of the definition used, prolines 12, 17, and 18 and 15-Gly are seen to stay at the surface most often. So does 11-Gly, but this residue never enters layer I; the interaction of its charged side chain with the surface is buffered

Figure 8. Comparison of the binding energy of the dipeptides (ϵij) and the sums of their constituent single amino acids (ϵi and ϵj) to the surface of cellulose (left) and gold (right).16 The scaling factors of the slope of the straight-line fits to the data (least-squares method) are 54% and 57%. 24412

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

0.092 ± 0.059 and 0.084 ± 0.048 in the helical and remaining segments, respectively. The results on f i displayed in Figure 13 allow us to estimate the fraction of the total adsorption time that the helical and remaining segments are in contact with the surface of the cellulose. We obtain 22% for the α-helix and 78% for the rest. A similar analysis done for gold yields the opposite result, 58% and 42%, respectively.

Table 3. Values of ϵ and σ for the Analogues and Their Corresponding Amino Acidsa analogue methane ammonium butylammonium methanol methanoate ethanoate guanidinium benzene water

ϵ (kJ/mol)

σ (nm)

amino acid

ϵ (kJ/mol)

σ (nm)

± ± ± ± ± ± ± ± ±

0.4 1.1 0.5 0.4 0.7 0.4 0.4 0.4 0.4

Ala Lys

6.7 ± 0.4 3.7 ± 0.4

0.4 0.7

Ser Asp

5.2 ± 0.4 3.2 ± 0.5

0.5 0.7

Arg Phe

6.0 ± 1.3 9.4 ± 1.1

0.5 0.5

2.5 0.2 1.8 2.4 0.5 2.2 5.1 7.8 1.1

0.3 0.3 0.4 0.3 0.3 0.3 0.3 0.4 0.3



FINAL COMMENTS Our results on protein adsorption are consistent with much of the affinity ranking of the single amino acids, but the ranking is not the only factor in understanding adsorption. The average value of 6.3 kJ/mol of the binding energies with the specificity of 28% for 20 single amino acids (see Table 1) does not seem to be sufficiently different from that of hydrophobic gold to make a significant difference in the adsorption pattern, yet substantially different binding modes are observed. We have demonstrated that binding energies of the amino acids in chains undergo an effective reduction in strength compared to the energies obtained for single amino acids. The shape of the protein and the flexibility of the surface also contribute to the adjustments in the affinity to bind. We have confirmed the major role of the aromatic residues in directing binding of proteins to the cellulose.29−31 Potentials of mean force can be used further to construct coarse-grained

a

For benchmarking, we also provide these parameters for the water molecule. The distances are measured between the center of mass of an analogue and the surface of cellulose. The analogues of entire side chains are in bold type, whereas the analogues of certain groups in the side chain are in lightface.

RMSF values for the non-hydrogen side chain atoms. We get 0.087 ± 0.061 and 0.094 ± 0.048 in the α-helix and in the rest of protein, respectively. The difference may have to do with the flat and rigid character of the surface of gold used in the studies. The rough and dynamic surface of the cellulose prefers the less ordered regions of the tryptophan cage. The RMSF values are

Figure 10. (Top) Plot of characteristic vertical positions in tryptophan cage during an adsorption event as derived through MD. The lines are explained in the text. (Bottom) Vertical positions of the lowest atoms of the individual amino acids corresponding to the same trajectory. 24413

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

allowed for a thorough statistical characterization of the adsorption events. This is hard to do for large proteins, for example, carbohydrate-binding modules that are involved in cellulose deconstruction. Nevertheless, the discovered validity of the affinity ranking of the single amino acids for tryptophan cage is expected to carry over to the larger proteins as well and should be of use in the interpretation of adsorption of these proteins. For large proteins, it would be interesting to explore the role of twists occurring in cellulose microfibrils.61 The biomolecules studied here are too small to be affected by the twists. There are several new physical insights obtained through our studies. One is that there are many protein adsorption modes to cellulose because of a moderate specificity in binding amino acids and the mildly hydrophilic character of the surface. These modes are driven by the aromatic side groups, which are not as important for adsorption to gold, if the gold surface is effectively hydrophobic with a comparable level of specificity. Another is that there seems to be a general systematic effective weakening of the binding strength of the amino acids as one transits from amino acids to proteins through dipeptides, tripeptides, etc.



Figure 11. Snapshot of adsorbed protein. The ions Na+ and Cl− are in the orange and lemon colors, respectively. The Cα atoms are enhanced by the black coats. The backbone of the protein is in ribbon representation. The water molecules are partially transparent.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b07118. Examples of near optimal conformations of the capped amino acids, dipeptides, tripeptides, analogues of amino acid side chains, and the molecule of water near the (100) surface of cellulose; potentials of the mean force for selected tripeptides and the resulting values of the binding energies; a comparison between the binding energies of the capped amino acids and of the corresponding analogues of their side chains; results of a similar analysis for dipeptides and tripeptides; and potentials of the mean force for arginine at two temperatures, 300 and 360 K (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Figure 12. Similar to the bottom panel of Figure 10, but the data shown here are obtained by combining two REMD-derived trajectories.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This research has been supported by the ERA-NET grant ERAIB (EIB.12.022) (FiberFuel) and the European Framework

models to study adsorption of large proteins to the surface of cellulose in implicit water. Our choice of the protein has

Figure 13. REMD-based results (the combination of the two runs) on the fraction of time that an individual amino acid of tryptophan cage bonds to the surface of cellulose. The fraction is relative to the total adsorption time of the protein. The darker colors (green and red) correspond to the amino acids staying in layer I, as defined in Figure 2, whereas the lighter colors (lime green and orange) correspond to the amino acids staying either in layer I or in layer II. The uniformly shaded bins (on the left side of each entry) are used when the positions of the amino acids are defined through their lowest atoms. The hatched bins (on the right side) are used for positions defined through the center of mass of the amino acids. 24414

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C

(23) Srinivas, G.; Cheng, X.; Smith, J. C. A Solvent-Free Coarse Grain Model for Crystalline and Amorphous Cellulose Fibrils. J. Chem. Theory Comput. 2011, 7, 2539−2548. (24) Matthews, J. F.; Skopec, C. E.; Mason, P. E.; Zuccato, P.; Torget, R. W.; Sugiyama, J.; Himmel, M. E.; Brady, J. W. Computer Simulations Studies of Microcrystalline Cellulose Iβ. Carbohydr. Res. 2006, 341, 138−152. (25) Bayer, E. A.; Lamed, R.; Himmel, M. E. The Potential of Cellulases and Cellulosomes for Cellulosic Waste Management. Curr. Opin. Biotechnol. 2007, 18, 237−245. (26) Bayer, E. A.; Belaich, J. P.; Shoham, Y.; Lamed, R. The Cellulosomes: Multi-Enzyme Machines for Degradation of Plant Cell Wall Polysaccharides. Annu. Rev. Microbiol. 2004, 58, 521−554. (27) Israelachvili, J. N., Ed. Intermolecular and Surface Forces, 2nd ed.; Academic Press: London, 1992; pp484−489. (28) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3, 300−313. (29) Payne, C. M.; Bomble, Y. J.; Taylor, C. B.; McCabe, C.; Himmel, M. E.; Crowley, M. F.; Beckham, G. T. Multiple Functions of Aromatic-Carbohydrate Interactions in a Processive Cellulase Examined with Molecular Simulation. J. Biol. Chem. 2011, 286, 41028−41035. (30) Tormo, J.; Lamed, R.; Chirino, A.; Morag, E.; Bayer, E. A.; Shoham, Y.; Steitz, T. A. Crystal Structure of Bacterial Family-III Cellulose-Binding Domain: A General Mechanism for Attachment to Cellulose. EMBO J. 1996, 15, 5739−5751. (31) Lehtio, J.; Sugiyama, J.; Gustavsson, M.; Fransson, L.; Linder, M.; Teeri, T. T. The Binding Specificity and Affinity Determinants of Family 1 and Family 3 Cellulose Binding Modules. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 484−489. (32) Caso, J. V.; Russo, L.; Palmieri, M.; Malgieri, G.; Galdiero, S.; Falanga, A.; Isernia, C.; Iacovino, R. Investigating the Inclusion Properties of Aromatic Amino Acids Complexing Beta-Cyclodextrins in Model Peptides. Amino Acids 2015, 47, 2215−2227. (33) Banik, A.; Saikia, M. D. Adsorptive Interaction of Chiral Amino Acids on β-Cyclodextrin Bonded to Silica Particles. J. Encapsulation Adsorpt. Sci. 2013, 3, 35−47. (34) Gomes, T. C. F.; Skaf, M. S. Cellulose-Builder: A Toolkit for Building Crystalline Structures of Cellulose. J. Comput. Chem. 2012, 33, 1338−1346. (35) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (36) Dragneva, N.; Floriano, W. B.; Stauffer, D.; Mawhinney, R. C.; Fanchini, G.; Rubel, O. Favorable Adsorption of Capped Amino Acids on Graphene Substrate Driven by Desolvation Effect. J. Chem. Phys. 2013, 139, 174711−174717. (37) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (38) Huang, J.; MacKerell, A. D. CHARMM36 All-Atom Additive Protein Force Field: Validation Based on Comparison to NMR Data. J. Comput. Chem. 2013, 34, 2135−2145. (39) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (40) Mackerell, A. D.; Feig, M.; Brooks, C. L. Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comput. Chem. 2004, 25, 1400−1415. (41) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E.; Mittal, J.; Feig, M.; Mackerell, A. D. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ, and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257−3273. (42) Guvench, O.; Greene, S. N.; Kamath, G.; Brady, J. W.; Venable, R. M.; Pastor, R. W.; Mackerell, A. D. Additive Empirical Force Field

Programme VII NMP grant 604530-2 (CellulosomePlus). It was cofinanced by the Polish Ministry of Science and Higher Education from the resources granted for the years 2014−2017 in support of international scientific projects. We thank Science Foundation Ireland (SFI) for computing time at the SFI/ Higher Education Authority Irish Center for High-End Computing (ICHEC).



REFERENCES

(1) Gray, J. J. The Interaction of Proteins with Solid Surfaces. Curr. Opin. Struct. Biol. 2004, 14, 110−115. (2) Nakanishi, K.; Sakiyama, T.; Imamura, K. On the Adsorption of Proteins on Solid Surfaces, a Common but Very Complicated Phenomeneon. J. Biosci. Bioeng. 2001, 91, 233−244. (3) Feng, L.; Andrade, J. D. Protein Adsorption on Low-Temperature Isotropic Carbon: I. Protein Conformational Change Probed by Differential Scanning Calorimetry. J. Biomed. Mater. Res. 1994, 28, 735−743. (4) Haynes, C. A.; Norde, W. Structures and Stabilities of Adsorbed Proteins. J. Colloid Interface Sci. 1995, 169, 313−328. (5) Whaley, S. R.; English, D. S.; Hu, E. L.; Barbara, P. F.; Belcher, A. M. Selection of Peptides with Semiconductor Binding Specificity for Directed Nanocrystal Assembly. Nature 2000, 405, 665−668. (6) Cuif, J. P.; Dauphin, Y.; Sorauf, J. E. Biominerals and Fossils Through Time; Cambridge University Press: Cambridge, 2011. (7) Smeets, P. J. M.; Cho, K. R.; Kempen, R. G. E.; Sommerdijk, N. A. J. M.; De Yoreo, J. J. Calcium Carbonate Nucleation Driven by Ion Binding in a Biomimetic Matrix Revealed by In Situ Electron Microscopy. Nat. Mater. 2015, 14, 394−399. (8) Olson, I. C.; Kozdon, R.; Valley, J. W.; Gilbert, P. Mollusk Shell Nacre Ultrastructure Correlates with Environmental Temperature and Pressure. J. Am. Chem. Soc. 2012, 134, 7351−7358. (9) Bachmann, M.; Goede, K.; Beck-Sickinger, A. G.; Grundmann, M.; Irback, A.; Janke, W. Microscopic Mechanisms of Specific Peptide Adhesion to Semiconductor Substrates. Angew. Chem., Int. Ed. 2010, 49, 9530−9533. (10) Walling, M. A.; Novak, J. A.; Shepard, J. R. E. Quantum Dots for Live Cell and In Vivo Imaging. Int. J. Mol. Sci. 2009, 10, 441. (11) Frasco, M. F.; Chaniotakis, N. Semiconductor Quantum Dots in Chemical Sensors and Biosensors. Sensors 2009, 9, 7266. (12) Lynch, I.; Dawson, K. A. Protein-Nanoparticle Interactions. Nano Today 2008, 3, 40−47. (13) Somorjai, G. A. Introduction to Surface Chemistry and Catalysis; John Wiley & Sons, Inc.: New York, 1994. (14) Nawrocki, G.; Cieplak, M. Amino Acids and Proteins at ZnOWater Interfaces in Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2013, 15, 13628−13636. (15) Nawrocki, G.; Cieplak, M. Interactions of Aqueous Amino Acids and Proteins with the (110) Surface of ZnS in Molecular Dynamics Simulations. J. Chem. Phys. 2014, 140, 095101−095111. (16) Nawrocki, G.; Cieplak, M. Aqueous Amino Acids and Proteins Near the Surface of Gold in Hydrophilic and Hydrophobic Force Fields. J. Phys. Chem. C 2014, 118, 12929−12943. (17) Starzyk, A.; Cieplak, M. Denaturation of Proteins Near Polar Surfaces. J. Chem. Phys. 2011, 135, 235103−235112. (18) Starzyk, A.; Cieplak, M. Proteins in the Electric Field Near the Surface of Mica. J. Chem. Phys. 2013, 139, 045102−045110. (19) Frey-Wyssling, A. The Fine Structure of Cellulose Microfibrils. Science 1954, 119, 80−82. (20) O'Sullivan, A. C. Cellulose: The Structure Slowly Unravels. Cellulose 1997, 4, 173−207. (21) Nishiyama, Y.; Langan, P.; Chanzy, H. Crystal Structure and Hydrogen-Bonding System in Cellulose Iβ from Synchroton X-ray and Neutron Fiber Diffraction. J. Am. Chem. Soc. 2002, 124, 9074−9082. (22) Shen, T.; Gnanakaran, S. The Stability of Cellulose: A Statistical Perspective from a Coarse-Grained Model of Hydrogen-Bond Networks. Biophys. J. 2009, 96, 3032−3040. 24415

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416

Article

The Journal of Physical Chemistry C for Hexopyranose Monosaccharides. J. Comput. Chem. 2008, 29, 2543−2564. (43) Guvench, O.; Hatcher, E. R.; Venable, R. M.; Pastor, R. W.; Mackerell, A. D. CHARMM Additive All-Atom Force Field for Glycosidic Linkages Between Hexopyranoses. J. Chem. Theory Comput. 2009, 5, 2353−2370. (44) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Potential. J. Chem. Phys. 1995, 103, 8577−8592. (45) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101− 014107. (46) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (47) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for FreeEnergy Calculations of Biomolecules. J. Comput. Chem. 1992, 13, 1011. (48) Lemkul, J. A.; Bevan, D. R. Assessing the Stability of Alzheimer’s Amyloid Protofibrils Using Molecular Dynamics. J. Phys. Chem. B 2010, 114, 1652−1660. (49) Wolf, M. G.; Jongejan, J. A.; Laman, J. D.; de Leeuw, S. W. Rapid Free Energy Calculation of Peptide Self-Assembly by REMD Umbrella Sampling. J. Phys. Chem. B 2008, 112, 13493−13498. (50) Efron, B. Bootstrap Methods: Another Look at the Jackknife. Ann. Stat. 1979, 7, 1−26. (51) de Ruiter, A.; Zagrovic, B. Absolute Binding-Free Energies between Standard RNA/DNA Nucleobases and Amino-Acid Sidechain Analogs in Different Environments. Nucleic Acids Res. 2015, 43, 708− 718. (52) Cieplak, M.; Koplik, J.; Banavar, J. R. Applications of Statistical Mechanics in Subcontinuum Fluid Dynamics. Phys. A 1999, 274, 281− 293. (53) Cieplak, M.; Koplik, J.; Bavanar, J. R. Molecular Dynamics of Flows in the Knudsen Regime. Phys. A 2000, 287, 153−160. (54) Bellesia, G.; Chundawat, S. P. S.; Langan, P.; Redondo, A.; Dale, B. E.; Gnanakaran, S. Coarse-Grained Model for the Interconversion between Native and Liquid Ammonia-Treated Crystalline Cellulose. J. Phys. Chem. B 2012, 116, 8031−8037. (55) de Ruijter, M. J.; Blake, T. D.; De Coninck, J. Dynamic Wetting Studied by Molecular Modeling Simulations of Droplet Spreading. Langmuir 1999, 15, 7836−7847. (56) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. On the Water-Carbon Interaction for Use in Molecular Dynamics Simulations of Graphite and Carbon Nanotubes. J. Phys. Chem. B 2003, 107, 1345−1352. (57) Bishop, C. A., Ed. Vacuum Deposition onto Webs, Films, and Foils, 2nd ed.;William Andrew Publishing: Norwich, NY, 2007; p 169. (58) Sheikh, S.; Blaszykowski, C.; Nolan, R.; Thompson, D.; Thompson, M. On the Hydration of Subnanometric Antifouling Organosilane Adlayers: A Molecular Dynamics Simulation. J. Colloid Interface Sci. 2015, 437, 197−204. (59) O’Mahony, S.; O’Dwyer, C.; Nijhuis, C. A.; Greer, J. C.; Quinn, A. J.; Thompson, D. Nanoscale Dynamics and Protein Adhesivity of Alkylamine Self-Assembled Monolayers on Graphene. Langmuir 2013, 29, 7271−7282. (60) Sikora, M.; Szymczak, P.; Thompson, D.; Cieplak, M. LinkerMediated Assembly of Gold Nanoparticles into Multimeric Motifs. Nanotechnology 2011, 22, 445601. (61) Hadden, J. A.; French, A. D.; Woods, R. J. Unraveling Cellulose Microfibrils: A Twisted Tale. Biopolymers 2013, 99, 746−756.

24416

DOI: 10.1021/acs.jpcc.5b07118 J. Phys. Chem. C 2015, 119, 24404−24416