The role of polymer coatings in lipid membrane penetration by graphene oxide dots

Abstract

Understanding the cell membrane penetration process of biomedical nanosystems and its dependence on nanomaterial properties and surface functionalization is crucial for the rational design of safe and efficient cellular internalization strategies. Computer simulations are powerful tools to evaluate the thermodynamic aspects of the process and to elucidate its underlying molecular mechanisms. In this work, the interaction between uncoated or polymer-coated graphene oxide (GO) dots and lipid bilayer models is investigated by coarse-grained (CG) molecular dynamics (MD) simulations. We first validate the coarse-grained model against all-atom MD simulations (AAMD). Then, we perform CGMD simulations and free energy calculations to assess the effect of the polymeric coating and of its features (grafting density, polymer end-group charge and polymer hydrophilic/hydrophobic character) on the interaction between GO dots of realistic size and lipid membranes. We find that the membrane penetration of GO dots is spontaneous when coated with a low-density polyethylene glycol (PEG) layer, while a high-density PEG coating prevents the penetration, and a mixed PEG/polyethylene (PE) coating excessively stabilizes the nanosystem in the inner membrane region. These findings will help to fine-tune how GO dots interact with cellular membranes.

---

1. Introduction

Graphene oxide (GO) is an extremely versatile material with environmental, electrochemical, sensing and biomedical applications.^1,2 GO is a 2-dimensional (2D) material composed of a layer of carbon atoms and oxygen-containing functional groups, including hydroxyl, epoxide, and carboxyl groups. The facile functionalization, high aqueous dispersibility and biocompatibility make it an appealing material for the design of various biomedical nanosystems. Indeed, GO can be used not only as a carrier for drug delivery, but also as a substrate for regenerative medicine and tissue engineering, and as a photoactive material for bioimaging and for photodynamic therapy.^3,4

Whenever GO is administered to the human body for any biomedical application, it must overcome many biological barriers, the most common of which is the cell membrane. In fact, cell membranes serve as barriers to the biodistribution of nanomaterials and small molecules.^5,6 In the case of nanomaterials carrying cargos to be released into the cytoplasm, a spontaneous membrane penetration is preferred, which can occur either by active cellular processes or by passive diffusion.⁷ A fundamental understanding of nanomaterial properties in aqueous solution is crucial for determining how they interact both with biological entities^8–12 and with cell membranes,¹³ i.e., whether they adhere to, penetrate, or cross the membrane and potentially oxidize and damage it.¹⁴ The most crucial parameters determining the outcome of nanoparticle (NP)/membrane interactions are (i) the particle size and shape, (ii) its electrostatic nature and (iii) its surface chemistry.¹⁵ In the case of GO, its structural features – thickness, lateral dimension, and surface functionalization – mainly govern its interaction with cell membranes and ultimately its suitability for the biomedical application of choice.¹⁶ In this work, we study the passive membrane penetration process and will not pursue internalization processes that require active cell engagements such as endocytosis.¹⁷ In particular, we focus our investigation on the passive translocation of GO across the plasma membrane, which is of particular interest for biomedical applications, both for therapy and toxicity studies. Other passive interaction events include the adhesion of GO to the outer surface of the plasma membrane or the incorporation of GO in between the leaflets of the plasma membrane.¹⁶

Insights into the nature and mechanism of the interaction of GO with model cell membranes have already been obtained both by experiments and by computational modeling.^18–28 At least four kinds of interactions are believed to be responsible for GO interaction with lipid bilayers: (i) the electrostatic attraction between deprotonated and negatively-charged edge carboxyl groups of GO and positively-charged choline groups of phospholipids, (ii) the electrostatic repulsion between GO carboxyl groups and lipid negatively-charged phosphate groups,^18,19 (iii) the hydrogen bonding of GO hydroxyl groups with lipid phosphate groups (might be mediated by water molecules), and (iv) the hydrophobic interaction of GO unoxidized areas with lipid aliphatic chains.²⁰ Moreover, it was found by MD simulations that GO nanosheets insert into lipid bilayers according to an inclined orientation^21,22,27 and induce pore formation,^23,24,27,29 favoring lipids extraction^23,25,26 and lipids scrambling.²⁷ Consequently, graphene and GO nanosheets can induce the degradation of the inner and outer cell membrane leaflets and reduce cell viability.²⁶

In fact, when used for biomedical applications, GO is usually coated with biocompatible polymeric chains to improve its dispersibility in biological solutions, to provide it with stealth properties and, finally, to enable the loading and/or conjugation of functional molecules, such as drugs, dyes, or active targeting ligands.^30–32 Among the available polymeric surface coatings, polyethylene glycol (PEG) is the most widely used one for GO.^31,33–39 PEG was proved both to limit the endocytotic uptake of GO–PEG by macrophages⁴⁰ and, conversely, to stimulate an immunological response when adhering to the cell surface.²⁸ Nevertheless, the non-specific passive interaction of polymer-coated GO with cell membranes remains poorly understood. In other words, it is not clear whether and how the polymer coating affects the passive membrane penetration process of GO nanostructures. Indeed, the surface chemistry of nanomaterials is of critical importance to determining their interaction with cell membranes.⁴¹ This is particularly relevant in the case of polymeric coatings with varying densities, charges, and hydrophobicity. For instance, it was found that PEG coating essentially blocked non-specific delivery of quantum dots into cells.⁴²

To fill this gap, in this work we assess the role played by the polymer coating of GO on its penetration into model lipid membranes by coarse-grained MD (CGMD) simulations. The latter ones are typically exploited to study the penetration into membranes or the interaction of NPs of different nature with lipids,^43–51 where simulation times of the order of μs and larger membrane models than the ones affordable by all-atom MD (AAMD) simulations are needed. This is achieved at the expense of a less detailed description of the system with respect to the atomic resolution of AAMD simulations, which, however, are spatially and temporally restricted to the observation of the early stages of NPs interaction with membranes. Apart from adopting standard MARTINI bead types,^21,22,24,52 few ad hoc coarse-grained models have been recently proposed for GO.^53–57 Among these, we select the GO FF parameters proposed by Wu et al.⁵⁴ to use in combination with the polarizable version of MARTINI2 FF.^58–61 We first validate the CG GO model for membrane penetration studies. Next, we build CG models of uncoated or polymer-coated 5 nm wide GO structures, called “dots” for their sub-10 nm size, and investigate each of them in terms of interaction with phospholipid membrane models, composed of 1-palmitoyl-2-oleoyl-glycero-3-phosphocholine (POPC) and cholesterol (CHOL, 10 mol%), under physiological conditions. POPC and CHOL are chosen as they are both natural and principal components of eukaryotic plasma membranes⁶² and the selected molar ratio is commonly used in liposome preparation and in similar computational studies.¹³ We evaluate the effect of five polymeric coatings on GO, which differ in the grafting density, in the charge of their end groups, and in the hydrophobic/hydrophilic character of the polymer chains. In fact, coatings of mixed polyethylene (PE) and PEG are generally adopted for gold NPs⁶³ to prevent protein adsorption and have been shown to favor TiO₂ NPs membrane penetration by MD simulations.¹³

2. Computational details

2.1 Systems and nomenclature

First, we introduce the nomenclature adopted in this work. The different types of polymer chains used for GO coating are depicted in Fig. 1 and listed in the following:

• Three PEG chains composed of 22 repeating units (PEG₂₂), which correspond to a molecular mass of approximately 1000 Da (experimentally referred to as PEG₁₀₀₀), and terminated with a neutral methyl group (PEG₂₂-CH₃), a negatively charged carboxyl group (PEG₂₂-COO⁻) or a positively charged amino group (PEG₂₂-NH₃⁺);

• A mixed PE-PEG chain composed of 11 repeating units of PE, followed by 11 repeating units of PEG and terminated with a neutral methyl group (PE₁₁-PEG₁₁-CH₃).

Fig. 1 Structural formulas and nomenclature of the polymer chains used to coat the GO dot in this work. The protonation state corresponds to the one at physiological pH.

Each polymer chain is referred to by its chemical acronym, i.e., PEG_n or PE_n-PEG_n, where n is the number of repeating units and “–” indicates covalent bonding. The chain name is followed by its terminal group, i.e., methyl (–CH₃), deprotonated carboxyl (–COO⁻) or protonated amine (–NH₃⁺). The amino group on the left side of the polymer chains in Fig. 1 is used for covalent conjugation with GO according to the two PEGylation strategies commonly used for GO: the ring opening of basal epoxide groups or the amidation of edge carboxyl groups (we refer the reader to ref. 64 for a more detailed description). The GO dot is named as GO and, when it is covalently functionalized with the polymer chains through their amino group, the resulting nanoconjugate is referred to as GO followed by the number and the name of the polymer chains, e.g., GO-4PEG₂₂-CH₃. The lipid membrane is named as MEMB and the systems composed of uncoated or coated GO interacting with it are indicated as MEMB followed by the corresponding GO nanomaterial, e.g., MEMB/GO or MEMB/GO-4PEG₂₂-CH₃, where “/” stands for non-covalent interaction.

