Coarse-grained molecular dynamics simulation for uptake of nanoparticles into a charged lipid vesicle dominated by electrostatic interactions
Naofumi Shimokawa Hiroaki Ito Yuji Higuchi


 


nshimo@jaist.ac.jp

School of Materials Science, Japan Advanced Institute of Science and Technology, Ishikawa 923-1292, Japan


[Present address: ]Department of Physics, Chiba University, Chiba 263-8522, Japan


Department of Mechanical Engineering, Graduate School of Engineering, Osaka University, Osaka 565-0871, Japan




Institute for Solid State Physics, University of Tokyo, Chiba, 227-8581, Japan





We use a coarse-grained molecular dynamics simulation to investigate the interaction between neutral or charged nanoparticles (NPs) and a vesicle consisting of neutral and negatively charged lipids.
We focus on the interaction strengths of hydrophilic and hydrophobic attraction and electrostatic interactions between a lipid molecule and an NP.
A neutral NP passes through the lipid membrane when the hydrophobic interaction is sufficiently strong.
As the valence of the positively charged NP increases, the membrane permeation speed of the NP is increased compared with the neutral NP and charged lipids are accumulated around the charged NP.
A charged NP with a high valence passes through the lipid membrane via a transient channel formed by charged lipids or transportlike endocytosis.
These permeation processes can be classified based on analyses of the density correlation function.
When the non-electrostatic interaction parameters are large enough, a negatively charged NP can be adsorbed on the membrane and a neutral lipid-rich region is formed directly below the NP.
The NP is spontaneously incorporated into the vesicle under various conditions and the incorporation is mediated by the membrane curvature.
We reveal how the NP's behavior depends on the NP valence, size, and the non-electrostatic interaction parameters.




                              
Coarse-grained molecular dynamics simulation 
 for uptake of nanoparticles
 
into a charged lipid vesicle dominated by electrostatic interactions
    Yuji Higuchi
    December 30, 2023
=================================================================================================================================================








§ INTRODUCTION


Phospholipids, which are typical amphiphilic molecules, self-assemble into a bilayer structure in the aqueous phase so as to protect their hydrophobic parts from the aqueous phase.
The bilayer structure is essentially the same as the structure of the plasma membrane, and thus phospholipid vesicles, which are spherical closed bilayers, are widely used as a simple model for cells.
Various physicochemical phenomena in lipid vesicles are deeply related to biocellular functions <cit.>.
In multicomponent lipid membranes, a compositional heterogeneous structure emerges spontaneously at lower temperatures <cit.>.
This phase separation has attracted great attention in relation to raft domains in biomembranes <cit.>.
Moreover, dynamic membrane deformation, such as endocytosis and exocytosis, can be reproduced in some vesicular systems <cit.>.


Biomembranes also contain negatively charged lipids in addition to zwitterionic lipids.
Negatively charged lipids, such as phosphatidylserine, phosphatidylglycerol, phosphatidylinositol, phosphatidic acid, and cardiolipin, play important roles in protein adsorption on membranes <cit.> and channel activity <cit.>.
In addition, several studies reported that negatively charged lipids greatly affect phase separation and membrane deformation in artificial lipid membranes.
For example, the suppression of the phase separation in unsaturated charged lipid-containing membranes was examined, and the addition of salt made the phase behavior similar to that of neutral membranes <cit.>.
This behavior was explained by theoretical and numerical studies <cit.>.
In addition, decreasing the temperature spontaneously changed the vesicle morphology from spherical to disk shape <cit.>.
Thus, electric charges on lipids play important roles in membrane stability, which is a crucial factor in the uptake of proteins and drugs into cells.


The interaction between nanoparticles (NPs) and the plasma membrane is important in NP cytotoxicity in the development of nanotechnology and drug delivery <cit.>.
Moreover, this type of interaction is also relevant to the uptake of viruses into cells <cit.>.
Many studies using living cells have been performed to reveal how NPs and nanomaterials are incorporated into cells and the cytotoxicity caused by NP incorporation.
Various materials, such as fullerene <cit.>, carbon nanotubes <cit.>, dendrimers <cit.>, gold NPs <cit.>, cobalt oxide particles <cit.>, and quantum dots <cit.>, have been investigated.
Studies using vesicles focus on the interaction between NPs and lipid membranes and the uptake behavior of NPs into vesicles <cit.>.
Changes in membrane surface properties <cit.>, membrane protrusion <cit.>, and emergence of lipid accumulation <cit.> have been observed by NP and protein adsorption on lipid membranes.
Furthermore, the uptake behavior of NPs into vesicles has also been studied; the membrane permeation and endocytic behavior of NPs or nanoemulsions <cit.> and the effects of phase-separated domains on NP permeation <cit.> have been examined experimentally.


One of the NP uptake behaviors is a passive process.
Passive processes are caused by the interaction between the NP and membrane without an external energy supply.
To understand NP uptake better, it is important to understand the passive process mechanism, which is the basic uptake mechanism  <cit.>.
Because it is difficult to observe the dynamic behavior of NP uptake on the nanometer scale experimentally, numerical simulation studies have been performed.
The hydrophobic interaction between the lipid membrane and the hydrophobic part of NPs, the specific interaction between the lipid membrane and ligands grafted onto the surface of NPs, and the electrostatic interactions between the polar groups on the lipid headgroups and charged NPs all play important roles in NP uptake.
In the studies of the hydrophobic interactions and ligand interactions, the size, shape, and surface properties of the NPs <cit.> and the cooperative behavior of multiple NPs <cit.> have been discussed.
In particular, the receptor-mediated endocytosis of NPs can be described by the competition between the curvature elastic energy of lipid membrane and the attractive energy between lipid membrane and NPs via ligand-receptor binding <cit.>.
The endocytic motion of larger NPs occurs easily, since the energy loss by curvature change becomes smaller and the energy gain by the attractive energy becomes larger.
However, the model is not sufficient because the lipid membrane is consisting of single component and the long-range interaction such as electrostatic interaction is not taken into account.


Although calculations that include the electrostatic interactions usually become complicated due to the long-range nature of electrostatic interactions, several studies have simulated NP uptake mediated by electrostatic interactions <cit.>.
In particular, the interaction between charged NPs and a mixed membrane consisting of neutral and charged lipids was investigated to mimic the cellular environment <cit.>.
Lipids with charge opposite to NPs are accumulated in the NP adsorption part <cit.>.
As a result, the large adsorption energy derived from electrostatic interaction is obtained, which overcomes the curvature energy loss, and endocytic movement occurs.
This picture is essentially the same as the receptor-mediated endocytosis. 
However, there are still few of these studies and most studies use a box-spanning flat membrane, where the asymmetry between the inner and outer space is not incorporated, and the vesicular environment with electrostatic interactions is usually not considered.
Furthermore, it is necessary to investigate how NP incorporation into a vesicle depends on interaction and NP size because the types of NPs and drugs are diverse.


We have previously conducted a coarse-grained molecular dynamics (CGMD) simulation to examine the dynamic behavior of binary lipid membranes consisting of neutral and charged lipids <cit.>.
In the present study, we investigate the dynamic behavior of an NP placed near the outer leaflet of a charged lipid membrane based on the previous CGMD simulation.
CGMD is suitable for describing the uptake of NPs into a vesicle while changing the interaction strength.
We systematically change the strengths of hydrophilic, hydrophobic, and electrostatic interactions between a lipid molecule and an NP.
The details of the CGMD simulation are introduced in Sec. <ref>. 
The interaction between a neutral NP and a charged membrane and the dynamic behavior of a neutral NP are presented in Sec. <ref>.
The interaction and the dynamic behavior of a charged NP are discussed in Sec. <ref>.
Finally, we discuss the implications of our work and compare it with other studies in Sec. <ref>.







§ METHODS



