Effect of Shape on the Entering of Graphene Quantum Dots into a Membrane: A Molecular Dynamics Simulation.

Graphene quantum dots (GQDs), a new quasi-zero-dimensional nanomaterial, have the advantages of a smaller transverse size, better biocompatibility, and lower toxicity. They have potential applications in biosensors, drug delivery, and biological imaging. Therefore, it is particularly important to understand the transport mechanism of the GQDs on the cell membrane. In particular, the effect of the GQD shapes on the translocation mechanism should be well understood. In this study, the permeation process of the GQDs with different shapes through a 1-palmitoyl-2-oleoylphosphatidylcholine membrane was studied using molecular dynamics. The results show that all small-sized GQDs with different shapes translocated through the lipid membrane at a nanosecond timescale. The GQDs tend to remain on the surface of the cell membrane; then, the corners of the GQDs spontaneously enter the cell membrane; and, finally, the entire GQDs enter the cell membrane and tend to stabilize in the middle of the cell membrane. Moreover, the GQDs do not induce notable damage to the cell membrane, indicating that they are less toxic to cells and can be used as a potential biomedical material.

Introduction

As the most promising new materials in the 21st century, nanomaterials have been widely used in the chemical industry, medical treatment, the aerospace industry, and other fields.[1−5] Graphene is called “black gold” because it is the thinnest and hardest nanomaterial and also has the highest electrical and thermal conductivities among the recently discovered nanomaterials. It is also considered to be a revolutionary material for the development of new technology in the 21st century.[6,7] Since the planar two-dimensional graphene material is composed of only sp2-hybridized carbon atoms, it exhibits unique electrical and mechanical properties.[8,9] Unlike for graphene, the radial size of graphene quantum dots (GQDs) is smaller than 100 nm.[10] The ultrasmall size and rich edge effects endow GQDs with performance characteristics different from those of graphene.[11−13] For example, GQDs have a larger specific surface area, higher activity, and stronger natural absorption.[14] The high specific surface area makes GQDs an ideal drug carrier in medicine, and their high activity and natural absorption make them better absorbed by cells.[15] Qin experimentally reported that GQDs can guide the delivery of ovarian drugs.[16] Fang et al. pointed out that GQDs inhibit the premature release of drugs, enhancing the efficacy of chemo-photothermal therapy and inhibiting tumor growth.[17]

However, in recent years, with the growing potential of GQDs in biological applications, their biosafety has attracted extensive attention.[18,19] Experiments by Chong et al found that GQDs had no obvious effect on mice.[20] Similarly, Tu et al found in experiments that GQDs with a large size can induce the degradation of Escherichia coli cell membranes, leading to bacterial apoptosis.[21] Investigations of GQDs with different surface modifications showed that compared to GQDs containing the −OH and −NH functional groups, the GQDs containing −COOH showed better biocompatibility and could be considered for biological applications.[22]

Although a large number of studies have experimentally proved the biosafety of GQDs, the mechanism of the translocation of GQDs through the lipid membrane needs to be explored much deeper. Molecular dynamics (MD) simulations have been widely used in the field of biomaterials and nanomaterials to study complex dynamic problems.[23−25] In particular, Wang et al. found that GQDs effectively inhibit the accumulation of insulin starch in the human body, indicating that GQDs have the potential to alleviate starch degeneration.[26] Chen et al. pointed out that the phospholipid extraction by graphene can activate integrin on the cell membrane.[27] Li et al. first pointed out that graphene could spontaneously pass through the cell membrane through edge protrusion or corner sites.[28] Puigpelat et al. used the MD method to study the spontaneous entry of graphene parallelograms of different sizes into the cell membrane and found that the graphene sheets tended to migrate to the disordered phase with less cholesterol content.[29] Our previous work indicated that the transport of GQDs on the cell membrane was affected by their size.[30] Li et al. used coarse-grained MD and all-atom steered MD simulations to find that the sharpest corner of rhombic graphene has the lowest energy barrier and is the most preferred entry.[28] In addition, Mohanty et al. could create nanostructures with predetermined shapes (square, rectangle, triangle, and ribbon) using nanocutting technology.[31] However, there is little research on whether the shape of GQDs affects the transport of graphene quantum dots across the cell membrane. In this study, we used MD simulations to study the transport mechanisms of GQDs with different shapes on the 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) membrane.