2.2 Building of GO coarse-grained models

All the GO and graphene CG models used in this work were built according to the GO model by Wu et al.,^54,65 which is herein referred to as “GO FF” and is compatible with the polarizable version of MARTINI2 FF (v2.2P, herein called MARTINI2P).^58–61 The GO FF mapping scheme is inspired by the standard MARTINI one and represented in Fig. S1 in the SI. Each bead represents four carbon atoms, plus all oxygen and hydrogen atoms bonded to them. The only exception is for carboxyl edge groups, whose CG bead represents only the carboxyl group, i.e., one carbon and two oxygen atoms. This mapping scheme, which retains the hexagonal structure of graphene and graphene derivatives, was inspired by the model proposed by Williams and Lísal.⁵³ However, they only define four bead types for GO, whereas GO FF model⁵⁴ includes nine bead types, which are represented in Fig. 2. The non-bonded interaction parameters for each of the nine bead types added to MARTINI2P by Wu et al. were taken from ref. 54 and can be found in the parameters file available in the SI. The bonded parameters (bonds, angles and dihedrals) were also defined according to the GO FF model, using the equilibrium and force constant values taken from ref. 54 and included in the parameters file in the SI. The topology of GO and graphene models was built in a bottom-up fashion from the corresponding AA models with the help of PyCGTOOL.⁶⁶ The parameter assignment was performed with an in-house python script.

Fig. 2 GO FF bead types from GO FF in ref. 54 used in this work.

The CG toy models of graphene and GO used in section 2.3 were also built from the corresponding AA models following the procedure described above and are represented Fig. S4 and S5.

The CG models of uncoated and polymer-coated GO dots of realistic size, used in section 2.4 and represented in Fig. 6, were prepared as described above and detailed in the following. The model of the uncoated GO dot was obtained as a flake-shaped 5 nm-wide GO structure, whose beads were assigned A, B, C, F, G, H or I type. Specifically, the number of beads for each bead type was chosen to reproduce the C/O ratio, the type and number of groups (hydroxyl, epoxide, and carboxyl groups) and the total mass of the AA random model, used in a previous work by some of us.⁶⁴ The beads were arranged in a random fashion as the aim of this work is to assess the impact of a polymeric coating, rather than the effects of oxidation heterogeneity, which requires an atomistic resolution to be properly investigated at this system size. The mass of each bead is set equal to the sum of the masses of the mapped atoms, yielding a total mass of 11 375 g mol⁻¹. Both the AA and the CG models of a random GO dot are represented in Fig. 3.

Fig. 3 AA (left) and CG (right) models of a random flake-shaped GO dot with a diameter of 5 nm, colored according to the atom type (color code on the left) or bead type (color code on the right), respectively. The AA model is the so-called random GO model in ref. 64.

Subsequently, the GO dot was coated with polymer chains, composed of only PEG units or mixed PEG and PE units, yielding five different nanoconjugates. In particular, the CG GO dot model was functionalized with: (i) 4 PEG₂₂-CH₃, (ii) 12 PEG₂₂-CH₃, (iii) 4 PEG₂₂-NH₃⁺, (iv) 4 PEG₂₂-COO⁻ or (v) 4 PE₁₁-PEG₁₁-CH₃ chains (structures in Fig. 1). All CG polymer structures, whose representation is shown in Fig. 4, were generated with polyply suite⁶⁷ and the MARTINI2 FF^{58–60,68,69} parameters were used. Each PEG or PE-PEG chain was covalently bonded to GO by setting a distance of 0.363 nm between the anchoring GO bead and the first bead of the polymer chain, which corresponds to the average value obtained from OPLS-AA AAMD simulations (Fig. S2), and applying a constraint on the distance with a force constant of 10 000 kJ mol⁻¹ nm⁻² during all the following CGMD simulations. For the edge functionalization sites, the GO anchoring bead was turned into a MARTINI P5 bead to reproduce the amide bond and its mass was changed accordingly; for the basal functionalization sites, the bead type remained unchanged, but the mass was modified accordingly. The functionalization sites were selected to be uniformly distributed across the GO surface and edges.

Fig. 4 CG bead representation of the polymer chains used to coat GO, according to MARTINI2 FF.^58–60,68 The methylene amino group on the right side of the structures belongs to the anchoring GO bead.

The A–I beads used to model GO are not part of the standard MARTINI force field, but were introduced by Wu et al.⁵⁴ and used with MARTINI2P. While most interactions with MARTINI beads were already defined by Wu et al., the cross-interaction parameters between PEG (EO bead) and GO (beads A–I) and PEG (EO bead) were missing from GO FF.⁵⁴ Therefore, they were defined by analogy with the ones between EO bead and MARTINI beads previously used to describe GO,^21,22,52,70 validated in section 2.3 and 3.1, and listed in Table S3. EO/SC4 cross-interaction parameters were used for apolar GO beads (EO/A and EO/B), EO/P1 for intermediately polar GO beads (EO/F and EO/G) and EO/P5 for more polar and charged GO beads (EO/C, EO/H and EO/I). The structure and topology files of the coarse-grained models of the nanoconjugates built in this work are available in the SI.

2.3 All-atom and coarse-grained MD simulations for validation

A set of AAMD simulations was performed to validate the CGMD models, whose computational details are found in the (section S1.1) and results and discussion in section 3.1 and in the (section S1.2).

To validate GO FF⁵⁴ for the description of graphene and GO solvation free energy in water, we compared the free energy profile for the translocation of graphene and GO toy models from water to vacuum at AA (CHARMM-AA FF) and CG (GO FF/MARTINI2P) level (details in section S1.1.1).

To validate GO FF⁵⁴ for the membrane interaction, we compared the free energy profile of graphene or GO toy models penetration through a POPC membrane computed at AA (CHARMM-AA FF) and CG (GO FF/MARTINI2P) level (details in section S1.1.2).

To evaluate the accuracy of MARTINI2 FF^58–60,68 in combination with polarizable MARTINI water⁶¹ for PEG and PE interaction with a POPC membrane, we compared the free energy profile of a PEG and a PE monomer penetration through a POPC membrane obtained at AA (CHARMM-AA FF) or CG (MARTINI v2.0, here called MARTINI2, or MARTINI2P) level (details in section S1.1.3).

To use GO FF⁵⁴ for PEG-coated GO dots, it was necessary to validate the cross-interaction parameters between GO bead types and PEG EO bead type in Table S3. Therefore, we performed a conformational and structural analysis of GO-4PEG₂₂-COO⁻ in solution at AA (OPLS-AA FF or CHARMM-AA FF) and at CG (GO FF/MARTINI2P) level (details in section S1.1.4).

2.4 Coarse-grained MD simulations

2.4.1 Polymer-coated GO in solution. Each of the polymer-coated GO dot models described in section 2.2 was solvated with the insane tool⁷¹ in a 13 × 13 × 13 nm³ cubic box filled with polarizable MARTINI water⁶¹ and positive (Na⁺) and negative (Cl⁻) ions to counterbalance the charge of the solute and mimic the physiological salt concentration of 0.15 M. After the energy minimization with the steepest descent algorithm, the system was heated to 310 K and equilibrated for 1 ns, then an NPT MD simulation was run up to 150 ns. The V-rescale⁷² thermostat with a coupling constant of 1.0 ps and the Parrinello–Rahman⁷³ barostat with a coupling constant of 12.0 ps were used to control temperature (310 K) and pressure (1 bar), respectively. The leap-frog⁷⁴ algorithm with a time step of 20.0 fs was used to integrate Newton's equations of motion. Long-range electrostatic interactions were handled with Particle Mesh Ewald (PME) method⁷⁵ with a cutoff distance of 1.2 nm and a dielectric constant of 2.5, while short-range repulsive and attractive interactions were treated by Lennard-Jones potential with a force switching function that ramps the energy smoothly between an inner cutoff of 0.9 nm and an outer cutoff of 1.2 nm. A 150 ns CGMD simulation was run for each system with the GROMACS⁷⁶ (version 2022) open-source code.

2.4.2 Lipid bilayer. The CG membrane model was designed using the CHARMM-GUI Membrane Builder module.^77,78 We prepared a 1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-choline/cholesterol (POPC/CHOL) lipid bilayer, with 10 mol% content of CHOL (2080 lipids in total, 25 × 25 nm²), solvated on both sides with polarizable MARTINI water.⁶¹ Then, standard minimization and equilibration steps were performed under semi-isotropic pressure coupling, using the MARTINI2P FF.^58–61 We used the same treatment of long-range interactions described above for nanoconjugates in water. All CGMD were performed by means of the GROMACS⁷⁶ (version 2021) open-source code.