Based on our previous CGMD simulation <cit.>, which was based on the CGMD simulation by Cooke et al. <cit.>, we calculate the behavior of charged lipid membranes.
This model describes the phase separation and spontaneous membrane morphological change of charged lipid membranes well <cit.>.
In the present study, we calculate the interaction between a charged lipid vesicle modeled by our previous study and an NP.



A lipid molecule consists of one hydrophilic head bead and two hydrophobic tail beads. 
For convenience, we refer to the hydrophilic bead, the hydrophobic bead in the center of lipid molecule, and the hydrophobic bead at the end of lipid molecule as a-bead, b-bead, and c-bead, respectively.
These beads are connected linearly through springs.
The lipid-lipid and NP-lipid interactions are indicated by superscript LL and NL, respectively.
The excluded-volume interaction between two beads separated by distance r is


    V^ (LL)_ re(r;b)=
                      4v[ ( b/r)^12 - ( b/r)^6 + 1/4],     r ≤ r_ c, 
    
                      0,     r>r_ c,


where r_ c=2^1/6b. 
v=k_ BT is the unit of energy, where k_ B is the Boltzmann constant and T is absolute temperature.
To form a stable bilayer, b_ head,head=b_ head,tail=0.95 σ and b_ tail,tail=σ are chosen <cit.>. 
Here, σ is the unit of length corresponding to the cross-sectional diameter of a single lipid molecule.
The stretching and bending potentials of bonds between connected beads are expressed as


    V_ bond(r)=1/2k_ bond(r-σ)^2


and


    V_ bend(θ)=1/2k_ bend(1-cosθ)^2,


where k_ bond=500v is the bonding strength of connected beads, k_ bend=60v is the bending stiffness of a lipid molecule, and θ is the angle between adjacent bond vectors.
The hydrophobic attractive interaction among hydrophobic beads is expressed as


    V^ (LL)_ at(r)=
                    -v,     r<r_ c, 
    
                    -v cos^2[π(r-r_ c)/2 w_ c],     r_ c≤ r ≤ r_ c+w_ c, 
    
                    0,     r>r_ c+w_ c,


where w_ c is the cut-off length for the attractive potential. 
The electrostatic repulsion between charged hydrophilic beads is described as the Debye-Hückel potential


    V^ (LL)_ el(r)=v ℓ_ B z_1z_2exp(-r/ℓ_ D)/r,


where ℓ_ B=σ is the Bjerrum length, z_1 and z_2 are the valencies of the interacting charged headgroups,
and ℓ_ D=σ√(ϵ k_ BT/n_0e^2) is the Debye–Hückel screening length, which is related to the bulk salt concentration, n_0, where ϵ and e are the dielectric constant of the solution and elementary charge, respectively. 
As a salt, we consider a symmetric monovalent salt, such as NaCl.
Because typical charged headgroups, such as phosphatidylglycerol, phosphatidylserine, phosphatidylinositol, and phosphatidic acid, have a monovalent charge, we set z_1=z_2=-1. 
We do not set a cutoff for this screened electrostatic interaction.
The modeling so far is the same as our previous study <cit.>.


Next, the interaction between an NP with diameter d and lipids is introduced.
The potential between an NP and lipid beads is expressed based on the Lennard–Jones potential,


    V^ (NL)_ LJ(r)=4v_ b[ ( b/r)^12 - ( b/r)^6].


We set v_ b=v_ nh for the hydrophilic beads and  v_ b=v_ nt for the hydrophobic beads.
The electrostatic interaction between a charged NP and charged beads is written as


    V^ (NL)_ el(r)=vℓ_ B z_1 z_ n/[1+(σ/2ℓ_ D)][1+(d/2ℓ_ D)]
    exp{-[r-(σ+d)/2]/ℓ_ D}/r,


where z_ n is the NP valence.
We consider that the electric charges are distributed on the surface of NP in this potential <cit.>.
Our simulation implicitly considers the same salt solution for exterior and interior of the vesicle and the same salt solution is also present inside the membrane in this modeling of electrostatic interaction.
It should be noted that we do not consider the presence of water molecules explicitly.
Instead, it is included as electrostatic interaction implicitly.
Therefore, in this calculation, volume change of water phase inside the vesicle and osmotic pressure are not taken into account.


Each bead and an NP obey the stochastic dynamics described by the Langevin equation


    md^2r_i/dt^2=-ηdr_i/dt+f^V_i+ξ_i,


where m=1 and η=1 are the mass and drag coefficients, respectively.
Forces f^V are obtained from the derivatives of the interaction potentials Eqs. (<ref>)–(<ref>) for the beads and Eqs. (<ref>) and (<ref>) for the NPs. 
Constant τ=ησ^2/v is chosen as the unit for the time scale, and the time increment is set at dt=7.5 × 10^-3τ which is sufficiently small to keep the stability of the calculation.
The Brownian force, ξ_i, satisfies the fluctuation-dissipation theorem


    < ξ_i(t)ξ_j(t')>=6vηδ_ijδ(t-t').




The calculated lipid membrane is a binary mixture of A-lipids with a monovalent electric charge at the head beads, corresponding to charged lipids, and neutral lipid (B-lipid)s. 
Charged lipid (A-lipid) is represented by red hydrophilic bead and green hydrophobic beads, while neutral lipid (B-lipid) is represented by yellow hydrophilic bead and cyan hydrophobic beads.
A lipid vesicle is composed of 2500 charged lipids (A-lipids) and 2500 neutral lipids (B-lipids), and the salt concentration is 100 mM.
The fraction of charged lipids in plasma membranes is roughly 20–30%.
On the other hand, various compositions have been used in the experiments with artificial lipid membranes to investigate the effect of electrostatic interaction <cit.>.
Hence, we calculate 50% charged lipid-containing membranes to see the effect of the electrostatic interaction clearly.
We constructed the vesicles with setting lipids randomly by hand as spheres as shown in Fig. S1(a) <cit.> and then started the simulations <cit.>.
Before we calculate the interaction between an NP and the lipid membrane, we perform the calculation until t=10 × 7500τ to relax the lipid membrane structure sufficiently <cit.>.
The cutoff lengths for the attractive interaction among hydrophobic beads are set to w_ c^ AA/σ=1.8, w_ c^ BB/σ=1.7, and w_ c^ AB/σ=1.575.
When the phase separation takes place at w_ c^ AB/σ=1.55, the attractive energy among hydrophobic beads V^ (LL)_ at gradually decreases and reaches the plateau after t=2 × 7500 τ as shown in Fig. S1(b) <cit.>.
Therefore, the membrane structure at t=10 × 7500 τ can be regarded as a sufficiently relaxed state, even if the phase separation does not occur at w_ c^ AB/σ=1.575.
More details of the relaxation process are shown in Fig. S1 <cit.>.
As the initial state of the membrane structure, we choose a sufficiently relaxed membrane at t=10 × 7500τ.
After an NP is put on the outer leaflet of a vesicle with zero velocity (we reset this time as t=0 × 7500τ), we start the calculation for t=1 × 7500τ.
As for NP size, it has been reported that NPs less than 3.1 nm can easily pass through a lipid membrane <cit.>.
Therefore, we use two sizes of NPs, d/σ=2 and 5 which roughly correspond to 1.4 nm and 3.5 nm, respectively.
We use v_ nh/k_ BT=1, 2, and 3, v_ nt/k_ BT=1, 2, and 3, and z_ n=0, +10, +50, +150, +200, and -50.
v_ nh and v_ nt are related to the phase equilibrium between adsorbed and dispersed states, where the solubility is roughly estimated as exp(-v^*/k_ BT) for v^* = v_ nh or v_ nt.
The all parameters are summarized in Table S1 <cit.>.
Calculations are performed at least three times for each condition to ensure data reproducibility.
We show a representative result in the calculations we performed.
The results shown are almost identical to the other calculation results which are shown in Figs. S6, S7, S8, S9, and S10 <cit.>.







§ RESULTS







 §.§ Neutral NP



