extfieldsofstudy:
Materials Science
Medicine

text:
Implicit-solvent dissipative particle dynamics force field based on a four-to-one coarse-grained mapping scheme

A new set of efficient solvent-free dissipative particle dynamics (DPD) force fields was developed for phospholipids and peptides. To enhance transferability, this model maps around four heavy atoms and their connected hydrogen atoms into a coarse-grained elementary bead based on functional group. The effective hybrid potential between any pair of beads is composed of a short-range repulsive soft-core potential that directly adopts the form of an explicit-solvent DPD model and a long-range attractive hydrophobic potential. The parameters of the attractive potentials for lipid molecules were obtained by fitting the explicit-solvent DPD simulation of one bead of any type in a water box, then finely tuning it until the bilayer membrane properties obtained in the explicit-solvent model were matched. These parameters were further extended to amino acids according to bead type. The structural and elastic properties of bilayer membranes, free energy profiles for a lipid flip-flop and amino acid analogues translocating across the membrane, and membrane pore formation induced by antimicrobial peptides obtained from this solvent-free DPD force field considerably agreed with the explicit-solvent DPD results. Importantly, the efficiency of this method is guaranteed to accelerate the assembly of vesicles composed of several thousand lipids by up to 50-fold, rendering the experimental liposome dynamics as well as membrane-peptide interactions feasible at accessible computational expense.