2.4.3 Uncoated or polymer-coated GO in a lipid bilayer. Following a 100 ns equilibration, each of the uncoated and the five polymer-coated GO dots was introduced in the center of the POPC/CHOL membrane by means of the gmx insert-molecules tool, which removed the lipid molecules overlapping with the nanoconjugates (Fig. S3 in the SI). Each of the six systems was minimized with the steepest descent algorithm and equilibrated for 10 ns at 310 K in a 25 × 25 × 25 nm³ cubic box with polarizable MARTINI water, counterbalancing ions and 0.15 M NaCl. GO beads were restrained to their minimized position with a force constant of 1000 kJ mol⁻¹ nm⁻² during the 10 ns equilibration to allow membrane equilibration first, before the nanoconjugate thermal equilibration. Then, the restrains were removed and a 1.5-μs-long production phase was run with the GROMACS (version 2021) open-source code.⁷⁶ We used the same treatment of long-range interactions described above for nanoconjugates in water and the system was kept at 310 K (V-rescale thermostat⁷² with a coupling constant of 1.0 ps) and 1 bar (semi-isotropic pressure coupling with the Parrinello–Rahman barostat,⁷³ coupling constant of 12.0 ps and compressibility of 3 × 10⁻⁴ bar⁻¹). Each unbiased MD simulation of the membrane/nanoconjugate systems was replicated three times with different initial velocity conditions to ensure statistical reliability and are herein referred to as (i), (ii) and (iii). The resulting data were averaged across all the three replicas, and standard deviations were calculated using the error propagation method. The input files and the starting-point structures used for unbiased simulations are available in the SI.

2.4.4. Free energy calculations. The free energy profile associated with the lipid bilayer penetration was calculated using the z-component of the distance between GO COM and the membrane COM, d, as a collective variable, with z being the normal direction to the bilayer (Fig. 5, left). This was done for MEMB/GO, MEMB/GO-4PEG₂₂-CH₃, MEMB/GO-12PEG₂₂-CH₃ and MEMB/GO-4PE₁₁-PEG₁₁-CH₃ systems.

Fig. 5 Schematic representation of the CVs used for US method: (left) the distance d, i.e., the z-component of the distance between the GO COM and the membrane COM to represent the distance of GO from the membrane and (right) the angle θ, between GO plane and z-axis to represent GO orientation in the membrane.

The free energy profile for different GO orientations within the lipid bilayer was obtained using the angle between the GO plane and z-axis as a collective variable (Fig. 5, right). This was done for MEMB/GO and MEMB/GO-4PEG₂₂-CH₃ systems, as representative uncoated or coated GO nanoconjugates.

As it was done in the aforementioned validation steps, a two-step protocol was employed in both cases: a steered MD (SMD) simulation and a series of umbrella sampling (US) window simulations. For SMD the same starting point of the unbiased simulations of the nanoconjugate in the membrane was used.

During the membrane penetration SMD run, the nanoconjugate was pulled out from the membrane by applying a harmonic pulling potential to the z-component of the distance between their COMs (pull-geometry = direction) with a force constant of 1000 kJ mol⁻¹ nm⁻² and a pulling rate of 0.001 nm ps⁻¹. We also imposed an orientational constraint (harmonic potential with force constant of 1000 kJ mol⁻¹ rad⁻²) on the GO dot to keep it in a perpendicular orientation with respect to the plane of the lipid bilayer. During GO orientation SMD run, a harmonic potential was applied to the angle between the GO diameter initially perpendicular to membrane plane (Fig. 5) and the z-axis (pull-geometry = angle-axis) with a force constant of 10 000 kJ mol⁻¹ rad⁻² and a pulling rate of 0.001 deg ps⁻¹. We also imposed a harmonic pulling potential to the z-component of the distance between the GO and membrane COMs with a force constant of 1000 kJ mol⁻¹ nm⁻² to keep the flake in the membrane. Positional restraints (force constant of 500 kJ mol⁻¹ nm⁻²) were imposed on lipid head groups during both SMD simulations, which were performed at 310 K (V-rescale thermostat,⁷² coupling constant of 1.0 ps) and 1 bar (semi-isotropic pressure coupling with the Berendsen barostat,⁷⁹ coupling constant of 5.0 ps and compressibility of 3 × 10⁻⁴ bar⁻¹). The same simulation setup described above for the nanoconjugates in the membrane was adopted. A box size of 25 × 25 × 45 nm³ was used.

US windows’ starting-point configurations were extracted from the SMD trajectories, at each 0.2 nm or 5° increment in the distance d and in the tilting angle θ, respectively (total of 51 windows for d and 19 for θ). In window simulations, d and θ were restrained to the starting-point value by applying a harmonic potential with a force constant of 1000 kJ mol⁻¹ nm⁻² and 3500 kJ mol⁻¹ rad⁻², respectively. During the orientation window simulations, we also imposed a harmonic pulling potential to the z-component of the distance between the GO and membrane COMs with a force constant of 1000 kJ mol⁻¹ nm⁻² to keep the flake in the membrane. Technically, to control GO orientation we also constrained the diameter, which is perpendicular to the one mentioned above, to always remain parallel to the membrane plane. For US simulations, the same MD settings of SMD simulations were used except for the Parrinello–Rahman barostat⁷³ instead of the Berendsen one. Each window simulation was run for 100 ns, and only the second half (50 ns) was considered for free energy profile calculation and analyses. The input files and starting-point structures used for the US simulations using d as a CV are available in the SI. The free energy profiles with their relative error were calculated from the count histograms using the WHAM method^80,81 and the bootstrap method implemented in the gmx wham tool.⁸²

2.5 Simulation analysis

Distances were computed through the gmx distance tool from GROMACS⁷⁶ or LOOS.⁸³

The contacts with water were counted with the gmx mindist GROMACS tool⁷⁶ and a cutoff of 0.6 nm for CGMD simulations, which was chosen to be equal to the distance corresponding to the minimum of the LJ potential for GO/water beads cross-interactions. Only the central W bead (the only bead with non-null van der Waals parameters) of polarizable MARTINI water was considered for contacts. The orientation of the GO flakes with respect to the normal of the lipid membrane was calculated with an in-house python code using the MDAnalysis python library^84,85 and was computed as the angle between the GO plane (defined by the positions of three particles belonging to GO edges) and the z-axis, which is perpendicular to the membrane plane.

The polymer radius of gyration, R_g, was calculated with gmx gyrate or gmx polystat GROMACS tools⁷⁶ and is defined by

where r_k is the position of each k-th heavy atom of the polymer chain, r_mean is the position of the center of mass of the chain and N is the total number of heavy atoms.

The number density profiles were obtained with gmx rdf with a bin size of 0.02 nm. The number density profiles along z-direction were obtained with the gmx density tool from GROMACS code,⁷⁶ subdividing the simulation box into parallelepiped-shaped slabs with a z-dimension of 0.1 nm, counting the number of particles for each bin and normalizing it by the volume of the slab. 2D density maps were generated by the gmx densmap tool from the GROMACS code.⁷⁶

The lipid order parameter, P₂, for CGMD simulations was calculated with the do-order-gmx5.py script from the MARTINI FF website.⁸⁶ P₂ is defined as:

P₂ = ½·〈3 cos² θ−1〉 (2)

where θ is the angle between the direction of the bond formed by two CG beads and the bilayer normal. The angular brackets represent the molecular and temporal ensemble average over all the lipids belonging to the bilayer and over the whole production phase.

Non-bonded interaction energies were calculated with gmx energy and summing van der Waals and Coulomb energy contributions.

Mean squared displacement (MSD) was calculated with the gmx msd tool from the GROMACS code.⁷⁶ The lateral diffusion coefficient for GO on the xy plane was calculated from the linear fitting of MSD, within the time intervals where a linear dependence with time t is observed, by the Einstein relation:

MSD = 2nDt (3)

where D is the diffusion coefficient, n is the dimensionality, which is equal to 2 for lateral diffusion. Since we are interested in relative trends, the finite size effect on computed diffusion constant was not investigated.^87,88

VMD was used for all structural representations.⁸⁹ The CG times reported in this work are simulation times and no conversion factor was used.

3. Results and discussion

This section is organized as follows: in section 3.1, we validate the CG models against AAMD simulations; in section 3.2, the molecular insights coming from unbiased CGMD simulations of uncoated or polymer-coated GO dots are reported and in section 3.3 the free energy profiles of GO dots translocation across the lipid bilayer are discussed.

3.1 Coarse-grained models validation

The preliminary step of this study is the validation of CG models against AAMD simulation data. We refer the interested reader to the section S1 for a detailed description of the validation procedure. In summary, we validate:

(i) GO FF model proposed by Wu et al.⁵⁴ for solvation free energies calculations by registering a close agreement between the AA and CG free energy profiles for the translocation of small graphene and GO toy models from water to vacuum (section S1.2.1);

(ii) CG GO parameters proposed by Wu et al.⁵⁴ in combination with polarizable MARTINI2 (MARTINI2P) ones for GO/membrane interactions by registering a fair agreement between the AA and CG free energy profiles for the translocation of small graphene and GO toy models across a POPC bilayer. The CG data are further validated against experimental POPC/water partition coefficients⁹⁰ (section S1.2.2);

