Biomimetic Phospholipid Membrane Organization on Graphene and Graphene Oxide Surfaces: A Molecular Dynamics Simulation Study

: Supported phospholipid membrane patches stabilized on graphene surfaces have shown potential in sensor device functionalization, including biosensors and biocatalysis. Lipid dip-pen nanolithography (L-DPN) is a method useful in generating supported membrane structures that maintain lipid functionality, such as exhibiting speci ﬁ c interactions with protein molecules. Here, we have integrated L-DPN, atomic force microscopy, and coarse-grained molecular dynamics simulation methods to characterize the molecular properties of supported lipid membranes (SLMs) on graphene and graphene oxide supports. We observed substantial di ﬀ erences in the topologies of the stabilized lipid structures depending on the nature of the surface (polar graphene oxide vs nonpolar graphene). Furthermore, the addition of water to SLM systems resulted in large-scale reorganization of the lipid structures, with measurable e ﬀ ects on lipid lateral mobility within the supported membranes. We also observed reduced lipid ordering within the supported structures relative to free-standing lipid bilayers, attributed to the strong hydrophobic interactions between the lipids and support. Together, our results provide insight into the molecular e ﬀ ects of graphene and graphene oxide surfaces on lipid bilayer membranes. This will be important in the design of these surfaces for applications such as biosensor devices.

D ip-pen nanolithography (DPN) 1 with phospholipids (L-DPN) 2 and the use of polymer pen lithography (PPL) 3 and other stamping techniques for the generation of lipid membranes on supports 4−6 have gained increasing interest in recent years. 7 Supported lipid bilayers have diverse applications in biomedical research, 8−11 sensor and device functionalization, 12−17 protein crystallization, 18 and in generating lipid sensor structures. 18−20 In L-DPN, the tip of an atomic force microscope (AFM) probe is covered with phospholipids and brought into close contact with a solid support, allowing the lipid ink to transfer to the substrate and self-assemble into stacks of membranes. 21 This method has the benefit of direct and precise spatial control during lipid deposition onto the support surface, with the ability to tailor the lipid mixture in the ink to a desired composition. Other methods of lipid deposition, including imprinting, have proved useful in elucidating the spreading behavior of lipid membranes on support surfaces of varying degrees of hydrophilicity and roughness and have indicated that different spreading mechanisms and velocities may occur depending on the nature of the surface. 22,23 An interesting difference when generating membrane patches by L-DPN compared with other methods (e.g., vesicle fusion, 24 micropipettes, 25,26 or spreading membranes in microfluidic systems 27,28 ) is that the fabrication of the lipid patches generally takes place in air (i.e., under ambient conditions). Therefore, the structures will usually be transferred into liquid only after the lithographic step is completed, as most applications take place in liquid phase. This suggests that the lipid structures may rearrange, and different structural models have been proposed for the molecular organization of the lipid membrane in air and water and on surfaces with varying hydrophilicity. 13,28−30 Computational methods such as molecular dynamics (MD) simulations have been used to study lipid membrane organization and dynamics on different support materials. 31−35 Coarse-grained (CG) simulations of lipid bilayers on hydrophilic model surfaces indicate that preformed lipid patches are more ordered on these surfaces compared to free-standing bilayers, affecting lateral lipid diffusion in the lipid leaflet that directly contacts the support and resulting in reduced lipid mobility. 35 Increasing surface roughness reduces lipid ordering and results in higher lipid mobility relative to smooth surfaces, confirming that lipid membrane organization is directly influenced by the underlying surface topology of the support. 35 Additionally, CG simulations modeling lipid selfassembly on supports of differing hydrophilicity suggest that lipid diffusion is more strongly affected, and reduced, on hydrophilic supports compared to on hydrophobic supports, due to stronger attractive forces. 36 The nature of the supporting surface is thus important. Recent studies have focused on parametrizing models of various support materials, for example, graphite, to more accurately represent interactions of organic molecules, such as long-chain alkanes, with the support. 37 These models were shown to reproduce phase transitions and molecular organization of organic molecules on the surface, in line with experimental data. 31,37 MD simulations are therefore useful in studying the molecular properties of lipid−surface interactions, allowing graphene surface in vacuum. The lipid choline, phosphate, glycerol, and carbon groups are shown as licorice representations colored in green, red, yellow, and gray, respectively. The graphene surface is shown as black van der Waals spheres. The Hyperballs program 41 was used for image generation. Partial density profiles were calculated for the last 25% of the simulation time for the inverted bilayer (C) and monolayer (D) simulations. (E) and (F) The time evolution of the average COM distance between the lipid phosphate head-groups and the graphene surface (dashed line) for the inverted bilayer and monolayer simulations in (A) and (B), respectively. See also Movie S1 and Movie S2.