In this section, we investigate the interaction between a neutral NP (z_ n=0) and a charged vesicle.
The non-electrostatic interaction parameters of the Lennard–Jones potential V^ (NL)_ LJ, v_ nh/k_ BT and v_ nt/k_ BT, are changed from
1 to 3, and the dynamic behavior of an NP with the diameter d/σ=5 or 2 is examined.
To ensure data reproducibility, the results obtained by the other two calculations are shown in Fig. S6 <cit.>.







Figure <ref> shows typical dynamic behavior of an NP for d/σ=5, v_ nh/k_ BT=2, and v_ nt/k_ BT=3, and the NP passes through the lipid membrane (Fig. <ref>(a)).
First, the NP is strongly adsorbed on the outer leaflet (snapshot at t=0.04 × 7500 τ in Fig. <ref>(a)), passes through the hydrophobic region of the lipid membrane  (snapshot at t=0.07 × 7500 τ in Fig. <ref>(a)), and is finally adsorbed on the inner leaflet (snapshots after t=0.18 × 7500 τ in Fig. <ref>(a)).
The Lennard–Jones potential between the NP and lipids, V^ (NL)_ LJ, is shown in Fig. <ref>(b).
Figure <ref>(c) shows the enlarged view from t=0 to 0.1× 7500τ in Fig. <ref>(b).
The potential V^ (NL)_ LJ is immediately decreased and becomes almost constant until t ∼ 0.05 × 7500 τ.
This step indicates the adsorption of the NP on the lipid membrane.
The sudden energy decrease by the NP adsorption is found in all three calculations as shown in Fig. S6 <cit.>.
Just after t ∼ 0.05 × 7500 τ, the potential V^ (NL)_ LJ is dramatically decreased again.
This process is the permeation of the NP through the lipid membrane.
This energy drop is also found in other two calculations at t=0.02 × 7500 τ and t=0.07 × 7500 τ (Fig. S6 <cit.>) which are almost the same time scale as t=0.05 × 7500 τ.
For NP passage through the bilayer membrane, lipids are pushed away by NP.
The displaced lipids return to fill the hole at t=0.18 × 7500 τ, in other words, the defect made by NP passage is completely fixed. 
The NP that has passed through the membrane is adsorbed on the inner leaflet.
The potential V^ (NL)_ LJ gradually decreases from t=0.18 × 7500 τ because the vesicle shape elongates to fit the NP (Fig. <ref>(a)).
The electrostatic interaction between charged lipids, V^ (LL)_ el, is not changed significantly (Fig. <ref>(d)).
Therefore, the lipid distribution in the membrane is not changed by the adsorption or the permeation of the NP.


We note that the time scale of permeation (t ≃ 0.1 × 7500 τ) is much shorter than that of membrane relaxation (t ≃ 2 × 7500 τ) as shown in Fig. S1 (b) <cit.>.
From Fig. S1 <cit.>, the attractive energy V^ (LL)_ at is decreased from -1.31 × 10^5 k_ BT to -1.33 × 10^5 k_ BT and the decrease is 2 × 10^3 k_ BT.
The energy difference per lipid is 2 × 10^3 k_ BT / 5000 = 0.4 k_ BT, which is smaller than v_ nh and v_ nt.
Therefore, during the permeation process, the interaction between NP and lipid is predominant as compared to lipid-lipid interaction, and the permeation time scale is almost determined only by the NP-lipid interaction. 


We focus on the process of the NP permeation through the lipid membrane.
At t=0.07 × 7500 τ, the NP passes through the lipid membrane.
In Fig. <ref>(e), the density correlation functions of A- and B-lipids around the NP are indicated by red and blue lines, respectively.
When the distance from the NP is defined as r_ n, these lines start from r_ n/σ=2.5 located on the surface of NP and go up to 0.6 immediately.
The values of density correlation functions keep around 0.6.
Therefore, these lines change in the same way, indicating that the A- and B-lipids are uniformly distributed around the NP.
Figure <ref>(f) shows the density correlation functions of hydrophilic a-beads (red line), hydrophobic b-beads located at the center of the lipid molecules (blue line), and hydrophobic c-beads located at the end of the lipid molecules (green line) around the NP.
There are more hydrophilic beads near the NP and we can see the peak value 0.7 at r_ n/σ=3.5.


Finally, the NP is adsorbed on the inner leaflet and induces the slight elongation of the vesicle.
The NP adsorbed to the inner leaflet tends to deform the lipid membrane and come into contact with the membrane of a larger area.
As a result, the vesicle is elongated.
This is caused by competition between adsorption energy and curvature energy to bend the lipid membrane.
The angle which determines the adsorption area is shown in Fig. S11 <cit.>, and it is 213.5 ± 1.0^∘ for z_ n=0, v_ nh/k_ BT=2, v_ nt/k_ BT=3.
In addition, the density correlation functions of A- and B-lipids in the inner and outer leaflets around the NP are shown in Figs. <ref>(g) and (h), respectively.
In the inner leaflet, there is a slight increase in neutral lipids (B-lipids) near the NP.
The inner leaflet near the NP has a large negative curvature to adsorb to the lipid membrane, and the charged headgroups are close to each other.
Consequently, it is difficult for charged lipids (A-lipids) to exist in this region due to the large electrostatic repulsion between charged lipids (A-lipids).
On the other hand, in the outer leaflet, there are some charged lipids (A-lipids) near the NP.
Because the outer leaflet has a positive curvature, the electrostatic repulsion is decreased, allowing charged lipids (A-lipids) to localize in this region.
However, the effect of lipid accumulation induced by the vesicle shape deformation is limited.


The NP is placed near the outer leaflet and spontaneously enters the interior of the vesicle.
Because the vesicle size is small, the NP on the inner leaflet can interact with many lipid beads without large membrane deformation; the curvature of NP fits the inner membrane easily.
The Lennard–Jones potential in the inner leaflet is smaller than that in the outer leaflet (Fig. <ref>(b)), which is the driving force that allows the NP to enter the vesicle spontaneously.






We summarize the dynamic behavior of a neutral NP with d/σ=5 in Fig. <ref>(a) and with d/σ=2 in Fig. <ref>(b).
Typical snapshots of the final state at t=1 × 7500 τ are shown in Fig. <ref>(c).
When the interaction parameter between the NP with d/σ=5 and the hydrophobic beads, v_ nt/k_ BT, is less than 2, the NP is adsorbed on the outer leaflet and does not enter the inside of the membrane or the vesicle (Fig. <ref>(a)).
At v_ nt/k_ BT=3, the NP can disturb the membrane structure and can remain in or pass through the lipid membrane.
On the other hand, a smaller NP with d/σ=2 is not adsorbed on the outer leaflet when the interaction parameters are too small (v_ nh/k_ BT=v_ nt/k_ BT=1).
Because the small NP interacts with only a few lipid beads, it desorbs from the outer leaflet.
Unlike the larger NP, the smaller NP can adsorb on the inner leaflet for v_ nh/k_ BT=3 and v_ nt/k_ BT=2.
Moreover, the smaller NP tends to localize in the hydrophobic region of the lipid membrane as it does not disturb the membrane structure greatly.
Some snapshots and density correlation functions for d/σ=2 are shown in Fig. S3 <cit.>.







 §.§ Charged NP







Next, we place a charged NP (z_ n≠ 0) on the outer leaflet.
First, the valence of an NP is set to z_ n=+50.
Typical dynamics for d/σ=5, v_ nh/k_ BT=2, and v_ nt/k_ BT=3 are shown in Fig. <ref>.
The parameters, except for the NP valence, are the same as for the calculation in Fig. <ref>.
To ensure data reproducibility, the results obtained by the other two calculations are shown in Fig. S7 <cit.>.
As shown in the snapshots in Fig. <ref>(a), the charged NP passes through the lipid membrane and is adsorbed on the inner leaflet.
This behavior is similar to that of a neutral NP (Fig. <ref>).
However, the membrane permeation speed for a charged NP is faster than that for a neutral NP.
For the charged NP, the NP is adsorbed on the inner leaflet at t=0.13 × 7500 τ, whereas the NP adsorption on the inner leaflet is observed from t=0.18 × 7500 τ for the neutral NP.
The Lennard–Jones potential between the NP and lipids, V^ (NL)_ LJ (black line in Fig. <ref>(b)), immediately decreases and the stepwise decrease observed in Fig. <ref>(b) does not occur.
In contrast to the neutral NP, the electrostatic interaction between the charged NP and the lipids, V^ (NL)_ el (red line in Fig. <ref>(b)), makes a strong contribution and the energy is significantly decreased.
The fast membrane permeation is attributed to this electrostatic energy gain. 