(iii) CG parameters proposed by Grunewald et al.⁶⁸ for PEG in combination with the polarizable MARTINI water for polymer/membrane interactions with an excellent agreement found between the AA and CG free energy profiles for the translocation of a PE and a PEG monomer across a POPC bilayer (section S1.2.3);

(iv) GO/PEG cross-interaction parameters, which were proposed in this work as not provided by Wu et al.,⁵⁴ with the structural and conformational analysis of the behavior of GO-4PEG₂₂-COO⁻ in solution (section S1.2.4).

3.2 Uncoated versus polymer-coated GO dots interaction with a POPC/CHOL bilayer

After validating the CG model against AAMD simulations data, here we investigate the interactions of GO-based nanoconjugates with a lipid bilayer composed of a binary mixture of POPC and CHOL (10 mol%), which are both natural and abundant components of eukaryotic cell membranes and commonly used for liposome synthesis. Fig. 6 comprehends the structures of all the nanosystems considered: (i) the uncoated GO dot and (ii–vi) the five polymer-coated GO dots.

Fig. 6 Top view of the uncoated GO dot model and top and side views of each of the polymer-coated GO dots, whose interaction with lipid bilayers is evaluated in this work. Color code: A beads in black, B beads in red, C beads in green, F beads in orange, G beads in sky blue, H beads in blue, I beads in magenta (see also Fig. 2 for reference), PEG beads are shown in purple, PE beads in orange and charged terminal beads, Q_d and Q_a, in blue and red, respectively. Each of the polymer-coated GO system was equilibrated in an aqueous solution 0.15 M NaCl at physiological temperature.

The uncoated GO dot model (i) was built as described in section 2.2. It is a 5 nm wide randomly oxidized GO flake bearing a C/O ratio of 2.5, which is commonly achieved by the most common GO synthesis methods,⁹¹ and 7 deprotonated carboxyl edge groups yielding a negative charge equal to −7 e. In terms of size, experimental studies have shown that smaller GO sheets are more readily taken up by cells,⁹² and since passive membrane penetration does not involve any activation pathways, it is expected only for nanometric structures. Highly oxidized, negatively charged GO sheets tend to repel each other, forming stable monomeric dispersions, as shown experimentally and theoretically.^53,93,94 Moreover, experimental evidence shows that GO size strongly affects aggregation, with nanometer-scale GO particles resulting in enhanced stability.⁹⁵ The nanoconjugates (ii–vi) were prepared by functionalizing GO with (ii) 4 PEG₂₂-CH₃, (iii) 12 PEG₂₂-CH₃, (iv) 4 PEG₂₂-NH₃⁺, (v) 4 PEG₂₂-COO⁻, or (vi) 4 PE₁₁-PEG₁₁-CH₃ chains, as detailed in section 2.2. The neutral coating of GO-4PEG₂₂-CH₃ nanoconjugate (ii) will serve as a reference system with respect to the denser coating of nanoconjugate (iii), the positively or negatively charged coating of nanoconjugates (iv) and (v), respectively, and the more hydrophobic coating of nanoconjugate (vi). The PEG content in nanoconjugates (ii) and (iv–v) correspond to 25 wt%, i.e., the PEG mass amounts to 25% of the total nanoconjugate mass, whereas for nanoconjugate (iii) the weight percentage is equal to 50 wt%. The wt.% of PEG over the total nanoconjugate mass is the often-reported quantity to express PEG content in experimental papers since it is directly assessed by Thermogravimetric Analysis (TGA).^{34,38,39,96–98}

A preliminary conformational characterization of the polymer chains, when covalently bonded to the GO dot, is reported in Table S5. According to both R_g and 〈h²〉^1/2 values, either the inclusion of a charged terminal group or the substitution of half of the PEG monomers with the more hydrophobic PE ones leads to the expansion of the polymeric chains. Instead, we do not register any effect of a higher coating density going from a PEG content of 25 wt% to that of 50 wt%. Moreover, we notice from the representations in Fig. 7, 8 and 9 that, due to the lack of atomic resolution in the CG model of GO, and the inability to account for carbon atom hybridization, our GO dot tends to favor flatter geometries relative to atomistic models⁶⁴ (see Fig. S10). This is a limitation of this CG model, and similar approximations have been adopted in earlier studies.^21,22

Fig. 7 Free energy profiles of GO reorientation in a POPC/CHOL (10 mol%) bilayer. GO z-alignment corresponds to the CV θ in Fig. 5. The standard deviations calculated with bootstrap method are shown with shadows. Color code: GO beads in red, POPC choline ones in blue, POPC phosphate ones in tan, POPC glycerol in pink, POPC tails in cyan and CHOL in fuchsia.

Fig. 8 Snapshots of the equilibrated systems at 1.5 μs (first column), number density profiles along the direction perpendicular to the membrane (second column), averaged over the last 0.5 μs of the 1.5 μs simulation of replica (i), for MEMB/GO, MEMB/GO-4PEG₂₂-CH₃ and MEMB/GO-12PEG₂₂-CH₃. Color code: GO beads in red, choline ones in blue, phosphate ones in tan, PEG ones in purple, and water in light blue. Red color is used for all GO beads and does not refer to the bead type. A surface representation is used for water beads within 0.6 nm of GO. Only POPC head groups are shown, water and ions are not shown for clarity. See Fig. S12 for replicas (ii) and (iii).

Fig. 9 Snapshots of the equilibrated systems at 1.5 μs (first column), number density profiles along the direction perpendicular to the membrane (second column), averaged over the last 0.5 μs of the 1.5 μs simulation of replica (i), for MEMB/GO-4PEG₂₂-NH₃⁺, MEMB/GO-4PEG₂₂-COO⁻ and MEMB/GO-4PE₁₁-PEG₁₁-CH₃. Color code: GO beads in red, choline ones in blue, phosphate ones in tan, PEG ones in purple, and water in light blue. Red color is used for all GO beads and does not refer to the bead type. A surface representation is used for water beads within 0.6 nm of GO. Only POPC head groups are shown, water and ions are not shown for clarity. See Fig. S13 for replicas (ii) and (iii).

Before moving to the focus of this section, it is useful for later discussion to briefly describe the structural properties of POPC/CHOL (10 mol%) lipid bilayers.⁹⁹ From the number density profiles of the different component of the equilibrated POPC/CHOL bilayer in Fig. S11A, it is possible to identify different regions: the most hydrophilic POPC head groups region, which is centered at a z-distance of about 2 nm from the membrane center; the intermediate glycerol region centered at 1.6 nm, where also CHOL hydroxyl groups lie at about 1.2 nm; the most hydrophobic central region, where lipid tails of the two opposing monolayers are found. Moreover, going from the bulk water region towards the membrane center, the water density begins to decrease at 2.5 nm and finally drops to 0 around 1.0 nm, which indicates that the head groups are solvated by water molecules, while the inner region is not. The distance between the peaks of POPC phosphate groups, which is a measure of membrane thickness, equals to ∼4 nm, as expected for a POPC/CHOL (10 mol%) lipid bilayer.^100,101 Also, the number density map in Fig. S11B shows a uniform lipid distribution in the xy plane, i.e., the membrane plane, and will be taken as a reference for the discussion on the influence of the nanoconjugates on membrane properties in section 3.2.1.

After having equilibrated the nanoconjugates and the membrane separately, each nanoconjugate in Fig. 6 was inserted into the lipid bilayer with the GO plane perpendicular to the membrane plane, as detailed in section 2.4.3. The perpendicular orientation was chosen on the basis of several considerations: (i) it represents the intermediate step of the membrane translocation process, which is the phenomenon of interest, (ii) it was shown that the starting configuration does not prevent the system from reaching the equilibrium orientation,²⁵ and (iii) this equilibrium orientation for small GO flakes in a lipid bilayer lies around a perpendicular one, according to both previous²⁷ and current free energy calculations.

Indeed, we computed the free energy profile for GO reorientation in a lipid bilayer, while constraining it to lie inside the membrane, both for MEMB/GO and for MEMB/GO-PEG₂₂-CH₃ systems, as representative systems for an uncoated or a coated GO (Fig. 7). To achieve this, we used the US method as detailed in section 2.4 and the angle θ as a collective variable, describing the GO alignment with the normal direction to the membrane (z), as depicted in Fig. 5. For both the systems, the perpendicular orientation of GO with respect to the membrane plane is the most stable one, with a free energy difference with the parallel orientation of ∼60 kcal mol⁻¹ and ∼30 kcal mol⁻¹ for GO and GO-PEG₂₂-CH₃, respectively. Therefore, we started all our CGMD simulations of a GO dot located in the membrane from the perpendicular orientation.

Three 1.5-μs CGMD simulations that are referred to as (i), (ii) and (iii) were run for each of the models in Fig. 6, as detailed in section 2.4. The final snapshot from replica (i) of each CGMD simulation is shown in Fig. 8 and 9.