ACS Nano
Article microscopic details to be characterized. In this study, we examine the organization of supported lipid bilayers as generated by L-DPN using CG molecular dynamics (CG-MD) simulations and correlate the results with experimental evidence. Our simulation results are in close agreement with AFM measurements of lipid layers on graphene and graphene oxide supports. The AFM experiments suggested that inverted bilayer configurations are formed on pristine graphene, while altered lipid structures can be formed on a graphene oxide surface. Our simulations reproduce this behavior, capturing the dynamic process underlying the altered lipid configurations. Furthermore, the simulations revealed effects on lipid ordering and diffusion within the bilayer structures relative to freestanding bilayers in water, depending on the support surface. Together, the combined experimental and simulation results provide an integrated study of the dynamical properties of lipid membranes on graphene-based supports.

RESULTS AND DISCUSSION
Previous studies have shown that the formation of stable lipid structures, such as multilayered lipid membranes, on graphene, silicon dioxide, and other substrates can be generated by L-DPN and are influenced by surface morphology and experimental conditions. 13,38−40 Here, we used CG-MD simulations to investigate the molecular details of lipid interactions with pristine graphene and with graphene oxide surfaces. The simulated systems were constructed to match the experimental conditions used during L-DPN and AFM studies as well as currently possible. These experiments were either performed for systems in air or systems immersed in liquid. 13 Experiments in air were typically performed at 20−30% relative humidity (RH). Even for the largest systems simulated (30 × 20 × 20 nm 3 ) this corresponds to only 2 to 3 water molecules, i.e. <1 CG water particle, within the simulation box. The experiments in air were, therefore, reasonably well approximated by simulations in vacuum.
Lipid Interactions with Graphene. L-DPN generated lipid structures on pristine graphene surfaces in air are thought to form inverted bilayer structures. 13,39 We investigated the stability of these inverted bilayer structures on a graphene surface in vacuum and in water using CG-MD. Two different graphene surface areas were modeled: a small (10 × 10 nm 2 ) and large (30 × 30 nm 2 ) surface. Preformed inverted bilayers consisting of 512 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) molecules (256 lipids per layer) were individually placed above the small and larger graphene surfaces and simulated for 0.5−1 μs. Simulations of the systems in vacuum (4 replicate simulations each of 4 μs duration) resulted in rapid adsorption of the lipid layers on the surface, mediated by strong hydrophobic interactions between the graphene surface and the lipid tails ( Figure 1A). The lipid structures maintained their inverted configurations ( Figure 1C), as shown by monitoring the center of mass (COM) distance between the lipid phosphate headgroup and the graphene surface ( Figure 1E). Visual inspection also indicated some degree of enhanced ordering of the lipid tails in the graphene-proximal layer. Polar interactions between the head-groups resulted in clustering of the lipid head-groups, stabilizing the inverted bilayer topology. Furthermore, CG-MD simulations of lipid monolayers with the small graphene surface (again, 4 replicate simulations each of 4 μs duration) resulted in spontaneous lipid reorganization to form an inverted bilayer structure in vacuum ( Figure 1B,D,F).  13 The bottom profile section reveals a bilayer membrane on top of a lipid monolayer (wetting layer). The wetting layer is similar to the layers observed on silicon dioxide, though thinner, probably due to reduced layer density. 28