Results and Discussion

Translocation Phenomenon

It was found in the simulations that all types of GQD shapes can permeate the membrane. First, the change of angle and the center of mass distance changes for the z coordinates (perpendicular to the cell membrane) between the GQDs and the POPC membrane were calculated, as shown in Figure c. At the beginning of the simulation, all the GQDs are located at 4.5 nm along the z coordinate. The small circular Cir-GQD37 and square Squ-GQD46 reach the surface of the POPC membrane in the shortest time but show two different transport modes. When Cir-GQD37 first came into contact with the cell membrane, the edges of GQDs were captured vertically by the cell membrane along the vertical direction and quickly through the cell membrane without a parallel attachment process. Unlike Cir-GQD37, when the square Squ-GQD46 reaches the surface of the POPC membrane, it is in a horizontal state with the surface of the membranes, and then GQDs adhered tightly to the cell membrane until 70 ns. As shown in Figure a,b, the snapshots of Cir-GQD37 and Squ-GQD46 at t = 4.5 ns were, respectively, taken. It can be seen that half of Cir-GQD37 was inside the POPC membrane at this time, while Squ-GQD46 was attached to the upper surface of the POPC membrane during this time. Similar to these two systems, rectangular Rec-GQD50 enter through the POPC membrane in a vertical state, remains on the surface of the membrane for a short time, and almost begins to penetrate through the membrane when it reaches the surface of the POPC membrane. For other systems, we found that the GQDs show parallel adsorption onto the membrane and remain on the cell membrane for a period of time. Although they cross the membrane in different ways, all GQDs can penetrate through the POPC membrane within 100 ns.

Figure 1

GQDs on the membrane at 4.5 ns MD simulation: (a) Cir-GQD37 and (b) Squ-GQD46. P (red) and N (yellow) atoms in the membrane are also shown in the VMD model. (c) Distance between the center of the mass of different GQDs and the center of the mass of the lipid membrane as a function of simulation time, where the green dotted line shows the top of the membrane.

To describe the state change of the GQDs entering the cell membrane in detail, we calculated the angle change between GQDs and the cell membrane in the simulation process. As shown in Figure , at the beginning of the simulation, θ = 90°, indicating that the GQDs are perpendicular to the membrane surface. In the simulation, the angle fluctuation was relatively large during the first 20 ns; this was due to the diffusion of GQDs in water. After the first 20 ns of the simulation, the change in the angle of all the GQDs can be divided into three different states. For circular Cir-GQD37 and rectangular Rec-GQD50, the angle range is mainly from 70° to 90°, combined with the change in the centroid coordinates, as shown in Figure , which indicates that the GQDs have entered the POPC membrane. However, the change in the angle of square Squ-GQD46 and parallel Par-GQD49 varies from 0° to 20°, combined with the change in the centroid coordinates shown in Figure and the observation of the motion state of the GQDs; at this time, Squ-GQD46 and Par-GQD49 are attached to the top of the membrane. Unlike for other simulation systems, for the Tri-GQD45, the angle and z coordination change dramatically before 40 ns because it still diffuses in the water, and it was also adsorbed onto the top of the membrane after 40 ns.

Figure 2

Angles between the GQDs with different shapes and the lipid membrane as a function of the simulation time.

Translocation Dynamics

