Role of surface ligands in nanoparticle permeation through a model membrane: a coarse-grained molecular dynamics simulations study

How nanoparticles interact with biological membranes is of significant importance in determining the toxicity of nanoparticles as well as their potential applications in phototherapy, imaging and gene/drug delivery. It has been shown that such interactions are often determined by nanoparticle physicochemical factors such as size, shape, hydrophobicity and surface charge density. Surface modification of the nanoparticle offers the possibility of creating site-specific carriers for both drug delivery and diagnostic purposes. In this work, we use coarse-grained molecular dynamic simulations to explore the permeation characteristics of ligand-coated nanoparticles through a model membrane. We compare permeation behaviors of ligand-coated nanoparticles with bare nanoparticles to provide insights into how the ligands affect the permeation process. A series of simulations is carried out to validate a coarse-grained model for nanoparticles and a lipid membrane system. The minimum driving force for nanoparticles to penetrate the membrane and the mechanism of nanoparticle–membrane interaction were investigated. The potential of the mean force profile, nanoparticle velocity profile, force profile and density profiles (planar and radial) were obtained to explore the nanoparticle permeation process. The structural properties of both nanoparticles and lipid membrane during the permeation, which are of considerable fundamental interest, are also studied in our work. The findings described in our work will lead to a better understanding of nanoparticle–lipid membrane interactions and cell cytotoxicity and help develop more efficient nanocarrier systems for intracellular delivery of therapeutics.