Both from the visual inspection of the last snapshots of CGMD trajectories (Fig. 8 and 9) and from the time evolution of the distance of GO COM from the membrane center in Fig. S14 and S15 it is evident that all nanoconjugates remain inside the membrane for the entire simulation time, except for GO-12PEG₂₂-CH₃, which bears the most dense PEG coating and leaves the membrane shortly after the simulation begins in all replicas (100–500 ns). On the contrary, uncoated GO and the other four nanoconjugates remain in the inner or head groups membrane region throughout the entire CGMD simulation. By comparing the number density profiles in Fig. 8, 9, S12 and S13 with the one of the isolated POPC/CHOL bilayer in Fig. S11, we do not observe any effect due to the presence of the nanoconjugate on the membrane as a whole, likely because the membrane model size (25 × 25 nm²) is substantially larger than the GO dimension (5 nm). However, as it can be noted in Fig. 8 and 9, in all the cases a notable amount of water penetrates the membrane forming a water pore, reaching even its most hydrophobic central region; moreover, some lipid head groups (blue and tan beads) are also relocated to the hydrophobic region of the bilayer.

We refer to section 3.2.1 for the effects of the nanoconjugates on membrane properties and integrity and to section 3.2.2 for a thorough discussion about the impact of the different polymer chains on GO interaction with the lipid bilayer.

3.2.1 Lipids ordering effects and water pore formation. Even though it is expected that the presence of a nanoconjugate may locally alter membrane properties, the type and extent of perturbation is not easily predictable. Lipids orientational ordering is evaluated by computing the lipids second-rank order parameter, P₂, which is defined in section 2.5. P₂ is a measure of the motional anisotropy of a particular bond and yields its average orientation.¹⁰² A complete alignment with the bilayer normal, i.e., perfect order, corresponds to P₂ = 1, a tendency of the bonds to align perpendicularly to the normal to the membrane to P₂ = −0.5, and a random lipid orientation to P₂ = 0.

Indeed, as shown in Fig. S16, the order parameter of the lipids closest to the nanoconjugates is consistently lower for both POPC tails when the nonconjugate is inserted into the bilayer than the value computed for the isolated membrane. We note a pronounced reduction in the order parameter for uncoated GO, followed by the PE-PEG and PEG coatings and, lastly, we do not observe any reduction in the case of the high-density PEG coating, as expected since the nanoconjugate is soon expelled from the membrane. The reduced ordering agrees well with the lipid re-orientation that we noticed around GO nanoconjugates in the previous section, and suggests that the lower coating densities and increased coating hydrophobicity induce greater lipid disorder. Nevertheless, given the large size of the membrane model with respect to the nanosystem dimensions, the effect registered in Fig. S16 is rather local and vanishes when averaging over all POPC lipids. Some perturbation might also affect the translational lipids order, on the membrane plane. However, we do not observe any increase in the lipids number density around any of the nanoconjugates nor a clustering of cholesterol (Fig. S17). We do note a decrease in the lipid lateral diffusion coefficients in Table S7 for both POPC and CHOL, when the membrane is perturbed by the nanoconjugates.

Regarding the penetration of water into the membrane, a water pore is formed in all the cases, independently of the presence and from the features of the polymer coating. In fact, in Fig. 10 we observe a non-zero and continuous value for the water density even in the central bilayer region. Water pore formation was reported both experimentally and computationally for uncoated GO nanostructures.^{23–25,27,29} Here, we report, to the best of our knowledge, the first observation of a pore formation due to polymer-coated GO dots. In particular, the largest water density at the membrane center is registered for the two GO flakes conjugated to PEG chains terminated with a charged group, especially COO⁻, whose negative charge adds up to GO dot negative charge. In the case of PE-PEG-CH₃ chains the lowest value is observed among coated GO dots, but still higher than the uncoated GO, since the polymer chains introduce a considerable degree of disorder in the membrane, as can be seen in Fig. S16 and S17, which may facilitate water access to the membrane inner region.

Fig. 10 Water number density for different systems. This is a focus of number density profiles of Fig. 8 and 9.

It is important to note that the present simulations rely on fixed-charge models, which do not account for possible changes in the protonation state of ionizable groups, e.g. carboxyl and amino groups, within the membrane. If the (de)protonation can be considered instantaneous, such changes in the protonation state could influence the formation of transient water pores. Nevertheless, a more accurate treatment would require constant pH simulations or other approaches capable of capturing also the protonation kinetics, which are beyond the scope of the present work.

3.2.2 Effect of polymer coating on GO interaction with a POPC/CHOL bilayer. Here, we analyze the effect of coating GO with polymers on its interaction with zwitterionic POPC membranes. In particular, using the 25 wt% PEG coating (GO-4PEG₂₂-CH₃) as a reference, we evaluate the influence of the following features of the coating chains: the terminal group charge and the hydrophobic/hydrophilic character. The effect of density will not be discussed here because it was seen that an increase in the PEG density makes GO leave the membrane.

First, from GO, PEG and PE number density profiles in Fig. 8, 9, S12 and S13 we notice that the PEG chains in GO-PEG nanoconjugates tend to escape from the membrane central region and locate at a z-distance between 2 and 4 nm from the bilayer center, whereas the PE portion of PE-PEG chains in GO-4PE₁₁-PEG₁₁-CH₃ is more likely to be found in the central region of the membrane because of its hydrophobic character. Moreover, we notice that PEG chains can assume either a symmetric (for MEMB/GO-4PEG₂₂-NH₃⁺ and MEMB/GO-4PEG₂₂-COO⁻ in Fig. 9) or an asymmetric distribution (for GO-4PEG₂₂-CH₃ in Fig. 8) around GO to maximize their interaction with the POPC hydrophilic head groups and the surrounding aqueous solution. Interestingly, in the presence of PEG or PE-PEG chains, GO adopts a slightly oblique orientation with respect to the membrane normal (10–20°), while it remains perpendicular to the membrane plane when uncoated (less than 5°) (Fig. S14, S15 and Table S6). Moreover, GO coating may influence its lateral diffusion in the membrane, as reflected by the computed lateral diffusion coefficients shown in Table S7. Indeed, the PEG coating slows down GO diffusion by halving the diffusion coefficient, independently of the terminal group charge. When PE-PEG chains are used, the interaction of the nanoconjugate with the membrane lipids, which is mainly driven by PE, is so strong that we did not observe a linear dependence of MSD on time and therefore could not report the diffusion coefficient within this time scale.

Lastly, by counting the number of contacts between the GO dot and water particles in Fig. S14 and S15 (average in Table S6), we find that GO is always, at least partially, solvated. Indeed, even when the GO dot is completely inserted in the lipid bilayer, the GO edges are located at the region of the highly solvated head groups, since the GO diameter (5 nm) is commensurate with the membrane thickness (∼4 nm). Nevertheless, water beads not only solvate the GO edges but also penetrate into the membrane (Fig. 8 and 9), leading to the eventual formation of a water pore, which was investigated in section 3.2.1.

In summary, the presence of a charged terminal group leads to a symmetric distribution of the PEG chains and to the formation of a larger water pore, as discussed in the previous section. An increase in the hydrophobic character of the polymer results in decreased mobility and solvation of the nanoconjugate in the lipid membrane, because of the increased interaction of the polymer chains with the inner region of the membrane, which is expected to help stabilizing the nanoconjugate in the membrane.

We next evaluate the effect of the GO coating and of the density and of the hydrophobic/hydrophilic character of the polymer coating on the membrane penetration process from a thermodynamic point of view in section 3.3 by free energy calculations.

3.3 Thermodynamics of uncoated versus polymer-coated GO dots penetration into a POPC/CHOL bilayer

Here, we calculate the Gibbs free energy profiles for membrane penetration of GO, GO-4PEG₂₂-CH₃, GO-12PEG₂₂-CH₃ and GO-4PE₁₁-PEG₁₁-CH₃ by US method. We have selected these four nanosystems, excluding the PEG coatings with positively and negatively charged end groups, based on our observations in unbiased CGMD simulations. Indeed, the behavior of GO-4PEG₂₂-NH₃⁺ and GO-4PEG₂₂-COO⁻ was similar to the correspondent neutral coating. Moreover, due to the presence of deprotonated carboxyl groups and polymer coatings, which are known to reduce the tendency of GO to aggregate,^36,94,103 we focus here on the membrane penetration behavior of individual GO nanoconjugates (see section S2 in the SI for a preliminary discussion on GO aggregation at AA level).

As detailed in section 2.4, we first performed SMD, during which the flake was constrained at a perpendicular orientation with respect to the membrane plane, because it is the most stable one for GO in the membrane as revealed in the free energy calculations presented above and it is also the strongly preferred one for membrane approaching, according to entropic considerations¹⁰⁴ and previous studies.¹⁰⁵ The distance d and force time evolution are shown in Fig. S19, where we notice that the lowest force is required to pull MEMB/GO-12PEG₂₂-CH₃ out from the membrane, as expected from its behavior in the CGMD unbiased simulations discussed in the previous section.