To clearly show the process of GQDs entering the POPC membrane and the movement state in the cell membrane, we take the parallelogram Par-GQD49 as an example to explore the processes. As shown in Figure d, the Par-GQD49 spontaneously adsorbed on the POPC membrane from the free state and adsorbed on the upper surface of the membrane, similar to the contact mode of triangular Tri-GQD45 and square Squ-GQD46 contacting the membrane. We selected the instantaneous snapshot in the phase where the GQDs entered the membrane, and the Par-GQD49 was still oriented horizontally on the membrane at 68.9 ns. Under the influence of the membrane, one corner of the Par-GQD49 shows flexible bending at 69.1 ns, and its corner has entered the hydrophobic region of the phospholipid membrane layer. With an increase in the simulation time, the angle of the Par-GQDs begins to change, and the carbon atoms of Par-GQD49 that did not enter the membrane at 71.3 ns bend at 45° with the membrane under the action of the force. After the simulation was carried out for 71.3 ns, the Par-GQD49 was completely perpendicular to the plane of the membrane, and half of the volume of Par-GQD49 was in the POPC membrane. Within the next 0.2 ns, the Par-GQD49 had fully penetrated the membrane.

Figure 3

(a) Initial positions of the membrane and GQDs. (b) Model of GQD intercalation into the membrane. (c) Interaction energy of GQDs with the membrane as a function of distance. (d) Center of mass coordinate in the z direction of GQDs (blue line) and POPC (red line) membrane, and the upper part of the figure shows the instantaneous snapshots of the GQDs at different time periods.

The GQDs can enter the membrane spontaneously, and the driving force is the interaction between the hydrophobic ester chain of the membrane and GQDs.[29] Similar to previous MD studies, the transport of the GQDs in the lipid membrane starts preferentially at the corner or at an irregular position, and the rough edges or the corners of the GQDs preferentially penetrate the membrane because of the action of the acyl chain in the membrane; after that the GQDs are rapidly pulled into the membrane.[28,32−34] In line with previous studies,[28] five variables were used to describe the process of membrane penetration: h (the height of the GQDs entering the membrane), hH (thickness of the head in the lipid monolayer), hT (thickness of the tail in the lipid monolayer), γH (energy density of the interaction between the GQDs and the head groups of lipids), and γT (energy density of the interaction between the GQDs and the tail groups of lipids).[28] For the Par-GQD49 entering the membrane, the energy change can be expressed as a piecewise function of the penetration depth h, as shown in eq , which is derived from the study of Li et al.[28]

The process of entering the membrane can be divided into two phases: the GQD cusp entering the head of the phospholipid (0 < h ≤ hH) and the cusp entering the tail of the phospholipid molecular layer (hH < h < hH + hT). The surface interaction energy densities γT and γH were estimated to be −7kBT and 7kBT nm–2,[35] respectively. For the hydrophilic and hydrophobic thickness of the membrane, the total calculation was carried out with the Gromacs 5.0 software, and the average value was taken at last. The hydrophilic hH and hydrophobic hT of the measured membrane were 0.5 and 1.5 nm, respectively; the thickness of the entire membrane was 4 nm; and the GQDs have a sharp angle of 60°. In the first regime, the energy change is E = 2h2γH tan α. Obviously, with an increase in GQDs entering the hydrophilic region of membrane, the energy required is increasing. The energy change in the second stage hH < h < hH + hT first increases and then decreases, and the peak energy appears at dE/dH = 0. This indicates that there is an energy barrier E = 2(1 – γH/γT)hH2γH tan α when the depth of the second stage is h = (1 – γH/γT)hH. The critical penetration depth h ≈ 1.0 nm and the critical energy E ≈ 4.04kBT can be obtained by introducing the formulawith the surface interaction energy densities γT and γH, which were estimated to be −7kBT and 7kBT nm–2. Using the known parameters, as shown in Figure c, eq is plotted as a function of the binding energy and distance. For h > 1 nm, the binding energy shows a downward trend, indicating that when the GQDs cross the E = 4.04kBT energy barrier, they enter the cell membrane spontaneously. When the angle α of GQDs approaches 90°, the energy E = 2(1 – γH/γT)hH2γH tan α value is infinite; hence, the energy needed to enter the membrane is also infinite. Obviously, GQDs cannot enter the cell membrane along 180°. Although the GQDs can enter the cell membrane along their small sharp corner according to the previous calculation theory and the insertion angle is the key factor affecting the insertion process, the size of GQDs should also be considered in the actual simulation. In our previous simulation process,[30] we found that the smaller GQDs easily enter the cell membrane. Herein, the Cir-GQD37 is the smallest size among all GQDs; thus, Cir-GQD37 first enters the cell membrane.

