Molecular Dynamics Study of Stability and Diffusion of Graphene-Based Drug Delivery Systems

Graphene, a two-dimensional nanomaterial with unique biomedical properties, has attracted great attention due to its potential applications in graphene-based drug delivery systems (DDS). In this work graphene sheets with various sizes and graphene oxide functionalized with polyethylene glycol (GO-PEG) are utilized as nanocarriers to load anticancer drug molecules including CE6, DOX, MTX, and SN38. We carried out molecular dynamics calculations to explore the energetic stabilities and diffusion behaviors of the complex systems with focuses on the effects of the sizes and functionalization of graphene sheets as well as the number and types of drug molecules. Our study shows that the binding of graphene-drug complex is favorable when the drug molecules and finite graphene sheets become comparable in sizes. The boundaries of finite sized graphene sheets restrict the movement of drug molecules. The double-side loading often slows down the diffusion of drug molecules compared with the single-side loading. The drug molecules bind more strongly with GO-PEG than with pristine graphene sheets, demonstrating the advantages of functionalization in improving the stability and biocompatibility of graphene-based DDS.


Introduction
The clinical use of various potent hydrophobic molecules, many of which are aromatic, is often hampered by their poor water solubility and low biocompatibility.Although water-soluble prodrugs may circumvent these problems, the efficacy of the drugs decreases.Nanomaterial-based drug carriers have become a hot spot of research at the interface between nanotechnology and biomedicine because the nanocarriers allow efficient loading, targeted delivery, and controlled release of drugs [1][2][3][4].Carbon nanotubes (CNTs) have attracted enormous interests in biomedicine as nanocarriers and much effort has been done in in vitro and in vivo biological applications recently [5][6][7][8].Graphene has emerged as a two-dimensional (2D) honeycomb lattice with interesting physical properties since it was synthesized in 2004 [9].The intensive research efforts are ongoing to investigate the potential applications of graphene in various fields including quantum physics, transparent conductors, nanoelectronic devices, nanocomposite materials, catalysis, energy research, and biomedicine [10][11][12].Graphene is a 2D monolayer nanomaterial composed of carbon atoms forming connected six-membered rings.Compared with CNTs, graphene sheets possess long edges and large accessible surfaces with delocalized  electrons [13,14], leading to excellent ability to immobilize a large number of substances including metals, drugs, biomolecules, and fluorescent molecules [15].Graphene and its derivatives, such as graphene oxide (GO), reduced graphene oxide (RGO), and GO-nanocomposites, with every atom and functional group exposed to environments, allow ultra-high drug loading efficiency and improve greatly the water solubility and targeting of drug molecules.Thus graphene-based materials may be potentially utilized in cancer treatments with encouraging therapeutic outcomes [16][17][18][19][20][21].
Tremendous efforts have been made on modification of graphene and its derivatives for various biomedical applications.Liu and coworkers [22] prepared PEGylated graphene 2 Journal of Nanomaterials oxide (NGO-PEG) sheets with ultra-small nanosizes (10-50 nm) by conjugating the amino groups of PEG with the carboxyl groups of GO.The NGO-PEG sheets exhibit high stability in physiological solutions and can absorb aromatic polymers such as SN38 via - interactions to improve the solubility of the carrier-drug systems.This research showed that the functionalized nanographene sheets are biocompatible without obvious toxicity and could be used potentially for drug delivery.Zhang and his coworkers [23] modified graphene oxides by sulfonic acid derivatives.The modified GO were used as nanocarriers to load two types of anticancer drugs (DOX and SN38) for targeted drug delivery.Other groups have used a series of macromolecular polymers to modify graphene oxide by covalent chemical reactions, including polyethyleneimine (PEI), chitosan (CS), poval (PVA), poly(N-isopropyl acrylamide) (PNIPAM), and amine-modified dextran (DEX).Another synthesis strategy is noncovalent functionalization via hydrophobic interactions, - interaction, or electrostatic binding [23][24][25][26][27][28][29][30][31].These experimental researches demonstrate that graphene can be used as a nanocarrier to load drug molecules and improve the solubility of carrier-drug systems effectively when functionalized with various hydrophilic molecules or polymers, implying potential applications in clinical treatments.Previous theoretical works on similar systems include the study of the interactions between biomolecules and graphene: Qin et al. [32] performed density functional theory (DFT) and molecular dynamics (MD) calculations to explore the interaction between graphene and L-leucine, showing that van der Waals interactions between the L-leucine and graphene play a dominant role in the absorption process.
Despite the experimental and theoretical studies mentioned above, the thermodynamics and kinetics of graphenedrug complex remain obscure at a molecular level.For noncovalent graphene-drug complexes, the major practical concerns are whether the drug molecules can be loaded onto the graphene sheets tightly.To design stable graphene-based drug delivery systems (DDS), therefore, it is necessary to investigate thermodynamics stability from binding energy point of view.Moreover, the diffusion behaviors of drug molecules on the graphene sheets can be used to measure the kinetic stability of DDS.In this work, we employed molecular dynamics method to study the interaction between graphene and drug molecules focusing on the various sizes of graphene sheets as well as the types, numbers, and loading modes of drugs.Specifically, we studied four types of anticancer drug molecules (CE6, DOX, MTX, and SN38, shown in Figure 1) absorbed onto graphene sheets with various sizes (10 Å, 20 Å, 30 Å, 40 Å, and 50 Å finite and periodic graphene sheets) and graphene oxide functionalized with aminemodified polyethylene glycol (GO-PEG).In addition to the absorption mechanism of drug molecules on graphene, we also investigated the binding strength and diffusion behavior of various graphene-based drug delivery systems.

