Molecular Dynamics Simulations of the Silica â Cell Membrane

Subscriber access provided by BUPMC - Bibliothèque Universitaire Pierre et Marie Curie

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Molecular Dynamics Simulations of the Silica – Cell Membrane Interaction: Insights on Biomineralization and Nanotoxicity Massimo Delle Piane, Sebastian Potthoff, C. Jeffrey Brinker, and Lucio Colombi Ciacchi J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b04537 • Publication Date (Web): 23 Aug 2018 Downloaded from http://pubs.acs.org on August 28, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Dynamics Simulations of the Silica – Cell Membrane Interaction: Insights on Biomineralization and Nanotoxicity Massimo Delle Piane†*, Sebastian Potthoff†, C. Jeffrey Brinker+, Lucio Colombi Ciacchi† †

Hybrid Materials and Interface Group, University of Bremen, Faculty of Production Engineering,

Bremen Center for Computational Material Science, Center for Environmental Research and Sustainable Technology (UFT) and MAPEX Center for Materials and Processes, Am Fallturm 1 28359, Bremen, Germany +

Advanced Materials Laboratory, Sandia National Laboratories, Albuquerque, NM 87106, United

States

1 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT The interaction of silica (SiO2) with biological systems is complex and contradictory. On the one hand, silica is at the basis of several biomineralization processes (e.g. in sponges). On the other hand, silica nanoparticles and dust may lead to silicosis and, at the cellular level, hemolysis. These toxic responses are strongly dependent on the silica polymorph and their root causes are still under debate. Both silica biomineralization and silica-induced nanotoxicity could be related to similar mechanisms of molecular recognition between the cellular membranes and the surface of the SiO2 particles. Based on this hypothesis, we employed classical molecular dynamics simulations, coupled to advanced-sampling techniques, to achieve an atomistic picture of the interactions between different types of silica nanoparticles and the membrane of erythrocytes. Our predicted free-energy profiles associated with membrane crossing give no evidence for segregation of nanoparticles at the membrane/water interface, irrespective of their Si nuclearity, structure and charge. The associated molecular trajectories, however, are suggestive of a possible direct translocation mechanism, in which silica nanoclusters elicit both local and large-scale effects on the membrane dynamics and stability. This gives hints on possible pathways for silica nanotoxicity based on nanoparticle-induced membrane perforation.

*

Corresponding author: [email protected] - +49 421 218 64579


Page 2 of 54

Page 3 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION The contact between inorganic materials and living tissues is of extraordinary relevance in biomedicine, pharmacology and biotechnology. Among the most important inorganic material classes, silica (silicon dioxide, SiO2) polymorphs play a key role, because they are chemically inert, structurally robust, easy to prepare and characterize.1 However, the interaction between silica and biological systems is complex, contradictory and to a large extent poorly understood. On the one hand, silica is at the basis of several biomineralization processes: organisms such as sponges or diatoms build up silica skeletons via the bio-templated polycondensation of silicic acid (Si(OH)4)2 dissolved in sea water into SiO2 structures.3 The accurate simulation of the silica polycondensation reactions and the generation of stable particles has been extensively studied in literature.4–6 Artificial, lab-made, biomineralization has also been described as a possible way to transfer cellular architectures into inorganic materials by means of silica biocomposites.7,8 On the other hand, inhalation of silica dust over extended periods of time leads to silicosis, a lung-related disease. At the cellular level, the contact of some types of silica with the membrane of erythrocytes is known to cause its rupture and cellular death (hemolysis).9–12 Since silica nanoparticles have been extensively proposed as carriers in nanomedicine,13 assessing their potential toxicity is of the utmost importance and a primary aim of fundamental research. Adding complexity to the issue, toxicity is known to be strongly dependent on the type of silica polymorph and its precise surface chemistry. For instance, the US FDA regards amorphous silica as GRAS – generally recognized as safe, while crystalline silica nanoparticles can be cytotoxic.10 However, so-called fumed silica (i.e. amorphous SiO2 nanoparticles produced by gas-phase pyrolysis) does induce hemolysis and is also cytotoxic to certain cellular strains. In contrast colloidal silicas such as Stöber silica14 and so-called mesoporous silica nanoparticles are much less harmful to eukaryotic cells and result in significantly less hemolysis than fumed silica.10,12 Notably, the sol-gel synthesis of Stöber silica nanoparticles happens via a monomer-addition process similar to naturally occurring silica biomineralization.14,15



The reasons for these markedly different behaviors are still under debate. Silica toxicity has been related to the potential of the particle to generate reactive oxygen species, but it has also been connected to the pattern of hydrogen bonded versus isolated silanols on the amorphous silica surface, in relation to silica-cellular membranes interactions..9,10,12 Interestingly, also the biomineralization capability of biomolecular templates is tightly related to a favorable match between the organic promoter and the growing mineral phase. The composition (in terms of ratio between the different phospholipid classes, cholesterol content, glycosylation pattern and exposed proteins) of the cell membrane16 most probably impacts on how a silica nanoparticle interacts with it.17 The work presented here is based on the hypothesis that both silica biomineralization and silicainduced rupture of erythrocyte membranes could be related to conceptually similar mechanisms of molecular recognition between the mutually interacting organic and inorganic phases. Due to the inherent complexity of biological systems, the experimental investigation of the details of these processes is unattainable through current methodologies. Experiments have been recently performed to study the contact between silica nanoparticles and model lipid membranes,17 but still no direct evidence of the mechanism of interaction has been found. Furthermore, no estimation of the energetics of the silica-membrane interface has been provided so far. This is especially true for the case of small silica clusters in the initial stage of Si(OH)4 polycondensation and biomineralization. Simulations, on the other hand, have been successfully demonstrated to be able to act as a “virtual microscope”, giving insights on the processes occurring at the organic/inorganic interface, with many examples regarding biomolecular interaction with silica surfaces.18–21 The interaction between inorganic nanoparticles and lipid membranes has been studied with increasing frequency by means of molecular dynamics (MD) simulations in the past years, but most investigations are performed at the coarse-grained level and employing highly simplified membrane models of homogeneous lipid composition.22–24 Very few studies are conducted with explicit inclusion of all atoms in the


Page 4 of 54

Page 5 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulation cells. Notable is a recent investigation of functionalized gold nanoparticles with phospholipid bilayers of asymmetric composition.25 Both naturally occurring2 (e.g. diatoms) and in-vitro7,8 silica biomineralization processes lead to the formation of linear, cyclic and branched chains, which further condense via the same general mechanism to a highly hydrophilic form of silica.26 For the present paper, we decided to focus on the first steps of these processes, by studying the interaction of small silica nanoclusters with the cell membrane. We aim at achieving an atomistic picture of the interactions taking place between different nanoclusters and the membrane of erythrocytes by means of all-atom molecular dynamics simulations, coupled to techniques that enable a prediction of free energy profiles across hybrid interfaces (well-tempered metadynamics), to gain a quantification of the energetics of the silicamembrane interaction. The ultimate goal is to shed some light on both the mechanism of silica biomineralization and the reasons behind the hemolytic power of some forms of silica.

COMPUTATIONAL DETAILS Silica models and force field Models for small silica clusters of increasing size and complexity were generated (Figure 1). They will be referred to using the nomenclature SIn, where n is the number of SiO4 units: SI1, SI2, SI6b and SI9b, where b highlights that the silica tetrahedra are arranged in branched fashion. The clusters were extracted from an amorphous silica bulk and the dangling bonds were saturated with hydrogens to form silanol groups. For the interactions of lipids and water with the silica particles, we employed the parameter set for silica by Butenuth et al.27, that, although originally validated for the water-silica interaction, has been proven reliable in describing silica-biomolecule interactions in several recent works.18,20,21,28 Notably, the partial atomic charges and Lennard-Jones terms in this force field do not depend on the silica cluster nuclearity and can be used without modification for single Si(OH)4 molecules or larger clusters, up to whole nanoparticles. The intramolecular Si−Si and Si−O interactions were described



starting from the parameters by Demiralp et al.29, as modified by Meißner et al.28, that describe these pairs with non-bonded potentials of the Morse type. Since the GROMACS code does not provide a simple implementation of a Morse potential for non-bonded interactions, a LJ function with the same dispersion parameters was employed instead. Since the different curvature in the minimum generates over/under-binding for close atoms, harmonic Si-O bond terms and Si-O-Si angle terms were added to compensate for these errors. The employed parameter set is reported in Table S1 of the Supporting Information. This excludes possible exchanges of H+ or OH- ions between the particles and their environment. Different protonation states, corresponding to the presence of net charges on the particles, need to be adjusted prior to each simulation. Structural properties for bulk amorphous SiO2 computed with our modified potential are in good agreement both with experimental data and simulations employing the original force field (Figure S1).

Figure 1. Model silica particles used in the present work. SIn models represent small silica clusters containing n SiO4 units, arranged in a linear or branched (b) fashion.

Membrane models and force field


Page 6 of 54