We show the density correlation functions of A- and B-lipids around the NP in Fig. <ref>(d) to investigate the process of the permeation through the lipid membrane (t=0.02 × 7500 τ).
Since the values of density correlation functions for charged lipid (A-lipid) and neutral lipid (B-lipid) are about 0.7 and 0.5 at r_ n/σ=3, respectively, the number of negatively charged lipids (A-lipids) is larger near the positively charged NP.
This may be caused by the electrostatic attraction between the NP and charged lipids (A-lipids).
Figure <ref>(e) shows the density correlation functions of a- to c-beads that make up the lipids around the NP.
There are slightly more hydrophilic beads near the NP (0.8 at r_ n/σ=3).
This effect is also small, because the indistinguishable profile is found in the neutral NP (z_ n=0) case as shown in Fig. <ref>(f).
The electrostatic attraction between the charged NP and charged lipids increases the membrane permeation speed; however, the permeation process is similar to that for the neutral NP.


We observe a crucial difference in the lipid distribution in the membrane between the charged and neutral NPs.
There is a charged lipid (A-lipid)-rich region in the inner leaflet near the charged NP (enlarged view in Fig. <ref>(a) and cross-sectional view in Fig. <ref>(h)).
Even in the outer leaflet, the charged lipids are accumulated in the protruding part of the vesicle (snapshot from the side view in Fig. <ref>(h)).
In Fig. <ref>(c), the electrostatic interaction between charged lipids, V^ (LL)_ el, increases over time.
This result also implies that the lipid distribution becomes heterogeneous.
The lipid accumulation is also indicated by the density correlation functions of A- and B-lipids in the inner and outer leaflets (Figs. <ref>(f) and (g)).
In both leaflets, the charged lipids (A-lipids) accumulate around the charged NP due to the electrostatic attraction between charged lipids and the charged NP.
Although the electrostatic repulsion between charged lipids is increased by the charged lipid accumulation, the system obtains a large energetic gain from the electrostatic attraction between the charged NP and charged lipids.
Because the energy gain shown in Fig. <ref>(b) overcomes the energy loss shown in Fig. <ref>(c), the charged lipids accumulate near the NP.
Moreover, the angle which determines the NP adsorption area is 273.8±13.3^∘ as shown in Fig. S11 <cit.>.
The angle becomes larger than that for neutral NP due to the contribution of the electrostatic interaction.


Neutral lipids (B-lipids) are abundant at r_ n/σ=7 in the inner leaflet (Fig. <ref>(f)), whereas there are fewer charged lipids (A-lipids).
Moreover, there are more neutral lipids (B-lipids) than charged lipids (A-lipids) around r_ n/σ=9 in the outer leaflet (Fig. <ref>(g)).
This region appears outside the peak of the charged lipids (A-lipids) because the charged lipids (A-lipids) accumulate near the NP and are depleted in the outer region.
In particular, the difference between the amounts of charged and neutral lipids is significant in the outer leaflet.
In the outer leaflet, the outside of the region where the NP is adsorbed to the membrane, that is, the part constricted by the protrusion of the vesicle, has a negative curvature, and it is difficult for charged lipids to localize in this region.
Thus, there is a region that is slightly neutral lipid-rich in the constricted part (side-view snapshot in Fig. <ref>(h)).


Figures <ref>(i) and (j) show the lipid distribution changes as functions of distance from the NP r_ n and time.
The difference in color indicates the difference in the values of density correlation functions between charged lipid (A-lipid) and neutral lipid (B-lipid).
Therefore, positive and negative values indicate abundance in charged lipids (A-lipids) and neutral lipids (B-lipids), respectively.
In the inner leaflet, a charged lipid (A-lipid)-rich region is formed around r_ n/σ=3 to 4 as shown in Fig. <ref>(i).
As the color in a charged lipid (A-lipid)-rich region gets darker over time, we can see that more charged lipids (A-lipids) are accumulated.
In addition, it is found that the charged lipid (A-lipid) is particularly depleted and the neutral lipid (B-lipid)-rich region is formed around r_ n/σ=6 to 8, which is outside the charged lipid (A-lipid)-rich region.
In the outer leaflet, a charged lipid (A-lipid)-rich region appears at r_ n/σ=3 to 4 at t<0.05 × 7500τ as shown in Fig. <ref>(j).
It can be seen that this region formed in a short time is a charged lipid (A-lipid)-rich region induced by NP adsorption to the outer leaflet.
The region enriched in charged lipids (A-lipids) moves away from the NP with the lapse of time.
This indicates that the NP passes through the lipid membrane and moves to the inner leaflet.
Therefore, the region rich in charged lipids (A-lipids) at r_ n/σ=6 to 7 can clearly be seen after t ≃ 0.4 × 7500τ.
Likewise, a neutral lipid (B-lipid)-rich region is formed along the outside of this charged lipid (A-lipid)-rich region.







We consider an increase in the NP valence for the charged NP with z_ n=+150.
In Fig. <ref>, we show the typical dynamic behavior of the charged NP with z_ n=+150, d/σ=5, v_ nh/k_ BT=2, and v_ nt/k_ BT=3.
The parameters, except for the valence of the NP, are the same as for the calculations in Figs. <ref> and <ref>.
To ensure data reproducibility, the results obtained by the other two calculations are shown in Fig. S8 <cit.>.
The charged NP permeates through the lipid membrane and is adsorbed on the inner leaflet (Fig. <ref>(a)).
This behavior is similar to the that of the neutral NP (Fig. <ref>) and that of the charged NP with z_ n=+50 (Fig. <ref>).
Because the valence of the NP (z_ n=+150) is larger than that in Fig. <ref> (z_ n=+50), the permeation speed is also faster.
The fast permeation is mediated by the large electrostatic energy gain by the attraction between the charged NP and charged lipids (Fig. <ref>(c)), although the change in the Lennard–Jones potential,  V^ (NL)_ LJ, in Fig. <ref>(b) is similar to that shown in Fig. <ref>(b).
In this case, there is also a charged lipid-rich region near the NP at t=1 × 7500 τ as shown in the enlarged view of Fig. <ref>(a).
In addition, the NP is completely wrapped by the lipid membrane at t=1 × 7500τ as shown in Fig. <ref>(a), and the adsorption angle becomes 360^∘ as shown in Fig. S11 <cit.>.


There is an interesting difference between the permeation processes for z_ n=+50 and +150.
The density correlation functions of A- and B-lipids near the NP at t=0.03 × 7500 τ are shown in Fig. <ref>(d) .
In the snapshot at t=0.03 × 7500 τ (Fig. <ref>(a)), the NP is just in the middle of the membrane.
In this step, the charged lipids (A-lipids) are accumulated near the charged NP (Fig. <ref>(d)) and are depleted outside the peak.
Moreover, the beads closest to the NP are the hydrophilic a-beads, followed by the hydrophobic b-beads, and then the hydrophobic c-beads (Fig. <ref>(e)), indicating that lipids are oriented with their hydrophilic heads toward the NP.
These results suggest that the charged lipids accumulate immediately near the NP, and the lipids are highly oriented toward the NP.
Therefore, a temporary channel covered with charged lipids is formed around the NP, the NP passes through the channel as shown in the enlarged view at t=0.03 × 7500τ of Fig. <ref>(a), 
and the channel disappears at t=0.12 × 7500 τ.