Modeling and Geometry Optimizations.
To study the size effects of graphene, we built graphene sheets with a variety of sizes: a periodic graphene sheet, five squared H-terminated graphene sheets with a side length, is 10 Å, 20 Å, 30 Å, 40 Å, and 50 Å, respectively.The periodic GO-PEG complex was constructed for studying the effects of surface functionalization.We built graphene-drug complex models initially by "adsorption locator" module implemented in Materials Studio [33] with Dreiding force field [34,35].Dreiding force field is a rule based classical force field and has good transferability for organic systems.We load multiple drug molecules ( = even) gradually on either single side or double sides of graphene sheets.In total we studied 420 graphene-drug complex models.Some examples of these atomic models are shown in Figure 2. All model structures were relaxed by geometry optimization using Dreiding force field.QEq charge assignments were adopted.The convergence criteria 2.0 × 10 −5 kcal/mol of energy and 0.001 kcal/mol/ Å of force were used.

MD Simulations.
To investigate the dynamic process of drug molecules on the graphene, we performed molecular dynamic (MD) simulations to elucidate the absorption behavior of drug molecules.On the basis of the geometry optimization, we first carried out the simulations at the constant volume and constant temperature (NVT) ensemble ( = 298 K).A total computation time of 1 ns and a time step of 1 fs were used to integrate Newton's equation of motion and the trajectories of the structures were saved every 0.1 ps.After the dynamics reach equilibrium, we used the last 500 ps data for the following analyses.
The binding energy (  ) of drug molecules on graphene is defined as follows: where (D  @G) is the average total energy of graphenedrug system,  D is the average total energy of a single drug molecule, (G) is the average total energy of graphene sheet, and  is the number of drug molecules absorbed onto graphene.Equation (1) shows that the binding energy of graphene and drug molecules consists of the interaction energy between graphene and drug molecules  int (D  @G), the interaction energy among drug molecules  int (D  ), and the deformation energy of graphene   (G) as well as drug molecules   (D  ).Average binding energy is an integrated representation of binding strength between drug molecules and graphene.In order to further understand the binding energy, we calculate the instantaneous interaction energy, instantaneous deformation energy, and instantaneous binding energy using  the last snapshot structures of NVT simulations.These energy components are defined as follows: where  1 int (D  @G) is the instantaneous interaction energy between graphene and drug molecules,  1 (D  @G) is the single-point energy of the last snapshot of graphene-drug system,  1 (D  ) is the single-point energy of all drug molecules as a whole, and  1 (G) is the single-point energy of graphene sheet.One has where  1  (D  ) is the instantaneous interaction energy and deformation energy of drug molecules,  1 (D  ) is the singlepoint energy of all drug molecules as a whole, and  D is the average energy of a single drug molecule.These data were taken from the last frame structures in the NVT MD simulations.One has       where  1  (G) is the instantaneous deformation energy of graphene,  1  (D  ) is the instantaneous deformation energy of drug molecules, and  1 (D) is the single-point energy of a single drug molecule from the last snapshot of graphene-drug system.
After NVT simulations reach equilibrium, we used the final frame structures of NVT simulations as the initial structures to run NVE (constant volume and constant energy) MD simulations.The total simulation time was 500 ps with a 1 fs time step.The MD trajectories were printed out every 100 fs.To study the thermal diffusion behavior of drug molecules on the surface of graphene, we calculated the diffusion coefficient of drug molecules on the surface of graphene as follows: where   () and   (0) stand for the position vector at times  and 0, respectively.⟨ ⟩ refers to the fact that data are averaged over the NVE ensemble.Moreover, we decomposed the total diffusion coefficient into its components along the direction of , , and .As we mainly concern the diffusion behavior of drug molecules on the surface of graphene, we chose the averaged diffusion coefficients between  and    components parallel to the graphene sheet surface, denoted as  , as a measure of the diffusion behavior of drug molecules.After the systems reach equilibrium in the NVE simulations, the trajectories in the last 250 ps were used for analyses.