The above-mentioned orientational constraint was then removed for window US simulations to let the flake assume the most favorable orientation at every penetration step. Nevertheless, consistently with our discussion in section 3.2, we observe that GO flakes maintain a quasi-perpendicular orientation with respect to the bilayer plane and gradually deviate from it until they reach the solution state, where they explore the entire orientational space. We ran 51 window simulations for each of the four systems and calculated the corresponding free energy profile from the count histograms in Fig. S20. The computed free energy profiles for the bilayer penetration of GO, GO-4PEG₂₂-CH₃, GO-12PEG₂₂-CH₃ and GO-4PE₁₁-PEG₁₁-CH₃ are shown in Fig. 11. Table 1 lists the predicted values of the free energy of penetration, ΔG_p, which is defined as:

ΔG_p = G_m(d) − G_w(d), (4)

where G_m and G_w are the absolute Gibbs free energies of the nanoconjugates in the center of the membrane and in the bulk water, respectively. A positive ΔG_p implies that the nanoconjugate is more stable in bulk water than in the membrane, thus preventing spontaneous translocation across the lipid membrane. Conversely, a negative ΔG_p indicates that the nanoconjugate exhibits enhanced thermodynamic stability when embedded within the membrane with respect to its solution state. A relatively low magnitude of ΔG_p is preferable for the nanoconjugate to spontaneously cross the membrane without being trapped within it. It is important to note that kinetic barriers associated with passive membrane translocation are not considered in this analysis.

Fig. 11 (Top) Gibbs free energy profiles of the membrane penetration process of GO, GO-4PEG₂₂-CH₃, GO-12PEG₂₂-CH₃ and GO-4PE₁₁-PEG₁₁-CH₃, with respect to the distance between GO and membrane COMs along z, i.e., CV d in Fig. 5. The standard deviations are calculated with bootstrap method and are shown with shadows. (Bottom) Representative figures of intermediate penetration stages. Color code: GO in red, POPC choline groups in blue, POPC phosphate ones in tan, POPC glycerol ones in pink, POPC tails in cyan and cholesterol beads in fuchsia.

Table 1 Free energy of penetration, ΔG_p, calculated as the difference between the Gibbs free energy value for the flake inside the lipid bilayer and the one for the flake in water. The standard deviation is calculated with bootstrap method

System	ΔGp [kcal mol^−1]
MEMB/GO	+4 (±1)
MEMB/GO-4PEG22-CH3	–13 (±1)
MEMB/GO-12PEG22-CH3	+40 (±1)
MEMB/GO-4PE11-PEG11-CH3	–26 (±1)

---

We note that the relative thermodynamic stability of the nanosystems in the membrane, from the least to the most stable, follows the order: GO-12PEG₂₂-CH₃ < GO < GO-4PEG₂₂-CH₃ < GO-4PE₁₁-PEG₁₁-CH₃. Even though for the uncoated GO, which bears an overall negative charge, one could expect a rather high energy cost for the membrane translocation due to its hydrophilic character, at least two considerations can explain the small value of +4 (±1) kcal mol⁻¹: GO charged edge groups are mainly located near the hydrophilic lipid head groups, since GO diameter is slightly larger than the membrane thickness, and POPC is a zwitterionic molecule thus, not completely repulsive to GO. Moreover, it must be considered that GO induces pore formation and that the lipid order is locally disturbed, so that even when GO COM is located within the inner membrane region, it is not surrounded by the highly hydrophobic lipid tail environment. Interestingly, for pure PEG coatings, we find that the stability of the nanoconjugate in the inner region of the membrane first increases by coating GO with a small number of PEG chains and then decreases dramatically by using higher PEG densities. The reasons for this non-monotonic trend in free energy values (GO-12PEG₂₂-CH₃ < GO < GO-4PEG₂₂-CH₃, from low to high stability in the core region of the bilayer) are rationalized in terms of the electrostatic and dispersion interactions between the components of the system (Table 2). Both at low and high densities, the PEG chains reach a similar conformation in terms of radius of gyration and end-to-end distance (see Table S8) and exhibit a similar non-bonded interactions pattern in Table 2. Nevertheless, we note that they face a stronger repulsive interaction against the GO flake with the denser coating relative to the low-density case. Moreover, the high-density coating shields the GO flake from forming favorable interactions with the lipids (Table 2). Finally, for the less dense coating, the conformational entropy cost for the four chains to adopt their equilibrium conformation is likely sufficiently small to be overcome by the favorable interactions between the nanoconjugate and the membrane. In contrast, for the denser PEG coating, wherein a three-times higher number of chains is conjugated to GO, there is a substantially higher conformational entropy cost for confining the chains within the membrane.

Table 2 Non-bonded interaction energies between pairs of system components with their relative standard deviation. The values are averaged over the last 50 ns of the window simulation centered on d = 0 nm, where d is CV in Fig. 5. The values for polymer chains are normalized by the number of chains. The polymer contribution is split into PE and PEG contributions for GO-4PE₁₁-PEG₁₁-CH₃

Energy [kcal mol^−1]	GO/membrane	Polymer/membrane	Polymer/GO	Polymer/polymer	Polymer/solution
GO	–745 (±7)	—	—	—	—
GO-4PEG22-CH3	–524 (±23)	–3.2 (±0.2)	48 (±1)	–5.51 (±0.01)	–440 (±2)
GO-12PEG22-CH3	–185 (±2)	–1.4 (±0.2)	54.3 (±0.1)	–5.69 (±0.01)	–451.5 (±0.2)
GO-4PE11-PEG11-CH3	–728 (±5)	PE: −66 (±1)	PE: 28.0 (±0.1)	PE: −0.4 (±0.1)	PE: −5.1 (±0.2)
		PEG: −1.02 (±0.07)	PEG: −0.03 (±0.01)	PEG: −2.470 (±0.008)	PEG: −226.2 (±0.2)
				PE/PEG: −0.21 (±0.01)

---

Lastly, as expected also from previous calculations,¹³ the membrane penetration of GO-4PE₁₁-PEG₁₁-CH₃ is a spontaneous process because of the hydrophobic character of PE portion of the polymer chains and the possibility for the remaining PEG part to still interact with the lipids head groups and water. These findings are in line with the unbiased MD simulations results in section 3.2, in which we observed that GO-12PEG₂₂-CH₃ was the only nanoconjugate to spontaneously leave the membrane. Moreover, Li et al.¹⁰⁶ performed CGMD simulations for model nanoparticles and found that the hydrophobic ones are thermodynamically stable around the middle of the hydrophobic region of the membrane, in contrast to the semi-hydrophilic ones’ preference for the adsorption on the surface of the bilayer.

Taken together, the thermodynamic data in this section highlight that not only the partial hydrophobic character of the PE-PEG coating but, surprisingly, also the low-density (25 wt%) PEG coating, have a stabilizing effect for GO in a lipid bilayer.

4. Conclusions

In this work, we have studied the interaction between negatively charged randomly oxidized GO dots, either uncoated or conjugated with polymer chains of different chemical composition or grafting density, and phospholipid membranes by coarse-grained MD simulations.

First, we validated the GO FF/MARTINI2P CG model for GO/membrane interaction, polymer/membrane interaction and PEGylated GO modeling, by systematic comparison against all-atom MD simulations and experimental data. We found a fair agreement between AA and CG data in terms of solvation free energy, membrane partitioning and PEG structural and conformational analysis. Then, from a preliminary free energy calculation, we found that the perpendicular orientation to the membrane plane is the most favorable one for both GO and PEGylated GO. Therefore, we ran replicated unbiased CGMD simulations of the uncoated GO dot and of the polymer-coated one perpendicularly inserted in a POPC/CHOL (10 mol%) zwitterionic bilayer. By changing the features of the polymer coating we aimed at assessing the effect of (i) the PEG density (25 wt% against 50 wt%), (ii) the terminal group charge (–CH₃, –NH₃⁺ or –COO⁻) and (iii) the hydrophobic/hydrophilic character of the polymer (PEG₂₂-CH₃ against PE₁₁-PEG₁₁-CH₃). By unbiased CGMD simulations we found that at the highest PEG density the nanoconjugate escapes from the inner region of the membrane towards the bulk water solution. The terminal group charge has little effects as long as the limited number of attached PEG chains allows them to escape from the hydrophobic part of the membrane. On the contrary, when the hydrophobic nature of the polymer increases going from pure PEG to mixed PE-PEG chains, the PE portion is stable in the inner part of the lipid bilayer.

Finally, we used umbrella sampling to estimate the Gibbs free energy differences associated with the membrane penetration process by the nanoconjugates, as a measure of its thermodynamic spontaneity. In agreement with unbiased CGMD simulations, the membrane translocation of the most densely (50 wt%) PEG-coated GO is not spontaneous, but associated with a free energy cost, while the PE-PEG-coated GO is stable within the inner region of the bilayer. Unexpectedly, the low-density (25 wt%) PEG coating favors membrane penetration, because of the low conformational entropy cost for the PEG chains to accommodate in the most hydrophilic membrane regions and in the surrounding aqueous phase.

Taken together, the inclusion of the amount of PE tested in this work stabilizes the nanosystem in the membrane inner region and prevents its translocation. In this respect, chains composed of a shorter PE portion and a longer PEG one (at least 40% PE and 60% PEG) might be preferable if the goal is to design a nanosystem capable of crossing the lipid membrane without being trapped inside. More interestingly, the PEG content (wt.%) plays a critical role in determining the spontaneity of the membrane penetration of coated GO.