Page 7 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Three different lipid bilayer models were assembled by using the freely accessible CHARMM-GUI Membrane Builder (www.charmm-gui.org):30–32 a) a homogeneous bilayer containing a single lipid species, 1-palmitoyl-2-oleoyl-sn-glycero-3phosphocholine (POPC), with identical leaflets, each containing 160 lipid molecules (Figure S2); b) a mixed bilayer with POPC and cholesterol (hereafter named POPC_CH), still symmetrical, but in which 32 POPC lipids were substituted with cholesterol molecules in both the inner and outer leaflet, resulting in 20% of all lipids of the membrane, consistent with the average cholesterol content of red blood cells33 (Figure S3); c) a heterogeneous lipid bilayer representing the typical composition of the erythrocyte membrane bilayer with asymmetric leaflets (hereafter named RBC, Figure S4). To reduce the complexity and the computational cost of our simulations, among the circa 160 different types of lipid species that constitute a real erythrocyte membrane,34 we resorted to considering only those species representing at least 2% in the typical composition of the bilayer.16,33–35 18 lipid species (belonging to the phosphatidylcholines,

PC,

phosphatidylserines,

PS,

phosphatidylethanolamines,

PE,

and

sphingomyelins, SM, classes) were selected, representing 64% of all lipids found in a real erythrocyte’s bilayer, and asymmetrically distributed among the inner and outer leaflets according to Virtanen et al.16 The resulting model consists of 128 lipids per leaflet plus 32 and 39 cholesterol molecules in the inner and outer leaflet, respectively. The final composition of the heterogeneous system can be found in Table S2 of the Supporting Information. In all models, the membrane lies on the xy plane, with the z axis normal to the lipid bilayer. The system size was chosen aiming at a surface area large enough to represent the erythrocyte bilayer and provide good statistics of the membrane behavior, and small enough to not include a membrane protein. This led to a square area with a side length of around 95 Å, with some variability between the models. Each system was built taking into account the physiological environment of an erythrocyte cell. This means that water molecules and 150 mM KCl were added to the simulation box. In order to guarantee neutrality for the heterogeneous RBC system, counter ions (K+) were



further added to balance the negative charge of the anionic lipids. The resulting systems included 90956, 72862, 88826 atoms for the POPC, POPC_CH, RBC models, respectively. The lipid bilayer was described using the CHARMM36 parameter set,36 that reproduces very accurately the experimental lipid properties and was recently found to perform well in the prediction of free energy barriers for the crossing of small molecules across lipid membranes.37 Furthermore, the CHARMM36 formulation allows to easily pair it with the silica force field by employing Lorentz−Berthelot combination rules to construct the Lennard-Jones (LJ) interactions between the atomic species.27 Water was described using the CHARMM-modified TIP3P water model (mTIP3P).38 Both the calculated area per lipid (i.e. the specific area that a lipid occupies) and the thickness of the lipid bilayers, averaged along the final 80 ns of production run (vide infra) relates well with the experiments39,40,41 (Table S3 of the Supporting Information) and correlate, as expected, with the cholesterol content. The good agreement is an expected feature of the well parameterized CHARMM36 lipid force field.36

Molecular Dynamics (MD) simulations The GROMACS code, in its 5.1.2 version, was used for all the MD simulations.42 For every membrane model, an energy minimization was first performed using the steepest descent algorithm with a tolerance of 1000 kJ mol-1 nm-1. Equilibration was then achieved using the multistep approach proposed by CHARMM-GUI:31 NVT dynamics is used for the first 4 equilibration steps (corresponding to a total of 100 ps) and NPT dynamics is used for the rest of the equilibration (5 more steps, or further 1 ns). Along the equilibration, restraints on both atomic positions and dihedrals of the lipid molecules are gradually decreased, so that at the end all atoms are completely free to move. The system temperature was set to 310 K to represent the physiological body temperature and the reference pressure was set to 1 bar, with compressibility set to 4.5 10-5 bar-1. The Berendsen thermostat and barostat were used for the equilibration, with time


Page 8 of 54

Page 9 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


constants of 1 and 5 ps, respectively. The time step was 1 fs. For the barostat, a semi-isotropic scaling behavior, treating the z axis (normal to the membrane) differently from the x and y axis, was chosen. Subsequent production MD runs for each system were performed for at least 100 ns. In this case a Nosé-Hoover thermostat, with a time constant of 1 ps, was used together with the ParrinelloRahman barostat, with a time constant of 5 ps, and semi-isotropic behavior. In production runs the time step was set to 2 fs. Silica models were also equilibrated in water at 310 K with a sequence of minimization-NVT-NPT steps. The starting structures for the free energy silica-membrane simulations (vide infra) were then built by inserting the silica clusters in the middle of the water phase above the membrane, removing all water molecules within 2 Å from any silica atom. Then, to thermalize the system, the multi-step equilibration was performed again, keeping the silica cluster in its position. Production runs (with or without the addition of the enhanced sampling methods described below) were then performed with the same parameters of the pure water-membrane systems for at least 300 ns. In all cases the non-bonded interaction cutoff was set to 1.2 nm and all hydrogen atoms in the system were constrained via the LINCS algorithm.43 The relative motion of the center of mass of the membrane with respect to the rest of the system was removed every 100 steps to ensure that the membrane does not drift out of the simulation box.



Figure 2. Schematics of how the different free energy differences of Table 1 (crossing barrier, ∆G‡, free energy of partition, ∆Gpart, free energy of adsorption, ∆Gads) have been computed. An exemplary free energy profile (from Figure 3c) is reported (red line), superimposed to the mass densities of both water (blue line) and lipids (black line) averaged along the corresponding simulation. Both free energy and density are plotted as a function of the collective variable CV, i.e. the fractional coordinate along the direction (z) normal to the membrane. The background is colored according to the different regions: water phase (cyan), where the water density is maximal, lipid phase (grey) where the water density is zero, interface (yellow) where the water and lipid densities crosses.

Well-tempered Metadynamics We resorted on methods developed to enhance the sampling of MD simulations, that make it possible to target free energies.44 In particular, our simulations were based on the well-tempered Metadynamics method (wt-metaD) introduced by Barducci et al.45 as implemented in PLUMED (version 2.2.3).46 Briefly, adaptive bias potentials were added during the course of the MD runs according to the metadynamics scheme47 in the well-tempered ensemble45 using a bias factor of 10. Gaussian hills with initial height of 1.0 kJ mol-1 were added every 0.6 ps. As a collective variable


Page 10 of 54

Page 11 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(CV) we chose the fractional spatial coordinate of the center of mass of the silica cluster in the direction normal to the membrane plane (z axis). The use of a fractional coordinate in place of a Cartesian one is required by the NPT ensemble: since the molecule can cross the box borders, a variable box length would otherwise cause discontinuities in the CV. The chosen CV (Figure 2) has values between 0 and 1 along the simulation box, with a value of 0.5 corresponding to the middle of the box, i.e. roughly in the middle of the lipid bilayer. The CV is periodic, spanning across the box limits. Value of the full width at half-maximum of the Gaussians was chosen to be 0.005, in CV units. Grid spacing was set to 0.001. Further details on the choice of the CV and on the convergence and analysis of the wt-metaD simulations will be provided, when required, in the Results and Discussion section of the paper.

RESULTS AND DISCUSSION Modeling the erythrocyte membrane As stated in the Introduction, focusing on the effect of membrane composition on the interactions with silica nanoparticles is of particular importance and, among the various cell types, red blood cells (erythrocytes) represent an obvious choice for a number of reasons, particularly their relevance to nano-toxicological studies.10,11,48 Yet, also limiting the study on an individual cell type, the complexity of the cell membrane poses enormous hurdles in its simulation. A priori, no factor determining the average membrane composition for a determined cell type (ratio between the different phospholipid classes, cholesterol content, glycosylation pattern and exposed proteins) can be ruled out as unimportant in the cell-silica interaction. However, due to the intrinsic limits of the chosen methodology, both in terms of size and time scale of the simulations, we decided here to focus exclusively on the lipid fraction of the erythrocyte cell membrane, leaving to further studies the investigation of the role played by the other components (proteins, in particular). Even so, the average lipid composition of an erythrocyte membrane includes at least 160 different types of lipid species,34 plus cholesterol,33 unevenly distributed among the inner and outer leaflets



of the bilayer.16 This is problematic for two reasons: firstly, including all the representative lipid species would require simulating a very large patch of cell membrane to reproduce the correct proportions, significantly increasing the computational cost of the simulations; secondly, more components included in the model mean a much harder interpretation of the results, i.e. obtaining the underlying structure-activity relationship. For this reason, we chose to follow a bottom-up approach, by gradually increasing the complexity of our cell membrane models, aiming at pinpointing the effect of specific components on the membrane-silica interaction. The building and features of these models are detailed in the Computational Details. The simplest model is a homogeneous bilayer containing a single lipid species, 1-palmitoyl-2oleoyl-sn-glycero-3-phosphocholine (POPC model, Figure 3a and Figure S2). Phosphatidylcholines (PC) represent the largest class of lipids in the erythrocyte lipid bilayer (~30% of the lipid fraction)33,35 and POPC in particular is one of the most abundant lipids in animal cell membranes.40,49 One advantage of dealing with mono-component homogeneous bilayers is the availability of experimental data to directly compare with, since this kind of membranes can be easily synthesized in laboratory.


Page 12 of 54