Introduction
Simulations of molecular dynamics (MD) play crucial roles in interpreting molecular interactions and assembly, particularly when they are unreachable experimentally because of limited time and length resolutions [1][2][3]. Examples encompass the assembly of a plethora of structures in soft matter and biological systems. In principle, all-atom (AA) MD captures the most full-scale information but provides access only to systems containing several hundred molecules on a nanosecond time scale because of the copious amount of computation [2,3]. For actual soft matter and biological systems such as biomembrane transformations (e.g., phase transition between gel and fluid, raft formation, microtube growth, vesiculation, and membrane fusion), the scale is far beyond the abilities of AAMD [4,5]. This hurdle between atomic resolved simulation and natural chemistry can be lessened by coarse-graining several atoms into one particle [4,5]. The traditional explicit-solvent coarse-grained (CG) method, however, continues to lack sufficient speed for solving the majority of membrane-associated processes. In these cases, the computational burden is related to the bulk solvent and becomes heavier when the membrane deforms to a great degree or the spatial size of the simulated system increases markedly [6,7]. In one example, in a CG simulation of a tether pulling from a membrane, the simulated system had as many as 4 million CG particles, of which more than 90% were CG water beads [6]. In another, a spherical vesicle of radius R was simulated. The fraction of its sheet-like membrane area (R 2 ) in the simulation box (L 3 ) was approximately R 2 /L 3 % 1/R, which decreased rapidly as a function of R [7]. By neglecting the solvent degrees of freedom, a significant enhancement in calculation efficiency can be achieved. With this consideration, various implicit-solvent force fields have been developed over the past decade. Milano and coworkers introduced a hybrid particle field MD scheme with compatibility across AA and/or particular CG models [8]. It was validated over various molecular models, including biomembranes [9]. The method combined a discrete representation of the molecules (AAMD or CGMD) with a continuous self-consistent field description of the inter-molecular interactions. A model such as this requires storage for all the particles' coordinates, and it loses a direct route to introduce implicit solvent. Wang and Deserno explicitly introduced a cohesive potential to mimic the hydrophobic attractions between the lipid tails in a one-component membrane in the fluid phase [10] and applied it to other similar, but simple, CG lipids by changing the topology of the lipid [11]. The model has poor transferability to complicated lipids. This transferability issue could be significantly solved by choosing a versatile and straightforward parameterization for the potential related to the implicit solvent. In this situation, the Newton's inversion method used by Lyubartsev et al. [12] or the force matching method employed by the Voth group [13,14] could be useful, but they remain insufficient at this moment because of the impractical procedures used to obtain potential and/or parameters, wherein the effective CG forces were extracted from the radial distribution functions or the effective non-bonded forces between selective CG interaction sites from a reference AA simulation.
Because chemical details matter in many cases (e.g., complicated peptide-membrane interactions, local membrane structural transformation induced by mini environmental stimulus, and the chemically heterogeneous Gram-negative bacteria membrane), a sufficiently high CG resolution is required. The CG Martini model provides a force field that considers transferability and chemical specificity [6,[15][16][17][18][19][20]. This model maps about four heavy atoms and their connected hydrogen atoms into a coarse-grained elementary bead. These beads are divided by charge type, polarizability, and hydrogen-bonding capacity. Embedded in the free-source, industry-wide, and constantly updated software GROMACS [21][22][23][24][25], Martini was developed to investigate a large variety of biological processes. The accuracy of Martini is guaranteed by the thermodynamic data calculations for chemical building blocks, usually representative small compounds, matched with experimental values [15][16][17][18][19][20]. In addition, a broad range of applications involving phospholipids [15], proteins [16,18], carbohydrates [17], glycolipids [19], DNA [20], etc. are realized with no reparameterization, benefitting from the extensively validated building blocks. For the sake of speed, the solvent-free version of Martini [26], nicknamed "Dry" Martini, was implemented based on the original Martini force field, inheriting the building block feature. The Dry Martini force field, however, continues to impose a severe computational requirement because of the "hard" Lennard-Jones (L-J) potential, which limits the time step used in Newton's equation of motion in the series of GROMACS-based force fields.
Compared to the hard sphere models in most CGMD force fields, including Martini, the dissipative particle dynamics (DPD) model maps atoms into soft beads that can overlap . These soft beads interact via short-ranged conservative, random, and dissipative forces. The soft potential increases the time step considerably, speeding up the simulation and making it easy for beads to escape local minimums formed by surrounding particles. In addition, DPD explicitly includes fluctuation and dissipation stemming from the neglected degrees of freedom. More importantly, the latest set of DPD parameters [45], updated by adopting the Martini-like four-to-one mapping scheme [15], has strengthened portability. It matches the correct water compressibility and the experimental Flory-Huggins χ parameter. Membrane properties simulated by this set of parameters are consistent with atomistic simulations [50][51][52] and experimental measurements [53][54][55][56].
A DPD model that includes an implicit solvent (Im-DPD) has been developed for studying ethanol-water mixtures [57], block copolymers [58,59], and liquid drops surrounded by a gas [60]. Silbermann et al. extracted effective pairwise additive CG potentials from AAMD simulations, but which lacked in an analytic form [57]. Panagiotopoulos's group introduced L-J type implicit potentials, which, however, were too "hard" [58,59]. Liu and Meakin combined an analytic soft-core smoothing function with DPD to simulate binary mixtures, Poiseuille flow, and liquid drop formation [60]. However, the implicit solvent potential has not been developed and applied to biomembrane systems. Only in recent years has a solventfree DPD model for lipid membranes been available [61,62]. Sevink and Fraaije rigorously derived the implicit solvent potential by systematically averaging the solvent degrees of freedom via a hybrid particle-continuum method. A long-range attractive hydrophobic potential with an analytic form was obtained. Three parameters in this additional potential were tuned to match the membranes and vesicles, as simulated by DPD with an explicit solvent (Ex-DPD) [31,38]. A 20-fold decrease in required calculational time was achieved by implicitly describing the bulk solvent. However, this implicit-solvent force field was based on a very simple CG model of lipid molecules and less transferable Ex-DPD force parameters [31,37,38]. Applying this to the Martini-like DPD beads with chemical specificity resulted in implicit-solvent force parameters that were difficult to tune using the methods of Ref. [61,62]. We have developed a more direct routine for constructing the Im-DPD potentials. In our method, two parameters were obtained by fitting the explicit-solvent force acting on a CG bead immersed in a pure water box, and one was finely tuned by matching the structure of a lipid bilayer membrane. The validity of this new solvent-free force parameter set was justified by simulating essential structural and elastic properties of bilayer membrane, partition energy of a lipid molecule or amino acid analogues into a bilayer membrane, and pore formation induced by antimicrobial peptides.

Theory of Im-DPD
In Ex-DPD simulation, for each pair of beads i and j separated by a distance r ij < r 0 , a conservative force a dissipative force and a random force are combined to describe their interactions. Here, r ij = r i -r j , r ij = |r ij |, ȓi j = r ij /r ij , and r 0 is the interaction cutoff. The vector v ij = v i -v j is the velocity difference between beads i and j. The DPD force parameter a ij (in units of k B T/r 0 , where k B is the Boltzmann constant and T is the temperature) represents the maximum repulsion strength. The parameter γ ij represents the friction coefficient (in units of ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B Tm 0 =r 0 2 p , where m 0 is the bead mass). The random force parameter z ij is a symmetrically and uniformly distributed random number.
To average the solvent degrees of freedom, Sevink and Fraaije [61] employed the hybrid particle/continuum method where the hybrid (free) energy for particle k and one monomeric (solvent) field ρ is given by The first term is the usual potential energy in a particle description. The second term is the usual free energy for the field description. Its form is shown in where Λ is the de Broglie thermal wavelength, and κ H is the Helfand compressibility. The final term is the coupling between the particles and the field, defined by A normalized Gaussian smoothing function is commonly adopted for the kernel function Here, the scalar spread σ and coupling parameter c k respectively describe the spatial extent and strength of interaction between the field and particle k. Assuming the solvent is always in a local equilibrium on a coarse-grained time scale, after a complicated derivation [62], the solvent-free potential between two particles i and j was obtained Here,c l ¼ s i c i , and s i is a scaling parameter. Besides ordinary particle-particle interactions, an additional attractive potential describing hydrophobic interactions is seen.
To describe the attractive potential, three variables must be parameterized for any type of CG bead, namely, the spreading width σ, the particle-field coupling strength c, and the scaling parameter s. Sevink and Fraaije constructed the parameters for a lipid bilayer membrane using a relatively simple CG lipid model wherein the CG lipid molecule took on a λ-shaped polymer H 3 (C 4 ) 2 composed of three head (H) particles and two hydrocarbon (C) tails of four particles each (see Fig 1A). They first explicitly obtained the field-particle coupling parameter c from existing DPD force parameters (a ww = 25) for water-water interaction and the Helfand compressibility (κ H = 4.6) via equation where a = a ww + Δa. The DPD repulsive force parameter a for the H 3 (C 4 ) 2 lipid used in Ref. [31] is given as {a WW , a WH , a WC } = {25, 35, 75}, offering {c H , c C } = {5. 55, 9.37}. The rescaled coupling parameterc (or s) and spread σ were then obtained by applying the iterative Boltzmann inversion procedure to determine the effective interaction potential after removing a fraction of the particles. For Δa = 0, the fittedc in the concentration range of interest for liposome simulations was roughly equal to 4.6, the same as was found for the unscaled value for the reference system, suggesting that the previously determined additional scaling factor s was not necessary. The spread (σ = 0.4) was the same across the concentration range. In practice, for Δa 6 ¼ 0, Sevink and Fraaije suggested setting s = 1 and manually tuning σ to match the macroscopic membrane properties. They found membrane integrity to be most sensitive to σ C , the spread of the lipid tail bead, but membrane stiffness was sensitive to σ H , the spread of the lipid head bead. They showed that the determined membrane properties (area per lipid and bending rigidity) for the model with an explicit solvent are best reproduced for {σ H , σ C } = {1.11, 1.14}. However, when we tried to apply this method to a more complicated CG lipid model, such as a Martini-like h-shaped lipid (cf. Fig 1B-1D), it was difficult to tune the parameters to reproduce membrane properties. We therefore developed an alternative procedure for creating Im-DPD force parameters.

Parameterization of Im-DPD potential for Martini-like CG lipids and peptides
We adopted a Martini-like CG model by mapping around four heavy atoms and their connected hydrogen atoms into one bead [47][48][49]. A lipid molecule was modeled as an h-shaped polymer connected by harmonic bonds and consisting of four hydrophilic head beads and two hydrophobic tails. The atomic representations and corresponding CG mappings of dimyristoyl phosphatidylcholine (DMPC), dimyristoyl phosphatidylglycerol (DMPG), and dioleoyl phosphocholine (DOPC) lipids are shown in Fig 1B-1D. The CG model for dipalmitoyl Implicit-solvent DPD force field phosphatidylcholine (DPPC) lipid was constructed by adding one more bead to each of the two tails in the CG DMPC model. An amino acid residue was represented by one backbone bead and one or more side-chain beads, as shown in Fig 2. Based on functional group, as in the Martini model, the DPD beads were sorted into charged (Q), polar (P), nonpolar (N), and apolar (C) types [15]. Each type was further divided into subtypes based on hydrogen donor capacity (d), hydrogen acceptor capacity (a), or no hydrogen bond forming capacity (0). This CG model retains the essential chemical specificities of the molecules. Based on this four-toone mapping scheme, we recently constructed a new set of Ex-DPD force fields [45], shown in Table 1. To reproduce the correct compressibility of water, a larger value, a WW = 100, was assigned to beads of the same type. The repulsion parameters between beads of different types Implicit-solvent DPD force field were obtained by matching to the Flory-Huggins χ parameter. Compared to the λ-shaped lipid model, this Martini-like DPD model allows for more types of CG beads. Therefore, an efficient and reliable approach is required for building an Im-DPD force field according to this CG mapping.
We first optimized parameters c and σ, determining the coupling between the solvent field and the particles. For the special condition where only one particle is immersed in a homogenous and equilibrated solvent box, the density of the entire particle-field system is approximately equal to the density of the field ρ, which is 3/r 0 3 in this model. Thus, the force acting on the particle by a slice of solvent can be obtained analytically using Eqs (6) and (7) and has the form Here, z is the distance between the particle and the center of the solvent slice, and z 2 -z 1 is the thickness of the slice. This force can also be obtained using Ex-DPD simulation, as a sum of the forces between the particle and all the water beads of the slice, within a cut-off of 1.0 r 0 . To calculate the force, Ex-DPD simulations governed by Eqs (1)-(3) were performed in a constant volume and constant temperature (NVT) ensemble and periodic boundary conditions using the velocity Verlet algorithm. Here, the reduced temperature T was set to 1, corresponding to 298 K. The time step was set to 0.
). The friction coefficient γ ij remained consistent at 4.5 for any pair of beads. In this four-to-one CG model, r 0 corresponds to 0.71 nm in physical units [28,45]. Each simulation was performed for up to 2×10 5 time steps. The final 10 5 time steps, containing 1000 equally separated frames, were used for subsequent analysis. The force acting on the immersed bead was calculated by averaging all forces between this target bead and water beads inside the slice of thickness 4r 0 . The force on a bead, where Δa = a iw -a ww = 30, as a function of distance z, is shown in Fig 3A. This curve was well fitted by Eq (10) and gave rise to the pair {c, σ}. The forces for beads along a series of a iw ranging from 90 to 130 were calculated to obtain the corresponding {c, σ}. As shown in Fig 3B and 3C, both increased almost linearly as a functions of Δa. This is represented by equations and s ¼ 0:5388 þ 9:564 Â 10 À 4 Da: ð12Þ From this fitting, we extracted the {c, σ} parameter set for all DPD bead types considered ( Table 2). As in Eq (9), we set the Helfand compressibility to κ H = 19.69 for this four-to-one CG model where a ww = 100, almost four times the compressibility κ H = 4.6 used in the Sevink and Fraaije model, where a ww = 25. The value of σ was only half of that used by Sevink and Fraaije. We propose that σ, when in the range of 0.54-0.57 r 0 , is reasonable because in traditional DPD simulation, the dimensionless interaction cutoff r 0 is only 1; σ cannot exceed this cutoff. The scaling parameter s was introduced to yield the equilibrium density when particles have assembled to a condensed phase. This implies that s must be related to the system's macroscopic properties. Therefore, we optimized s by first simulating the self-assembly process of  DMPC vesicles then examining the structural and elastic properties of a planar bilayer membrane. A planar membrane was prepared by placing N lipid = 1152 DMPC lipids into a grid spanning the x-y plane with their head groups facing outward and tail groups facing the membrane center. The three box dimensions (L x = L y = L z ) were determined by L x = (0.5N lipid a prj ) 0.5 , where a prj is the projected area per lipid at initial setup. To simulate vesicle assembly, a planar membrane inside a large box of 1.5L x ×1.5L y ×1.5L z was prepared. A box this large can provide sufficient room for phospholipids to escape their periodic images, allowing the membrane to curl up along the edges with vesicle growth. The integration algorithm in Im-DPD is same as in Ex-DPD. The implicit-solvent forces governing the motion of the DPD beads are described by and The Im-DPD forces were pairwise as in Ex-DPD and conserved the overall momentum, in sharp contrast to typical Brownian Dynamics. Specifically, the cutoff r c im was equal to 2.5 r 0 .
A carefully tuned γ ij im was assigned a value of 0.8 to maintain the temperature at 1.0 k B T. The random number z im ij differs from z ij in usual DPD forces. In a lipid molecule or peptide, the molecular skeleton was maintained via harmonic potentials, including 2-body bonding interactions and 3-body angular constraints. These are described by and where K 2 and K 3 are the spring constants of the bond strength and angular stiffness, respectively, and L 0 and θ 0 are the equilibrium bond length and angular degree, respectively. The equilibrium CG bond lengths and angles and related force constants can be obtained by fitting the bond distributions into AAMD simulations [63]. Tables 3 and 4 list these constants for the typical lipid molecules and peptides.
To maintain the secondary structures of polypeptides, dissociable Morse potentials were applied to mimic hydrogen bond interactions between beads of a peptide skeleton that is helical or sheet-like. This is described by For a α-helical structure, MB 1-3 Morse bonds were introduced between every two backbone beads separated by two normal harmonic bonds, and MB [1][2][3][4][5] Morse bonds were introduced between beads separated by four bonds. The equilibrium distance r M and depth of the potential well K M were set to 0.60 r 0 and 3.0 k B T, respectively, for bonds MB 1-3 and 0.78 r 0 and 6.0 k B T, respectively, for bonds MB [1][2][3][4][5] . In both cases, the width of the potential well α was set to 6.4 r 0 , and the cutoff for the Morse potential was 2.0 r 0 . To stabilize a β-sheet structure, a Morse potential with K M = 6.0 k B T, r M = 0.5 r 0 , and α = 6.4 r 0 for every two hydrogen-bonded beads was applied. Simultaneously, a harmonic angle potential θ 0 = 180˚was applied for every two backbone beads separated by two bonds. If both beads of a pair were charged, electrostatic interactions were calculated empirically [32] Here, r c = r 0 = 0.71 nm, representing the length unit. The electrostatic smearing radius r e was set to 1.6 r 0 . At 298 K, the coupling constant G ¼ e 2 k B T 0 r r c ¼ 65:78 [32], where 0 is the dielectric constant of a vacuum, and r is the relative permittivity of water at room temperature. Because the water was implicit, the relative permittivity constant r was only 15.0 [26] (in Ex-DPD, it was 78.3 for water).
The parameter s was modulated so that the assembled structures of lipids could be examined. Initially, s was made equal for both head beads and tail beads. As illustrated in Fig 4, when s was small, the lipids were unable to assemble, but when it was large, tail particles  densely packed inside the core of a hard spherical shape. Only when s ranged between 0.7 and 0.8 could the lipids self-assemble into a vesicle structure. However, the area per lipid of a tensionless planar bilayer a prj 0 became 1.091 r 0 2 (detailed definitions of a tensionless bilayer and its area are forthcoming), much less than the 1.355 r 0 2 obtained in the Ex-DPD simulation [45]. To correctly reproduce the membrane area, s was finely tuned for the lipid head and tail separately. For simplicity, DPD beads were sorted only into hydrophilic (where a iw is close to a ww ) and hydrophobic (where a iw is much larger than a ww ) types while parameterizing s. After a series of simulations, a well-ordered vesicle and a tensionless bilayer with a prj 0 close to 1.355 r 0 2 ( Fig 5A) were found when {s H , s C } ranged between {0.46-0.54, 0.730-0.734}. To refine the parameters more precisely, the elastic response of a membrane to stress (i.e., the membrane tension as a function of area per lipid) was investigated. As shown in Fig 5B-5D . These values for s were assigned to amino acid polypeptide beads based on their hydrophilicity and hydrophobicity ( Table 2).