In conclusion, coarse-grained simulations provide valuable insights into the membrane penetration process of different graphene-based nanoconjugates but come with inherent limitations. The loss of atomic-level resolution hinders the understanding of fine structural details. Additionally, the reduction in degrees of freedom smooths free energy landscapes, and the acceleration of kinetics complicates the interpretation of dynamic behavior. Despite these drawbacks, coarse-grained models remain essential tools for evaluating the relative stability of nanoconjugates within membrane environments versus the bulk water phase at a reasonable computational cost.

These findings are helpful to the experimental community as they suggest that the spontaneity of coated GO nanosystems in the lipid membrane can be fine-tuned by adjusting the PEG content of the attached polymer chains, e.g., by increasing their grafting density. Ultimately, it is important to understand how to optimize the surface coating of nanoconjugates to meet the treatment and toxicity requirements of the biomedical application of interest.

Conflicts of interest

The authors declare no conflict of interest.

Data availability

The data supporting this article have been included as part of the SI.

Supplementary Information file contains a detailed description of coarse-grained model validation, supplementary figures and tables, and a brief discussion on graphene oxide aggregation. Supplementary data include starting-point structures and input files for the coarse grained simulations described in the main text. See DOI: https://doi.org/10.1039/d5nr00838g

Acknowledgements

The authors are grateful to Lorenzo Ferraro for his technical support and to Zeke A. Piskulich and Reilly Osadchey Brown from Boston University for useful discussions. The research leading to these results has received funding from the European Union – NextGenerationEU through the Italian Ministry of University and Research under PNRR – M4C2-I1.3 Project PE_00000019 “HEAL ITALIA” to Prof. Cristiana Di Valentin CUP H43C22000830006 of the University of Milano Bicocca. GF and QC acknowledge the support by the NSF Center for Sustainable Nanotechnology under Grant No. CHE-2001611 for partially supporting the visit of GF to the Cui group.

References