Page 13 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. (a) Homogeneous 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane model. Phosphatidylcholines in blue, ions in pink and water as a transparent surface. Simulation box in blue. (b-f) Free energy profiles as a function of the collective variable CV, i.e. the fractional coordinate of the center of mass of the molecule along the direction normal to the membrane, from well-tempered metadynamics (wtmetaD) simulations of water and the SIn clusters of Figure 1 across the POPC membrane model. The evolution of the free energy profiles during the wt-metaD simulations is reported according to the corresponding color scale. In all cases the zero of free energy corresponds to the minimum of the final profile.



Figure 4. (a) Mixed POPC + 20% cholesterol (POPC_CH) membrane model. Color code as in Figure 3. Cholesterol molecules in yellow. (b-f) Free energy profiles as a function of the collective variable from wtmetaD simulations of water and the SIn clusters of Figure 1 across the POPC_CH membrane model. Details as in Figure 3.

The subsequent step was to replace a fraction of the phospholipids with cholesterol, matching the average cholesterol content (~20%) of erythrocytes33 (POPC_CH model, Figure 4a and Figure S3). This was motivated by the important role that cholesterol is known to play in determining the chemical-physical properties of the cell membrane, in particular its flexibility and overall mechanical stability.50,51 Furthermore, recent coarse-grained based simulations have found that the structural reorganization of the membrane after uptake of anionic nanoparticles strongly depends on the presence and concentration of cholesterol.24 It must be stated here that this POPC_CH model, albeit informative, is not representative of a real system: indeed, it is experimentally known that


Page 14 of 54

Page 15 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


phase separation should be observed for a POPC/cholesterol binary system below 40% cholesterol concentration.52 The size and time scale of our all-atom simulations, however, are not large enough to show this behavior and the system remains perfectly stable during our production runs. We therefore believe that it should not impact the validity of our results.

Figure 5. (a) Heterogeneous red blood cell-like (RBC) membrane model. Color code as in Figure 4. Phosphatidylcholines in blue, phosphatidylserines in

red, phosphatidylethanolamines in green,

sphingomyelins in magenta. (b-f) Free energy profiles as a function of the collective variable from wt-metaD simulations of water and the SIn clusters of Figure 1 across the POPC_CH membrane model. Details as in Figure 3.

Finally, the third model (RBC model, Figure 5a and Figure S4) aims at representing the lipidic complexity of the erythrocyte membrane, by including, in the correct proportions, all phospholipids



Page 16 of 54

representing at least 2% in the typical composition of the bilayer.16,33–35 Figure 5a shows that lipid classes are realistically distributed asymmetrically between the leaflets. Of particular notice is the fact that phosphatidylserines, that expose a net negative charge at physiological pH, only occur on the inner leaflet of the membrane bilayer. This results in an asymmetric distribution of charges, leading to an electrostatic potential difference between the inner and outer surface of the membrane model. The second-rank deuterium lipid order parameter, SCD, is a metric of the order of the acyl chains of the phospholipids and describes the orientation of the lipid tails with respect to the bilayer normal. It is defined as:

= 〈3 − 1〉 .

Eq. 1

The brackets denote that it is averaged over all lipids of the system and over the time of the simulation. θ is the angle between the bilayer normal and the vector joining two Cn-1 and Cn+1 atoms of the acyl chain. A value of 1 corresponds to perfect alignment with the bilayer normal and 0 to a random orientation of molecules. Computational prediction of this parameter is particularly valuable since it can also be gathered experimentally by 2H NMR spectroscopy. The value of SCD along the palmitic acyl chains is reported in Figure 6 for both the POPC and POPC_CH models (full and empty red circles). Although in both cases the values describe a system in the so-called liquid disordered phase,34 the bilayer with the added cholesterol shows a higher order of the acyl chains: this underlines the known fact that the cholesterol rigid structure stiffens the carbons in its surroundings.24,50 As the tails of the lipids can move more freely, the order parameter decreases towards the end of the chain in both cases. SCD values for the RBC model (data not shown) matches those obtained for POPC_CH, confirming that the stiffening effect is indeed due to the added cholesterol.


Page 17 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Deuterium order parameters (SCD) of each atom of the palmitic acyl tail of POPC, for both the homogeneous POPC (full symbols) and the mixed POPC + 20% cholesterol (empty symbols) membrane models. In both cases, results are reported both in absence of silica (red circles) and when a silica particle is at the center of the lipid phase, in this case comparing lipids close to the particle (at less than 15Å from any silica atom, orange triangles) and lipids far from the particle (at more than 20Å from any silica atom, blue squares). Error in the SCD calculation is visualized as a shaded area around the curves.

Modeling the silica–membrane interaction With the aim of giving atomistic insights on silica biomineralization as well as nanotoxicity, a similar bottom-up approach was followed. Our interest here lies not in the formation mechanisms or the energetics of the condensation process, but only on the interactions between the condensates and the cell membrane. To this purpose we extracted possible nanocluster structures from an amorphous silica bulk, obtaining the SIn models shown in Figure 1. Most of the results have been obtained using uncharged models, with fully protonated silanol, as representative of the monomers and clusters in silification processes, that typically occur around pH 3 (surface silanols have an average pKa of about 7.5853).7 Surface charge can be relevant, however, for the hemolytic behavior of silica, that take place at neutral pH and the issue of particle charge will be addressed later in the paper.



By comparing results obtained with the silica models of Figure 1 we can: a) determine if and at which point of the silica nucleation process the condensate favors the interaction with the membrane surface with respect to the overlying water phase; b) study the effect of cluster size and structure on the lipid bilayer, particularly when particle internalization occurs. To obtain such information, knowledge of how the free energy of the system varies with respect to the cluster location (in the water phase, at the water-lipid interface, within the lipid bilayer) is required. From the obtained free-energy profile one can then extract evidence of a preferential segregation of the cluster at the membrane, i.e. if a free energy minimum is present at the interface, or quantify the free energy barrier for the cluster translocation across the membrane. Such free energy profile can only be obtained from a complete sampling of the configurational space of the system. Standard MD at a finite temperature T is not able to provide such sampling in a reasonable amount of computational time, due to the high energy barriers (that is, larger than a few kBT) that the system must cross.44 As detailed in the Computational Details, we resorted on well-tempered metadynamics (wt-metaD) simulations45 to achieve the required sampling, with the fractional coordinate of the center of mass of the cluster along the direction (z) normal to the membrane as the collective variable (CV). Metadynamics, with the correct set up, has been proven to be reliable for the computation of the transfer free energy of small solutes into model lipid membranes.54 Figure 2 shows the range of the CV superimposed to a sample free energy profile. For clarity, the mass densities of water and lipids are also reported. In each system we can define the CV region where water density is zero as the lipid phase. Being occupied by the acyl tails, it is fully hydrophobic. The CV region where water density corresponds to the TIP3P water density in liquid bulk (ρwat ≈ 0.998 g/cm3) is the water phase. Since the membrane is fluid and the position of the phospholipid polar heads changes continuously during the simulation time, the water-membrane interface region has a considerable thickness, and can be defined as the range of CV in which the two mass densities cross each other (see Fig. 2).


Page 18 of 54

Page 19 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Silica-membrane free energy profiles Figures 3-5 (b-f) report the obtained free energy profiles across the three POPC, POPC_CH and RBC membrane models for a water molecule, chosen as a reference (b), and for the four SIn silica nanoclusters (c-f). The graphs show the evolution of the profiles along the wt-metaD simulations, color coded according to the simulation time, displaying the gradual deposition of the bias potential in the CV space, slowly converging on the underlying free energy profile. Energies are provided, in kJ mol-1, relative to the minimum of the last profile. Simulations were run for at least 300 ns. However, many of them were run longer, since the profiles did not appear satisfyingly converged at that point. In the case of SI9b in interaction with the POPC_CH bilayer, for example, the final simulation time reached ~600 ns. It can be noticed that, apart from water, the region corresponding to the lipid phase is explored by the system mostly at later times of the simulation, since the free energy basin in the water phase needs to be filled to let the particle overcome the free energy barrier that gives access to the inside of the bilayer. Indeed, as expected for highly hydrophilic molecules like our silica nanoclusters, the region of minimal free energy lies completely within the water phase, where the profiles are almost perfectly flat. On the other hand, for the same reason, free energy reaches a maximum in the lipid phase, for all molecules and membranes. There appears to be, however, a marked difference between smaller molecules (water, SI1 and SI2), for which the free energy profile is rather flat in that region, and larger and more complex ones (SI6b, SI9b), that show for all membranes a “peaklike” shape of the free energy across the bilayer. One noticeable feature shared by all profiles is the clear absence of any minimum at the water-lipid interface. This means that no preferential segregation of the silica clusters from solution to the interface (or, in other words, no adsorption at the membrane surface) is expected under the heresimulated conditions. Minima at this interface have been observed in simulations of similar systems, both for the case of small molecules55 and of non-silica nanoparticles,24 albeit with great variety depending on the chemical nature of the solute. Indeed, recent experimental results on the



Page 20 of 54

interaction between silica nanoparticles and model membranes have already suggested that electrostatic interactions between phospholipids and silanol groups are not the driving force for membrane disruption.17