When the valence of a charged NP is increased to z_ n=+200, we observe a different permeation process through the lipid membrane.
The typical dynamic behavior of a charged NP with z_ n=+200, d/σ=5, v_ nh/k_ BT=2, and v_ nt/k_ BT=3 is shown in Fig. <ref>.
The parameters, except for the valence of the NP, are the same for the calculations in Figs. <ref>, <ref>, and <ref>.
To ensure data reproducibility, the results obtained by the other two calculations are shown in Fig. S9 <cit.>.
As shown in the snapshots in Fig. <ref>(a), the NP is adsorbed on the outer leaflet immediately, buried in the membrane at t=0.02 × 7500 τ, and covered with the lipids at t=0.03 × 7500 τ.
At t=0.1 × 7500 τ, the NP is completely wrapped with the lipid membrane and the adsorption angle becomes 360^∘ as shown in Fig. S11 <cit.>.
After t=0.1 × 7500 τ, the vesicle shape is elongated to fit the NP, and the bypass bilayer, which is under the NP, disappears at t=0.4 × 7500 τ.
Finally, the NP is adsorbed on the inner leaflet of the membrane and the charged lipid-rich region is observed near the NP at t=1 × 7500 τ as shown in Fig. <ref>(a).
Because the NP is covered with the lipid membrane and is transported to the interior of the vesicle, we call this process transportlike endocytosis.







The Lennard–Jones potential between the NP and lipids, V^ (NL)_ LJ, in Fig. <ref>(b) is similar to those in Figs. <ref>(b) and <ref>(b).
There is a significant decrease in the electrostatic interaction between the NP and charged lipids, V^ (NL)_ el, in Fig. <ref>(c).
Although the electrostatic interaction, V^ (NL)_ el, for z_ n=+150 also decreases dramatically, almost no transportlike endocytosis is observed.
This is because the energy gain by the electrostatic interaction, V^ (NL)_ el, does not compensate for the energy loss induced by the large membrane deformation for endocytic movement.


Here, we focus on the process of transportlike endocytosis.
From the density correlation functions of the A- and B-lipids near the NP at t=0.1 × 7500 τ in Fig. <ref>(d), the charged lipids (A-lipids) are accumulated near the NP by the strong electrostatic attraction.
This result is similar to that for the transient channel formation (Fig. <ref>).
Figure <ref>(e) shows the complicated peak appearance, with the hydrophilic a-beads beads closest to the NP, followed by the hydrophobic b-beads, hydrophobic c-beads, hydrophobic b-beads, and hydrophilic a-beads.
This peak pattern indicates bilayer formation around the NP, and the NP is covered with the lipid bilayer at t=0.1 × 7500 τ (Fig.  <ref>(a)).







Finally, we show the dynamic behavior of the charged NP with z_ n=-50, which has the same sign as the charged lipids.
Figure <ref> shows the typical dynamic behavior of the charged NP with z_ n=-50, d/σ=5, v_ nh/k_ BT=2, and v_ nt/k_ BT=3.
The parameters, except for the valence of the NP, are the same as for the calculations in Figs. <ref>, <ref>, <ref>, and <ref>.
To ensure data reproducibility, the results obtained by the other two calculations are shown in Fig. S10 <cit.>.
The NP is adsorbed on the outer leaflet and continues to adhere until t ∼ 0.14 × 7500 τ (Fig. <ref>(a)).
After t=0.14 × 7500 τ, the NP starts to bury into the membranes and is localized in the hydrophobic region of the membrane at t=0.19 × 7500 τ (Fig. <ref>(a)).
The Lennard–Jones potential between the NP and the lipids, V^ (NL)_ LJ, and the electrostatic interaction between the NP and the lipids, V^ (NL)_ el, are shown in Fig. <ref>(b).
The Lennard–Jones potential, V^ (NL)_ LJ, decreases quickly and the value remains almost constant until t ∼ 0.14 × 7500 τ.
This step indicates the adsorption of the NP on the outer leaflet.
The Lennard–Jones potential decreases again from t=0.14 to 0.2 × 7500 τ.
This process corresponds to the incorporation of the NP into the hydrophobic core of the membrane.
Because the charges of the NP and charged lipids (A-lipids) are the same sign, the electrostatic interaction, V^ (NL)_ el, has a positive value.
This repulsion tends to prevent the NP adsorption on the leaflets and the NP permeation through the membrane.


The density correlation functions at t=1 × 7500 τ are shown in Fig. <ref>(d) and (e).
There are many neutral lipids near the NP due to the electrostatic repulsion between the NP and the charged lipids (Fig. <ref>(d)).
Although there is a slight increase in the number of hydrophilic beads near the NP (Fig. <ref>(e)), the effect is small and the result is the same as that obtained for the neutral NP (Fig. <ref>(f)).
Therefore, the NP is in the hydrophobic region of the lipid membrane and the orientation of the neutral lipids near the NP is not affected by the NP, and behaviors such as transient channel formation and transportlike endocytosis are not observed.


While the NP is adsorbed on the outer leaflet (t=0.03 - 0.14 × 7500 τ), a neutral lipid (B-lipid)-rich region forms just under the NP (Fig. <ref>(h)).
This is because the charged lipids (A-lipids) near the NP are repelled from the region around the NP by the electrostatic repulsion and the neutral lipids (B-lipids) are present instead.
The concentration of charged lipids (A-lipids) in the outer region just under the NP hardly changes, even when the neutral lipid (B-lipid)-rich region is formed; thus, the electrostatic interaction between charged lipids (A-lipids), V^ (LL)_ el, in Fig. <ref>(c) remains almost constant.
The density correlation functions in the inner and outer leaflets are shown in Figs. <ref>(f) and (g).
In both leaflets, there are peaks for the neutral lipids near the NP.
Because the NP is adsorbed on the outer leaflet, the difference between the charged and neutral lipids is greater in the outer leaflet due to the large electrostatic repulsion.


We calculated the lipid distribution changes as functions of distance from the NP r_ n and time as shown in Figs. <ref>(i) and (j).
As in Figs. <ref>(i) and (j), positive and negative values indicate charged lipid (A-lipid)-rich and neutral lipid (B-lipid)-rich regions, respectively.
Figure <ref>(i) shows the lipid distribution change in the inner leaflet.
At t<0.15 × 7500τ, a region rich in neutral lipids (B-lipids) appears at r_ n/σ=6 to 8.
This formation is caused by the adsorption of NP to the outer leaflet.
The region rich in neutral lipids (B-lipids) suddenly moves to r_ n/σ=3 to 5 just before t=0.2 × 7500τ.
This indicates that the NP adsorbed on the membrane surface entered into the membrane.
From the snapshot of Fig. <ref>(a), it can be seen that the NP localizes in the lipid membrane at t=0.19 × 7500τ.
Because the NP keeps the localization in the lipid membrane, the neutral lipid (B-lipid)-rich region does not migrate greatly after t=0.2 × 7500τ.
The lipid distribution change in the outer leaflet is shown in Fig. <ref>(j).
A neutral lipid (B-lipid)-rich region is formed at r_ n/σ=3 to 5, and this region does not migrate during the calculation.
The neutral lipid (B-lipid)-rich region is formed by the NP adsorption to the outer leaflet or the NP localization in the lipid membrane.
The distances between NP and the outer leaflet for these two situations are almost the same.
Therefore, it is considered that the position of the neutral lipid (B-lipid)-rich region is not changed.