Size Effects of Graphene Sheets on Graphene-Drug Binding.
To investigate the influences of sizes of graphene sheets on the binding strength of graphene-drug complexes, we load drug molecules on the squared graphene sheets with a variety of finite side lengths, for example, 10 Å, 20 Å, 30 Å, 40 Å, and 50 Å, as well as periodic squared sheets with a unit cell 50 × 50 × 20 Å modeling infinite large graphene sheets.The four types of drug molecules, namely, CE6, DOX, MTX, and SN38, are loaded, respectively, on either one side (D  @Gf  ) or two sides [(D  ) (D  )@Gf  ] of graphene sheets, where D  represent  drug molecules and Gf  or Gp  are finite or periodic graphene sheets with a length of  Å.To focus on the graphene-drug bindings and avoid the interactions among drug molecules, we first load one drug molecule on one side of graphene sheets denoted as single-side loading mode (D 1 @Gf  ), or one drug molecule on each side of graphene sheets denoted as double-side loading mode [(D 1 ) (D 1 )@Gf  ].To evaluate the binding strength of graphene-drug complex, we discuss the results of binding energies and their components including interaction energies and deformation energies, respectively, as follows.
3.1.1.Average Binding Energy.The average binding energy   as defined in (1) represents the binding strength of graphenedrug complex with respect to the corresponding isolated components.We calculated the average binding energies   as functions of graphene sheet sizes for the four types of drug molecules (Figure 3).The results show that the average   of graphene-drug complexes vary depending on the sizes of graphene sheets.Specifically, the average   maximizes for 20-30 Å graphene sheets in both single-side and double-side loading modes.When finite graphene sheets are larger or smaller than 20-30 Å, the average binding energies become smaller and comparable to the   of periodic graphene sheets.These results suggest that the appropriate sized graphene sheets or coverage density (20-30 Å per molecule) helps stabilize the binding of graphene and drug molecules.
Comparisons of   between single-side and double-side loadings show that single-side bindings are stronger than double-side bindings for graphene sheets of size 20-30 Å, while becoming comparable for graphene sheets of size 40-50 Å.For periodic graphene sheets, the double-side binding is stronger than the single-side binding except for CE6.The binding strength between a single drug molecule and periodic graphene decreases in the order of CE6 > SN38 > DOX > MTX in the single-side loading mode, while the binding strength decreases in the order of MTX > SN38 > DOX > CE6 in the double-side loading mode.