Article
These results corroborate AFM measurements of L-DPN deposited lipid layers on pristine graphene surfaces, 13,39 revealing that phospholipids form a flat inverted bilayer of uniform thickness (Figure 2A) on the hydrophobic support in vacuum (simulation) or air (experiment). The measured thickness of the inverted bilayer structure from the simulations compares well to the height measurements from the AFM experiments, corresponding to ∼4 nm. This can be seen from the density peaks of lipid tail groups as shown in the partial density profiles in Figure 1B,D. This matches the average height of ∼4−4.5 nm of the inverted bilayer as measured by AFM ( Figure 2A).
However, the inverted bilayer topology was less stable (in terms of maintenance of the inverted bilayer structure after interactions with the support surface) in CG-MD simulations with the large graphene surface (30 × 30 nm 2 ) compared to the small graphene surface (10 × 10 nm 2 ). In simulations with the large graphene surface, the strong hydrophobic interactions between the lipid tails and graphene beads dominated lipid organization, resulting in disassembly of the layer and spreading of lipids across the available surface area, aligning with the surface in a lateral fashion ( Figure S2A). The same behavior was also observed in simulations of preformed regular bilayers with the large graphene surface ( Figure S2C). This contrasts with our initial observations of stable inverted bilayer configurations during simulations with the small graphene surface area, which did not disassemble. The difference between these simulations is the overall surface area of the graphene model. The preformed inverted lipid structures are roughly 10.5 × 10.5 nm 2 in x,y dimensions (prior to simulation) and thus occupy a slightly larger surface area than the small graphene surface (10 × 10 nm 2 ). Therefore, the larger density of lipids on this small graphene surface may affect the stability of the observed lipid configuration compared to the large graphene surface, resulting in maintenance of the inverted structure. Experimentally, inverted bilayer patches generated by L-DPN are stable on more extended graphene surfaces (20 × 20 μm 2 ), for which lipid spreading is thought to be a self-limiting process. 39 Consequently, the difference in the lipid density on the large graphene surface areas within the AFM experiment compared the simulations may be a deciding factor in the observed disassembly of the inverted bilayer structure (512 lipids) on the large graphene surface in the simulation. Interestingly, simulation of larger inverted bilayers (2110 lipids) exhibited higher stability and maintained the inverted configuration to a larger extent compared to the smaller bilayer systems (512 lipids) ( Figure S2E). However, lipid tail entanglement and lateral interactions with the surface, as well as headgroup clustering, resulted in deviation from the more "ideal" inverted bilayer configuration observed in simulations of the smaller bilayers on the small graphene surface in vacuum (10 × 10 nm 2 ; Figure 1). Importantly, the increased stability of larger lipid layers suggests that lipid density on the graphene surface is another determinant of the topology of a supported membrane.
Lipid Interactions with Graphene Oxide. Surface polarity is known to affect the molecular properties of supported lipid bilayers, with effects on topology, 29,30

ACS Nano
Article diffusion, 42 and lipid order. 36,43 To investigate these effects, L-DPN was also performed on graphene oxide substrates in air. AFM height measurements of the deposited lipids on the hydrophilic substrate indicated that the lipids organized into a "1.5 bilayer" configuration ( Figure 2B), as has also been suggested for lipid structures on silicon dioxide surfaces. 13 Two distinct lipid layers could be distinguished: a "wetting" layer composed of phospholipids with their head-groups oriented toward the hydrophilic support surface, and a second inverted bilayer that formed on top of this wetting layer. The same lipid organization is seen in CG-MD simulations (of duration 5 μs) of (non-inverted) bilayers with graphene oxide in vacuum. The lipid bilayer interacted with the surface, undergoing substantial reorganization within ∼0.2 μs to form a 1.5 bilayer on top of the oxidized surface ( Figure 3A). Simulation (also of duration 5 μs) of a preformed inverted bilayer configuration positioned above the graphene oxide surface converged to a similar 1.5 bilayer structure ( Figure 3B), indicating that the outcome was robust to the initial bilayer (inverted vs non-inverted) model used. The significant rearrangements facilitating formation of . Partial density profiles were calculated for the last 25% of the simulation time shown for the graphene (E, G) and graphene oxide (F, H) systems. The same color scheme is used as in Figure 1. The graphene oxide surface is shown in black (carbons) and red (oxygens). The graphene surface is shown in black.

ACS Nano
Article the 1.5 bilayer structure are thought to be driven by polar headgroup interactions with the hydrophilic graphene oxide surface. Partial density profiles calculated for the last 25% of CG-MD simulation time show very similar peaks in lipid headgroup and tail densities for both the inverted and noninverted bilayer systems ( Figure 3C,D), confirming that the differing starting structures converged to similar end structures. These results support the interpretation of a 1.5 bilayer configuration by AFM height measurements for supported lipid membranes on graphene oxide, suggesting that this arrangement is the preferred molecular state of lipid layers on polar oxide surfaces in air.
In contrast, the 1.5 lipid bilayer configuration (which exposes the hydrophobic tails of the lipid molecules at the "uppermost" layer) was not observed in CG-MD simulations of the lipidgraphene oxide systems in water. Instead, preformed lipid bilayers reorganized to form bicelle-like configurations on both small and larger graphene oxide surfaces in water ( Figure 4E− H). The rearrangement of inverted bilayers was driven by lipid headgroup interactions with the underlying support, either directly or through bridging water particles, as well as by interactions with the surrounding solvent, resulting in a stable bicelle structure which does not subsequently move as a whole on the graphene oxide surface. Conversely, the addition of water to lipid-covered graphene systems in the CG-MD simulations resulted in destabilization of the previously formed inverted bilayers (formed in vacuum) and triggered disassembly of the structures ( Figure 4A−D). Subsequently, the lipids spread over the available surface area with head-groups oriented toward the surrounding solvent molecules. Importantly, similar behavior is observed for lipid-covered graphene transferred to an aqueous environment by AFM measurements.
Two observations can be made from the combined simulation and experimental data. First, the addition of aqueous solvent can result in different lipid configurations on the same surface, for example, bicelle-like formations on graphene oxide in an aqueous environment compared to a 1.5 bilayer configurations in air (or vacuum). 13,21,28 Second, the polarity difference resulting from the oxygen-containing functional groups in graphene oxide compared to the hydrophobic surface presented by pristine graphene resulted in stabilization of different lipid configurations, for example, 1.5 layers vs inverted bilayers. These observations imply that the configuration of the supported membrane can be tuned as a function of system solvation as well as underlying surface polarity, resulting in relatively different end structures.
Direct observation of lipid reorganization on graphene and silicon dioxide in liquid by AFM has been reported by a previous study. 13 These observations necessitated the use of bovine serum albumin (BSA), which binds to the surface, to

