Molecular modeling of interaction between lipid monolayer and graphene nanosheets: implications for pulmonary nanotoxicity and pulmonary drug delivery

Abstract

Understanding how nanoparticles interact with the pulmonary surfactant monolayer (PSM) is of great importance for safe applications in biomedicine and for evaluation of both health and environment impacts. Here, by performing molecular dynamics simulations, we propose a possible origin of the pulmonary nanotoxicity of graphene-based nanoparticles that comes from a rigidifying effect of graphene nanosheets (GNs) on PSM. This, in reality, indicates that once captured by the PSM, inhaled GNs are hard to be removed from the PSM partially because the expiration or PSM compression is locally restrained, possibly leading to GN accumulation on the PSM. The local rigidifying effect, which is enhanced as multiple GNs approach each other, is found to be dependent on the GN hydrophobicity. In the expiration or PSM compression process, the hydrophilic GN keeps adhering to the monolayer–air interface, while the hydrophobic GN tends to be hosted in the hydrophobic interior and internalize into the PSM via self-rotation. Besides the spontaneous internalization via PSM compression, our pulling simulations indicate that both pulmonary internalization and externalization of GNs can be accomplished by direct translocation across the PSM. The effect of GN hydrophobicity on the direct PSM translocation is well supported by the free energy analysis. This work will help our understanding of pulmonary nanotoxicity of GNs and provide useful guidelines for molecular design of GN-based pulmonary drug delivery materials.

---

1 Introduction

In recent years, bionanotechnology has opened new avenues for both clinical diagnostics and therapeutics.¹ In particular, nanoparticles (NPs) with different physicochemical properties have been designed for drug delivery and bio-sensing.^2–4 Nevertheless, accumulating evidence suggests that NPs may cause adverse health effects during or after their internalization.^5,6 Therefore, understanding how NPs interact with different bio-interfaces is of quite importance for their safe applications in broad biomedical area.⁷

In the last decade, interactions of NPs with lipid membranes have been intensively studied both experimentally^8–17 and theoretically.^18–46 In general, both pathway and mechanism of NP–membrane interaction are affected by the NP properties, like NP size,^11,18,19,42 NP shape,^{8,20,21,24,25,28,31,32,35} surface charge,¹⁰ hydrophobicity,^13,36 and ligand arrangement,^16,23,37 but also membrane properties, like membrane tension,⁴² membrane asymmetry⁴⁵ and transmembrane potential.³⁸

Besides traditional internalization pathways, the human lung is easily accessible for NPs due to its specific properties, like thin epithelial barrier, large surface area, abundant underlying vasculature and low acidity.⁴⁷ The alveolar epithelial cells, where the gas exchange takes place, are covered by a thin pulmonary surfactant monolayer (PSM) exposed to the air. The PSM is composed of approximately 90% lipids, with dipalmitoylphosphatidylcholine (DPPC) being the most abundant single lipid component for most mammalian species, and 10% proteins. The main physiological function of PSM is to ‘switch off’ the tension of the air–water interface. On a molecular level, the reduction of surface tension is mediated by reversible morphology transformations from monolayer to bilayer and multilayer. Both pathways and molecular mechanisms of the PSM morphology transformation have been intensively investigated by both molecular dynamics simulations and experiments.^48–51

During respiration, external NPs can be captured by the PSM and disrupt its function, thus leading to enhanced susceptibility to pathological conditions, such as pneumonia and asthma. As the PSM is the first barrier against NP internalization, the specific NP–PSM interaction directly determines both the extent of pulmonary nanotoxicity and subsequent fate of NPs. On the positive side, the study of NP–PSM interactions will inform novel guidelines for the design of targeted pulmonary drug delivery materials with both high efficiency and low toxicity. Experimentally, adverse biophysical effects of NPs on natural PSM have been investigated by Zuo and his coworkers.^14,52,53 Some other research groups have performed coarse-grained molecular dynamics simulations to explore the interaction mechanism between spherical NPs and model PSM.^39,40 Nevertheless, the question of how to balance the pulmonary nanotoxicity and pulmonary drug delivery is still open. Therefore, it is in urgent need to further study the pathway and mechanism by which the inhaled NPs interact with the PSM.

In this work, our focus is mainly on the interaction of model PSM with graphene nanosheets (GNs). The choice of GN is based on its extraordinary electrical, mechanical, and thermal properties,^54,55 and its great potential for various bio-applications.^{23,33,56–60} In general, the GNs are hydrophobic in nature while hydrophilic GNs are usually obtained via surface functionalization. At present, GNs have been developed in many different forms in terms of shapes, sizes, hydrophobicities, chemical modifications and other characteristics that can produce dramatically different results when studied biologically. Therefore, the detailed interaction between PSM and GNs with different properties is of equal importance for both evaluating their pulmonary nanotoxicity and developing the pulmonary drug delivery materials. Here, the PSM is simply represented by a pure DPPC monolayer. Both mechanical perturbation of PSM by GN adsorption and pulmonary internalization/externalization mechanisms will be discussed.

It should be noted that the PSM is a mixture, with DPPC the major component but there are several other relevant components, like POPG, cholesterol and surfactant proteins (including SP-A, SP-B, SP-C and SP-D). Although the DPPC is the major component of PSM responsible for the surface tension reduction, the effect of other components is equally important. For example, the addition of POPG has a fluidizing effect on PSM and thus assists the morphology transformation.⁶¹ As an important component of PSM, cholesterol was found to modulate the phase behavior of PSM. Accordingly, the local mechanical property and the resultant morphology transformation will be changed.⁶² Among the surfactant proteins, SP-A and SP-D are large hydrophilic glycoproteins that are members of the Ca²⁺-dependent carbohydrate-binding collectin family. They can only weakly associate with the PSM but play a primary role in pulmonary host defense.⁶³ Comparatively, both SP-B and SP-C are hydrophobic and thus participate directly in the morphology transformation process of the PSM.⁶⁴ Therefore, it will be much more complex when we take all the PSM components into account to study the GN–PSM interaction. Rather than that, we choose to first understand the interaction between GNs and pure DPPC monolayers. The effect of each other component will be discussed individually in our future work.

2 Model and methods

2.1 Coarse-grained models

Compared with atomistic simulation models, the coarse-grained (CG) models, which map a group of atoms into one interaction site, allow simulations to be carried out with a larger length scale and longer time scale.^65–67 Specifically, one of the most widely used CG force field is the Martini force field.⁶⁵ It has been widely used to simulate biomolecules, such as lipids,⁶⁸ biopolymers⁶⁹ and proteins,⁷⁰ but also to simulate NPs.⁷¹ The demonstrated ability of the Martini force field to reproduce biological phenomena makes it an attractive tool to perform CG simulations of GN–PSM interactions. The Martini force field is based on a four-to-one mapping, i.e., on average four heavy atoms are represented by a single interaction site. Only four main types of CG particles, including polar (P), nonpolar (N), apolar (C), and charged (Q), are considered. Each particle type has a number of subtypes, which allow for a more accurate representation of chemical properties of the underlying atomic structure.

The CG DPPC and water force field parameters were downloaded in details from http://md.chem.rug.nl/cgmartini/. The honeycomb structure of GN is reduced to a triangular lattice of CG beads, where every three carbons in the all-atom GN are modeled as one particle. The angular force constants and equilibrium angles are given by K_angle = 700 kcal mol⁻¹ rad⁻² and θ₀ = 60°, respectively. The dihedral angle force constants, multiplicity and potential minima are defined as K_d = 3.1 kcal mol⁻¹, n = 2, and θ = 180°, respectively. The specific GN CG model comes from the previous simulation work of Titov et al., and has been successfully used to simulate the GN-membrane interaction.⁷² Besides, both C1 and Nda particle types are used to represent a GN particle with hydrophobic and hydrophilic properties, respectively (C1 and Nda are parameters in Martini force field).⁴⁰ The GN edge length, l_g is fixed to 7.7 nm and two different initial orientations, including parallel and perpendicular with respect to the PSM, are compared in the simulations.

2.2 Simulation method