Instantaneous Binding Energy, Interaction Energy, and
Deformation Energy.Instantaneous binding energies ( 1  ), calculated based on the last MD snapshots, show similar trends as the binding energies (  ) averaged over the whole MD trajectories; all four types of drug molecules bind most strongly with the graphene sheets of 20-30 Å (Figure 4).To understand the origin of such trends, we decompose the binding energy into its components including the interaction energy ( 1 int ) between graphene and drug as well as   the deformation energies ( 1  ) of graphene and drug, respectively (equation ( 1)).It is clearly seen that the deformation energies of graphene sheets dominate the binding energies, leading to their very similar trends.The graphene sheets of 20-30 Å exhibit the smallest deformation energies when one drug molecule is loaded possibly because the sizes of graphene sheets and drug molecules match well at this range as shown in Figure 1.On the other hand, each drug molecule deforms similarly from an energy point of view independent of the sizes of graphene sheets.Finally, the interaction energy between the graphene and drug barely changes in spite of the varied sizes and deformation of graphene sheets.Combining all analyses above, we conclude that the size match between finite graphene sheets and drug molecules is critical to determining the deformation of graphene sheets and in turn their binding strength.

Size Effects of Graphene Sheets on Diffusion of Drug
Molecules on Graphene Sheets.The in-plane diffusion coefficients,    defined in (6), measure the ease of diffusion behavior of drug molecules on the surface of graphene sheets from a kinetics point of view.We calculated the in-plane diffusion coefficients as functions of the sizes of graphene sheets when one drug molecule is loaded in either single-side or double-side modes (Figure 5).Our calculations show that the diffusion coefficients of the single molecule are negligible (<0.2 Å2 /s) when loaded on the graphene sheet of 10-20 Å    in both single-side and double-side modes, indicating that the movement of drug molecule is restricted within the small sized graphene sheets.The diffusion coefficients of drug molecules increase gradually as the graphene sheets become larger and approach towards the corresponding values on the periodic graphene sheets.This trend indicates that larger graphene sheets provide more free spaces for faster diffusion.
The diffusion coefficients of drug molecules loaded on the periodic graphene sheets in the double-side modes are always smaller than those in the single-side modes except for CE6.This means that double-side loading often slows down the diffusions of drug molecules compared with single-side loading probably due to the interactions between the drug molecules separated by the graphene sheets.The diffusion coefficients of single molecules loaded on the periodic graphene sheets increase in the order of CE6 < DOX < MTX < SN38 in the single-side mode, while the diffusion coefficient increases in the order of DOX < CE6 < MTX < SN38 in the double-side mode.

Effects of the Number and Types of Drug Molecules on
the Graphene-Drug Binding and Diffusion 3.2.1.Average Binding Energy.The average binding energy   defined in (1) for loading multiple drug molecules measures the binding strength between multiple drugs and graphene as well as those among drug molecules normalized by the number of drug molecules.The diverse configurations of    drug molecules and graphene make   fluctuate around 10-50 kcal/mol for graphene sheets of 50 Å while fluctuating around 20-70 kcal/mol for periodic graphene sheets independent of the number and types of drug molecules (Figure 6).