As shown in Figure , to observe the stable state of GQDs more clearly, we selected the state of the GQDs in the membrane when the simulation was carried out for 100 ns. Although all of the GQDs permeate into the membrane, they do not remain in the center of the membrane but rather are stable at the higher positions in the membrane. The lowest energy position in the membrane is located at both ends of the POPC membrane, explaining why the GQDs remain at a higher position in the POPC membrane after entering the membrane, which is consistent with our previous results.[30] By observing the final stable state of the GQDs, they tend to be parallel to the main chain of the POPC molecule, which is consistent with the angle measurement shown in Figure . To further explore the interaction mechanism between the GQDs and the POPC membrane, we calculated the interaction force between the GQDs and the cell membrane. As shown in Figure , the van der Waals (vdW) force and electrostatic interaction force were calculated to explore the interaction between the GQDs and the POPC membrane. Analysis and comparison of these two interaction forces show that the vdW interaction accounts for a large part of the binding energy, while the contribution of the electrostatic interaction is relatively small. To more intuitively show the magnitude of the two interactions in the stable state, Table lists the average values of these two interactions for the simulation times from 80 to 100 ns in detail. vdW forces comprise more than 95% of the total binding energy, indicating that the electrostatic interactions in the simulation process are weak, and the main driving force is the vdW force. However, considering that Squ-GQD46 has not fully reached a stable state, the calculated results have a certain deviation. In Figure a, the line with zero potential indicates that there is no contact between the GQDs and the POPC membrane. It is observed that the GQDs and the POPC membrane have a stage that maintains an interaction energy of 200 kJ/mol, indicating that the GQDs are adsorbed on the membrane. In the process of entering the POPC membrane, the vdW force between the GQDs and the POPC membrane increases rapidly to reach a steady state. Because of the different shapes and sizes of the GQDs, the energy levels after the stabilization are also different, and the vdW interaction energy is maintained at 600–800 kJ/mol; further observation of Table shows that the magnitude of the binding force is greatly affected by the size of the GQDs, whereas the shape of the GQDs has almost no effect on the binding force.

Figure 4

Snapshots at ultimate time GQDs with different shapes in the membrane after 100 ns MD simulation: (a) protocell membrane, (b) Cir-GQD37, (c) Par-GQD49, (d) Rec-GQD50, (e) Squ-GQD46, and (f) Tri-GQD45. The GQD and POPC are shown by the VMD model. GQDs are represented by a red VMD model, and the P atoms in the membrane are represented by a blue VMD model. Water molecules are hidden for clearer display.

Figure 5

Interaction energy between the GQDs and the POPC membrane. (a) vdW interaction energy. (b) Electrostatic interaction energy.

Table 1

Details of the Binding Interaction between the GQDs and POPC Membranes

system	total (kJ/mol)	vdW (kJ/mol)	electrostatic (ele) (kJ/mol)	vdW/total (%)	ele/total (%)
Cir-GQD37	–685.46	–627.89	–30.57	95.37	4.63
Par-GQD49	–819.82	–790.36	–29.45	96.41	3.59
Rec-GQD50	–854.64	–825.10	–29.53	96.54	3.46
Squ-GQD46	–707.67	–631.92	–75.75	89.30	10.7
Tri-GQD45	–756.15	–739.21	–16.94	97.78	2.22

Table 2

Details of the Simulations Performed in This Study

