2.1 Elemental metals and alloys

Protein–metal surface interactions can be studied with experimental techniques such as atomic force microscopy (AFM) (Binnig et al., Reference Binnig, Quate and Gerber1986), surface plasmon resonance (SPR) (Jönsson et al. Reference Jönsson, Fägerstam, Ivarsson, Johnsson, Karlsson, Lundh, Löfås, Persson, Roos and Rönnberg1991) and localized SPR (Stuart et al. Reference Stuart, Haes, Yonzon, Hicks and Van Duyne2005). Due to their chemical inertness and unique optical properties (Jain et al. Reference Jain, Huang, El-Sayed and El-Sayed2008), the noble metals, gold and silver, are the most commonly used metals employed as probes or sensors in these techniques. Along with well-known applications of metal surfaces, such as biosensors and implants (Liu et al. Reference Liu, Chu and Ding2004), metal surfaces are also used in bioelectronics as electrodes as they allow controlled exchange of electrons with metalloproteins immobilized on them (Alessandrini et al. Reference Alessandrini, Salerno, Frabboni and Facci2005; Andolfi et al. Reference Andolfi, Bruce, Cannistraro, Canters, Davis, Hill, Crozier, Verbeet, Wrathmell and Astier2004). The interactions of proteins with bare metal surfaces have been the subject of many computational studies, which cover a wide range of different elemental and alloy surfaces, such as Cu(100) (Chen & Wang, Reference Chen, Wang, Huang and Gao2010), Au(111) (Bizzarri, Reference Bizzarri2006; Hoefling et al. Reference Hoefling, Monti, Corni and Gottschalk2011; Siwko & Corni, Reference Siwko and Corni2013; Venkat et al. Reference Venkat, Corni and Di Felice2007; Zanetti-Polzi et al. Reference Zanetti-Polzi, Daidone, Bortolotti and Corni2014), Au(100) (Hagiwara et al. Reference Hagiwara, Sakiyama and Watanabe2009), Au nanoparticle (Todorova et al. Reference Todorova, Chiappini, Mager, Simona, Patel, Stevens and Yarovsky2014), Fe (Zhang et al. Reference Zhang, Chen, Feng, Ding, Zhou and Jia2009b), Ni (Yang & Zhao, Reference Yang and Zhao2006), Pd (Coppage et al. Reference Coppage, Slocik, Briggs, Frenkel, Heinz, Naik and Knecht2011), Pt (Kantarci et al. Reference Kantarci, Tamerler, Sarikaya, Haliloglu and Doruker2005), Ag (Aliaga et al. Reference Aliaga, Ahumada, Sepúlveda, Gomez-Jeria, Garrido, Weiss-López and Campos-Vallette2011; Ghosh et al. Reference Ghosh, Jana and Guchhait2012) and steel (Imamura et al. Reference Imamura, Kawasaki, Awadzu, Sakiyama and Nakanishi2003).

2.2 Oxides and minerals

Metal surfaces (excluding noble metals such as gold, platinum and palladium) are oxidized when exposed to water or air, forming metal oxides that are very common in the Earth's crust. Due to their good mechanical stability, catalytic properties and biocompatibility (Andreescu et al. Reference Andreescu, Ornatska, Erlichman, Estevez, Leiter and Matijević2012), metal oxides and minerals are used in a wide range of applications that include fabrication of biomaterials (Whaley et al. Reference Whaley, English, Hu, Barbara and Belcher2000), cellular delivery of drugs and biomolecules (Kievit & Zhang, Reference Kievit and Zhang2011; Xu et al. Reference Xu, Zeng, Lu and Yu2006), tissue engineering (Shin et al. Reference Shin, Jo and Mikos2003) and proteomics (Sugiyama et al. Reference Sugiyama, Masuda, Shinoda, Nakamura, Tomita and Ishihama2007). Computational studies to investigate the interactions of oxide and mineral surfaces with proteins or peptides have mostly been carried out for different forms of titanium dioxide, such as rutile and anatase (Carravetta et al. Reference Carravetta, Monti and Zhang2009; Kang et al. Reference Kang, Li, Tu, Wang and Ågren2010; Köppen et al. Reference Köppen, Bronkalla and Langel2008; Monti, Reference Monti2007; Monti et al. Reference Monti, Carravetta, Battocchio, Iucci and Polzonetti2008; Sun et al. Reference Sun, Han, Lindgren, Shen and Laaksonen2014a; Wu et al. Reference Wu, Chen, Skelton, Cummings and Zheng2013), silicon dioxides (Chen et al. Reference Chen, Su, Neoh and Choe2009a; Nonella & Seeger, Reference Nonella and Seeger2008; Patwardhan et al. Reference Patwardhan, Emami, Berry, Jones, Naik, Deschaume, Heinz and Perry2012; Rimola et al. Reference Rimola, Sodupe and Ugliengo2009; Tosaka et al. Reference Tosaka, Yamamoto, Ohdomari and Watanabe2010), calcite (Wierzbicki et al. Reference Wierzbicki, Sikes, Madura and Drake1994) and mica (Kang et al. Reference Kang, Huynh, Xia, Zhang, Fang, Wei and Zhou2013).

2.3 Self-assembled monolayers

SAMs are thin films that coat surfaces by spontaneous adsorption of organic molecules that form ordered molecular assemblies. Typically, in a SAM, molecules are chemisorbed onto a surface substrate through their reactive head groups and are, therefore, very stable. SAMs can be categorized into two groups according to their head-group type: thiol-based and silate-based (Schreiber, Reference Schreiber2004). The head-group is attached to a tail group that, by following a fast adsorption phase (seconds), undergoes a slow reorganization (hours), in which interactions with other tail groups increase and the packing is improved (Love et al. Reference Love, Estroff, Kriebel, Nuzzo and Whitesides2005). The tail group can be functionalized with small chemical groups or large molecules, such as peptides. These, together with the length of the tail group, allow the physico-chemical interfacial properties, in particular, hydrophobicity and ionization to be adjusted according to the desired application. Alkanethiols are the most common type of SAM. They have the chemical formula S-(CH₂)_n-R, where R stands for a functional group, such as -CH₃, -COOH, -NH₂ or -OH.