We summarize the dynamic behavior of charged NPs with d/σ=5 and 2 in Fig. <ref>, and the valences and the diameters of the NPs (z_ n, d/σ) are (+10, 5) in (a), (+50, 5) in (b), (+150, 5) in (c), (+200, 5) in (d), (-50, 5) in (e), (+10, 2) in (f), (+50, 2) in (g), (+150, 2) in (h), (+200, 2) in (i), and (-50, 2) in (j). 
Some snapshots and density correlation functions for d/σ=2 are shown in Fig. S3 <cit.>.
For the neutral NP, roughly v_ nt/k_ BT=3 is needed to localize in the membrane or permeate through the membrane (Figs. <ref>(a) and (b)).
The electrostatic attraction between the NP and lipids decreases the critical hydrophobic interaction strength for permeation through the membrane or localization in the membrane as v_ nt/k_ BT=2 (Figs. <ref>(a) and (f)).
For z_ n=+10, no charged lipid (A-lipid)-rich region is formed just under the charged NP, because the electrostatic repulsion between the charged lipids when the lipids are accumulated is larger than the electrostatic attraction between the NP and lipids due to the low valence of the NP.
For z_ n=+10, the smaller NP with d/σ=2 is more likely to be present in the membrane, similar to the neutral NP.
This similarity arises because the small NP does not disturb the ordering of the lipids when it is present in the membrane.


When the valence of the NP is z_ n=+50, the NP passes through the hydrophobic region of the membrane and is adsorbed on the inner leaflet, regardless of the values of v_ nh and v_ nt (Figs. <ref>(b) and (g)).
Although z_ n=+50 is too low to form a transient channel, it is high enough for the NP to pass through the membrane.
In Figs. <ref>(c) and (h), permeation through the membrane occurs via the transient channel for z_ n=+150, except for d/σ=5 and v_ nh/k_ BT=v_ nt/k_ BT=1.
At d/σ=5 and v_ nh/k_ BT=v_ nt/k_ BT=1, transportlike endocytosis is the most common dynamic behavior.
When z_ n is increased to +200, transportlike endocytosis occurs under many conditions for d/σ=5 (Fig. <ref>(d)).
Interestingly, transportlike endocytosis is not observed for z_ n=+200 and d/σ=2 (Fig. <ref>(i)).
This is because the lipid membrane must be curved sharply to cover the small NP with the lipid membrane in endocytosis, and the electrostatic attraction between the NP and lipids cannot compensate for the energy loss arising from this large deformation.


When v_ b is small, especially v_ nh, it is difficult for the negatively charged NP with z_ n=-50 to be adsorbed on the outer leaflet by overcoming the electrostatic repulsion between the NP and charged lipids for d/σ=5 (Fig. <ref>(e)).
If the Lennard–Jones potential is strong enough, the negatively charged NP can be adsorbed on the outer leaflet.
For large v_ nt, the NP can exist in the hydrophobic region of the lipid membrane through the favorable attraction between the NP and hydrophobic beads.
On the other hand, the smaller NP with d/σ=2 cannot be adsorbed on the outer leaflet (Fig. <ref>(j)).
Because the number of lipid beads interacting with the NP is less for the smaller NP, a sufficient attractive interaction is not achieved and the NP desorbs from the outer leaflet.








§ DISCUSSION



Several points merit further discussion.
Some experimental studies have reported that cationic dendrimers are internalized into cells via endocytosis <cit.>.
As the charge density of the dendrimers increased, the dendrimers were more easily internalized into cell <cit.>.
This result is good agreement with our results (Fig. <ref>), and endocytosis is frequently found for the large NP (d/σ=5) and high charge density NP (z_ n=+200) in our simulation.


In evaluating the relationship between experiments and our results, it is important to consider the behavior of multiple NPs and their cooperative motion.
Membrane protrusion was observed experimentally when NPs were added to vesicles  <cit.>.
This behavior was explained by the adsorption of NPs on the outer leaflet.
However, we observed membrane protrusion induced by NP adsorption on the inner leaflet.
In our calculations, NPs are relatively large since the vesicle size is small.
Therefore, the vesicle may be just stretched to fit the NP.
Also in the experiments, it is difficult to identify the location of NPs by microscopy due to the optical resolution.
It is important to calculate the interaction between a larger vesicle and multiple NPs and directly compare experimental results and simulation results in near future.
Moreover, domain formation induced by protein adsorption has been observed experimentally <cit.>.
Small domain formation has also been reported, and this is consistent with our findings and other simulation work <cit.>.
Scheve et al. reported that numerous small domains are formed by steric repulsion between adsorbed proteins <cit.>.
Calculations with multiple NPs may produce similar domain formation, which would be confirmed in future work.
Moreover, electrostatic repulsion as well as the steric repulsion between NPs could form positional ordering arrangements, such as two-dimensional crystals <cit.>.


Liu et al. investigated the interaction between charged NPs and vesicles consisting of zwitterionic lipids by CGMD simulations with a MARTINI force field <cit.>.
Although the net charge of the hydrophilic headgroup is neutral, they considered a headgroup with a dipole, such as phosphatidylcholine.
Charged NPs mainly interact with the positive charge on the headgroup, which is located near the membrane surface.
Consequently, the behavior of NPs depends on their charge density, ρ (e nm^-2); localization in the hydrophobic core is observed for ρ=-1 and +1, adsorption is observed for ρ=-3.58, adsorption or penetration is observed for ρ=-7.51, and desorption is observed for ρ=+7.15.
Based on the charge density of NPs, our NPs are ρ=+0.26 for (z_ n,d/σ) = (+10,5),  ρ=+1.3 for (z_ n,d/σ) = (+50,5), ρ=+3.9 for (z_ n,d/σ) = (+150,5),  ρ=+5.2 for (z_ n,d/σ) = (+200,5),  ρ=-1.3 for (z_ n,d/σ) = (-50,5),  ρ=+1.6 for (z_ n,d/σ) = (+10,2), ρ=+8.1 for (z_ n,d/σ) = (+50,2), ρ=+24.4 for (z_ n,d/σ) = (+150,2), ρ=+32.5 for (z_ n,d/σ) = (+200,2), and ρ=-8.1 for (z_ n,d/σ) = (-50,2).
These charge densities are summarized in Table S2 <cit.>. 
The behavior of the NPs summarized in Fig. <ref> is mainly consistent with that reported by Liu et al..
The essential differences between our study and Liu et al. are that we consider charged membranes (net charge of membrane is not zero), a more coarse-grained lipid molecule,
and 5000 lipid molecules, compared with 1182 used by Liu et al..
Although we consider a simple lipid molecule consisting of three beads connected linearly by a spring, this model is sufficient to consider the behavior of NPs with a lipid membrane qualitatively.
Our model reproduces the dynamic behavior of NPs with a charged lipid membrane for a wide range of interaction parameters, and NP valences and sizes.


An NP on an outer leaflet spontaneously passes through the lipid membrane and is adsorbed on the inner leaflet (Fig. <ref>) because the attractive energy between the NP and lipids is lower when the NP is adsorbed on the inner leaflet (Fig. <ref>(b)).
This difference in attractive energy between the NP on the outer and inner leaflets arises from the lipid membrane curvature.
To evaluate the effect of the membrane curvature on the NP uptake behavior, we perform the same calculations using a smaller vesicle consisting of 1000 lipid molecules. 
The vesicle consists of 500 charged lipids (A-lipids) and 500 neutral lipids (B-lipids) and the symmetric monovalent salt concentration is 100 mM.
We also perform the calculation until t = 10 × 7500τ to relax the lipid membrane structure sufficiently with w_ c^ AA/σ=1.8, w_ c^ BB/σ=1.7, and w_ c^ AB/σ=1.575.
The NP behavior is summarized in Fig. S2 <cit.>.
Some snapshots and density correlation functions are shown in Figs. S4 and S5 <cit.>.
The neutral NP and charged NPs with z_ n=+10 and -50 pass more easily through the membrane of the smaller vesicle (Figs. <ref>, <ref>, and S2).
For z_ n=+50 and +150, there is no large difference in NP behavior between the larger and smaller vesicles because the electrostatic interaction is dominant and the effect of the membrane curvature is negligible.
However, the transportlike endocytosis, which is seen for z_ n=+200, is not observed for the smaller vesicle.
Because of the large membrane curvature, a large membrane deformation is needed to cover the NP surface for the smaller vesicle.
Because this deformation costs a large amount of energy, transportlike endocytosis is suppressed.