system	shapes of GQDs	number of atoms	simulation time (ns)	GQDs (H:C)	GQDs:POPC
Cir-GQD37	circular	125 364	100	24:96	1:256
Par-GQD49	parallelogram	125 358	100	30:126	1:256
Rec-GQD50	rectangle	125 346	100	32:130	1:256
Squ-GQD46	square	125 352	100	30:120	1:256
Tri-GQD45	triangle	125 334	100	30:118	1:256
Par-GQD169	parallelogram	126 407	100	54:390	1:256
Squ-GQD218	square	125 297	100	62:496	1:256

Effect of GQD on the Membrane

To investigate the effect of the GQD shape on the membrane, we calculated the order parameters of the lipid membrane in the simulation process. The structure of the lipid membrane in the liquid is very flexible; therefore, we introduce an order parameter that represents the orientation and sensitivity of the molecular arrangement of the cell membrane that in turn can be related to the configuration entropy. The order parameters were calculated using the method described in the literature.[36,37] The order parameter formula of our previous calculation was used to perform the calculations for the molecules of the POPC membrane. The final formula is given by[38]where, in eq , m represents the atomic mass of the ith atom in the jth lipid molecule and r is the distance vector between the ith atom in the jth lipid molecule and the centroid of the jth lipid molecule. α and β indicate two Cartesian (xyz) directions, and rα is the α direction vector of r. The long axis of a single phospholipid molecule was defined as the eigenvector α. Using the eigenvector a of the major axis of molecules, all molecules can be represented by diagonalized order tensors, as shown in eqs and 5, where Nm is the number of the lipid molecules in the POPC bilayer membrane, and λmax is the order parameter Sorder.

As shown in Figure , for better comparison, a set of POPC membranes without the GQD system was also studied as a control experiment. The deuterium order parameter of the initial POPC membrane was approximately 0.72. As the simulation progresses, we find that this value does not change significantly and remains between 0.7 and 0.8; moreover, it can be observed that with the entry of the GQDs, the order parameters do not change significantly. Additionally, compared with the order parameters of the films without GQDs, the change in the structure of the POPC membrane is relatively small during the entire simulation phase. This indicates that the influence of GQDs with different shapes and sizes on the POPC membrane was relatively small.

Figure 6

Order parameter of the membrane as a function of simulation time in all of the simulations.

To better understand the effect of GQDs on the POPC membrane atoms, we calculated the density of water, GQDs, and P atoms in the z direction. As shown in Figure S1 in the Supporting Information, the black line represents the density of the P atom and the distance between the peaks indicates the thickness of the POPC membrane. It is observed that the P atom and water densities in the membrane did not change significantly. Water molecules and P atoms are distributed symmetrically on both sides of the membrane, indicating that the GQDs do not destroy the structure of the POPC membrane, further confirming the previous calculation results for the Sorder of the membrane. We also present the density distribution of the GQD atoms that can be observed in the main molecules of the GQDs at 6–8 nm along the z direction, which is the higher position in the POPC membrane; the obtained results are consistent with the results of the previous calculations.

Dynamic Simulation of Larger GQDs

To explore the effect of the quantum dot size on the cell membrane, we constructed larger diamond-shaped and square GQDs, named Par-GQD125 and Squ-GQD218. Similar to the previous simulation process, we performed MD simulations for 100 ns. The simulation results are shown in Figure . The GQDs completely entered the membrane at 100 ns. Because the GQDs in the parallelogram are larger than the thickness of the POPC membrane, the position of the upper end of the GQDs entering the membrane is exposed to the outside. The square GQDs were completely inserted into the membrane when the simulation was completed. Because these GQDs are larger than those studied previously, there is no phase of the GQDs spreading to the membrane in the simulation. It is important to note that, consistent with the previous results, the first corner of the GQD enters the POPC membrane and then whole GQD enters the membrane.

Figure 7

Snapshots of Par-GQD125 and Squ-GQD218 entering the cell membrane. (a and b) Initial structures and (c and d) structures at 100 ns.

Computational Details

System Setup