Bare metals and oxide surfaces are prone to non-specific adsorption of proteins and other organic molecules. These adsorption processes may result in undesirable agglomeration of adsorbates on the surfaces. The self-assembly of organic molecules on a metal surface in a SAM creates a physical barrier between the surface and the adsorbates, acting as an electrical insulator and passivating the surface atoms (Love et al. Reference Love, Estroff, Kriebel, Nuzzo and Whitesides2005). The properties of SAMs have been reviewed in (Chaki & Vijayamohanan, Reference Chaki and Vijayamohanan2002; Gooding et al. Reference Gooding, Mearns, Yang and Liu2003; Love et al. Reference Love, Estroff, Kriebel, Nuzzo and Whitesides2005; Senaratne et al. Reference Senaratne, Andruzzi and Ober2005). Modeling of protein–SAM interactions has been reported, mostly for alkanethiol SAMs, in (Alvarez-Paggi et al. Reference Alvarez-Paggi, Martín, DeBiase, Hildebrandt, Martí and Murgida2010; Hsu et al. Reference Hsu, Sheu and Tsay2008; O'Mahony et al. Reference O'Mahony, O'Dwyer, Nijhuis, Greer, Quinn and Thompson2013; Sun et al. Reference Sun, Welsh and Latour2005; Utesch et al. Reference Utesch, Millo, Castro, Hildebrandt, Zebger and Mroginski2013; Wang et al. Reference Wang, Zhao, Zhao, Wang, Yang, Yu and Zheng2010a, Reference Wang, Zhao, Yu, Zhao, Li and Zhengb; Xie et al. Reference Xie, Liu and Zhou2012) and peptide–SAM interactions have been modeled by Nowinski et al. (Reference Nowinski, Sun, White, Keefe and Jiang2012).

2.4 Polymers

Polymer surfaces have attracted much attention in the nanotechnology field due to their mechanical stability, low cost and their wide applicability (Nie & Kumacheva, Reference Nie and Kumacheva2008).

Particularly synthetic polymer-based biomaterials, due to their non-fouling properties, are currently being investigated intensively for applications in controlled drug delivery (Hoffman, Reference Hoffman2008), in highly sensitive biosensors (Anker et al. Reference Anker, Hall, Lyandres, Shah, Zhao and Van Duyne2008), and in bioelectronics (Senaratne et al. Reference Senaratne, Andruzzi and Ober2005). Nanostructured polymer materials used for bio-related and medicinal research include electrostatic polymer brushes, micelles, layer-by-layer deposition and thin films (Stuart et al. Reference Stuart, Huck, Genzer, Müller, Ober, Stamm, Sukhorukov, Szleifer, Tsukruk, Urban, Winnik, Zauscher, Luzinov and Minko2010). Polymer brushes are surface modifiers that share many properties in common with SAMs (Senaratne et al. Reference Senaratne, Andruzzi and Ober2005). They are prepared by grafting polymers of the same or varying kinds on surfaces forming homopolymer and mixed brushes, respectively. Polymeric micelles are formed through self-assembly of amphiphilic copolymers and typically have a diameter of size 30–50 nm (Otsuka et al. Reference Otsuka, Nagasaki and Kataoka2003; Stuart et al. Reference Stuart, Huck, Genzer, Müller, Ober, Stamm, Sukhorukov, Szleifer, Tsukruk, Urban, Winnik, Zauscher, Luzinov and Minko2010).

These micelles are particularly important due to their lower critical micelle concentration (CMC), higher stability and slower rate of dissociation than surfactant micelles. These properties have allowed polymeric micelles to act as effective cancer treatment tools with high drug deposition at the target site (Otsuka et al. Reference Otsuka, Nagasaki and Kataoka2003). Polymer surfaces and their applications have been reviewed by (Barbey et al. Reference Barbey, Lavanant, Paripovic, Schüwer, Sugnaux, Tugulu and Klok2009; Kim et al. Reference Kim, Khang and Lee2008; Nie et al. Reference Nie, Petukhova and Kumacheva2010; Otsuka et al. Reference Otsuka, Nagasaki and Kataoka2003; Senaratne et al. Reference Senaratne, Andruzzi and Ober2005; Stuart et al. Reference Stuart, Huck, Genzer, Müller, Ober, Stamm, Sukhorukov, Szleifer, Tsukruk, Urban, Winnik, Zauscher, Luzinov and Minko2010). The limited number of computational studies of protein-polymer surfaces to date have been carried out for polymer types including polystyrene, polyethylene and polydimethylsiloxane (Boughton et al. Reference Boughton, Andricioaei and Chen2010; Jeyachandran et al. Reference Jeyachandran, Mielczarski, Rai and Mielczarski2009; Liu et al. Reference Liu, Wu, Feng, Shao and Cai2012; Lu et al. Reference Lu, Lee and Park1992; O'Brien et al. Reference O'Brien, Stuart, Bruce and Latour2008; Raffaini & Ganazzoli, Reference Raffaini and Ganazzoli2007; Zhang et al. Reference Zhang, Zhao and Sun2009a; Wei et al. Reference Wei, Carignano and Szleifer2011).

2.5 Carbon allotropes

Pure carbon may exist in a number of different allotropes. Owing to their unique thermal, electrical, chemical and mechanical properties, carbon-based nanomaterials have been the subject of numerous applications in analytical chemistry (Scida et al. Reference Scida, Stege, Haby, Messina and García2011). These applications mostly focus on carbon nanomaterials with sp²-carbon bonding, such as fullerenes, carbon nanotubes (CNT), graphene and graphite. This is due to the extremely high surface areas of fullerenes and CNTs relative to their size, which makes them suitable for design as highly efficient drug carriers. Furthermore, the excellent electrical properties of graphene and CNTs make them suitable for biosensor applications (Liu & Liang, Reference Liu and Liang2012). In all the four materials, the carbon atoms make three chemical bonds with other carbons in the surface-plane with delocalized π electron clouds in the direction perpendicular to the surface (Scida et al. Reference Scida, Stege, Haby, Messina and García2011). This configuration makes the mutual van der Waals interactions between CNTs very strong and hence leads them to be very hydrophobic (Guldi et al. Reference Guldi, Rahman, Sgobba and Ehli2006). To alter the hydrophobicity, surface modifications with surface defects and polar groups have been suggested, but these affect the stability of the materials as well as their mechanical and electrical properties (Scida et al. Reference Scida, Stege, Haby, Messina and García2011). Computational studies of protein–carbon surface interactions have mostly focused on graphene/graphite (Mereghetti & Wade, Reference Mereghetti and Wade2011; Mücksch & Urbassek, Reference Mücksch and Urbassek2011; Raffaini & Ganazzoli, Reference Raffaini and Ganazzoli2010; Kang et al. Reference Kang, Huynh, Xia, Zhang, Fang, Wei and Zhou2013; Sun et al. Reference Sun, Feng, Hou and Li2014b; Yu et al. Reference Yu, Wang, Lin, Zhao, Zhao and Zheng2012b), CNT (Balamurugan et al. Reference Balamurugan, Gopalakrishnan, Raman and Subramanian2010; Chen et al. Reference Chen, Wang, Liu, Wu, Kang, Moore and Gubbins2009b; Tallury & Pasquinelli, Reference Tallury and Pasquinelli2010; Wang et al. Reference Wang, Humphreys, Chung, Delduco, Lustig, Wang, Parker, Rizzo, Subramoney, Chiang and Jagota2003) and fullerenes (Durdagi et al. Reference Durdagi, Mavromoustakos, Chronakis and Papadopoulos2008; Kraszewski et al. Reference Kraszewski, Tarek, Treptow and Ramseyer2010; Noon et al. Reference Noon, Kong and Ma2002).

3.1 Ionization and hydration

SAMs and metal oxide surfaces may be protonated and/or hydroxylated to varying degrees depending on the environmental conditions and on the material itself, e.g. pH, material shape and size. Many theoretical studies on SAMs do not take into account the ionization states of the functional groups of SAMs. However, a recent study (Utesch et al. Reference Utesch, Millo, Castro, Hildebrandt, Zebger and Mroginski2013) reporting a titration curve of an amino-terminated alkanethiol SAM showed that the level of ionization was very sensitive at pH values around 6 (with 16 ± 5% protonation at pH 7 and 52 ± 9% at pH 6). A positive correlation between the ionization level of a SAM and the strength of adsorption of charged proteins was found (Zhou et al. Reference Zhou, Zheng and Jiang2004; Utesch et al. Reference Utesch, Millo, Castro, Hildebrandt, Zebger and Mroginski2013). Therefore, it is crucial to compute and model the surface ionization to obtain reliable results compatible with experiments. The protonation states of SAMs in modeling studies are mostly represented by random assignment of protonated and deprotonated groups (Alvarez-Paggi et al. Reference Alvarez-Paggi, Castro, Tórtora, Castro, Radi and Murgida2013; Utesch et al. Reference Utesch, Sezer, Weidinger and Mroginski2012; Zhou et al. Reference Zhou, Zheng and Jiang2004). In several studies, they are represented either by a uniform distribution of small partial charges (Sun et al. Reference Sun, Welsh and Latour2005) or by large partial charges being assigned to functional groups in neutral surfaces (Wang et al. Reference Wang, Zhao, Yu, Zhao, Li and Zheng2010b).

As with SAMs, SiO₂ and TiO₂ surfaces have different levels of ionization of the surface groups depending on the pH of the environment, and this has been shown to determine selectivity of the adsorption of proteins (Patwardhan et al. Reference Patwardhan, Emami, Berry, Jones, Naik, Deschaume, Heinz and Perry2012). Studies have shown that the concentration of silanol groups (Si-OH) and the degree of their ionization define the hydrophobicity of silicon dioxide (silica) surfaces and govern their adsorption properties and thus also, the behavior of silica-based materials in processes such as biomolecular adhesion and biomineralization (Sumper & Brunner, Reference Sumper and Brunner2008; Voskerician et al. Reference Voskerician, Shive, Shawgo, von Recum, Anderson, Cima and Langer2003). In computational studies of protein–oxide surface interactions, ionization may be represented explicitly (Friedrichs et al. Reference Friedrichs, Köppen and Langel2013; Köppen & Langel, Reference Köppen and Langel2010; Patwardhan et al. Reference Patwardhan, Emami, Berry, Jones, Naik, Deschaume, Heinz and Perry2012; Tosaka et al. Reference Tosaka, Yamamoto, Ohdomari and Watanabe2010), or by assigning uniform partial charges to selected surface atoms (Kubiak-Ossowska & Mulheran, Reference Kubiak-Ossowska and Mulheran2012). Köppen & Langel (Reference Köppen and Langel2010) showed that the adhesion energies of peptides on titanium dioxide surfaces are sensitive to the values of the partial charges of the surface hydroxyl groups. Therefore, care must be taken to ensure a reliable parameterization of the ionization charges.

Conventional simulations of proteins usually neglect changes in the protonation states of ionizable groups as a function of time. However, an accurate simulation of a protein-surface system may require a more sophisticated approach, such as constant pH simulation, to treat the variation of the protonation states of residues in the interfacial region, as well as of the surface interaction sites. Although constant pH simulation techniques have been applied to various molecular systems, including small chemotherapeutic drug molecules binding to nanodiamonds (Adnan et al. Reference Adnan, Lam, Chen, Lee, Schaffer, Barnard, Schatz, Ho and Liu2011), they have not, as far as we are aware, been applied to simulations of peptide/protein–surface interactions.

The most important issue for modeling the adsorption of biomolecules on oxide surfaces is the treatment of the properties of the surface hydration shell. Due to dissociation, the hydroxyl groups and hydrogen atoms form bonds with unsaturated surface metal (e.g. Ti) and O atoms, respectively. The degree of water dissociation defines the physical properties of the surface and strongly affects the binding properties of biological molecules. Kang et al. (Reference Kang, Li, Tu, Wang and Ågren2010) investigated the role of the water in the adsorption process on rutile surfaces and observed that human serum albumin (HSA) adsorbs on a rutile surface modified with -OH groups more strongly than on an unmodified rutile surface. Hydrogen bond analysis in the same study showed that the bonds formed between the structured hydroxyl groups on the modified surface reduced the possibility of hydrogen bond formation between the surface and the water molecules, hence making it easier for the protein to adsorb onto the modified surface.

It should be noted that water dissociation is reversible on all oxides (Henderson, Reference Henderson2002). Therefore, the concentration and position of hydroxyl groups may change with time. Even though relying on the average distribution of hydroxyl groups on a particular surface may be sufficient, it may not always give accurate results, especially if the protein adsorption kinetics are determined by surface hydration kinetics. This poses a problem for conventional modeling techniques for biomolecules, which ideally have to capture the dynamics of dissociative adsorption and associative desorption of water molecules on oxide surfaces.

3.2 Polarization

The electrostatic potential of solutes and solvent molecules induces an attractive polarization of metal surfaces. Metal polarization is negligible for neutral adsorbates that do not have large dipoles and, since proteins usually have a relatively small total charge, it is often neglected in simulations (Braun et al. Reference Braun, Sarikaya and Schulten2002). Early studies of metal surfaces also showed that the effect of pure water on polarization of the metal is often negligible (Barabino & Marchesi, Reference Barabino, Cesare-Gavotti and Marchesi1984; Shelley et al. Reference Shelley, Patey, Bérard and Torrie1997; Spohr, Reference Spohr1995) and that charge-induced polarization does not cause any change in the structure of the water on the surface (Feng et al. Reference Feng, Pandey, Berry, Farmer, Naik and Heinz2011). However, for the adsorption of biomolecules with considerable dipole moments, including some peptides and proteins, the surface polarization contributes to the binding energy and influences the binding mode. Although studies have shown that, for water-metal surface systems, the energy due to polarization is less than 10% of the total binding energy (Feng et al. Reference Feng, Pandey, Berry, Farmer, Naik and Heinz2011; Neves et al. Reference Neves, Motheo, Fartaria and Silva Fernandes2007; Siepmann & Sprik, Reference Siepmann and Sprik1995; Vila Verde et al. Reference Vertegel, Siegel and Dordick2009, Reference Vigmond, Iwakura, Mizutani and Katsura2011), for proteins, in particular, on Au(111) surfaces, the contribution from polarization has been estimated to be about 10–20% of the total binding energy (Heinz et al. Reference Heinz, Jha, Luettmer-Strathmann, Farmer and Naik2011). Further, on surfaces such as Au(100), where the van der Waals attraction is weaker, polarization was found to tune the adsorption of proteins (Heinz et al. Reference Heinz, Farmer, Pandey, Slocik, Patnaik, Pachter and Naik2009) and act as a major contributor to the adsorption of highly charged peptides (Heinz et al. Reference Heinz, Jha, Luettmer-Strathmann, Farmer and Naik2011). For the simulation of amino acid adsorption on Au(111), polarization effects were also found to be important for reproducing experimental binding propensities (Hoefling et al. Reference Hoefling, Monti, Corni and Gottschalk2011).

The polarization of surfaces of other types can also play an important role in the interactions of surfaces with their environments. Schyman & Jorgensen (Reference Schyman and Jorgensen2013) showed that while a non-polarizable FF is adequate for the description of interactions between water and small hydrocarbons, such as benzene (C₆H₆) and coronene (C₂₄H₁₂), a polarizable FF is required for CNTs and fullerenes in order to reproduce interaction energy values obtained from density functional theory (DFT) calculations. Therefore, the effect of polarization should be carefully considered in computer simulations of biomolecules, in particular, of charged proteins, with metal and carbon surfaces.

3.3 Reconstruction

The atomic structure of the surface of a material generally differs from that of its interior because of differences in the forces acting on the atoms in the vicinity of the surface. The type and the degree of reconstruction of a surface is determined by environmental conditions, such as temperature and pressure (Somorjai & Li, Reference Somorjai and Li2011), as well as the structure of the material and may be affected by adsorbing molecules. Ideally, reconstruction of surfaces has to be taken into account in the simulation of adsorption (Ghiringhelli et al. Reference Ghiringhelli, Hess, van der Vegt and Delle Site2008): an enormous task in most cases. On the other hand, it has been reported in several studies that the reconstruction of some surfaces in certain conditions is negligible (Feng et al. Reference Feng, Pandey, Berry, Farmer, Naik and Heinz2011; Iori et al. Reference Iori, Corni and Di Felice2008; Raffaini & Ganazzoli, Reference Raffaini and Ganazzoli2012; Wright et al. Reference Wright, Rodger, Corni and Walsh2013a). However, other experimental and computational studies show that large scale surface reconstruction may take place after adsorption of small molecules (Eralp et al. Reference Eralp, Shavorskiy and Held2011; Gibbs et al. Reference Gibbs, Ocko, Zehner and Mochrie1990; Lal et al. Reference Lal, Plummer, Richmond and Smith2004, Reference Lal, Plummer and Smith2006). Therefore, care must be taken in modeling if a major reconstruction of the surface takes place upon adsorption.

3.4 Topography

The topographic characteristics of a surface or nanoparticle at the micro- to nanometer-scale are important determinants of protein adsorption properties, such as binding affinities and surface saturation values (Fenoglio et al. Reference Fenoglio, Fubini, Ghibaudi and Turci2011; Gagner et al. Reference Gagner, Lopez, Dordick and Siegel2011, Reference Gagner, Shrivastava, Qian, Dordick and Siegel2012; Roach et al. Reference Roach, Farrar and Perry2006). The topography of a surface can be characterized by its exposed crystal planes, its roughness and its defects (due to locally varying chemical composition or the crystalline structure of the surface), as well as kinks, edges and steps that occur during the growth of a crystal. Studies have shown that the same type of protein may have different adsorption energies, varying from highly favorable to highly unfavorable, for surfaces with different lattice structures but the same chemical composition (Heinz et al. Reference Heinz, Farmer, Pandey, Slocik, Patnaik, Pachter and Naik2009; Oren et al. Reference Oren, Tamerler and Sarikaya2005). Although modeling of the intrinsic crystal structure is straightforward, it is important to consider that different crystal planes of a material (which can be kinetically or thermodynamically favored) may be exposed during a surface adsorption process. For example, during the growth of a nanoparticle, different surface planes of the same particle can display different structures at the same time, e.g. (100) and (111) (Korzeniewski et al. Reference Korzeniewski, Climent, Feliu, Bard and Zoski2011).

Similarly, protein binding properties depend on the material structure. Rechendorff et al. (Reference Rechendorff, Hovgaard, Foss, Zhdanov and Besenbacher2006) showed that adsorption of fibrinogen on a tantalum surface can be induced by up to 70% by increasing the surface (root mean squared) roughness from 2·0 to 32·9 nm. The increase was much greater than the increase in the surface area due to the surface roughness of around 20%. On the other hand, Rechendorff et al. (Reference Rechendorff, Hovgaard, Foss, Zhdanov and Besenbacher2006) found that the adsorption of the more globular protein, bovine serum albumin, to tantalum induced by the surface roughness was similar to the increase in the surface area, thus demonstrating a selective effect of the material structure on the adsorption processes of different proteins. Finally, using molecular dynamics (MD) simulations, Nada (Reference Nada2014) investigated the interactions of an aspartic acid with step edges and kinks, as well as flat regions of a calcite crystal surface. They showed that aspartic acid binds preferentially to an acute step edge and not to an obtuse one, due to the ordered structure of water near the obtuse step edge that prevented the aspartic acid from binding strongly (Nada, Reference Nada2014). In summary, these studies demonstrate that the topography of a surface may determine its binding characteristics to proteins and peptides. Although, many topographic features can be neglected in modeling studies when comparing to experiments with well-characterized surfaces, the characteristics of the surface topography should be taken into account to realistically model and simulate the interactions of proteins with the materials used in real-life applications.

3.5 Morphology

The size of a nanoparticle, and therefore the curvature of its surface, has a strong impact on both the physico-chemical characteristics of the nanoparticle itself, such as surface polarizability, surface charge (Chiu et al. Reference Chiu, Moore, Shinoda and Nielsen2009) and isoelectric point (Suttiponparnit et al. Reference Suttiponparnit, Jiang, Sahu, Suvachittanont, Charinpanitkul and Biswas2010), and the properties of the layer of coating molecules which can determine its acidity (Wang et al. Reference Wang, Nap, Lagzi, Kowalczyk, Han, Grzybowski and Szleifer2011). The curvature of a nanoparticle may change the properties of its hydration shell as well. The solvation free energy of a nanoparticle becomes more favorable as the surface curvature decreases (and the size of the particle increases) owing to the increased water–particle interaction energy, which in turn determines its polarity (Chiu et al. Reference Chiu, Moore, Shinoda and Nielsen2009).

The effect of surface curvature on protein adsorption properties, i.e. the kinetics, thermodynamics and structural stability of adsorbed proteins, increases with decreasing particle size (Lacerda et al. Reference Lacerda, Park, Meuse, Pristinski, Becker, Karim and Douglas2010), but also depends on the adsorbate. Moreover, a strong difference in binding to spherical nanoparticles and nanorods is often observed (Gagner et al. Reference Gagner, Lopez, Dordick and Siegel2011). A structural adaption to the curvature of the nanoparticle surface upon adsorption may lead to a loss of enzymatic activity of some proteins (Wu & Narsimhan, Reference Wu and Narsimhan2008), or it may lead to significant changes in secondary or tertiary structures of self-assembling peptides (Shaw et al. Reference Shaw, Middleton, Volk and Lévy2012) or proteins (Tavanti et al. Reference Tavanti, Pedone and Menziani2015; Yang et al. Reference Yang, Johnson, Wu, Woods, George and Murphy2013) adsorbing onto the surfaces. It has been shown that while surface curvature may help to retain the tertiary structure of some proteins with globular structures adsorbed on small nanoparticles (Vertegel et al. Reference Vila Verde, Acres and Maranas2004; Lundqvist et al. Reference Lundqvist, Sethson and Jonsson2004), it can also cause significant loss in the secondary structure of a protein upon adsorption (Gagner et al. Reference Gagner, Lopez, Dordick and Siegel2011).

The effects of the surface curvature of nanoparticles can be neglected in computational studies when the binding of relatively small proteins to large nanoparticles is studied. However, the size and curvature of a nanoparticle may play a major role in the adsorption patterns of proteins not only due to geometric adaptation of the protein to a nanoparticle of similar size, but also due to changes in the physico-chemical properties of the nanoparticle itself. DFT calculations of amino acid adsorption on a single-walled carbon nanotube (SWCNT) showed that glycine adsorbs more strongly on a nanotube (3,3) than on a flat graphite surface, whereas phenylalanine adsorbs more strongly on a flat graphite surface, and the amino acids cysteine and histidine, showed no significant change in their adsorption energies (Roman et al. Reference Roman, Diño, Nakanishi and Kasai2006). Binding free energy calculations of amyloidogenic apoC-II(60–70) peptide on fullerene, CNT and graphene using MD simulations showed that the binding affinity was weakest for the fullerene and strongest for the graphene due to reduced efficiency of π-stacking interactions between the aromatic side chains of the peptide and the fullerene and CNT arising from the increased surface curvature (Todorova et al. Reference Todorova, Makarucha, Hine, Mostofi and Yarovsky2013). Xie et al. (Reference Xie, Luo, Lin, Xi, Yang and Wei2014) calculated the adsorption free energies of Alzheimer's β-amyloid peptide fragments (Aβ) to two different fullerene nanoparticle systems, C₁₈₀ and three C₆₀ (3C₆₀). They found tighter binding of the peptides on the larger C₁₈₀.

Raffaini & Ganazzoli (Reference Raffaini and Ganazzoli2013) showed that the binding strength of a protein to a carbon nanotube depends on the morphology and that the interaction energy between the protein and the nanoparticle was larger for the concave interior surface of the nanotube than for the convex outer surface. In a similar study, Chen et al. (Reference Chen, Wang, Liu, Wu, Kang, Moore and Gubbins2009b) showed that the length and the diameter of CNTs affect the energetics of the interactions with a peptide drug. While longer CNTs provide more space to trap the peptide inside the tube, a smaller diameter increases the interaction energy.

The surface morphology not only alters the binding characteristics of a peptide or protein to a surface but also the interactions between biomolecules near a surface or a nanoparticle, thus leading to intermolecular structural changes. Li et al. (Reference Li, Luo, Derreumaux and Wei2011) performed 3 separate sets of simulations of Aβ(16–22) peptides that abnormally self-assemble into β-rich aggregates: amorphous peptide in solution, amorphous peptide with SWCNT, and prefibrillar peptide with SWCNT. They observed that without the CNTs, the amorphous peptides form β-sheet structures. On the other hand, simulations with SWCNTs showed that the amorphous peptides tend to form disordered coils, whereas the β-sheet structures formed by the prefibrillar peptides were destabilized due to the interactions of the peptides with the CNTs.

Finally, noteworthy to mention is that a change in the morphology of a surface/particle will lead to changes in its physico-chemical properties. In a study by Emami et al. (Reference Emami, Puddu, Berry, Varshney, Patwardhan, Perry and Heinz2014b), it was shown that the size of the silica nanoparticles determines their surface ionization levels. The larger nanoparticles have higher ionization and therefore bind more strongly to peptides with high net positive charge. Further, Baier et al. (Reference Baier, Blumenstein, Preusker, Jeurgens, Welzel, Do, Pleiss and Bill2014) recently investigated the binding free energy profile of 12-mer peptides on polar (001) and non-polar (100) ZnO surfaces. Employing an enhanced sampling approach (see Section 10.2), the authors showed that there is a positive correlation between the ZnO particle sizes obtained in experiments in the presence of the peptides and the calculated affinities of the peptides for the ZnO surface. Their result showed that the selective adsorption of a peptide can impact the growth of certain nanocrystals. Therefore, one should consider modeling, not only the size and curvature of the material, but also other properties that depend on the morphology, in particular for nanoparticles.

6.1 Interaction potentials

The classical potential energy functions employed in the all-atom molecular mechanics force fields (MM FFs) for biomolecules, such as AMBER (Cornell et al. Reference Cornell, Cieplak, Bayly, Gould, Merz, Ferguson, Spellmeyer, Fox, Caldwell and Kollman1995), CHARMM (Brooks et al. Reference Brooks, Bruccoleri, Olafson, States, Swaminathan and Karplus1983), GROMOS (Oostenbrink et al. Reference Oostenbrink, Villa, Mark and Van Gunsteren2004), and OPLS-AA (Jorgensen et al. Reference Jorgensen, Maxwell and Tirado-Rives1996), are widely used and thoroughly evaluated for simulation of biomolecules in aqueous solution. Most commonly used biomolecular FFs are expressed as a sum of pairwise interaction terms that represent changes of the (i) chemical bond lengths and bond angles of a molecule as harmonic spring functions; and (ii) torsions as periodic functions (dihedral angles, or torsional rotation of atoms around a central bond); and (iii) non-bonded electrostatic and van der Waals intra-  and inter-molecular interactions:

$${E_{{\rm Total}}} = {E_{{\rm bond}}} + {E_{{\rm bond}}}_{\_{\rm angle}} + {E_{{\rm torsion}}} + {E_{{\rm electrostatic}}} + {E_{{\rm vdW}}}$$

Electrostatic interactions between the atoms in the system are approximated by Coulomb's law with fixed point charges assigned to each atom, whereas van der Waals interactions are typically described by the Lennard–Jones (12–6) potential:

$${U_{{\rm LJ}}} = 4{\varepsilon _{ij}}\left[ {{{\left( {\displaystyle{{{\sigma _{ij}}} \over r}} \right)}^{12}} - {{\left( {\displaystyle{{{\sigma _{ij}}} \over r}} \right)}^6}} \right]$$

where the two terms represent repulsive and attractive interactions, respectively, and the parameters, and σ are expressed as a combination of parameters of atoms i and j. If this form of the FF is suitable for the solid surface to be studied, it can also be applied for simulations of protein adsorption.

The FFs derived specifically for simulation of inorganic materials are generally aimed at reproducing the structural and physico-chemical properties of the bulk or surface of a material in vacuum and often have a functional form that is different from that of conventional biomolecular force-fields. Usually, monatomic systems with close-packed structures can be reasonably well described by two-body interaction potentials, whereas a three-body potential must be used for semiconductors (Vashishta et al. Reference Vashishta, Kalia, Rino and Ebbsjö1990). Additional angle-dependent terms can be employed to preserve the directionality of bonds, e.g. in the crystal packing arrangement of a bulk solid (Cruz-Chu et al. Reference Cruz-Chu, Aksimentiev and Schulten2006). Furthermore, a Buckingham potential of the form Ae ^−Br −(C/r ⁶), where A, B and C are force field parameters, is often used instead of a Lennard–Jones potential to better describe the repulsion between ions in inorganic materials (see, for example, Van Beest & Kramer, Reference Van Beest, Kramer and van Santen1990). Moreover, a balanced parameterization of a simple pairwise interaction energy model commonly used in standard biomolecular FFs is not straightforward and is not always possible for the water-material interface. For instance, a combination of two-body and three-body non-bonded potentials was employed in a silica force-field to preserve the tetrahedral structure of the silica glass and to reproduce accurately silica surfaces, pores and surface wettability (Cruz-Chu et al. Reference Cruz-Chu, Aksimentiev and Schulten2006). Another limitation of a simple pair-wise function is that it is based on a fixed point charge approximation for Coulomb's interaction, which omits the interactions between higher order multipoles and polarization effects arising from the electrostatic field of the environment (a non-polarizable force-field). The partial atomic charges are usually assigned for a molecule in an aqueous environment and, thus, the polarization effect of the aqueous surroundings is implicitly accounted for. Non-polarizable FFs have been extensively validated over the last decades and currently offer a robust description of the equilibrium properties of biomolecules in aqueous solution. However, polarization effects of the surface (such as metal or carbon surfaces, as discussed earlier) may contribute notably to the binding energy of an adsorbate, as well as to the structure of the surface itself.

Polarizable FFs for biological molecules are currently under active development and validation (Halgren & Damm, Reference Halgren and Damm2001). However, apart from their notably higher computational cost, polarizable FFs for biomolecules are not designed to reproduce the polarization properties of inorganic surfaces. In Section 7, we will consider several models adapted for simulation of surface polarization effects.

Despite all the limitations mentioned above, the early simulations of protein and peptides on solid surfaces have typically been performed using standard biomolecular FFs. In the last decade, however, much effort has been invested to evaluate and adapt biomolecular FFs to the simulation of interfaces between biomolecules in water and specific inorganic materials, and several new approaches for developing general bio-inorganic FFs have been reported. For example, a dual-scale FF that includes two sets of partial atomic charges, optimized for each phase separately, has been developed (Biswas et al. Reference Biswas, Vellore, Yancey, Kucukkal, Collier, Brooks, Stuart and Latour2012). In the recently reported INTERFACE FF (Heinz et al. Reference Heinz, Lin, Kishore Mishra and Emami2013), partial atomic charges and van der Waals parameters have been optimized using a non-polarizable fixed-charge model for a large number of different materials based on the properties of the solid surfaces, e.g. crystal unit parameters, surface tension, water contact angle and interfacial tension, and adsorption energy of selected peptides. Some applications of these new models for specific surfaces will be discussed in Section 7.

6.2 Solvation models

The properties of the hydration shell of a solid substrate, in particular variations of physical properties, such as the density, free energy and dielectric constant of water in a hydration shell, and chemical properties, such as surface reconstruction and ionization, often govern adsorption processes. Therefore, a rigorous treatment of water, whether explicit or implicit, is required to account for the effects of water on adsorption. FF parameters for water molecules are usually parameterized for the bulk water state and thus do not take into account the specific features of the surface or protein hydration shell. Moreover, the properties of a solid surface hydration shell differ from those of small solutes and biomolecules in that, generally, more water–water hydrogen bonds are broken and the water is forced to build an extended ordered surface layer. In general, there is still uncertainty as to whether standard explicit water models are able to adequately reproduce adsorption phenomena for proteins and surfaces. Studies of water properties on a wide range of surfaces have been reviewed by Henderson (Reference Henderson2002).

The microscopic properties of the water molecules on a surface can be parameterized using DFT ab initio computations, which have been widely used for predicting the adsorption energies of water molecules as well as water dissociation properties on different surfaces. The AIMD technique mentioned earlier is a powerful tool for exploring water properties in a monolayer or even in several layers (e.g. on gold, see (Velasco-Velez et al. Reference Velasco-Velez, Wu, Pascal, Wan, Guo, Prendergast and Salmeron2014)). However, it was shown that the spatial and temporal limitations of DFT simulations may prevent accurate modeling of water structural transitions from the hydration shell to bulk, which is important for the modeling of protein adsorption (Große Holthaus et al. Reference Große Holthaus, Köppen, Frauenheim and Colombi Ciacchi2012). Accordingly, an all-atom MM representation of water has to be used for exploring the adsorption processes of molecules in solution, which usually occur on a much longer time-scale. Likewise, an all-atom MM model of the solvent can be replaced by a continuum solvent model to enable longer simulations of large proteins or protein solutions for investigating their adsorption properties.

Implicit solvation is a common way to reduce the computational cost of simulations that arises in large part from calculating the interactions between the explicit water molecules. In this approach, the contribution of the solvent to the binding energies between the solutes is represented with a mean field model, and is usually decomposed into two major contributions: electrostatic contributions due to charge interactions and the change of dielectric polarization of their surroundings, and non-polar contributions due to the change in entropy and the dispersion energy upon removing water from the solute surface. The non-polar desolvation energy is commonly estimated as a linear function of the solvent accessible surface area. The most common approaches for computing the continuum electrostatic contribution to the binding free energy is to solve the Poisson–Boltzmann (PB) equation or the Generalized Born (GB) model, which provides an approximation to the PB method. The PB method is often used with a rigid molecule approximation because it is relatively computationally expensive to solve the PB equation for every macromolecular configuration during a simulation. For applications of these methods to free energy calculations of protein-surface systems, see Section 10.

The implicit solvent models traditionally used for molecules in solution do not, per se, account for the unique properties of the hydration shell and the interactions of the water in direct contact with the inorganic surface. However, they can be parameterized to account for the hydration shell characteristics of individual materials. An example of such a continuum solvent model is ProMetCS, which was derived for protein-gold (111) interfaces (Kokh et al. Reference Kokh, Corni, Winn, Hoefling, Gottschalk and Wade2010). In ProMetCS, the effects of the hydration shell distortion at the protein-surface interface were accounted for by the addition of an analytical function parameterized to reproduce the potential of mean force (PMF) of a probe ion/atom on the gold surface obtained from explicit water MD simulations. The effects arising from the partial replacement of the metal hydration shell by a protein adsorption site are described by a free energy penalty, proportional to the protein-surface contact area, compensating the Lennard–Jones attraction to a large extent. The PMF of an atom on the gold surface computed using the ProMetCS FF reproduces the profile of the PMF function obtained in explicit solvent MD simulations (Fig. 7).

Fig. 7. Modeling of protein-Au(111) interactions using a continuum solvent model parameterized by comparison with explicit solvent MD simulations. Solid line: The potential of mean force (PMF) for a test atom as a function of atom-gold surface distance, as obtained from MD simulations in explicit water solvent; squares: the corresponding LJ potential; dashed line: their difference, associated with the desolvation energy; dotted line: PMF computed using the protein-metal continuum solvent ProMetCS model (which includes both LJ and metal hydrophobic desolvation energies, see text). Reprinted with permission from (Kokh et al. Reference Kokh, Corni, Winn, Hoefling, Gottschalk and Wade2010). Copyright (2010) American Chemical Society.

7.1 Metal surfaces

The metal polarization in the presence of the charges of a solute, as discussed earlier, has a significant effect on adsorption and, therefore, efforts have been put into integrating polarization effects into FF models. The simplest way to introduce the induced polarization of the metal surface in a computational model is using a classical image-charge approximation for a charge and a zero-potential metal surface in a continuum dielectric medium with an interaction energy expressed by Coulomb's law between the charge and its image of the opposite sign (see Fig. 8a ). This approximation was shown to give an accurate description for the interaction energy of an Al(111) surface with a charge when the separation distance was more than 2·5 Å (Finnis & Finnis, Reference Finnis1991), while at closer distances to the surface, the image charge model overestimated the interaction potential (Finnis & Finnis, Reference Finnis1991; Smith et al. Reference Smith, Chen and Weinert1989). The choice of the image plane position was discussed by Heinz et al. (Reference Heinz, Jha, Luettmer-Strathmann, Farmer and Naik 2011 ) who suggested that the image plane may deviate from the jellium edge (which corresponds to a plane half a lattice spacing above the first layer of surface atoms). The position of this image plane is further away from the surface atoms than the jellium edge and its deviation is represented by a system dependent offset, δ, which is determined by the response of the electron density of the metal (Smith et al. Reference Smith, Chen and Weinert1989). Heinz et al. (Reference Heinz, Jha, Luettmer-Strathmann, Farmer and Naik2011) exploited this scheme to estimate polarization effects a posteriori from non-polarizable fully atomistic MD simulations of a water-Au surface, as well as a water-peptide-gold system.

Fig. 8. Models of polarizable gold surfaces. (a) The method of image charges. Charges q _i and q _j induce polarization charges inside the metal shown as −q _i and −q _j respectively. (b) Classical polarization models. Each gold atom of the surface is assigned a dipole with either variable (Drude oscillator) (left) or fixed (rigid rod that is free to rotate) (right) moment. Reproduced with permission from (Iori & Corni, Reference Iori and Corni2008). Copyright (2008) John Wiley and Sons.

Although implementation of the image-charge approximation in an all-atom MD simulation is straightforward, it is impractical for large systems due to the increased computational load that scales with the number of interacting particles. In an alternative approach, which can be incorporated into any commonly used MD energy function, virtual dipoles or rods that can adjust their position in response to the external electrostatic field are introduced on all the surface atoms. Each virtual dipole is constrained at one end to a real surface atom and, depending on the model, can either change the magnitude of its dipole moment (Drude oscillator) or its orientation (rigid rods, see Fig. 8b ) and, hereby screen the external electrostatic field. A model with virtual rigid rods for simulation of the polarization of a metal surface was proposed and implemented by Iori & Corni (Reference Iori and Corni2008) and used in the GolP family of FFs optimized for Au surfaces (Iori et al. Reference Iori, Di Felice, Molinari and Corni2009; Wright et al. Reference Wright, Rodger, Corni and Walsh2013a, Reference Wright, Rodger, Walsh and Cornib), as well as for Ag and graphene.

The first optimized sets of Lennard–Jones parameters to include metal–water and metal–amino acid interactions was developed by Ghiringhelli et al. (Reference Ghiringhelli, Hess, van der Vegt and Delle Site2008) on the basis of DFT calculations, and by Heinz et al. (Reference Heinz, Vaia, Farmer and Naik2008) and Vila Verde et al. (Reference Vertegel, Siegel and Dordick2009) to reproduce Au(111) hydrophilicity. In the GolP FF (Iori et al. Reference Iori, Di Felice, Molinari and Corni2009), which is based on the OPLS FF, additional parameters describing the interactions of biomolecular groups with the gold surface are parameterized from QM calculations and experimental studies of adsorption energy. To reproduce the binding energy and orientation of small molecular fragments, the authors included and optimized a set of additional parameters that describe the van der Waals interactions with gold as well as stronger, chemical-like bonds between aromatic groups and gold atoms in the form of a Lennard–Jones function. This parameterization can be directly used in standard MM force-fields. The GolP FF has been applied in MD simulations of the adsorption of amino acids, as well as several proteins, on gold (Brancolini et al. Reference Brancolini, Kokh, Calzolai, Wade and Corni2012; Cohavi et al. Reference Cole, Payne, Csányi, Mark Spearing and Colombi Ciacchi2011; Hoefling et al. Reference Hoefling, Iori, Corni and Gottschalk2010a, Reference Hoefling, Iori, Corni and Gottschalkb, Reference Hoefling, Monti, Corni and Gottschalk2011; Kokh et al. Reference Kokh, Corni, Winn, Hoefling, Gottschalk and Wade2010). Furthermore, the ProMetCS (Kokh et al. Reference Kokh, Corni, Winn, Hoefling, Gottschalk and Wade2010) continuum solvent model describes protein–surface interactions in atomic detail using the GolP Lennard–Jones interaction parameters together with an image charge model for protein-metal surface electrostatic interactions.

Following the same strategy, a new force-field, GolP-CHARMM (Wright et al. Reference Wright, Rodger, Corni and Walsh2013a, Reference Wright, Rodger, Walsh and Corni2013b), has recently been developed on the basis of the CHARMM FF for studying the interaction of proteins with Au(111) and Au(100) surfaces. In GolP-CHARMM, special attention has been paid to the correct description of the properties of the water molecules adsorbed on gold, with the properties being compared with ab initio MD simulations (Cicero et al. Reference Clementi, Nymeyer and Onuchic2011). The new FF improves the quality of the simulation of the hydration shell, reproducing the increased tendency of interfacial water molecules to donate hydrogen bonds to other water molecules compared with those in the bulk water and reproducing the energetics of water molecules on a gold surface.

Heinz et al. (Reference Heinz, Vaia, Farmer and Naik2008) developed Lennard–Jones parameters for face-centered cubic metals based on experimentally determined densities and surface tensions at 298 K under atmospheric pressure for use in modeling of mixed interfaces. They provide two sets of parameters for metal atoms (Ag, Al, Au, Cu, Ni, Pb, Pd, Pt); one set suited for 12–6 Lennard–Jones functions to be used in AMBER, CHARMM, CVFF and OPLS-AA force-fields, and another for 9–6 Lennard–Jones functions to be used in the COMPASS (Sun, Reference Sun1998) and PCFF force-fields. Penna et al. (Reference Penna, Mijajlovic and Biggs2014) used surface models and the Lennard–Jones parameters of Heinz et al. (Reference Heinz, Vaia, Farmer and Naik2008) with a CHARMM FF (MacKerell et al. Reference MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, Joseph-McCarthy, Kuchnir, Kuczera, Lau, Mattos, Michnick, Ngo, Nguyen, Prodhom, Reiher, Roux, Schlenkrich, Smith, Stote, Straub, Watanabe, Wiórkiewicz-Kuczera, Yin and Karplus1998) to investigate the interactions of peptides with Au and Pt surfaces.

In addition to polarization, surface charges should be taken into account when modeling a surface. Bizzarri (Reference Bizzarri2006) investigated the dynamics and electron transfer (ET) properties of the bacterial ET protein, azurin, anchored to neutral, positively and negatively charged gold surfaces using classical MD simulations. They calculated the ET rate between the azurin copper atom and the sulfur atom of the cysteine surface anchor using classical Marcus theory (Marcus & Sutin, Reference Marcus and Sutin1985). The computed ET rate was highest on the positively charged surface and lowest on the neutral surface. However, it was not clear to what extent the relatively larger structural changes observed in the molecule anchored on the positively charged surface play a role in this trend.

7.2 Titanium oxide surfaces

Because of its high stability and resistance to corrosion, titanium is the material of choice for many medical and technical applications and it is, therefore, one of the most widely used materials in computational modeling. Crystals of TiO₂ (titania) appear in three forms: rutile, anatase and brookite. The rutile form, in particular its interaction with water, has been the most studied by computer modeling. An oxide layer, predominantly of amorphous structure, is formed at the titanium-water interface (Schneider & Colombi Ciacchi, Reference Schneider and Colombi Ciacchi2010). The properties of this layer depend strongly on the physico-chemical conditions and the temperature regime of adsorption.

There is controversy between theory and experiment regarding the extent to which water dissociates on TiO_x; simulations often show a larger degree of dissociation than experiment (Diebold, Reference Diebold2003; Huang et al. Reference Huang, Gubbins, Li and Lu2014; Sun et al. Reference Sun, Liu, Selloni, Lu and Smith2010). It is believed that for rutile (110), dissociation takes place predominantly on defects, while water can also dissociate to some extent on a defect-free (100) surface, and effectively dissociates on a (001) surface. Accordingly, the density distribution of water strongly varies for different TiO₂ surfaces from two sharp maxima on a rutile surface to a smooth distribution on hydrophobic titania (Huang et al. Reference Huang, Gubbins, Li and Lu2014). The mixture of hydroxyl groups and molecularly adsorbed water may form a partially ordered solvation layer that influences interactions of peptides with the surface (Li et al. Reference Li, Monti and Carravetta2012).

The behavior of water on different types of TiO₂ surfaces, including partial dissociation of water, was investigated using the reactive FF ReaxFF (Huang et al. Reference Huang, Gubbins, Li and Lu2014; Kim et al. Reference Kim, Kumar, Persson, Sofo, van Duin and Kubicki2013) and AIMD (Carravetta et al. Reference Carravetta, Monti and Zhang2009; Schneider & Colombi Ciacchi, Reference Schneider and Colombi Ciacchi2011). More details on theoretical studies of titania–water interactions can be found in the review by Sun et al. (Reference Sun, Liu, Selloni, Lu and Smith2010).

As mentioned before, water dissociation on oxide surfaces is reversible and cannot be simply described by classical MM force-fields. Bandura & Kubicki (Reference Bandura and Kubicki2003) pointed out that the two-body interaction potential is not suitable for capturing the correct directionality of hydroxyl groups on a TiO₂ surface and that this problem can be solved by using a polarizable FF or by introducing bond stretching and angle bending terms for Ti-O-H. They used the Matsui & Akaogi (Reference Matsui and Akaogi1991) model which incorporated additional bond-stretching and angle-bending terms for O-Ti-H into the Buckingham potential for the atoms of the (110) surface of rutile. The FF was later refined by Předota et al. (Reference Předota, Bandura, Cummings, Kubicki, Wesolowski, Chialvo and Machesky2004) and then revised for anatase (101), (112) and rutile (110) surfaces. In the FF developed by Borodin et al. (Reference Borodin, Smith, Bandyopadhyaya and Byutner2003), an exponential term for modeling short-range repulsion was added to the standard Lennard–Jones function together with an additional polarization term, and the FF was optimized against the results of QM calculations for polyethylenoxide on small TiO₅H₉ clusters and then applied to TiO₂ surfaces.

Since the use of AIMD and reactive FFs is too time-consuming for simulations of large systems and/or long time-scales, TiO_x parameters have been incorporated into standard MM FFs employed in biomolecular simulations. A classical FF with a Finnis–Sinclairtype many-body potential for the surface coupled to a Buckingham potential for the short-range repulsive interactions was developed by Schneider and Colombi Ciacchi and applied to simulating peptide interactions on a Ti/TiO_x interface (Schneider & Colombi Ciacchi, Reference Schneider and Colombi Ciacchi2010). The parameters of the AMBER FF were optimized for TiO₂ interactions with collagen peptides in a study by Köppen & Langel (Reference Köppen and Langel2010). The ratio of protonated, hydroxylated and stoichiometric units on the surface was adapted to the physiological pH value and retained during MD simulations. The charges in the bulk of TiO₂ were fixed, but different charges on the surface oxygen and hydrogen atoms were tested. The binding of Glu and Lys to TiO₂ was found to be mediated by hydroxyl groups on the surface, with the adsorption energy strongly dependent on the charge model used.

Friedrichs & Langel (Reference Friedrichs and Langel2014) re-parameterized the Matsui–Akaogi potential to the standard Lennard–Jones form for a description of the interactions between the rutile surface atoms and used their model to simulate adsorption of peptides on the rutile surface. Although this enabled standard MD software to be used, the simplicity of the form of the potential gave rise to some deviation of physical properties, such as crystal structure parameters and permittivity, from those observed in experiments and simulations with the original potential.

The adsorption of different peptides onto a rutile surface, TiO₂(110), was simulated in a series of studies by Monti and colleagues of the Ala-Glu and Ala-Lys dipeptides (Carravetta & Monti, Reference Carravetta and Monti2006) and (Monti et al. Reference Monti, Carravetta, Battocchio, Iucci and Polzonetti2008), and β-sheet conformations of the EAK16 and RAD16 oligopeptides (Carravetta et al. Reference Carravetta, Monti and Zhang2009). The FF for biomolecule-surface and water–surface interactions used in these studies was optimized against data from QM calculations for a small titania cluster and experimental data (Carravetta & Monti, Reference Carravetta and Monti2006). The TIP3P water model and a combination of a Buckingham potential for the interaction between the surface and surface-water atoms and a Lennard–Jones potential for most of the rest of the atoms were employed. Peptide binding to the surface was found to preferentially occur through the backbone atoms; particularly carboxyl oxygens and the C-terminal carboxy group on the peptides interacted with exposed Ti atoms. Water molecules were also observed to mediate peptide-surface binding through hydrogen bonds. Since the dissociation of water was not taken into account, the contribution of hydroxyl groups to peptide binding was not analyzed. Finally, the adsorption of glycine (Li et al. Reference Li, Monti and Carravetta2012), and cysteine (Li et al. Reference Li, Monti, Ågren and Carravetta2014) onto TiO₂ was explored using the reactive force-field, Reaxff (Van Duin et al. Reference Van Duin, Dasgupta, Lorant and Goddard2001), extended to treat the interaction of solid inorganic substrates with biomolecules. It is interesting that these adsorbing molecules, especially glycine, tend to form self-assembled clusters and that only weak adsorption was observed. Additionally, the proton transfer reactivity of the amino acid is observed to be enhanced by the presence of the surface.

It should be noted that MM simulations of proteins on reactive surfaces like titanium oxide should be performed with caution, because the presence of defects may change the surface properties dramatically, and thus affect the adsorption behavior. In a study by Utesch et al. (Reference Utesch, Daminelli and Mroginski2011), the BMP-2 protein interacted unexpectedly weakly with the TiO₂ surface, because it was unable to penetrate the highly-ordered water layer formed on the defectless surface.

7.3 Silicon oxide surfaces

A number of FFs have been developed for MD simulation of bulk silicon oxides and their interfaces with water due to the wide range of technological applications of silica-based materials (see reviews by Butenuth et al. (Reference Butenuth, Moras, Schneider, Koleini, Köppen, Meißner, Wright, Walsh and Colombi Ciacchi2012) and by Rimola et al. (Reference Rimola, Costa, Sodupe, Lambert and Ugliengo2013)). There are, however, some challenges in building a FF for silicon oxide surfaces. In particular, for modeling the tetragonal bulk structure of SiO₂ (silica), many-body potentials are often used, which are generally not compatible with common biomolecular force-fields. At the SiO_x-water interface, water molecules dissociate and saturate dangling Si and O groups, which leads to formation of different types of exposed functional groups: hydrophilic silanols (Si-OH), and hydrophobic siloxanes (with oxygen buried, Si-O-Si) and silanes (with silica atoms saturated by hydrogen). Silanol groups can make hydrogen bonds to each other or two silanol hydroxyl groups can bind to one Si atom. Since they are hydrophilic in character, these groups form hydrogen bonds with surface water molecules. Crystalline silica has a dense ordered network of geminal hydrogen-bonded silanol groups leading to ordered water layers (Notman & Walsh, Reference Notman and Walsh2009) (see Fig. 9), while amorphous silica has isolated geminal and hydrogen-bonded silanol groups resulting in only locally ordered regions of water molecules and hydroxyl groups (Aarts et al. Reference Aarts, Pipino, Hoefnagels, Kessels and van de Sanden2005).

Fig. 9. Lateral density profiles of the first layer of structured water on the (101-0), (0001) and (011-1) surfaces of silica. Reprinted with permission from (Notman & Walsh, Reference Notman and Walsh2009). Copyright (2009) American Chemical Society. Positions of the water molecules are dominated by the positions of the silanol hydroxyl groups and to a lesser extent by the underlying crystal structure.

Silica has an isoelectric point of about 2–3 (Parks, Reference Parks1965), and, at neutral pH, between 5% and 20% of silanol groups are deprotonated, leading to a negative surface charge. Experimental studies of peptide adsorption on silica nanoparticles in aqueous solution (Patwardhan et al. Reference Patwardhan, Emami, Berry, Jones, Naik, Deschaume, Heinz and Perry2012; Puddu & Perry, Reference Puddu and Perry2012) showed that binding was dominated by electrostatic interactions between negatively charged siloxide groups and ammonium groups at N-termini, and Lys and Arg containing peptides and, to a lesser extent, hydrogen bonds between negatively charged silanol and siloxide groups and polar groups in Ser, His and Asp-containing peptides. Binding of hydrophobic peptides was observed at low pH (Puddu & Perry, Reference Puddu and Perry2012). A FF for MD simulations of a deprotonated negatively charged SiO₂ surface in solution at neutral pH was reported by Butenuth et al. (Reference Butenuth, Moras, Schneider, Koleini, Köppen, Meißner, Wright, Walsh and Colombi Ciacchi2012).

Recently, several FFs for silica compatible with biomolecular FFs have been developed. A FF based on the CHARMM FF was parameterized by Lopes et al. (Reference Lopes, Murashov, Tazi, Demchuk and MacKerell2006) to obtain a structural and dynamic representation of water in the vicinity of neutral quartz crystalline surfaces using MD simulations. This FF was used for a set of peptides known to be strong and weak binders to the quartz (100) surface and showed that Pro, Trp and Leu were the main residues forming close contacts with the surface (Oren et al. Reference Oren, Notman, Kim, Evans, Walsh, Samudrala, Tamerler and Sarikaya2010).

The GLASSFF_2·01 FF is compatible with CHARMM and with TIP3P water parameters and was developed by Cruz-Chu et al. (Reference Cruz-Chu, Aksimentiev and Schulten2006) for simulation of amorphous silica surfaces and nanopores. This FF employs standard two-body Lennard–Jones and Coulomb potentials as well as a three-body directional term for modeling the tetrahedral arrangement of bulk silica. To construct amorphous silica, the authors applied annealing cycles in MD simulations, which initiated surface reconstruction from crystal to amorphous silica accompanied by formation of dangling atoms (oxygens with less than two bonds and silicons with less than four bonds) that were required for modeling the wetting properties of silica correctly.

A comparison of FFs for the prediction of the physical and chemical properties of bulk silica and its aqueous interface such as atomic charge, bond length and angles, density of silanol groups on the surface and degree of ionization of silanol groups can be found in the papers by Skelton et al. (Reference Skelton, Fenter, Kubicki, Wesolowski and Cummings2011) and by Emami et al. (Reference Emami, Puddu, Berry, Varshney, Patwardhan, Perry and Heinz2014a). In the latter study, a data set with more than 20 different silica surface models that cover most important types of silica surface chemistry at different surface ionization states was assembled. Then, atomic parameters for the different chemical groups in silica were derived to reproduce a large set of chemical and physical properties of bulk and surface silica as well as the water contact angle, heat of water immersion and water adsorption isotherms observed in experiments. The parameters were integrated into several FFs (AMBER, CHARMM, COMPASS, INTERFACE (Heinz et al. Reference Heinz, Lin, Kishore Mishra and Emami2013)) compatible with the TIP3P and SPC water models. In particular, the CHARMM-INTERFACE FF (Heinz et al. Reference Heinz, Lin, Kishore Mishra and Emami2013) was applied to simulation of peptide adsorption on silica nanoparticles of different size at different pH values (Emami et al. Reference Emami, Puddu, Berry, Varshney, Patwardhan, Perry and Heinz2014b), showing good agreement of nanoparticle coverage with experimental data (Puddu & Perry, Reference Puddu and Perry2012).

The Dual-FF  (Biswas et al. Reference Biswas, Vellore, Yancey, Kucukkal, Collier, Brooks, Stuart and Latour2012), (discussed in more detail in Sec.7·5) was optimized to reproduce the experimental binding energies of TGTG-X-GTGT peptides, with X = N, D, G, K, F,T, W and V, on hydroxylated quartz (100) and glass surfaces in MD simulations (Snyder et al. Reference Snyder, Abramyan, Yancey, Thyparambil, Wei, Stuart and Latour2012). In the Dual-FF, CHARMM parameters were used for biomolecules, while interfacial parameters, such as partial charges and Lennard–Jones parameters for silica–peptide and silica–water interactions, were tuned to obtain the best agreement with peptide adsorption energies.

7.4 Mineral surfaces

The main difficulty in simulating mineral-biomolecule interfaces in aqueous solution is that, while both minerals and biomolecules may have atom types in common, e.g. oxygen, the FF for bonded terms are different in minerals and biomolecules. Specifically, minerals are ionic systems and are best represented by ionic bonds based on Coulomb electrostatics and Pauli repulsion terms (Hauptmann et al. Reference Hauptmann, Dufner, Brickmann, Kast and Berry2003). In contrast, biomolecules are covalently bonded. As for oxide surfaces, water may dissociate on a mineral surface leaving adsorbed hydroxyl groups at charged group sites on the surface with a partial occupancy of about 50% of the charged sites (De Leeuw, Reference De Leeuw2010). This process cannot be modeled in a standard MM FF, but the presence of hydroxyl groups/hydrogens on the surface may strongly affect the adsorption properties of biomolecules.

A number of FFs for modeling the bulk and surface properties of minerals have been developed (see (Mishra et al. Reference Mishra, Flatt and Heinz2013) and references therein), and several non-reactive FFs for the simulation of mineral-water interfaces have been reported. A general FF, CLAYFF, for modeling the interface of a hydrated multicomponent mineral with aqueous solution was proposed by Cygan et al. (Reference Cygan, Liang and Kalinichev2004), where charges on oxygen atoms and hydroxyl groups vary according to their structural environment. This FF includes a flexible SPC-based model for water and hydroxyl groups and treats most interatomic interactions as non-bonded with Lennard–Jones and Coulomb terms.

A methodology for the generation of a general bio-compatible FF for minerals was proposed by Freeman et al. (Reference Freeman, Harding, Cooke, Elliott, Lardge and Duffy2007). In this work, existing FFs and models were used for the various components of the system (including the TIP3P model for water), while new parameters were derived for cross-species interaction terms, for example, for the interaction of a hydroxyl group with a Ca²⁺ ion. This methodology was used for simulation of small organic molecules on calcite and magnesite surfaces (Freeman et al. Reference Freeman, Asteriadis, Yang and Harding2009) by combining the AMBER FF for organic molecules with inter-species terms proposed by Freeman et al. (Reference Freeman, Harding, Cooke, Elliott, Lardge and Duffy2007); and for simulation of the binding of ovocleidin-17 protein to calcium carbonate nanoparticles (Freeman et al. Reference Freeman, Harding, Quigley and Rodger2010) and surfaces (Freeman et al. Reference Freeman, Harding, Quigley and Rodger2011). Furthermore, Katti et al. reported parameters for modeling clay minerals (NaSi₁₆(Al₆FeMg)O₄₀(OH)₈) in the CHARMM FF (Katti et al. Reference Katti, Schmidt, Ghosh and Katti2005b) and amino-acid adsorption to these clay minerals (Katti et al. Reference Katti, Ghosh, Schmidt and Katti2005a).

Hydroxyapatite (HAP, Ca₁₀ (PO₄)₆ (OH)₂) is the main mineral component of bones and teeth and thus a promising material for application in bone replacement. A FF for HAP has been developed with a Born–Mayer–Huggins potential with an exponential repulsion term similar to that of Buckingham potential instead of the standard Lennard–Jones plus Coulomb expression for the non-bonded interactions. The FF parameters were derived by fitting the QM electrostatic field in a 3D crystal environment and experimental crystal parameters (Hauptmann et al. Reference Hauptmann, Dufner, Brickmann, Kast and Berry2003). This FF was subsequently re-parameterized for monoclinic hydroxyapatite using an energy function consistent with common biomolecular FFs (Bhowmik et al. Reference Bhowmik, Katti and Katti2007). Dong et al. (Reference Dong, Wang, Wu and Pan2007) adapted these HAP parameters for the standard biomolecule FF for modeling the adsorption dynamics of the BMP-2 protein on a HAP (001) surface. In other studies, the Lennard–Jones and Coulomb parameters proposed by Hauptmann et al. (Reference Hauptmann, Dufner, Brickmann, Kast and Berry2003) were combined with CHARMM parameters for the protein using the standard mixing rule in order to study the adsorption of fibronectin on HAP (001) (Shen et al. Reference Shen, Wu, Wang and Pan2008) and with OPLS-AA parameters for simulating the HAP-water interface (Zahn & Hochrein, Reference Zahn and Hochrein2003).

7.5 Self-assembled monolayer surfaces

The applicability of several commonly used biomolecular FFs for the simulation of peptide adsorption on hydrophobic and hydrophilic SAMs has been systematically evaluated by Latour and colleagues for CHARMM19 (Sun & Latour, Reference Sun and Latour2006; Sun et al. Reference Sun, Dominy and Latour2007; Vellore et al. Reference Vellore, Yancey, Collier, Latour and Stuart2010), and CHARMM22, OPLS-AA, and AMBER94 (Collier et al. Reference Colombi Ciacchi and Payne2012) FFs. In the latter study by Collier et al. (Reference Colombi Ciacchi and Payne2012), the FF parameters and partial charges of the SAMs (-OH and -COOH) were assigned by analogy to amino acids with similar functional groups. The authors compared experimental findings with their computed free energies of adsorption and qualitative behavior of peptides on different SAM surfaces such as the change in conformation upon adsorption and the orientation on the surface. This study (Collier et al. Reference Colombi Ciacchi and Payne2012) demonstrated that although some FFs perform reasonably well (the best match to experiment was obtained in simulations using CHARMM22 and AMBER94 (Cornell et al. Reference Cornell, Cieplak, Bayly, Gould, Merz, Ferguson, Spellmeyer, Fox, Caldwell and Kollman1995)), none of the FFs capture the specific interaction properties of the SAM-water interface. In particular, systematic overestimation of the binding strength of hydrophobic peptides and underestimation for negatively charged peptides was observed. On the other hand, Collier et al. (Reference Colombi Ciacchi and Payne2012) also noted that altering FF parameters to reproduce the properties of adsorbed peptides unavoidably led to alteration of peptide behavior in solution. These results suggested that for accurate simulations of peptide adsorption a new FF strategy was required.

A new approach, the dual-FF, was proposed by Biswas et al. (Reference Biswas, Vellore, Yancey, Kucukkal, Collier, Brooks, Stuart and Latour2012). In this FF, different sets of non-bonded parameters, i.e. atomic partial charges and parameters of the Lennard–Jones potential, were used to represent intra-phase (i.e. between peptides or between water molecules) and inter-phase (peptide–SAM) interactions simultaneously in simulations. Furthermore, the non-polarizable TIP3P (Jorgensen et al. Reference Jorgensen, Chandrasekhar, Madura, Impey and Klein1983) water model was reparameterized for water-SAM interfaces based on the electrostatic properties of water molecules in the interface computed using the polarizable TIP4P-FQ (Rick et al. Reference Rick, Stuart and Berne1994) FF. In addition, extended experimental benchmarks of binding free energies of peptides with different sequences on various functionalized SAM surface interfaces by Wei & Latour (Reference Wei and Latour2009, Reference Wei and Latour2010) were used for optimization of the parameters for amino-acids.

Utesch et al. (Reference Utesch, Sezer, Weidinger and Mroginski2012) studied the adsorption of sulfite oxidase (SO) on a functionalized SAM surface under different ionic strength conditions using the CHARMM32 FF (MacKerell et al. Reference MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, Joseph-McCarthy, Kuchnir, Kuczera, Lau, Mattos, Michnick, Ngo, Nguyen, Prodhom, Reiher, Roux, Schlenkrich, Smith, Stote, Straub, Watanabe, Wiórkiewicz-Kuczera, Yin and Karplus1998) for parameters of both the protein and the SAM surface. The catalytic activity of SO was earlier experimentally observed to be dependent on the flexibility of the tether connecting the molybdenum cofactor (Moco) domain and the cytochrome b5 domain of SO. In line with experimental findings, their results indicated that the adsorption of the enzyme's cytochrome b5 domain is inhibited at high ionic strength (750 mM), whereas under much lower ionic strength conditions (100 mM), stable interactions with the surface take place, leading to a loss of flexibility of the tether.

Effects of water on adsorption on SAMs have been reported in several studies. Residence times (τ _F) and self-diffusion coefficients (D_F) of interface water molecules on several different functionalized SAM surfaces of type S(CH₂)-F were inspected by Wang et al. (Reference Wang, Zhao, Zhao, Wang, Yang, Yu and Zheng2010a, Reference Wang, Zhao, Yu, Zhao, Li and Zhengb) and found to be in the order: τ _COOH > τ _NH2 > τ _OH > τ _CH3 > τ _bulk and D_COOH < D_OH < D_NH2 < D_CH3 < D_bulk. This ranking was found to be directly correlated with the number of hydrogen bonds formed between the interfacial water and the modified SAM surfaces. Further, it was shown that the structure and the dynamics of the water determine the adsorption behavior of amyloid-β peptide on the modified SAM surfaces and that the energy barrier for the adsorption of the peptide is lower on hydrophobic CH3-SAM than on the hydrophilic OH-SAM.

The mobility of interfacial water on the zwitterionic sulfobetaine (SBT)-SAM surface was investigated in a separate study by Xie et al. (Reference Xie, Liu and Zhou2012). They found that the lower density packing and the higher flexibility of the SBT-SAM allowed water molecules to penetrate into the surface, thereby increasing the number of interactions of SBT-SAM with the water molecules and decreasing the mobility of the interface water molecules significantly compared with that of pure bulk water (D_SBT < D_bulk/20). The mean force-distance profile showed that the affinity of the SBT-SAM surface towards neuromedin-B peptide was smaller than for the other two SAM surfaces (CH₃-, OH-SAM). The observations support the experimental results showing the antifouling properties of SBT-SAM surfaces (Xie et al. Reference Xie, Liu and Zhou2012).

7.6 sp²-Carbon surfaces

The ability of a biomolecular FF to predict the interaction properties of a single water molecule described by a standard TIP3P or TIP4P model, and several ions on CNT and fullerene surfaces was explored by Schyman & Jorgensen (Reference Schyman and Jorgensen2013). Adsorption energies, as well as the binding site on the surface and the orientation of the water molecule obtained from the simulations using MM FFs, were compared with those of DFT calculations. The non-polarizable OPLS-AA FF was found to be adequate for the description of water properties on benzene (C₆H₆) and coronene (C₂₄H₁₂). However, with increasing surface size above (C₅₄H₁₈) and especially for C₆₀, the simulations became inaccurate, indicating that surface polarization played an important role. Indeed, a polarizable biomolecular FF with induced dipoles on all non-hydrogen atoms (OPLS-AAP) was shown to yield reasonably good agreement with QM calculations for the interaction properties of water molecules and ions.

Another polarizable FF, AMOEBA (Ren & Ponder, Reference Ren and Ponder2002, Reference Ren and Ponder2003), has been extended for the description of the non-bonded interactions between peptides and CNT (De Miranda Tomasio & Walsh, Reference De Miranda Tomásio and Walsh2007) and graphite (Walsh, Reference Walsh2008), and exploited for simulations of peptide binding on CNT (Tomásio et al. Reference Tomásio and Walsh2009). These studies demonstrated that surface polarization caused a strong interaction between the carbon surface and aromatic rings of the side chains of the residues, particularly of tryptophan, placed in the surface plane, which caused strong binding of peptides rich in aromatic residues. The binding of peptides on CNT was also shown to be affected by the interplay between the ability of the aromatic rings to align on a non-flat surface and the total peptide-surface contact area that is increased due to surface curvature. The importance of π-stacking between aromatic side chains and carbon surfaces was also shown in a study of peptide adsorption on a single-layer graphene substrate by Akdim et al. (Reference Akdim, Pachter, Kim, Naik, Walsh, Trohalaki, Hong, Kuang and Farmer2013) using the polarizable AMOEBABIO-09 FF implemented in TINKER 5·1. (Saint Louis, Washington. Software Tools for Molecular Design.). Although these studies demonstrate that polarizable MM FFs are able to provide reasonably accurate descriptions of the interaction properties of peptides on carbon substrates, such simulations are still computationally very demanding for the use of protein adsorption simulations. A viable approach is the inclusion of polarizability in graphene simulations via the rigid rod model, which was employed by Hughes et al. (Reference Hughes, Tomásio and Walsh2014) to calculate the adsorption free energies of single amino acids on a aqueous graphene surface.

The role of the water in peptide binding to a graphene surface was explored by Camden et al. (Reference Camden, Barr and Berry2013), who simulated GXG tripeptides for 20 amino acid residues (X) with an explicit water model. In the TEAM FF used in these simulations, the carbon atoms of the graphene surface were not polarizable, and thus important polarization-driven interactions between aromatic groups and the surface were missing. However, this study demonstrated that desolvation effects may play a very important role in peptide binding to a graphene surface. The dynamics of water on carbon-based hydrophilic graphite-COOH and hydrophobic graphite-CH₃ surfaces was studied by Li et al. (Reference Li, Liu, Li, Ye, Chen, Fang, Wu and Zhou2005) using the standard OPLS-AA FF. They observed that the diffusion of water was slowed dramatically in the vicinity of the surfaces, and this effect extended up to 15 Å from the surface. The water diffusion coefficient within a 3 Å distance from a hydrophilic surface (-COOH) was 4 times smaller than that of a hydrophobic surface (-CH₃).

Implicit solvent models offer an alternative way to explore biomolecule-carbon interfaces in aqueous solution. It can be expected that the hydration shell on carbon surfaces has a much less defined structure than on hydrophilic surfaces and thus, is better suited for a continuum model. Indeed, no strong influence of structural water on the conformation of a test peptide was observed in evaluation of implicit and explicit water models (Walther et al. Reference Walther, Jaffe, Halicioglu and Koumoutsakos2001). In an implicit solvent model, however, one has to account for hydrophobic desolvation effects, which may promote binding of biomolecules to carbon surfaces, explicitly. This was demonstrated by Mereghetti & Wade (Reference Mereghetti and Wade2011) in a study of hydrophobin adsorption on graphite.

Raffaini and Ganazzoli reported simulations of protein adsorption on graphite surfaces in a series of studies using the consistent valence FF (CVFF) (Dauber-Osguthorpe et al. Reference Dauber-Osguthorpe, Roberts, Osguthorpe, Wolff, Genest and Hagler1988) with a Morse potential applied to bonded interactions (Raffaini & Ganazzoli, Reference Raffaini and Ganazzoli2003, Reference Raffaini and Ganazzoli2004a, Reference Raffaini and Ganazzolib, Reference Raffaini and Ganazzoli2006, Reference Raffaini and Ganazzoli2007, Reference Raffaini and Ganazzoli2010). In their simulations, many proteins were observed to denature partially or completely upon adsorption to the hydrophobic surface.

Raffaini & Ganazzoli (Reference Raffaini and Ganazzoli2010) also carried out kinetic calculations for lysozyme spreading on a graphite surface by fitting the time evolution of the interaction energies, protein-surface center of mass distance and components of the radius of gyration of the protein to a function consisting of an exponential term and a stretched exponential term (Kohlrausch function):where γ stands for an arbitrary quantity that decreases with time (t), t ₀ and t ₁ are relaxation times, and A, B, C and δ are fitting parameters. In one simulation, they found that, upon adsorption, lysozyme underwent a liquid-like spreading on the surface suggested by a very fast initial decrease of distance between the center of mass of the protein and the surface with t ₀ = 28 ps and by loss of most of the secondary structure. The initial kinetic step was followed by a longer stage with a t ₁ ~ 1 ns in which all the secondary structure was lost. The authors point out that the fast kinetics of spreading, observed in the study may be attributable to the implicit solvent model used in the simulations.

9.1 Molecular dynamics

Among the different sampling methods available, classical MD is by far the most widely used technique to investigate the characteristics of protein–surface interactions. Examples include simulations for structural refinement of docked complexes (Aliaga et al. Reference Aliaga, Ahumada, Sepúlveda, Gomez-Jeria, Garrido, Weiss-López and Campos-Vallette2011; Alvarez-Paggi et al. Reference Alvarez-Paggi, Castro, Tórtora, Castro, Radi and Murgida2013; Brancolini et al. Reference Brancolini, Kokh, Calzolai, Wade and Corni2012, Reference Brancolini, Corazza, Vuano, Fogolari, Mimmi, Bellotti, Stoppini, Corni and Esposito2015; Imamura et al. Reference Imamura, Kawasaki, Nagayasu, Sakiyama and Nakanishi2007), and for the investigation of the binding orientations of proteins on surfaces (Alvarez-Paggi et al. Reference Alvarez-Paggi, Martín, DeBiase, Hildebrandt, Martí and Murgida2010; Boughton et al. Reference Boughton, Andricioaei and Chen2010; Coppage et al. Reference Coppage, Slocik, Ramezani-Dakhel, Bedford, Heinz, Naik and Knecht2013); the kinetic mechanisms of adsorption (Raffaini & Ganazzoli, Reference Raffaini and Ganazzoli2010); the ET pathways and properties of adsorbed proteins (Bizzarri, Reference Bizzarri2006; Siwko & Corni, Reference Siwko and Corni2013; Utesch et al. Reference Utesch, Sezer, Weidinger and Mroginski2012; Zanetti-Polzi et al. Reference Zanetti-Polzi, Daidone, Bortolotti and Corni2014; Zhou et al. Reference Zhou, Zheng and Jiang2004); the effects of pH (Emami et al. Reference Emami, Puddu, Berry, Varshney, Patwardhan, Perry and Heinz2014b; Imamura et al. Reference Imamura, Kawasaki, Awadzu, Sakiyama and Nakanishi2003; Tosaka et al. Reference Tosaka, Yamamoto, Ohdomari and Watanabe2010; Utesch et al. Reference Utesch, Millo, Castro, Hildebrandt, Zebger and Mroginski2013) and ionic strength (Bizzarri, Reference Bizzarri2006) on adsorption; the role of ions in mediating adsorption (Wu et al. Reference Wu, Chen, Skelton, Cummings and Zheng2013); and structural and energetic aspects of adsorption of proteins on surfaces (Apicella et al. Reference Apicella, Soncini, Deriu, Natalello, Bonanomi, Dellasega, Tortora, Regonesi and Casari2013; Jose & Sengupta, Reference Jose and Sengupta2013; Hoefling et al. Reference Hoefling, Monti, Corni and Gottschalk2011; Hung et al. Reference Hung, Mwenifumbo, Mager, Kuna, Stellacci, Yarovsky and Stevens2011; Kubiak-Ossowska & Mulheran, Reference Kubiak-Ossowska and Mulheran2010a, Reference Kubiak-Ossowska and Mulheranb; O'Mahony et al. Reference O'Mahony, O'Dwyer, Nijhuis, Greer, Quinn and Thompson2013; Vila Verde et al. Reference Vertegel, Siegel and Dordick2009, Reference Vigmond, Iwakura, Mizutani and Katsura2011; Wang et al. Reference Wang, Zhao, Zhao, Wang, Yang, Yu and Zheng2010a, Reference Wang, Zhao, Yu, Zhao, Li and Zhengb; Yu et al. Reference Yu, Becker and Carri2012a; Steckbeck et al. Reference Steckbeck, Schneider, Wittig, Rischka, Grunwald and Colombi Ciacchi2014; Sun et al. Reference Sun, Han, Lindgren, Shen and Laaksonen2014a; Sun et al. Reference Sun, Feng, Hou and Li2014b). Further, conventional MD can be used to simulate physical perturbations, such as mechanical or electrical forces exerted on molecules in experiments. Examples include voltage-biased MD simulations of peptides on a gold surface to validate the signals resulting from a surface enhanced Raman spectroscopy (SERS) experiment (Chen et al. Reference Chen, Cruz-Chu, Woodard, Gartia, Schulten and Liu2012). Mimicking the experimental conditions helps to understand how experiments work and how their results should be interpreted.

Classical MD simulations are very helpful in understanding many molecular phenomena at atomic detail. However, the capabilities and limitations of the MD method, as of other sampling methods, must be known in order to set up a simulation and interpret the results obtained properly. When comparing MD results with the thermodynamic and kinetic properties of the same system obtained from experiments, it should be remembered that a trend observed in an MD simulation may not always agree with the expected average properties. A system kinetically trapped in an energy well, may stay in its metastable state for the rest of the simulation and appear as if in equilibrium (Wei et al. Reference Wei, Carignano and Szleifer2011). In fact, a single simulation corresponds to a single trajectory of a system, and it may not always ensure an adequate sampling. Drawing conclusions from a statistically inadequate sampling will thus often lead to incorrect interpretations. To improve the statistical significance of their sampling, Penna et al. (Reference Penna, Mijajlovic and Biggs2014) chose a direct approach and performed more than 240 explicit solvent MD simulations for the investigation of adsorption mechanisms of a platinum-binding peptide and a gold-binding peptide on neutral Pt(111) and Au(111) surfaces, respectively. From these simulations, they obtained statistics on many binding characteristics, such as anchoring events between each of the peptide residues and the surfaces and the distribution of transition times from the anchoring phase to the initial lockdown phase.

Another important sampling problem is the simulation time required to reproduce an adsorption event on a molecular scale. Even though a simulation that captures the dynamics of a system for several nanoseconds is helpful for structural refinement of molecules or for understanding the initial stages of a process, many molecular events take place over much longer time scales. In their paper, Wei et al. (Reference Wei, Carignano and Szleifer2011) stressed the necessity of very long simulations to be able to investigate a complete adsorption process on a surface in detail. To bypass the energy barriers that normally require long simulation times to overcome, Emami et al. (Reference Emami, Puddu, Berry, Varshney, Patwardhan, Perry and Heinz2014b) employed 1 ns-long annealing simulations at 500 K prior to their production simulations at 298 K. From these classical MD simulations, they investigated the adsorption of three peptides differing in net charge to silica surfaces under four different pH values (with changing ionization levels of the silica surface depending on the pH value). The authors showed that the times spent by each of the peptides on the silica surfaces in the simulations could be used as a relative measure of adsorption strength as a function of pH (see Fig. 10). Even though the trends in adsorption level with pH for each of the peptides from the simulations and the experiments are in good agreement, the time spent on the surface by each of the peptides fails to fully reproduce the difference in relative adsorption levels between the negatively charged and overall neutral peptides measured experimentally.

Fig. 10. Adsorption of three different peptides (positively charged, neutral and negatively charged) on silica nanoparticles as a function of pH. (a) Amounts of adsorbed peptides measured in experiments (Puddu & Perry, Reference Puddu and Perry2012). (b) Times spent by the peptides on the surface calculated from the MD simulations. Reprinted with permission from (Emami et al. Reference Emami, Puddu, Berry, Varshney, Patwardhan, Perry and Heinz2014b). Copyright (2014) American Chemical Society.

To improve the sampling characteristics of conventional MD simulations, replica exchange MD (REMD) (Sugita et al. Reference Sugita and Okamoto1999) was developed based on a MC simulation method called parallel tempering (Hansmann, Reference Hansmann and Okamoto1997) (see the following section). In conventional MD simulations, structures may become trapped in local energy minima on the potential energy landscape, hence, sampling of configurational phase space is limited. To overcome high energy barriers, the REMD method makes use of sampling of the configurational space with independent MD simulations (i.e. replicas) of a system at different replica conditions, e.g. at different temperatures or with a different Hamiltonian. By exchanging the replicas, configurations that are trapped in local energy minima can explore other parts of the energy landscape and an improved Boltzmann-weighted sampling can be achieved. REMD methods are useful to sample a large set of configurations in different potential wells separated by high energy barriers, which otherwise is not possible in a classical MD simulation. Indeed, temperature-based REMD (T-REMD) simulations were used in several studies to accelerate the sampling of peptide–surface interactions (Corni et al. Reference Corni, Hnilova, Tamerler and Sarikaya2013; Li et al. Reference Li, Luo, Derreumaux and Wei2011; Oren et al. Reference Oren, Notman, Kim, Evans, Walsh, Samudrala, Tamerler and Sarikaya2010) and of interactions of basic fibroblast growth factor (bFGF), a small protein, with a hydroxyapatite (001) surface (Liao & Zhou, Reference Liao and Zhou2014). Liao and Zhou observed that while the protein displaces the surface hydration shell and binds tightly to the surface in the T-REMD simulations (five replicas in the range from 310 K to 2500 K), bFGF did not contact the surface directly in classical MD simulations at 310 K (Liao & Zhou, Reference Liao and Zhou2014). However, the T-REMD method is usually not applicable to large systems, e.g. a large protein and a surface with explicit solvent, due to the number of replicas required which scales as O(f ^1/2), with f being the system's total number of degrees of freedom (Fukunishi et al. Reference Fukunishi, Watanabe and Takada2002; Wright & Walsh, Reference Wright and Walsh2013).

Typically, simulation of a protein-surface system in atomic detail requires solvation in tens of thousands of water molecules. The large number of water molecules becomes problematic because the water–water interaction energy term (E _ww) dominates the energy terms due to interactions between the solute molecules and between the solute and water molecules, thereby demanding more replicas in a temperature-based REMD simulation (Huang et al. Reference Huang, Hagen, Kim, Friesner, Zhou and Berne2007). To overcome the poor scaling of the T-REMD method with the size of the system, replica exchange with solute tempering (REST) was proposed (Liu et al. Reference Liu, Kim, Friesner and Berne2005). In this approach, the potential energy function is tailored such that the E _ww term is eliminated from the acceptance probability.

This increases the chances of acceptance of exchanges between replicas compared with regular TREMD simulations. Wright & Walsh (Reference Wright and Walsh2013) investigated the transferability of the REST method to peptide–surface interactions by using a quartz binding protein and a fully hydroxylated α-quartz system as a benchmark for the REST and the T-REMD simulations. The results showed that the REST approach was 50% less expensive in terms of CPU time than the T-REMD method using the same simulation parameters. Furthermore, up to about 80% of the CPU time could be saved with the REST method, owing to the smaller number of replicas used in the simulations. Using the REST parameters optimized for peptide–surface interaction simulations (Wright & Walsh, Reference Wright and Walsh2013), Tang et al. (Reference Tang, Palafox-Hernandez, Law, E. Hughes, Swihart, Prasad, Knecht and Walsh2013) performed REST simulations for numerous gold binding peptides to obtain an accurate ensemble of peptide conformations adsorbed on a gold surface.

The REMD and REST methods sample the phase space with a canonical distribution (fixed number of atoms, volume and temperature) and therefore fall short in exploring the very low probability states. Although investigating the dynamics of a system does not necessarily require an extensive exploration of these states, it may lead to erroneous results, particularly when calculating free energies. This problem can be greater for systems with peptides strongly interacting with surfaces (Wright & Walsh, Reference Wright and Walsh2013). To obtain more exhaustive sampling of the phase space, replica exchange methods can be used in conjunction with methods employing biasing potentials such as metadynamics (Bussi et al. Reference Bussi, Gervasio, Laio and Parrinello2006a).

9.2 Monte Carlo methods

MC methods are stochastic techniques used to solve problems that obey certain probability functions. MC methods have long been used in molecular simulations as an alternative to MD, owing to their efficiency in sampling of conformational ensembles with a Boltzmann distribution. From these ensembles, the geometric averages and thermodynamic properties of a system can be approximated (Schlick, Reference Schlick2002). The MC methods applied to biomolecular systems have been reviewed (see Hansmann & Okamoto, Reference Hansmann1999; Taylor et al. Reference Taylor, Jewsbury and Essex2002), and information on MC methods can be found elsewhere (Landau & Binder, Reference Landau and Binder2009).

Commonly employed for atomistic and CG simulations of protein-surface systems, Metropolis sampling (Metropolis et al. Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953) is a simple and yet very efficient Markov chain Monte Carlo (MCMC) scheme. The Metropolis scheme is used to create a random walk on a set of points with an importance-weight based on a probability distribution (Frenkel & Smit, Reference Frenkel and Smit2001). MC simulations with the Metropolis algorithm have been used to explore the adsorbent conformations and protein denaturation on the surfaces (Anand et al. Reference Anand, Sharma, Kumar and Belfort2009; Castells et al. Reference Castells, Yang and Van Tassel2002; Euston & Naser, Reference Euston and Naser2005; Skepö, Reference Skepö2008; Skepö et al. Reference Skepö, Linse and Arnebrant2006; Zhdanov & Kasemo, Reference Zhdanov and Kasemo2000a), to investigate ordering of the proteins adsorbed on the surface (Zhdanov & Kasemo, Reference Zhdanov and Kasemo2000b), and to calculate the binding free energy of adsorption (Lund et al. Reference Lund, Åkesson and Jönsson2005).

Al-Mekhnaqi et al. (Reference Al-Mekhnaqi, Mayeed and Newaz2009) applied MC simulations with an energy minimization scheme to predict the conformations of a HSA subdomain A on a graphite surface with transition probabilities determined using the Metropolis method. The protein structure was modeled with a united atom approximation, and the structures obtained during MC simulations were minimized using a direct search scheme, which is a technique used for solving optimization problems without constraints. The simulations started with the peptide chain fully extended. In each step, one dihedral angle was chosen randomly and given a random value. After an energy minimization step, a new step was performed and this conformation was accepted with the Metropolis criterion. Simulations were ended after a predefined number of perturbations was achieved. Using this method, they were able to reproduce very similar conformations of the albumin subdomain adsorbed to a graphite surface to those obtained by all-atom molecular simulations (Raffaini & Ganazzoli, Reference Raffaini and Ganazzoli2003).

The Metropolis MC method can fail in sampling the configuration space of a system if the system becomes trapped in a low-energy state at low temperatures. Therefore, typically, many independent Metropolis MC simulations of a system are carried out. Alternatively, an energy-based MC method, parallel tempering Monte Carlo (PTMC), (Swendsen & Wang, Reference Swendsen and Wang1986) can be used. It is the basis for the REMD method. Every time a certain number of MC moves have been performed, the configurations are exchanged between a replica and its neighboring replicas that have the closest temperature values. The exchange operation is ruled by the Boltzmann distribution and done with the Metropolis acceptance probability. The PTMC method has been applied to protein–surface interactions in several studies, including prediction of the binding orientations of proteins (Xie et al. Reference Xie, Zhou and Jiang2010), and investigation of thermal and mechanical properties of proteins tethered on surfaces using a hybrid PTMC method with MD (Knotts et al. Reference Knotts, Rathore and de Pablo2005).

Xie et al. (Reference Xie, Zhou and Jiang2010) compared the performance of the PTMC method with that of conventional serial Metropolis MC simulations of lysozyme on charged surfaces using the same parameters. Figure 11a shows the distribution of the orientations of the lysozyme on a negatively charged surface obtained from the different MC methods. A comparison of these distributions with the adsorption energy landscape (Fig. 11b ) shows that with the PTMC method, the two lowest energy minima were found. In the four serial MC simulations, on the other hand, each simulation revealed a different energy minimum: two local energy minima in addition to the two lower minima found with PTMC.

Fig. 11. (a) Distribution of orientations of lysozyme on a negatively charged surface obtained by a single parallel tempering MC (labeled as p₀) simulation and by four separate conventional serial MC (labeled as s1-s4) simulations. Orientations are represented by the cosine of angle θ, the angle between the unit vector along the dipole moment of lysozyme and the unit vector normal to the surface. (b) Energy landscape of interaction of lysozyme with the surface. Energy minima corresponding to the most visited orientations in the MC simulations are indicated with arrows. Reprinted with permission from (Xie et al. Reference Xie, Zhou and Jiang2010). Copyright (2010), AIP Publishing LLC.

In two other studies, the same group applied the PTMC method to obtain the most probable orientations of fibronectin on a HAP surface (Liao et al. Reference Liao, Xie and Zhou2014) and of protein G on -COOH and -NH₂ functionalized SAMs (Liu et al. Reference Liu, Liao and Zhou2013). These orientations were used as starting orientations for MD simulations to investigate the adsorption mechanisms of the proteins.

Finally, basin-hopping Monte Carlo (BHMC) (Li & Scheraga, Reference Li and Scheraga1987; Wales & Doye, Reference Wales and Doye1997) is worth mentioning, since it is one of the most reliable MC methods for searching the lowest energy configurations of molecules (Iwamatsu & Okabe, Reference Iwamatsu and Okabe2004). Although the BHMC method has been successfully applied to conformational sampling of peptides in implicit solvent (Strodel & Wales, Reference Strodel and Wales2008), to the best of our knowledge, this method has yet to be employed for peptides or proteins near inorganic surfaces.

9.3 Brownian dynamics

The diffusive dynamics of a Brownian particle suspended in a solution can be modeled by the BD technique. BD is commonly applied for simulation of diffusion-driven binding processes leading to the formation of ‘diffusional encounter complexes’. Further information on details of BD methods, and their applications can be found in other reviews (Allison et al. Reference Allison, Northrup and McCammon1986; Gabdoulline & Wade, Reference Gabdoulline and Wade2002; Madura et al. Reference Madura, Davist, Gilson, Wade, Luty, McCammon, Lipkowitz and Boyd1994).

BD has been used in many studies of adsorption of proteins on surfaces with CG (Ravichandran & Talbot, Reference Ravichandran and Talbot2000) and all-atom (Ravichandran et al. Reference Ravichandran, Madura and Talbot2001) representations of the proteins. Simplified BD simulation methods designed for CG models have been proposed (De la Torre et al. Reference De la Torre, Hernández Cifre, Ortega, Schmidt, Fernandes, Pérez Sánchez and Pamies2009; Gorba & Helms, Reference Gorba and Helms2003) and one (Gorba & Helms, Reference Gorba and Helms2003) was tested for the diffusional dynamics of cytochrome c on a charged surface. As an alternative to CG models, several BD studies of protein adsorption to different surfaces have been performed with rigid, atomic detail models of the solutes (Brancolini et al. Reference Brancolini, Kokh, Calzolai, Wade and Corni2012; Mereghetti & Wade, Reference Mereghetti and Wade2011; Kokh et al. Reference Kokh, Corni, Winn, Hoefling, Gottschalk and Wade2010). Mereghetti and Wade (Mereghetti & Wade, Reference Mereghetti and Wade2011) applied BD to the simulation of hydrophobin proteins on a graphite surface. They explained the high affinity of the protein for the graphite surface from the average proximity and orientation of hydrophobin on the surface with its hydrophobic face towards the surface obtained from the simulations.

Brancolini and coworkers employed BD simulations for docking of ubiquitin (Brancolini et al. Reference Brancolini, Kokh, Calzolai, Wade and Corni2012) and of human β2-microglobulin (Brancolini et al. Reference Brancolini, Toroz and Corni2014, Reference Brancolini, Corazza, Vuano, Fogolari, Mimmi, Bellotti, Stoppini, Corni and Esposito2015; Brancolini et al. Reference Brancolini, Corazza, Vuano, Fogolari, Mimmi, Bellotti, Stoppini, Corni and Esposito2015) on bare, citrate covered and thiol-protected gold nanoparticles using the ProMetCS model (Kokh et al. Reference Kokh, Corni, Winn, Hoefling, Gottschalk and Wade2010). The most abundant orientations of ubiquitin on the gold obtained from BD simulations were used as starting orientations for fully flexible MD simulations. By using a combination of BD and MD methods, they were able to reduce the number of possible initial orientations for MD simulations, evaluate the stability of each orientation with MD simulations and observe the dynamics of adsorption in detail.

Overall, BD is an effective method for simulation of biomolecules as well as interactions of molecules with surfaces. An implicit representation of the solvent and a rigid-body treatment of the molecules help to reduce the time required for calculations, allowing thousands of simulations to be performed within a few hours with current computer technology. However, while the rigid-body model may be appropriate for globular proteins, it is less suitable for peptides and unstructured proteins. For these, flexible CG models of the solute can be employed in BD simulations.

10.1 Equilibrium methods

A commonly used equilibrium sampling technique to obtain a free energy profile of a system is umbrella sampling (Torrie & Valleau, Reference Torrie and Valleau1974, Reference Torrie and Valleau1977). In this method, a reaction coordinate is defined that connects two thermodynamic states for which the free energy difference is computed. The reaction coordinate used is usually selected based on a geometric entity, e.g. distance, angle, etc. The reaction coordinate is then divided into distinct windows in which a bias (umbrella) potential is applied to keep the system near that coordinate point by means of restraints. For each window, a separate simulation of the system is performed for sampling around the corresponding coordinate point. Then, the simulations are combined with a reweighting procedure since they are performed in biased ensembles. Two reweighting procedures are weighted histogram analysis method (WHAM) (Kumar et al. Reference Kumar, Rosenberg, Bouzida, Swendsen and Kollman1992; Souaille & Roux, Reference Souaille and Roux2001) and umbrella integration (Kästner & Thiel, Reference Kästner and Thiel2005). Using the final probability distribution of the configurations, the PMF and hence, the free energy profile of the system, is obtained.

The umbrella sampling technique has been successfully employed for the calculation of the free energy of adsorption of a model surfactant on hydrophobic and hydrophilic silica surfaces (Xu et al. Reference Xu, Yang and Yang2008), nanoparticles on phospholipid membranes (Li & Gu, Reference Li and Gu2010), ions on hydrophobic surfaces (Horinek et al. Reference Horinek, Serr, Bonthuis, Boström, Kunz and Netz2008) and amino acids on a ZnO surface (Nawrocki & Cieplak, Reference Nawrocki and Cieplak2013). Umbrella sampling combined with WHAM has also been applied to peptide–surface interactions. Examples include adsorption of various tripeptides to a CH₃-SAM (Sun et al. Reference Sun, Dominy and Latour2007), and an RGD peptide to a titanium oxide surface (Schneider & Colombi Ciacchi, Reference Schneider and Colombi Ciacchi2010). Although, as in these two examples, the distance of the peptide from the surface is typically used as the reaction coordinate for sampling, other definitions can also be employed. For instance, Boughton et al. (Reference Boughton, Andricioaei and Chen2010) used the angular orientation of the α-helix of magainin 2, an α-helical amphiphilic peptide, on a polystyrene surface as a reaction coordinate to compute its adsorption free energy by umbrella sampling combined with WHAM.

Umbrella sampling, even though it restricts the spatial sampling necessary to calculate free energy differences, can still be subject to conformational sampling problems. Wang et al. (Reference Wang, Stuart and Latour2008) therefore combined umbrella sampling with a biased REMD technique. In this method, an initial umbrella sampling of the system is employed to obtain a rough estimate of the PMF profile of the system. The negative of the PMF is then used as a biasing potential for REMD simulations. A final PMF profile is constructed from a probability-ratio method analysis for the REMD simulations. Comparison of the estimated PMF profile of a Na⁺ ion on a carboxylic acid functionalized SAM surface with the theoretical profile, showed that the biased REMD gives much better agreement with the theoretical profile compared with those obtained from conventional REMD and umbrella sampling simulations. The method was applied in separate studies to adsorption of the G₄-K-G₄ peptide on a polylactic acid surface for the investigation of binding free energy profiles (O'Brien et al. Reference O'Brien, Stuart, Bruce and Latour2008), for the assessment of the protein FF parameters to be used in protein-surface simulations (Vellore et al. Reference Vellore, Yancey, Collier, Latour and Stuart2010), and for investigation of the effect of pressure on the calculation of protein-surface binding free energies from MD simulations (Yancey et al. Reference Yancey, Vellore, Collier, Stuart and Latour2010).

Another common approach used in molecular simulations to calculate the free energy landscape as well as investigate new reaction pathways is metadynamics (Bussi et al. Reference Bussi, Laio and Parrinello2006b; Laio & Gervasio, Reference Laio and Gervasio2008; Laio & Parrinello, Reference Laio and Parrinello2002). In this approach, the dynamics of the system is steered by thermodynamic forces and compensated by a history-dependent biased potential (Laio & Parrinello, Reference Laio and Parrinello2002). Owing to the bias potential, the frequency at which a state is visited during a simulation typically decreases linearly with its free energy (Bussi et al. Reference Bussi, Gervasio, Laio and Parrinello2006a). Accurate free energy profiles can be obtained from metadynamics with an appropriate reweighting scheme (Bonomi et al. Reference Bonomi, Barducci and Parrinello2009) to recover an unbiased probability distribution (Meissner et al. Reference Meißner, Schneider, Schiffels and Colombi Ciacchi2014). To date, metadynamics has been successfully applied for calculation of the free energy of binding of a formate anion to a TiO₂ surface (Mori et al. Reference Mori, Hamers, Pedersen and Cui2013) and of amino acids to silver and gold surfaces (Palafox-Hernandez et al. Reference Palafox-Hernandez, Tang, Hughes, Li, Swihart, Prasad, Walsh and Knecht2014), and for calculation of the free energy landscapes of alanine dipeptide in its free and gold surface bound states from 2μs-long simulations (Bellucci & Corni, Reference Bellucci and Corni2014).

The combination of metadynamics with other enhanced sampling techniques allows a better exploration of low probability states, hence increasing the accuracy of the free energy calculations (Bussi et al. Reference Bussi, Gervasio, Laio and Parrinello2006a). Qin & Buehler (Reference Qin and Buehler2014) investigated the adsorption free energy of mussel adhesion proteins onto a silica surface using metadynamics starting from configurations obtained from REMD simulations.

Schneider & Colombi Ciacchi (Reference Schneider and Colombi Ciacchi2012) have successfully applied a hybrid REST and metadynamics to study the adsorption of small peptides on Si and Ti surfaces. With only four replicas at temperatures of 300, 350, 400 and 450 K, they were able to obtain agreement with the experimental adsorption free energy of the RKLPDA peptide on a Ti surface (experimental: 38·0 ± 8 kJ mol⁻¹, computed: 38·6 ± 3·9 kJ mol⁻¹). The same approach was employed recently by Meissner et al. (Reference Meißner, Schneider, Schiffels and Colombi Ciacchi2014) to estimate the circular dichroism (CD) spectra of a peptide in its free state and bound to a silica surface. CD spectroscopy is a useful technique to monitor secondary structure content in biomolecules, in particular, α-helix content in peptide structures. Meissner et al. calculated the CD ellipticity values of conformations from simulation snapshots and used these values as external collective values in their reweighting procedure. The estimated fractional helicity values of the free and adsorbed peptides obtained from the simulations show good agreement with experiments. Even though these results are very promising, application of this method to protein adsorption requires careful selection of parameters, e.g. adequate number of replicas, appropriate selection of temperatures.

Thermodynamic integration is another approach used in free energy calculations of protein–surface interactions. This method, like umbrella sampling, requires a predefined reaction path between the initial and the final states. Juffer et al. (Reference Juffer, Argos and de Vlieg1996) applied thermodynamic integration to compute the adsorption free energies of the enzyme cutinase and its variants to a charged surface in the presence of explicit ions. Hoefling et al. (Reference Hoefling, Iori, Corni and Gottschalk2010a) applied the method to obtain the PMF profiles of amino acids adsorbing on an Au(111) surface. Further, Schneider & Colombi Ciacchi (Reference Schneider and Colombi Ciacchi2010) employed thermodynamic integration, along with the umbrella sampling/WHAM method, to compare the free energy profiles of adsorption of RGD peptides on oxidized Ti obtained by the two methods. Their results showed a perfect agreement between the results of the two methods. Finally, Friddle et al. (Reference Friddle, Battle, Trubetskoy, Tao, Salter, Moradian-Oldak, De Yoreo and Wierzbicki2011) used the adaptive biasing force simulation method, which is a technique based on thermodynamic integration to obtain free energy profiles using a biasing force, for the adsorption of a 12-mer C-terminal fragment of the amelogenin protein on various different crystal terminations of two different hydroxyapatite surfaces ((100) and (001)). By complementing AFM measurements with the free energy calculations, Friddle et al. (Reference Friddle, Battle, Trubetskoy, Tao, Salter, Moradian-Oldak, De Yoreo and Wierzbicki2011) were able to predict the crystal termination of the hydroxyapatite surface used in experiments and identify the interactions governing the adsorption of the amelogenin peptide to it.

10.2 Non-equilibrium methods

Observing adsorption/desorption events of large molecules using conventional MD methods requires very long simulation times. These limitations can be overcome by steered molecular dynamics (SMD) simulations (Utesch et al. Reference Utesch, Daminelli and Mroginski2011), which is the most commonly used non-equilibrium simulation technique for biomolecular systems. In a typical SMD simulation, a protein or a peptide is pulled with a non-physical external force along a predefined reaction coordinate with a constant velocity or with a constant force, thereby accelerating adsorption or desorption events. Further, SMD allows simulation of systems under mechanical stress such as stretching, shearing and bending, and can hence be used to predict the effects of external disturbances on protein-surface systems (Hamdi et al. Reference Hamdi, Ferreira, Sharma and Mavroidis2008). These disturbances may lead to uncertainties in experimental measurements. Particularly, SMD simulations conducted to investigate the effect of the tip of an AFM showed that pushing a molecule towards a surface with an AFM tip will bias the results by increasing adhesion, as the force exerted on the molecule leads it to spread more on the surface (Horinek et al. Reference Horinek, Serr, Bonthuis, Boström, Kunz and Netz2008; Mücksch & Urbassek, Reference Mücksch and Urbassek2011). With SMD simulations, the adsorption mechanisms of peptides and proteins on different types of surfaces have been investigated (Alvarez-Paggi et al. Reference Alvarez-Paggi, Martín, DeBiase, Hildebrandt, Martí and Murgida2010; Dong et al. Reference Dong, Wang, Wu and Pan2007; Emami et al. Reference Emami, Puddu, Berry, Varshney, Patwardhan, Perry and Heinz2014b; Friedrichs et al. Reference Friedrichs, Köppen and Langel2013; Hamdi et al. Reference Hamdi, Ferreira, Sharma and Mavroidis2008; Shen et al. Reference Shen, Wu, Wang and Pan2008; Utesch et al. Reference Utesch, Daminelli and Mroginski2011; Yang & Zhao, Reference Yang and Zhao2007) and compared with experimental measurements (Schneider & Colombi Ciacchi, Reference Schneider and Colombi Ciacchi2010, Reference Schneider and Colombi Ciacchi2012). In addition to accelerating the sampling of certain events, the SMD method is used to calculate free energy differences. It was shown by Jarzynski (Reference Jarzynski1997) that a free energy difference can be obtained from a non-equilibrium process following the equation: e ^{− βΔF}  = 〈e ^−βW 〉 where β = 1/(k _B T) (k _B is the Boltzmann constant and T is the absolute temperature) and ∆F and W stand for the free energy difference and the non-reversible work done along a reaction coordinate, respectively. The average is taken on all (in principle) the realizations of the non-equilibrium process that connect the initial and final states of interest. Based on this equation, a free energy difference can be obtained using a population of SMD simulations that start from an equilibrium ensemble. Binding free energies of peptides and proteins to various different surfaces computed using the Jarzynski equation have been reported in a number of papers (Chen et al. Reference Chen, Wang, Liu, Wu, Kang, Moore and Gubbins2009b; Kang et al. Reference Kang, Liu, Wang, Shen, Wu and Guan2009; Mijajlovic et al. Reference Mijajlovic, Penna and Biggs2013). Baier et al. (Reference Baier, Blumenstein, Preusker, Jeurgens, Welzel, Do, Pleiss and Bill2014) investigated the adsorption free energies of ZnO-binding peptides on ZnO surfaces. To this end, they applied a hybrid steered MD with umbrella sampling method. In this approach, the peptides were steered towards the surface with a constant force applied to the peptides. After an equilibration period without any external forces applied upon adsorption, the peptides were pulled back from the ZnO surface. Selected conformers were then simulated for another 5 ns with a harmonic potential applied to their centers of mass. The PMF profiles were obtained using the WHAM method. Combined SMD simulations with umbrella sampling were also used in a study to calculate the adsorption free energy of a cationic peptide on a silica surface (Emami et al. Reference Emami, Puddu, Berry, Varshney, Patwardhan, Perry and Heinz2014b).