Energetics of the silica-membrane interaction So far, we have provided a qualitative assessment of the free energy profiles of Figures 3-5. From these profiles, however, a more quantitative description of the silica-membrane interaction can be obtained, by computing the free energy differences associated with the involved processes. We will focus here on three measures of the interaction: the crossing barrier, ∆G‡, the free energy of partition, ∆Gpart, and the free energy of adsorption, ∆Gads. The meaning of these energies is visually represented in Figure 2, while the corresponding values, averaged through the last 50 ns of each wtmetaD run, are reported in Table 1. The associated error is computed as the standard deviation during this time. The simplest approach is to compute the free energy difference between the minimum (in water) and the maximum (in the lipid phase) of free energy, evaluating the energy barrier, ∆G‡, that the silica cluster must cross in order to pass the lipid bilayer. This measure, however, is extremely sensitive to both the choice of the CV and the convergence behavior of the metaD simulation. One way to mitigate the variability in the quantitative analysis of the silica-membrane interaction is to focus on the whole free energy profile instead than looking at individual points. We can first compute the probabilities ρwat and ρlip of finding the cluster in the two water and lipid phases, respectively, by Boltzmann integration of the one-dimensional free energy profile G(CV) along the collective variable CV (see Fig. 2):

= +,

"$%&,'%( 1 , −, "$%&,')*

"-).,'%( 1 = +,, −+,, "-).,')*

! "

! "

#

#


Eq. 2

Eq. 3

Page 21 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The free energy of partition between the two phases can be computed as:

∆0,1 = −23 4 ln 7

+, 8

Eq. 4

with a positive value meaning that the molecule prefers the water phase to the lipid phase. We can use the same approach to estimate the free energies of particle adsorption at the membranewater interface, by focusing on the free energy of adsorption (∆Gads):

∆09: = −23 4 ln 7

8

Eq. 5

with

=

"-).,')*, 1 +,, −, "$%&,'%(

! "

#

Eq. 2

For the POPC and POPC_CH cases, the two faces of the bilayers were considered equivalent, while for the RBC membrane model the inner and outer surfaces were treated separately, as reported in Table 1.

∆G‡

∆Gpart

∆Gads

H2 O

SI1

SI2

SI6b

SI9b

POPC

28.6 ± 1.3

33.3 ± 7.4

31.5 ± 1.0

82.1 ± 0.7

85.2 ± 0.8

POPC_CH

27.5 ± 1.4

40.6 ± 1.5

57.7 ± 2.1

74.2 ± 1.5

124.6 ± 0.6

RBC

32.5 ± 1.6

32.5 ± 1.6

54.7 ± 0.8

77.1 ± 0.8

106.0 ± 21.5

POPC

19.0 ± 0.9

19.0 ± 4.0

18.6 ± 1.1

33.9 ± 1.3

42.8 ± 0.5

POPC_CH

24.6 ± 1.4

37.2 ± 1.7

26.8 ± 1.1

29.2 ± 0.8

39.7 ± 1.0

RBC

24.9 ± 0.9

27.1 ± 1.6

33.8 ± 1.5

41.3 ± 1.2

50.7 ± 8.1

POPC

-0.3 ± 0.3

-0.3 ± 0.6

0.6 ± 0.4

1.2 ± 0.8

2.1 ± 0.5

POPC_CH

0.9 ± 0.6

1.5 ± 0.5

1.3 ± 0.4

3.9 ± 0.9

-0.3 ± 0.4

RBC (i)

3.6 ± 0.4

3.8 ± 0.5

2.9 ± 0.8

3.4 ± 0.7

4.5 ± 0.9

RBC (o)

2.6 ± 0.5

1.9 ± 0.6

1.5 ± 0.4

1.2 ± 0.6

2.2 ± 0.1

Table 1. Energetics (in kJ mol-1) of the silica/membrane system. Results are reported for the SIn clusters of Figure 1 and for water as a reference, for all three membrane models (Figures 2-4). ∆G‡ is the membrane crossing barrier; ∆Gpart is the free energy of partition between the water and lipid phases; ∆Gads is the free energy of adsorption at the water-lipid interface. In the latter case, results are given separately for the inner



(i) and outer (o) layers of the RBC membrane model, due to their heterogeneous composition. Figure 2 visualizes the meaning of these energies on an exemplary free energy profile. All values have been obtained by averaging through the last 50 ns of each wt-metaD simulation and include an error computed as the standard deviation in this range.

Membrane crossing barrier (∆G‡) As a reference, the crossing barrier (∆G‡) for water was computed as 28.6 kJ mol-1 for the POPC membrane. It remains mostly unchanged for the POPC_CH case, with just a moderate increase for the crossing of the erythrocyte RBC model membrane. These values compare well with a simulation recently reported in literature (26.3 kJ mol-1) for water crossing a homogeneous POPC membrane, with a similar setup but employing a different free energy estimation method (Orthogonal Space Tempering).56 These data are also in good agreement with the experimental water to hexadecane transfer free energy (24.9 kJ mol-1, from the Minnesota Solvation Database). Looking at the barriers for the silica clusters, for the POPC membrane, SI1 and SI2 show similar values (~30 kJ mol-1) that are only slightly larger than water. On the other hand, SI6b and SI9b display much larger free energy barriers (~80 kJ mol-1), suggesting that the increase in size, as expected, further hinders the crossing of the membrane by the particle. Adding cholesterol to the membrane causes mixed results. While in some cases (SI1 and SI6b) the addition does not impact the crossing barrier, in other cases (SI2 and SI9b) they increase significantly (from 85.2 to 124.6 kJ mol-1 in the SI9b case). The result is that for the POPC_CH model a clear trend of ∆G‡ with size can be observed (H2O < SI1 < SI2 < SI6b < SI9b). The same trend with size is obtained for the RBC model, suggesting a similar impact of cholesterol on the crossing barrier. However, ∆G‡ values for the silica clusters are generally lower than the POPC_CH case, in some instances being even lower than the pure POPC bilayer. One tentative explanation is the different proportion of saturated vs unsaturated lipid tails (cf. RBC model composition in Table S2) that might affect the lipid mobility upon particle insertion, i.e. the crossing could be less hindered in the RBC case, despite the cholesterol stiffening effect. Some observed non-monotonous trends in the results suggest,


Page 22 of 54

Page 23 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


however, that this measure of the interaction might be overly affected by technical convergence issues in the free energy profiles. Indeed, as can be noticed from Figures 3-5, the maximum of free energy, in the middle of the membrane, often shows abrupt changes also at later stages of the wtmetaD simulations, due to the pathologically poor sampling of this region. This is particularly evident for larger particles, for which most of the sampling of the lipid phase occurs only after several hundreds of ns. In turn, this behavior of the profiles causes rapid changes in the computed ∆G‡ (Figure S5 in the Supporting Information) during the simulation time. Free energy of partition (∆Gpart) Since we are now averaging along a larger region of the CV space, ∆Gpart shows a smoother convergence with time than ∆G‡ (Figure S6 in the Supporting Information). Looking at the ∆Gpart values in Table 1, they are in all cases, as can be expected, lower than the corresponding ∆G‡ values, particularly for the larger particles, whose “peak-shaped” free energy profiles make more evident the difference between the two measures of the silica-lipids energetics. The ∆Gpart values are all positive, reflecting the hydrophilic character of the investigated molecules. Indeed, starting from one frame with SI2 inside the lipid phase extracted from the wt-metaD trajectory, the molecule was expelled from the membrane in a few ps of unconstrained MD, returning and then remaining in the water phase. Interestingly, differences between the membrane models are also smoothed out by this approach. This was anticipated, since the partition free energy captures mainly the difference in solubility between the two phases and this is not supposed to be largely affected by the membrane composition. Nevertheless, the presence of cholesterol, both in the POPC_CH and RBC cases, produces generally larger ∆Gpart values for most of the molecules, suggesting that cholesterol further decreases the solubility of the particles in the lipid phase, since its stiffening effect on the acyl tails reduces the membrane ability to quickly adapt to the internalized molecule. Differences among particle sizes seem also reduced, with, in some cases (c.f. SI1 vs SI6b in the POPC_CH case) the silicic acid monomer exhibiting similar or even larger ∆Gpart values than bigger clusters. As will be described afterwards, internalization of larger clusters involves the formation of water-filled holes



in the membrane, with the molecules keeping contact with water also when in the middle of the bilayer. This does not happen for SI1 and sometimes SI2, that are usually completely surrounded by the hydrophobic lipid tails at CV = 0.5. As already evidenced in simulations on aminoacids across lipid bilayers,57 the balance between the energetics of maintaining hydration and deforming the membrane is delicate and can strongly effect the partition free energy. Apart from the silicic acid molecules, the largest ∆Gpart values were found for the RBC membrane, with in this case a clear trend with particle size. The presence of lipids with unsaturated acyl tails (less mobile) in this model possibly enhances the stiffening effect of cholesterol. Free energy of adsorption (∆Gads) Finally, the evolution of the ∆Gads values during the simulation time (Figure S7 in the Supporting Information) shows a higher degree of convergence, due to the good sampling of the considered CV regions. The ∆Gads values in Table 1 quantify what was already evident from an inspection of the free energy profiles of Figures 3-5: no preferential interaction between the silica clusters and the membrane surface can be detected. The obtained ∆Gads values are either insignificantly small or slightly positive, in this case suggesting a marginal preference by the silica particles for the interaction with the bulk water phase. Small differences can be observed, for the RBC case, when separating the values for the inner and outer layers of the membrane, since the heterogeneous composition makes the two surfaces inequivalent: adsorption at the outer layer appears less unfavorable for all silica models. This might be due to the presence, on the inner leaflet, of head groups exposing a net negative charge, although the small absolute numbers for ∆Gads do not allow this observation to be conclusive.