Instantaneous Binding Energy, Interaction Energy, and
Deformation Energy of Graphene-Drug Complexes D  @Gf 50 .Instantaneous binding energies ( 1  ) as functions of the number of drug molecules () were evaluated based on the last MD snapshots and shown in Figure 7.  1   of all four types of drugs decrease monotonically as  increases.Specifically, the instantaneous binding energies are relatively small when 1 or 2 drug molecules are loaded on the graphene sheets of 50 Å, while increasing quickly and approaching towards constant values.Since these results are obtained based on the independent systems with a variety of number and types of drug molecules, the trends of instantaneous binding energies are statistically meaningful.
To understand the trend of the instantaneous binding energies, we decompose  1  further into the interaction energy ( 1 int ) between graphene and drugs, deformation energies ( 1  ) of graphene and deformation-interaction energies of drugs ( 1  ) according to (2).It is found that the deformation energies of graphene exhibit very similar trends as the instantaneous binding energies, while the graphenedrug interaction energies and the deformation-interaction energies of drugs barely change with the number of drug  molecules.These results demonstrate that the deformation of graphene sheet determines the binding strength of graphenedrug complex when drug coverage density is low, for example, 1-2 molecules/(50 × 50 Å2 ).The larger drug coverage leads to almost invariant instantaneous deformation and binding energies, indicating that the stability of the graphene-drug complex is independent of the number of densely loaded drug molecules.The typical atomic structures of graphenedrug complexes exhibit the H-bonding and electrostatic interactions among drug molecules and vdW interactions between drug molecules and graphene sheets as shown in Figure 2.

Instantaneous Binding Energy, Interaction Energy, and
Deformation Energy of Graphene-Drug Complexes D  @Gp 50 .When multiple drug molecules are loaded onto the periodic graphene sheets, neither the instantaneous binding energies nor their components change significantly with the number of the drug molecules (Figure 8).These results are different from the case of finite graphene sheets where the deformation energies of finite graphene sheets vary depending on the number of drug molecules.Apparently, it is hard for the periodic graphene sheets to deform permanently due to binding with drug molecules besides their thermal vibrations.

Instantaneous Deformation Energies and Interaction
Energies of Drug Molecules.To examine the effects of deformation energies and interaction energies of multiple drug molecules separately, we decompose further the deformation-interaction energies ( 1  ) into the deformation ( 1  ) and interaction ( 1 int ) energies (Figure 9).The results show that the deformation energies increase slightly except for SN38, while the interaction energies become larger as the number of drug molecules increases.This means that the drug molecules deform themselves more significantly in order to interact with each other more strongly as the number of drug molecules increases, leading to the total deformationinteraction energies being barely changed.

Diffusion Coefficients of Multiple Drug Molecules.
The average diffusion coefficients of multiple drug molecules on the surface of graphene sheets were calculated as the functions of the number of the drug molecules (Figure 10).The results show that, with the increase of the number of drug molecules, the diffusion coefficients of drug molecules on the surface of graphene fall quickly in both singleside and double-side loading modes.When the number of drug molecules reaches 6-8/(50 × 50 Å2 ), the diffusion coefficients become small, meaning that the drug molecules are locked.The diffusion coefficients of drug molecules loaded on the double sides of graphene are lower than those loaded on the single side, suggesting the dragging effect of molecular motion caused by the interactions between the drug molecules separated by the graphene sheets.Moreover, the drug molecules diffuse more quickly on the surface of periodic graphene than on the finite sized graphene sheets, indicating that the limited sizes of graphene sheets slow down the movement of drug molecules probably due to the boundary effects.

Loading Drug Molecules onto PEGylated
Graphene Oxide (GO-PEG)   coefficients (Figure 11).Compared with pristine graphene sheets, the absolute values of the binding energies increase significantly when loading drug molecules onto the GO-PEG sheets.These mean that the oxidized functional groups and PEG chains on the surface of graphene have strong attractive interactions with drug molecules.When loading one drug molecule, the binding energies of drug molecule and GO-PEG complexes decrease in the order of DOX > SN38 > MTX > CE6.When the number of drug molecules increases, the binding energies decrease and the differences of binding energies among various types of drug molecules become smaller, probably because of the fact that multiple drug molecules distort the flexible PEG chains more significantly causing energy penalty from the deformation of PEG chains.
The typical atomic structures of MD snapshots of D  @GO-PEG complex are shown in Figure 12.We found that the arrangement of drug molecules is scattered, indicating that the drug molecules tend to bind with the PEG chains rather than aggregate by themselves.The in-plane diffusion coefficients (Figure 11(b)) show that, compared with the drugs loaded on the pristine graphene sheets, the diffusion coefficients of drug molecules loaded on the GO-PEG sheets decrease significantly (<0.1 Å2 /s) and remain unchanged with the increase of the number of drug molecules.This suggests that the drug molecules on the GO-PEG sheets are almost immobile.Considering both the large binding energies and small diffusion coefficients discussed above, GO-PEG-drug complexes are obviously more stable than the graphene-drug complexes.This can be understood by the fact that there are mainly vdW interactions between drug molecules and graphene in D  @G, while additional electrostatic and H-bonding interactions in D  @GO-PEG contribute to the stronger bindings.

