Size-dependent Interaction of Nanoparticles with Nonionic Bilayers† Sang Young Noh,∗a and Rebecca Notmanb Understanding the mechanism of transit of a nanoparticle (NP) through a biomimetic bilayer has been at the forfront of research for the design of efficient drug-delivery mechanisms, nanotechnology and biomedicine. Establishing a consistent picture of how the transit mechanism depends on the physiochemical property of a NP is ritical to understanding what approach may be the most effective for nanomedicine design. In this study, using molecular simulation techniques, we have analyzed the key properties of a NP that may affect the mechanism of transit - the effect of size and hydrophobicity. By using a continuum model of a NP based on the Hamaker potential, we have created NP of tunable hydrophobic properties. The effect of hydrophilic, hydrophobic, and mixed properties of the NP is analyzed against a biomimetic bilayer - we show that this model can illustrate three distinct properties - where the hydrophilic type shows rupture of the bilayer, the hydrophobic type showing a entrapment of the NP around the hydrophobic tailgroups of the bilayer, and the mixed type showing a distinct, direct translocation type mechanism. Increasing the NP size shows different effects for each type of NP, and hence, may provide insight into the design of NPs with these types of mechanisms involved. 1 Introduction The self-assembly of nanoparticles (NPs) and amphiphilic macromolecules offers a powerful route to generate functional soft materials with controllable structure and properties. At the most basic level, polymeric vesicles are liquid-containing sacs and have a compartmentalised enclosed volume shielded from the outside liquid environment by a bilayer commonly composed of amphiphilic molecules. The tendency of surfactants to aggregate into large vesicular structures allows polymeric amphiphiles to be used as a model to emulate properties of more complex biomembranes. As a result, numerous studies with representations of complex bilayers have been attempted - biomembranes contain numerous heterogeneous components and factors such as the presence of membrane proteins, varying rigidity of the lipids and varying charge density near the headgroups, to name a few 1 . These factors control the curvature, diffusive properties a Department of Chemistry, University of Warwick, Coventry, United Kingdom. E-mail: s.y.noh@warwick.ac.uk b Department of Chemistry, University of Warwick, Coventry, United Kingdom. E-mail: r.notman@warwick.ac.uk † Electronic Supplementary Information (ESI) available: [details of any supplementary information available should be included here]. See DOI: 00.0000/00000000. ‡ Additional footnotes to the title and authors can be included e.g. ‘Present address:’ or ‘These authors contributed equally to this work’ as above using the symbols: ‡, §, and ¶. Please place the appropriate symbol next to the author’s name and include a \footnotetext entry in the the correct place in the list. and mechanisms of translocation through a bilayer. Given the need to allow selective permeation of molecular species, the additional mechanisms biolipid systems have evolved involve complex mechanisms which allow passive and active diffusion, controlled by a concentration gradient or an energy-input driven process respectively. Despite the complexity of the systems involved, recent successes in particle delivery, such as the development of inhalable insulin 2 , amongst others, has highlighted the possibilities of NPs as novel agents in clinical applications. Difficulties remain, however, in determining which component-bycomponent interactions are responsible for phenomena ranging from the direct translocation of a NP to the rupture of the bilayer. Questions remain with regards to the mechanism of how foreign particles can interact near the bilayer interior - what properties dictate the mechanism of transit? Passive translocation processes for example, do not require an external energy input, which would imply that the surface functionalization of the NP primarily dictates the mechanism of membrane translocation 3 . The purpose of this work is to provide an overview of the mechanism of transit of a NP with respect to it’s surface composition. To address this challenge, we have investigated the mechanisms by which NPs of varying size and hydrophobicity interact with and cross biomimetic polymer membranes. To this end, we have simulated a polyoxyethylene-glycol (PEG-PEO) surfactant bilayer; this surfactant provides a large variety of phase transitions depending on the temperature, concentration and on 1–7 | 1 structural properties, and computationally intensive coulombic potentials are not taken into account, which alleviates the MD computation. We have developed nonspecific model NP-bilayer systems which comprised a coarse-grained poly(ethylene)6-blockpoly(ethylene oxide)2 (C12 E2 ) bilayer in water and a generic NP of 1.0 or 2.0 nm radius (R) with hydrophilic, hydrophobic or mixed character (Hence, nine different systems to analyze in total). Similar continuum style models, targeting the efficiency through the use of generic parameters have been studied before using DMD/LME simulations of hydrophilic models 4 . However, using similar continuum models to tune the hydrophobicity has not been attempted before. We show that by adjusting this continuum model of a NP, we can reproduce many of the properties of hydrophilic/hydrophobic NPs and also show that NPs of mixed hydrophobicity can undergo a direct translocation. From our simulations, we have identified three distinct, size-dependent mechanisms - where the hydrophilic NP forms a pore like deformation of varying sizes, depending on the radius of the NP. For the mixed character type, we show that this NP models a direct penetration through the bilayer, with minimal disruption in the bilayer, coupled with a flattened free energy profile. With the hydrophobic NP type, we show that the hydrophobic groups within the bilayer wraps around the NP to maximise the hydrophobic interactions its hydrophobic components. 2 OA EO EO EO CM CM CM CT2 EO OA CT2 CM CM CM Fig. 1 Coarse-grained mapping of the surfactant C12 E2 used in this study. The EO and OA groups constitute the hydrophilc headgroups, while the CM and CT2 beads represent the hydrophobic tailgroups Simulation Methodology 2.1 Coarse-grained model We have useed the Shinoda-DeVane-Klein (SDK) CG model 5,6 which was recently expanded to the SPICA model 7 . The model has been parameterized against thermodynamic properties (densities, interfacial tensions, transfer free energies) and has been applied successfully by itself or as a model to study a range of soft matter systems 8–10 , and recently extended to include representations of amino acid residues and rigid molecules such as cholesterol. A brief description of their coarse-graining approach follows here - within this model 3-5 heavy atoms are represented by a single interaction site or "bead" (one water bead represents three water molecules). The system in this study consisted of a single NP, a non-ionic surfactant bilayer consisting of 2000 poly(oxyethylene) (denoted C12 E2 ) surfactant molecules and 57312 water beads. The surfactant molecule was described using four CG bead types: OA (-CH2 OH) and EO (CH2 OCH2 -) which represent the hydrophilic head group and CM (-CH2 CH2 CH2 -) and CT2 (CH3 CH2 CH2 -), which represent the hydrophobic tail groups. A schematic of the atomistic to CG mapping is given in Figure 1. The non-bonded interactions between the surfactant molecules and between the surfactant molecules and water were described using the Lennard-Jones-like formulae ULJ12−4 (r) = ULJ9−6 (r) = 2| √   3 3ε  σ 12  σ 4 − 2 r r (1a)   27ε  σ 9  σ 6 − 4 r r (1b) 1–7 Fig. 2 The side and top views of the polymeric bilayer we will be using for this simulation. The dimensions of the expanded bilayer (a) are 16 × 16 nm2 , while the thickness of the bilayer is approximately 3.6 nm where ε is the well depth (the numerical factors in the prefactor are chosen such that the minimum value of the potential is −ε ) and U(σ ) = 0, following a standard Lennard-Jones model. The first (LJ12−4 ) form is only used for interactions involving water (W) beads. The intra-molecular interactions were described using a molecular mechanics force field: Ubond = kℓ (ℓ − ℓ0 )2 (2a) Uangle = kθ (θ − θ0 )2 (2b) where kℓ and kθ are the bond stretching and bending force constants and ℓ0 and θ0 are the equilibrium bond lengths and angles respectively. Parameters for both the bonded and non-bonded interactions are taken from parameters published by Shinoda et al 11 and are tabulated in supplementary information for completeness. Rather than explicitly model the NP as a rigid collection of interaction sites, the particles are modelled as spheres composed of uniformly distributed interaction sites i.e. the particle was treated as a surface-to-molecule potential, which is the ap- proach taken by Chiu et al, which itself is based upon Hamaker’s model for modelling the non-bonding potential of surfaces 12,13 . Using Hamaker’s model allows us to concentrate on the bulk property of the NP rather than its individual constituents, and since we are primarily interested in the effect of the surface properties of the NP, it fits this study. To represent the hydrophilic and hydrophobic NPs, the interaction sites were taken to be of type W and CT2 respectively. The interactions between the NP and CG beads were derived from integrating the interaction potential between a CG bead and an interaction site within the NP over the particle’s volume. Assuming that the interaction between a CG bead and single interaction site can be described through the modified LJ functions above (Equations 1a and 1b) the interaction between a NP and a CG bead may be written as: √ 2 3πρεσ 12 R3 5R6 + 45R4 r2 + 63R2 r4 + 15r6 12−4 UNP (r, R) = 15 (r2 − R2 )9 − √   3πρεσ 4 R 3 3πρεσ 4 r−R + ln 2r r+R r 2 − R2 9−6 (r, R) = 9πρεσ 9 R3 UNP (3a) 3R4 + 42R2 r4 + 35r4 9πρεσ 6 R3 − 2 (3b) r − R2 35r(r2 − R2 )6 where ρ is the density of the NP and R is the NP radius. The density of the NP was taken to be 1.0 kg m−3 (based on the density of water). Full derivations of these formulae are given in the supplementary information. The parameters ε and σ are the Van der Waals parameters for the W (hydrophilic) or CM (hydrophobic) CG beads. As well as purely hydrophilic and hydrophobic NPs we also consider a NP of mixed hydrophobicity; in that case the interaction potential between the NP and solvent bead is taken as a simple interpolation of the hydrophilic and hydrophobic potentials:  phile phobe phile(r)=(x)UNP (r)+(1−x)UNP   :x=1   UNP NPproperty = phobe(r)=(x)U phile (r)+(1−x)U phobe NP NP UNP   phile phobe   U mixed(r)=(x)UNP (r)+(1−x)UNP NP :x=0 : x = 0.5 where 0 ≤ x ≤ 1. The plots of the interaction potentials of each NP are provided in the supplementary information. Two R values were used for creating the NPs - R = 1.0 and 2.0 nm. The 1.0 nm radius NP represents a NP that is comparable to the tailgroup (CM and CT2) length of the amphiphile, while the 2.0 nm radius NP represents a NP with a width that is comparable to the thickness of the bilayer. The schematic for each type of NP is shown in Figure 3. Initial configuration of the bilayer All simulations were performed using the LAMMPS simulation package 14 . To generate the initial configuration of the bilayer, 3500 CG-beads of the C12 E2 were used to arrange a bilayer in a simulation box ranging from dimensions 80 Å × 80 Å × 200 Å with the rest of the box being filled with 14328 water beads. This system was equilibrated using a NVT ensemble for 1 ns simulation time with increasing timesteps to ensure the overlaps be- (a) Hydrophilic (b) Mixed (c) Hydrophobic Fig. 3 Color Schematic for each NP. Two values of radiuses (R) for each NP was designed - R = 1.0 and 2.0 nm tween molecules were resolved and to ensure the system did not fluctuate substantially. To reduce the likelihood of finite size effects, the simulation box was replicated once in each of the x and y directions, meaning that the simulation box was expanded from 80Å × 80 Å × 200 Å to 160 Å × 160Å × 200Å. 14000 CG beads of the C12 E2 were used to construct the larger bilayer, which equals to 2000 molecules. This new system was also subjected to 1 ns of NVT equilibration in the same manner as the smaller bilayer. To insert the NP, the water beads were removed, the simulation box dimension in the direction normal to the bilayer (z-direction) was set to 200 Å and the NP was inserted into the system 100 Å above the centre of mass of the bilayer followed by re-solvation. Simulation setup for NPT production runs The system was simulated in the NPT -ensemble, with temperature and pressure set to 303 K and 1 atm respectively. The temperature and pressure were controlled using a Nosé-Hoover thermostat and barostat 15–17 , both with relaxation times of 0.2 ps. The equations of motion were integrated using the rRESPA multipletimestep algorithm 18 with a 2.0 fs inner (bonded) and 10.0 fs (non-bonded) outer timesteps. Non-bonded interactions were truncated at rcut = 15 Å. Periodic Boundary Conditions (PBCs) were set in the x and y directions. A production run of 200 ns was run to ensure that the surfactant bilayer had reached consistency in it’s properties with past experimental and simulation results. 3 Free energy calculations - Umbrella Sampling (US) In order to get a complete picture of NP translocation energetics, we have implemented the Umbrella sampling method (US) 19 . A spring constant of k = 50.0 kcal mol−1 was used to anchor the NP along the bilayer normal from 0 Å to 30 Å in 1 Å intervals, resulting in 30 windows for each nanoparticle. The coordinate distributions of each window was checked for overlap. The implementation of the weighted-histogram analysis method (WHAM) program developed by Grossfield 20–22 was used to extract the final potential mean force (PMF) of the trajectory. The probability distribution of the umbrella histograms are shown in the supporting information. 4 Results 4.1 Simulation Analysis The density profile of the surfactant (Figure 4) shows that the average thickness of each surfactant leaflet is approximately 17-18 1–7 | 3 NP type Hydrophilic Hydrophilic Mixed Mixed Hydrophobic Hydrophobic R (nm) 1.0 2.0 1.0 2.0 1.0 2.0 Number of simulations 20 20 20 20 20 20 (a) Hydrophilic 20 z = 20 Å ∆G (kcal mol−1 ) Index 1 2 3 4 5 6 z = 10 Å R = 1 nm R = 2 nm 10 0 -10 -20 z = 0 Å -30 0 5 10 15 20 Table 1 The catalogue of US simulations undertaken in this work (b) Mixed 20 R = 1 nm R = 2 nm ∆G (kcal mol−1 ) z = 20 Å z = 10 Å 10 0 -10 -20 z = 0 Å -30 0 5 10 15 20 (c) Hydrophobic 20 10 ∆G (kcal mol−1 ) z = 20 Å z = 10 Å Fig. 4 The density profile of each CG beads in the bilayer. The EO and OA beads represents those within the headgroup region of the surfactant, while the CM and CT2 beads represent those composing the hydrophobic tailgroups. W represents the bulk water above and below the bilayer on the bilayer normal Å, which amounts to a total thickness of 36 Å. This is approximated from the bilayer centre to the peak of the EO bead. This thickness is consistent with those shown in previous simulation and experimental data 23–26 . Additional studies looking at the mean square displacement (MSD) and area per lipid (APL) are shown in in the supplementary information. 4.2 Analysis of the NP trajectories and the free energy profiles Figure 5(a) shows the trajectory snapshots for the R = 1.0 nm and R = 2.0 nm simulations, with the corresponding free energy profile. In the R = 1.0 nm case, we see a free energy trough at approximately z = 10 Å, where we had observed the adsorption of the NP onto the surfactant surface. From 10 - 0 Å distance, we see a levelling off of the free energy which corresponds to the NP inducing curvature in the surfactant bilayer. In the R = 2.0 nm case, we see a significant increase in free energy between z = 15 - 12.5 Å distance of approximately 13.0 kcal mol−1 , which corresponds to the significant bending of the bilayer and pore formation around the NP as it nears the center of of the bilayer. Figure 5(b) shows the trajectory snapshots for the R = 1.0 nm and R = 2.0 nm simulations, along with the corresponding free energy profile. In the R = 1.0 nm case, we see a slight energetic trough at approximately z = 12.5 Å, which is followed by a decrease in the free energy from z = 12.5 - 0 Å, with an overall free energy change of -25.0 kcal mol−1 . In the R = 2.0 nm case, we see the absence of this free energy trough, where instead we see a slight 4| 1–7 R = 1 nm R = 2 nm 0 -10 -20 -30 -40 -50 z = 0 Å -60 0 5 10 15 20 Distance from bilayer center (nm) R = 1.0 nm R = 2.0 nm Fig. 5 Snapshots of the US simulations for the hydrophilic, mixed and hydrophobic NPs at z = 0, 10 and 20 Å in the bilayer normal for the R = 1.0 and 2.0 nm, with the corresponding free energy profiles shown in (a), (b) and (c) respectively increase of 5.0 kcal mol−1 , at z = 10 Å, but otherwise showing a ‘flattened’ free energy surface, which is in great contrast to the hydrophilic NP example. Figure 5(c) shows the trajectory snapshots for the R = 1.0 nm and R = 2.0 nm simulations, with the corresponding free energy profile. In this case, we see a consistent wrapping of the NPs into the hydrophobic tailgroups of the bilayer, which is demonstrated by the similar free energy profiles, where we see an overall change in the range of -40.0 - -50.0 kcal mol−1 for both NPs as the it appoaches the bilayer interior. 5 Discussion Due to the nature of the construction of the NP, there is no direct experimental equivalent that we can compare the PMF to in a quantitative manner. However, by comparing the properties of the system with respect to the NP position, comparisons with published literature can be made. Similar approaches at fine-graining the NP property has been done; using discontinuous molecular dynamics (DMD), NP studies done by Curtis et al 4 illustrated a radius limit for the hydrophilic NP to change from a bilayer ‘embedding’ phase (where the bilayer curves around the bilayer) to a bilayer ‘wrapping’ phase (where the bilayer completely wraps around the bilayer). Here, the range of radius measured in our study (R = 1.0 - 2.0 nm) falls within the ‘embedding’ phase, which is consistent with the DMD study, which is consistent with the phenomena observed for the hydrophilic and mixed Hydrophilic ∆G (kcal mol−1) 0 -100 -200 -300 -400 -500 R = 1.0 nm R = 2.0 nm Mixed 0 R = 1.0 nm R = 2.0 nm -50 -100 -150 -200 Hydrophobic 0 -50 -100 -150 -200 -250 -300 R = 1.0 nm R = 2.0 nm 0 2 4 Time (ns) 6 Fig. 6 The total interaction energy of the surfactants to each of the NPs. (a), (b) and (c) show the hydrophilic, mixed and hydrophobic examples respectively, for the R = 1.0 and 2.0 nm NPs. The interaction energy was computed for the z = 0 nm US simulation, where the NPs are held in the center of the bilayer. NPs. To distinguish and analyze the mechanistic difference between each NP, the magnitude of change in the free energy profile and the shape of the PMF will be the main factors of consideration in this discussion. This is consistent with simulations of hydrophilic NPs done by Rocha et al 27 and Gu 28 , amongst others, where a higher charge density on the NP (which corresponded to a greater hydrophilicity) correlated with a stronger adsorption on the bilayer surface. Simulations of hydrophilic NPs by Gu indicated that there is competition between the favourable electrostatic interaction between the charged groups in the bilayer and the positive/negatively charged NP, and the unfavourable energy associated with bending the bilayer. Simulations of acylated and un-acylated dendrimers through a lipid bilayer 29 show that the unacylated dendrimer acts analogous to the hydrophobic NP, as shown by the insertion of the dendrimer into the tailgroups; the acylated version shows adsorption on the surface of a 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) bilayer, analogous to the hydrophilic NPs simulated. The model we have used does not include charged components within the surfactant groups nor implement surface charges on the NP surface. However, the association of the hydrophilic headgroups around the hydrophilic/mixed NPs follows the pattern of previous studies. This ‘flattening’ effect is consistent with recent studies demonstrated by Su et al 30 , who suggested that, at mixed hydrophobicity, the NP is distributed between the head and tail groups, and in-between the bilayer leaflet layers, all accessible through a flattening of the free energy profile across the bilayer normal. Su also suggested that the barrier entailing the tail end of the flattened PMF (shown in 0 - 5 Å in the R = 1.0 nm mixed NP profiles) is associated with the entropic cost of inducing order within the bilayer tailgroups. In that case, the barrier was shown to be within 10’s of kT, which is clearly smaller than the barrier observed in this example, (which is in the range of 50 - 100 kcal mol−1 ). It may be the case that the NP model used in this example overestimates the barrier. Simulations with intemediate type dendrimers by Hwankyu et al 29 demonstrated that the mixed dendrimers adsorbed onto the bilayer surface (which is a ‘acylated’ trait) with strands of the dendrimer penetrating inside the bilayer, which suggests that such NPs may allow the mixed property to ‘latch-on’ to the hydrophilic surface, and subsequently penetrating through the bilayer normal through increased hydrophobic contacts. Another mechanism of insertion of mixed type NPs was demonstrated using striated NPs as simulated by Yinfeng et al 31 , who demonstrated that the rotational dynamics of hydrphobic/hydrophilic striped NPs can rotate while inside the bilayer to minimise the interaction with the hydrophobic interior, which again, leads to a ‘flattening’ of the free energy profile. In other examples, Zhang et al 32 performed DPD membrane translocation simulations of NPs with fully hydrophilic (WNP), hydrophobic (ONP) and a random mixture of hydrophilic and hydrophobic points to represent a NP of mixed hydrophobicity (RNP). Here, the resistance force (the combined repulsion force from the hydrophilic headgroups and hydrophobic tailgroups of the bilayer) acting against the RNP was minimal - it has been hypothesised that the randomly distributed hydrophilic forces counteract against the attractive hydrophobic forces near the tailgroup region. Furthermore, it may be the case that larger aggregations will undergo cooperative mechanisms of translocation, which has not been covered with model and will be a subject for a future work. This observation was followed by a flattening of the PMF near the membrane region, which is entirely consistent with the observation in this study in the case of R = 1.0 nm mixed NPs. Another example of this type of direct penetration through membranes with mixed hydrophobicity comes with striped NPs, as designed from monolayer-protected NPs 33,34 , which suggested that homogenous patterning of hydrophilic/hydrophobic beads on a NP surface enables the passive translocation of NPs through a lipid bilayer. In contrast to the smaller examples, the R = 2.0 nm hydrophilic and mixed NPs clearly shows the formation of holes or ‘doors’ on approach to the bilayer surface, as demonstrated in the US simulations. The formation of doors has been speculated upon in the past; MD simulations by Lin et al 35 illustrated that the interaction of cationic alkanethiol-ligand AuNPs with the transmembrane potential as the driver of the movement of the NP, with the formation of similar doors on approach to the bilayer surface - it is thought that this mechanism is responsible for the rapid intake of cationic NPs and cell-penetrating peptides (CPPs) 36 . Here, as we have not included the ionic potential as a factor in the experiment, an alternate mechanism for direct penetration in a non-ionic case is necessary. We see that for R = 1.0 nm, the NP does not rupture the bilayer when placed in the center of the bilayer at z = 0 Å. With the R = 2.0 nm example, we see that both the hydrophilic and mixed NPs induce deformation of the bilayer and subsequent formation of ‘doors’ on approach to the bilayer. However, the shape of the bilayer around the rupture differs; in the hydrophilic example, the surfactants around 1–7 | 5 Hydrophilic 10 8 6 4 2 0 0 1 2 3 4 5 6 7 8 9 10 6 7 8 9 10 Mixed 7 6 5 4 3 2 1 0 0 1 2 3 4 5 Hydrophobic 2 1.5 1 Conflicts of interest 0.5 0 There are no conflicts to declare. 0 1 2 3 4 5 6 7 8 9 10 Fig. 7 Water density around the R = 1.0 and 2.0 nm hydrophilic, mixed and hydrophoic NP, shown in a radius of 10 nm around each NP at z = 0 Å, illustrated with the RDF density profile of water around each NP. the NP become ‘folded’, while in the mixed example we observe a clear hole around the NP lateral plane, which indicates that the energetically favourable interactions with the hydrophobic tailgroups clearly affects rupture geometry. To ascertain the total force contribution on the NP surface, we have computed the total interaction energy between the NP and the bilayer. This is shown in Figure 6. The hydrophilic and mixed examples show similar trends, but we see a diminished interaction energy with the surfactants with the mixed NP types, where we see a interaction energy of of approximately -100 kcal mol−1 , compared to the limit of -300 kcal mol−1 we observe with the hydrophilic example, as would be expected from our model. When comparing the hydrophobic example with the hydrophilic/mixed example, we see a lower fluctuation in the interaction energy of the hydrophobic NPs, as they are embedded within the tailgroup (CM and CT2) regions of the bilayer, wheras the hydrophilic/mixed NPs show tenuous interaction with the hydrophilic headgroups, as well as the surrounding water solvent, as shown with the water density profiles around each NPs, shown in Figure 7. 6 Conclusion We have adapted the Hamaker-like continuum potentials to construct 3 different types of NPs; a hydrophilic, hydrophobic and mixed particle of two different sizes; R = 1.0 and 2.0 nm, and analyzed the translocation of each NP through a neutral C12 E2 bilayer using US simulations with the bilayer normal as the reaction coordinate. By analyzing the free energy profile and surfactant configuration around the NP, we have attempted to characterize the mechanism for each hydrophobicity. We see that the hydrophilic NP, at R = 1.0 nm, the favourable interaction between the hydrophilic headgroups of the surfactants and the NP surface 6| overcomes the energetic cost of forming the curvature around the NP, while for the R = 2.0 nm, we observe that the NP forms a ‘hole’ upon approach to the surfactant surface. In the mixed NP case, we observe a flattening of the PMF in the surfactant region in the cases of the R = 1.0 nm, while with the R = 2.0 nm case, we see the formation of a hole on approach to the surfactant surface. Hence, we demonstrate that the flattening of the free energy profile accross the membrane bilayer is a general property for semi-hydrophobic NPs, but may still be subject to size-dependent effects. While numerous challenges to the clinical applications of NPs remain, this study shows that by tuning the hydrophobicity of the NP surface, it may be possible to design targeted NP delivery mechanisms that can strictly bypass the bilayer - as the lifetime within the bilayer may be short, the concerns regarding nanotoxicity and complications regarding accumulation of NPs within the interior of the membrane. 1–7 Acknowledgements We would like the acknowledge the Center for Scientific Computing (CSC) at the University of Warwick for computing resources required for this work. Notes and references 1 N. Lewinski, V. Colvin and R. Drezek, Small, 2008, 4, 26–49. 2 M. M. Bailey, E. M. Gorman, E. J. Munson and C. J. Berkland, Langmuir, 2008, 23, 13614 – 13620. 3 F. Liu, D. Wu, R. D. Kamm and K. Chen, Biochimica et Biophysica Acta, 2013, 1828, 1667 – 1673. 4 E. M. Curtis, A. H. Bahrami, T. R. Weikl and C. K. Hall, Nanoscale, 2015, 7, 10799 – 10808. 5 S. Bandyopadhyay, M. Tarek, M. L. Lynch and M. L. Klein, Langmuir, 2000, 16, 942 – 946. 6 W. Shinoda, R. DeVane and M. L. Klein, Soft Matter, 2008, 4, 2454–2462. 7 S. Seo and W. Shinoda, J. Chem. Theory Comput., 2019, 15, 26–49. 8 D. Alemani, F. Collu, M. Cascalla and M. D. Peraro, The Journal of Chemical Theory and Computation, 2010, 6, 315 – 324. 9 R. Devane, M. L. Klein, C. Chiu, S. O. Nielson, W. Shinoda and P. B. Moore, J. Phys. Chem. B, 2002, 296, 6386 – 6393. 10 D. N. LeBard, B. G. Levine, P. Mertmann, S. A. Barr, A. Jusufi, S. Sanders, M. L. Klein and A. Z. Panagiotopoulos, Soft Matter, 2012, 8, 2385 – 2397. 11 W. Shinoda, R. DeVane and M. L. Klein, Molecular Simulation, 2007, 33, 27–28. 12 C. Chi-Cheng, P. B. Moore, W. Shinoda and S. O. Nielson, The Journal of Chemical Physics, 2009, 131, 244706. 13 H. C. Hamaker, Physica, 1937, 4, 1058 – 1072. 14 S. Plimpton, Journal of Computational Physics, 1995, 117, 1– 19. 15 S. Nose, The Journal of Chemical Physics, 1984, 81, 511. 16 M. E. Tuckermann, C. J. Mundy and G. J. Martyna, Journal of Chemical Physics, 2001, 115, 1678. 17 M. E. Tuckermann, C. H. Mundy and G. J. Martyna, Europhys . Lett, 1999, 45, 149. 18 M. Tuckermann, B. J. Berne and G. J. Martyna, J. Chem. Phys, 1992, 97, 1990 – 2001. 19 G. M. Torrie and J. P. Valleau, Chem Phys Lett, 1974, 28, 578– 581. 20 A. Grossfield, WHAM: the weighted histogram analysis method, http://membrane.urmc.rochester.edu/ content/wham. 21 S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem, 1992, 13, 1011 – 1021. 22 M. Sousaile and B. Roux, Comput. Phys. Commun, 2001, 135, 40 – 57. 23 W. Shinoda and M. L. Klein, Soft Matter, 2008, 4, 2454–2462. 24 S. Bandyopadhyay and J. Chanda, Langmuir, 2003, 19, 10443–10448. 25 S. Bandyopadhyay and J. Chanda, J. Phys. Chem. B, 2006, 110, 10443–10448. 26 S. S. Fumari and G. Rapp, J. Phys. Chem. B, 1997, 101, 732 – 27 28 29 30 31 32 33 34 35 36 739. E. L. da Rocha, G. F. Caramori and C. R. Rambo, Phys. Chem. Chem. Phys, 2013, 15, 2282 – 2290. Y. Li and N. Gu, Journal of Physical Chemistry B, 2010, 114, 2749 – 2754. L. Hwankyu and R. G. Larson, Journal of Physical Chemistry B, 2006, 110, 18204 – 18211. C. Su, H. Merlitz, H. Rabbel and J. Sommer, Phys. Chem. Lett, 2017, 8, 4069 – 4076. L. Yinfeng, L. Xuejin, L. Zhongha and H. Gao, Nanoscale, 2012, 4, 3768. H. Zhang, Q. Ji, C. Huang, S. Zhang, B. Yuan, K. Yang and Y. Ma, Scientific reports, 2015, 5, 10525. A. Verma, O. Uzun, Y. Hu, Y. Hu, H. Han, N. Watson, S. Chen, D. J. Irvine and F. Stellachi, Nature Materials, 2008, 7, 588 – 595. P. Gkeka, L. Sarkisov and P. Angelikopoulos, Journal of Physical Chemistry Letters, 2013, 4, 1907 – 1912. J. Lin and A. Katz, ACS Nano, 2013, 7, 10799 – 10808. P. Nativo, I. A. Prior and M. Brust, ACS Nano, 2008, 2, 1639 – 1644. 1–7 | 7