Details of the silica-membrane interaction By extracting and analyzing individual trajectory frames from the wt-metaD simulations, corresponding to different configurations of the systems (e.g. particle at the interface vs particle


Page 24 of 54

Page 25 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


inside the lipid phase), precious information can be gained regarding the atomistic details of the silica-membrane interaction, that might help putting the energy values in Table 1 in more context.

Figure 7. (a) Distribution of the number of H-bonds between the silanols (hydroxyl) groups of the SI6b silica cluster and the phosphate groups of the POPC phospholipids in the homogeneous POPC model. Data are computed tacking in to account the frames of the wt-metaD simulations in which the silica particle is in the middle of the water-lipid interface. (b) Comparison between the distributions of the total number of Hbonds made by the SI6b silica cluster when at the interface with the POPC model (red) and when in the water phase (blue). Lines are Gaussian fittings to help the comparison.

Silica at the membrane/water interface As already stated, no driving force for preferential segregation of the silica clusters at the membrane/water interface has been observed. All particles, decidedly hydrophilic due to the high number of silanol (OH) groups that they expose (that can act as both strong donors and acceptors of hydrogen bonds19,58), interact stably with the water molecules when they are in the water phase. However, the computed ∆Gads values (Table 1) indicate that only a negligibly small energy penalty is associated with the clusters moving from the bulk water region to the interface region. This is because the polar functionalities exposed by the phospholipids heads, particularly the phosphate groups, can act as new sites of hydrogen-bond interaction, as shown in Figure 7.a for the exemplary case of SI6b interacting with the POPC membrane. On average, the total number of hydrogen bonds



formed by the silica clusters is slightly smaller when they are at the interface compared to when they are in the water phase (Figure 7.b for the SI6b@POPC case), which may explain the very small positive adsorption energies reported in Table 1. Similar analyses for the particles at the interface with the POPC_CH membrane model show that cholesterol participates in the interactions with silica through its polar hydroxyl functionality, that it exposes at the membrane surface (data not shown). As regard the RBC case, one question is whether the lipid composition affects the interaction. We found that the distribution of the silica hydrogen bonds among the lipid classes represented in the model follows accurately the relative ratios of the lipids within the bilayer, suggesting that no preferential interaction exists with any of the lipid species. Silica inside the membrane Of great interest, particularly regarding possible pathways leading to silica-induced nanotoxicity, is what happens when the silica nanoclusters cross the water/membrane interface, move into the lipid bilayer and eventually emerge on the opposite side of the membrane. Indeed, if signs of a destabilization of the membrane structure can be spotted in our simulations, these could in turn provide some atomistic rationalization to the red blood cell membrane rupture (hemolysis) observed in nanotoxicity tests. For this purpose, we extracted from the wt-metaD simulations all the frames in which the silica cluster lies inside the lipid phase. From a visual inspection (as shown by the snapshots reported in Figure 8.a, for the POPC case), one noticeable common feature is the formation of “holes” filled with water molecules that connect the cluster within the membrane to the external solvent region. These holes have increasing diameter following the size of the silica clusters. In all cases, the simulations show that their formation is preceded by an interaction of the silica cluster with the lipid heads. The fluidity of the hydrocarbon tails then allows the particle to retain a large amount of its hydration shell while entering the lipid membrane region. Reported here for the first time for silica clusters, similar water-filled membrane defects have been indeed


Page 26 of 54

Page 27 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


described in literature in the case of membrane permeation by several polar and charged species, such as ions and aminoacids.57,59

Figure 8. (a) Snapshots displaying silica clusters (depicted as vdW spheres) of increasing size inside the lipid bilayer. To show the formation of “holes” in the membrane, the lipids are displayed as a transparent isosurface and water molecules (ball and stick models) close to the lipids are also included. (b) Radial pair distribution function (ga,b) between oxygens of silica and water averaged along all frames in which the silica particle (SI6b, in this example) is inside the POPC lipid bilayer. (c) Magnification of the SI6b case of Figure 8.a, to better show the co-translocation of water molecules inside the membrane.

The events leading to eventual membrane crossing are closely related with the capacity of the solvent to re-catch the particles on the opposite side of the membrane. This is due to intrinsic membrane/solvent fluctuations combined with the relatively low energy barriers associated with the penetration of water molecules deep into the lipid region (see Figs. 3-5). Only in the case of the smallest clusters, SI1 and SI2, does the hole close over the molecule causing its full internalization before it is expelled again toward the exterior, after renewed water-particle contact (Figure 9). This



behavior is in line with the rather flat free-energy peaks in the profiles associated with the crossing (Figs. 3-5). For the larger clusters, SI6b and SI9b, the hole remains fully open even when the particle is in the very middle of the lipid bilayer (Figure 8.a). The transition to the opposite side is the result of a sharp instability (pointed free-energy peaks), in which the particle is wet by partial solvent holes either on one or the other side of the membrane.

Figure 9. Snapshots displaying the crossing of the POPC lipid bilayer by the SI1 silica cluster. The water and lipid phases are visualized as transparent isosurfaces (cyan and white, respectively), while the silica cluster and water molecules close (< 5Å) to it are reported in ball and stick models. The corresponding value of the CV is also provided.

The extent of water co-translocation can be observed by looking at the O(sil)–O(wat) radial pair distribution function, computed when the particle is in the lipid phase (Figure 8.b). The clear peaks characteristics of hydrogen bond interactions between the two components are present, confirming that water is internalized together with the particle. It must be noted, however, that in no case these observed hole features evolved in full pores spanning across the bilayer. Whether this may eventually happen for larger nanoparticles or for nanoparticle aggregates characteristic of fumed silica nanopowder is an interesting question that


Page 28 of 54

Page 29 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


shall be answered in follow-up works. The formation of transmembrane water-filled pores would in fact have important implications for silica-induced hemolysis. The here-observed temporary silica internalization within the lipid bilayer appears to have impact on the structure of the whole simulated membrane patch. One way to visualize this is by measuring the membrane deformation during one internalization event. Deformation is defined as the interpolated location of the head-groups after translating the chosen leaflet's center at z = 0.60 Therefore, deformation takes on both negative and positive values. Irrespective of the leaflet, deformation is positive when it expands the membrane, and negative otherwise. By plotting the time-averaged deformation on the xy plane, for both leaflets, a deformation map can be generated, as the one reported in Figure 10 for one penetration of the POPC membrane by the SI9b particle. Although the maximum size of this particle is about 1 nm when fully extended, its effect on the bilayer is not restricted to its local environment. A deformation is measured on the whole plane of the membrane, in the form of an inflection around the newly formed hole (green area in Figure 10, left), with the neighboring parts elevating in response. Interestingly, the same deformation map appears in reverse on the opposite leaflet (Figure 10, right) with a bulge in correspondence of the particle, showing that the particle-induced deformation spans through the entire model system, involving both lipid layers. The effects of particle penetration can be also measured on a more local scale. As already stated, the deuterium lipid order parameter, SCD, is a measure of the mobility of the lipid acyl chains, reflecting the order of the bilayer, and, as reported in Figure 6, this can be used to quantify the stiffening effect of cholesterol on the membrane. Coarse-grained MD simulations of anionic nanoparticles across lipid bilayers have recently shown a local disordering of the lipids induced by the particle.24 We investigated this effect in our simulations. The results are reported, for the SI9b case inside the POPC and POPC_CH bilayers, in Figure 6. When limiting the measure of the order parameter to the lipids in the vicinity of the particle (closer than 15 Å), a decrease with respect to the pristine membranes can be noticed. This indicates a local disordering induced by the nanocluster, more



pronounced, as expected, on the middle section of the acyl chains. When we move further away from the particle, we recover the original SCD values in the POPC case, while still falling slight below the originals for the POPC_CH model, suggesting that the particle-induced disordering might involve a larger area when cholesterol is present in the membrane.

Figure 10. Time-averaged deformation maps for the top and bottom leaflets, reported for one event in which the SI9b cluster is inside the POPC membrane. Deformation is defined as the interpolated location of the head-groups after translating the chosen leaflet's center at z = 0. Therefore, deformation takes on both negative and positive values. Irrespective of the leaflet, deformation is positive when it expands the membrane, and negative otherwise. The particle location corresponds to the area of maximum negative deformation of the top leaflet, that is the side from which it entered the membrane.

Effect of particle charge The hydroxyl groups (silanols) exposed by a silica particle are moderately acidic, with an average pKa of about 7.53 However, assignment of a pKa value to specific silanol is particularly difficult and the number of negatively charged sites on a silica surface is strongly dependent on the pH of the solution and the ionic strength.19 This results in very different measured values for the surface charge density of silica particles, with respect to different synthesis methods (sol-gel vs. fumed silica, for example) and to heat treatment.61,62 Furthermore, no data exist regarding the expected charge density for silica nanoclusters of the size scale investigated in this work. In addition, aiming