Introduction
Biological membranes are one of the major structural elements of cells, and play a key role as a selective barrier and substrate for many proteins that facilitate transport and signaling processes. These selective permeable membranes define the boundary and maintain the essential intracellular environment of the cell. Transport of nanoparticles across biological membranes is of significance in separation, bio-imaging and drug delivery systems. Recently, the available experimental techniques have been reported [1][2][3] to synthesize nanoparticles with well-controlled size, geometry and surface coatings, opening the door to use these nanoparticles for such applications effectively. Designing multifunctional nanoparticles for biomedical application requires a fundamental understanding of the interactions between nanoparticles and biological membranes, which remain largely unknown. It has been shown that such interactions are often determined by general physicochemical factors such as particle size, shape, hydrophobicity, surface charge density and the characteristics of the environment such as type of membrane and the interaction with other biological entities present in the system [4][5][6]. It has also been demonstrated in a few cases that the details of the molecular structure of the ligands and their distribution on the surface can have a profound effect on the mechanism by which the nanoparticle is transported across the membrane or incorporated into the cell [7,8]. The ability of these engineered nanoparticles, which are expected to be used increasingly in industry, to enter and be transported within biological bodies in ways that larger particles cannot, may have concomitant adverse toxicity effects [9][10][11]. Therefore, there is a clear need for physical insights into the questions regarding the permeation process of nanoparticles across biological membranes, which can provide an understanding of membrane structural changes during the permeation in general, and insight into how the nature of the interactions between lipid membrane and chemical species on the surface of the core particle in particular determines the details of membrane penetration by nanoparticles. Investigation of nanoparticle biocompatibility and toxicity has been a growing interest in addressing the impact of nanotechnology on human health and the life environment [12][13][14][15]. Recently, efforts have been focused on the interactions between nanomaterials and lipid membranes, experimentally [16][17][18][19] and theoretically [20][21][22][23][24][25][26].
Molecular dynamics (MD) is a powerful tool that can provide structural and dynamic details of the permeation process not readily available experimentally. Coarse-grained (CG) models where small groups of atoms are treated as single beads provide a promising method to study large biomolecular systems [27]. Marrink and co-workers [28] recently developed a coarse-grained force-field called the MARTINI force field for simulation of lipids and surfactants, and extended to amino acids and proteins [29]. The MARTINI force field has been shown to reproduce semi-quantitatively fundamental structural and thermodynamic properties of lipid bilayers and proteins.
Some molecular dynamics simulation studies have been conducted on the interactions of nanomaterials and biological membranes and a number of possible mechanisms have been reported recently. Those nanomaterials include fullerenes [22,[30][31][32], carbon nanotubes [33,34], bare nanoparticles [23,35,36] and functionalized nanoparticles [25,37]. The results from their initial studies are very promising. These theoretical studies of the translocation of nanomaterials have examined various aspects of the permeation process, including altered membrane thickness around the embedded nanomaterials, the phase transformations of lipid bilayers, the structural properties of the bilayer during permeation such as average order parameters of the tail, area per lipid and the thermodynamics of adsorption and penetration for nanomaterials. There is also an investigation of the effect of nanoparticle shape on its permeation through lipid bilayers using dissipative particle dynamics [24], in which the minimum driving forces required for various shaped nanoparticles (ellipsoids, cylinders, pushpin shapes) to translocate across the lipid bilayer was investigated. In summary, computer simulations in this field are at an early stage and many efforts are still underway to investigate the effects of nanostructured materials in biological environments.
In our previous work, we used a coarse-grained model to simulate the permeation of small molecules such as O 2 and CO 2 through pure DPPC lipid bilayers [38] and DPPC lipid bilayers with embedded Outer Membrane Protein A [39], with permeability results in satisfactory agreement with experimental results. We have also investigated the transport characteristics of nanocrystals (nanoparticles without surface ligands) across a model lipid membrane and the membrane response in terms of the structural and mechanical properties of the lipid membrane under the perturbation of nanocrystals of various sizes [40]. In the present work, we extend our studies to functionalized nanoparticles with surface ligands and explore how those ligands vary the details of the permeation process. Uncharged hydrophobic ligands of different lengths were investigated and an external force was added as the driving force to aid nanoparticles across the lipid bilayer. The external force applied on the nanoparticle mimics the force a nanoparticle can experience during transport in biological systems as a result of collisions and external magnetic fields [41], or optical forces acting on a gold nanoparticle under laser irradiation at the plasmon resonance [42]. In some experiments, penetration of nanoparticles across the bilayer has been aided by an external force, for example by using a nanoinjector [43]. The nanoparticle-lipid interaction was investigated in our simulations, in which the potential of the mean force profile, the minimum forces needed to penetrate the model membrane were obtained and the permeation characteristics for nanoparticles with and without ligands were examined. We also studied the structural properties of both the nanoparticles and lipid bilayers during the penetration, including local and bulk properties, which are of considerable fundamental interest. The findings described in the present work will lead to a better understanding of the permeation by ligand-coated nanoparticles and help in developing efficient nanocarrier systems for intracellular delivery of therapeutics, as well as understanding and predicting the cytotoxicity of nanoparticles. We plan to extend these studies to ligands that have hydrophilic sections to investigate possible nanoparticle penetration without the use of external forces.

Coarse-grained models
A coarse-grained model allows us to extend the space and time scales of simulations compared with the all-atoms model. The MARTINI force field [28] is one of the widely used Coarse-Grained (CG) models in MD simulations. More details about the MARTINI CG force field can be found in the literature [28,29].
Here, we will only introduce the model briefly.
In the MARTINI CG model, all particle pairs (in the Martini force field) i and j at distance r ij interact via a Lennard-Jones (L-J) potential: The well depth " ij depends on the interacting particle types and values range from " ij ¼ 5.6 kJ mol -1 for interactions between strong polar groups to " ij ¼ 2.0 kJ mol -1 for interactions between polar and apolar groups, mimicking the hydrophobic effect. The effective size of particles is governed by the L-J parameter ¼ 0.47 nm for all normal particle types, except that for the interaction between charged (Q type) and most apolar types (C1 and C2), the range of repulsion is extended by setting ¼ 0.62 nm. In addition to the L-J interaction, charged groups interact via a shifted Coulombic potential function: In the simulations described here, the non-bonded interactions are cut off at r cut ¼ 1.2 nm. The L-J potential is shifted from r shift ¼ 0.9 nm to 1.2 nm and the electrostatic potential is shifted from r shift ¼ 0.0 nm to 1.2 nm following a standard shift function [44]. The bonds are described by a harmonic potential V bond (R), and a cosine-type harmonic potential V angle ðÞ is used for bond angles,

Simulation of nanoparticles in solution
There is increasing interest in using functionalized gold nanoparticles as attractive nanomaterials for biological and biomedical applications, such as bio-imaging, single molecule tracking, drug delivery and diagnostic purposes [3,[45][46][47]. For example, gold nanoparticles can be engineered to accumulate preferentially in tumor cells using properly functionalized targeting ligands, thus providing an effective tool for cancer diagnosis and therapy [48]. These experiments motivate our work to use gold nanoparticles as our model nanoparticles. As in our previous study, the structure of the gold nanocrystal (nanoparticle without ligands (AuNP_bare)) is obtained by cutting a nearly spherical nanocrystal out of a bulk face-centered-cubic (FCC) structure gold lattice, with a diameter of 2.1 nm. For the present study, we attach ligands to the surface of such a 2.1 nm gold nanocrystal by the following procedure: we first put the nanocrystal in the center of a 12.0Â12.0Â12.0 nm 3 simulation box and fill the box with butanethiol ligands in large excess compared with the number required to form a compact monolayer. A cycled annealing simulation procedure is then followed to condense the ligands onto the surface of the nanocrystal, which is a method also reported by Luedtke et al. [49] using atomistic simulations. The temperature was subsequently increased from 200 to 500 K and then cycled between these temperatures to allow desorption of excess ligand molecules and exploration of stable binding sites. The final number of equilibrated butanethiol chains on the gold core is 87, yielding a thiolate per surface gold atom coverage of 48.3% (a surface density of 6.28 ligands nm -2 ), which is in the range of experimental coverage measurements, up to 52-57% for 2.1 nm diameter alkanethiolate gold nanoparticles [50]. As shown in Figure 1(a), after the annealing process, the surface sulfur atoms of the nanoparticles are uniformly distributed and the distances between sulfur atoms are in the range 0.44-0.51 nm (Figure 1(b)). For studies using variable ligand lengths, the butanethiol ligands (AuNP_SL) are replaced by R¼(CH 2 ) 8 and R¼(CH 2 ) 12 to form gold nanoparticles with neutral hydrophobic ligands of medium length (AuNP_ML) and long length (AuNP_LL), shown in Figure 1(c). The 2.1 nm diameter gold nanoparticle used here is smaller than typically used in biomedical applications. A wide range of sizes of Au nanoparticles have been used experimentally for drug delivery and as imaging agents. A size range of 1-100 nm or more is easily achieved [51]. For example, PEG-coated AuNPs (4 nm and 100 nm) have been administered intravenously to mice [52]. Pan et al. have investigated the size dependence of the cell toxicity of water-soluble gold nanoparticles ranging from 0.8 to 15 nm in diameter in four cell lines and find that all prove most sensitive to gold particles 1.4 nm in size [53]. Hainfeld [54] used Au NPs 1.9 nm in size for imaging in mice. Thus, our 2.1 nm nanoparticles, although small, are of practically relevant size. Gold nanoparticles used for biomedical applications (gene and drug delivery) are typically larger (20-100 nm) because they are conjugated with various biomolecules or drugs and these sizes are known to cross cell membranes efficiently using other mechanisms such as endocytosis, which is not studied in our work [6,19,55].
In our simulations, the atomistic to coarse-grained mapping strategy is as follows. The gold and sulfur atoms are mapped 1:1 and are fixed rigid and the residues of alkyl chains are 4:1 mapped and flexible (Figure 1(d)). The interaction between the gold nanoparticle and lipid molecules is described by L-J potentials. Various potential parameters have been used previously for gold atoms, from All-Atoms [56][57][58][59] to Coarse-Grained (atomistic structure, but modeling gold atoms as C-class [25] and P-class [60] using MARTINI force fields). Here, we assign MARTINI C5-type interaction sites for gold atoms, N0-type for sulfur atoms and C1-type for alkyl chains. The classes of interaction sites and their accompanying potential parameters have previously been tested and verified against atomistic simulations by Marrink et al. [28] and we have tested them in our simulation systems against experimental data for lipid membranes [39]. The non-bond and bond parameters for gold nanoparticles are shown in Table 1. For the crossinteractions between gold nanoparticles and lipid/ water sites, we use the standard Lorentz-Berthelot mixing rules [61] as a starting point in this study.
A range of properties of ligand-coated gold nanoparticles was examined to validate the effectiveness of the model we have adopted for the gold nanoparticles. A 10.8Â10.8Â10.8 nm 3 simulation box is created with explicit solvent and a gold nanoparticle placed at the center. First, we calculate the nanoparticle radius of gyration by letting each gold nanoparticle dissolve in hexane solution at 300 K. Two C1 beads having a 0.47 nm harmonic bond between them represent the CG hexane. The particle-solvent system was equilibrated for 20 ns with a time step of 10 fs and the radius of gyration is averaged every 400 ps. Second, we calculate the nanoparticle diffusion coefficient by letting each nanoparticle insert into a CDCl 3 solution. A C4 bead represents the CG CDCl 3 and the simulation setup is identical to the above case. The diffusion coefficient is obtained from the long-time slopes of mean square displacements. The nanoparticle radius of gyration and the diffusion coefficient for AuNP_ML and AuNP_LL nanoparticles, which are shown in Table 2, are consistent with Lin's coarse-grained simulation report [25]. We also notice that the radius of gyration of our nanoparticles is slightly higher than the experimental measurements, concomitant with a smaller diffusion coefficient in CDCl 3 solution compared with the experimental value [62]. These results are consistent with each other because a larger gyration radius will induce more friction in the solution, thus having a retarding effect on the diffusion. In summary, our simulation results are in reasonable agreement with available experimental results, which validates the effectiveness of the coarse-grained model we are using for the ligandcoated gold nanoparticles.

Simulation of lipid membrane with nanoparticles
We perform molecular dynamics simulations for a lipid membrane system with one nanoparticle. The model membrane in the present work consists of a DPPC (C 16 ) lipid bilayer. The lipid membrane/water system consists of 512 DPPC molecules and 11828 CG waters (four water molecules in each CG water) in a 12.6Â12.8Â16.4 nm 3 simulation box. We have previously shown that this bilayer membrane self-assembles from the isotropic solution of lipids in CG simulations [39] and that the properties of the self-assembled bilayers are in good agreement with experimental measurements [40], which validates the effectiveness of the coarse-grained model we are using for the lipid bilayers.
After allowing the equilibration of the lipid-water system (50 ns), we introduced one gold nanoparticle with flexible ligands into our simulation system, as shown in Figure 2. To make room for the nanoparticle, roughly 800 CG water molecules were removed. All the simulations were performed using the LAMMPS simulation package [63]. A Langevin thermostat [64] was applied in the NVT ensemble to maintain the desired temperature (323 K). To ensure stability, we used a time step of 10 fs. A typical simulation takes about 0.75 h per ns on an Intel Core2Quad Note: a MARTINI classification. Figure 2. Side and top view of the simulation system for investigating the transport of a nanoparticle across the DPPC lipid membrane (yellow represents the gold core, cyan the surface sulfur atoms, green the hydrophobic chains, blue the choline group, orange the phosphate group, magenta the glycerol group, red-orange the acyl chain tail group, white dots are water molecules).
CPU system. An external driving force was applied to aid the permeation of the nanoparticle across the lipid membrane. We note that the results we report here are for a single nanoparticle permeating the lipid membrane. The results observed may significantly change if several nanoparticles were permeating the membrane simultaneously.
3. Results and discussion 3.1. Permeation characteristics for bare and ligand-coated nanoparticles We first examine the minimum driving force needed for the nanoparticles to permeate through the lipid membrane. The minimum driving force for crossing the first and second layers of the lipid membrane are shown in Figure 3. In our simulations we have defined the minimum force as the smallest force required to permeate the membrane in 100 ns. Our minimum force for nanoparticles permeating across the first layer is in the range of 175-225 pN, and 350-550 pN for permeating both layers. Typical forces applied to single cells for AFM imaging, which are not large enough to cause cell rupture, are in the range 50-1200 pN [65]. For example, in the Vakarelski experiments, AFM 20-25 nm tips applied loads of only 100-200 pN [66]. The external forces used in our simulations are of the same order of magnitude. To permeate the first layer, the force needed for a nanoparticle with ligands (AuNP_SL, AuNP_ML and AuNP_LL) is smaller than for a nanoparticle without ligands (AuNP_bare). Compared with the bare nanoparticle, the ligands introduce more disruption in the first layer as the ligand-coated nanoparticles get close to the surface of the first layer under the same external forces, making it easier to open up the lipids to accept the nanoparticles into the bilayer. We shall see the details more clearly in the snapshots shown in Figure 5 below. For this type of nanoparticle in the range of ligand lengths studied, we also find the force needed for ligand-coated nanoparticles decreases a little bit with increasing ligand length on permeating the first layer. For the permeation of both first and second layers, the required force increases as the ligand length increases. Due to the interaction between hydrophobic ligands and hydrophobic lipid tails, a larger force is needed for ligand-coated nanoparticles to permeate the second layer, resulting in the larger force needed for the permeation. Therefore, the longer the length of the ligands, the larger the force needed for penetration. Since we are using the same gold core for all the simulations, the effective nanoparticle size is increasing with increasing ligand length. In a previous bare nanoparticle study [40], we also found the larger the nanocrystal, the larger the force need for permeation.
We also obtain and compare the potential of mean force (PMF) profiles for bare and ligand-coated nanoparticles. We constrain the nanoparticles and move the nanoparticles at a constant velocity of 0.1025 m s -1 through the simulation box. This nanoparticle velocity used here is much larger than experimentally obtainable (which would require simulations almost 400 times longer) [67]. However, we believe our simulations still represent the permeation process realistically. This is demonstrated by the recovery of the lipid layer between two permeation cycles, indicating no permanent damage to the membrane even at these high velocities [40]. Results obtained in other simulations using similar velocities also appear reasonable [34,68]. The force on the nanoparticles is measured and integrated over the positions to obtain the PMF profile. This standard integration method was also described in Ref. [69]. The PMF profile of nanoparticles as a function of position along the z-direction (direction of nanoparticle motion), calculated from the water phase to the center of the lipid phase, is shown in Figure 4. Our results show that the ligand-coated nanoparticles as well as the bare nanoparticle show a strong preference for the center of the membrane, however there are energy barriers at the water-lipid interface as the nanoparticle moves closer to the lipid membrane, preventing the nanoparticle from moving into the membrane spontaneously. Therefore, an external force is necessary for nanoparticles to penetrate into the membrane initially, and once the nanoparticles have permeated the head group they would move towards the center of the membrane. We note that if several nanoparticles were to attempt this simultaneously, it is conceivable that one would have enough energy to overcome this barrier spontaneously without the need for an artificial external force. The PMF profiles of our hydrophobic nanoparticles show agreement with those of fullerenes [22,30,32] and nanoparticles [25] having a hydrophobic nature in previous studies. The longer the ligand length, the larger the energy barrier encountered at the interface and the lower PMF is found in the interior of the lipid membrane. If Figure 3 is compared with Figure 4, there appears to be an apparent inconsistency. While the minimum force decreases with the length of the ligand (F min LL 5 F min ML 5 F min SL ), the energy barrier at the interface shows the opposite trend. The minimum force is a point property, while the PMF is an integral property. When the ligands are longer, the region of interaction between the nanoparticle and the interface increases, which leads to a larger energy barrier, even though the minimum force needed is lower for the longer ligands.
Another interesting observation relates to the position of the peak of the energy barriers. For a bare nanoparticle the barrier is at about 0.7 nm inside the interface. The bare nanoparticle touches the interface at a center-of-mass distance about 1.0 nm away from the surface. At that point the nanoparticle faces resistance from the lipid bilayer until it is about 0.7 nm into the interface. At this point the nanoparticle is essentially sucked into the center of the bilayer as the head groups begin to realign and the tails push the nanoparticle into the bilayer center. Snapshots of the permeation of a bare nanoparticle through the first layer of a lipid membrane are shown in Figure 5. For a nanoparticle with ligands, the ligands play a crucial role in opening up the lipid bilayer interface, so the energy barrier peak roughly corresponds to the nanoparticle center-of-mass being at the interface. The position of the peak is slightly inside the interface for the shortest ligands, since they cause the least disruption inside the lipid interface. Once the ligands are inside the bilayer they can easily open the lipid tails and move to the center of the bilayer. Figure 5 also shows a nanoparticle with ligands as it moves in the bilayer.
To compare the dynamic characteristics of nanoparticles with and without ligands during the permeation process, we plot the velocity and force profile along the direction of motion of the nanoparticles (z-direction). We consider AuNP_bare, AuNP_SL, AuNP_ML and AuNP_LL nanoparticles permeating under the same driving force (600 pN) to study how the ligands change the permeation dynamics. We obtain the force profile by letting the nanoparticles permeate at a constant velocity (0.41 m s -1 ) through the lipid membrane. Typical velocity profiles and force profiles    Figure 5).
for the permeation are shown in Figures 6(a) and (b). As we trace out the changes of the velocity and force profiles during the permeation process, our understanding of the molecular events is aided by examining the changes in the distribution of the lipid molecules in front of the nanoparticle (in the direction of its motion) at the same time. For this purpose, we obtained the lipid membrane xy-plane density profiles and the radial density profiles. These were obtained by calculating the number density of lipid molecules in front of the nanoparticles (in the direction of motion) and are shown in Figures 7(a) and (b), respectively. We discuss the profiles shown in Figure 6 enlightened by the lipid molecule distributions shown in Figure 7. For all nanoparticles, the velocity decreases when the nanoparticles approach the head groups of the lipid membrane and the resistance increases gradually in this region. Therefore, an external force would be necessary to permeate the membrane initially. In the entry region (2.5-4 nm in the z-axis in Figure 6), in order to permeate the first layer of the lipid membrane, the nanoparticle has to compress the first layer and push the head groups apart in order to make room for its cross-sectional area. Therefore, the velocity of the nanoparticle decreases in this region and the resistance from the head group increases continuously in this region accordingly.
As the nanoparticles move deeper into the membrane center (4-5.5 nm in the z-axis in Figure 6), the velocity of the nanoparticles slows down further. The deformation of the first layer gradually induces the deformation of the second layer, maintaining a pore in the direction of motion of the nanoparticle. Therefore, the size of the effective pore created in the membrane remains the same as the entry region (Figure 7(f)). Thus, similar velocity and force curves for AuNP_bare, AuNP_SL, AuNP_ML and AuNP_LL nanoparticles are obtained in this region.
After the bare and ligand-coated nanoparticles attach to the second layer and move through the exit region (5.5-7.5 nm in the z-axis in Figure 6), the xy-plane density profiles are not as informative as the radial density profiles. At the entrance, within the bilayer, and at the exit, the xy-plane density profiles look essentially the same in all regions. In contrast, the radial density profiles clearly exhibit systematic differences between these regions and also show differences that arise from different ligand lengths. As seen in our previous report [40], when the bare nanoparticle has totally crossed the first layer, the first layer starts to recover, the second layer is to be compressed and the tails of the second layer separate from each other, dragging the head group apart to form a pore in the second layer even before the nanoparticle arrives there. Accordingly, a speed-up in the velocity profile ( Figure 6(a)) and a significant drop in the force profile ( Figure 6(b)) are observed as the AuNP_bare nanoparticle passes through the second layer. The xy-plane lipid density profile (Figure 7(g)) shows recovery of the lipid membrane as the bare gold nanoparticle moves through the exit region. The peak at 1.7 nm in the radial lipid density profile (Figure 7(h)) drops significantly because the bare gold nanoparticle pushes the lipid molecules away and the number of lipid molecules in front of the nanoparticle decreases, resulting in a reduction of the effective pore size created in the lipid membrane in this region (Figure 7(i)).
In contrast, the ligand-coated gold nanoparticles (AuNP_SL, AuNP_ML and AuNP_LL) move more slowly in the exit region and exhibit a minimum in the velocity profile (Figure 6(a)) after these particles have crossed the second layer. Due to the attractive interactions between lipid tails and ligands, the resistance continues to rise for these nanoparticles and achieves a maximum after they cross the second layer ( Figure 6(b)). As shown in the xy-plane lipid density profile (Figure 7(g)), compared with the bare nanoparticle, the local decrease in the lipid density profile does not obviously recover, which means that the recovery of the lipid membrane is slower after permeation for the ligand-coated gold nanoparticle compared with that for bare gold particles. In the radial density profile (Figure 7(h)), we find a peak at a distance 1.7 nm from the center of the gold core and the magnitude of the peak remains the same, so the effective pore still exists in the exit region for ligand-coated nanoparticles (Figure 7(i)). The entanglement of ligands and lipid molecules causes the gold nanoparticles with ligands to move more slowly in this region compared with the bare gold particles. Thus, once the nanoparticles enter the exit region of the bilayer, the ligand-coated nanoparticles exhibit behavior different from the bare nanoparticle. Let us now examine more closely the effective pore created by the passage of the nanoparticle. We see that the density profiles in the xy-plane show similar behavior for all nanoparticles. The peak in the radial density profile is located at about 1.7 nm away from the center of the gold core for all three nanoparticles, which means similar effective pores are created in the lipid membrane by bare and ligand-coated nanoparticles (Figure 7(c)) in this region. We also observed that some water molecules entered the lipid bilayer region during the permeation of nanoparticles, in contrast to the water density profile for the unperturbed membrane. This occurs because the penetration of nanoparticles creates a pore in the lipid bilayer and allows some water molecules into the lipid region. This is further evidence of the disruption of the integrity of the lipid bilayer by the perturbation of a nanoparticle.
In Figure 8 we show the change in the effective pore size in the lipid membrane by obtaining the radial lipid density profile along the direction of motion of a nanoparticle permeating the second layer, in the region starting from the center of the lipid membrane (at 5.0 nm in Figure 6) to the head group of the second layer (at 7.0 nm in Figure 6). The density of lipid molecules around the bare nanoparticles (Figure 8(a)) drops significantly during its penetration of the second layer, indicating that the effective pore size decreases in this region. On the other hand, for nanoparticles with ligands (Figure 8(b)), no significant differences are found for the density of lipid molecules around the nanoparticles, indicating that the effective size of the pore is maintained in this region. The radial density profiles for ligands and lipids are shown together in Figure 8(c), illustrating that some lipid molecules move inside the nanoparticle ligand region, resulting in entanglement of lipid molecules with ligands. Because this occurs for ligand-coated nanoparticles, the effective size of the pore resulting from the boundary lipid molecules does not change significantly in the exit region during permeation. The velocity remains slow and the force is larger, as seen in Figure 6(a).
Furthermore, we studied the orientation of boundary lipid molecules during the nanoparticle permeation process. We define a vector pointing from the lipid head group to the bead at the end of the lipid tail and calculate the tilt angle between this vector and the z-axis (normal to the surface of the lipid membrane). The lipid molecular orientation profiles for the permeation of the nanoparticles are shown in Figures 9(a) and (b), for the entry region and the exit region, respectively. In the entry region, boundary lipid molecules have similar tilt angles when accepting both bare and ligand-coated nanoparticles. The angle for the majority of lipid molecules is in the range of 20-60 , while around 20% of the lipid molecules have smaller or larger angles. On the other hand, in the exit region, the lipid molecule orientations for AuNP_bare and AuNP_LL nanoparticles show different trends. For the AuNP_bare nanoparticle, most lipid molecules are observed to tilt in the range of 0-40 . Smaller numbers of lipid molecules are tilted by larger angles. For the bare gold nanoparticle, the lipid molecules in the exit region tilt less than in the entry region. This occurs because of the decreased dynamic space for lipid molecules in this region, resulting from the push back by the bare nanoparticle. In contrast, for the ligand-coated nanoparticle, the distribution of tilt angles is wider and more uniform, compared with the bare gold nanoparticle. Around 10% of lipid molecules (c) Radial lipid membrane and ligand density profiles for AuNP_LL nanoparticle (7 nm in Figure 5).
have large tilt angles (that is above 90 ), showing a tilt towards the membrane surface in this region. The tilting of lipid molecules maintains the effective size of the pore during nanoparticle permeation of the second layer. The interpretation of Figure 9(b) can also be seen in snapshots (c) and (d). We see from Figure 9(c) that the bare gold nanoparticle creates a hole before the nanoparticle reaches the head groups of the lipid membrane, and Figure 9(d) shows that the tails of lipid molecules move towards the ligand region and become entangled with the ligands, which results in the large tilt angle seen in Figure 9(b).

Internal order and structural properties of the lipid membrane and nanoparticles during permeation
It is known experimentally that the penetration of nanoparticles can affect the stability and the mechanical strength of lipid membranes [19]. Microscopy experiments have examined the formation of nanoscale holes caused by nanoparticles in model membranes.
Chen et al. [18] observed dendrimer nanoparticles making 3 nm diameter holes in living cell membranes. In our simulations, we observe structural changes of the lipid bilayer at a more fundamental level, including lipid membrane bulk and local properties. We observe the tail segment order parameter, the average tail length of the lipid molecule and the thickness of the lipid membrane, from which we can quantify the structural changes of lipid membranes during the permeation of different nanoparticles. We can also observe whether the structural properties of lipid membranes can recover or not after nanoparticle permeation, within our simulation period. We obtain the bulk and local structural properties, including the thickness of the lipid membrane, the tail length of the lipid molecule and the tail segment order parameter when bare and ligand-coated nanoparticles are in the lipid membrane region. We call it a bulk property because the value is an average over all lipids in the simulation, not only those lipid molecules close to the nanoparticle. For the local property, the value is an average over only those lipid molecules at the boundary of the nanoparticle. Figure 10 shows that the lipid membrane exhibits similar bulk properties during the permeation of bare and ligand-coated nanoparticles. The membrane thickness is shown for the pure membrane and under the perturbation of a bare or ligand-coated nanoparticle in Figure 10(a). The thickness of the lipid membrane starts to increase once the nanoparticles permeate the surface of the lipid membrane, compared with its equilibrium undisturbed condition. After the nanoparticles move into the membrane, the thickness of the membrane increases significantly from its equilibrium value. Longer ligands induce more perturbation in the lipid membrane, resulting in a greater increase in the thickness. We observe that the thickness of the lipid membrane recovers after perturbation by all three nanoparticles. This is expected since lipid membranes have been observed to self-assemble even from an initial random configuration of lipids in solution. Furthermore, we obtain the bulk and local average tail lengths of the lipid molecules by calculating lipid tail end-to-end distances. Accompanying the increase in the thickness of the lipid membrane, the tail length of the lipid in the bulk also increases during the permeation of nanoparticles, which is shown in Figure 10(b). We observe that the longer the ligand, the greater the increase in the membrane thickness and tail length it will induce. Moreover, the average length of the tails is shorter for those lipids closer to the nanoparticles because of compression by the nanoparticle during passage. Similar behavior is observed for both the bare nanoparticle and the ligand-coated nanoparticle.
We also calculated the order parameter of the tail segment of the lipid membrane to characterize the internal order of the lipid membrane when the nanoparticle is in the lipid membrane region. The order parameter in the bulk is shown in Figure 10(c). Overall, the bulk order parameters change only slightly upon insertion of the nanoparticles; minor structural changes are also observed in the work of Wong-Ekkabut et al. on the simulation of C 60 insertion into a bilayer [22]. The bulk order parameter of the tail segment is slightly larger than that for the unperturbed lipid membrane. We believe that this phenomenon is due to the decrease of the dynamic space for the lipid membrane after the insertion of the nanoparticle, which induces lower mobility and a smaller extent of isotropic averaging for lipid molecules of the entire lipid bilayer.
Furthermore, we calculated the order parameter for those lipid molecules 2.0 nm, 4.0 nm and 6.0 nm away from the center of the nanoparticles. It can be seen from Figure 11 that the local order parameter (bond 1) for bare and ligand-coated nanoparticles shows similar behavior when these particles move into the membrane and stay inside the membrane (Figures 11(a) and (b)), which means that the ligands do not have a significant effect on the local order parameter in these two regions. In the exit region (Figure 11(c)), the tails of the boundary lipids are more ordered during the permeation of a bare gold nanoparticle since the dynamic space for local lipids is more compressed. However, for ligand-coated nanoparticles, local lipid molecules maintain their order parameter by changing their orientation. The order parameters for bonds 2 and 3 show the same trends reported in Figure 11 and hence are not shown here.
In our simulations, we also observe the structural changes of the ligand-coated nanoparticle. We calculate the size of the nanoparticles, which is the average distance from the gold center to the last bead of the nanoparticle (Figure 12(a)). As for the order parameter  of the lipid membrane, we obtain the order parameter for the ligands to study how the structure of the nanoparticle itself changes during permeation (Figure 12(b)). It is observed that the size of the nanoparticle increases when it moves into the lipid membrane region, which accompanies the stretching of surface ligands when the particle is in this region. The hydrophobic property of nanoparticles makes it energetically favorable to stay inside the membrane. At the same time, due to the decrease of the dynamic space, the ligands become more ordered, compared with their equilibrium value in the aqueous phase. When the nanoparticles move out of the membrane, the size of the nanoparticle reaches a maximum. This occurs because the drag force from lipid molecules makes the ligands stretch more, which may significantly disturb the structure of the ligands of the nanoparticle, making the nanoparticle less compact. We observe that the ligands become less ordered in this region. Furthermore, we notice that the effective size increase of the nanoparticles is smaller for those with longer ligands. Also, the bond closer to the surface of the gold core (bond 1) is more ordered than the remote ligand bond (bond 2), which indicates that, during the permeation process, a nanoparticle with short ligands is more rigid than a nanoparticle with long ligands.

Conclusions
In summary, we have investigated the translocation of nanoparticles across a lipid membrane using coarsegrained molecular dynamics simulations. We compare the different permeation behaviors for nanoparticles without surface ligands and functionalized    nanoparticles with different ligand lengths and study how the surface ligands change the details of the permeation process. The velocity profile of the nanoparticles during the translocation across the lipid membrane is obtained in the present work, as well as the force profile arising from the interaction between the nanoparticles and the lipid molecules. We find that the minimum force for penetrating the first layer is smaller while the minimum force for permeating both layers is larger for functionalized nanoparticles with neutral hydrophobic surface ligands. The effective size of the pore created by a nanoparticle in the lipid membrane remains the same as the nanoparticle moves inside the membrane. This effective pore size is essentially independent of the length of the ligand within our studied range. The effective pore size decreases when bare nanoparticles (having no ligands) move out of the membrane, while remaining the same for nanoparticles having surface ligands. Accordingly, the tilt angle for lipid molecules close to the nanoparticle is similar for all nanoparticles in the entry region. We observe larger lipid tilt angles in the local lipid molecules when functionalized nanoparticles are in the exit region, which is not observed for bare nanoparticles. Observation of the lipid membrane xy-plane density profile shows that the membrane recovery has a longer timescale after penetration by functionalized nanoparticles than for bare nanoparticles. The thickness of the lipid membrane and the lipid tail length recover after the permeation is completed, which points to the elastic nature of the lipid bilayer membrane. We find that the lipid tail segment of bulk lipids changes only slightly during the permeation for all nanoparticles. In the exit region, the tails of the surrounding lipids are more ordered during the permeation of bare nanoparticles than for nanoparticles with ligands. Finally, we studied the structural properties of functionalized nanoparticles themselves. We found that the size of the nanoparticles and their ligand order parameters increase when they move into the membrane. The size of the nanoparticles reaches a maximum value in the exit region, while the ligand order parameter decreases significantly in this region. All findings are consistent with the ability of lipid membranes to recover after penetration by nanoparticles, not only for those nanoparticles without strong or specific interactions with lipid molecules, but also for nanoparticles having strong interactions with lipid molecules. These results are for neutral ligands uniformly distributed on the gold surface. We plan to investigate nanoparticles functionalized with other types of ligands, including charged and amphiphilic ligands, which may have different permeation characteristics. We also plan to investigate penetration by nanoparticles with distributions of non-uniform surface ligands. The findings of our work will lead to a better understanding of the mechanisms for the translocation of nanoparticles across lipid membranes, may help to develop efficient nanocarrier systems for intracellular delivery of therapeutics, and may also help us gain a better understanding of the mechanisms for cytotoxicity of some types of nanoparticles.