ACS Nano
Article block excessive lipid spreading and stabilize the membrane patches in aqueous solution. In this study, direct AFM measurements of the lipid layer structures on graphene oxide in liquid have been hindered by the blocking layers of BSA needed to stabilize the small membrane patches during scanning. Given that specific binding interactions in L-DPN applications are usually conveyed by functionalization at the lipids headgroup, it was inferred that the lipid patches undergo reorganization on graphene oxide, exposing the (otherwise buried functional groups) to the liquid phase, similar to the lipid configurations identified from the CG simulations of the graphene oxide systems in water. 13,44 Lipid Layer Properties on Graphene and Graphene Oxide: Lipid Order Parameters. So far, we have provided an overview of the topological differences generated in lipid organization on surfaces of varying polarity and system solvation. How does this affect the dynamic properties of the lipids within the observed structures? Calculation of lipid order parameters (S) as a function of simulation time provides a metric for lipid rearrangements on a support surface (see Methods). Briefly, the lipid order parameter is calculated by estimating the angle between the bonds connecting the lipid particles and the z-axis of the simulation box (i.e., the normal to the bilayer). Thus, alignment of the lipid particles with the zaxis yields an order parameter of S = 1, whereas a value of S = 0 corresponds to randomly orientated lipids. Order parameters were calculated for the phosphatidylcholine lipids of selected systems and compared to those for a free-standing bilayer (i.e., no supporting surface) simulated in water (Table S1).
Systems that underwent substantial lipid reorganization showed major changes in the lipid order parameters over time, particularly for the simulations starting from an inverted bilayer on the graphene oxide surface in vacuum ( Figure 5). In this system, the average order parameter values change significantly during the first 2 μs of simulation time, particularly for the phosphate-glycerol bond and consecutive bonds for both sn-1 and sn-2 lipid tails up to the third CG particle. The changes represent altered configurations of the lipids as the molecules reorganized to form the 1.5 bilayer structure. Specifically, the average order parameters for particular lipid bonds decrease from around 0.3 to ∼0.1 during the first microsecond of the simulation time, indicating a random arrangement for many lipids during formation of an inverted bilayer on top of the monolayer (wetting layer). Toward the end of the rearrangement (∼2 μs), the average order parameters, for example, the phosphate-glycerol bond, return to 0.3, suggesting lipid alignment with the z-axis following the transition ( Figure 5). Other order parameter values, such as the C2−C3 bond, however, did not completely recover their initial value, which may be attributed to a degree of bilayer distortion brought about by micelle-like lipids within the bilayer formed on top of the wetting layer ( Figure 3A). This value (∼0.1) is reflected in those for subsequent bonds (e.g., C3−C4) in both sn-1 and sn-2 chains, characterizing the random arrangement of lipids tails within the 1.5 bilayer structure on graphene oxide.
In general, the average lipid tail order parameters for all systems suggest that the tails were less ordered on both graphene and graphene oxide surfaces compared to a freestanding bilayer in water ( Figure 6; Table S1). For example, an average order parameter of 0.18 is obtained for lipid tails in an inverted bilayer-graphene simulation in vacuum (small surface area), compared to 0.36 for lipids within the free-standing bilayer. However, inspection of the time evolution of the order parameters from the inverted bilayer-graphene simulation in vacuum reveals a slight gain in order over the last 0.5 μs of simulation time, particularly for the sn-2 chain of the lipids (e.g., C1B−C5B) ( Figure S3A). This transition indicates the increased alignment of the lipid tails, as the inverted bilayer configuration becomes more stable on the pristine graphene surface. Furthermore, simulations of different starting lipid structures with graphene oxide resulted in similar order parameter values, particularly toward the end of the simulations, further indicating convergence to the 1.5 bilayer organization on this surface in vacuum ( Figure S3C and D).
The overall values compare well with lipid tail order parameters calculated from CG simulations of related systems. 36 The overall reduced lipid order on the graphene (oxide) supports could be a result of lateral interactions between the lipid tails and the surface (partial alignment of the lipid tails with the surface), as is observed in the pristine graphene systems (due to strong hydrophobic forces), and/or headgroup interactions with the graphene oxide surface, disrupting lipid order in the layer.
Lipid Layer Properties on Graphene and Graphene Oxide: Lipid Diffusion. Characterization of the lateral diffusion of lipids within supported membrane systems is fundamental to understanding the dynamic properties of phospholipids within these membranes. 45−48 Experimental studies report both linear and anomalous diffusion regimes of lipids depending on the systems and increasingly highlight the importance of the time and length scale over which the diffusion data are collected when making a distinction between anomalous and linear diffusion. 42,46,47,49 To characterize lateral diffusion of lipids within our simulated supported bilayer systems, both linear and anomalous diffusion models were considered. Mean square displacements (MSDs) were collected by tracking the displacement of lipids (see Methods). The resultant MSDs were then fitted with either a linear 50 or an anomalous 51 diffusion equation.
This analysis revealed differences in the lateral displacement of lipids on surfaces relative to those free-standing bilayers. The diffusion of lipids within supported membranes, on either graphene or graphene oxide, was slower in the presence of the support compared to free-standing bilayers (Figures 7A and  S5). 52 The presence of water appeared to result in increased lipid diffusion on the graphene oxide surface (Figures 7B and

ACS Nano
Article S6), suggesting the water layers could act as a lubricant for lipid diffusion within the layer. Indeed, bridging water particles between the lipids and the support surface were often observed within the simulations with this surface. Similar effects have also been observed in experimental systems of supported lipid membranes on hydrophilic surfaces in an aqueous environment. 46 Furthermore, fitting MSD data using the anomalous diffusion equation produced α values of 0.7−0.9 for all the supported membrane systems, suggesting that the lipids within the supported membranes may exhibit different diffusion regimes as a consequence of their interactions with the surface ( Figure  8). In particular, lipids exhibited subdiffusion (i.e. α < 1) for all supported membrane simulations.
Importantly, the displacement of lipids within the freestanding bilayer was broadly consistent with a linear diffusion model, producing a value of D lin = 3.8 × 10 −7 cm 2 /s (±0.5). This value should be scaled by a factor of 4 to account for the reduced degrees of freedom within CG simulations 53 in order to compare with experimental estimates, resulting in D lin = 0.98 × 10 −7 cm 2 /s, which compares well to atomistic simulations of DOPC bilayers in water (1.5 × 10 −7 cm 2 /s). 54,55 However, it does remain challenging to compare diffusion data for membranes on surfaces with experimental values, or even other simulation studies, given the dependence on both the time scale over which diffusion is measured and the length scale of the system. 42,45,47,49 Nonetheless, we can reasonably conclude from the MSD data and fits across all of the supported membrane systems (see Figures S5−6) that in general D is lower for lipids in the supported bilayers relative to a free-standing bilayer, with the largest reduction (more than 2fold, with α = 0.7) for a membrane on pristine graphene or on graphene oxide in vacuum.