As shown in Fig. 1, the simulation setup includes a water slab bounded by air on both sides with two symmetric monolayers at each air–water interface. In our system, the air space is represented by vacuum and the PSM is simply represented by a pure DPPC monolayer, which is composed of 2 × 4096 DPPC molecules. The number of water beads is fixed to 383 744. It should be noted that the similar bi-monolayer model has been used by other research groups to investigate both structural and dynamic properties of PSM.^39,40,48,49

Fig. 1 Schematic representation of the simulation setup. (a) bi-monolayer PSM model, (b) DPPC lipid model, (c) hydrophilic (purple) and hydrophobic (yellow) GN model, with the edge length of 7.7 nm.

For all simulations, a cutoff of 1.2 nm was used for van der Waals interactions, and the Lennard-Jones potential was smoothly shifted to zero between 0.9 nm and 1.2 nm to reduce the cutoff noise. For electrostatic interactions, the coulombic potential, with a cutoff of 1.2 nm, was smoothly shifted to zero from 0 to 1.2 nm. Lipids, water, and GN were coupled separately to Berendsen heat baths at T = 310 K, with a coupling constant τ = 1 ps. The monolayer compression was simulated using semi-isotropic pressure coupling. The system compressibility was set to 5 × 10⁻⁵ bar⁻¹ in the lateral direction and zero in the normal direction. The latter ensures constant size of the box normal to the membrane. The neighbor list for non-bonded interactions was updated every 10 steps. In order to accomplish the translocation of GNs across the monolayer, an external spring force is exerted on the center of mass (COM) of the GN. The pulling rate and force constant are fixed to 0.0002 nm ps⁻¹ and 1000 kJ mol⁻¹ nm⁻², respectively. Snapshots of the simulation system in this paper were rendered by VMD.⁷³ All simulations were performed using GROMACS 4.5.3.⁷⁴

3 Results and discussion

3.1 Rigidifying effect of GNs on the PSM

First, to understand the mechanical perturbation of PSM by GN adsorption, we calculate the average surface tension under different lateral pressures (Fig. 2). Four different systems are considered here. They are the pure DPPC monolayer, DPPC monolayer with one hydrophilic GN, DPPC monolayer with one hydrophobic GN, and DPPC monolayer with four hydrophilic GNs. In order to improve the statistics, 5 independent simulations are performed under each lateral pressure. Each simulation is performed for 80 ns and the surface tension is averaged for the last 40 ns. In order to verify the reliability of the shorter simulation time, we choose two typical examples and extend the simulations to 1.0 μs.

Fig. 2 The average surface tension as a function of lateral pressure. Four different systems are considered, including the pure DPPC monolayer (black), DPPC monolayer with one hydrophilic GN (red), DPPC monolayer with one hydrophobic GN (blue), and DPPC monolayer with four hydrophilic GNs (cyan).

For pure DPPC monolayer, we gradually increase the lateral pressure from P = −5 bar to as large as P = 5 bar. We find that the surface tension first undergoes a linear decrease as the increase of lateral pressure (−5 bar < P < 3 bar). As the lateral pressure exceeds a critical value (P = 3 bar), the monolayer surface tension keeps nearly unchanged. Accordingly, a morphology transformation from monolayer to bilayer occurs (Fig. S1†).^48,49 Next, we placed a GN, either hydrophilic or hydrophobic, on the monolayer–air interface and performed simulations under different lateral pressures. As shown in Fig. 2, the biophysical function of DPPC monolayer is in some extent restrained by the GN adsorption, that is, in the presence of GN a larger lateral pressure is required to induce the monolayer collapse compared with the situation in the absence of GN. Comparatively, the effect of hydrophobic GN is less obvious than that of hydrophilic GN, especially under the lateral pressure of monolayer buckling (P ≥ 2 bar). Moreover, when we place four hydrophilic GNs on the PSM, the surface tension keeps decreasing even when the lateral pressure is as large as 5 bar. This implies that the GN–PSM interaction is not only dependent on the GN hydrophobicity but also affected by the GN concentration.

It should be noted that the effect of GN adsorption on surface tension of PSM is less distinguishable under negative lateral pressure than that under positive lateral pressure (Fig. 2). We interpret that the limited effect is on one hand because of the much smaller GN size (49 nm²) compared with that of the PSM (1950 nm² under zero surface tension). On the other hand, the effect of GN adsorption on PSM is mainly reflected by the increased local bending rigidity. Under negative lateral pressure, the PSM tends to undergo expansion while no deformation is induced. Accordingly, the effect of increased bending rigidity on the tension reduction of PSM is negligible. Under positive lateral pressure, contrarily, the PSM has the tendency to undergo the buckling even collapse to release the surface energy. As the GN adsorption can locally rigidify the PSM and thus restrain the morphology transformation, the effect on surface tension reduction is accordingly enhanced. Recently, Barnoud et al. performed MD simulations to show that the lipid monolayer collapse is promoted by C₆₀ fullerene.⁷⁵ They interpret that the promoted monolayer collapse is ascribed to the lower bending rigidity of monolayer in the presence of C₆₀. The opposite effect of NP adsorption on the monolayer collapse indicates that the detailed interaction between NPs and lipid monolayer is strongly affected by the NP properties, like size, shape, and hydrophobicity.

Next, we monitored the interaction between DPPC monolayer and GN with different hydrophobicity and different concentration. According to Fig. 2, the effect of GN adsorption on the tension reduction is much more distinguishable under positive lateral pressure. Thus, we first fixed the lateral pressure to P = 3 bar, under which the monolayer–bilayer transformation occurs in the absence of GN (Fig. 2 and S1†). When a hydrophilic GN is adsorbed on the monolayer surface, however, the monolayer–bilayer transformation is strongly restrained (Fig. 3). In this case the hydrophilic GN keeps adhering on the monolayer surface and thus rigidifies the local region of the monolayer where it covers. The rigidifying effect is clearly illustrated by the typical snapshots, in which two obvious monolayer deformations occurs far from the GN, while the local monolayer beneath and adjacent to the GN keeps flat during the simulation (Fig. 3). As a result, the subsequent monolayer–bilayer transformation is restrained by the enhanced bending rigidity.

Fig. 3 Restraining effect of hydrophilic GN on morphology transformation of DPPC monolayer. The lateral pressure is fixed to P = 3 bar, under which a monolayer–bilayer transformation occurs in the absence of GN (Fig. S1†). Two membrane deformations are labeled by white arrows.

It should be noted that the time scale of microseconds is not uncommon for the simulations of lipid monolayer. For example, Baoukina et al. recently performed 10–25 μs simulations to investigate the collapse of heterogeneous lipid monolayer.⁷⁶ The time scale of microseconds is mainly because the lipid demixing of heterogeneous monolayer generally needs much longer simulation time. Comparatively, the simulations of morphology transformation of lipid monolayers without lipid demixing are much cheaper.^48,49 Here, the PSM is simply modeled by the pure DPPC monolayer, thus the shorter time scale is suitable for the simulation of GN–PSM interaction. Nevertheless, we believe that longer simulation time can make the system to relax and thus further minimize the free energy. Therefore, we extend the above simulation to 1.0 μs, during which the hydrophilic GN keeps adhering on the PSM surface and no obvious morphology transformation is observed (Fig. S2†).

To verify the rigidifying effect of hydrophilic GN on PSM, we next increase the local concentration of GNs adsorbed on the PSM. To this end, four hydrophilic GNs are placed on the monolayer surface. As in the case having one GN, four hydrophilic GNs keep adhering on the PSM surface and no monolayer–bilayer transformation is observed at P = 3 bar (Fig. 4A). As we increase the lateral pressure to P = 5 bar, the surface tension keeps unchanged (Fig. 2) and accordingly the monolayer–bilayer transformation occurs in the presence of one GN. When four close-located GNs adsorb on the PSM, however, no morphology transformation is observed (Fig. 4C). We then extend the simulation to 1.0 μs and find that the four hydrophilic GNs keep adhering on the monolayer surface and thus restrain the subsequent morphology transformation (Fig. S3†). This indicates that the local rigidifying effect of GNs on PSM is dose dependent and is further enhanced as multiple GNs approach each other.