Page 30 of 54

Page 31 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


at investigating both biomineralization and nanotoxicity, the large range of pHs involved in these processes further complicates the issue. For example, the artificial biomineralization of cells performed in Kaehr et al. was performed at pH 3, where the silica clusters are generally uncharged,8 that is very far from the usual pH of blood, that lies around 7.63 For this reason, the choice was made to perform most of the simulations without simulating the deprotonation of silica silanols, as reported in the structures of Figure 1. Surface charge is known, however, to have a big impact on the adsorption behavior of silica.21 For comparison, one wt-metaD simulation was performed for the POPC case with a negatively charged SI6b2- cluster, in which two silanols were deprotonated (Figure S8 in the Supporting Information). The simulation failed to converge after 500 ns, lacking a sampling of the inner part of the membrane. This implies that the energy barrier to cross the watermembrane interface (that cannot be estimated due to the under-sampling) is much higher than the neutral particle. This is could be partly explained with repulsive interactions between the negatively charged silanols on the silica cluster and the negatively charged phosphate groups on the lipid heads (particularly when the membrane exhibits a net negative charge, such as in the RBC case, not tested here), but it can also be related to the energy cost of inserting a stably hydrated charged group into the hydrophobic membrane core.57 Moreover, also in this case the free-energy profile did not show any local minima in correspondence of the membrane/water interface, ruling out spontaneous interface segregation of negatively charged clusters.

CONCLUSIONS We have here presented a first investigation, by means of all-atom molecular dynamics simulations augmented with enhanced sampling techniques, of the interaction of silica with the lipid phase of cell membranes. By using simple silica nanoclusters of different size and topology, we were able to pursue both a deeper understanding of the molecular recognition processes at the basis of biomineralization and give a hint for a possible pathway of nanotoxicity induced by small particles of colloidal synthesis. The lipid phase was simulated targeting a realistic representation of the lipid



composition of a red blood cell membrane, given its relevance to nano-toxicology, with particular interest in the effect of cholesterol. Well-tempered metadynamics allowed to overcome the time-scale limits of standard all-atom MD and to sample the full configurational space of the silica nanoclusters (in both the water and lipid phases and at the interface), providing free energy profiles of the molecules across the lipid bilayer. These profiles allowed us, in turn, to predict the energetics of the silica-membrane interaction, in the form of crossing barrier, ∆G‡, free energy of partition, ∆Gpart, between the water and lipid phases, and free energy of adsorption, ∆Gads, at the water-lipid interface, also investigating their dependence on particle size, particle charge and membrane composition. These results showed that silica nanoclusters are highly hydrated hydrophilic molecules that must overcome high barriers to cross the water-lipid interface already at nuclearities of a few Si atoms. However, the barriers encountered by individual water or Si(OH)4 molecules are of nearly equal height, strongly suggesting that either can cross the RBC membrane with roughly the same probability. This explains why intracellular structures are easily infiltrated by silicic acid and undergo silicification just as extracellular components do.7,8 Moreover, the energy barriers are markedly peaked in the central zone of the lipid region, but present smooth and shallow tails extending towards the membrane/water interfaces. This means that the clusters are able to spontaneously enter deep into the region occupied by the lipids heads, as in fact confirmed by our structural analyses. There, they form a tight network of hydrogen bonds with the negatively charged phosphate groups. Especially under acidic conditions, these may be thus actively involved in polycondensation reactions, leading their incorporation in a growing silica network, similarly as phosphate groups in phosphosilicate glasses.64 This putative process can explain the repeatedly observed silicification of phospholipid membrane structures even in the absence of a net driving force for segregation of small silica clusters at the membrane/water interface. Alternatively, segregation may take place at higher Si nuclearities that the ones


Page 32 of 54

Page 33 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


considered in this work; or membrane components other than lipids (e.g. proteins) may act as preferential adsorption sites for silicic acid moieties. Cholesterol was found to be an import actor in governing the interaction, through its stiffening action on the lipid tails, opposing the penetration of the membrane by the silica particles. When the particles reach the inner part of the lipid bilayers, however, they cause both local and global effects on the membrane structure. The former appears in the form of an induced disordering on the lipid tails within 1-2 nm from the particles; the latter in the form of water-filled holes deforming the entire membrane sheet. These conclusions, if confirmed for larger nanoparticles, could be in principle correlated to possible pathways of silica-induced nanotoxicity (hemolysis) based on local perforation and thus uncontrolled permeabilization of the cell membrane. From a methodological point of view, our simulations demonstrate that an all-atom approach is not only feasible, but also strictly necessary to study the interaction of small nanoparticles with the lipid membrane, an investigation that has usually been performed through less detailed coarse-grained simulations.22–24 Limitations remain, however, particularly in the form of the relatively short achievable time scales, that hinders the convergence of the free energy simulations and might not allow to observe longer-lasting effects of the particles on the lipid bilayer. Further work is needed to extend these simulations to larger nanoparticles models with varied morphologies, taking into account in particular both Stöber-like and fumed-like silica. Finally, other components of the cell membrane should be included in the simulations, with particular reference to membrane proteins, that may play an important role in nanoparticle/cell interactions, especially when nanoparticle incorporation occurs via endocytosis rather than passive translocation.

ACKNOWLEDGEMENTS Models have been visualized and manipulated by VMD.65 This work was supported by the Deutsche Forschungsgemeinschaft under grant CO 1043/11-1. Computational resources were provided by the North German Supercomputing Alliance (HLRN), project hbc00020. CJB was



supported by the Sandia National Laboratories LDRD program. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

ASSOCIATED CONTENT Supporting Information The Electronic Supporting Information contains: force field parameters for the silica nanoclusters and validation; composition of the RBC membrane model; images of the three membrane models; area per lipids and membrane thickness for the three membrane model; evolution of the free energy values of Table 1 along the wt-metaD simulations; free energy profile as a function of the collective variable from wt-metaD simulations of water and a charged silica nanocluster.

BIBLIOGRAPHY (1)

Iler, R. K. The Chemistry of Silica; 2nd ed.; John Wiley & Sons, Inc.: New York, Chichester, Brisbane, Toronto, 1979.

(2)

Hildebrand, M. Diatoms, Biomineralization Processes, and Genomics. Chem. Rev. 2008, 108, 4855–4874.

(3)

Ning, R. Y. Reactive Silica in Natural Waters — A Review. Desalin. Water Treat. 2010, 21, 79–86.

(4)

Cuko, A.; Macià, A.; Calatayud, M.; Bromley, S. T. Global Optimisation of Hydroxylated


Page 34 of 54

Page 35 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Silica Clusters: A Cascade Monte Carlo Basin Hopping Approach. Comput. Theor. Chem. 2017, 1102, 38–43. (5)

Trinh, T. T.; Jansen, A. P. J.; van Santen, R. A. Mechanism of Oligomerization Reactions of Silica. J. Phys. Chem. B 2006, 110, 23099–23106.

(6)

Bhattacharya, S.; Kieffer, J. Molecular Dynamics Simulation Study of Growth Regimes during Polycondensation of Silicic Acid: From Silica Nanoparticles to Porous Gels. J. Phys. Chem. C 2008, 112, 1764–1771.

(7)

Townson, J. L.; Lin, Y.-S.; Chou, S. S.; Awad, Y. H.; Coker, E. N.; Brinker, C. J.; Kaehr, B. Synthetic

Fossilization

of

Soft

Biological

Tissues

and

Their

Shape-Preserving

Transformation into Silica or Electron-Conductive Replicas. Nat. Commun. 2014, 5, 5665. (8)

Kaehr, B.; Townson, J. L.; Kalinich, R. M.; Awad, Y. H.; Swartzentruber, B. S.; Dunphy, D. R.; Brinker, C. J. Cellular Complexity Captured in Durable Silica Biocomposites. Proc. Natl. Acad. Sci. 2012, 109, 17336–17341.

(9)

Pavan, C.; Fubini, B. Unveiling the Variability of “Quartz Hazard” in Light of Recent Toxicological Findings. Chem. Res. Toxicol. 2017, 30, 469–485.

(10)

Pavan, C.; Tomatis, M.; Ghiazza, M.; Rabolli, V.; Bolis, V.; Lison, D.; Fubini, B. In Search of the Chemical Basis of the Hemolytic Potential of Silicas. Chem. Res. Toxicol. 2013, 26, 1188–1198.

(11)

Shi, J.; Hedberg, Y.; Lundin, M.; Odnevall Wallinder, I.; Karlsson, H. L.; Möller, L. Hemolytic Properties of Synthetic Nano- and Porous Silica Particles: The Effect of Surface Properties and the Protection by the Plasma Corona. Acta Biomater. 2012, 8, 3478–3490.

(12)

Zhang, H.; Dunphy, D. R.; Jiang, X.; Meng, H.; Sun, B.; Tarn, D.; Xue, M.; Wang, X.; Lin, S.; Ji, Z.; et al. Processing Pathway Dependence of Amorphous Silica Nanoparticle Toxicity: Colloidal vs Pyrolytic. J. Am. Chem. Soc. 2012, 134, 15790–15804.



