Carbon nanotubes (CNTs), which are entirely composed of sp2 carbons, can be viewed as graphite sheets that have been rolled into seamless cylinders. Since their discovery in 1991 , CNTs have become the “superstar” in the field of nanoscience and nanotechnology because of their unique structural, mechanical, and electronic properties . Various promising biomedical applications have also been proposed, such as drug design , drug delivery , tumor therapy , tissue engineering , DNA recognition , and biosensor design . It was predicted that the global market for CNTs could grow to $2 billion by 2014 .
On the other hand, with the fast advance of the nanotechnology and its widespread applications in various industries, there is a growing concern on the biosafety of these nanomaterials to human health [10-17]. Human skin, lungs, and gastrointestinal tract are in constant contact with the environment. Despite the skin is an effective barrier to foreign substances, the lungs and gastrointestinal tract are vulnerable. Recent experiments showed that unrefined CNTs might be aerosolized and release fine particles into air . When the concentration of CNTs in air is sufficiently high, they can be inhaled and then migrate into other parts of body . These inhaled CNTs can enter cells and accumulate in cytoplasms [20, 21], which may lead to the lung insult [22-24], immunologic toxicity , and adverse cardiovascular effects (on mice) . Injections and implants, which are required for some biomedical applications, are the other possible routes of the exposure of human cells to nanoparticles. The experiments on mice showed that the tissue exposing directly to the injected CNTs may lead to asbestos-like pathogenic behaviors, including inflammation and the formation of lesions known as granulomas . In addition, the injecting CNTs into the tail vein of male mice may cause the reversible testis damage . Therefore, there is an urgent need for the understanding of nanotoxicity and the development of biocompatible nanomaterials.
Recently, the interactions between nanoparticles, such as CNTs, and biological molecules, such as proteins, have received great attention due to the aforementioned biosafety concerns. Proteins are the functional units of life. Studies on the interaction of CNTs and proteins may provide a key to understand the basic questions in the nanotoxicology and nanopharmacology. In 1999, Mioskowski and coworkers observed that the protein streptavidin can grow into helical crystals on the surface of CNTs in experiment . In 2003, Park et al. showed experimentally that the CNT is effective in blocking some biological membrane ion channels due to its structural match in size and shape . Using atomic force microscope (AFM) and Fourier transform infrared spectroscopy (FTIR), Karajanagi et al. investigated the conformational changes of two enzymes, α-chymotrypsin (CT) and soybean peroxidase (SBP), after adsorbing onto the SWCNT . They argued that the distribution of the hydrophobic residues on the surface of the enzymes may determine the protein conformational changes when adsorbed onto the CNT. The important role of hydrophobic residues during the interactions between proteins and CNTs was also found by the experiments of Goldberg-Oppenheimer and Regev , in which they explored the conformational changes of bovine serum albumin (BSA) bound with CNTs via Cryogenic temperature transmission electron microscopy (Cryo-TEM). And the experiments on the interactions between small peptides and CNTs also showed that the tryptophan residues play a key role in the binding of the peptides on CNTs [33, 34].
In addition to these studies mentioned above, various experimental techniques have also been developed to shed light on this challenging problem [35-37]. However, the mechanism, particularly the effect of CNTs on the structure and function of biomolecules which may lead to the loss of the original function of proteins, is far from being understood. The inherent difficulties involved in these complex settings can limit their applications due to too many factors involved (“side effects”) and experimental results may vary for the same proteins interacting with the same nanoparticles (such as same chemical composition, same shape, length, and aggregation property, etc.) under same conditions (such as same concentration, temperature, and pressure, etc). The computational simulations, on the other hand, might mimic and model this kind of interactions and avoid some complicated side effects. Molecular dynamics (MD) is one of these computational techniques, which is widely used in the studies of biomolecules [38-55], and nanoscale systems [56-60]. It is found effective in providing insights to the interactions between proteins and CNTs. For example, a recent MD simulation revealed some stepwise conformational changes of the sub-domain of human serum albumin (HSA) due to its adsorption onto the CNT surfaces .
In this chapter, we will review some of our recent works on the interactions between single-wall carbon nanotubes (SWCNTs) and proteins, and their effects on protein structure and function, using large scale molecular dynamics simulations. It is widely accepted that specific interactions between proteins and other molecules, often based on the structures of active sites, play a vital role in the function of proteins. Thus, if a SWCNT can occupy or disrupt the active site of a protein, it will very likely affect the function of the protein. First, we investigated the probability of the SWCNT disrupting the active site by using a signaling and regulatory protein WW domain (a relatively small protein) as an example. It was found that the SWCNT could plug into the hydrophobic cores of WW domains. Two key residues forming the functional scaffold in the native structure were separated by the SWCNT. This plugging of the SWCNT seriously broke the structure of the active site and reduced the possibility of the WW domain in binding with its partners. Second, we investigated the three-way binding competition among a SWCNT, protein SH3 domain (another small signaling and regulatory protein domain), and a proline rich motif ligand. It was found that the SWCNT had a very high probability occupying the active site of the SH3 domain, which prevented the ligand from binding to the active site. In these two cases, the effects of the SWCNT on the proteins would lead to the loss of the original function of the proteins. Finally, the interactions between the SWCNT and human blood serum proteins, including BSA (bovine serum albumin), γ-Ig (gamma-globulin), Tf (transferring), and BFg (fibrinogen), were also investigated by MD simulations, to confirm that the mechanisms identified with these small protein domains, such as WW domains and SH3 domains, also apply to much larger proteins. It is found that the hydrophobic interactions between the SWCNTs and nonpolar residues, particularly those aromatic ones through so-called π−π stacking, play a key role in the interaction between the SWCNTs and proteins.
2.1. Disruption of protein active sites
Proteins are molecular machines, building blocks, and arms of living cells. The enormous variety of protein functions is based on their high specificity in partnering with other molecules, with which they interact -- a relationship that resembles a key and lock. This specific relationship demands a fairly unique spatial structure of the protein, in which the hydrophobic residues are buried in the core (named hydrophobic core). The pristine CNTs are extremely hydrophobic. Thus it is possible for the pristine CNTs to plug into the hydrophobic core due to the strong affinity between the CNT and the hydrophobic residues in the core. This plugging will lead to the loss of the original function of the protein, and thus is a potential mechanism for the toxicity of the CNTs (and maybe other hydrophobic nanoparticles).
We take a WW domain interacting with various sized SWCNTs as an example to illustrate this idea . The WW domain is a protein module in signaling and regulatory proteins as the functional module to identify and bind the proline-rich motifs (PRMs) of its binding partner [63-67]. They exist as a triple stranded anti-parallel β-sheet structure [68-70]. Two highly conserved aromatic residues, a tryptophan residue in the third β-strand and an aromatic residue in the second β-strand, are together in the native structure to form the scaffold to bind the proline residue in the PRM. For the WW domain used in our simulation, which is a module of the human Yes-associated protein , the two key residues are W39 at the third β-strand and Y28 in the second β-strand. The SWCNTs used in the simulation are armchair of (m, n), where m = n = (4, 5, 6), corresponding to the tube diameters of (5.38 Å, 6.73 Å, 8.08 Å), respectively. The SWCNT and protein are initially well separated, with a distance between the geometric centre of protein domain and the SWCNT of 15.0 Å. Each system was solvated with TIP3P model water. We performed 36 independent simulations for the YAP65 WW domain and these three sizes of SWCNTs with different initial velocities, each with 200 ns (see more detail in Computational Methods).
Interaction of YAP65 WW domain with SWCNT
In our simulations, we found several cases that the SWCNTs plugged into the YAP65 WW domain between the second and third β-strands, as one representative structure shown in Fig. 1(a). Fig. 1(b) shows the superposition of the structure of the WW domain in the complex and in the native state. It was found that the main change of the protein conformation in the complex was its third β-strand. It unfolds into a loop and wraps on the SWCNT. Most of the contacts between the second and third β-strand are broken. Essentially, these conformational changes of the WW domain are pre-requirements for the formation of this complex of the protein domain and the SWCNT. In the complex, nine residues in the second and third β-strands contact with the SWCNT. Most of them are hydrophobic residues, such as W39. This implies that the hydrophobic interaction, the dominant force for protein folding , is the key to this phenomenon, in which the SWCNT plugged into the hydrophobic core of the WW protein domain to form a stable protein-SWCNT complex.
Fig. 2(a) displays some representative snapshots of these complexes at different times to show the plugging process of the SWCNT. The interface area between the protein domain and the SWCNT (denoted by S, shown in Fig. 2(b)) is used to illustrate this process. Here S is defined as half of the difference between the solvent accessible surface area of the complex and the sum of solvent accessible surface areas of the protein and the SWCNT . Thus at t = 0, S is small (~50 Å2) since the YAP65 WW domain and the SWCNT are initially separated. The YAP65 WW domain and the SWCNT approached each other very quickly. As a result, S rose to about 200 Å2 within only 1 ns. This indicated that the SWCNT was adsorbed onto the protein surface (see the snapshot of t = 1 ns in Fig. 2(a)). In the next 45 ns, S gradually increased from ~200 Å2 to ~250 Å2, with significant fluctuations along the way. It was found that the contacting surface region of the domain kept changing during the process of adsorption. That is, the SWCNT was constantly seeking for a more stable binding site. Around t = 45 ns, there is a negative impulse in the interface area S. A careful examination on the trajectory showed that the SWCNT began to plug into the second and third β-strands and the negative impulse was a result of the opening of the second and the third β-strands before the SWCNT was “swallowed” by the protein. S reached its maximal value of ~300 Å2 at around t = 80 ns, meaning that the SWCNT was finally successful in plugging into the protein domain. In these structures, the SWCNT was basically wrapped by the second and the third β-strands, forming a complex of the WW domain with the SWCNT. From t = 80 ns up to 200 ns, there were only minor fluctuations of S at ~300 Å2, showing the stability of this protein-SWCNT complex (data after 100 ns were not shown since the plugging was mainly done in the first 100 ns).
In order to investigate the conformational change of the YAP65 WW domain, we computed the RMSD from its native structure in the MD trajectory, as shown in Fig. 2(c). As a comparison, we had also performed simulations for the YAP65 WW domain only (without any SWCNT or ligand). It was found that without the disruption of the SWCNT, the RMSD of the WW domain remained at a value of around 2.0 Å. From the RMSD shown in Fig. 2(c), the conformational change of the YAP65 WW domain can be described as four stages. In the first stage, from 0 to ~45 ns, the RMSD maintained at roughly 2.0 Å. That is, there was little impact of the SWCNT on the conformation of the YAP65 WW domain despite its adsorption on the SWCNT surface. In the second stage, from ~45 to ~70 ns, the RMSD increased gradually from ~2.0 Å to ~3.5 Å, which was a notable conformational change and corresponded to the opening of the second and third β-strands due to the insertion of the SWCNT. The third stage is from ~70 to ~90 ns. In this stage, a small plateau for the RMSD, around 3.0 Å, appeared. It was due to the partial recovery of the contacts between the second and third β-strands. Around 90 ns, there was a significant jump in the RMSD, from ~3.0 to ~4.0 Å, which initiates the fourth stage. Detailed studies showed that there was a change in both the orientation and binding site of the indole ring of W39 residue on the SWCNT during this jumping process. From t = ~95 ns, the RMSD kept roughly constant at ~4.0 Å, indicating the final complex of the YAP65 WW domain and SWCNT was quite stable. This relatively large RMSD indicates a significant conformational alteration of the YAP65 WW domain caused by the insertion of the SWCNT.
The W39 and Y28 are the key residues for the function of the YAP65 WW domain. The interaction between the SWCNT and these two key residues were also investigated. Fig. 2(d) shows the distances between these two residues, and distances between each residue with the SWCNT as a function of time, denoted by dW39,Y28, dW39,SWCNT, and dY28,SWCNT, respectively. Corresponding to the adsorption of the SWCNT on the YAP65 WW, there is a drop of dW39,SWCNT and dY28,SWCNT in the first ns. And then dW39,Y28, dW39,SWCNT, and dY28,SWCNT all kept almost as a constant at ~6.0 Å with small fluctuations until t = 45 ns, showing that the two key residues were “confined” in their native state, and the initial adsorption of the SWCNT does not impact much on their structures and orientations. In the period of t = 45 - 50 ns, there were an impulse jump in dW39,Y28 to a value around 9 Å and a drop in dY28,SWCNT to a valley around 4 Å. These changes indicate the detachment of W39 with Y28 and the approach of Y28 to the SWCNT. However, dW39,SWCNT decreased sharply to < 4 Å at the same time and remained nearly unchanged for the rest of the simulation times, which indicates a strong and favorable interaction between W39 and the SWCNT. It is interesting to note that this favorable tryptophan-CNT interaction has also been observed in recent experiments for CNT-peptide interactions [33, 34]. Therefore, even though we do not have direct experimental evidence at this point for this plugging of SWCNTs, it might not be totally impossible due to the strong and favorable interactions between CNTs and hydrophobic residues with aromatic rings such as tryptophan residues. It was very interesting to find that both dY28,W39 and dY28,SWCNT returned to their previous values (before t = 45 ns) at ~50 ns. A careful examination of the trajectory movie showed that there was a rotation of the SWCNT orientation, with the SWCNT firmly seizing the W39 in the period of t = 50 -70 ns. After t = 70 ns, dY28,SWCNT decreased to a small value again (~4 Å), and the dW39,Y28 turned to a significantly larger value (8-12 Å). This indicates a significantly larger separation of the two active-site residues, which corresponds to a wide opening of the second and third β-strands, indicating a potential loss of its function.
Stability and hydrophobic interaction
We had computed various interaction energies of the system as functions of the interface area between YAP65 WW domain and SWCNT (shown in Fig. 3). For the interaction energy of the total system (using the initial state as reference, see Fig. 3(a)), there is a basin at about S = 280 Å2, with the interaction energy about 30.0 kcal/mol lower than the reference state, the state where YAP65 WW domain and the SWCNT are well separated. From Fig. 2(b), we know that S ~ 280 Å2 corresponds to the state in which SWCNT plugs into the YAP65 WW domain. Therefore, the complex with SWCNT inserted into the core of YAP65 WW domain is the energetically more favored state.
Fig. 3(b-d) show the interaction energies of the solute (YAP65 WW domain and SWCNT), the solvent (the water molecules), and between solute and solvent, as functions of the interface area, respectively. As shown in these figures, the internal interaction energies of the solute and the solvent (see fig. 3(b and c)) change in the opposite direction of the interaction energy between the solute and the solvent (see fig. 3(d)). Clearly, in our current system, the internal interaction energies of the solute and the solvent dominate the total energy of the system, which both have an energetic basin around S=280 Å2. This implies that the solute and the solvent (water) are favored to interact with themselves and not favored to interact with each other, i.e., a separation of the two is energetically more favored, which is best described as so-called hydrophobic interactions. It is the hydrophobic interaction between the YAP65 WW, SWCNT and water that makes the plugging of SWCNT into YAP65 WW domain.
Binding site and structural flexibility
In each case of a SWCNT plugging into the YAP65 WW domain, we have found that the YAP65 WW domain displays the same binding site after relentless search. It was shown above that the hydrophobic interaction dominated the binding and plugging of the SWCNT. However, there are also hydrophobic residues, including a tryptophan residue, on the first β-strand of the YAP65 WW. Thus the question was why the SWCNT plugged into the gap between the second and third β-strand instead of that between the first and second β-strand. The flexibility of different part of the YAP65 WW domain may be a factor of selection of the plugging site. To illustrate this idea, Fig. 4 shows the root mean square fluctuations (RMSF) of the average of non-hydrogen atoms of each residue for the YAP65 WW domain (the simulation was performed without any SWCNT or ligand). Generally, the loops and turns have larger fluctuations. The notable point of this figure is the difference between the three β-strands. Clearly, the fluctuations of the third β-strand residues are significantly higher than those of the other two β-strands. This indicated that the third β-strand was more flexible than the other two β-strands. That is, the gap between the second and third β-strand is easier to be opened, and thus the SWCNT is easier to penetrate it.
Effect on protein function
As shown above, the SWCNT plugged into the second and third β-stands and formed a stable complex. In the protein-SWCNT complex, two active residues of YAP65 WW domain (W39 and Y28) are separated, thus the primary functional unit, the scaffold, for binding with prolines in the proline rich motifs (PRM) is disrupted. In general, the SWCNT insertion of the binding site is much more disruptive than the typical surface adsorption, which may result in the loss of the original function of the WW domain. To demonstrate this, we performed simulations for the system with an additional PRM ligand (sequence GTPPPPYTVG), which binds to both the YAP65 WW domain in its native state , and to the protein-SWCNT complex. In the PRM binding simulations, the PRM was initially placed at 25 Å away from the center of the protein-SWCNT complex. There were 7045 TIP3P water molecules and one Cl- (which was used to neutralize the system), and all other MD setups and parameters were similar to the previous simulations. Ten independent runs were performed for each system in NVT ensemble at 298 K.
Our simulations showed that the protein-SWNT complex was stable, and the PRM always bound to the SWCNT instead of the binding site of the YAP65 WW, as shown in one representative structure in Fig. 5(a). As for comparison, we also showed the native structure with the PRM in Fig. 5(b). Clearly, the binding site (labeled by the green stick in these two figures) had been blocked by the SWCNT in the YAP65 WW domain-SWCNT complex. These findings indicate that the insertion of the SWCNT disrupted and blocked the binding site of the YAP65 WW domain, and disabled its function of identifying and binding with the PRM. To validate this prediction from our in silico approach, in vitro and/or in vivo experiments are highly needed. Interestingly, there is recent evidence from both experiment and theory showing that some antibody from mouse immune repertoire can “absorb” the C60 nanoparticle through similar hydrophobic interactions [75, 76].
2.2. Competitive binding with ligands to receptors
The scheme of protein functioning can be roughly described into three steps: binding, transforming, and releasing. The binding, especially the specific binding, between proteins and ligands is the primary step for the protein functioning. Again, hydrophobic interactions often play a vital role in the binding of proteins with ligands. If pristine CNTs (which are extremely hydrophobic) pass through cell membranes and possess stronger binding affinity than ligands with hydrophobic binding sites of some key proteins, these CNTs may hold their hydrophobic binding sites, obstruct their binding processes, and disrupt the functions of these proteins. This competitive binding may be another potential mechanism of toxicity of hydrophobic nanoparticles, including CNTs. All the aforementioned studies focus on the direct interactions of nanoparticles with proteins. The effect of nanoscale particles on the ligand-receptor binding, i.e., a three way nanoparticle-ligand-protein competitive binding, is much less studied.
We take a proline-rich ligand (PPPVPPRR) and its binding module, the SH3 domain, together with a pristine SWCNT as an example to illustrate the idea. The SH3 domain is one of the most abundant protein interaction domains , and usually found in the signaling and regulatory proteins as the functional module to identify and bind PRM of their binding partners [78-84]. Fig. 6(a) shows the native structure of the SH3 domain binding with the ligand. It exists as a characteristic β-barrel fold which consists of five β-strands arranged as two tightly packed anti-parallel β-sheets . The residues, including F8, W36, P50, and Y53, are the key residues to bind ligands [78, 86]; and the hydrophobic interaction is the dominant force for the binding process . The SWCNTs are armchair of (m, n), where m = n = 3, corresponding to the tube diameters of 4.04 Å with a length of 19.54 Å. The ligand, SWCNT and protein are initially well separated from each other. The initial distances of the geometric centers, of both the SWCNT and the ligand, from that of the SH3 domain are set at 30 Å, but in different orientations. The distance between the geometric centers of the SWCNT and ligand is about 42 Å. They were solvated with TIP3P model water as well. We performed 20 independent simulations for the SH3 domain, the PRM ligand, and the SWCNT with different initial velocities, each with 200 ns (see more detail in Computational Methods).
Interaction of SH3 domain, ligand and SWCNT
We find that the SWCNT occupies the binding pocket and contacts directly with the key residues of the SH3 domain in 13 of the 20 simulation trajectories. In this way, the SWCNT prevents the ligand from forming the native contacts with the SH3 domain binding pocket.
Fig. 6(b) shows a typical structure for the case that SWCNT occupies the binding pocket of SH3 and the ligand binds to the SWCNT. The SWCNT directly contacts with nine residues of the SH3 domain in the pocket, including the key residues W36, P50, and Y53. In the other 12 cases, the SWCNT contacts with at least one of the key residues of the binding pocket of the SH3 domain. And in 11 of those 13 trajectories, the ligand binds partly to the SWCNT and partly to the surface of the SH3 domain (near the SWCNT). In other two cases, the ligand only binds to the surface of the SH3 domain (but not the active site) and does not bind to the SWCNT. Even though the other seven trajectories do not show a direct binding of the SWCNT with the SH3 domain binding pocket, the SWCNT does display some tendency of the binding. For example, in three of the seven trajectories, the SWCNT is getting close to the key W37 residue in the binding pocket; and in the other two trajectories, it is getting close to the key Y03 residue. The incomplete binding might be caused by the fact that our simulation lengths, 200+ ns, are still not long enough.
The binding process is illustrated by representative snapshots as well as the interface areas between the SH3 domain, the ligand, and the SWCNT, using a typical trajectory (see Fig. 7), in which the SWCNT first binds onto the SH3 and then the ligand binds onto the SWCNT.
To illustrate the binding process of the ligand, we also show the interface areas of the ligand with the SH3 domain (protein-ligand), and the ligand with the SWCNT (ligand-SWCNT), denoted by SPL and SLC, respectively (see Fig. 7(b)). At t = 0, SPL = SLC = 0 since the ligand is separated from both the SH3 domain and the SWCNT initially. In the first 5 ns, SPL increases to about 200 Å2 very quickly, and maintains that value until t = ~ 90 ns. At the snapshot of t = 40 ns, the ligand is adsorbed onto the SH3 domain (but not at the binding pocket). On the other hand, SLC still keeps at about 0. The most notable event happened at about 90 ns. After that, SPL decreases quickly from 200 Å2 to a very small value, and later recovers back to ~150 Å2. Meanwhile, SLC increases to about 100 Å2. At t = 89 ns, the ligand has temporarily separated from the SH3 domain, diffuses to the other side of the SH3 domain, and then binds to both the SWCNT and the SH3 domain (the structure is shown in the snapshot at 100 ns). From t = 90 ns to t = 200 ns, both SPL and SLC remain at these values with normal fluctuations, indicating a stable three-way bound structure among the SWCNT, the ligand, and the SH3 domain.
From the interface areas, we showed the main picture of the binding process of the SWCNT, the ligand and the SH3 domain. To give a more detailed picture about the binding process at the residue level, we had computed the contact number of the SWCNT and the ligand with the residues of SH3 domain. Here we used the protein residue as a unit base in accounting (see Fig. 8), and a contact is counted if the distance between a non-hydrogen atom in the
bound object and a non-hydrogen atom in a residue is less than 5 Å. Generally speaking, the closer the object is to the protein residue, the larger the contact number is between the protein and the object, and the stronger the object binds to the protein residue. Similarly, from the contacts of the SWCNT on the SH3 domain, we found that there were three stages in the binding process. In the initial stage, from 0 to ~3 ns, the SH3 domain and the SWCNT were separated; in the second stage, from ~3 to ~35 ns, the SWCNT was adsorbed onto four residues of the SH3 domain; in the third stage, from ~35 ns to the end of the simulation, the SWCNT bound closely with many residues (up to 9) of the SH3 domain, including the key residues W36, P50, and Y53.
Furthermore, the contacts between the SWCNT and the ligand are found to be quite stable, i.e., these contacts remain in contact once they appear (see Fig. 8). On the contrary, the contacts between the ligand and the protein (both ligand on protein and protein on ligand) vary with time, indicating that the relative positions between the ligand and the SH3 domain change with time much more often. Further, the great changes of the ligand-protein contacts happen at t = ~88 ns (Fig. 8), corresponding to the previous observation that the ligand temporarily separates from the SH3 domain, then diffuses to the other side of the SH3 domain, and finally re-adsorbs onto both the SWCNT and the SH3 domain. Conformational fluctuations of the ligand on the SWCNT and SH3 domain are observed after the final adsorption. For example, from 90 ns to 150 ns, the ligand binds both the SH3 domain (on the residues 12-13, 16, and 36) and the SWCNT by the C-terminal residues; while from 150 ns to 200 ns, the ligand binds the SH3 domain (on the residues 9~11, and 53) by the N-terminal and the SWCNT by C-terminal (see also Fig. 8).Free energy landscape
Fig. 9(a) shows the free energy landscape as a function of two reaction coordinates: the interface area between the SH3 domain and the SWCNT (SPC), and the minimal distance between the SWCNT and the key residues of the SH3 domain (including F08, W36, P50, and Y53), denoted by DKC. Clearly, the SWCNT will prevent the regular binding of the ligand to the SH3 domain if it binds to one of the key residues. Thus, the reaction coordinate DKC could be a good reaction coordinate to describe the binding free energy landscape. To obtain a reasonable result, we sampled the structures with a time interval of 20 ps in the last 180 ns of all the 20 trajectories. The global minimum in this landscape is found at SPC ~ 220 Å2 and DKC < 5.0 Å (see Fig. 9(a)) with a binding free energy of -6.08 kcal/mol, which corresponds
to the state where the SWCNT binds to the SH3 domain binding pocket. Two local minima, in which SPC ~ 150 Å2 & DKC ~ 8.0 Å, and SPC ~ 180 Å2 & DKC ~ 18.0 Å, are also observed. Detail studies show that these two states correspond to the SWCNT binding to the side (the region around the residue W37) and the opposite (the region around the residue Y03) of the binding pocket of the SH3 domain, respectively.
The binding of the ligand to the SH3 domain is also described by the free energy landscape as a function of the two equivalent reaction coordinates: the interface area between the ligand and the SH3 domain (SPL), and the distance between the ligand and the key residues of the SH3 domain (DKC) (shown in Fig. 9(b)). The lowest binding affinity for the ligand is -5.51 kcal/mol, which is 0.57 kcal/mol less than the SWCNT shown above. This explains why the PRM ligand loses the competition to the SWCNT in the binding to the SH3 domain. Interestingly, the binding free energy landscape also shows that there are two basins, with the free energy of -5.51 kcal/mol (the lowest for the ligand binding) and -5.16 kcal/mol, respectively. Both states have similar values of SPL ~ 180 Å2, i.e., adsorbing onto the SH3 domain with similar contacting surface areas. The difference lies in the DKL, i.e., the distance to the key residues of the SH3 domain binding pocket. The deeper one has a very small value of DKL ~ 5.0 Å, which indicates that the ligand is around the binding pocket. The other basin has a larger value of DKL ~ 12.0 Å, indicating that the ligand is far away the binding pocket. We note that the same DKL (~12 Å) does not mean that there is only one binding site for the ligand in this case. A detailed study shows that the ligand binds to a broad region of the SH3 domain with residues 11-19, 23-25, and 43-47 involved.
π−π stacking interaction
From Fig. 9(a), there are three basins in the free energy landscape of the SWCNT binding on the SWCNT. When the SWCNT was in the global minimum (binding on the binding pocket), it contacted with at least one of the aromatic residues, including F08, W36 and Y53. The upper panel of Fig. 10 shows the typical binding mode of the SWCNT contacting with these three aromatic residues, in which their side chains of are almost parallel with the SWCNT. For the two local minima, the SWCNT contacted with the Y03 and W37, respectively (shown in the lower panel of Fig. 10). Overall, the SWCNT always contacts with aromatic residues when it binds onto the SH3 domain. In the SH3 domain native structure, there are a total of 47 residues exposed to environment (solvent, ligand, or SWCNT), but only ten of them are hydrophobic and five of them aromatic residues. This indicates the SWCNT is favored to contact with the hydrophobic and aromatic residues (more below).
In order to describe quantitatively the degree by which the SWCNT favors to contact with the hydrophobic and aromatic residues, we estimated the p-values of distribution of the hydrophobic and aromatic residues in contact with the SWCNT based on the χ2 hypothesis  (see Fig. 11). It was found that the SWCNT is remarkably favored in binding with the hydrophobic residues of the SH3 domain. The p-value for the SWCNT in contacting with hydrophobic residues (Ile, Val, Leu, Met, Ala, Trp, and Phe, following the hydrophobicity classifications of this NCBI website: http://www.ncbi.nlm.nih.gov/Class/Structure/aa/aa_explorer.cgi), aromatic residues (His, Tyr, Phe, and Trp), and the hydrophobic aromatic residues (Phe and Trp; here we do not include Tyr residue, even though many people might classify Tyr to be hydrophobic aromatic residue as well) were 6.51 x 10-3, 2.42 x 10-12, and 2.00 x 10-19, respectively. This indicates that the interactions between the SWCNT and the hydrophobic residues, particularly the aromatic residues (π−π stacking interactions), play an important role in the binding of SWCNT and proteins. We note that favorable interactions of the carbon nanotube with the hydrophobic and aromatic residues had also been observed in recent carbon nanotube-peptide experiments [33, 34, 89]. That is, these interactions are independent of the sizes of proteins and CNTs. Therefore, even though the SWCNT used in our simulation is quite small, the observation of the high probability to occupy the binding pocket of the protein by the SWCNT is extendable to the hydrophobic nanoscale particles of larger lengths, which is expected to be observed with more extensive computations with the development of the supercomputers as well as advanced simulation techniques. It should also be noted that the π−π stacking interactions might be underestimated in standard force fields due to the lack of polarization, and more advanced techniques such as quantum mechanics simulations might be needed to fully catch the effect; however, a recent study shows that classical simulations can be sometimes even better to capture this π−π stacking interactions than the quantum mechanics simulations due to the facts that the classical force fields are often directly fitting from experimental data, while the quantum mechanics simulations, on the other hand, might suffer from the limited size and boundary conditions .
Effect on the binding event
It is noteworthy that a previous MD simulation has shown that the ligand could recognize the binding pocket and form the native binding mode with the SH3 domain, in the absence of nanoscale particles . The long-range electrostatic interactions are found to play an important role in the ligand recognition of the binding pocket, through its two arginine residues on the C-terminal . The PRM ligand’s R7 and R8 residues have favorable long-range electrostatic interactions with SH3-domain’s acidic residues, D14, E16, D17, and E33 on RT and n-sCr loops, near one end of the binding pocket, which provides guidance to an accelerated ligand binding mode search. This electrostatically accelerated association of proteins was also previously observed in experiments [91, 92]. And studies had confirmed that the arginine residue is conserved in the PRM to this SH3 domain [78, 79]. The stability of the final binding, on the other hand, is driven by the structural match and hydrophobic interactions between the ligand (proline residues and V4) and the SH3 domain (F8, F10, W36, P50, P52, and Y53), in addition to the salt-bridges formed by R7/R8 (see Fig. 6(a)). In our current complex system with the presence of SWCNT, the SWCNT wins the competition over the ligand for the binding pocket even without any guidance from the long-range electrostatic interactions (since SWCNT is uncharged in our simulations), indicating a very strong hydrophobic interaction between the SWCNT and the SH3 domain. Interestingly, the ligand still manages to locate the area of the active site and bind with both the SWCNT and SH3 domain (though no longer in the binding pocket), which implies the PRM ligand still maintains some recognition capability due to the long-range electrostatic interactions through its two C-terminal arginine residues.
The structural match between the ligand and active site stabilized the native binding mode. However, it raised the demand for an extensive conformational space mapping even with the electrostatic acceleration. And the hydrophilic residues in the ligand (the two terminal arginines) often shield the hydrophobic ones in the water environment to “inhibit” their exposure to the active site. On the other hand, the hydrophobic interactions between the SWCNT and the hydrophobic residues of the binding pocket, particularly those aromatic ones, are strong and nonspecific. This makes the SWCNT more straightforward to be adsorbed onto the binding pocket. Thus it takes only about 35 ns for the SWCNT to reach a stable binding state at the active site, while it takes ~130 ns for the ligand to come near the SH3 domain in the previous MD simulations . This is also consistent with the above free energy landscape analysis, which reveals that the SWCNT is more favorable than the PRM ligand to be bound to the active site of the SH3 domain. That is, in the competition for the binding pocket of the SH3 domain, the SWCNT has advantages over the PRM ligand in both kinetics and thermodynamics. The SWCNT essentially occupies the binding pocket of the SH3 domain and interrupts its native binding with the PRM ligand.
2.3. Interaction with human blood serum proteins
All the above simulations were done with relatively short CNTs (~20 Å) and small proteins (SH3 domain has 57 amino acids and WW 26 amino acids) due to the enormous computational resources needed. In order to see whether longer SWCNTs and larger proteins might exhibit similar behaviors, we have recently performed simulations with IBM BlueGene supercomputers on much larger systems – the SWCNTs interacting with human blood serum proteins BSA (bovine serum albumin), γ-Ig (gamma-globulin), Tf (transferring), and BFg (fibrinogen), also in an attempt to compare with experiments directly. The total number of atoms in the solvated systems is up to 1.06 million. Figure 12 shows one example of the system setup with BFg dimmer interacting with SWCNT in solution. Since this is still an ongoing research, only a very brief description of the current results will be given. More simulation results and detailed comparisons with experiments will be reported elsewhere .
To elucidate the mechanism of protein adsorption by the SWCNT, up to 150 ns molecular dynamics simulations have been performed for each protein-SWCNT complex in explicit solvent. In our experiments, including AFM images, fluorescence spectroscopy, and SDS-PAGE, we found that these serum proteins display a competitive binding/adsorption with the SWCNTs, with BFg showing the largest adsorption capacity. In addition, the adsorption capacity displays a general trend of BFg > γ-Ig > Tf > BSA .
As mentioned above, the aromatic residues, tryptophan (Trp), tyrosine (Tyr) and phenylalanine (Phe), could produce π-π stacking interactions between aromatic rings of amino acids and six-member rings of the SWCNT. As a first analysis, the numbers of these three kinds of amino-acid residues on the surface of initial conformation for each protein molecule were counted. The approximate numbers of (Trp, Tyr, Phe) for BFg, γ-Ig, Tf and BSA were in sequence of (14, 47, 27); (2, 25, 10); (1, 13, 15); and (1, 9, 10); respectively. This implies that the aromatic residues on the surface of these blood serum proteins particularly consist of Tyr and Phe, and lack of Trp, except for BFg which has 14 surface Trp residues. Most of the Trp residues are wrapped in the core of proteins, and play a lesser role in these serum protein-SWCNT complex adsorptions, which seems to be different from the above WW domain and SH3 domain cases where Trp residues play a significant role.
The snapshots of MD simulation results describe the preferred binding sites on proteins for SWCNT. If we highlight the residues within 5 Å distance from the surface of SWCNT, the aromatic residues Tyr and Phe stand out as the key residues involved in the adsorption. Interestingly, only Tyr and Phe were observed in the adsorption on the surface of SWCNT (at least in the first 150 ns simulations, longer runs might reveal larger protein conformational changes and exposure of Trp residues to the SWCNT surface as well), and Trp residues made little contribution for adsorption in the present MD simulations, which are still very short in time scale as compared to experiments (5 min to 1 hour). Nevertheless, our simulations show that there is a positive correlation between the total numbers of Trp, Tyr, and Phe residues and the binding surface areas, thus the adsorption content, of the proteins, respectively.
In order to illuminate the above positive correlation relationship more clearly, we calculated the contact residue number and contact surface area of various proteins, in adsorption of SWCNT versus time. We used the average contact residue number (ACRN) of the protein to quantify its adsorption ability with SWCNT. The ACRN for BFg during 140 to 150 ns is 94. The ACRN for γ-Ig after its saturation during 70 to 150 ns is 28. The ACRN for Tf after its saturation during 80 to 150 ns is 13. The ACRN for BSA after its saturation during 65 to 150 ns is 12. The relationship of ACRN of various protein complex systems shows an order of BFg > γ-Ig > Tf > BSA, which is in good agreement with the experimental findings  (and also compatible with the above simple analysis using aromatic hydrophobic residues). Similarly, we computed the average contact surface area (ACSA) for these blood serum proteins (see Fig. 13). The ACSA for BFg, γ-Ig, Tf and BSA during the responding time ranges are 4360 Å2, 1240 Å 2, 580 Å 2 and 470 Å 2 respectively, which again shows an order of BFg > γ-Ig > Tf > BSA, in good agreement with the experimental findings .
The most abundant Tyr and Phe residues on the surface of BFg molecule can be attributed to the highest protein quantities binding to SWCNT. Interestingly, the π-cation interactions are also found to play a role, though smaller, in the BFg-SWCNT binding. The average Arg and Lys residues in contact with SWCNT (ACRN-Arg, ACRN-Lys) are found to be about ~1.3 and ~3.5, respectively (out of the total ~94 ACRN above). The other three proteins γ-Ig, Tf, and BSA show similar results, with the (ACRN-Arg, ACRN-Lys) having values of (0.0, 0.2), (1.4, 3.2), (0.1, 2.4), respectively. These findings are consistent with a previous report by Kagan and coauthors, where they found that the interactions between human myeloperoxidase (hMPO) and the carboxylated nanotubes involved both arginine (π–π-cation interaction) and tyrosine (π–π interaction) residues . Our current simulations show that the π–π stacking interactions are the driving forces for the competitive binding of human serum proteins onto SWCNTs.
Finally, it should be noted that carbon nanotubes often exist in much longer lengths in in vitro experiments as compared to our simulation lengths, which are typically a few to a few tens of nanometers. However, there are recent studies showing that the enzyme horseradish peroxidase (HRP) can biodegrade SWCNTs into very small fragments . Therefore short SWCNTs may exist in the living body and studying SWCNTs with relatively short lengths might be of significant importance to human health as well.
We had investigated the effect of the CNTs on the conformation and function of proteins by using the large scale molecular dynamics simulations. WW domain, SH3 domains were used as examples to illustrate our ideas that the CNTs can affect the function of the proteins by either disrupting the structure of active site or shielding the active site from competing ligands. In the simulation of the SWCNT and the WW domain, it was found that the SWCNT can plug into the hydrophobic core of WW domain to form a stable protein-SWCNT complex. This insertion of the SWCNT broke the scaffold which used to bind the proline rich motifs (PRM), and thus reduced the possibility of the direct binding between the PRM and the WW domain. In the case of SH3 domain, we studied a three-way binding competition among the SWCNT, the PRM ligand, and the SH3 domain, and found that the SWCNT had a very high probability occupying the binding pocket of the ligand in the SH3 domain, with about 0.6 kcal/mol more favorable binding affinity than the original PRM ligand. Moreover, in most of the simulation trajectories, the ligand will be adsorbed onto the SWCNT, which may further reduce the possibility of the correct binding for the ligand onto the active site of the SH3 domain. These adverse effects of the SWCNTs on proteins, including both the disruption of the active sites and the competitive binding with the ligands, might seriously damage the original functions of proteins, suggesting the potential nanotoxicity of the SWCNT. The interactions between the SWCNTs and the hydrophobic residues, particularly the aromatic residues (π−π stacking interactions), are found to play a key role in both the insertion of the SWCNT and the competitive binding with ligands. Moreover, we also simulated much larger systems, with human serum proteins interacting with SWCNTs, and found a competitive binding among all these serum proteins, with the adsorption capacity in an order of BFg > γ-Ig > Tf > BSA, in good agreement with experimental findings. Again, the π–π stacking interactions are found to be the driving forces for the competitive binding of human serum proteins onto SWCNTs. These findings might have shed some light to the general mechanism of the interactions between hydrophobic nanoparticles and proteins.
4. Computational method
Modeling WW domain and SWCNT
The protein YAP65 WW domain was prepared from the Protein Data Bank (PDB code: 1JMQ, truncated to include residues 15-40), and modeled by AMBER03 force field . The lengths of all those SWCNTs are 19.54 Å. The carbon atoms of the SWCNTs were modeled as uncharged Lennard–Jones particles with a cross-section of σcc = 3.40 Å, and a depth of the potential well of εcc= 0.36 kJ/mol [56, 57]. The interactions between these carbon atoms of SWCNT and other atoms were generated by the AMBER03 force field . The combined systems were then solvated in a rhombic dodecahedral periodic box with the distance between the solutes and box boundary at least 8 Å. The numbers of water molecules were 2694, 2709 and 2735 for the system with the SWCNT of m = 4, 5, and 6, respectively. And a Cl- was added into solution to neutralize the system.
Modeling SH3 domain, Ligand and SWCNT
The SH3 domain and the ligand were prepared from the Protein Data Bank (PDB code: 1CKB , residues 134-190 for the SH3 domain, and residue 1-8 for the ligand), and modeled by AMBER03 force field . The carbon atoms of the SWCNT were also modeled as uncharged Lennard–Jones particles with a cross-section of σcc = 3.40 Å, and a depth of the potential well of εcc = 0.36 kJ/mol [56, 57]. Two three-way (SWCNT + ligand + protein) binding complex systems were set up. In the first system, the initial separations of the geometric centre of the SWCNT and the ligand from that of the SH3 domain were both set as 30 Å, a distance long enough to avoid the starting point bias. The initial orientations of the SWCNT and the RPM ligand versus the SH3 domain were set at different directions. The resulting complex was then solvated in a rhombic dodecahedral periodic box, with the distance between the solutes and the boundary of the box at least 10 Å. The final complex system size is 12,378 atoms. For the second complex system, a similar procedure was followed, but with the SWCNT and ligand swapped in space (thus in the opposite orientations when compared to the first system) for effective sampling. The final system size is 13,367 atoms. The TIP3P  water model was used for the salvation, and three Na+ were added into solution to neutralize each system.
WW domain and SH3 domain MD simulations
The MD simulations were performed by using the Gromacs package 4.0 . In the simulations, the covalent bonds involving H atoms were constrained by the LINCS algorithm , allowed a time step of 2 fs. The long range electrostatic interactions were treated with the particle-mesh Ewald method (PME) with a grid spacing of 1.2 Å . The cutoff for the van der Waals interaction was set to 10 Å. After energy minimization, all the systems were equilibrated by MD simulations for 200 ps at a constant pressure of 1 bar and temperature of 298 K using Berendsen coupling . Then the production simulations were performed in the NVT ensemble at 298 K.
Human serum proteins interacting with SWCNT
The human blood serum proteins simulated in this work were taken from the Protein Data Bank and Swiss-Model Repository as BFg with PDB ID: 1DEQ, γ-Ig with PDB ID: 3HR5, Tf with PDB ID: 2HAV, and BSA with SWISS-MODEL ID: 432779d395a52bfc9f6574bc3e98afcd_1. All SWCNTs used in the molecular dynamics (MD) simulations for the adsorption of CNT for proteins BFg, Ig, Tf, and BSA were armchair CNT (14, 14) with diameter as around 2.0 nm, which were consistent with the experimental ones. The geometrical coordinate parameters of CNT were generated by using Nanotube Modeler software (www.jcrystal.com/products/wincnt/).
MD simulations were performed by NAMD  of version 2.6. The various blood proteins were modeled by CHARMM 32b1 force field . Carbon atoms of SWCNTs were assumed to be type CA with Lennard-Jones parameters kcal/mol, Å. The interaction between these carbon atoms of SWCNTs and other atoms were generated by CHARMM 32b1 force field. The various blood proteins were well-separated from SWCNT with the minimum distances up to 40.0 Å initially. All systems were solvated in periodic TIP3P modeled water box, and neutralized by adding sodium and chloride ions with 0.2 M salt concentration as physiological condition, with system sizes up to 1,058,598 atoms (see Fig. 12 for one example of BFg).
Simulations were performed at constant temperature of 310 K, and pressure of 1 atm, with the PME method for long-range electrostatic interactions and time step of 2 fs. The cutoff for the Van der Waals interaction was set to 12 Å. All systems were simulated in the NPT ensemble for more than 150 ns. Visual Molecular Dynamics (VMD 1.8.7.) graphics viewer software was utilized to illustrate the adsorption of proteins and SWCNT in different representations.
We thank Dr. Peng Xiu, Dr. Chunlei Wang, Dr. Payel Das, Dr. Seung-gu Kang, and Prof. Bruce Berne for helpful discussions. This research is supported in part by grants from NNSFC (10825520), NBRPC (973 Program: 2007CB936000 and 2007CB814800), Shanghai Leading Academic Discipline Project (B111), and Shanghai Supercomputer Center of China. RZ acknowledges the support from the IBM BlueGene Science Program.