Fig. 4 Rigidifying effect of 4 hydrophilic GNs on local DPPC monolayer. (A–C) show typical snapshots under different pressure (A: 3 bar; B: 4 bar; C: 5 bar), (D) shows time evolution of area per lipid under lateral pressure of P = 3 bar in the presence of 1 and 4 GNs; (E) shows time evolution of area per lipid under lateral pressure of P = 5 bar in the presence of 1 and 4 GNs.

Similar with the NP concentration, we speculate that the rigidifying effect of GN on PSM is also affected by the GN size. Here, we construct three GNs with different sizes, including L = 5.5 nm, L = 11 nm, and L = 16.5 nm, respectively. After a number of simulations under different pressures, we calculate the surface tension of PSM as functions of average lipid area (Fig. S4A†). As a result, we can determine the area compressibility modulus K_A from the linear relation,⁷⁷

where represents the mechanical tension of the PSM adsorbed with GN, A is the lipid area and A₀ is the lipid area under zero mechanical tension.

According to Fig. S4A,† the surface tension of PSM increases linearly but with different gradient as functions of the average lipid area. We then calculate the area compressibility modulus to K_A = 291.3 mN m⁻¹, 356.3 mN m⁻¹, and 439.8 mN m⁻¹ for PSMs adsorbed with GNs with the size of L = 5.5 nm, 11 nm, and 16.5 nm, respectively (Fig. S4B†). It should be noted that the area compressibility modulus for pure lipid monolayer was determined to K_A = 290 mN m⁻¹.⁴⁸ Therefore, it is evident that the rigidifying effect of GN on PSM can be enhanced by increasing the GN size.

Compared with the GN size and local GN concentration, increasing the GN thickness has the similar but different effect. Noting that the production of multilayered graphene is much easier than that of monolayered graphene, we next explore the effect of GN thickness on the interaction with PSM. To this end, we construct three multilayered GNs with the same edge length (L = 7.7 nm) but different thicknesses (N_layer = 2, 4, 6). The lateral pressure is fixed to P = 3 bar during the simulations. Our simulations suggest that increasing the GN thickness has two aspects of effect. One is to increase the inner bending rigidity of GN. Similar with that of GN concentration, the PSM compression and deformation is further restrained (Fig. S5†). The other effect is to reduce the anisotropy of GN. According to the typical snapshots, GN with larger thickness (N_layer = 6) can rotate more easily on the PSM surface (Fig. S5C†). Accordingly, the PSM undergoes a severe local deformation to wrap the GN. As the interaction of PSM with hydrophilic GN is less favorable than that with hydrophobic GN, we speculate that the effect of GN thickness may be more significant for hydrophobic GNs.

Now we turn to the detailed interaction pathway for the case of hydrophilic GN at P = 5 bar. As the monolayer is locally rigidified, the monolayer buckling has to occur far from the site where the GN locates (Fig. 5A, t = 10 ns, white arrow). Contrarily, the local monolayer in the vicinity of GN is relatively flat until the monolayer collapse occurs. As the simulation proceeds, one of the buckles develops into a collapse and the monolayer–bilayer transformation occurs (Fig. 5A, t = 55 ns, 65 ns). It is noteworthy that the hydrophilic GN keeps adhering on the monolayer–air interface during the compression and transformation process. The evolution of angle between GN and membrane normal is shown in Fig. 5B, which further illustrates the interaction pathway of hydrophilic GN with the monolayer. In order to energetically explain the interaction mechanism, we monitored the interaction energy between GN and DPPC molecules (Fig. 5C). After a drastic decrease, the interaction energy keeps fluctuating around the value of −3700 kJ mol⁻¹.

Fig. 5 Monolayer buckling and collapse in the presence of hydrophilic GN. (A) Shows the typical snapshots, (B) shows the time evolution of the angle between GN and membrane normal, and (C) shows the time evolution of the interaction energy between GN and DPPC molecules. The lateral pressure is fixed to P = 5 bar.

To understand the interaction of hydrophobic GN with DPPC monolayer, we fixed the lateral pressure to P = 5 bar, under which the monolayer–bilayer transformation occurs in the presence of one hydrophilic GN (Fig. 6). According to the typical snapshots and time evolution of angle between GN and membrane normal, the hydrophobic GN first adheres on the monolayer surface and undergoes a slight tilt during the monolayer buckling process (Fig. 6). Similar to the case of hydrophilic GN, the adsorption of hydrophobic GN rigidifies the local monolayer. As a result, the monolayer buckling tends to occur far from the GN. As the simulation proceeds, the local monolayer–bilayer transformation occurs at one of the buckling sites. Meanwhile, the hydrophobic GN is gradually encapsulated in the hydrophobic interior as the collapse proceeds. The specific embedding process is accompanied by an obvious GN rotation, which is illustrated in both typical snapshots (Fig. 6A, t = 60 ns, 70 ns) and time evolution of the angle between GN and membrane normal (Fig. 6B, t = 60 ns). Moreover, the GN embedding and GN rotation are accompanied by a drastic decrease of the GN–DPPC interaction energy (Fig. 6C). This is because the hosting of GN in the hydrophobic interior maximizes the contact area between hydrophobic GN and hydrophobic lipid tails. Therefore, the specific PSM internalization of hydrophobic GN via the monolayer–bilayer transformation is energetically favorable.

Fig. 6 Monolayer buckling and collapse in the presence of hydrophobic GN. (A) Shows the typical snapshots, (B) shows the time evolution of the angle between GN and membrane normal, and (C) shows the time evolution of the interaction energy between GN and DPPC molecules. The lateral pressure is fixed to P = 5 bar.

Here we should note that the cut-off is used to treat the electrostatic interaction. The choice of cut-off based electrostatics is based on the fact that the GN is electroneutral and thus no electrostatic interaction exists between GN with PSM. Moreover, many other simulation works have used the cut-off based electrostatic interaction to investigate the dynamic properties of pure lipid monolayer.^48,49,78 We thus speculate that the effect of using cut-off electrostatic interaction on the GN–PSM interaction can be very limited. Nevertheless, it has been shown in another work that the PAMAM dendrimer can induce pore in DMPC membrane when the PME was used to treat the electrostatic interaction. However, the membrane pore was not observed when authors used the cut-off based electrostatics.⁷⁹ Therefore, it is important to explore how much the results are influenced by use of cut-off based electrostatics. To this end, we use the PME to treat the electrostatic interaction and repeat the above simulation of interaction between hydrophobic GN and PSM under lateral pressure of P = 5 bar. According to the typical snapshots, the monolayer–bilayer transformation occurs at about t = 58 ns, after which the hydrophobic GN is wrapped by the deformed monolayer (Fig. S6†). This is similar with that using the cut-off based electrostatics interaction. Therefore, it can be speculated that the influence of simulation results by using cut-off based electrostatic interaction is negligible in this work.

In fact, recent experiments showed that hydrophobic NPs have minor influence on the lateral pressure-area isotherm, but do alter the morphology of the condensed domains in the DPPC monolayer.⁸⁰ This is consistent with our present simulations, in which the restraining effect of hydrophilic GN on PSM compression is more obvious than that of hydrophobic GN. Moreover, we explain that the minor effect for hydrophobic GNs is ascribed to the detailed interaction pathway. Specifically, the hydrophobic GNs tend to be wrapped by the deformed PSM and thus internalize into the hydrophobic interior. This accordingly reduces the restraining effect on the PSM compression. Comparatively, the hydrophilic GNs keep adhering on the PSM surface during the compression process. As a result, the subsequent morphology transformation can be effectively restrained.

3.2 Pulmonary internalization of GNs via PSM translocation

Due to the extraordinary electrical, mechanical, and thermal properties,^54,55 GN-based materials have become prevalent in a variety of biomedical applications.^56,57 For example, the large surface area of GN allows high-density biofunctionalization or drug loading.⁸¹ In the above section, we showed that hydrophilic GN stays adhering on the monolayer–air interface, while the hydrophobic GN can be hosted in the hydrophobic interior during the expiration or PSM compression process. This implies a potential pulmonary internalization pathway for hydrophobic GNs. In the inspiration or PSM expansion process, however, the unique pathway for GN internalization is the direct translocation across the monolayer. Therefore, we next fixed the lateral pressure to P = −15 bar and performed Pulling MD simulations.