To study the transport mechanism of GQDs with different shapes and sizes through the cell membranes, we built a complete zigzag graphene sheet using the visual molecular dynamics (VMD) software[39] and obtained different shapes of graphene sheets by cutting the complete sheet with Materials Studio (BIOVIA). As shown in Figure , to ensure that the sizes of the different shapes of GQDs are as close as possible to each other, the number of C atoms in cut graphene sheets is kept as close as possible to the same value for all GQDs with the same size. The edges of all of the GQDs were saturated with hydrogen atoms, and the GQDs were labeled based on the total number of their benzene rings. The size of the GQDs is defined according to the number of carbon rings and is preceded by an abbreviation for shape, e.g., Cir-GQD37, Par-GQD49, Rec-GQD50, Squ-GQD46, and Tri-GQD45. All of the GQD structures were optimized using the Gaussian 09 software (Gaussian, Inc.),[40] and the optimized GQD structure was used as the initial structure for the MD simulation. Using a bilayer membrane containing 256 POPC molecules (156 phospholipids per layer), the upper and lower layers were symmetrically distributed. The structure of the POPC membrane was optimized by an isothermal−isobaric (NPT) simulation carried out for 100 ns. After equilibration, the POPC membrane was used as the initialization structure. The initialized GQDs were placed 0.2 nm above the POPC membrane. Then, the system was dissolved in TIP4P water in a simulation box with the dimensions of 9.40 × 9.30 × 14.28 nm3. The instantaneous images of the MD simulation results were obtained using the VMD software.[39]

Figure 8

GQDs with different shapes: Cir-GQD37, Par-GQD49, RecGQD50, Squ-GQD50, Squ-GQD46, and Tri-GQD45.

MD Simulations

All simulations were performed using the open source Gromacs 5.0 software package, using the force-field parameters of Charmm36. As in our previous work,[30] the outermost hydrogen atoms and the carbon atoms connected to the hydrogen atoms carry charges of +0.115e and −0.115e, respectively, and the remaining carbon atoms are neutral. The harmonic angle, harmonic potential, and dihedral angle of the GQDs were all derived from Cohen-Tanugi and Grossman[41] Prior to the MD simulation, the energy was minimized for all systems, and the equilibrium ensemble was used to ensure that the simulation system reached a stable state. The semi-isotropic Parrinello–Rahman constant pressure method was used to control the pressure, and the temperature was maintained at 310 K using the Berendsen thermostat. The bond length and bond angle of the H atom are limited by the LINCS algorithm in MD simulations.[30] The transport of small-sized GQDs in the POPC membrane was simulated using the NPT simulation for 100 ns. To explore the influence of the size of a single shape on the cell membrane, we increased the size of parallelogram and square GQDs, named Par-GQD169 and Squ-GQD218, to explore whether the large GQDs can enter the cell membrane at 100 ns. Same as the previous simulation, the simulation step size was 2 fs, and the trajectory file was saved every 2 fs.

Conclusions

In summary, we explored the translocation of GQDs with different shapes and sizes in the POPC membrane by MD simulations. In the simulations carried out for 100 ns, all of the GQDs penetrated the membrane regardless of their shape. During the simulation process, the GQDs were first adsorbed parallel to the membrane. Exploiting their flexibility, the GQDs bend their structure and permeate the POPC membrane. After entering the membrane, because of the vdW and Coulomb interactions, small GQDs tend to be vertically dispersed on the side of the membrane without being adsorbed onto the center of the membrane, and large GQDs completely entered the membrane. Individual, small GQDs induce relatively little mechanical damage to the membrane, and large and angled GQDs can also enter the membrane, damaging the membrane to some degrees. In our simulation, all the atoms are in a completely free state, and the cell membrane is also considered as a flat membrane. However, in the actual body, change in the curvature of the cell membrane is considered to be a dynamic change. Therefore, whether the change in the curvature of the cell membrane will affect the process of GQDs entering into the cell membrane will be further discussed.