1. D. Chen, H. Feng and J. Li, Chem. Rev., 2012, 112, 6027–6053.
2. N. Panwar, A. M. Soehartono, K. K. Chan, S. Zeng, G. Xu, J. Qu, P. Coquet, K.-T. Yong and X. Chen, Chem. Rev., 2019, 119, 9559–9656.
3. N. Karki, H. Tiwari, C. Tewari, A. Rana, N. Pandey, S. Basak and N. G. Sahoo, J. Mater. Chem. B, 2020, 8, 8116–8148.
4. S. S. An, S.-Y. Wu and J. Hulme, Int. J. Nanomed., 2015, 10, 9.
5. S. Wilhelm, A. J. Tavares, Q. Dai, S. Ohta, J. Audet, H. F. Dvorak and W. C. W. Chan, Nat. Rev. Mater., 2016, 1, 16014.
6. P. Siani, E. Donadoni, L. Ferraro, F. Re and C. Di Valentin, Biochim. Biophys. Acta, Biomembr., 2022, 1864, 183763.
7. H. Meng, W. Leong, K. W. Leong, C. Chen and Y. Zhao, Biomaterials, 2018, 174, 41–53.
8. H. Liu, P. Siani, E. Bianchetti, J. Zhao and C. Di Valentin, Nanoscale, 2021, 13, 9293–9302.
9. P. Siani, G. Frigerio, E. Donadoni and C. Di Valentin, J. Phys. Chem. C, 2023, 127, 9236–9247.
10. S. Motta, P. Siani, A. Levy and C. Di Valentin, Nanoscale, 2021, 13, 13000–13013.
11. S. Motta, P. Siani, E. Donadoni, G. Frigerio, L. Bonati and C. Di Valentin, Nanoscale, 2023, 15, 7909–7919.
12. G. Frigerio, E. Donadoni, P. Siani, J. Vertemara, S. Motta, L. Bonati, L. De Gioia and C. Di Valentin, Nanoscale, 2024, 16, 4063–4081.
13. E. Donadoni, P. Siani, G. Frigerio, C. Milani, Q. Cui and C. Di Valentin, Nanoscale, 2024, 16, 9108–9122.
14. R. M. De Souza, P. Siani, T. F. Schmidt, R. Itri and L. G. Dias, J. Phys. Chem. B, 2017, 121, 8512–8522.
15. X. Zhang, G. Ma and W. Wei, NPG Asia Mater., 2021, 13, 52.
16. K. Kostarelos and K. S. Novoselov, Science, 2014, 344, 261–263.
17. Z. Tu, G. Guday, M. Adeli and R. Haag, Adv. Mater., 2018, 30, 1706709.
18. R. Frost, G. E. Jönsson, D. Chakarov, S. Svedhem and B. Kasemo, Nano Lett., 2012, 12, 3356–3362.
19. S. Li, A. J. Stein, A. Kruger and R. M. Leblanc, J. Phys. Chem. C, 2013, 117, 16150–16158.
20. L. Wu, L. Zeng and X. Jiang, J. Am. Chem. Soc., 2015, 137, 10052–10055.
21. Z. Chen, X. Lu, J. Liu, X. Tian, W. Li, K. Yang and B. Yuan, J. Phys. Chem. B, 2021, 125, 3589–3597.
22. J. Wang, Y. Wei, X. Shi and H. Gao, RSC Adv., 2013, 3, 15776.
23. X. Liu and K. L. Chen, Langmuir, 2015, 31, 12076–12086.
24. Q. Hu, B. Jiao, X. Shi, R. P. Valle, Y. Y. Zuo and G. Hu, Nanoscale, 2015, 7, 18025–18029.
25. J. Chen, G. Zhou, L. Chen, Y. Wang, X. Wang and S. Zeng, J. Phys. Chem. C, 2016, 120, 6225–6231.
26. Y. Tu, M. Lv, P. Xiu, T. Huynh, M. Zhang, M. Castelli, Z. Liu, Q. Huang, C. Fan, H. Fang and R. Zhou, Nat. Nanotechnol., 2013, 8, 594–601.
27. X. Zhu, C. Huang, N. Li, X. Ma, Z. Li and J. Fan, J. Colloid Interface Sci., 2023, 637, 112–122.
28. N. Luo, J. K. Weber, S. Wang, B. Luan, H. Yue, X. Xi, J. Du, Z. Yang, W. Wei, R. Zhou and G. Ma, Nat. Commun., 2017, 8, 14537.
29. P. Chen, H. Yue, X. Zhai, Z. Huang, G.-H. Ma, W. Wei and L.-T. Yan, Sci. Adv., 2019, 5, eaaw3192.
30. G. Duan, S. Kang, X. Tian, J. A. Garate, L. Zhao, C. Ge and R. Zhou, Nanoscale, 2015, 7, 15214–15224.
31. S. F. Kiew, L. V. Kiew, H. B. Lee, T. Imae and L. Y. Chung, J. Controlled Release, 2016, 226, 217–228.
32. S. Ghosh and K. Chatterjee, Int J Nanomedicine, 2020, 15, 5991–6006.
33. O. Akhavan and E. Ghaderi, Small, 2013, 9, 3593–3601.
34. E. Bidram, A. Sulistio, H. J. Cho, A. Amini, T. Harris, A. Zarrabi, G. Qiao, A. Stewart and D. E. Dunstan, ACS Appl. Mater. Interfaces, 2018, 10, 43523–43532.
35. X. Pei, Z. Zhu, Z. Gan, J. Chen, X. Zhang, X. Cheng, Q. Wan and J. Wang, Sci. Rep., 2020, 10, 1–15.
36. X. Sun, Z. Liu, K. Welsher, J. T. Robinson, A. Goodwin, S. Zaric and H. Dai, Nano Res., 2008, 1, 203–212.
37. W. Miao, G. Shim, S. Lee, S. Lee, Y. S. Choe and Y. K. Oh, Biomaterials, 2013, 34, 3402–3410.
38. Z. Xu, S. Wang, Y. Li, M. Wang, P. Shi and X. Huang, ACS Appl. Mater. Interfaces, 2014, 6, 17268–17276.
39. Z. Xu, S. Zhu, M. Wang, Y. Li, P. Shi and X. Huang, ACS Appl. Mater. Interfaces, 2015, 7, 1355–1363.
40. N. Luo, D. Ni, H. Yue, W. Wei and G. Ma, ACS Appl. Mater. Interfaces, 2015, 7, 5239–5247.
41. A. M. Alkilany, L. Zhu, H. Weller, A. Mews, W. J. Parak, M. Barz and N. Feliu, Adv. Drug Delivery Rev., 2019, 143, 22–36.
42. Y. Xiao, S. P. Forry, X. Gao, R. D. Holbrook, W. G. Telford and A. Tona, J. Nanobiotechnol., 2010, 8, 13.
43. L. Monticelli, J. Chem. Theory Comput., 2012, 8, 1370–1378.
44. F. Simonelli, D. Bochicchio, R. Ferrando and G. Rossi, J. Phys. Chem. Lett., 2015, 6, 3175–3179.
45. S. Salassi, F. Simonelli, D. Bochicchio, R. Ferrando and G. Rossi, J. Phys. Chem. C, 2017, 121, 10927–10935.
46. S. Salassi, E. Canepa, R. Ferrando and G. Rossi, RSC Adv., 2019, 9, 13992–13997.
47. M. Das, U. Dahal, O. Mesele, D. Liang and Q. Cui, J. Phys. Chem. B, 2019, 123, 10547–10561.
48. T. Lunnoo, J. Assawakhajornsak and T. Puangmali, J. Phys. Chem. C, 2019, 123, 3801–3810.
49. T. Lunnoo, J. Assawakhajornsak, S. Ruangchai and T. Puangmali, J. Phys. Chem. B, 2020, 124, 1898–1908.
50. E. Canepa, D. Bochicchio, M. Gasbarri, D. Odino, C. Canale, R. Ferrando, F. Canepa, F. Stellacci, G. Rossi, S. Dante and A. Relini, J. Phys. Chem. Lett., 2021, 12, 8583–8590.
51. M. Ivanov and A. P. Lyubartsev, J. Comput. Chem., 2024, 45, 1364–1379.
52. N. Willems, A. Urtizberea, A. F. Verre, M. Iliut, M. Lelimousin, M. Hirtz, A. Vijayaraghavan and M. S. P. Sansom, ACS Nano, 2017, 11, 1613–1625.
53. C. D. Williams and M. Lísal, 2D Mater., 2020, 7, 025025.
54. Y. Wu, J. Yang, A. van Teijlingen, A. Berardo, I. Corridori, J. Feng, J. Xu, M. Titirici, J. C. Rodriguez-Cabello, N. M. Pugno, J. Sun, W. Wang, T. Tuttle and A. Mata, Adv. Funct. Mater., 2022, 32, 2205802.
55. L. Liu, Q. Ma, J. Cao, Y. Gao, S. Han, Y. Liang, T. Zhang, Y. Song and Y. Sun, Cancer Nanotechnol., 2021, 12, 1–31.
56. L. Ruiz, W. Xia, Z. Meng and S. Keten, Carbon, 2015, 82, 103–115.
57. Z. Meng, R. A. Soler-Crespo, W. Xia, W. Gao, L. Ruiz, H. D. Espinosa and S. Keten, Carbon, 2017, 117, 476–487.
58. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. De Vries, J. Phys. Chem. B, 2007, 111, 7812–7824.
59. L. Monticelli, S. K. Kandasamy, X. Periole, R. G. Larson, D. P. Tieleman and S. J. Marrink, J. Chem. Theory Comput., 2008, 4, 819–834.
60. D. H. De Jong, G. Singh, W. F. D. Bennett, C. Arnarez, T. A. Wassenaar, L. V. Schäfer, X. Periole, D. P. Tieleman and S. J. Marrink, J. Chem. Theory Comput., 2013, 9, 687–697.
61. S. O. Yesylevskyy, L. V. Schäfer, D. Sengupta and S. J. Marrink, PLoS Comput. Biol., 2010, 6, e1000810.
62. G. van Meer, D. R. Voelker and G. W. Feigenson, Nat. Rev. Mol. Cell Biol., 2008, 9, 112–124.
63. T. A. Larson, P. P. Joshi and K. Sokolov, ACS Nano, 2012, 6, 9182–9190.
64. G. Frigerio, S. Motta, P. Siani, E. Donadoni and C. Di Valentin, J. Controlled Release, 2025, 379, 344–362.
65. A. van Teijlingen, M. C. Smith and T. Tuttle, Acc. Chem. Res., 2023, 56, 644–654.
66. J. A. Graham, J. W. Essex and S. Khalid, J. Chem. Inf. Model., 2017, 57, 650–656.
67. F. Grünewald, R. Alessandri, P. C. Kroon, L. Monticelli, P. C. T. Souza and S. J. Marrink, Nat. Commun., 2022, 13, 68.
68. F. Grunewald, G. Rossi, A. H. de Vries, S. J. Marrink and L. Monticelli, J. Phys. Chem. B, 2018, 122, 7436–7449.
69. E. Panizon, D. Bochicchio, L. Monticelli and G. Rossi, J. Phys. Chem. B, 2015, 119, 8209–8216.
70. A. V. Titov, P. Král and R. Pearson, ACS Nano, 2010, 4, 229–234.
71. https://github.com/Tsjerk/Insane.
72. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
73. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
74. R. W. Hockney, S. P. Goel and J. W. Eastwood, J. Comput. Phys., 1974, 14, 148–158.
75. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
76. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
77. E. L. Wu, X. Cheng, S. Jo, H. Rui, K. C. Song, E. M. Dávila-Contreras, Y. Qi, J. Lee, V. Monje-Galvan, R. M. Venable, J. B. Klauda and W. Im, J. Comput. Chem., 2014, 35, 1997–2004.
78. S. Jo, J. B. Lim, J. B. Klauda and W. Im, Biophys. J., 2009, 97, 50–58.
79. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
80. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021.
81. B. Roux, Comput. Phys. Commun., 1995, 91, 275–282.
82. J. S. Hub, B. L. De Groot and D. Van Der Spoel, J. Chem. Theory Comput., 2010, 6, 3713–3720.
83. T. D. Romo, N. Leioatts and A. Grossfield, J. Comput. Chem., 2014, 35, 2305–2318.
84. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.
85. R. Gowers, M. Linke, J. Barnoud, T. Reddy, M. Melo, S. Seyler, J. Domański, D. Dotson, S. Buchoux, I. Kenney and O. Beckstein, in Proceedings of the 15th Python in Science Conference, SciPy, 2016, pp. 98–105.
86. https://cgmartini.nl/docs/downloads/tools/other-tools.html.
87. M. Vögele and G. Hummer, J. Phys. Chem. B, 2016, 120, 8722–8732.
88. M. Linke, J. Köfinger and G. Hummer, J. Phys. Chem. Lett., 2018, 9, 2874–2878.
89. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
90. S. A. Van Der Heijden and M. T. O. Jonker, Environ. Sci. Technol., 2009, 43, 8854–8859.
91. D. C. Marcano, D. V. Kosynkin, J. M. Berlin, A. Sinitskii, Z. Sun, A. Slesarev, L. B. Alemany, W. Lu and J. M. Tour, ACS Nano, 2010, 4, 4806–4814.
92. J. Ma, R. Liu, X. Wang, Q. Liu, Y. Chen, R. P. Valle, Y. Y. Zuo, T. Xia and S. Liu, ACS Nano, 2015, 9, 10498–10515.
93. A. J. Paulista Neto and E. E. Fileti, J. Phys. Chem. B, 2018, 122, 2578–2586.
94. Y. Jin, Z. Xu, Y. Guo and X. Yang, J. Phys. Chem. C, 2018, 122, 4063–4072.
95. T. Szabo, P. Maroni and I. Szilagyi, Carbon, 2020, 160, 145–155.
96. M. Xu, J. Zhu, F. Wang, Y. Xiong, Y. Wu, Q. Wang, J. Weng, Z. Zhang, W. Chen and S. Liu, ACS Nano, 2016, 10, 3267–3281.
97. T.-R. Jan, H.-Y. Wu, K.-J. Lin, P.-Y. Wang, Y.-J. Lu, C.-C. M. Ma, H.-W. Yang and C.-W. Lin, Int. J. Nanomed., 2014, 4257.
98. L. Zhang, Z. Wang, Z. Lu, H. Shen, J. Huang, Q. Zhao, M. Liu, N. He and Z. Zhang, J. Mater. Chem. B, 2013, 1, 749–755.
99. E. Awoonor-Williams and C. N. Rowley, Biochim. Biophys. Acta, Biomembr., 2016, 1858, 1672–1687.
100. B. N. Olsen, A. A. Bielska, T. Lee, M. D. Daily, D. F. Covey, P. H. Schlesinger, N. A. Baker and D. S. Ory, Biophys. J., 2013, 105, 1838–1847.
101. F. Favela-Rosales, A. Galván-Hernández, J. Hernández-Cobos, N. Kobayashi, M. D. Carbajal-Tinoco, S. Nakabayashi and I. Ortega-Blake, Biophys. Chem., 2020, 257, 106275.
102. J.-P. Douliez, A. Leonard and E. J. Dufourc, Biophys. J., 1995, 68, 1727–1739.
103. C. J. Shih, S. Lin, R. Sharma, M. S. Strano and D. Blankschtein, Langmuir, 2012, 28, 235–241.
104. F. Ahmadpoor, G. Zou and H. Gao, Mech. Mater., 2022, 174, 104414.
105. Z. Li, X. Zhu, J. Li, J. Zhong, J. Zhang and J. Fan, Nanoscale, 2022, 14, 5384–5391.
106. Y. Li, X. Chen and N. Gu, J. Phys. Chem. B, 2008, 112, 16647–16653.

---