It should be noted that the pulling simulations of PSM translocation of GNs are based on the fact that the inhaled NPs have nonzero speed and thus first collide with the PSM after their inhalation. We speculate that the PSM translocation can be accomplished once the membrane resistance is overcome. In reality, the triggering force needed for PSM translocation of NPs is derived from the increased NP translocation speed, which is generally determined by the inhalation speed. It was measured that the average inspiratory reserve volume is about 3.0 litres for adults. Besides, the average human respiratory rate is 12–20 breaths per minute. Therefore, it normally takes 1 second to get about 500–1000 cm³ of air into the lung. As the volume flow rate is Q = velocity × area, the inhalation velocity is thus Q/area. Considering the terminal bronchiole is about 0.5 mm in diameter and the average number is about 60 000, the total cross area of the terminal bronchiole is about 0.012 m². Therefore, the inhalation velocity is estimated to about 0.04–0.08 m s⁻¹.

In order to accomplish the GN translocation, an external spring force is exerted on the center of mass (COM) of the GN. Both the COM pulling energy and the angle between GN and membrane normal is monitored to thermodynamically and dynamically elucidate the translocation process. First, the translocation of hydrophobic GN is divided into two stages (Fig. 7 and Video S1†). In the first stage, the GN keeps adhering on the monolayer surface and no translocation occurs. Meanwhile, the monolayer undergoes a severe caving due to pushing of the GN (Fig. 7A, 96 ns). In the second stage, the GN translocation is initiated by piercing of one corner across the monolayer (Fig. 7A, 108 ns). Subsequently, the translocation is completed via a spontaneous GN rotation (Fig. 7A and B). After translocation completes, a number of DPPC molecules is extracted out of the monolayer and adsorb on the GN surface (Fig. 7A, 120 ns). In fact, graphene was shown to have antibacterial activity on Escherichia coli.^82,83 Recent MD simulations revealed that graphene can penetrate into and extract large amounts of lipids from the cell membranes.³³ The destructive extraction offers a novel mechanism for molecular basis of graphene's cytotoxicity. This is similar with our present case, in which a large number of DPPC molecules are extracted by the PSM translocation of GNs. Besides, it has been suggested experimentally that large structural disruptions of DPPC monolayer can make some DPPC molecules detach from the air–water interface, which may further affect PS homeostasis.⁸⁴ Therefore, we speculate that the structure and biofunction can be disrupted in some extent by the penetration of hydrophobic GNs.

Fig. 7 Pulmonary internalization of hydrophobic GN via translocation across the monolayer. The GN is placed initially in the air phase and is parallel with the monolayer. (A) Shows the typical snapshots, (B) shows the time evolution of angle between GN and membrane normal, and (C) shows the time evolution of COM pulling energy.

It is noteworthy that the GN self-rotation is followed by a sudden decrease of COM pulling energy. This energetically highlights the importance of GN rotation in its translocation process (Fig. 7C). In fact, the promoting effect of spontaneous nanoparticle rotation on cellular internalization is generous and has been proved by other simulation works.^{21,23–25,85,86} Nevertheless, the high COM pulling energy indicates that the pulmonary internalization of hydrophobic GNs via translocation across the PSM is rather difficult, at least when they are initially parallel with the monolayer surface. As the high COM pulling energy is mainly ascribed to the severe monolayer deformation by GN pushing, we next set the GN keep standing during the translocation process. Indeed, both typical snapshots and evolution of COM pulling energy indicate that the translocation of GN across the monolayer is highly facilitated by decreasing the contact area between GN and monolayer (Fig. S7†).

We note that the translocation speed is fixed to 0.0002 nm ps⁻¹ (0.2 m s⁻¹), which is much higher than the normal inhalation speed (0.04–0.08 m s⁻¹). Considering that in some cases one need to intentionally increase or decrease the inhalation speed, it is significant to explore the effect of GN speed on the PSM translocation. To this end, we decrease the pulling rate to 0.00002 nm ps⁻¹ to repeat the PSM translocation of hydrophobic GN. According to the typical snapshots and time evolution of COM pulling energy, the extent of membrane resistance is strikingly reduced by decreasing the translocation speed (Fig. S8†). In more detail, the PSM first undergoes a deformation to resist the GN translocation. Subsequently, the PSM translocation is initiated by flipping of the GN (Fig. S8A,† t = 430 ns). We note that at higher translocation speed (0.0002 nm ps⁻¹), a number of lipid molecules are extracted from the monolayer and adsorb on the GN surface (Fig. 7). Comparatively at lower speed, the GN flipping and translocation is accompanied by the local morphology transformation of the PSM. This is because the decreased translocation speed can leave more time for the PSM to relax and thus the GN is encapsulated in the hydrophobic interior and forms the GN-sandwiched structure. Finally, the GN is detached from the PSM and the COM pulling energy further decreases to zero (Fig. S8B†).

Comparatively, the translocation of hydrophilic GNs across the PSM is different from that of hydrophobic ones. First, when the hydrophilic GN is parallel with the monolayer, the translocation is rather difficult and is accompanied by the monolayer rupture (Fig. S9†). This is ascribed to the large contact area and unfavorable interaction between hydrophilic GN and hydrophobic monolayer surface. However, when the hydrophilic GN is initially perpendicular with the monolayer, the translocation is highly improved (Fig. S10†). It should be noted that a number of lipids are extracted from the monolayer by the hydrophilic GN. This is mainly due to the attraction between hydrophilic bead of lipids and hydrophilic GN. Nevertheless, we analyze that the lipid extraction is energetically unfavorable because the hydrophobic tail of lipid is exposed to the water (Fig. S10,† t = 80 ns). Therefore, we speculate that the extracted lipid may return to the monolayer if we decrease the translocation speed. Besides the initial orientation of GNs, the PSM translocation of hydrophilic GNs is also affected by the translocation speed. As we decrease the translocation speed from 0.0002 nm ps⁻¹ to 0.00002 nm ps⁻¹, no monolayer rupture is observed during the translocation process. Besides, the translocation of hydrophilic GN is initiated by piercing of the corner across the PSM, which is similar with that of hydrophobic GNs (Fig. S11†). Moreover, the strikingly decreased COM pulling energy indicates that the translocation of hydrophilic GN can also be accomplished without PSM damage by decreasing the translocation speed.

3.3 Pulmonary externalization of GNs via PSM translocation

After pulmonary internalization, one may ask whether or how the internalized GNs are exhausted from the lung. To answer this question, we placed one hydrophobic GN in the water phase and performed pulling simulations. According to the typical snapshots, the GN first adheres on the monolayer–water interface (Fig. 8 and Video S2†). The upward pushing force exerted on the monolayer results in a severe monolayer deformation (Fig. 8). Accordingly, the COM pulling energy keeps increasing until the translocation starts (Fig. 8B). At about 35 ns, the GN translocation is initiated by a transient GN flipping (Fig. 8A and C), which is accompanied by a sudden decrease of the COM pulling energy (Fig. 8B). After translocation, the exhausted GN keeps adhering on the hydrophobic monolayer surface and keeps pulling the monolayer upwards. As a result, the COM pulling energy continues to increase after translocation. Comparatively, the translocation is much easier when the GN is initially perpendicular with the monolayer (Fig. S12 and Video S3†). After the tip-first translocation (Fig. S12A,† 12 ns), the GN undergoes a rotation (Fig. S12A,† 14–40 ns) and keeps adhering on the hydrophobic monolayer surface (Fig. S12A,† 80 ns). Therefore, once internalized into the PSM, the reverse translocation of hydrophobic GNs across the monolayer is easy but the detachment from the monolayer surface is rather difficult.

Fig. 8 Pulmonary externalization of hydrophobic GN by translocation across the monolayer. The GN is placed initially in the water phase and is parallel with the monolayer. (A) Shows the typical snapshots, (B) shows the time evolution of COM pulling energy, and (C) shows the time evolution of the angle between GN and membrane normal.