Conclusions
In this work, molecular dynamic simulations were performed to investigate the energetic stabilities and diffusion behaviors of the graphene-drug complexes.We focus on the influences of the size of graphene sheets, the number and types of drug molecules, and the loading modes.Our simulations show that the binding strength of graphene-drug complex is mainly determined by the deformation of finite graphene sheets.When the areas that drug molecules occupy have comparable sizes as graphene sheets, for example, 20-30 Å/molecule, the deformation of graphene sheets is minimized and the graphene-drug bindings are the strongest.The average binding strength per molecule fluctuates between 10 and 70 kcal/mol that is not sensitive to the number and types of drug molecule as well as their loading modes.If the density of drug coverage is low and graphene sheets are relatively large, the binding strength is mainly determined by the interaction of graphene-drug.The limited sizes of graphene sheets restrict the movement of drug molecules.Multiple drug molecules may form clusters that slow down the diffusion on graphene sheets.Diffusion in the doubleside loading mode is often slower than that in the singleside loading mode.Compared with pristine graphene sheets, graphene oxide functionalized with PEG chains has stronger bindings with drug molecules so that the drug molecules are essentially immobilized.These results give physical insights into the stability and dynamics of graphene-drug complexes, helpful for designing novel graphene-based drug delivery systems.

Figure 3 :
Figure 3: Average binding energies (  ) of graphene-drug complexes when one drug molecule is loaded onto the graphene sheets with sizes of 10-50 Å in (a) single-side and (b) double-side modes.The dashed lines represent   when one drug molecule is loaded onto the periodic graphene sheets in (a) single-side and (b) double-side modes.

Figure 5 :
Figure 5: Average in-plane diffusion coefficients of one drug molecule loaded on the graphene sheets of 10-50 Å in (a) single-side and (b) double-side modes.The dashed lines represent the average in-plane diffusion coefficients of one drug molecule loaded on periodic graphene in (a) single-side and (b) double-side modes.The diffusion coefficients of four types of drug molecules CE6, DOX, MTX, and SN38 are shown in the figures, respectively.

Figure 6 :
Figure 6: Average binding energies (  ) of graphene-drug complexes when loading various numbers and types of drug molecules onto ((a) and (b)) the finite graphene sheets of 50 Å and ((c) and (d)) the periodic graphene sheets in ((a) and (c)) single-side and ((b) and (d)) double-side modes.

Figure 9 :
Figure 9: (a) Instantaneous deformation energies and (b) instantaneous interaction energies of different types and numbers of drug molecules loaded onto the periodic graphene sheets.

Figure 10 :
Figure 10: In-plane diffusion coefficients (   ) of various types and numbers of drug molecules loaded on ((a) and (b)) the finite graphene sheets of 50 Å and ((c) and (d)) the periodic graphene sheets in ((a) and (c)) single-side modes and ((b) and (d)) double-side modes.

(Figure 11 :
Figure 11: (a) Average binding energies and (b) diffusion coefficients of D  @(GO-PEG) complex with various types and numbers of drug molecules in single-side loading modes.

Figure 12 :
Figure12:(a) Atomic structure of GO-PEG complex.The GO-PEG is divided into four parts decorated by four types of functional groups (i) hydroxy, (ii) carboxide, (iii) epoxy group, and (iv) carboxyl, respectively.(b) Typical atomic structures of graphene oxide-drug complex CE6 8 @GO-PEG (top view and side view).