In this study, the Debye-Hückel approximation is used as the electrostatic interaction, and it will be appropriate for low electric potential and/or high salt concentration.
For NPs with high valence, it is known that such a mean-field theory does not treat ion-ion correlation correctly.
Therefore, the calculation including the presence of ions explicitly will be important for future work.
Moreover, there is a hydrophobic core in the center of lipid membrane.
Although the electrostatic interaction across the lipid membrane is considered as one of the important interactions <cit.>,
it is almost decoupled due to low dielectric constant of the hydrophobic part.
Harries and his coworkers proposed a theoretical model for describing the endocytosis of charged globular protein <cit.>.
In this model, they considered that the outer leaflet of the lipid bilayer membrane wrapping an NP does not electrostatically interact with the NP.
Since we did not consider the low dielectric nature of hydrophobic core, our simulation overestimates the electrostatic interaction.
However, our CGMD simulation qualitatively consistent with the theoretical model, where the endocytosis occurs for charged protein with higher valence <cit.>.
For further quantitative comparison of our results with experimental and theoretical studies, we should propose a model including the ion-ion correlation and the low dielectric nature of hydrophobic core in future.


In the present study, we examine the vesicle configuration and the effect of the membrane curvature is different from that of the flat membrane.
For example, NPs are adsorbed onto inner leaflet under various conditions (Figs. <ref>, <ref>, <ref>, and <ref>) because the curvature of NPs and that of inner leaflet fit, and thus the NPs can interact with more lipid beads in the inner leaflet without large membrane deformation.
Therefore, this uptake mechanism is mediated by the membrane curvature.
Since the vesicular diameter is on the order of 10 nm and it is much smaller than the cell size, which is on the micrometer scale. It is difficult to directly extrapolate the obtained results to the situation of much larger vesicles or cells.
However, the cell is not completely spherical, and has various local curvatures.
NP uptake mediated by membrane curvature could occur at the curved parts in the cell.
Moreover, Harries et al. indicated that the inner leaflet composition of the lipid membrane wrapping an NP should be satisfied with the condition of charge matching with NP <cit.>.
In our calculation, the number of lipids which can contact with NP is about 100 for d/σ=5.
Although the condition of charge matching is z_ n≃ +100, we can see the endocytosis from z_ n=+200.
We consider two effects to explain this difference.
First, the theoretical model does not take into account the intermediate state between spherical protein adsorbing state and protein wrapped state.
Since the intermediate state is generally in a very high energy state, it cannot be discussed just by comparing the energies of the two stable states.
Second, it is considered that NPs in our calculations are too small to collect the sufficient number of charged lipids on the NP surface, 
and it can be difficult to gain sufficient electrostatic energy that overcomes the deformation energy for endocytosis.
Based on these points, it is also important to check the behaviors of larger vesicles, which needs to increase the speed and scale of calculation.
In future work, we will perform similar calculations in vesicles of different sizes to reveal more detailed dependence of the NP uptake behavior on the membrane curvature as well as endocytosis.







§ CONCLUSION



In this paper, we report the interaction between charged NPs and negatively charged binary lipid vesicles based on a CGMD simulation.
For a neutral NP, the NP can be incorporated into the membrane hydrophobic region and the interior of a vesicle when the hydrophobic interaction parameter, v_ nt, is large enough.
A positively charged NP can permeate through the membrane owing to the electrostatic interaction, even at smaller v_ nt.
Moreover, at higher valences, membrane permeation can occur easily and the NP is transported via a transient charged lipid channel or endocytosis.
We can classify the permeation processes as passing through the hydrophobic core, transient channel formation, and transportlike endocytosis based on the density correlation function.
When the NP is small (d/σ=2), transportlike endocytosis does not occur because a large membrane deformation is necessary.
A negatively charged NP can be adsorbed onto the membrane for the large hydrophilic parameter, v_ nh, region.
The spontaneous uptake into the vesicle is caused by membrane curvature.
Using the CGMD model, we investigated the NP behavior with a combination of various parameters.
Another interesting phenomenon is that the NP delivered into a cell via endocytosis is released from the endosome <cit.>.
Based on our model, the release of NP due to the change in the electrostatic interaction can be examined by considering the difference in ionic strengths between inside and outside of the cell.
In future, the cooperative behavior of multiple NPs, the NP uptake behavior for a phase-separated membrane, and the NP uptake in larger and smaller vesicles will be investigated to discuss the mechanism of NP uptake into cells and vesicles in more detail.







Most of the calculations were performed using the parallel computer “SGI UV3000” provided by the Center for Information Science at JAIST. 
We acknowledge support from the Bilateral Joint Research Project (Japan-Slovenia) from the Japan Society for the Promotion of Science (JSPS). 
N.S. acknowledges support from a Grant-in-Aid for Scientific Research (C) 
(Grant No. JP17K05610) from JSPS. 
H.I. acknowledges support from a Grant-in-Aid for JSPS Research Fellow (Grant No. JP18J00259). 








0

KA
S. Komura and D. Andelman,
Adv. Colloid Interface Sci., 208, 34 (2014).

Keller1
S. L. Veatch and S. L. Keller,
Phys. Rev. Lett., 89, 268101 (2002).

Keller2
S. L. Veatch and S. L. Keller,
Biophys. J., 85, 3074 (2003).

raft
K. Simons and E. Ikonen,
Nature, 387, 569 (1997).

Simons
K. Simons and J. L. Sampaio,
Cold Spring Harbor Perspect. Biol., 3, a004697 (2011).

Hamada
T. Hamada, R. Sugimoto, M. C. Vestergaard, T. Nagasaki, and M. Takagi,
J. Am. Chem. Soc., 132, 10528 (2010).

Ben-Tal
N. Ben-Tal, B. Honig, C. Miller, and S. McLaughlin,
Biophys. J., 73, 1717 (1997).

Hille
B. Hille, E. J. Dickson, M. Kruse, O. Vivas, and B.-C. Suh,
Biochim. Biophys. Acta, 1851, 844 (2015).

Shimokawa1
N. Shimokawa, M. Hishida, H. Seto, and K. Yoshikawa,
Chem. Phys. Lett., 496, 59 (2010).

Dimova1
C. C. Vequi-Suplicy, K. A. Riske, R. L. Knorr, and R. Dimova,
Biochim. Biophys. Acta, 1798, 1338 (2010).

Himeno1
H. Himeno, N. Shimokawa, S. Komura, D. Andelman, T. Hamada, and M. Takagi,
Soft Matter, 10, 7959 (2014).

Dimova2
B. Kubsch, T. Robinson, R. Lipowsky, and R. Dimova,
Biophys. J., 110, 2581 (2016).

May
S. May, D. Harries, and A. Ben-Shaul,
Phys. Rev. Lett., 89, 268102 (2002).

Shimokawa2
N. Shimokawa, S. Komura, and D. Andelman,
Phys. Rev. E, 84, 031919 (2011).

Shimokawa3
N. Shimokawa, H. Himeno, T. Hamada, M. Takagi, S. Komura, and D. Andelman,
J. Phys. Chem. B, 120, 6358 (2016).

Ito
H. Ito, Y. Higuchi, and N. Shimokawa
Phys. Rev. E, 94, 042611 (2016).

Himeno2
H. Himeno, H. Ito, Y. Higuchi, T. Hamada, N. Shimokawa, and M. Takagi,
Phys. Rev. E, 92, 062713 (2015).

Albanese
A. Albanese, P. S. Tang, and W. C. W. Chan,
Annu. Rev. Biomed. Eng., 14, 1 (2012).

Jaggers
R. W. Jaggers and S. A. F. Bon,
Soft Matter, 14, 6949 (2018).

virus
S. Mukherjee, R. N. Ghosh, and F. R. Maxfield,
Physiol. Rev., 77, 759 (1997).