For hydrophilic GNs, the pulmonary externalization via PSM translocation is relatively easier, at least when the GN is initially perpendicular with the monolayer. According to the typical snapshots and time evolution of COM pulling energy, the translocation process is stepwise (Fig. S13 and Video S4†). In the first stage, the COM pulling energy keeps increasing until the whole monolayer thickness is crossed by the hydrophilic GN (Fig. S13A,† 30 ns; Fig. S9B†). The increase of COM pulling energy is ascribed to the unfavorable interaction between hydrophilic GN and hydrophobic lipid tails. After that, the translocation is much easier as the COM pulling energy starts to decrease. It is noteworthy that no obvious membrane deformation is observed during the simulation. However, when the GN is initially parallel with the monolayer, the translocation is strongly restrained by a severe membrane deformation (Fig. S14A†). Moreover, an obvious monolayer pore is generated by the large upward pushing force (Fig. S14B†).

3.4 Free energy analysis

To understand the effect GN hydrophobicity on the PSM translocation, we then calculate the potential of mean force (PMF) to provide an estimate of the free energy change as the GN penetrates across the PSM. The detailed procedure is as follows. First, the harmonic potential with a force constant of 1000 kJ mol⁻¹ nm⁻² was applied to the GN. Consequently, the GN is pulled to penetrate across the PSM in a constant rate (0.0001 nm ps⁻¹). For simplicity, both hydrophobic and hydrophilic GNs are constrained to keep parallel with the PSM normal during the translocation process. After the pulling simulation completes, a number of frames are extracted from the trajectory that correspond to the desired COM of the GN. Then the umbrella sampling simulation was performed on each configuration to restrain it within a window corresponding to the chosen COM position. Finally, we used the Weighted Histogram Analysis Method (WHAM) to extract the PMF from the biased distribution.⁸⁷

For the hydrophilic GN, the PSM translocation can be divided into three stages. The first stage is featured by the short increase in the free energy (Fig. 9A). This short membrane resistance is ascribed to the unfavorable interaction between hydrophilic GN and hydrophobic lipid tails (0 → 1). Once the whole monolayer thickness is crossed by the GN, the PMF shows a long decrease in the free energy (Fig. 9A). In this stage, the favorable interaction between hydrophilic GN and water dominates the PSM translocation. In addition, the monolayer undergoes an upward local deformation in order to maximize the contact number between GN and water molecules (2). After the whole GN is submerged in the water phase, the local monolayer deformation gradually disappears and the PSM slowly decreases to zero (3).

Fig. 9 Potential of mean force (PMF) for the PSM translocation of both hydrophilic (a) and hydrophobic (b and c) GNs. Two different initial orientations for hydrophobic GNs are considered. One is perpendicular with the PSM (b) and the other is parallel with the PSM (c).

For the hydrophobic GN, the corresponding PMF profile first shows a short decrease in the free energy and reaches a local minimum as the GN enters the lipid monolayer (Fig. 9B, 0 → 1). The decreasing PMF indicates that the interaction between hydrophobic GN and hydrophobic lipid tails is favorable and thus promotes the insertion process. However, after the whole monolayer thickness is crossed (Z < 12.5 nm), the PSM turns to increase as more and more GN beads get contacted with the hydrophilic lipid heads and water molecules. In order to protect the hydrophobic GN from water, a number of lipid molecules are extracted from the monolayer and get attached on the GN surface. Meanwhile, the lipid monolayer undergoes a severe damage induced by the GN insertion (2). This is energetically unfavorable and thus the PMF undergoes a continuous increase until more than half of the GN is below the monolayer–water interface (1 → 2). As the GN further moves downwards, more GN surface is wrapped by the lipid and the monolayer damage is gradually repaired. Consequently, the PMF profile shows a long drastic decrease in the free energy (2 → 3). This indicates the preferential location of GN below the lipid monolayer and completed wrapped by the lipid molecules (3). However, the PMF shows a drastic increase of the free energy as the wrapped GN further moves downwards. We analyze that the increased membrane resistance is mainly ascribed to the severe monolayer deformation induced by the GN (4). After the neck connecting the GN and monolayer breaks (5), the extent of membrane resistance is effectively reduced.

We have to note that the perpendicular orientation of GNs is less statistically likely in reality due to the extremely small contact area with the PSM. Therefore, it promotes us to further calculate the free energy change of PSM translocation of GNs which are initially parallel with the PSM. According to Fig. 9C, the PMF keeps increasing during the whole translocation process. This clearly indicates that PSM translocation of GN is very energetically unfavorable when the GN is parallel with the PSM. The typical snapshots suggest that the high energy barrier comes from the severe membrane deformation (Fig. 9C). Nevertheless, the PMF increases with a lower gradient once the PSM translocation is initiated by the GN flipping (Fig. 9C, 3). This is quite similar with our previous pulling simulations, in which the pulling energy drops after the GN flipping takes place (Fig. 7 and S8†).

3.5 Implications to pulmonary nanotoxicity and pulmonary drug delivery

It has been known that PSM reduces the surface tension of the air–water interface by reversible morphology transformation and thus can accommodate the large size change of alveoli. During respiration, either harmful or beneficial NPs can be inhaled. As the PSM forms the first alveolar capillary barrier against inhaled NPs, both examination of pulmonary nanotoxicity and development of pulmonary drug delivery should consider the NP–PSM interaction. However, we must admit that it is quite difficult to systematically probe and visualize the interaction under various conditions, due to available experimental technologies. Computer simulations, on the contrary, may offer the underlying mechanism because it allows the investigation on a molecular level. Besides, it can provide distinct insights from the microscopic to mesoscopic view into the interactions in well-controlled conditions. Recently, Ding et al. have presented the recent theoretical and computational progress on the NPs–membrane interaction under different biological processes.⁸⁸ Besides, the links between simulations and experiments are also connected and well discussed. Here, both equilibrium and non-equilibrium MD simulations are performed to investigate the interaction mechanism between GNs and DPPC monolayer. We think our simulation results have implications for both understanding of the pulmonary nanotoxicity and designing of NP-based pulmonary drug delivery.

First, our equilibrium simulations suggest that the inhaled GNs have adverse effects on the biophysical function of PSM. The pulmonary nanotoxicity of GNs is in some extent derived from the rigidifying effect of GN adsorption on the PSM. As one of the results, the respiration will be strongly hampered due to the restrained PSM compression. Although the rigidifying effect of NPs on the PSM has been previously reported experimentally,¹⁴ we speculate that the rigidifying effect is more obvious for GNs because of the two-dimensional shape. On the other hand, once captured by the PSM, the GNs especially hydrophobic GNs are difficult to be removed from the PSM. This may lead to their accumulation on the PSM surface and thus further intensify the adverse effect. In fact, several animal studies have verified the pulmonary accumulation of NPs by determining their quantity in total lung homogenate preparations following their ingestion.^89–91 Consequently, the likelihood of free NPs being respired is reduced and thus the pulmonary nanotoxicity will be further amplified.⁹²

Here, we should note that the rigidifying effect of GN on PSM is strongly affected by the inner bending rigidity of GN. For example, GNs with higher bending rigidity usually have stronger rigidifying effect on PSM, while flexible GNs are expected to undergo deformations during the PSM buckling process. Experiments on bulk graphite give ≈1 Tpa for the in-plane Young's modulus,⁹³ while Young's modulus in graphene varies between 0.5 and 1.0 TPa.^94,95 The GN CG model used in this work comes from the previous simulations by Titov et al. By adjusting the bond and angular constant, they obtained the Young's modulus to 928 Gpa, which is similar with the experiment result.⁹⁵ It also should be noted that the bending rigidity of graphenes can be tuned via a number of strategies, like introducing defects or modulating the functional groups.⁹⁶ Therefore, exploring how bending rigidity of GNs affects the interaction with PSM will help us to design drug delivery materials with both high efficiency and low toxicity.

The pulmonary nanotoxicity of GNs is also reflected by the destructive extraction of lipids from the PSM. Our non-equilibrium simulations suggest that the PSM translocation of hydrophobic GNs is generally accompanied by the lipid extraction. This is speculated because the lipid extraction from the air–water interface can induce the structural destruction of PSM and finally affects the PSM homeostasis.⁸⁴ Besides, GNs have been shown to have antibacterial activity on Escherichia coli.^82,83 Recent MD simulations explained that the GNs can penetrate into cell membranes and extract a number of lipids from the bilayer.³³ This offers a novel mechanism for molecular basis of GN's cytotoxicity. Therefore, we suggest that the intentional inhalation of NPs should avoid the lipid extraction from the air–water interface.