CONCLUSIONS
The dynamic organization of lipid membranes on graphene supports has been investigated using both experimental and computational techniques. As has been observed in a number of other studies of supported lipid membranes, the dynamic properties of the lipid configurations presented here were affected by the polar/apolar nature of the surface as well as system hydration. 32,[34][35][36]56 AFM measurements suggested the formation of inverted bilayer topology on pristine graphene, while a 1.5 bilayer structure was observed on graphene oxide in air ( Figure 2). 13 CG simulations of these systems supported initial observations of the differing lipid topologies on the different surfaces. Encouragingly, different preformed starting lipid structures spontaneously rearranged to converge to the same end structure, as observed by the AFM studies. For example, regular bilayers formed 1.5 bilayers on graphene oxide, while single lipid monolayers formed inverted bilayer structures on pristine graphene, converging to the observed end structures despite different starting configurations (Figures 1 and 3). Furthermore, the CG simulations suggested that stability of the lipid organization was related to lipid density on the surface. While small inverted bilayers were observed to disassemble on large graphene surfaces due to strong hydrophobic interactions

ACS Nano
Article with the surface, larger bilayers were more stable ( Figure S2). Importantly, initial AFM characterization of phospholipid patches on pristine graphene in air indicated that lipids are more mobile on this surface than on hydrophilic surfaces, such as silicon dioxide, also attributing to the strong hydrophobic interactions between lipid hydrocarbon tails and graphene. 13,39 Spreading of a very small membrane patch on a large graphene surface is, thus, perhaps not unexpected.
System hydration was an important factor affecting lipid stability and organization on support surfaces. Specifically, lipid organization on graphene oxide is affected by the addition of water, resulting in formation of stable bicelle-like structures on graphene oxide surfaces, contrasting the preferred 1.5 bilayer lipid organization observed on the same surface in vacuum (simulation) or in air (AFM). Given the apparent high stability of bicelle structures on graphene oxide in water, patterning of hydrophilic surfaces in a discrete manner should be possible in aqueous support systems. A recent study investigating lipid localization across lipid monolayer−bilayer junctions exploited precisely this topological difference in lipid layering, which results from the different energies of interaction of lipids with surfaces of varying hydrophobicity. 57 Thus, lipid monolayer− bilayer junctions were established by patterning glass surfaces with hydrophobic molecules, resulting in formation of lipid monolayers on the hydrophobic sections and lipid bilayers on the surrounding hydrophilic (glass) sections. 57 Lipid bilayer formation on graphene oxide in solution has also been observed by Okamoto et al., who used AFM measurements to propose that both single and double lipid membranes are thermally stable on graphene oxide supports after vesicle fusion with the surface. 30 Furthermore, a recent study investigating lipid assembly on oxidized graphene surfaces confirmed that planar lipid bilayers were stable on this support. 58 The same study also observed that, in contrast with lipid bilayer formation on oxidized graphene, lipid monolayers were formed on a pristine graphene surface in solution, in which lipid tails interacted with the underlying graphene support. This is similar to what is suggested by our CG-MD simulations of preformed lipid bilayers on larger graphene surfaces, which disassembled into monolayers in the presence of water, again suggesting that the hydrophobic interactions between the lipid carbon tails and graphene surface dominate lipid organization on this support. A recent functional study of L-DPN generated membranes on pristine graphene observed the same lipid behavior upon the addition of buffer solution, in which fluorescence was used to verify the solvent-exposed orientation of the lipid head-groups with respect to surface. 39 The study also indicated that lipid spreading on pristine graphene is self-limiting and is more favorable than lipid interactions with silicon dioxide surfaces, attributing to the contrast in hydrophilicity between the two supports. 39 A similar trend in lipid spreading behavior has also been observed on other support surfaces, including glass (hydrophilic) and n-octadecyltrichlorosilane (hydrophobic). 59 Using ellipsometry to determine lipid film heights, it was suggested that lipids preferentially spread as monolayers on hydrophobic surfaces, while lipid bilayers form on more hydrophilic surfaces. 59, 60 We also observed that lipids may undergo dynamic changes in their configurations on either support surface, as is evidenced by the time evolution of lipid tail order parameters. This is particularly true for simulations of regular (i.e., non-inverted) lipid bilayers on the small graphene oxide surface, in which complete restructuring of the lipid layers occurred to form the apparently more favorable 1.5 bilayer configuration ( Figure 5). Furthermore, characterization of lipid mobility within the supported membrane structures suggested that lipid diffusion was anomalous, exhibiting subdiffusion for most systems. This is expected given that interactions between the lipid and the support alter MSD values over the longer time sampling windows. Tero et al. also identified subdiffusive behavior of DOPC molecules within supported lipid membranes on hydrophilic titanium oxide and related this to the atomic topology of the support surface. 42 Furthermore, other experimental studies of lipid diffusion in supported membranes have reported both linear and anomalous diffusion regimes and related this to lipid interactions with the support surface, surface topology of the support, membrane composition, and topology (e.g. planar bilayer vs vesicles) as well as system hydration. 42,46,47,49 The distinction between diffusion regimes is, however, also related to the temporal and spatial scales at which diffusion is measured. 45 It would, therefore, be of interest in the future to extend our simulations to substantially larger length scales and longer time scales, enabling more direct comparison with experimental data for similar systems. 61 The use of CG resolution simulations of the lipid membranegraphene systems has allowed us to characterize supported membrane properties for longer (i.e. several microsecond) simulation times than would be readily available, for example, all atom simulations. A number of previous studies have indicated the utility of CG simulation methods for studying lipid behavior on surfaces. 31,34−36 It is important to note that the mapping scheme used to model the CG graphene and graphene oxide surfaces (2 atoms:1 CG particle) differs from the 4:1 mapping scheme used for the lipid molecules. However, the Martini force field was parametrized to remain internally consistent and compatible between parametrization efforts (e.g., comparing Martini v2.1 and Martini v2.2). 37,53 Consequently, the parameters used to model the interactions between the lipid molecules and graphene/graphene oxide surfaces are largely independent of the mapping schemes employed to generate the graphene and lipid models. This is because the interactions between the different components are modeled using an interaction matrix, with different values assigned to interactions of each bead type. 53,54 The original parametrization of the graphite surface on which the graphene (oxide) model in this study is based involved reproducing thermodynamic values and adsorption behavior of long-chain alkanes on the graphite support. 37 The nonbonding interaction parameters between the different beads within the Martini models were adjusted to reproduce this behavior, so that the different mapping schemes did not affect the overall behavior of the systems.
Full atomistic simulations would be of interest to capture in more detail the interactions between lipid molecules and the support surface in more detail. Importantly, the electronic structure of graphene and the polarizability of interacting molecules, which may play key roles in biomolecular adsorption and interactions, 62−64 are not captured by the CG model. With the development of polarizable force fields 62,65,66 and further characterization of the unique properties of graphene supports, 63 future simulations encompassing these aspects, while starting from the CG models described in the current study, would provide enhanced insight into the interaction forces underlying different supported membrane configurations. It would be of interest to explore the free energy landscape underlying the different lipid configurations observed on these support surfaces, as this would be likely to provide further