(13)

Barbé, C.; Bartlett, J.; Kong, L.; Finnie, K.; Lin, H. Q.; Larkin, M.; Calleja, S.; Bush, A.; Calleja, G. Silica Particles: A Novel Drug-Delivery System. Adv. Mat. 2004, 16, 1959–1966.

(14)

Stöber, W.; Fink, A.; Bohn, E. Controlled Growth of Monodisperse Silica Spheres in the Micron Size Range. J. Colloid Interface Sci. 1968, 26, 62–69.

(15)

McIntosh, G. J. Theoretical Investigations into the Nucleation of Silica Growth in Basic Solution Part I - Ab Initio Studies of the Formation of Trimers and Tetramers. Phys. Chem. Chem. Phys. 2013, 15, 3155–3172.

(16)

Virtanen, J. A.; Cheng, K. H.; Somerharju, P. Phospholipid Composition of the Mammalian Red Cell Membrane Can Be Rationalized by a Superlattice Model. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 4964–4969.

(17)

Kettiger, H.; Québatte, G.; Perrone, B.; Huwyler, J. Interactions between Silica Nanoparticles and Phospholipid Membranes. Biochim. Biophys. Acta - Biomembr. 2016, 1858, 2163–2170.

(18)

Derr, L.; Hildebrand, N.; Köppen, S.; Kunze, S.; Treccani, L.; Dringen, R.; Rezwan, K.; Colombi Ciacchi, L. Physisorption of α-Chymotrypsin on SiO2 and TiO2: A Comparative Study via Experiments and Molecular Dynamics Simulations. Biointerphases 2016, 11, 011007.

(19)

Rimola, A.; Costa, D.; Sodupe, M.; Lambert, J.-F.; Ugliengo, P. Silica Surface Features and Their Role in the Adsorption of Biomolecules: Computational Modeling and Experiments. Chem. Rev. 2013, 113, 4216–4313.

(20)

Meißner, R. H.; Wei, G.; Ciacchi, L. C. Estimation of the Free Energy of Adsorption of a Polypeptide on Amorphous SiO2 from Molecular Dynamics Simulations and Force Spectroscopy Experiments. Soft Matter 2015, 11, 6254–6265.

(21)

Hildebrand, N.; Köppen, S.; Derr, L.; Li, K.; Koleini, M.; Rezwan, K.; Colombi Ciacchi, L. Adsorption Orientation and Binding Motifs of Lysozyme and Chymotrypsin on Amorphous


Page 36 of 54

Page 37 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Silica. J. Phys. Chem. C 2015, 119, 7295–7307. (22)

Salassi, S.; Simonelli, F.; Bochicchio, D.; Ferrando, R.; Rossi, G. Au Nanoparticles in Lipid Bilayers: A Comparison between Atomistic and Coarse-Grained Models. J. Phys. Chem. C 2017, 121, 10927–10935.

(23)

Rossi, G.; Monticelli, L. Gold Nanoparticles in Model Biological Membranes: A Computational Perspective. Biochim. Biophys. Acta - Biomembr. 2016, 1858, 2380–2389.

(24)

Gkeka, P.; Angelikopoulos, P.; Sarkisov, L.; Cournia, Z. Membrane Partitioning of Anionic, Ligand-Coated Nanoparticles Is Accompanied by Ligand Snorkeling, Local Disordering, and Cholesterol Depletion. PLOS Comput. Biol. 2014, 10, e1003917.

(25)

Heikkilä, E.; Martinez-Seara, H.; Gurtovenko, A. A.; Vattulainen, I.; Akola, J. Atomistic Simulations of Anionic Au144(SR)60 Nanoparticles Interacting with Asymmetric Model Lipid Membranes. Biochim. Biophys. Acta - Biomembr., 2014, 1838, 2852–2860.

(26)

Kim, J. M.; Chang, S. M.; Kong, S. M.; Kim, K.-S.; Kim, J.; Kim, W.-S. Control of Hydroxyl Group Content in Silica Particle Synthesized by the Sol-Precipitation Process. Ceram. Int. 2009, 35, 1015–1019.

(27)

Butenuth, A.; Moras, G.; Schneider, J.; Koleini, M.; Köppen, S.; Meißner, R.; Wright, L. B.; Walsh, T. R.; Ciacchi, L. C. Ab Initio Derived Force-Field Parameters for Molecular Dynamics Simulations of Deprotonated Amorphous-SiO2/Water Interfaces. Phys. status solidi 2012, 249, 292–305.

(28)

Meißner, R. H.; Schneider, J.; Schiffels, P.; Colombi Ciacchi, L. Computational Prediction of Circular Dichroism Spectra and Quantification of Helicity Loss upon Peptide Adsorption on Silica. Langmuir 2014, 30, 3487–3494.

(29)

Demiralp, E.; Çağin, T.; Goddard, W. A. Morse Stretch Potential Charge Equilibrium Force Field for Ceramics: Application to the Quartz-Stishovite Phase Transition and to Silica Glass.



Phys. Rev. Lett. 1999, 82, 1708–1711. (30)

Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A Web-Based Graphical User Interface for CHARMM. J. Comput. Chem. 2008, 29, 1859–1865.

(31)

Wu, E. L.; Cheng, X.; Jo, S.; Rui, H.; Song, K. C.; Dávila-Contreras, E. M.; Qi, Y.; Lee, J.; Monje-Galvan, V.; Venable, R. M.; et al. CHARMM-GUI Membrane Builder toward Realistic Biological Membrane Simulations. J. Comput. Chem. 2014, 35, 1997–2004.

(32)

Jo, S.; Lim, J. B.; Klauda, J. B.; Im, W. CHARMM-GUI Membrane Builder for Mixed Bilayers and Its Application to Yeast Membranes. Biophys. J. 2009, 97, 50–58.

(33)

Dougherty, R. M.; Galli, C.; Ferro-Luzzi, A.; Iacono, J. M. Lipid and Phospholipid Fatty Acid Composition of Plasma, Red Blood Cells, and Platelets and How They Are Affected by Dietary Lipids: A Study of Normal Subjects from Italy, Finland, and the USA. Am. J. Clin. Nutr. 1987, 45, 443–455.

(34)

Myher, J. J.; Kuksis, A.; Pind, S. Molecular Species of Glycerophospholipids and Sphingomyelins of Human Erythrocytes: Improved Method of Analysis. Lipids 1989, 24, 396–407.

(35)

Dodge, J. T.; Phillips, G. B. Composition of Phospholipids and of Phospholipid Fatty Acids and Aldehydes in Human Red Cells. J. Lipid Res. 1967, 8, 667–675.

(36)

Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; MondragonRamirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of the CHARMM AllAtom Additive Force Field for Lipids: Validation on Six Lipid Types. J. Phys. Chem. B 2010, 114, 7830–7843.

(37)

Paloncýová, M.; Fabre, G.; DeVane, R. H.; Trouillas, P.; Berka, K.; Otyepka, M. Benchmarking of Force Fields for Molecule–Membrane Interactions. J. Chem. Theory Comput. 2014, 10, 4143–4151.


Page 38 of 54

Page 39 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(38)

MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586–3616.

(39)

Kučerka, N.; Nieh, M.-P.; Katsaras, J. Fluid Phase Lipid Areas and Bilayer Thicknesses of Commonly Used Phosphatidylcholines as a Function of Temperature. Biochim. Biophys. Acta - Biomembr. 2011, 1808, 2761–2771.

(40)

Kučerka, N.; Tristram-Nagle, S.; Nagle, J. F. Structure of Fully Hydrated Fluid Phase Lipid Bilayers with Monounsaturated Chains. J. Membr. Biol. 2006, 208, 193–202.

(41)

McCaughan, L.; Krimm, S. X-Ray and Neutron Scattering Density Profiles of the Intact Human Red Blood Cell Membrane. Science 1980, 207, 1481 LP-1483.

(42)

Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1–2, 19–25.

(43)

Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463–1472.

(44)

Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing Important Fluctuations: Rare Events and Metadynamics from a Conceptual Viewpoint. Annu. Rev. Phys. Chem. 2016, 67, 159–184.

(45)

Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 20603.

(46)

Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; et al. PLUMED: A Portable Plugin for FreeEnergy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961– 1972.



(47)

Page 40 of 54

Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. 2002, 99, 12562–12566.

(48)

Dobrovolskaia, M. A.; Clogston, J. D.; Neun, B. W.; Hall, J. B.; Patri, A. K.; McNeil, S. E. Method for Analysis of Nanoparticle Hemolytic Properties In Vitro. Nano Lett. 2008, 8, 2180–2187.

(49)

Tattrie, N. H.; Bennett, J. R.; Cyr, R. Maximum and Minimum Values for Lecithin Classes from Various Biological Sources. Can. J. Biochem. 1968, 46, 819–824.

(50)

Fantini, J.; Garmy, N.; Mahfoud, R.; Yahi, N. Lipid Rafts: Structure, Function and Role in HIV, Alzheimer’s and Prion Diseases. Expert Rev. Mol. Med. 2002, 4, 1–22.

(51)