Nowadays, the concern about potential diagnostic and therapeutic applications of NPs in respiratory and systemic diseases has been growing.^97–99 As the PSM forms the first barrier of pulmonary internalization of NPs, the detailed NP–PSM interaction will not only determines the extent of pulmonary nanotoxicity but also affects the subsequent fate of inhaled NPs. Here we distinguish two possible PSM internalization mechanisms of GNs. On one hand, the PSM internalization of GNs can be a result of the morphology transformation. In this scenario, the GN is wrapped by the deformed PSM and delivered to the hydrophobic interior of PSM. Our simulations suggest that this pathway is feasible only for the hydrophobic GNs (Fig. 6), while hydrophilic GNs keep adhering on the PSM surface during PSM compression process (Fig. 5). The other scenario is that the internalization is due to direct translocation across the PSM, although the structure of PSM can be disrupted more or less. In reality, the inhaled NPs have nonzero speed, the value of which is mainly determined by the inhalation speed. Thus, the inhaled NPs first collide with the PSM and penetrate across it once the barrier is overcome.

4 Conclusions

In this paper, we have presented a systematical investigation on interactions of graphene nanosheets (GNs) with pulmonary surfactant monolayer (PSM). First, our simulations suggest that the pulmonary nanotoxicity of GNs can in some extent be reflected by rigidifying locally a region of the monolayer larger than the domain area they covered. This, in reality, indicates that once a number of GNs are captured by the PSM, the expiration or PSM compression will be strongly restrained. Comparatively, the restraining effect is more obvious for GNs with hydrophilic property. In particular, hydrophilic GNs stay adhering on the monolayer–air interface during the expiration process, while the hydrophobic GNs can be encapsulated in the hydrophobic interior and form the GN-sandwiched structure. The specific interaction pathway implies a potential pulmonary internalization pathway for the hydrophobic GNs.

In the inspiration or PSM expansion process, the unique pulmonary internalization pathway for GNs is the direct translocation across the PSM. Our pulling simulations suggest that the pulmonary internalization of GNs is strongly dependent on both the GN hydrophobicity and the initial orientation with respect to the monolayer. In general, the translocation of GNs with an initial orientation perpendicular with the monolayer is much more favorable than that parallel with the monolayer. After pulmonary internalization, the externalization via translocation is quite different for hydrophobic and hydrophilic GNs. Moreover, after translocation, the detachment of hydrophobic GNs from the monolayer surface is comparatively difficult, while for hydrophilic GNs the detachment is almost spontaneous. The present investigation of interaction between GN and model PSM should provide fundamental information needed to help understand pulmonary nanotoxicity and imply potential applications in pulmonary drug delivery.

Acknowledgements

This work was financially supported by National Natural Science foundation of China (no. 21303269, 21103230, and 21273287) and the Natural Science Foundation for Distinguished Young Scholar of Shandong Province (no. JQ201008). T. YUE is grateful for the support of Natural Science foundation of Shandong Province (no. ZR2013BQ029) and the Qingdao Science and Technology Project (no. 13-1-4-235-jch). The authors thank the National Super-computing Center in Jinan and Shenzhen for providing computer time.

References