Applications of Im-DPD
This new Im-DPD force field was applied to investigate the structural and elastic properties of three types of bilayer membranes, the free energy cost of moving a target molecule (a lipid molecule) or a molecular fragment (an amino acid analogue) into a bilayer membrane, and  Implicit-solvent DPD force field antimicrobial peptide (AMP)-induced membrane pore formation. The efficiency of this method was examined by simulating the self-assembly process of large vesicles.

Bilayer properties
The bilayer's tension as a response to stress was investigated first, then the membrane thicknesses, projected lipid headgroup areas, group density profiles, lipid order parameters, and bending rigidities of tensionless DMPC, DPPC, and DOPC bilayers were calculated. To evaluate membrane tension as a function of the bilayer's stretched area, planar membranes with fixed numbers of lipids but different initial projected areas per lipid a prj (i.e., a different lateral area L x ×L y ) were simulated. In equilibrium, the tension ∑ of the bilayer was estimated via an integrated stress profile [31], where P x , P y , and P z are the pressure components in three directions. A membrane with ∑ = 0 Implicit-solvent DPD force field is said to be in the tensionless state. The corresponding projected area per lipid in the tensionless state is a prj 0 . The membrane thickness L mem was obtained using the average distance between the choline groups in the upper and lower leaflets. The density distribution of relevant functional groups such as lipid head groups and lipid tail groups, perpendicular to the membrane (along the z-axis), was averaged along the x and y axes. The orientation order of the alkyl tails S chain is defined in where θ is the angle between the normal bilayer plane and the orientation of the vector along the hydrocarbon chain. The square bracket symbolizes an ensemble average.  Table 5. Extreme stretching caused the membrane to rupture at a critical tension ∑ rup , also listed in Table 5.
Once the membrane ruptured, its tension fell rapidly. Even though the values for tensionless a prj 0 were different for DMPC, DOPC, and DPPC (T = 1.2) bilayers, their tension curves were similar and increased monotonically with membrane area before rupture, indicating they were in fluid phases. In contrast, the membrane tension for a DPPC bilayer simulated at T = 1.0 Implicit-solvent DPD force field exhibited different elastic responses to stress. Near the tensionless state, tension increased quickly as the bilayer was stretched, as seen by a sharp slope in the curve. This indicated that the bilayer was in a gel phase. When the membrane area was stretched more than 7%, tension began to decrease. Once it was stretched close to 13%, tension increased again. Corresponding snapshots show that in the abnormal region, a local fluid phase (see Fig 7B) appeared in the bulk gel phase. The lipids in this local fluid phase were less compact and relatively easier to stretch than those in the gel phase; thus, the tension curve of the mixed fluid-gel phase membrane decreased when membrane area increased. The phase transition process quickly finished upon sufficient stretch, after which the DPPC bilayer behaved like a normal fluid-phase membrane. This direct membrane rupture upon stretching for fluid phase bilayers and phase transition under stress for gel phase bilayers were in good agreement with Ex-DPD results [45]. More detailed snapshots and density distribution profiles of tensionless DMPC, DOPC, DPPC (T = 1.0), and DPPC (T = 1.2) bilayers are shown in Figs 8 and 9, respectively. As observed in Ex-DPD simulations [45], in the fluid phase, the lipid hydrophobic tails randomly spread out, and the terminal beads from opposite leaflets can touch each other. On the other hand, in the gel phase, the alkyl tails arranged like a hairbrush, and two leaflets greatly separated. The two distinguishable phases were further demonstrated by their order parameters of alkyl tails S chain as shown in Table 5. The gel phase lipids arranged orderly and therefore had larger values for S chain (as high as 0.9), but in the fluid phase, the tails beads were able to swing more freely and thus had slightly lower values for S chain (less than 0.7).
A noticeable difference between Im-DPD and Ex-DPD models was the considerable vacuum (indicated by a minimum in the center of the tail bead distribution profile in Fig 9) arising at the centers of the membranes in Im-DPD simulations. In contrast, in Ex-DPD simulations [45], only shallow minima were observed for the gel phase DPPC (T = 1.0) bilayer. The tail bead distribution profiles for fluid phase DMPC, DOPC, and DPPC (T = 1.2) membranes were smooth. Based on these phenomena, the implicit-solvent bilayers were slightly thicker and had larger areas a prj 0 than explicit-solvent bilayers, implying that implicit-solvent bilayers are less compact. This vacuum was frequently observed in atomistic membrane [50,51,[64][65][66] and coarse-grained membrane [9,10,13,26,67]. The lack of vacuum in Ex-DPD is due to the short-range purely repulsive forces for every two interacting beads. It brings the Implicit-solvent DPD force field lipids to assemble into compact structures to avoid the penetration of water. However, longrange L-J interaction in AAMD and other CGMD methods as well as a solvation term explicitly included in Im-DPD take into account the hydrophobic attractions, allowing less compact lipid packing. In this way, Im-DPD even compensates the imperfection of Ex-DPD. Elasticity is another important property of lipid bilayers that influences cell growth, cell division, cell apoptosis, and so on. The fluctuation spectra method was used to identify the bending modulus of a membrane. In this method, the instantaneous positions of the choline groups were chosen, thus defining the two surfaces of the bilayer. The local height of the membrane u(x,y) is defined by where Z 1 (x, y) and Z 2 (x, y) are the heights of the upper and lower surfaces, respectively. This bilayer was subsequently mapped to a coarser discrete grid with spacing h (h is typically assigned a value slightly larger than the membrane thickness). The discrete height u i,j represents the average value of u(x, y) belonging to lattice (i, j). The Fourier expansion of the bilayer height is defined as where q is the wave vector. The Fourier coefficients u(q) are given by In the lattice model, r = hn, n = n x i + n y j, and q ¼ 2pn , and À L y 2h n y < L y 2h . Radially averaging |u(q)| 2 over q = |q| combined with time averaging over 1000 frames yields the fluctuation spectrum This fluctuation spectrum is related to the membrane tension S and bending rigidity κ via A larger membrane, where L x = 2L y = L z % 64r 0 , was prepared to guarantee sufficient membrane fluctuation. The membrane was then relaxed for up to 2×10 5 time steps in the NVT ensemble. The final 1000 coordinates were collected for fluctuation analysis. A typical fluctuation spectrum for a tensionless DMPC bilayer is shown in Fig 10. By fitting Eq (26), S = 0.1   Table 5, were consistent with Ex-DPD results [45] and experimental data [53][54][55].

Potential of mean force for flip-flop of a lipid molecule
A smaller bilayer containing 512 DMPC (or DOPC or DPPC) lipids was prepared to calculate the potential of mean force (PMF) for lipid flip-flop. By using the umbrella sampling method [68], a harmonic biasing potential was employed to restrain a lipid molecule (of the same type used for simulated membranes) at 61 continuous heights within the membrane, producing 61 serial simulation windows. The spacing between adjacent windows was set at 0.15 r 0 . The harmonic force with a constant of 110 k B T acted on the phosphate head bead of the lipid. After 10 5 time steps of equilibration for each simulation, the weighted histogram analysis method [69] was utilized to produce PMF profiles from the subsequent 10 5 time steps. The PMF profiles obtained from both Ex-DPD and Im-DPD are shown in Fig 11. Similar tendencies are seen: a free energy minimum appears near the head group where the confined Implicit-solvent DPD force field lipid molecule is in equilibrium, a steep slope is seen as the head group moves into the bilayer center, and a flat plateau or a sharp peak is seen at the bilayer center. As one can see, the PMF profiles exhibited wide plateaus in the hydrophobic core region in Im-DPD and sharp peaks in Ex-DPD. Actually, both shapes were observed in the AAMD [70] and CG Martini simulations [51].
The energy differences for flip-flopping a DMPC, DOPC, and DPPC in Ex-DPD simulations are 16, 17, and 34 kJ/mol, respectively. These increase to 39, 34, and 47 kJ/mol, respectively, in Im-DPD simulations. Dozens of available free energy barriers from AAMD simulations [52,66,70], CG Martini simulations [51], as well as experimental excess chemical potentials [71] were summarized in Table 6. Besides DMPC, DOPC, and DPPC, the energy cost for DLPC flip-flop was also listed. Compared to DMPC lipid, a DLPC lipid has shorter tails with less two carbons in each tail, while they both have the same CG structures in either DPD or CG Martini method. The overall trend of the free energy barriers can be summarized as Ex-DPD<Im-DPD<AAMD%Experiment<CGMD. It has been verified that, the PMF had a strong dependence on truncation radius [70], the treatment of electrostatic interactions [51], or even the force fields [70]. Taking the Berger DLPC lipid bilayer as an example, the increase of the cut-off from 0.9 nm to 1.4 nm brought about an increase of the energy barrier from 17 kJ/mol to 50 kJ/mol [70]. Another instance is the Martini membrane, wherein the normal water model estimated a free energy of 90 kJ/mol for DPPC flip-flop, which value increased to 160 kJ/mol in the polarized big multiple water model [51]. Even two AAMD lipids, e.g. Berger type and Charmm36 type, had a discrepancy as high as 40-50 kJ/mol [70]. Therefore, the differences between both DPD barriers and experiment references as well as AAMD estimates were within the realistic range. Both DPD underestimated the PMF barrier because they Implicit-solvent DPD force field employed soft-core short-range interaction. The Im-DPD energy barrier was slightly greater than that of Ex-DPD because of the thicker Im-DPD membrane.
Specific snapshots for confined DMPC and DPPC lipids at various depths were taken from the simulations so that their conformations in fluid and gel phases could be compared. As shown in Fig 12, at membrane center (0.00 r 0 ), the DMPC lipid became trapped in the vacuum region, and its tails were primarily oriented parallel to the membrane, but they could also interact with other lipids in both leaflets and have a variety of orientations. When pulled into equilibrium (3.00 r 0 ), an elongated conformation was gradually adopted, resembling surrounding lipids. Forcing the lipid out of the membrane environment (7.95 r 0 ) resulted in more random orientations. The DPPC lipid ( Fig 12B) adopted a similar conformation distribution near the membrane surface (3.60 r 0 ) and within the bulk implicit solvent (9.90 r 0 ), while distinctions were seen near the membrane center (0.00 r 0 ). One tail of the DPPC lipid was parallel and another tail was perpendicular to the normal membrane. In contrast, in Ex-DPD simulation, the two tails of DPPC splay to insert into two opposing leaflets separately. Once again, the larger vacuum in the bilayer center in Im-DPD allowed at least one tail to remain there.

PMFs for CG amino acid beads across membranes
To verify the validity of the Im-DPD force field for polypeptides, the free energy cost for transferring side chain beads of amino acids across a DOPC bilayer membrane was investigated. For simplicity, 18 amino acid side chains (alanine and glycine lack side-chains) were sorted into 8 groups according to charge, hydrophilicity, and polarizability. They were the SER group (serine, threonine, and cysteine), VAL group (valine, leucine, isoleucine, methionine, and proline), ASN group (asparagine and glutamine), ASP group (aspartic acid and glutamic acid), Implicit-solvent DPD force field LYS group (lysine and arginine), TRP group (tryptophan and tyrosine), PHE group (phenylalanine), and HIS group (histidine). Fig 13 presents the PMFs obtained from both Ex-DPD and Im-DPD simulations. The available PMF data from AAMD simulations are shown for comparison [64,65]. To account for differences in the thickness of the membrane, the thickness of Ex-DPD and Im-DPD membrane was enlarged or shrunk to be consistent with AAMD membrane. Energy costs for transferring amino acid analogues from the bulk solvent to the membrane center are given in Table 7. The PMF profiles for polar uncharged ASN and SER groups (Fig 13A and 13B) showed energy peaks in the center of the bilayer and energy troughs in the bulk solvent, with a difference of less than 10 kJ/mol (e.g., 31.41-23.85 = 7.56 kJ/mol), providing a good agreement between atomistic and both DPD force fields. Using AAMD simulation, a free energy barrier appeared in the proximity of the lipid head group, implying the amphipathic nature of ASN and SER side-chains. However, in the coarse-grained method, ASN and SER side-chains were modeled simply as hydrophilic beads, suppressing any energy barrier.  The PMFs for charged residues (LYS and ASP) showed relatively large differences between AAMD and DPD simulations (Fig 13C and 13D). The free-energy penalties for burying the residues inside the membrane by the atomistic force field were much larger than those found using either Ex-DPD or Im-DPD force field. The CG bead model carrying a single charge underestimated the hydrophilicity of a charged residue.
As shown in Fig 13E and 13F, the aromatic residues of TRP and PHE entered into the hydrophobic membrane center more easily in DPD simulations. In the DPD model, TRP and PHE side-chains were represent by 2 to 3 C-type beads, which are prone to be compatible with lipid tails. However, a barrier in the center of membrane in the AAMD simulation indicated the difficulty in maintaining a large ring residue even though it is hydrophobic. The CG model does not reproduce this characteristic.
In the process of transferring a VAL into the bilayer center, all three models were thermodynamically admissible (Fig 13G). Although VAL has a PMF very similar to those of TRP and PHE, because all are primarily composed of C-type beads, it differs in that it is smaller, comprising only one CG bead. As a result, VAL had fewer PMF discrepancies between the DPD and AA models.
In the PMF profile of HIS, atomistic references were unavailable. In DPD representations, its side-chain has a ring structure linked by P 1 -, Q d -, and C-type particles, making it is amphiphilic. As a result, the PMF curve had a minimum at the hydrophilic-hydrophobic interface of the bilayer.
Overall, the free energy of transferring an amino acid side-chain across a lipid bilayer obtained using the Im-DPD model was satisfactory, reproducing the Ex-DPD results well. However, a large discrepancy as high as 25 kJ/mol was obtained for the aromatic ring residue, e.g. TRP, PHE, and HIS, while 6 to 13 kJ/mol was for the rest, in terms of the difference between two DPD methods. We speculate this was of the same reason as between Martini and Dry Martini [26]. No matter in Martini or DPD force field, the aromatic beads partly overlap. Thereupon, the implicit description for water was overestimated, unwantedly improving the hydrophilicity of aromatic rings, also bringing larger energy costs. As mentioned already, the PMF data had strong dependence on the force fields and simulation conditions. For most amino acids, the barrier differences obtained from Ex/Im DPD simulations were within the realistic range compared to AAMD estimations. Nevertheless, room for improvement remains. For example, a more precise scaling factor s could be introduced for specific CG beads.

Antimicrobial peptide-induced membrane pore formation
Antimicrobial peptides are small and electropositive peptides with antimicrobial activities against many pathogenic microorganisms such as fungi, bacteria, etc. DPD simulation at the mesoscopic level provide a better understanding of the relationship between AMP structures and their actions. For example, we recently employed an Ex-DPD method to simulate the interactions between membranes and a couple AMPs with various secondary structures [49]. In this work, we investigated membrane pore formations induced by magainin and melittin as a practical application of Im-DPD.
Similar to the work of Ref. [49], the secondary structures of magainin and melittin were maintained as helixes to resemble their conformations upon binding to the membrane surface. A planar bilayer composed of 1682 zwitterionic DMPC and anionic DMPG lipids (DMPC:DMPG = 7:3) was used to mimic a bacterial membrane. A series of simulations were performed at peptide/lipid (P/L) molar ratios of 0.5%, 1%, 2%, 3%, 4%, and 5%. Counterions were loaded to neutralize the system. Results were consistent with those found using Ex-DPD simulation [49]. As shown in Fig 14, when the magainin concentration was low, the absorbed AMPs induced membrane buckling to a small degree. When more peptides were bound onto the membrane surface, the buckling was considerable, and the lipids were strongly disordered. To relax the peptide-induced compression (or tension), one or two magainin peptides were inserted into the membrane to form a small pore. Additional peptides and lipid head groups were added to enlarge the pore. Melittin peptides also induced pore formation in membranes with sufficient absorptions (Fig 15). Because the charge distribution along the helix was more discrete, magainin induced larger pores composed of 5 to 6 peptides. On the other hand, the charge on melittin was located primarily near the C-terminus, inducing the formation of smaller pores (each of 2-4 peptides). The shape and size of the peptide-induced pore and the critical peptide concentration necessary to induce pore formation in Im-DPD were good reproductions of those found using Ex-DPD [49] and AAMD [72,73], as listed in Table 8. Implicit-solvent DPD force field

Efficiency of Im-DPD
To demonstrate the efficiency of Im-DPD, we simulated vesicle formation starting from planar bilayers composed of 1152, 2048, 3025, and 4096 lipids. The time evolution of small vesicle formation is shown in Fig 16. A planar pre-assembled bilayer quickly transformed into a circular plate-like bicelle to minimize its circumference, followed by bending, wrapping up, and closing into a vesicle. We noted several detached free lipids 'flying' in the simulation box ( Fig 16B). This is a normal phenomenon in implicit-solvent CG simulations [62,77]. These lipids followed a ballistic trajectory because of the low local particle density [62]. The bilayer of 1152 lipids required 8.89 CPU hours to assemble a vesicle using Im-DPD; it required 95 CPU hours using Ex-DPD. Around 90 percent of the simulation time was saved when the solvent degrees of freedom was omitted. More remarkably, for larger vesicles composed of 2048, 3025, or 4096 lipids, a 20-fold to 50-fold increase in speed was acquired (Fig 17), demonstrating the high efficiency of our Im-DPD force field. Implicit-solvent DPD force field

Limitations of Im-DPD
Traditional AAMD and CGMD simulations of lipid membranes employ the Lennard-Jones potential, which has a repulsive r -12 part to approximate the short-distance repulsion, as well as an attractive r -6 part to describe the van-der-Waals interaction. The Ex-DPD non-bonded https://doi.org/10.1371/journal.pone.0198049.g015 Table 8. Membrane pores induced by magainin and melittin: The threshold pore forming peptide concentration C pep , the pore model, the number of peptides inside a pore N pep , and the inner diameter D pore . DT denotes the disordered toroidal model [74], B denotes the barrel-stave model [75], and C denotes the carpet model [76].  Implicit-solvent DPD force field interaction is purely repulsive, rendering too compact membranes, which is relieved in Im-DPD by introducing an extra long-range attractive term, as evidenced from the density profile and the free energy barriers (Fig 9, Tables 6 and 7). However, the mismatching between DPD and AAMD still exists. To further match them, Im-DPD parameters need to be obtained from AAMD, rather than from Ex-DPD simulations. Perhaps a thorough solution would be the development of a brand new Ex-DPD, that strikes the proper balance between the soft-bead DPD potential and the hard-core Lennard-Jones potential. Another limitation arises from the criteria of parameterization, e.g., the elastic property for DPD, while the partition free energy for Martini. As a matter of fact, the preferred properties are well reproduced at the expense of weakly sacrificing other characters. For instance, the DMPC bending modulus estimated by Martini is two to three times (1.1~1.7 ×10 −19 J) [26] larger than experimental value (0.56 ×10 −19 J) [54] or DPD (0.4~0.5 ×10 −19 J). On the contrary, DPD is a weak reproduction of AAMD, in comparison to the rather good correspondence between Martini and AAMD, in terms of the flip-flop energy consumption (Tables 6 and 7).

Pore Properties
A common limitation in CGMD/DPD method is the unrealistic treatment of the conformational change of a (poly)peptide. Because of the losing atomistic information, a peptide cannot spontaneously and correctly rearrange its backbone to fit different circumstances, e.g., the membrane surface or the water environment. In the current version of DPD force field, a Morse potential is applied to mimic fixed secondary structures, as the dihedral potential is employed in Martini. Therefore, the considerable change of the peptide structure is not handled reasonably. To alleviate this problem, one could consider the elastic network model [78], or more precisely the polarized CG model [79].

Conclusion
An efficient and reliable solvent-free DPD force field for phospholipids and polypeptides was developed based on a popular Martini CG mapping rule to improve transferability. Implicit DPD parameters were obtained from explicit DPD simulations combined with complementary refinements that yielded better performance. The validity of this Im-DPD force field was justified by investigating structural and elastic properties of lipid bilayer membranes, free energy profiles for lipid flip-flop and amino acid analogues translocating across the membrane, and pore formation induced by the AMPs magainin and melittin. Overall results showed good agreement with Ex-DPD simulations, atomistic MD simulations, and experimental measurements, providing a versatile and reliable application. Most significantly, when used to simulate a large vesicle composed of thousands of lipids, nearly 99 percent of the computation time was saved by implicitly processing the bulk solvent. Therefore, we expect a wide range of applications for Im-DPD that will bridge the gap between the atomistic and mesoscopic levels.


text:
Clustering and separation of hydrophobic nanoparticles in lipid bilayer explained by membrane mechanics

Small hydrophobic gold nanoparticles with diameter lower than the membrane thickness can form clusters or uniformly distribute within the hydrophobic core of the bilayer. The coexistence of two stable phases (clustered and dispersed) indicates the energy barrier between nanoparticles. We calculated the distance dependence of the membrane-mediated interaction between two adjacent nanoparticles. In our model we consider two deformation modes: the monolayer bending and the hydroxycarbon chain stretching. Existence of an energy barrier between the clustered and the separated state of nanoparticles was predicted. Variation analysis of the membrane mechanical parameters revealed that the energy barrier between two membrane embedded nanoparticles is mainly the consequence of the bending deformation and not change of the thickness of the bilayer in the vicinity of nanoparticles. It is shown, that the forces between the nanoparticles embedded in the biological membrane could be either attractive or repulsive, depending on the mutual distance between them.

that the method of loading of 2 nm-sized gold NPs into lipid vesicles affects the NPs distribution within the biological membrane 34 . If the NPs and phosphatidylcholine lipids are coextruded, the NPs form a dense monolayer in the hydrophobic core of the vesicle membrane (Fig. 1A). In contrast, vesicles which are formed by extrusion and then dialyzed in the presence of NPs dispersed with detergent contain a membrane region with NPs and a membrane region without NPs (Fig. 1B). In these vesicles, so-called the Janus-like vesicles, the NPs aggregate and form clusters in the membrane regions that are rich with NPs.
The two stable states of NPs in the vesicle membrane may be explained by the existence of an energy barrier between aggregated and separated NPs mediated by a biomembrane 34 . The aim of this study is to test this hypothesis by estimation of the energy of a biomembrane induced by redistribution of NPs.
It was shown in the past that considering the membrane as continuum elastic medium may well explain interactions between membrane inclusions like transmembrane proteins 24,35,36 . Two key modes of membrane deformation caused by hydrophobic inclusion are the hydrophobic mismatch causing deformation of hydroxycarbon chains (stretching/compression) and the membrane bending of both membrane lipid bilayers 10 . Within this study, two rigid NPs of diameter r separated by the distance d are considered (Fig. 2). The local deformation of the membrane by the inclusions is then analysed by the variation of the membrane elastic energy which is increased by intercalation of NPs into hydrophobic moiety of membrane 28 . The local equilibrium shape of the membrane is determined as a shape with minimal elastic energy of the membrane in the deformed state. The variations of the intrinsic curvature C 0 , the bending constant κ b , and the compression-expansion constant κ c are performed to  Geometrical model of the lipid bilayer with two embedded NPs of radius r at distance d. The thickness of the phospholipid hydroxycarbon tail region in the unperturbed planar lipid bilayer is denoted by ξ 0 . NPs induce changes in the local curvature of the bilayer leaflets C and the change of the membrane thickness due to deformation of hydroxycarbon tails λ = ξ/ξ 0 .
SCiEnTiFiC REPORTS | (2018) 8:10810 | DOI:10.1038/s41598-018-28965-y study the effect of membrane properties on the elastic energy of membrane with the NP inclusions. The range of the membrane elastic parameters adopted in simulations is shown in Table 1. Figure 3 shows the dependence of the calculated membrane elastic energy on the distance between two NPs and the corresponding membrane elastic energy for various values of NP diameter. It can be seen in Fig. 3 that the condensed state of NPs is energetically more favorable than the separated state. Small nanoparticles require less elastic energy to embed into the lipid bilayer, in agreement with Wi et al. 10 .

Results
Further, our simulations confirms the hypothesis of Rasch et al. 34 , that there exists an elastic energy barrier between disperse and aggregate state of NPs embedded in the biological membrane (Fig. 3).
The bending contribution to the energy per lipid molecules is higher than the stretching contribution of hydroxycarbon tails (Fig. 4). Figure 5 shows the height of the energy barrier in condensed (Fig. 5A) and separated (Fig. 5B)    membrane bending rigidity κ b has the largest effect on the magnitude of the energy barrier. The higher the values of intrinsic curvature C 0 , membrane bending κ b and compression modulus κ c , the higher the energy barrier that the nanoparticles must overcome to separate (Fig. 5A) or to cluster themselves (Fig. 5B). The height of the elastic energy barrier is generally lower for transition from separated to clustered phase than in the opposite directionfrom clustered into separated state.

Discussion
It was experimentally observed that the lateral distribution of hydrophobic NPs in a lipid membrane vesicle depends on the method of loading of NPs into the lipid vesicle 34 . If the NPs are incorporated into the membrane of the already existing lipid vesicle by the dialysis process, they form NP clusters in the membrane indicating that the NP clustering is energetically favorable 37 . However, if the membrane is self-assembled with NPs, it was observed that NPs are distributed uniformly in the membrane 34 . Our results could explain the aggregation of NPs in the lipid bilayer prepared by dialysis observed by Bonnaud et al. 37 and Rasch et al. 34 . Insertion of the first NP disturbs the membrane locally. It is energetically more favorable for the next intercalated particles to be clustered with the other membrane-inserted particle (Fig. 5A) than to be dispersed. If the clustered nanoparticles are separated, e.g. by small thermal fluctuations, the membrane-mediated interaction attract them together. However, if the membrane is self-assembled together with NPs, the NPs are embedded in the membrane at a distance larger than the one corresponding to the membrane elastic energy peak in Fig. 3. In this case, membrane-mediated interactions cause repulsive forces between NPs (Fig. 3) and NPs remain unclustered.
In general, the membrane-mediated interactions may induce either attractive or repulsive forces between the biomembrane inclusions [23][24][25][26][27]29 . Considering the gradient of membrane elastic energy in Figs 3 and 4, we may conclude, that both cases are possible for small hydrophobic NPs embedded in hydrophobic core of membrane, depending on the mutual distance between NPs.
In agreement with our model, the recently published coarse-grained molecular dynamics study 9 also predicted the existence of energy barrier between two anionic ligand-coated NPs embedded in a cholesterol-free lipid membrane. The peak energy of the energy barrier was observed at roughly 8 nm distance between the NPs of diameter 3 nm. It corresponds to the distance d = 5 nm as defined in Fig. 2. Our simulations predict energy barrier peak at the distance d equal to 4 nm. However, the previous study 9 was focused on NPs coated with amphiphilic ligands and they are therefore not fully included in the hydrophobic core of the lipid bilayer. Also, the particles in this study 9 were negatively charged (ammonic), therefore they repel each other electrostatively.
Rasch et al., 2010 hypothesized that the energy barrier appears mainly due to the phospholipid tail stretching deformation. Our simulations show that instead of the deformation of membrane hydrophobic core, bending of the monolayer is the decisive factor influencing the membrane elastic energy barrier between two neighbouring membrane embedded particles (Fig. 4). Increase of the bending stiffness of the phosholipid monolayers considerably increases the bilayer energy associated with NP insertion in the membrane and the height of the energy barrier between the condensed (Fig. 5A) and the dispersed (Fig. 5B) state of NPs. The present study employs the number of assumptions to intentionally keep the model simple. Our main goal was to qualitatively illustrate the membrane elasticity-driven interactions between NPs in the biomembrane. We have adopted a continuum approach although the studied length scale is at molecular dimensions. Describing the biomembrane as elastic continuum is certainly an approximation. However, studies of protein channels 38 show that hydrophobic matching and membrane-mediated interactions can be well understood also within simple elasticity theory.
The model could be further upgraded by considering the lipid tilt 24,36,39,40 . May et al. 36 showed that the lipid tilt modulus has two energy contributions. The first one is included in the present model as the stretching of the hydrocarbon chains and reflects the loss of chain conformational freedom. The second contribution that was not considered in this work is the entropic contribution resulting from the constraints imposed by the tilt deformation on the fluctuations of the hydroxycarbon chain. We may assume that the entropic contribution is less pronounced in the case of a fully embedded NP than in the case of conformational restrictions that flexible hydroxycarbon chains experience in the vicinity of a rigid transmembrane inclusion 24 .
The model assumes that there are no constraints on the distribution of phospholipids between NPs and planar biomembrane. This assumption corresponds to strong hydrophobic interaction between the NP and hydroxycarbon chains of membrane lipids. In the study of Rasch et al., 2010, such strong hydrophobic interaction is mediated by dodecanethiol coating of gold NPs In the case of weaker interaction, the larger NPs might induce membrane poration 33 , that is not considered within this study.
Our model also does not explicitly take into account the increase of the disorder of lipids in the vicinity of the NP 25 . This effect might be important in the lipid bilayer containing cholesterol, where the cholesterol content was shown to decrease near NP 25 . It was indicated that purely hydrophobic NPs accumulate at the interface between the cholesterol-rich, ordered domains and cholesterol-lean, liquid domains to reduce the net interfacial free energy 9 . The interfacial energy minimization was also proposed for protein self-assembly mediated by disordered domains around protein inclusions within ordered lipid bilayers 41 . The presence of cholesterol may also contribute to the depletion forces 24 between the NPs. These additional attractive forces may be responsible for NP clustering observed in the presence of cholesterol in the membrane 9 .
Our simulations were performed for phosphatidylcholine lipids with negative intrinsic curvature C 0 ( Table 1). The effect of variation of the intrinsic curvature on the membrane energy is minor as contributions of the positive and the negative curvature regions to elastic membrane energy are canceled (Fig. 2). Therefore, we may expect the existence of energy barrier between the condensed and the separated state of NPs also for membranes composed of the lipids with positive intrinsic curvature.
Within this study a membrane composed of single lipid sort was considered. The lipid mixture with various intrinsic curvatures 42 could influence NP-biomembrane interaction. We may expect also lipid sorting at the highly curved NP surface 43 . In this case, the Gaussian energy term should be included in equation (1) and entropic contribution [44][45][46] due to non-homogenous lateral distribution of lipids should be considered.
The presented study is based on two-dimensional analysis by taking the second principal curvature in the membrane to be zero. In three-dimensions, the non-zero second principal curvature should be taken into account. The Gaussian bending energy could also be considered constant in three-dimensions if the studied membrane area is the part of closed vesicle and the presence of NPs will not influence the membrane shape outside the studied region. Considering the second principal curvature in three dimensions would likely change the values of the elastic energy, but it would not influence the main conclusions about the role of the energy barrier between the NPs in biomembrane, preventing either their clustering or disassembling of already formed clusters.

Methods
The NP embedded within the lipid bilayer hydrophobic region induces perturbation in bilayer thickness and bends both bilayer leaflets (Fig. 2). The local shape of the lipid bilayer wrapping around the NP is governed by the interplay between the lipid tails' stretching energy G c and the bending energy of both monolayers G b . The local bending energy can be expressed as deviation of the membrane local curvature C from initrinsic membrane curvature C 0 47 , while the transversal stretching energy is considered as the deformation of the hydrophobic core 35,48 . Within this study, the problem of NP interaction is reduced to two-dimensions, as shown in Fig. 2. Two rigid NPs of diameter r separated by the distance d are considered. The average energy per lipid molecule can be expressed as where κ b and κ c are the bending and the compression-expansion constants of the lipid layers, respectively, ξ is the thickness of the deformed monolayer in the vicinity of the NP, ξ 0 is the length of the lipid tails of an unperturbed monoloayer, and a 0 is the area per lipid mocule. The integration is performed over the contour length l of the monolayers at the hydrophobic/water interface. The optimal shape is determined by the minimum of the average elastic energy per lipid in the presence of nanoparticle g with respect to its energy in the planar lipid bilayer g planar (Δg = g − g planar ). The energy per lipid molecule in planar bilayer can be expressed by taking C = 0 and ξ = ξ 0 in Eq. (1).
The configuration of the system corresponding to its minimal elastic energy was obtained by using the custom written optimization program employing basin hopping method 49 . The membrane contour is discretized and each step of the random perturbation is followed by the local minimization of the membrane elastic energy using the sequential quadratic programming 50 . The simulated annealing acceptance test based on standard Monte Carlo was used 51 . To reduce computational demands, one quarter of the geometry is solved considering symmetry about bilayer midplane and midplane between NPs (Fig. 2). The predictions of our model were verified by comparison with the results of an analytical model of Sub Wi et al. 10 for a single NP embedded in the membrane. The membrane element size was chosen on the basis of the mesh convergence test. The lipid parameters chosen for the simulation are shown in Table 1. Simulations are based on default values of parameters, if not stated otherwise. The range of the NP size is taken after Gopalakrishnan et al. 52 , who observed no NPs larger than 4 nm embedded in the lipid bilayer. To study the effect of membrane properties on the elastic energy of membrane with the NP inclusion, the variations of the intrinsic curvature C 0 , the bending κ b , and the compression-expansion constants κ c are performed. The range of the membrane elastic parameters adopted in simmulations is shown in Table 1.
Data availability. All data generated during this study are included in this manuscript. The code of the computational framework in Python is available from the corresponding author on reasonable request.


extfieldsofstudy:
Materials Science
Physics
Chemistry
Medicine

text:
Polarizable Molecular Simulations Reveal How Silicon-containing Functional Groups Govern the Desalination Mechanism in Nanoporous Graphene

We report a molecular dynamics (MD) simulation study of reverse osmosis desalination using nanoporous monolayer graphene passivated by SiH$_2$ and Si(OH)$_2$ functional groups. A highly accurate and detailed polarizable molecular mechanics force field model was developed for simulating graphene nanopores of various sizes and geometries. The simulated water fluxes and ion rejection percentages are explained using detailed atomistic mechanisms derived from analysis of the simulation trajectories. Our main findings are: (1) The Si(OH)$_2$ pores possess superior ion rejection rates due to selective electrostatic repulsion of Cl$^-$ ions, but Na$^+$ ions are attracted to the pore and block water transfer. (2) By contrast, the SiH$_2$ pores operate via a steric mechanism that excludes ions based on the size and flexibility of their hydration layers. (3) In absence of ions, water flux is directly proportional to the solvent accessible area within the pore; however, simulated fluxes are lower than those inferred from recent experimental work. We also provide some hypotheses that could resolve the differences between simulation and experiment.

Introduction
Seawater desalination plays a crucial role in the research for overcoming water scarcity, one of the most serious challenges humans are facing. 1,2 Multiple approaches have been investigated, including multi-stage flash, 3 membrane distillation, 3 mechanical vapor compression, 4 and reverse osmosis (RO). 5 Among them, reverse osmosis technology is by far the most promising method, as it requires the least amount of energy input to purify the same amount of water.
Conventional RO membranes made of composite materials, like polymers, has improved over decades to reach an acceptable energy efficiency; 6 this has led the growth of RO plants all over the world, but improvements in membrane permeability have been limited since the 1990s. 7 Therefore, people seek new forms of materials, such as carbon nanotubes, 8 and nanoporous graphene 9,10 that hold promise for desalination applications.
Single-layer nanoporous graphene is a promising desalination membrane due to its unmatched thinness, good mechanical strength, and chemical stability. Experimental studies on nanoporous graphene involve various treatments to create nanopores in graphene, including electrical pulse, 11 ion bombardment, 12 and O 2 plasma etching. 10 Ion passage through the prepared nanoporous graphene was observed to be highly selective in several studies, [11][12][13] with cation/anion selectivity ratio of over 100; ion selectivity is positively correlated with salt rejection. 11 Recently, one experimental study has shown nanoporous graphene can act as a desalination membrane with surprisingly high water flux and nearly perfect salt rejection. 10 The reported 10 6 g m −2 s −1 water flux (3000 water molecules per pore per ps) under only 0.17 bar driving pressure is several orders of magnitude higher than conventional membranes characterized by fluxes of about 12 g m −2 s −1 under 83 bar. 6 Scanning transmission electron microscopy (STEM) images suggested that the pores produced by O 2 plasma were around 1 nm in diameter. Lee et al. 14 proposed that graphene nanopores could be stabilized by silicon based on STEM imaging and density functional calculations; these results suggest that the edges of nanopores produced by oxygen plasma may also contain silicon. These studies have contributed essential insights into the nanopore structure; however, the molecular mechanism of water transport and ion rejection responsible for the high performance remain unknown, creating an opportunity for theory and simulation to provide detailed insights into the desalination process in nanoporous graphene.
Several published theoretical studies have focused on the possibility of nanoporous graphene as a desalination membrane. In 2012, Cohen-Tanugi and Grossman investigated the use of hydrogen-or hydroxyl-passivated graphene nanopores to separate NaCl ions from water. 9 Their molecular dynamics (MD) study suggested a maximum diameter of 5.5Å for maintaining a high percentage of salt rejection. They also followed up the studies with reduced graphene oxide 15 and multilayer porous graphene. 16 Strong and Eaves investigated the hydrodynamics of the water transferring through porous graphene 17 and suggested that the non-equilibrium dynamics of the water transfer is controlled by two competing mechanismstranslocation mechanism and evaporation-condensation mechanism. Ebrahimi 18 modeled the nanoporous graphene passivated by Si atoms and found that the membrane curvature reduces the water flux. The relation between pore ion selectivity and surface charge was studied by Zhao et al., 19 who predicted that negative charges at pore edges could impede the passage of Cl − . Very recently, Li et al. 20 found that the ion selectivity of multi-layer graphene nanopores may decrease if the surface charge becomes too strong due to counterions entering the pore and causing a charge inversion. These studies all used fixed-charge models that do not account for the environmental dependence of molecular dipole moments, raising the possibility that a more detailed physical model is needed to explain the water transfer phenomenon at the interface of the nanopore.
The inclusion of polarizability greatly improves on the physical descriptions of molecular mechanics (MM) force field models by explicitly modeling induced polarization from external electric fields or intermolecular interactions. 21 In particular, the AMOEBA model includes polarizable point dipoles and also incorporates atomic multipole moments up through the quadrupole, 22-25 providing a more detailed description of molecular electrostatics compared to other widely implemented polarizable force fields such as the Drude oscillator 26,27 and the fluctuating-charge models. 28,29 As a result, AMOEBA models can reach quantitative agreement against gas-phase ab initio interaction energies at various distances and orientations, and produce reliable predictions when the models are used to simulate condensed-phase phenomena. [30][31][32][33][34] Because molecular polarization depends on changes in the electrostatic environment, we expect a polarizable model to yield a more accurate description of water transfer and ion rejection as water molecules and ions move from the bulk and through the nanopore.
In this paper we describe an accurate "AMOEBA-type" polarizable force field for siliconpassivated nanoporous graphene, and its application in MD simulations on model nanopores that reveal their performance in desalination. We find that hydrophobic nanopores, passivated by SiH 2 groups (also denoted as Si-H), has superior water transfer rates compared to hydrophilic Si(OH) 2 nanopores (also denoted as Si-OH) ; on the other hand, hydrophilic Si-OH nanopores selectively block the negatively charged Cl − ions resulting in an improved salt rejection. We also provide new atomistic insights into the physical mechanisms that govern water transfer rates and salt rejection percentages on the atomic scale; one of our novel observations is that ions attracted by hydrophilic functional groups may block the pores, thereby significantly reducing water transfer rates. However, under the same assumed pore density as in the experimental work, our calculated transfer rates are still two orders of magnitude lower than experiment, 10 suggesting there still exist differences between the experimental result and our molecular picture.
The rest of the paper is organized as follows. First, we describe the development of our polarizable nanopore model and force field, the setup of desalination simulations, and the simulation analysis methodology. Next, we describe the optimized parameters and simulation results, including our simulated water transfer and ion rejection rates and discussions of the molecular mechanisms. Finally, we provide some hypotheses that may explain the remaining differences between simulation and experiment.

Polarizable Model
The general functional form of the AMOEBA polarizable force field can be divided into four terms: The first term describe valence interactions which include bond stretching, angle bending, and torsional rotations; they may also include out-of-plane bending terms and higher-order couplings. The last three terms describe the non-bonded interactions, including the van der Waals (vdW) force and the electrostatic contributions from permanent and induced multipoles. Two types of local molecular frames used for defining the permanent dipole and quadrupole moments are shown in the Supporting Information. ( Figure S1 and Table S6).
We constructed an atomistic model of nanoporous graphene passivated by Si atoms in the form of five-membered rings according to previous STEM imaging and density functional studies. 14 Figure 1 shows the molecular models used in QM calculations to develop of force field parameters for the Si-H pore. We constructed two medium-sized molecules, C 26 H 16 Si 2 ( Figure 1a) and C 24 H 14 Si 2 (Figure 1b), to characterize both the Si-C-C-C-Si and the Si-C-C-Si structures in the complete pores. We also created a small model incorporating only one Si atom (SI Figure S2), but it was not included in the final data set because it lacks the interaction between Si passivating groups. For carbon atoms, five different atom types, namely Ca1, Ca2, Cb, Cc, Cs, were created based on the topology and their distance to the pore center, as shown in different colors. Figure 1c illustrates how a complete Si 12 H 24 pore can be built by assembling the two types of medium-sized models. This large model contains 120 carbon, 12 silicon, and 60 hydrogen atoms.  Table 1 summarizes the reference data set for the Si-H pore model. In order to generate data describing the intramolecular forces, we carried out ab initio MD simulations at the B3LYP/6-31G* level of theory 35 using the TeraChem software package. 36,37 Under the simulated temperature of 1500 K, 300 snapshots were generated with an 0.5 fs time interval, from which single-point energies and gradients were computed at the DF-MP2/cc-pVTZ level of theory using the Psi4 software package. 38 We then characterized the non-bonded interactions by computing interaction energies at various model geometries including some The water-"large pore" model contains 195 atoms and 4567 basis functions with def2-TZVP basis set, which limited the number of calculations that we could afford to only four. We included them mainly as a sanity check to ensure the results did not deviate significantly from the medium models. The entire procedure was repeated for the Si-OH pore models, in which all of the SiH 2 terminating groups were replaced by Si(OH) 2 groups. Total 2247 A. DF-MP2/cc-pVTZ on geometries from B3LYP/6-31G* AIMD simulations. B. DLPNO-CCSD(T)/def2-TZVP on geometries optimized at B3LYP/6-311G* level of theory. The graphene carbon atoms more than four bonds away from Si atoms were described using one extra atom type, Ca0, which was not included in these models. The interaction parameters of Ca0 were fitted to coronene-water binding energies (Supporting Information Section S3). Although the modeling of graphene with finite polycyclic aromatic hydrocarbons is imperfect due to the closing of the HOMO-LUMO gap and emerging radical character at the graphene edges, 43 Lazar et al. have shown that for seven small organic molecules, experimental adsorption enthalpies on graphene can be well predicted by ab initio adsorption enthalpies on coronene. 44 To ensure transferability from the finite-sized coronene model systems to the continuous graphene model, we fixed the charge of Ha atom type (C-H) and both charge and dipole parameters for the Ca0 atom type to be zero.
The total charge of the large model was balanced to zero in our parameterization. In keeping with our assumption that the nanoporous graphene should be overall neutral in the simulations, the simulation setup of various pores included a uniform neutralizing charge to all "pore" atoms with a magnitude that never exceeded 0.0001 elementary charge.
Following the generation and organization of ab initio data, the force field parameters were optimized using ForceBalance. 45 An initial set of bonded parameters and vdW parameters were adopted from the carbon and hydrogen parameters of benzene and methane in the AMOEBA09 force field for organic molecules. 24 For the electrostatic interactions, we followed the protocol utilizing the GDMA program, 46 described in the supporting information of the AMOEBA09 paper 24 to generate the initial guess for the charge, dipole, and quadrupole moments. The AMOEBA03 water and ion model parameters were employed in our parameterization and simulation. 22,47   All desalination simulations used a Langevin thermostat with a collision frequency of 0.1 ps −1 to maintain a temperature of 348.15 K. The particle-mesh Ewald (PME) method 22 with cut-off distance at 0.9 nm was adopted when evaluating non-bonded interactions. The systems are first equilibrated for 5 ns with step length of 1 fs; a Monte Carlo barostat with 1 atm pressure was used to adjust the size of the periodic box once every 25 MD steps. All three dimensions of the periodic box were allowed to change individually to relax internal tensions within the graphene model. To prevent ions from leaking during this equilibration period, a virtual blocking force is added that occupies the space of the pore and only repels the ions (Supporting Information Section S5).

Simulation Setup
Following equilibration, the volume is held constant and an external force is applied to the graphene wall to simulate the driving pressure, pushing water from the feed side to the permeate side. Simulation data is collected for 10 ns and repeated ten times from different equilibrated geometries to obtain averages and standard error estimates, for a total of 100 ns per system. Each of the seven Si-H and seven Si-OH pores were simulated at four driving pressures (100, 200, 500 and 1000 bar) for a total of 5.6 µs of simulation data in the main study. We also performed several simulations with the graphene removed and/or ions removed to address specific questions, described in later sections.

Trajectory Analysis
To gain further insights on the water and ion transfer process, we developed a trajectoryanalyzing tool called the "buffered-region method". Our analysis method first estimates the z-coordinate and thickness of the nanoporous graphene sheet. The z-coordinate is computed as the averaged geometric center of C atoms, and the thickness is estimated as the distance where the radial distribution function (RDF) between pore and water atoms reaches the chosen threshold of 0.3. As shown in Figure 5a, the membrane thickness turns out to be 0.28 nm from center to surface, though this is rather insensitive to the chosen RDF threshold.
The simulation box is divided into three sections: feed, middle, and permeate, and the events for each individual water molecule crossing the boundary of two sections are recorded. As illustrated in Figure 5b, a forward transfer event is counted when a water molecule crosses from the feed to middle region, followed by a crossing from the middle to permeate region; the reverse ordering is used to count backward transfer events. Any other crossing sequences, such as feed to middle then back to feed, are ignored. The total exchange rate is defined as the sum of forward and backward transfer rates, while the net transfer rate is the difference between forward and backward transfer rates. The analysis tool also allows us to access the characteristics of individual transfer events including transfer times and changes in the hydrogen bonding environment. Our analysis of the water transfer and ion rejection mechanism necessitate two complementary definitions of pore area. Based on the total water exchange, we calculate the effective pore area as: where A pore;eff is proportional to the total water exchange N pore , and A 0 /N 0 is a normalization factor calculated from water exchange across an unimpeded cross section of the bulk solution.
We also define the accessible pore area by integrating the region inside the pore with a nonzero water density; (Figure 6) in contrast to the former definition, the accessible pore area does not depend on dynamical properties such as the total exchange rate. In practice, the accessible pore area is converged to within 1% after collecting 100 ns of the simulation data. All molecular dynamics simulations were performed using the OpenMM package. 48 The three-dimensional rendering were generated using VMD. 49 Trajectory analyses were carried out in Python using the MDTraj library. 50 3 Results and Discussion

Optimized Force Field Parameters
The final force field with optimized parameters accurately reproduces the QM interaction energy data. For example, Figure 7 shows the SiH-Cl The partial charges and dipole moments of the optimized pore models are shown in

Water Transfer
The net water flux as a function of pore size and driving pressure is plotted in Figure 9a; an increase in net transfer proportional to applied pressure is observed as expected. The effect of functional groups is significant, as the flux for each size of the Si-OH pore is at most half of the corresponding Si-H pore.

Pore Area Comparison
SiH, Accessible SiH, Effective SiOH, Accessible SiOH, Effective SiOH, No Ions, Effective Figure 10: Comparison between accessible pore areas and effective pore areas; the latter is computed using total water exchange under 100 bar driving pressure.
The first hypothesis is tested directly by computing the accessible and effective pore areas. Figure 6 illustartes how the accessible pore area is computed for the SiH 2 -9 pore; the water density pattern is hexagonal, indicating that three SiH 2 groups interact less closely with water than the other six. The effective pore area is computed from the total water exchange as in Eq. 2, and the comparison of two area measurements is shown in Figure 10.
Despite the fact that water density is not equally distributed inside the pores, we found that the accessible pore area for each Si-H pore matches very well with the effective pore area. This agreement indicates that pore permeability is proportional to the accessible pore area, consistent with the macroscopic description of maximum flow rate in the limit of short pipe lengths predicted by Bernoulli's principle, Q max = A 2∆P/ρ. 51 For Si-OH pores, the accessible pore areas are smaller than Si-H by around 5-10Å 2 , but are not enough to rationalize their much smaller effective pore areas; this indicates the larger size of OH functional groups (Hypothesis #1) cannot fully account for the observed difference between the pores.
After carefully inspecting the simulated trajectories, we found that the Si(OH) 2 -9 pore is occupied by three or more Na + ions. (SI section 6) These ions become immobile as they balance the negative partial charges from oxygen atoms on the Si-OH edge groups, and as a result the water transfer through the pore is greatly slowed down. The average transfer time is 100 ps through the Si(OH) 2 -9 pore, compared to 10 ps through the SiH 2 -9 pore. To confirm the effect of ion blockage, we ran a separate set of simulations with the ions removed.
Without the ions blocking the pore, the water flow rates are much higher, and the resulting effective pore area for Si-OH closely matches the accessible pore area. This rules out the kinetic barrier hypothesis (#2) and confirms the ion blockage hypothesis (#3). The interaction between water and pore edges can be further understood by quantifying how the hydrogen bonding network changes as a water molecule passes through the pore. Figure 11 shows the number of hydrogen bonds (N hyd ) each water molecule is involved in, calculated using the hydrogen bond criterion of Baker and Hubbard. 52 We computed a reference value of < N hyd >= 3.6 for the average number of hydrogen bonds in bulk water simulated with the AMOEBA force field at 348 K; the KDE shows a clear peak at n=4 reflecting the underlying histogram. When the analysis is narrowed down to only the molecules undergoing transfer events, the distribution of N hyd is shifted towards smaller numbers. For the smallest and largest Si-H pore, < N hyd >= 1.7 and 3.2 respectively; in particular < N hyd >= 2.88 for SiH 2 -9. The reduced hydrogen bonding indicates that water molecules are partially desolvated as they are transferred through the pore, with < N hyd > decreasing by 0.5-2 depending on pore size. For the Si(OH) 2 -9 pore, the OH group pore edges can serve both as hydrogen bond donor and acceptor, and we expected the transferring water to make more hydrogen bonds than in the corresponding SiH 2 -9 pore; however, the distribution is highly similar to the SiH 2 -9 pore with < N hyd > = 2.83. The reduction in hydrogen bonding is again caused by the Na + ions filling the pore; when ions were removed from the system, the water molecules transferring through the Si-OH pores have < N hyd > = 3.76, even a little higher than bulk water. Additionally, we noticed that the population for zero hydrogen bonds is negligible, suggesting that the evaporation-condensation mechanism 17 is not playing a significant role in these simulations.
We investigated whether water transfer events are correlated in time by examining the time windows in which a water molecule occupies the middle region for all transfer events.
For the SiH 2 -9 pore we found that 70% of water transfer events overlap with other water transfer events in time. However, 70% of the entire trajectory contains at least one water transfer event, which means if the water transfer events are distributed randomly with no correlation, there is still a 70% chance that any one of them would overlap with another.
Therefore, we conclude that existing water transfer events did not have a significant effect on the possibility of the occurrence of another event.  Figure 12: Comparison of the Na + and Cl − rejection rates for Si-H and Si-OH pores, under 100 bar driving pressure (left) and 500 bar driving pressure (right). Negative rejection rates were rectified to zero. Extremely small number of water transfer prevents us from calculating a statistically meaningful rejection rate for Si(OH) 2 -6 pores.

Ion Rejection
The ion rejection rates are calculated by comparing the ion concentration of transferred water molecules and ions relative to the initial feed solution as: where R = 100% corresponds to zero ion leakage, and R = 0% means the transferred molecules have the same or higher ion concentration than in the feed; only transfer events in the forward direction are counted. We find significant differences in ion rejection between positively and negatively charged ions and a strong dependence on the functional groups on the pore edges. The Si-H pores tend to reject Na + ions with a higher percentage than the Cl − ions, as shown in Figure 12; by contrast, the Si-OH pores selectively blocks the Cl − ions and does not block Na + ions at all. These differences can be explained by two distinct mechanisms of ion rejection and the characteristics of Na + and Cl − hydration layers.  Figure S135) shows that a broad peak at low values of the acute angle vanishes for ions inside the pore, indicating that water oxygen atoms in the Cl − hydration layer can undergo significant conformational changes. Figure 13 shows the solvation structure of ions found inside a SiH-9 nanopore and provides intuition for how a more flexible hydration layer around Cl − enables it to permeate the pore more effectively.
(a) (b) Figure 14: Two snapshots of a Na + ion blocking a SiOH-9 nanopore. Left: The Na + hydration layer easily makes hydrogen bonds with the SiOH groups. Right: The SiOH groups can also replace water molecules in the Na + hydration layer.
By contrast, the Si-OH pores repel ions primarily by electrostatic interactions (also called the charge rejection mechanism); the OH groups on the pore edge create a negatively electrostatic potential that attracts cations while exclusively blocking anions. Figure 14 shows two snapshots of a Na + ion occupying the interior of a Si(OH) 2 -9 nanopore; Na + may retain its first hydration shell (left panel), or Si-OH groups may take the place of water molecules and coordinate to Na + directly (right panel). This is further substantiated by a large drop in the Na + -O coordination number from 5.7 in the bulk to 4.0 inside the Si-OH pore, which does not occur for Si-H pores (SI Table S16). The ability of Si-OH groups to substitute for water in the first hydration shell implies that the accessible pore area for Na + inside Si-OH pores is greater than in Si-H pores, and explains how Na + may pass through even the smallest Si(OH) 2 -6 pore with a diameter less than 6Å. By contrast, when a Cl − ion permeates the Si-OH pore (which occurs only rarely), we find there are no strong interactions with the pore edges; a snapshot of this system is shown in SI Figure S136.
Because the permeate solution cannot develop a macroscopic net charge, blocking one type of ion will effectively prevent the large-scale leakage of counter-ions as well; thus we should use the Na + rejection rates for the Si-H pores, and Cl − rejection rates for the Si-OH pores to characterize their overall salt rejection performance. The ion rejection rates for each pore size under different driving pressures are compared in the left and right panel of Figure 12; the ion rejection rates of Si-H pores are nearly 100% for small pores and starts to decrease for SiH 2 -9 and larger pores, in agreement with the size rejection mechanism.
On the other hand, Si-OH pores have nearly 100% ion rejection rates even for the relatively large Si(OH) 2 -11 whose diameter exceeds the hydration shell of Cl − . These results differ from previous studies 9,58 in which the authors concluded a lower salt rejection for hydroxylated (C-OH) pores than hydrogenated (C-H) pores; the discrepancy could be attributed to differences in the pore structures used in the simulations, levels of detail in the force fields, and/or simulation analysis approaches. Ion rejection also decreases with an increase of driving pressure; if we try to maximize the water flux while holding ion rejection at a constant value, we find that smaller pores with higher driving pressures have superior performance to larger pores with lower driving pressure. For example, water flux through SiH 2 -9 under 1000 bar (83 water molecules per ns) is higher than SiH 2 -11 under 200 bar (39 water molecules per ns), yet both have ion rejection rates above 90%.

Discussions
In our simulations, the largest Si-H and Si-OH pores with close to 100% ion rejection rates are SiH 2 -9 and Si(OH) 2 -11. Although Si(OH) 2 -11 has a larger accessible pore area, its average net water transfer rate (8.5 water molecules per ns) is smaller than SiH 2 -9 (12.2 water molecules per ns) under 200 bar of driving pressure; the second value is consistent with a previous study that estimates 10.9 water molecules per ns through C-H pore with 5.5 A diameter. 59 Although the Si-OH pores suffer from cations entering and blocking the pore, the potential transfer rate can be as high as 32.9 water molecules per ns if the ion blockage could be avoided. In experiments, ion blockage was suspected as the main reason causing the degrading of synthesized nanoporous graphene, but smaller ion concentrations were used (6 mM/mol) and the performance degradation happens on the time scale of 24 hours. Another possible route to reduce ion blockage is to reduce the negative charge density by replacing a fraction of OH groups with H or other less charged functional groups; ideally a good balance could be achieved that does not diminish the good Cl − rejection performance.
In the experimental study of Ref., 10 the authors concluded on a massive 3000 water molecules ns −1 pore −1 , more than two orders of magnitude higher than our simulated results.
The large difference might be explained by the large area of nanoporous graphene in direct contact with water. In the experimental setup, the nanoporous graphene prepared by oxygen plasma was suspended over a SiN wafer, in which a 5 µm diameter hole was made for water to pass through, and the hole area was considered as the total nanoporous graphene area that was operational in the desalination experiment. However, the area of the nanoporous graphene in direct contact with water is hundreds of times larger than the area of the 5 µm diameter hole in the SiN wafer. Therefore, if water could be transported laterally between the graphene and SiN wafer, the nanoporous graphene area responsible for the water transfer may be larger than the hole area. This hypothesis can be tested by expanding the 5 µm hole to 10 µm and even larger, and see if the water transfer rate approaches a limit instead of increasing proportionally to the hole area.
The same paper reports a separate experiment where the authors calculate the flux of forward osmosis (water entering the feed side) to be 70 g m −2 s −1 atm −1 . 10 Based on this number and assuming a pore density of 1 per 100 nm 2 , we calculated the initial water transfer rate under 50 atm of osmotic pressure to be 10 molecules per ns per pore. This is in close agreement with our simulated forward osmosis transfer rate with 1 mol/L NaCl solution and no driving pressure, which is 6 water molecules per ns per pore for the SiH-9 pore. The level of agreement between forward osmosis simulations and experiment suggests that the microscopic description of the pore and ion/water transfer is reasonable, and the remarkable desalination performance may originate from macroscopic features of the experimental configuration.
Another possibility for explaining the higher water transfer numbers in experiment is that our optimal pore size could be too small. We estimated that the optimal pores for high flux and salt rejection should have a diameter of about 6Å, which is consistent with earlier theoretical work 9 using non-polarizable force fields that predicted high flux and salt rejection with pores approximately 5.5Å in diameter). It is possible that long-range charge transfer in the graphene layer could allow nanopores to accumulate significant amounts of net charge, thereby allowing more ions to be rejected via the charge rejection mechanism (though Ref. 20 predicts an inverse effect when charges are too high).

Conclusions
In the present work, we investigated the molecular mechanism of water desalination using nanoporous graphene passivated by Si-H and Si-OH groups. Our simulation used a newly developed AMOEBA-type polarizable force field that reproduces high-level ab initio quantum chemistry data with good accuracy. Reverse osmosis simulations were carried out with various pore sizes and driving pressures to measure net water transfer and ion rejection rates; the effective pore area of Si-H pores based on water transfer rates is commensurate with the solvent-accessible pore area, while for Si-OH pores the water transfer is impeded due to Na + ions occupying the pore interior. Size-based and charge-based ion rejection mechanisms were discussed for the Si-H and Si-OH types of pores respectively. Our simulations and analyses revealed complex relationships between water, ions, and pore edges at the nanoporous graphene interface, suggesting that intermediate values of pore-ion interaction strengths may provide a route toward optimizing this system. While our models provide an accurate and physically realistic description of intermolecular interactions between pore edges and water/ions, they do not account for long-range charge transfer and local deformations of the graphene monolayer that may also impact desalination performance. Inclusion of these additional effects would further improve on the modeling of nanoporous graphene and further advance our understanding of this promising system.