Moon, S.; Yan, R.; Kenny, S. J.; Shyu, Y.; Xiang, L.; Li, W.; Xu, K. Spectrally Resolved, Functional Super-Resolution Microscopy Reveals Nanoscale Compositional Heterogeneity in Live-Cell Membranes. J. Am. Chem. Soc. 2017, 139, 10944–10947.

(52)

de

Almeida,

R.

F.

M.;

Sphingomyelin/Phosphatidylcholine/Cholesterol

Fedorov, Phase

A.; Diagram:

Prieto,

M.

Boundaries

and

Composition of Lipid Rafts. Biophys. J. 2003, 85, 2406–2416. (53)

Sauer, J.; Ugliengo, P.; Garrone, E.; Saunders, V. R. Theoretical-Study of Van-Der-Waals Complexes at Surface Sites in Comparison with the Experiment. Chem. Rev. 1994, 94, 2095– 2160.

(54)

Bochicchio, D.; Panizon, E.; Ferrando, R.; Monticelli, L.; Rossi, G. Calculating the Free Energy of Transfer of Small Solutes into a Model Lipid Membrane: Comparison between Metadynamics and Umbrella Sampling. J. Chem. Phys. 2015, 143, 144108.

(55)

Menichetti, R.; Kanekal, K. H.; Kremer, K.; Bereau, T. In Silico Screening of DrugMembrane Thermodynamics Reveals Linear Relations between Bulk Partitioning and the Potential of Mean Force. J. Chem. Phys. 2017, 147, 125101.


Page 41 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(56)

Lv, C.; Aitchison, E. W.; Wu, D.; Zheng, L.; Cheng, X.; Yang, W. Comparative Exploration of Hydrogen Sulfide and Water Transmembrane Free Energy Surfaces via Orthogonal Space Tempering Free Energy Sampling. J. Comput. Chem. 2016, 37, 567–574.

(57)

MacCallum, J. L.; Bennett, W. F. D.; Tieleman, D. P. Distribution of Amino Acids in a Lipid Bilayer from Computer Simulations. Biophys. J. 2008, 94, 3393–3404.

(58)

Delle Piane, M.; Corno, M.; Ugliengo, P. Ab Initio Modeling of Hydrogen Bond Interaction at Silica Surfaces With Focus on Silica/Drugs Systems. In Modelling and Simulation in the Science of Micro- and Meso-Porous Materials 1st Edition; Catlow, C. R. A.; van Speybroeck, V.; van Santen, R., Eds.; Elsevier, Amsterdam, 2018; pp. 297–328.

(59)

Vorobyov, I.; Olson, T. E.; Kim, J. H.; Koeppe II, R. E.; Andersen, O. S.; Allen, T. W. IonInduced Defect Permeation of Lipid Membranes. Biophys. J. 2014, 106, 586–597.

(60)

Guixà-González, R.; Rodriguez-Espigares, I.; Ramírez-Anguita, J. M.; Carrió-Gaspar, P.; Martinez-Seara, H.; Giorgino, T.; Selent, J. MEMBPLUGIN: Studying Membrane Complexity in VMD. Bioinformatics 2014, 30, 1478–1480.

(61)

Kobayashi, M.; Skarba, M.; Galletto, P.; Cakara, D.; Borkovec, M. Effects of Heat Treatment on the Aggregation and Charging of Stöber-Type Silica. J. Colloid Interface Sci. 2005, 292, 139–147.

(62)

Sonnefeld, J.; Göbel, A.; Vogelsberger, W. Surface Charge Density on Spherical Silica Particles in Aqueous Alkali Chloride Solutions. Colloid Polym. Sci. 1995, 273, 926–931.

(63)

Kellum, J. A. Determinants of Blood PH in Health and Disease. Crit. Care 2000, 4, 6–14.

(64)

Plotnichenko, V. G.; Sokolov, V. O.; Koltashev, V. V; Dianov, E. M. On the Structure of Phosphosilicate Glasses. J. Non. Cryst. Solids 2002, 306, 209–226.

(65)

Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph.



1996, 14, 33–38.


Page 42 of 54

Page 43 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC GRAPHIC



Figure 1. Model silica particles used in the present work. SIn models represent small silica clusters containing n SiO4 units, arranged in a linear or branched (b) fashion. 165x127mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 44 of 54

Page 45 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Schematics of how the different free energy differences of Table 1 (crossing barrier, ∆G‡, free energy of partition, ∆Gpart, free energy of adsorption, ∆Gads) have been computed. An exemplary free energy profile (from Figure 3c) is reported (red line), superimposed to the mass densities of both water (blue line) and lipids (black line) averaged along the corresponding simulation. Both free energy and density are plotted as a function of the collective variable CV, i.e. the fractional coordinate along the direction (z) normal to the membrane. The background is colored according to the different regions: water phase (cyan), where the water density is maximal, lipid phase (grey) where the water density is zero, interface (yellow) where the water and lipid densities crosses. 356x292mm (300 x 300 DPI)



Figure 3. (a) Homogeneous 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane model. Phosphatidylcholines in blue, ions in pink and water as a transparent surface. Simulation box in blue. (b-f) Free energy profiles as a function of the collective variable CV, i.e. the fractional coordinate of the center of mass of the molecule along the direction normal to the membrane, from well-tempered metadynamics (wtmetaD) simulations of water and the SIn clusters of Figure 1 across the POPC membrane model. The evolution of the free energy profiles during the wt-metaD simulations is reported according to the corresponding color scale. In all cases the zero of free energy corresponds to the minimum of the final profile. 178x179mm (300 x 300 DPI)


Page 46 of 54

Page 47 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. (a) Mixed POPC + 20% cholesterol (POPC_CH) membrane model. Color code as in Figure 3. Cholesterol molecules in yellow. (b-f) Free energy profiles as a function of the collective variable from wtmetaD simulations of water and the SIn clusters of Figure 1 across the POPC_CH membrane model. Details as in Figure 3. 177x177mm (300 x 300 DPI)



Figure 5. (a) Heterogeneous red blood cell-like (RBC) membrane model. Color code as in Figure 4. Phosphatidylcholines in blue, phosphatidylserines in red, phosphatidylethanolamines in green, sphingomyelins in magenta. (b-f) Free energy profiles as a function of the collective variable from wt-metaD simulations of water and the SIn clusters of Figure 1 across the POPC_CH membrane model. Details as in Figure 3. 180x178mm (300 x 300 DPI)


Page 48 of 54

Page 49 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Deuterium order parameters (SCD) of each atom of the palmitic acyl tail of POPC, for both the homogeneous POPC (full symbols) and the mixed POPC + 20% cholesterol (empty symbols) membrane models. In both cases, results are reported both in absence of silica (red circles) and when a silica particle is at the center of the lipid phase, in this case comparing lipids close to the particle (at less than 15Å from any silica atom, orange triangles) and lipids far from the particle (at more than 20Å from any silica atom, blue squares). Error in the SCD calculation is visualized as a shaded area around the curves. 163x131mm (300 x 300 DPI)



Figure 7. (a) Distribution of the number of H-bonds between the silanols (hydroxyl) groups of the SI6b silica cluster and the phosphate groups of the POPC phospholipids in the homogeneous POPC model. Data are computed tacking in to account the frames of the wt-metaD simulations in which the silica particle is in the middle of the water-lipid interface. (b) Comparison between the distributions of the total number of H-bonds made by the SI6b silica cluster when at the interface with the POPC model (red) and when in the water phase (blue). Lines are Gaussian fittings to help the comparison. 176x83mm (300 x 300 DPI)


Page 50 of 54

Page 51 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8. (a) Snapshots displaying silica clusters (depicted as vdW spheres) of increasing size inside the lipid bilayer. To show the formation of “holes” in the membrane, the lipids are displayed as a transparent isosurface and water molecules (ball and stick models) close to the lipids are also included. (b) Radial pair distribution function (ga,b) between oxygens of silica and water averaged along all frames in which the silica particle (SI6b, in this example) is inside the POPC lipid bilayer. (c) Magnification of the SI6b case of Figure 8.a, to better show the co-translocation of water molecules inside the membrane. 152x161mm (300 x 300 DPI)



Figure 9. Snapshots displaying the crossing of the POPC lipid bilayer by the SI1 silica cluster. The water and lipid phases are visualized as transparent isosurfaces (cyan and white, respectively), while the silica cluster and water molecules close (< 5Å) to it are reported in ball and stick models. The corresponding value of the CV is also provided. 342x159mm (300 x 300 DPI)


Page 52 of 54

Page 53 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 10. Time-averaged deformation maps for the top and bottom leaflets, reported for one event in which the SI9b cluster is inside the POPC membrane. Deformation is defined as the interpolated location of the head-groups after translating the chosen leaflet's center at z = 0. Therefore, deformation takes on both negative and positive values. Irrespective of the leaflet, deformation is positive when it expands the membrane, and negative otherwise. The particle location corresponds to the area of maximum negative deformation of the top leaflet, that is the side from which it entered the membrane. 395x178mm (300 x 300 DPI)



Table of Content 101x49mm (144 x 144 DPI)


Page 54 of 54

Molecular Dynamics Simulations of the Silica â Cell Membrane

Recommend Documents

Molecular Dynamics Simulations of the Silica â Cell Membrane