1. T. L. Doane and C. Burda, The unique role of nanoparticles in nanomedicine: imaging, drug delivery and therapy, Chem. Soc. Rev., 2012, 41, 2885–2911.
2. D. Bechet, P. Couleaud, C. Frochot, M. L. Viriot, F. Guillemin and M. Barberi-Heyob, Nanoparticles as vehicles for delivery of photodynamic therapy agents, Trends Biotechnol., 2008, 26, 612–621.
3. Y. W. C. Cao, R. C. Jin and C. A. Mirkin, Nanoparticles with roman spectroscopic fingerprints for DNA and RNA detection, Science, 2002, 297, 1536–1540.
4. T. Wei, J. Liu, H. Ma, Q. Cheng, Y. Huang, J. Zhao, S. Huo, X. Xue, Z. Liang and X. J. Liang, Functionalized nanoscale micelles improve the drug delivery for cancer in vitro and in vivo, Nano Lett., 2013, 13, 2528–2534.
5. P. R. Leroueil, S. Hong, A. Mecke, J. R. Baker, B. G. Orr and M. M. B. Holl, Nanoparticle interaction with biological membranes: does nanotechnology present a janus face?, Acc. Chem. Res., 2007, 40, 335–342.
6. C. J. Murphy, A. M. Gole, J. W. Stone, P. N. Sisco, A. M. Alkilany, E. C. Goldsmith and S. C. Baxter, Gold nanoparticles in biology: beyond toxicity to cellular imaging, Acc. Chem. Res., 2008, 41, 1721–1730.
7. A. E. Nel, L. Mädler, D. Velegol, T. Xia, E. M. V. Hoek, P. Somasundaran, F. Klaessig, V. Castranova and M. Thompson, Understanding biophysicochemical interactions at the nano-bio interface, Nat. Mater., 2009, 8, 543–557.
8. S. E. A. Gratton, P. A. Ropp, P. D. Pohlhaus, J. C. Luft, V. J. Madden, M. E. Napier and J. M. DeSimone, The effect of particle design on cellular internalization pathways, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 11613–11618.
9. Y. Roiter, M. Ornatska, A. R. Rammohan, J. Balakrishnan, D. R. Heine and S. Minko, Interaction of nanoparticles with lipid membrane, Nano Lett., 2008, 8, 941–944.
10. E. C. Cho, J. Xie, P. A. Wurm and Y. Xia, Understanding the role of surface charges in cellular adsorption versus internalization by selectively removing gold nanoparticles on the cell surface with a I2/KI etchant, Nano Lett., 2009, 3, 1080–1084.
11. W. Jing, B. Y. S. Kim, J. T. Rutka and W. C. W. Chan, Nanoparticle-mediated cellular response is size-dependent, Nat. Nanotechnol., 2008, 3, 145–150.
12. A. Albanese and W. C. W. Chan, Effect of gold nanoparticle aggregation on cell uptake and toxicity, ACS Nano, 2011, 5, 5478–5489.
13. A. Verma and F. Stellacci, Effect of surface properties on nanoparticle–cell interactions, Small, 2010, 6, 12–21.
14. Q. Fan, Y. E. Wang, X. Zhao, J. S. C. Loo and Y. Y. Zuo, Adverse biophysical effects of hydroxyapatite nanoparticles on natural pulmonary surfactant, ACS Nano, 2011, 5, 6410–6416.
15. A. K. Sachan, R. K. Harishchandra, C. Bantz, M. Maskos, R. Reichelt and H. Galla, High-resolution investigation of nanoparticle interaction with a model pulmonary surfactant monolayer, ACS Nano, 2012, 6, 1677–1687.
16. A. Verma, O. Uzun, Y. Hu, Y. Hu, H. Han, N. Watson, S. Chen, D. J. Irvine and F. Stellacci, Surface-structure-regulated cell-membrane penetration by monolayer-protected nanoparticles, Nat. Mater., 2008, 7, 588–595.
17. R. Frost, G. E. Jönsson, D. Chakarov, S. Svedhem and B. Kasemo, Nano Lett., 2012, 12, 3356–3362.
18. S. Zhang, J. Li, G. Lykotrafitis, G. Bao and S. Suresh, Size-dependent endocytosis of nanoparticles, Adv. Mater., 2009, 21, 419–424.
19. H. Gao, W. Shi and L. B. Freund, Mechanics of receptor-mediated endocytosis, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 9469–9474.
20. X. Shi, A. Bussche, R. H. Hurt, A. B. Kane and H. Gao, Cell entry of one-dimensional nanomaterials occurs by tip recognition and rotation, Nat. Nanotechnol., 2011, 6, 714–719.
21. Y. Li, H. Yuan, A. Bussche, M. Creighton, R. H. Hurt, A. B. Kane and H. Gao, Graphene microsheets enter cells through spontaneous membrane penetration at edge asperities and corner sites, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 12295–12300.
22. X. Yi, X. Shi and H. Gao, Cellular uptake of elastic nanoparticles, Phys. Rev. Lett., 2011, 107, 098101.
23. Y. Li, X. Li, Z. Li and H. Gao, Surface-structure-regulated penetration of nanoparticles across a cell membrane, Nanoscale, 2012, 4, 3768–3775.
24. K. Yang, B. Yuan and Y. Ma, Influence of geometric nanoparticle's rotation on cellular internalization process, Nanoscale, 2013, 5, 7998–8006.
25. K. Yang and Y. Ma, Computer simulation of the translocation of nanoparticles with different shapes across a lipid bilayer, Nat. Nanotechnol., 2010, 5, 579–583.
26. A. Alexeev, W. E. Uspal and A. C. Balazs, Harnessing janus nanoparticles to create controllable pores in membranes, ACS Nano, 2008, 2, 1117–1122.
27. M. Dutt, O. Kuksenok, M. J. Nayhouse, S. R. Little and A. C. Balazs, Modeling the self-assembly of lipids and nanotubes in solution: forming vesicles and bicelles with trans-membrane nanotube channels, ACS Nano, 2011, 5, 4769–4782.
28. R. Vácha, F. J. Martinez-Veracoechea and D. Frenkel, Receptor-mediated endocytosis of nanoparticles of various shapes, Nano Lett., 2011, 11, 5391–5395.
29. B. J. Reynwar, G. Illya, V. A. Harmandaris, M. M. Müller, K. Kremer and M. Deserno, Aggregation and vesiculation of membrane proteins by curvature mediated interactions, Nature, 2007, 447, 461–464.
30. M. Deserno and W. M. Gelbart, Adhesion and wrapping in colloid-vesicle complexes, J. Phys. Chem. B, 2002, 106, 5543–5552.
31. C. Huang, Y. Zhang, H. Yuan, H. Gao and S. Zhang, Role of nanoparticle geometry in endocytosis: laying down to stand up, Nano Lett., 2013, 13, 4546–4550.
32. S. Dasgupta, T. Auth and G. Gompper, Shape and orientation matter for the cellular uptake of nonspherical particles, Nano Lett., 2014, 14, 687–693.
33. Y. Tu, M. Lv, P. Xiu, T. Huynh, M. Zhang, M. Castelli, Z. Liu, Q. Huang, C. Fan, H. Fang and R. Zhou, Destructive extraction of phospholipids from Escherichia coli membranes by graphene nanosheets, Nat. Nanotechnol., 2013, 8, 594–601.
34. J. Wong-Ekkabut, S. Baoukina, W. Triampo, I. Tang, D. P. Tieleman and L. Monticelli, Computer simulation study of fullerene translocation through lipid membranes, Nat. Nanotechnol., 2008, 3, 363–368.
35. X. Xie, A. Xu, M. R. Angle, N. Tayebi, P. Verma and N. A. Melosh, Mechanical model of vertical nanowire cell penetration, Nano Lett., 2013, 13, 6002–6008.
36. H. M. Ding and Y. Q. Ma, Designing nanoparticle translocation through membranes by computer simulations, ACS Nano, 2012, 6, 1230–1238.
37. H. M. Ding and Y. Q. Ma, Role of physicochemical properties of coating ligands in receptor-mediated endocytosis of nanoparticles, Biomaterials, 2012, 33, 5798.
38. J. Lin and A. Alexander-Katz, Cell membranes open “doors” for cationic nanoparticles/biomolecules: insights into uptake kinetics, ACS Nano, 2013, 7, 10799–10808.
39. G. Hu, B. Jiao, X. Shi, R. P. Valle, Q. Fan and Y. Y. Zuo, Physicochemical properties of nanoparticles regulate translocation across pulmonary surfactant monolayer and formation of lipoprotein corona, ACS Nano, 2013, 7, 10525–10533.
40. X. Lin, T. Bai, Y. Y. Zuo and N. Gu, Promote potential applications of nanoparticles as respiratory drug carrier: insights from molecular dynamics simulations, Nanoscale, 2014, 6, 2759–2767.
41. N. Arai, K. Yasuoka and X. C. Zeng, A vesicle cell under collision with a janus or homogeneous nanoparticle: translocation dynamics and late-stage morphology, Nanoscale, 2013, 5, 9089–9100.
42. R. Lipowsky and H.-G. Döbereiner, Vesicles in contact with nanoparticles and colloids, Europhys. Lett., 1998, 43, 219–225.
43. M. Raatz, R. Lipowsky and T. R. Weikl, Cooperative wrapping of nanoparticles by membrane tubes, Soft Matter, 2014, 10, 3570–3577.
44. J. C. Shillcock and R. Lipowsky, Tension-induced fusion of bilayer membranes and vesicles, Nat. Mater., 2005, 4, 225–228.
45. T. Yue, X. Zhang and F. Huang, Membrane monolayer protrusion mediates a new nanoparticle wrapping pathway, Soft Matter, 2014, 10, 2024–2034.
46. A. D. Phan, T. X. Hoang, T. Phan and L. M. Woods, Repulsive interactions of a lipid membrane with graphene in composite materials, J. Chem. Phys., 2013, 139, 184703.
47. H. S. Choi, Y. Ashitate, J. H. Lee, S. H. Kim, A. Matsui, N. Insin, M. G. Bawendi, M. Semmler-Behnke, J. V. Frangioni and A. Tsuda, Rapid translocation of nanoparticles from the lung airspaces to the body, Nat. Biotechnol., 2010, 28, 1300–1303.
48. S. Baoukina, L. Monticelli, H. J. Risselada, S. J. Marrink and D. P. Tieleman, The molecular mechanism of lipid monolayer collapse, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 10803–10808.
49. S. Baoukina, L. Monticelli, M. Amrein and D. P. Tieleman, The molecular mechanism of monolayer–bilayer transformations of lung surfactant from molecular dynamics simulations, Biophys. J., 2007, 93, 3775–3782.
50. M. M. Lipp, K. Y. C. Lee, D. Y. Takamoto, J. A. Zasadzinski and A. J. Warning, Coexistence of buckled and flat monolayers, Phys. Rev. Lett., 1998, 81, 1650–1653.
51. E. Hatta, Sequential collapse transitions in a Langmuir monolayer, Langmuir, 2004, 20, 4059–4063.
52. Y. E. Wang, H. Zhang, Q. Fan, C. R. Neal and Y. Y. Zuo, Biophysical interaction between corticosteroids and natural surfactant preparation: implications to pulmonary drug delivery using surfactant as a carrier, Soft Matter, 2012, 8, 504–511.
53. R. P. Valle, C. L. Huang, J. S. C. Loo and Y. Y. Zuo, Increasing hydrophobicity of nanoparticles intensifies lung surfactant film inhibition and particle retention, ACS Sustainable Chem. Eng., 2014, 2, 1574–1580.
54. C. Lee, X. Wei, J. W. Kysar and J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene, Science, 2008, 321, 385–388.
55. K. S. Novoselov, Nobel Lecture: graphene: materials in the flatland, Rev. Mod. Phys., 2011, 83, 837–849.
56. P. K. Ang, M. Jaiswal, C. H. Y. X. Lim, Y. Wang, J. Sankaran, A. Li, C. T. Lim, T. Wohland, O. Barbaros and K. P. Loh, A bioelectronic platform using a graphene-lipid bilayer interface, ACS Nano, 2010, 4, 7387–7394.
57. C. Wang, J. Li, C. Amatore, Y. Chen, H. Jiang and X. Wang, Gold nanoclusters and graphene nanocomposites for drug delivery and imaging of cancer cells, Angew. Chem., Int. Ed., 2011, 50, 11644–11648.
58. K. Yang, L. Feng, X. Shi and Z. Liu, Nano-graphene in biomedicine: theranostic applications, Chem. Soc. Rev., 2013, 42, 530–547.
59. L. Yan, F. Zhao, S. Li, Z. Hu and Y. Zhao, Low-toxic and safe nanomaterials by surface-chemical design, carbon nanotubes, fullerenes, metallofullerenes, and graphenes, Nanoscale, 2011, 3, 362–382.
60. R. Guo, J. Mao and L. Yan, Computer simulation of cell entry of graphene nanosheet, Biomaterials, 2013, 34, 4296–4301.
61. R. Veldhuizen, K. Nag, S. Orgeig and F. Possmayer, The role of lipids in pulmonary surfactant, Biochim. Biophys. Acta, 1998, 1408, 90–108.
62. S. Baoukina, E. Mendez-Villuendas and D. P. Tieleman, Molecular view of phase coexistence in lipid monolayer, J. Am. Chem. Soc., 2012, 134, 17543–17553.
63. F. X. McCormack and J. A. Whitsett, The pulmonary collectins, SP-A and SP-D, Orchestrate innate immunity in the lung, J. Clin. Invest., 2002, 109, 707–712.
64. J. Pérez-Gil, Structure of pulmonary surfactant membranes and films: the role of proteins and lipid-protein interactions, Biochim. Biophys. Acta, 2008, 1778, 1676–1695.
65. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. de Vries, The MARTINI force field: coarse grained model for biomolecular simulations, J. Phys. Chem. B, 2007, 111, 7812–7824.
66. L. Y. Lu and G. A. Voth, Systematic coarse-graining of a multicomponent lipid bilayer, J. Phys. Chem. B, 2009, 113, 1501–1510.
67. Z. J. Wang and M. Deserno, A systematically coarse-grained solvent-free model for quantitative phospholipid bilayer simulations, J. Phys. Chem. B, 2010, 114, 11207–11220.
68. S. Baoukina, S. J. Marrink and D. P. Tieleman, Molecular structure of membrane tethers, Biophys. J., 2012, 102, 1866–1871.
69. G. Rossi, L. Monticellli, S. R. Puisto, I. Vattulainen and T. Ala-Nissila, Coarse-graining polymers with the MARTINI force-field: polystyrene as a benchmark case, Soft Matter, 2011, 7, 698–708.
70. M. Baaden and S. J. Marrink, Coarse-grain modeling of protein-protein interactions, Curr. Opin. Struct. Biol., 2013, 23, 878–886.
71. S. Baoukina, L. Monticelli and D. P. Tieleman, Interaction of pristine and functionalized carbon nanotubes with lipid membranes, J. Phys. Chem. B, 2013, 117, 12113–12123.
72. A. V. Titov, P. Král and R. Pearson, Sandwiched graphene–membrane superstructures, ACS Nano, 2010, 4, 229–234.
73. W. Humphrey, A. Dalke and K. Schulten, VMD: visual molecular dynamics, J. Mol. Graphics, 1996, 14, 33–38.
74. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation, J. Chem. Theory Comput., 2008, 4, 435–447.
75. J. Barnoud, L. Urbini and L. Monticelli, C60 fullerene promotes lung monolayer collapse, J. R. Soc., Interface, 2015, 12, 20140931.
76. S. Baoukina, D. Rozmanov, E. Mendez-Villuendas and D. P. Tieleman, The mechanism of collapse of heterogeneous lipid monolayers, Biophys. J., 2014, 107, 1136–1145.
77. B. Rózycki and R. Lipowsky, Spontaneous curvature of bilayer membranes from molecular simulations: Asymmetric lipid densities and asymmetric adsorption, J. Chem. Phys., 2015, 142, 054101.
78. S. Baoukina, L. Monticelli, S. J. Marrink and D. P. Tieleman, Pressure-area isotherm of a lipid monolayer from molecular dynamics simulations, Langmuir, 2007, 23, 12617–12623.
79. H. Lee and R. G. Larson, Coarse-grained molecular dynamics studies of the concentration and size dependence of fifth- and seventh-generation PAMAM dendrimers on pore formation in DMPC bilayer, J. Phys. Chem. B, 2008, 112, 7778–7784.
80. S. Tatur and A. Badia, Influence of hydrophobic alkylated gold nanoparticles on the phase behavior of monolayers of DPPC and clinical lung surfactant, Langmuir, 2012, 28, 628–639.
81. L. Zhang, J. Xia, Q. Zhao, L. Liu and Z. Zhang, Functional graphene oxide as a nanocarrier for controlled loading and targeted delivery of mixed anticancer drugs, Small, 2010, 6, 537–544.
82. W. Hu, C. Peng, W. Luo, M. Lv, X. Li, D. Li, Q. Huang and C. Fan, Graphene-based antibacterial paper, ACS Nano, 2010, 4, 4317–4323.
83. K. Krishnamoorthy, M. Veerapandian, L. Zhang, K. Yun and S. J. Kim, Antibacterial efficiency of graphene nanosheets against pathogenic bacteria via lipid peroxidation, J. Phys. Chem. C, 2012, 116, 17280–17287.
84. J. Perez-Gil and T. E. Weaver, Pulmonary surfactant pathophysiology: current models and open questions, Physiology, 2010, 25, 132–141.
85. Y. Li, T. Yue, K. Yang and X. Zhang, Molecular modeling of the relationship between nanoparticle shape anisotropy and endocytosis kinetics, Biomaterials, 2012, 33, 4965–4973.
86. T. Yue, X. Wang, F. Huang and X. Zhang, An unusual pathway for the membrane wrapping of rodlike nanoparticles and the orientation- and membrane wrapping-dependent nanoparticle interaction, Nanoscale, 2013, 5, 9888–9896.
87. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method, J. Comput. Chem., 1992, 13, 1011–1021.
88. H. M. Ding and Y.-Q. Ma, Theoretical and computational investigations of nanoparticle–biomembrane interactions in cellular delivery, Small, 2014, 11, 1055–1071.
89. J. J. Li, S. Muralikrishnan, C. T. Ng, L. Y. Yung and B. H. Bay, Nanoparticle-induced pulmonary toxicity, Exp. Biol. Med., 2010, 235, 1025–1033.
90. E. Fabian, R. Landsiedel, L. Ma-Hock, K. Wiench, W. Wohlleben and B. van Ravenzwaay, Tissue distribution and toxicity of intravenously administered titanium dioxide nanoparticles in rats, Arch. Toxicol., 2008, 82, 151–157.
91. J. Wang, G. Zhou, C. Chen, H. Yu, T. Wang, Y. Ma, G. Jia, Y. Gao, B. Li, J. Sun, Y. Li, F. Jiao, Y. Zhao and Z. Chai, Acute toxicity and biodistribution of different sized titanium dioxide particles in mice after oral administration, Toxicol. Lett., 2007, 168, 176–185.
92. S. T. Stern and S. E. McNeil, Nanotechnology safety concerns revisited, Toxicol. Sci., 2008, 101, 4–21.
93. O. L. Blakslee, D. G. Proctor, E. J. Seldin, G. B. Spence and T. Weng, Elastic constants of compression-annealed pyrolytic graphite, J. Appl. Phys., 1970, 41, 3373–3382.
94. I. W. Frank, D. M. Tanenbaum, A. M. van der Zande and P. L. McEuen, Mechanical properties of suspended graphene sheets, J. Vac. Sci. Technol., B: Microelectron. Nanometer Struct.--Process., Meas., Phenom., 2007, 25, 2558–2561.
95. C. Lee, X. Wei, J. W. Kysar and J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene, Science, 2008, 321, 385–388.
96. H. C. Schniepp, K. N. Kudin, J.-L. Li, R. K. Prud'homme, R. Car, D. A. Saville and I. A. Aksay, Bending properties of single functionalized graphene sheets probed by atomic force microscopy, ACS Nano, 2008, 2, 2577–2584.
97. J. C. Sung, B. L. Pulliam and D. A. Edwards, Nanoparticles for drug delivery to the lungs, Trends Biotechnol., 2007, 25, 563–570.
98. C. Y. Dombu and D. Betbeder, Airway delivery of peptides and proteins using nanoparticles, Biomaterials, 2013, 34, 516–525.
99. D. B. Buxton, Nanotechnology in the diagnosis and management of heart, lung and blood diseases, Expert Rev. Mol. Diagn., 2007, 7, 149–160.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra04922a

---