Rouse
J. G. Rouse, J. Yang, A. R. Barron, N. A. Monteiro-Riviere,
Toxicol. in Vitro, 20, 1313 (2006).

Shvedova
A. Shvedova, V. Castranova, E. Kisin, D. Schwegler-Berry, A. Murray, V. Gandelsman, A. Maynard, and P. Baron,
J. Toxicol. Environ. Health Part A, 66, 1909 (2003).

Seib
F. P. Seib, A. T. Jones, and R. Duncan,
J. Control. Release, 117, 291 (2007).

Albertazzi
L. Albertazzi, M. Serresi, A. Albanese, and F. Beltram,
Mol. Pharmaceutics, 7, 680 (2010).

Chithrani
B. D. Chithrani, A. A. Ghazani, and W. C. W. Chan,
Nano Lett., 6, 662 (2006).

Bossi
E. Bossi, D. Zanella, R. Gornati, and G. Bernardini,
Sci. Rep., 6, 22254 (2016).

Ruan
G. Ruan, A. Agrawal, A. I. Marcus, and S. Nie,
J. Am. Chem. Soc., 129, 14759 (2007).

Chen
K. L. Chen and G. D. Bothun,
Environ. Sci. Technol., 48, 873 (2014).

Granick
B. Wang, L. Zhang, S. C. Bae, and S. Granick,
Proc. Natl. Acad. Sci. U.S.A., 105, 18171 (2008).

Wei
X. Wei, W. Jiang, J. Yu, L. Ding, J. Hu, and G. Jiang,
J. Hazard. Mater., 287, 217 (2015).

Li
S. Li and N. Malmstadt,
Soft Matter, 9, 4969 (2013).

Hinderliter
A. Hinderliter, R. F. F. Almeida, C. E. Creutz, and R. L. Biltonen,
Biochemistry, 40, 4181 (2001).

McLaughlin
S. McLaughlin and D. Murray,
Nature, 438, 605 (2005).

Scheve
C. S. Scheve, P. A. Gonzales, N. Momin, and J. C. Stachowiak,
J. Am. Chem. Soc., 135, 1185 (2013).

Banerji
S. K. Banerji and M. A. Hayes,
Langmuir, 23, 3305 (2007).

Bihan
O. Le Bihan, P. Bonnafous, L. Marak, T. Bickel, S. Trépout, S. Mornet, F. De Haas, H. Talbot, J.-C. Taveau, and O. Lambert,
J. Struct. Biol., 168, 419 (2009).

Miyake
M. Miyake, Y. Kakizawa, N. Tobori, M. Kurioka, N. Tabuchi, R. Kon, N. Shimokawa, Y. Tsujino, and M. Takagi,
Colloids Surf. B: Biointerfaces, 169, 444 (2018).

Takeuchi
S. Ichikawa, N. Shimokawa, M. Takagi, Y. Kitayama, and T. Takeuchi,
Chem. Commun., 54, 4557 (2018).

Vacha
R. Vácha, F. J. Martinez-Veracoechea, and D. Frenkel,
Nano Lett., 11, 5391 (2011).

Gompper
S. Dasgupta, T. Auth, and G. Gompper,
Nano Lett., 14, 687 (2014).

Lehn1
R. C. Van Lehn, P. U. Atukorale, R. P. Carney, Y-S. Yang, F. Stellacci, D. J. Irvine, and A. Alexander-Katz,
Nano Lett., 13, 4060 (2013).

Lehn2
R. C. Van Lehn and A. Alexander-Katz,
J. Phys. Chem. A, 118, 5848 (2014).

Ye_Li
Y. Li, X. Zhang, and D. Cao,
Soft Matter, 10, 6844 (2014).

Yue
T. Yue and X. Zhang,
ACS Nano, 6, 3196 (2012).

Liu
L. Liu, J. Zhang, Z. Zhao, Z. Mao, N. Liu, Y. Zhang, and Q. H. Liu,
Phys. Chem. Chem. Phys., 18, 31946 (2016).

Yang_Li
Y. Li,
Phys. Chem. Chem. Phys., 16, 12818 (2014).

Li_Gu
Y. Li and N. Gu,
J. Phys. Chem. B, 114, 2749 (2010).

Nangia
S. Nangia and R. Sureshkumar,
Langmuir, 28, 17666 (2012).

Cooke
I. R. Cooke, K. Kremer, and M. Deserno,
Phys. Rev. E, 72, 011506 (2005).

Sakaue
A. A. Zinchenko, T. Sakaue, S. Araki, K. Yoshikawa, and D. Baigl,
J. Phys. Chem. B, 111, 3019 (2007).

SI
See Supplementary Information.

Srivastava
S. Srivastava, D. Nykypanchuk, M. Fukuto, and O. Gang,
ACS Nano, 8, 9857 (2014).

Baciu
C. L. Baciu and S. May,
J. Phys. Condens. Matter, 16, S2455 (2004).

Wagner
A. J. Wagner and S. May,
Eur. Biophys. J., 36, 293 (2007).

Harries
D. Harries, A. Ben-Shaul, and I. Szleifer,
J. Phys. Chem. B, 108, 1491 (2004).

Vacha2
R. Vácha, F. J. Martinez-Veracoechea, and D. Frenkel,
ACS Nano, 6, 10598 (2012).








§ SUPPORTING INFORMATION




 §.§ Parameters of calculations






 §.§ Surface charge densities of NPs








 §.§ Relaxation of membrane structure


Figure S1 (a) shows snapshots of the membrane structure in the binary 
charged lipid mixture composed of 2500 charged lipids (A-lipids) (red or dark) and 
2500 neutral lipids (B-lipids) (yellow or light), where the cutoff lengths for the 
attractive interactions among hydrophobic beads were set to 
w_ c^ AA/σ=1.8, w_ c^ BB/σ=1.7, and
w_ c^ AB/σ=1.55.
The membrane separates into A-lipid-rich and B-lipid-rich phases at 
t=1 × 7500 τ, and there is a clearer phase-separated structure 
at t=10 × 7500 τ. 
The change in the attractive potential, V_at^(LL), is shown
 in Fig. S1(b). 
As the phase separation progresses, the potential decreases substantially 
until  t ≃ 2 × 7500 τ. 
At this time, the vesicle may reach the equilibrium state because the 
potential does not change significantly after t ≃ 2 × 7500 τ. 
Figure S1(c) shows the density correlation functions between lipid species 
at t=10 × 7500 τ. 
The density autocorrelation function of A-lipids (red line) intersects with 
the density cross-correlation function between the A- and B-lipids (green line). 
The spatial distance at this intersection is the typical domain size, and this 
result also indicates that phase separation occurs.
When the cutoff length, w_ c^ AB/σ, is changed to 
w_ c^ AB/σ=1.575, no phase separation occurs (Fig. S1(d)). 
No decrease in the attractive potential, V_at^(LL), is visible 
immediately after the start of the calculation in Fig. S1(e), and the density 
autocorrelation function of the A-lipids (red line) and the density 
cross-correlation function between the A- and B-lipids (green line) do not 
intersect at t=10 × 7500 τ in Fig. S1(f). 
Therefore, the phase separation does not occur for 
w_ c^ AB/σ=1.575. 
V_at^(LL) continuously increases during the calculation in 
Fig. S1(e) because the lipids gradually escape from the vesicle. 
We regard the structure at t=10 × 7500 τ as the static states of 
the preassembled vesicles; such states observed experimentally have been 
treated as practically static because of the much slower time scale of the 
complete melting of the structure. 
Our results are essentially the same as those reported in a previous study <cit.>.











 §.§ Calculation for smaller vesicle









 §.§ Snapshots and density correlation functions



















 §.§ Data reproducibility

Although we performed three calculations at least for each condition to ensure data reproducibility, 
we show a representative result in the calculations we performed in main text.
Here, we show the two other calculation results, and these two calculations will be referred to as 2nd calculation and 3rd calculation for convenience.


























 §.§ Adsorption angle