ACS Nano
Article insights into the origin of stability of the lipid structures observed for each of the surfaces. Characterization of the energetic differences between the different lipid structures formed on pristine graphene and the graphene oxide depending on system solvation might also be expected to provide insights into the link between surface wettability and the eventual lipid structure formed. However, accurate estimation of free energy landscapes for lipid/graphene interactions would be challenging in terms of achieving adequate sampling of the system to reach convergence. In a previous study, we have shown that both the choice of a suitable collective variable and adequate sampling are key parameters in obtaining a converged free energy landscape even for relatively simple systems, such as single protein−lipid interactions. 67 In conclusion, our combined simulation and experimental results highlight the effects of both surface polarity and the solvent environment on phospholipid membrane organization when interacting with a support surface. We demonstrate that lipids can form multilayered structures such as 1.5 bilayers on hydrophilic surfaces, such as graphene oxide, and can spontaneously rearrange to form a preferred topology even when starting from different structures (e.g., regular bilayers). Hydrophobic surfaces such as pristine graphene interact highly favorably with lipids, affecting lipid bilayer membrane stability on this surface. Characterizing these interactions will provide important insight into the applications of graphene/graphene oxide within biotechnology, including sensor devices and drug delivery systems. 68,69 METHODS Coarse-Grained Graphene (Oxide) Models. CG simulations were performed using the Martini force field. 53,70 The Martini model for graphite 37 was used to construct the graphene and graphene oxide model surfaces (single layer; Figure S1). Briefly, the model involved a 2:1 mapping scheme of atomistic carbon atoms into hexagonally packed graphene beads (SG4) making up the sheet. The published parametrization of this model was based on reproducing the adsorption and topological behavior of long-chain alkane molecules on graphite, suggesting the model is suitable to study the behavior of lipids on graphene surfaces. 37 The CG graphene model was used to build the graphene oxide model. To do this, the carbon beads (SG4) comprising the pristine graphene surface were randomly substituted by oxygen beads (SP1), which represent either a C−O−C (epoxy) group or a C−O−H group, 71 eventually reproducing the oxygen content of the graphene oxide substrates used in the experiments (2.82% C−O− C: 49.3% C−O−H based on X-ray photoelectron spectroscopy data) ( Figure S1).
Simulation Details. The GROMACS 4.6 (www.gromacs.org) simulation suite 72 and the Martini 2.2 force field 53,54,71 were used to perform all CG simulations. Lipid bilayers consisting of DOPC molecules were constructed by self-assembly simulations; this lipid composition reflects the main lipids employed in the experiments and related studies. 13,39 Lipid molecules were randomly inserted into cubic box of dimensions 10 × 10 × 5 nm 3 (small bilayers) or 20 × 20 × 5 nm 3 (large bilayer); the systems were then energy minimized using the steepest descent algorithm for 10,000 steps. The z-dimension of the box was then extended to 15 or 30 nm. The system was subsequently solvated and simulated for 100 ns to allow bilayer self-assembly to occur. The inverted bilayers that mimic lipid structures seen in previous AFM studies 13,39 were constructed by placing two lipid monolayers on top of each other, with the lipid head-groups pointing toward each other, using the editconf tool from GROMACS. These lipid configurations were then placed above either graphene or graphene oxide surfaces.
The polarizable Martini water model 73 was employed in simulations containing water. A round of equilibration simulations were performed for the solvated systems. This involved performing 10,000 steps of equilibration using incremental simulation time steps of: 1, 2, 5, 10, and 20 fs. The equilibrated systems were used to initiate production simulations performed at constant temperature and pressure (NPT ensemble). Temperature was regulated using the Berendsen thermostat 74 and a coupling constant of 0.3 ps at 298 K. Pressure was regulated using the Berendsen barostat 74 with a coupling constant of 3.0 ps, applying anisotropic pressure coupling using a compressibility of 0.5 × 10 −5 bar −1 in x and y and 3.0 × 10 −5 bar −1 in the z-direction. Simulations of systems in vacuum were performed at constant temperature and volume (NVT ensemble). The lipids, water, and graphene were coupled to separate external baths. Non-bonding interactions were modeled using shift functions; both LJ and Coulombic interactions were evaluated within a 1.2 nm cutoff and shifted within a 0.9 nm cutoff distance. These parameters reflect those applied in parametrization studies of the CG graphene model. 37 Lipid Order Parameter Analysis. Second rank order parameters were calculated for all of the lipids in the select CG systems ( Figures  S3 and S4). These were calculated according to the S n cc second rank order equation: where θ is the angle between the bond connecting lipid bead particles (B n − B n−1 ) connecting beads within the lipid and the z-axis of the simulation box, normal to the bilayer. The angle brackets denote the mean angle calculated for all of the lipids in the system (ensemble average). The second rank term refers to the second-order Legendre polynomial used to describe the order parameter. 75 This has been used in many experimental and computational studies of similar bilayer systems. 55,76−78 Diffusion Analysis. The diffusion of lipids within the simulated systems was analyzed using documented open-source code (http://dx. doi.org/10.5281/zenodo.11827). This code employs an algorithm which calculates the MSD (MSD) values of individual lipid centroids over a range of time sampling windows, including: 1, 2, 5, 10, 20, 50, 70, 100, and 200 ns. The diffusion coefficients are then calculated by fitting the MSD vs time data using a linear diffusion equation 50 or an anomalous diffusion equation. 51 The linear approximation uses a leastsquares first-degree fit of the data, whereas the anomalous approximation uses a nonlinear least-squares fitting to a two-parameter equation of the form: where D α is the fractional diffusion coefficient, measured in units of length 2 time −α . For the linear diffusion fit, standard deviations of the diffusion coefficients were estimated as the difference of the slopes from the first and second halves of the MSD vs time data. This is the approach used by the GROMACS tools function g_msd. For the anomalous diffusion fit, the scaling exponent (α) was estimated as the slope of the log MSD vs log time data. The standard deviations of both parameters (D α and α) were calculated from the square root of the diagonal of the covariance matrix from the anomalous fit. This code has been used by previous simulation studies of lipid diffusion in virus particles. 79,80 Diffusion constants are reported without correction for the reduced degrees of freedom resulting from applying the Martini CG model. A simple scaling factor of 4 has previously been applied to compare diffusion of lipid and water particles in CG systems to experimental measurements. 54 Lipid Dip-Pen Nanolithography (L-DPN). L-DPN was performed on graphene and graphene oxide samples. Lithography took place in a DPN5000 system (Nanoink, U.S.A.) with a single cantilever probe (A-Type, ACST, U.S.A.). The cantilever was coated with 1,2dioleoyl-sn-glycero-3-phosphocholine (DOPC) using a microfluidic inkwell (ACST, U.S.A.) at high humidity (∼70%). The tip was then conditioned by writing some sacrificial patterns to strip off excessive ink. Membranes used for AFM imaging were written at 25% humidity.
Atomic Force Microscopy (AFM). AFM imaging was performed on a Dimension Icon setup (Bruker) in tapping mode. NSC15

ACS Nano
Article cantilevers (MikroMasch) with a nominal force constant of 46 N/m and a resonance frequency of 325 kHz were used.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsnano.6b07352.
Additional details of the simulated graphene and graphene oxide models, as well as lipid diffusion and order parameter analyses (PDF) Movie S1 corresponds to the simulation in Figure 1A (AVI) Movie S2 corresponds to the simulation in Figure 1B (AVI) Movie S3 corresponds to the simulation in Figure 3A (AVI) Movie S4 corresponds to the simulation in Figure 3B (AVI)