Interplay between ligand mobility and nanoparticle geometry during cellular uptake of PEGylated liposomes and bicelles

Abstract

We explore the cellular uptake process of PEGylated liposomes and bicelles by investigating their membrane wrapping process using large-scale molecular dynamics simulations. We find that due to the mobility of ligands on the liposome/bicelle, the membrane wrapping process of a PEGylated liposome/bicelle can be divided into two stages, whose transition is determined by a critical wrapping fraction f_c; it is reached when all the ligands are exhausted and bound to receptors within the cell membrane. Before this critical scenario is approached, the grafted polyethylene glycol (PEG) polymers aggregate together within the membrane–wrapped region of the liposome/bicelle, driven by ligand–receptor binding. For wrapping fractions f > f_c, membrane wrapping cannot proceed unless a compressive membrane tension is provided. By systematically varying the membrane tension and PEG molar ratio, we establish phase diagrams about wrapping states for both PEGylated liposomes and bicelles. According to these diagrams, we find that the absolute value of the compressive membrane tension required by a fully wrapped PEGylated bicelle is smaller than that of the PEGylated liposome, indicating that the PEGylated bicelle is easily internalized by cells. Further theoretical analysis reveals that compared to a liposome, the flatter surface at the top of a bicelle makes it energetically more favored beyond the critical wrapping fraction f_c. Our simulations confirm that the interplay between ligand mobility and NP geometry can significantly change the understanding about the influence of NP geometry on the membrane wrapping process. It can help us to better understand the cellular uptake process of the PEGylated liposome/bicelle and to improve the design of lipid-like NPs for drug delivery.

---

1. Introduction

Nanoparticle (NP) based cancer nanomedicine aims to improve the diagnosis and treatment of diseases by engineering NPs to specifically identify or deliver drugs to tumor sites.^1–3 After being injected into the human body, NPs need to circulate in blood flow,⁴ penetrate leaky vessel walls under the enhanced permeability and retention (EPR) effect,⁵ diffuse through the extracellular matrix of tumor stroma,^6,7 and finally enter into tumor cells. During this process, there are many biological barriers confronted by NPs. For instance, serum proteins in the blood flow can detect and absorb NPs on their surfaces, and thus act as an indicator for macrophage cells to clear these NPs.⁸ The cell membrane is an additional physical barrier to inhibit the NPs from entering the cytoplasm of tumor cells.⁹ Therefore, compositions and surface modifications of NPs are of great importance to NP delivery efficiency.^10–14

Lipid-like NPs stand themselves out among various NP candidates.^11,15–17 For instance, a liposome that consists of a fluid or gel state lipid bilayer shell and an aqueous core is one of the first studied NPs as a drug carrier.^18–20 The phospholipid surface of a liposome makes it biologically inert, weakly immunogenic, and less toxic compared to other NP formulations.^15,21 Specifically, Doxil (PEGylated liposome-encapsulated doxorubicin) is the first US Food and Drug Administration approved NP-based delivery formulation for clinical applications.¹⁵ By grafting polyethylene glycol (PEG) polymers on their surface, liposomes demonstrate a prolonged blood circulation time, which significantly improves the drug accumulation within tumor sites through the EPR effect.²² By loading doxorubicin into the PEGylated liposomes, Doxil demonstrates a higher efficacy for cancer therapy than freely administered doxorubicin molecules.¹⁵ Nevertheless, liposomes still have some limitations, such as morphology instability,^18,23 low delivery efficiency²⁴ and difficulty in size control.^18,25

A bicelle that is composed of a single disc-like lipid bilayer has attracted a lot of attention due to its unique properties and its biocompatibility which it shares with the liposome.^26–28 For instance, bicelles are found to be able to penetrate through the narrow intercellular spaces of the stratum corneum and show a promising platform for dermal applications.²⁶ Furthermore, compared with liposomes, bicelles are found to have a better chance to be internalized by tumor cells. For example, Wang et al.²⁹ compared the cellular uptake of bicelles and liposomes with the modification of the octa-arginine (R8) sequence in four different cell lines, including MM-231, human breast cancer cell MCF-7, human umbilical vein endothelial cell HUVEC and murine macrophage cell RAW264.7. They found that bicelles are more efficiently internalized by all these cells. In particular, the number of bicelles internalized by MCF-7 cells is 2.5 times larger than that of liposomes. Additionally, Aresh et al.^30,31 also found that the cellular uptake efficiency of bicelles is significantly higher compared with that of liposomes, based on three different human cancer cell lines of CCRF-CEM, KB and OVCAR-8.

Despite these important and promising applications in nanomedicine, knowledge about the interactions between cells and a liposome/bicelle is still quite limited. One of the key reasons is the possible interplay between NP elasticity, NP geometry, and the ligand mobility during the cellular uptake process. First, a fluid or gel state liposome/bicelle is soft and has the ability to deform itself.^32–36 Second, a liposome is spherical, while a bicelle is disc-like. These two shapes have dramatically different curvature distributions on their surfaces, leading to different internalization kinetics.³⁷ Third, ligands tethered on PEGylated lipids can freely diffuse on the liposome/bicelle surface.^38,39 Each of these factors can influence the cellular uptake process of a PEGylated liposome/bicelle. For instance, due to the deformation of elastic NPs, a soft NP is less energetically favorably wrapped by the membrane than its rigid counterpart.^40,41 Compared with a spherical NP, a disc-like NP is more difficult to be fully wrapped by the cell membrane due to its highly curved surface edge.^42–44 Moreover, the mobility of ligands would induce the aggregation of PEG polymers and ligands, which suppresses the membrane wrapping.^38,39 However, it is still not clear how the NP elasticity, geometry and ligand mobility can play together to influence the whole membrane wrapping process. In particular, the question on why a bicelle is more efficiently internalized by cancer cells remains to be answered.

To fill the knowledge gap between experimental results and our current understanding, we explore the cellular uptake of a PEGylated liposome and bicelle through large-scale molecular dynamics simulations. In our simulations, the cellular uptake process is mimicked by a membrane wrapping process initiated by ligands on grafted PEG terminals and receptors in the cell membrane (cf.Fig. 1). To obtain a full picture for the membrane wrapping process of the PEGylated liposome/bicelle, the membrane tension and PEG molar ratio are systematically varied. In this work, we define the wrapping fraction f as the ratio between the membrane's wrapped surface area and the total surface area of the liposome/bicelle. As the ligands are mobile on the surface, the ligand–receptor binding fraction is defined separately. We find that due to the mobility of ligands on the liposome/bicelle, the membrane wrapping process can be divided into two stages, separated by a critical wrapping fraction f_c. This f_c is defined as the wrapping fraction f when all the ligands are exhausted and bound to receptors. As long as f < f_c PEG polymers aggregate within the membrane–wrapped region of the liposome/bicelle, driven by ligand–receptor binding, while for f > f_c the membrane wrapping cannot proceed unless a compressive membrane tension is provided. From the phase diagrams about different wrapping states, we find that the absolute value of the compressive membrane tension boundary for the PEGylated bicelle is smaller than that of the PEGylated liposome, indicating that the PEGylated bicelle is easily internalized. We further proceed to analyze the free energy changes due to PEG polymers, NP elastic deformations, and membrane bending. The free energy analysis reveals that the energy barrier induced by PEG polymers is much larger than that of the NP elastic deformation. These energy barriers can be overcome by ligand–receptor binding when f < f_c. Otherwise, the membrane wrapping is driven by a compressive membrane tension. Our theory confirms that due to its flatter surface at the top, a bicelle is easily wrapped compared to a liposome when f > f_c. Our simulations reveal the interplay between the mobility of PEG polymers and NP geometry and demonstrate that under certain situations a disc-like NP is energetically more favored than a spherical NP. These findings help us better understand the cellular uptake process of the PEGylated liposome/bicelle, and can thus become useful in the design of lipid-like NPs for drug delivery.

Fig. 1 Molecular structures of computational models: (A) PEGylated bicelle, (B) PEGylated liposome, and (C) lipid and PEG polymer models for both PEGylated bicelle and liposome. Lipid heads and tails in a liposome/bicelle are colored in light blue and ice blue, respectively. The PEG polymers are colored in blue. The ligands (targeting moieties) conjugated on the distal ends of PEG polymers are represented by red beads. (D) Lipid membrane with over-expressed receptors. Lipid heads and tails in the membrane are colored in green and gray, respectively. The molecular structure of receptors is the same as that of lipids in simulations. The bead colored in tan on the receptor head is the active site to specifically bind to a ligand. As highlighted in (A) and (B), the maximum radii of the liposome and bicelle are the same, sharing a value of 11r₀.

2. Computational model and methods

All coarse-grained molecular simulations performed in this work are based on the dissipative particle dynamics (DPD) method.^45,46 The basic interacting sites in DPD simulations are represented by soft beads. Between each pair of DPD beads, effective two-body interactions consist of three major forces:^45,46 a conservative force F^C, a random force F^R and a dissipative force F^D. Specifically, the conservative force between beads i and j is F^C_ij = a_ijω(r_ij)e_ij, where r_ij denotes the spatial distance between the two beads i and j, and e_ij is the unit vector pointing from i to j; a_ij represents the maximum repulsion force strength. The weighting factor ω(r_ij) is a normalized distribution function as ω(r_ij) = 1 − r_ij/r₀ for r_ij ≤ r₀, while ω(r_ij) = 0 for r_ij > r₀. Here r₀ is the cutoff distance for pairwise interactions. The random forces are specified by , where α represents a normal distributed Gaussian random number with zero mean and unit variance, Δt = 0.01τ () denotes the integration time step, and β_ij is a bead friction coefficient. The dissipative force is given by F^D_ij = −β_ijω²(r_ij)(e_ij·v_ij)e_ij, where v_ij is the relative velocity vector between beads i and j. All pair-wise interactions a_ij between different types of beads are listed in Table S1 of the ESI.†

Two different lipid models are adopted in our simulations to represent the lipid molecules in a liposome/bicelle and the membrane, respectively. A DPD lipid model that mimics 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) is utilized to assemble the liposome and bicelle. In this model, the head group of each lipid molecule is represented by three linearly connected hydrophilic beads, while each of the two tails is represented by 5 hydrophobic beads,^47,48cf. Fig. 1(C). Adjacent two beads are connected by a harmonic spring potential U_s1 = K_s1(r_ij − r_s1)². The stiffness of the head and tail groups is controlled by a harmonic bending potential applied on the adjacent three beads U_θ1 = K_θ1(ϑ − ϑ₀₁)², where ϑ denotes a bending angle. Due to the flexible head group in this DPPC lipid model, the energy penalty to form the bilayer edge in the bicelle is relatively small, facilitating the formation of a larger sized bicelle. The line tension of the edge in this model is around λ = 1.4k_BT/r₀. Because of the favorable energy state, a bicelle will transform into a liposome only above a critical size. It is impossible in the simulation to create a liposome and bicelle with an identical lipid number. Following the idea in experiments,^29,31 the liposome and bicelle are assumed to have the same size, around R = 11r₀ (as denoted in Fig. 1), in our simulations. Our liposome and bicelle have 800 and 392 lipid molecules, respectively. A DPD lipid model that mimics 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) is utilized to form the planar membrane. In this lipid model, two lipid tails (carrying four tail beads each) are connected with two head beads, respectively, cf. Fig. 1(D). The head group contains three head beads. Adjacent beads in each lipid molecule are connected by a harmonic spring potential U_s2 = K_s2(r_ij − r_s2)². The stiffness of the lipid tails is guaranteed by the bending potential U_ϑ2 = K_ϑ2(1 − cos ϑ). Using this DMPC lipid model, the tension of a planar bilayer is linearly related to the lipid molecular area.^49,50 The stretch modulus of the membrane can be obtained from the related slope (cf. ESI Fig. S1†) as K_A = 17.42k_BT/r²₀. The bending rigidity of the membrane is given by^51,52κ_m = K_Ad²_hh/48 with κ_m ≈ 6k_BT, which is within the range of 5–50k_BT obtained in the experiments.^53,54 The size of the planar membrane bilayer used to investigate the membrane wrapping is (70 × 70)r²₀. Please refer to the ESI† for details about the interaction potential parameters and the calibration of mechanical properties.

A hydrophilic PEG polymer in our DPD simulations is modeled by a linear chain consisting of coarse-grained monomers. The PEG monomers are linearly connected by the harmonic bond potential U_s3 = K_s3(r_ij − r_s3)², with spring stiffness K_s3 = 2111.3k_BT/r²₀ and equilibrium distance r_s3 = 0.4125r₀. The semi-flexibility of the PEG polymer is taken into account by adding the bending potential U_ϑ3 = K_ϑ3(cos ϑ − cos ϑ₀₃)², with bending stiffness K_θ3 = 16.4946k_BT, and equilibrium bending angle θ₀₃ = 130° between each three consecutive monomers. Such a DPD PEG model could correctly reproduce the conformation of a PEG polymer in water, including the radius of gyration and end-to-end distance, as shown in our previous studies.^37,47 To describe the PEGylated lipid, one end of the PEG polymer is bonded to the lipid head bead through a harmonic bond potential. In addition, the monomers at the free end of PEG polymers are defined to act as targeting moieties (ligands) (Fig. 1). The polymerization degree of PEG polymers in our simulation is set as N = 20 (representing a molecular weight of around 660 Da), falling within the typical range of 500–3000 Da in experiments.^8,55,56 Four different sets of PEG polymer molar ratios of 40, 50, 60, and 70 mol% will be investigated. The corresponding numbers of PEG polymers n_c on liposomes are 320, 400, 480 and 560, respectively. Those n_c on bicelles are 156, 196, 235 and 274, respectively.

To mimic the ligand–receptor interaction, we assume that 50% of lipid molecules in the planar bilayer act as receptors. In this way, the receptor diffusion will not be a factor that limits the efficiency of membrane wrapping in our simulations.^57,58 Receptors in the planar membrane follow the same configuration as a lipid, cf. Fig. 1(D), with the head bead acting as an active site to interact with the ligand. The ligand–receptor interaction follows a modified Lennard–Jones potential^37,47 as U_ij = 4ε_ligand[(σ_b/r_ij)¹² − (σ_b/r_ij)⁶] − U_cut, when r_ij ≤ r_cut and U_ij = 0 otherwise. Here, r_cut = r₀ for a short-range attractive interaction and U_cut = 4ε_ligand[(σ_b/r₀)¹² − (σ_b/r₀)⁶]. Due to the difference between this ligand–receptor cutoff distance and bilayer thickness, ligands on the NPs can only interact with the receptors on the outer leaflet. The equilibrium distance is fixed by using σ_b = 0.624r₀. Additionally, the repulsive force is limited to 25k_BT/r₀. We use ε_ligand = 12k_BT. The single ligand–receptor binding energy is then around 6.8k_BT.⁴⁷ This pair-wise interaction setting between ligands and receptors is a commonly used strategy in simulations to speed up the membrane wrapping process.^14,32,33,36 Different from the valence-limited interactions,⁵⁹ the pair-wise potential can lead to multivalent ligand–receptor interactions. The maximum bound receptors for a ligand is 6 in our simulations. Such a multivalent interaction is also experimentally possible by engineering antibodies.⁶⁰

Due to the elasticity of the lipid membrane, the tension of membranes in cells can be adjusted, changing from 0.01 to 10 mN m⁻¹.⁶¹ The N-varied DPD method is applied during the membrane wrapping process to ensure a constant tension of the planar membrane.^{37,47,62–64} In practice, boundaries of the lipid bilayer are treated as a lipid reservoir for addition and removal of lipids. If the lipid number per unit area is larger (or smaller) than a target density ρ, lipid molecules will be deleted (or inserted) into this boundary region to maintain a constant lipid number density. Meanwhile, the corresponding number of water molecules will be randomly inserted into (or deleted from) the simulation box to ensure a constant bead density of 3.0/r³₀ in the simulation box. The target density ρ is taken based on the relationship between the membrane tension and lipid area given in ESI Fig. S1.† By using the N-varied DPD protocol, the lipid density in the membrane can easily be controlled to maintain the membrane's lateral tension during the membrane wrapping process.

The physical length corresponding to our simulation unit is obtained by comparing the membrane thickness in simulations d_HH = 4r₀ to the thickness of a real membrane, d_HH ≈ 3.53 nm,⁶⁵ indicating r₀ = 0.9 nm. The experimental lipid lateral diffusion coefficient of DMPC is D ≃ 5 μm² s⁻¹.⁶⁶ In our simulations, the lipid lateral diffusion coefficient D_lipid ≃ 7.3 × 10⁻²r²₀/τ is calculated by averaging the values under different membrane tensions. Hereafter, we deduce the physical time scale τ = 11.8 ns. Note that these spatial and temporal mappings are only used to approximate the length and time scales of all DPD simulations, which are different from the real length and time scales in all-atom molecular dynamics simulations.

3. Results and discussion

3.1. Membrane wrapping of the PEGylated liposome and bicelle

PEG polymers aggregate during the membrane wrapping process. We firstly investigate the membrane wrapping process of the PEGylated liposome and bicelle under a membrane tension of −0.74k_BT/r²₀, where 1k_BT/r²₀ corresponds to 5 mN m⁻¹. The PEG polymer molar ratio in both the liposome and bicelle is 60% mol. As given in Fig. 2, the PEGylated liposome and bicelle are initially placed above the membrane at a distance of 5r₀. The PEG polymers are almost evenly distributed on the liposome and bicelle surfaces at t = 0. Due to the attraction between the ligand and receptor, the PEGylated liposome quickly adheres on the membrane. At the same time, the membrane starts to bend and partially wrap the liposome. Furthermore, the PEG polymers aggregate within the wrapped region of the liposome because of their mobility.³⁸ The aggregation is more pronounced at t = 118 μs in Fig. 2(A), which results in a PEG/ligand free region at the top of the yet unwrapped part of liposome. As investigated in our previous work,³⁸ the ligand depletion on top of a NP will suppress the membrane wrapping on the PEGylated liposome compared to a rigid NP. At t = 472 μs, the membrane begins to protrude and spread over the top part of the liposome. At t = 944 μs, the liposome is fully wrapped (f = 1). Associated with this full wrapping, the PEG polymers rearrange and redistribute over the entire liposome surface. As shown in the lower panel of Fig. 2(A), the liposome adjusts itself to deform along with the evolution of the wrapped state.

Fig. 2 Membrane wrapping process of the PEGylated liposome and bicelle. (A) The snapshots in the upper panel represent the membrane wrapping process of the PEGylated liposome. The snapshots in the lower panel show the corresponding morphology change of liposome. (B) The snapshots in the upper panel represent the membrane wrapping process of the PEGylated bicelle. The snapshots in the lower panel show the corresponding morphology change of the bicelle. Water beads are not shown for clarity. The membrane tension is maintained at −0.74k_BT/r²₀. The PEG molar ratio of both the PEGylated liposome and bicelle is 60%. The scale bar indicates 20 nm.

The wrapping process of the PEGylated bicelle shows a similar configurational pathway in terms of the PEG polymers’ aggregation and rearrangement in Fig. 2(B). Half of the bicelle is quickly wrapped by the membrane at t = 59 μs. PEG polymers start to aggregate within the wrapped region. At t = 118 μs, with more PEG polymers aggregated, the bicelle elongates into a strip-like shape. Compared to a disc-like geometry, this strip-like shape has a longer edge, in which it produces a large curvature. This large curvature at the edge can release the increased interactive energy between aggregated PEG polymers.⁶⁷ At t = 236 μs, the bicelle slightly rotates to release the membrane energy. At the same time, lipids in the membrane start to protrude. Following the membrane protruding and PEG polymer rearrangement, the bicelle recovers its initially disc-like shape when fully wrapped at t = 336 μs. It is also interesting to note that during the wrapping processes for both the liposome and bicelle, the ligands arrange in an ordered strip-like pattern. This specific arrangement of ligands should be caused by the competition between the PEG free energy and ligand–receptor binding energy.

To explore more details behind the membrane wrapping process, we proceed to calculate both the ligand–receptor binding fraction and wrapping fraction f. The binding fraction is defined as the ratio of the bonded ligand number to total ligand number in the liposome/bicelle. As given in Fig. 3(A) and (B), due to the mobility of ligands on PEG polymers, the increment of the binding fraction is much faster than that of f in both the PEGylated liposome and the bicelle. As highlighted by the dashed lines, it is worth noticing that all the ligands in the liposome and bicelle are bound to receptors in the membrane and exhausted at the critical wrapping ratio f_c ≈ 0.7. It means that no driving force can be provided by ligands at the later wrapping state (f > f_c). To fully wrap the liposome/bicelle, the key question is what will drive the membrane wrapping after ligands have been exhausted.

Fig. 3 Detailed information about the membrane wrapping process of the PEGylated liposome and bicelle. (A and B) The comparison of the ligand–receptor binding ratio and membrane wrapping fraction for the liposome and bicelle, respectively. (C) The evolution of asphericity for both the liposome and bicelle. The asphericity is defined by r²_gz − 0.5(r²_gx + r²_gy), where r²_gx, r²_gy and r²_gz are the principal moments of the gyration tensor (r²_gx ≤ r²_gy ≤ r²_gz) for a NP. (D) The evolution of average end-to-end distance for PEG polymers in both the PEGylated liposome and bicelle.

We further calculate the asphericity of the liposome/bicelle, and the end-to-end distance of grafted PEG polymers. The asphericity of the liposome/bicelle, defined by r²_gz − 0.5(r²_gx + r²_gy), is plotted against f in Fig. 3(C). Here, r²_gx, r²_gy and r²_gz are the ordered principal moments of the NP gyration tensor (r²_gx ≤ r²_gy ≤ r²_gz). A large asphericity value indicates a pronounced anisotropic shape.⁶⁸ As we can see in Fig. 3(C), the variation of asphericity is more pronounced for the bicelle. In particular, near the wrapping fraction of f ≈ 0.75, the asphericity of the bicelle has increased by a factor 2 compared with its initial value. This large variation corresponds to the strip-like shape of the bicelle at t = 118 μs, cf. Fig. 2(B), when all the PEG polymers are aggregated within the wrapped region. The average end-to-end distance R_ee = 〈R²_ee〉^1/2 as given in Fig. 3(D) is another indicator of PEG polymer configurations. The R_ee for both the liposome and bicelle increases with increasing wrapping fraction below f ≈ 0.7, since PEG polymers still tend to aggregate within the wrapping region. Afterwards, their R_ee values change slightly with the redistribution of PEG polymers.

PEGylated liposome and bicelle cannot be fully wrapped under positive membrane tension. To explore possible different states for both the PEGylated liposome and bicelle during the membrane wrapping process, we next investigate the situation with a slightly positive membrane tension 0.09k_BT/r²₀. Additionally, a series of different PEG molar ratios, ranging from 40% mol to 70% mol, are considered. As we can see in Fig. 4, under the positive membrane tension, all of the PEGylated NPs cannot be fully wrapped by the membrane, and remain in a partial wrapped state. Furthermore, due to their mobility, all PEG polymers are aggregated within the wrapped region of the liposome/bicelle, leaving a ligand-free region at the top. Specifically, the morphology of the liposome at the trapped state is highly dependent on the PEG molar ratio as shown in Fig. 4(A). At 40% mol and 50% mol, the liposomes deform into ellipsoidal shapes. The membrane wrapping process of 50% mol is given in ESI Fig. S3.† At 60% mol, it is interesting to find that lipids at one side of the liposome's contact edge form a tubular shape. This tubular shape should be caused by the continuously aggregating PEG polymers, which results in highly packed PEG polymer (and linked lipids) states in the contact edge (cf. ESI Fig. S3†). The formation of the tubular shape is favorable as it releases the increased steric interaction for both PEG polymers and lipids. More importantly, the liposome with 70% mol PEG polymers ruptures during the wrapping process and deforms into a strip-like shape at the end of simulation. The rupture of the liposome should be induced by a decrease of the water storage space. Details about the rupture process of the PEGylated liposome during the membrane wrapping process are available in the ESI.†

Fig. 4 Membrane wrapping states of the PEGylated (A) liposome and (B) bicelle with different PEG molar ratios under a membrane tension of 0.09k_BT/r²₀. All the simulations are run long enough without further wrapping state changes.

Concerning the PEGylated bicelle, its configuration within the trapped state is also dependent on the PEG molar ratio. At 40% mol and 50% mol, the trapped bicelle maintains its disc-like shape, while it finally deforms into the above-mentioned strip-like shape at 60% mol and 70% mol. The membrane wrapping processes for bicelles with 50% mol and 70% mol PEG polymers are given in ESI Fig. S4.†

Phase diagrams reveal a compressive (negative) membrane tension boundary. In view of the different wrapping states of the PEGylated liposome/bicelle, to come up with a complete picture, we proceed to investigate the influence of both the membrane tension and PEG molar ratio. The membrane tension is varied from negative tension of −0.74k_BT/r²₀ to positive tension of 1.0k_BT/r²₀, beyond which the membrane has large chances to ripple or rupture in the simulation. A phase diagram is obtained for both the PEGylated liposome and bicelle as shown in Fig. 5. For the PEGylated liposome, according to Fig. 5(A), there are three main regions in the phase diagram: (1) when the PEG molar ratio is high, the liposome ruptures. This upper boundary about the PEG molar ratio limits the maximum number of decorated ligands. (2) If the PEG molar ratio is not large enough or the membrane tension is positive, the liposome is only partially wrapped by, and trapped inside the membrane. (3) A fully wrapped state for the PEGylated liposome exists only for the situation with a high PEG molar ratio and compressive (negative) membrane tension. For the phase diagram of PEGylated bicelles shown in Fig. 5(B), there are two main regions: (1) when the membrane tension is positive and the PEG molar ratio is small, the PEGylated bicelle is partially wrapped by, and trapped inside the membrane. (2) If the membrane tension is sufficiently small and the PEG molar ratio is high, the PEGylated bicelle is fully wrapped by the membrane. Additionally, both the PEGylated liposome and bicelle have extra minor regions about tubular and strip-like shapes, respectively.

Fig. 5 Phase diagrams for membrane wrapping of the PEGylated (A) liposome and (B) bicelle under the influence of the membrane lateral tension and PEG molar ratio. The red circles represent the cases of being fully membrane wrapped. The green triangles represent the rupture of PEGylated liposomes. The blue squares are the cases of being partially membrane wrapped. The half open blue squares in the PEGylated liposome represent the cases that the liposome is protruded into a tubular shape. The open squares in the PEGylated bicelle are the cases that the bicelle is elongated to a strip-like shape.

Comparing these two phase diagrams, we find that the PEGylated bicelle has a larger chance to be fully wrapped. First, the possibility for the PEGylated liposome to rupture sets an upper limit boundary to the PEG polymer molar ratio, which does not exist for the PEGylated bicelle. Additionally, the minimum required PEG molar ratio for the PEGylated liposome is larger than that of the PEGylated bicelle. More importantly, it is interesting to find that the full wrapping (f = 1) of the PEGylated liposome and bicelle can only occur when the effective membrane tension is compressive. Due to the mobility of PEG polymers, the ligands on the liposome/bicelle are exhausted before its fully wrapped state is reached. It is the compressive tension that makes it energetically favorable and drives the membrane wrapping after the ligands are exhausted. Furthermore, the PEGylated bicelle has a larger membrane tension boundary compared to that of the PEGylated liposome, i.e., the absolute value of the compressive membrane tension boundary for the PEGylated bicelle is smaller than that of the PEGylated liposome. This suggests that a PEGylated bicelle is more easily fully wrapped by the membrane than its liposomic counterpart. The key question about the membrane tension boundary is why the PEGylated bicelle is more favorable. In living cells, the required compressive membrane tension can be produced through active mechanisms.⁶⁹ For instance, the contraction of actomyosin can produce a compressive stress on the cytoskeletal network, resulting in a locally compressive membrane tension.^70,71 It might also suggest that for the cellular uptake of the PEGylated liposome/bicelle, other biological mechanisms should be involved apart from the receptor-mediated membrane wrapping.

Considering the importance of ligand exhausting, we calculate and compare the critical wrapping fraction f_c for a PEGylated liposome and bicelle. When f > f_c, any further membrane wrapping is driven by the compressive membrane tension alone. As shown in Fig. 6, the f_c is increasing with the increment of the PEG molar ratio. At the same PEG molar ratio, the PEGylated liposome and bicelle share a similar f_c value. Moreover, within the studied range of the PEG molar ratio, all the f_c values are larger than 0.68.

Fig. 6 Critical wrapping fraction against the PEG molar ratio for the PEGylated liposome and bicelle. The critical wrapping fraction is defined by the wrapping fraction for the ligand–receptor binding fraction approaching 1.0. The membrane tension for the cases listed in the figure is −0.74k_BT/r²₀.

3.2. Free energy analysis of the membrane wrapping process

To gain deep insights into why a PEGylated bicelle is easier to be wrapped than a PEGylated liposome under otherwise similar conditions, we perform a free energy analysis for the membrane wrapping process. During this process, free energy changes are composed of four major parts: (1) the ligand–receptor binding energy ΔF_LR that provides the driving force for membrane wrapping. (2) The energy change caused by PEG polymer aggregation and configurational changes, ΔF_PEG. (3) The energy change induced by the deformation of a soft liposome or bicelle, ΔF_NP. (4) The energy change associated with the membrane bending, ΔF_m. In our simulations, the individual binding strength between the ligand and receptor is about 6.8k_BT. Therefore, ΔF_LR has a linear relationship with the ligand–receptor binding fraction. In the following parts, we will analyze these free energy changes. Note that the possible free energy changes induced by the membrane fluctuation and translational entropy loss of receptors (and solvent molecules) are small compared with the above free energy changes.^72,73 We ignore their contributions to simplify the free energy analysis.

Aggregation of PEG polymers leads to a large energy barrier. We firstly investigate the free energy change of PEG polymers. Along with the aggregation and rearrangement of PEG polymers, their free energy change consists of three parts:^37,47,74 (1) elastic energy change ΔF_el, associated with the stretching or compressing of polymer chains, in analogy to elastic springs. (2) Interaction free energy change ΔF_int, caused by the change of interaction strength between PEG polymers. (3) Translational free energy change ΔF_trans, reflecting the translational entropy loss of polymers. Therefore, the total free energy change of PEG polymers can be expressed as

ΔF_PEG = ΔF_el + ΔF_int + ΔF_trans. (1)

Following the approaches in our previous studies, each part of ΔF_PEG can be estimated by feeding self-consistent mean field (SCF) theory with the information about PEG polymer configurations and local volume fractions from DPD simulations.^{37,38,47,75,76} Based on the SCF theory, the mean elastic energy F_el per chain is linearly proportional to the mean squared end-to-end distance 〈R²_ee〉, F_el/n_ck_BT = 3〈R²_ee〉/2R²₀, where R₀ represents the equilibrium span of an unconstrained PEG polymer with polymerization degree N, and n_c denotes the number of PEG chains. We employ R²₀ ∝ 〈R²_ee〉_w using the available N-dependent 〈R²_ee〉_w values for a single PEGchain in water.⁴⁷ The mean interaction free energy of each chain F_int is quantified through the spatially inhomogeneous volume fraction φ(r) of PEG monomers. , where f_mix(φ) = (φ² + φ³)/v is a mixing free energy density, and v = 0.0633 nm² is the excluded volume for a PEG monomer. The integral extends over the whole PEG–populated volume. The translational free energy ΔF_trans is directly related to the variable distribution of mobile PEG tethering points on the NP surface. It is estimated by the volume fraction profile as . Note that the intrinsic relaxation time of a PEG polymer with polymerization degree N = 20 is smaller than 100 ps,⁷⁷ while a time step Δt in our DPD simulations corresponds to 118 ps. It indicates that PEG polymers are able to relax themselves on a time scale that is short compared with the total simulation time, the equilibrium SCF approach can be considered applicable, and we can extract ΔF_PEG in the course of simulation time.

Taking the membrane wrapping processes in Fig. 2 for example, the local volume fraction distributions φ(r) of PEG polymers change dramatically during their aggregation and redistribution, as quantified in Fig. 7. Under the conditions of 60% mol PEG polymers and −0.74k_BT/r²₀ membrane tension, the evolution of φ(r) can be divided into two different stages. Before the exhausting of ligands (t = 472 μs and t = 118 μs for the PEGylated liposome and bicelle, respectively), the PEG volume fraction in the wrapped part of the liposome/bicelle significantly increases due to the aggregation of PEG polymers and at the expense of a dramatically reduced φ(r) value, associated with the depletion of PEG polymers, in the unwrapped part. This highly inhomogeneous PEG volume distribution manifests itself in coexisting highly compacted and reduced occupied regions for PEG polymers, disfavored by F_mix. Along with the spreading of the membrane over the top part of the liposome/bicelle, the PEG volume fraction starts to recover its initial homogeneous distributed state. Note that if the liposome/bicelle remains at the trapped state as for the positive membrane tension cases, the PEG volume fraction distribution is maintained at the highly inhomogeneous state as shown in ESI Fig. S5.†

Fig. 7 (A and B) Cross-sectional views of PEG polymer volume fraction distribution during the membrane wrapping process: (A) PEGylated liposome and (B) PEGylated bicelle.

Combining the above PEG volume fraction distribution and end-to-end distance of PEG polymers in Fig. 3, we can estimate the free energy change of each of the three contributions to ΔF_PEG as shown in Fig. 8. The dashed lines denote the location of the f_c values. For the PEGylated liposome, ΔF_PEG dramatically increases along with the ligand–receptor binding and PEG aggregation as long as f < f_c. The maximum value of ΔF_PEG is around 1000k_BT for the liposome. The changes of elastic energy ΔF_el and interaction energy ΔF_int are seen to dominate the PEG energy increment. Compared with ΔF_el and ΔF_int, the translational entropy change of PEG polymers ΔF_trans is near zero and can be neglected. With recovering to a homogeneous PEG distribution, the ΔF_PEG decreases. The PEG free energy changes of the PEGylated bicelle follow a similar trend, while they are overall smaller due to the smaller number of PEG polymers in our bicelle; the maximum energy barrier is at ΔF_PEG ≈ 600k_BT for the bicelle. Note that for the trapped state of the liposome/bicelle, the ΔF_PEG will maintain at its maximum value due to the end-to-end distance increment and highly inhomogeneous PEG volume fraction distribution as given in ESI Fig. S5.† In short, the PEG polymers produce a large energy barrier during the membrane wrapping process for both the PEGylated liposome and bicelle.

Fig. 8 Free energy changes of PEG polymers during the membrane wrapping process for the PEGylated (A) liposome and (B) bicelle, respectively.

Nanoparticles’ elastic deformation energy is small. We proceed to analyze the energy change ΔF_NP accompanying the deformations of the liposome and bicelle. We consider the membrane wrapping processes in Fig. 2 as representative examples to estimate the energy penalty due to NP deformation. For the deformed liposome as given in Fig. 2(A), its surface area is usually assumed to be time-independent. Additionally, the energy change caused by an osmotic pressure variation is negligible because of the relaxed initial state of the liposome.⁴¹ Therefore, the major contribution to ΔF_NP for the liposome should be the bending energy associated with curvature variation. This bending energy can be estimated by fitting the shape of the liposome by an ellipsoid as we proposed in our previous work.³⁸ Based on this fitted shape, the bending energy of the liposome can be obtained by following Helfrich's expression, , where κ_lip is the assumed spatially homogeneous bending rigidity of the liposome, and c¹_lip and c²_lip are the principal curvatures of the fitted ellipsoidal surface. For a bicelle, its total deformation energy consists of bending energy and edge energy.^78,79 The bending energy change of the bicelle is negligible, due to the small variation of curvature during the entire membrane wrapping process, cf. ESI Fig. S6.† To estimate the energy change from edge length variation, we fit the bicelle by an ellipse. The energy change of the bicelle can then be obtained as E_bic = Lλ, where L is the perimeter of the fitted ellipse, and λ is the line tension of the bicelle. The estimated free energy changes ΔF_NP for both the liposome and bicelle are given as a function of f in Fig. 9. The energy variation of the liposome is smaller than 25k_BT during the whole wrapping process, consistent with the observed small asphericity (Fig. 3). For the bicelle, its energy increases more significantly with its elongation. Its maximum value is around 70k_BT. However, compared with the energy change ΔF_PEG due to PEG polymers, F_NP is much smaller for both the liposome and bicelle.

Fig. 9 Free energy changes from elastic deformations of the liposome and bicelle versus wrapping fraction f.

Critical membrane tension is crucial for full wrapping. As we mentioned, the PEG polymers produce a large free energy barrier during the membrane wrapping process. This barrier can be overcome below the f_c by consuming the ligand–receptor binding. However, beyond the f_c, ligands on the liposome/bicelle are exhausted. The membrane wrapping cannot proceed unless any other driving force is provided. As revealed by the phase diagrams, a compressive membrane tension is required for a successful wrapping. Therefore, the membrane energy change plays a crucial role for f > f_c. Here, we provide a theoretical analysis about membrane energy change to reveal and compare the required compressive membrane tensions for the liposome and bicelle.

The membrane energy E_m consists of two parts: (1) membrane bending energy E_mbend and (2) membrane tension energy E_mtens,

E_m = E_mbend + E_mtens. (2)

To estimate E_m, we assume that an uncorrugated membrane wraps around a solid NP. The geometry of the solid NP is taken based on the shape of the liposome and bicelle in Fig. 2. As illustrated in Fig. 10(A), three different shapes are considered. The first one is a spherical NP to represent the liposome. The radius of the spherical NP is taken as 14.5r₀, based on the radius of the liposome 11r₀ and the thickness 3.5r₀ of the PEG polymer shell. The second one is an oblate NP to mimic the undeformed bicelle. Its principle radii are R₁ = R₂ = 14.5r₀ and R₃ = 6r₀, obtained from the principle radius of the bicelle and thickness of the PEG polymer shell. The third one is a biaxial ellipsoidal shape to represent the elongated bicelle, with the principle radii R₁ = 19r₀, R₂ = 10r₀ and R₃ = 4r₀. These values are taken from the fitted ellipsoidal shapes of the entire PEGylated bicelle during the wrapping process (cf. ESI Fig. S7†). To simplify the analysis, we only consider the wrapped (red) part of NPs, which can effectively catch the main membrane energy penalty during the wrapping process.^36,80,81 With these parameters at hand, the contributions to the membrane energy (Eq. 2) at a certain wrapping fraction f can be expressed as

E_mtens(f) = σΔS, (4)

where κ_m is the bending rigidity of the membrane; c¹_m and c²_m are its spatially dependent principal curvatures. The integration is taken over the wrapped region of the NP; σ is the membrane tension and ΔS is the excess membrane area caused by membrane bending. The details about the expressions of curvatures and excess membrane area are given in the ESI.† We firstly calculate the membrane bending energy E_mbend for these three NPs. As given in Fig. 10(B), the E_mbend is always positive and monotonically increasing with increasing f due to the large curvature of NPs. Additionally, near the critical wrapping fraction (f ≈ 0.7), the E_mbend of the bicelle is larger than that of the liposome, consistent with our previous observations.³⁷ However, the derivative of E_mbend with respect to f is smaller for the bicelle than for the liposome. Note that, according to our previous discussions, the membrane wrapping is driven by the compressive membrane tension for f > f_c. To gain an energetically favorable state by membrane wrapping, the derivative of the membrane energy with respect to f should be smaller than zero, dE_m(f)/df ≤ 0. Based on the derivative of the membrane energy, we can thus conclude that σ(f) ≤ σ_critical must hold, where the critical membrane tension is given by

Fig. 10 Membrane energy change. (A) Illustration of the membrane wrapping for NPs with different geometries at f = 0.85. The red part of NPs denotes the wrapped region. The unwrapped part of NPs is colored by blue. (B) Membrane bending energy versus f. The inserted figure is the derivative of membrane bending energy versus f. (C) Critical tension versus f for different NP geometries. The inserted figure is the amplified plot for the regime 0.7 < f < 1.0. The critical tension σ_critical is defined by the minimum required membrane tension that can drive the membrane wrapping, when the wrapping is only driven by compressive tension.

The corresponding numerically evaluated σ_critical for all three shapes is given in Fig. 10(B). Comparing the critical membrane tension for the liposome and bicelle, it is interesting to find that there are two regions. For wrapping fractions f < 0.7, the σ_critical for the bicelle is smaller than that of the liposome because of its highly curved edges. However, for f > 0.7, the σ_critical for the bicelle is larger than that of the liposome, which indicates that it is easier to wrap the bicelle after f = 0.7. After f = 0.7, both the σ_critical values for the liposome/bicelle are around σ_critical = −0.1k_BT/r²₀. Interestingly, this value is consistent with the membrane tension boundary we found in the phase diagrams (Fig. 5). Moreover, the transition wrapping fraction of f = 0.7 is in good agreement with the critical wrapping fraction f_c in Fig. 6. These agreements further prove that our membrane energy analysis can estimate the key energy barrier for the bilayer wrapping on a NP. In short, our theoretical analysis explains why the bicelle is energetically more favorable than the liposome after the exhausting of ligands.

The free energy analysis above reveals that the membrane tension required by the disc-like bicelle is larger than that of the liposome. It suggests that under the driving of compressive membrane tension beyond the f_c, the disc-like bicelle is energetically more favorable than the spherical liposome. Our result is different from the previous conclusion that driven by the ligand–receptor binding, a disc-like NP is energetically difficult to be fully wrapped compared to a spherical NP.^42–44 Our results indicate that the interplay between ligand/PEG mobility and NP geometry can significantly change the existing picture about the influence of NP geometry on the membrane wrapping process. To confirm the interplay between ligand/PEG mobility and NP geometry, we further analyze the critical membrane tension for spherical and oblate NPs with different aspect ratios at f > 0.5. The surface areas of spherical and oblate NPs are considered identical. As given in ESI Fig. S8,† all the oblate NPs with aspect ratio >1 require larger membrane tensions than spherical NPs. Furthermore, the wrapping fraction, beyond which an oblate NP is energetically favorable, is decreasing with increasing aspect ratio. This fact indicates that driven by the compressive membrane tension, an oblate NP with a larger aspect ratio is more easily fully wrapped in the later membrane wrapping stage (at f > f_c) due to its flat surface on top.

4. Conclusions

In this work, we performed large scale DPD simulations to understand the membrane wrapping processes of the PEGylated liposome and bicelle. We find that during the wrapping process, characterized by the wrapping fraction f, PEG polymers on the liposome and bicelle aggregate within their wrapped region due to the mobility of polymers and ligand–receptor binding. This aggregation promotes the ligand–receptor binding in the early membrane wrapping state. But the quick ligand and receptor binding results in a critical wrapping fraction f_c, after which the ligands are exhausted. According to the f_c, the entire membrane wrapping process for the PEGylated liposome/bicelle can be divided into two different stages: (1) as long as f < f_c, the membrane wrapping is driven by ligand–receptor binding. (2) For f > f_c, no driving force can be provided by ligand–receptor binding. Membrane wrapping cannot proceed unless a compressive membrane tension is provided. Furthermore, by systematically varying the molar ratio of PEGpolymers and membrane tension, we find that the PEGylated liposome is overall more difficult to be fully wrapped than the PEGylated bicelle because of two major reasons: (1) the possibility of rupture of the liposome at a high PEG molar ratio sets the upper limit of the ligand number. Such rupture did not occur for the bicelle under all PEG molar ratios studied in our simulations. (2) The absolute value of the compressive membrane tension boundary of the PEGylated bicelle is smaller than that of the PEGylated liposome, which indicates that the bicelle is the one that is easier to be fully wrapped.

Our free energy analysis revealed that PEG polymer aggregation leads to a large free energy barrier, while the energy barrier caused by NP deformation is relatively small and negligible. By analyzing the membrane energy, we find that the absolute value of the compressive membrane tension required by a disc-like bicelle is smaller than that of a spherical liposome, which suggests that the disc-like bicelle is energetically more favorable than the spherical liposome at f > f_c, where a compressive membrane tension is required to provide the driving force. The compressive membrane tension in living cells can be produced through active mechanisms.^70,71,82 Our results confirm that the interplay between ligand mobility and NP geometry can significantly change our understanding about the influence of NP geometry on the membrane wrapping process. This work can also help understand the cellular uptake process of the PEGylated liposome and bicelle, which might improve the design of new lipid-like drug delivery platforms.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Z. S., H. Y. and Y. L. are grateful to support from the Department of Mechanical Engineering at the University of Connecticut. This research benefited in part from the computational resources and staff contributions provided for Booth Engineering Center for Advanced Technology (BECAT) at the University of Connecticut. Part of this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which was supported by National Science Foundation grant number ACI-1053575. The contribution by M. K. was promoted by the Swiss National Science Foundation through grants 200021_156106 and 200021L_185052. This work was partially supported by a fellowship grant from GE's Industrial Solutions Business Unit under a GE-UConn partnership agreement. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of Industrial Solutions or UConn.

References

1. B. Y. Kim, J. T. Rutka and W. C. Chan, N. Engl. J. Med., 2010, 363, 2434–2443.
2. M. E. Davis, J. E. Zuckerman, C. H. J. Choi, D. Seligson, A. Tolcher, C. A. Alabi, Y. Yen, J. D. Heidel and A. Ribas, Nature, 2010, 464, 1067.
3. S. Nie, Nanomed., 2010, 5, 523–528.
4. H. Ye, Z. Shen, L. Yu, M. Wei and Y. Li, Proc. R. Soc. London, Ser. A, 2018, 474, 20170845.
5. R. A. Petros and J. M. DeSimone, Nat. Rev. Drug Discovery, 2010, 9, 615.
6. E. A. Sykes, Q. Dai, C. D. Sarsons, J. Chen, J. V. Rocheleau, D. M. Hwang, G. Zheng, D. T. Cramb, K. D. Rinker and W. C. Chan, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, E1142–E1151.
7. Q. Dai, S. Wilhelm, D. Ding, A. M. Syed, S. Sindhwani, Y. Zhang, Y. Y. Chen, P. MacMillan and W. C. Chan, ACS Nano, 2018, 12, 8423–8435.
8. C. D. Walkey, J. B. Olsen, H. Guo, A. Emili and W. C. W. Chan, J. Am. Chem. Soc., 2012, 134, 2139–2147.
9. S. Zhang, H. Gao and G. Bao, ACS Nano, 2015, 9, 8655–8671.
10. Z. Shen, M.-P. Nieh and Y. Li, Polymers, 2016, 8, 83.
11. A. Albanese, P. S. Tang and W. C. Chan, Annu. Rev. Biomed. Eng., 2012, 14, 1–16.
12. L. Y. Chou, K. Ming and W. C. Chan, Chem. Soc. Rev., 2011, 40, 233–245.
13. Y. Li, Y. Lian, L. T. Zhang, S. M. Aldousari, H. S. Hedia, S. A. Asiri and W. K. Liu, Interface Focus, 2016, 6, 20150086.
14. H.-M. Ding and Y.-Q. Ma, Biomaterials, 2012, 33, 5798–5802.
15. Y. C. Barenholz, J. Controlled Release, 2012, 160, 117–134.
16. A. C. Eifler and C. S. Thaxton, Biomed. Nanotechn., Springer, 2011, pp. 325–338.
17. M. Wang, K. Alberti, S. Sun, C. L. Arellano and Q. Xu, Angew. Chem., Int. Ed., 2014, 53, 2893–2898.
18. A. Sharma and U. S. Sharma, Int. J. Pharm., 1997, 154, 123–140.
19. T. Lian and R. J. Ho, J. Pharm. Sci., 2001, 90, 667–680.
20. Z. Shen, A. Fisher, W. K. Liu and Y. Li, Engineering of Biomaterials for Drug Delivery Systems, Elsevier, 2018, pp. 1–26.
21. M. L. Immordino, F. Dosio and L. Cattel, Int. J. Nanomed., 2006, 1, 297.
22. S. M. Moghimi and J. Szebeni, Prog. Lipid Res., 2003, 42, 463–478.
23. P. Marmottant, T. Biben and S. Hilgenfeldt, Proc. R. Soc. London, Ser. A, 2008, 1781–1800.
24. G. H. Petersen, S. K. Alzghari, W. Chee, S. S. Sankari and N. M. La-Beck, J. Controlled Release, 2016, 232, 255–264.
25. F. Olson, C. Hunt, F. Szoka, W. Vail and D. Papahadjopoulos, Biochim. Biophys. Acta, Biomembr., 1979, 557, 9–23.
26. L. Barbosa-Barros, G. Rodríguez, C. Barba, M. Cócera, L. Rubio, J. Estelrich, C. López-Iglesias, A. de la Maza and O. López, Small, 2012, 8, 807–818.
27. Z. Dai, S. Hameed and P. Bhattarai, Front. Chem., 2018, 6, 127.
28. U. H. Dürr, R. Soong and A. Ramamoorthy, Prog. Nucl. Magn. Reson. Spectrosc., 2013, 69, 1.
29. X. Wang, L. Lin, R. Liu, M. Chen, B. Chen, B. He, B. He, X. Liang, W. Dai and H. Zhang, et al. , Adv. Funct. Mater., 2017, 27, 1700406.
30. W. Aresh, Ph.D. thesis, University of Connecticut, U.S.A., 2016.
31. W. Aresh, Y. Liu, J. Sine, D. Thayer, A. Puri, Y. Huang, Y. Wang and M.-P. Nieh, J. Biomed. Nanotechnol., 2016, 12, 1852–1863.
32. T. Yue and X. Zhang, Soft Matter, 2013, 9, 559–569.
33. V. Schubertová, F. J. Martinez-Veracoechea and R. Vácha, Soft Matter, 2015, 11, 2726–2730.
34. T. Yue, Y. Xu, M. Sun, X. Zhang and F. Huang, Phys. Chem. Chem. Phys., 2016, 18, 1082–1091.
35. S. Li, Z. Luo, Y. Xu, H. Ren, L. Deng, X. Zhang, F. Huang and T. Yue, Biochim. Biophys. Acta, Biomembr., 2017, 1859, 2096–2105.
36. R. Guo, J. Mao and L.-T. Yan, ACS Nano, 2013, 7, 10646–10653.
37. Y. Li, M. Kröger and W. K. Liu, Nanoscale, 2015, 7, 16631–16646.
38. Z. Shen, H. Ye, M. Kröger and Y. Li, Nanoscale, 2018, 10, 4545–4560.
39. L. Di Michele, P. K. Jana and B. M. Mognetti, Phys. Rev. E, 2018, 98, 032406.
40. Z. Shen, H. Ye and Y. Li, Phys. Chem. Chem. Phys., 2018, 20, 16372–16385.
41. X. Yi, X. Shi and H. Gao, Phys. Rev. Lett., 2011, 107, 098101.
42. L. Chen, S. Xiao, H. Zhu, L. Wang and H. Liang, Soft Matter, 2016, 12, 2632–2641.
43. A. H. Bahrami, Soft Matter, 2013, 9, 8642–8646.
44. S. Dasgupta, T. Auth and G. Gompper, Soft Matter, 2013, 9, 5473–5482.
45. P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett., 1992, 19, 155.
46. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423.
47. Y. Li, M. Kröger and W. K. Liu, Biomaterials, 2014, 35, 8467–8478.
48. R. D. Groot and K. Rabone, Biophys. J., 2001, 81, 725–736.
49. A. Grafmüller, J. Shillcock and R. Lipowsky, Biophys. J., 2009, 96, 2658–2675.
50. A. Grafmüller, J. Shillcock and R. Lipowsky, Phys. Rev. Lett., 2007, 98, 218101.
51. M. Mutz and W. Helfrich, J. Phys., 1990, 51, 991–1001.
52. R. Goetz, G. Gompper and R. Lipowsky, Phys. Rev. Lett., 1999, 82, 221.
53. R. S. Gracià, N. Bezlyepkina, R. L. Knorr, R. Lipowsky and R. Dimova, Soft Matter, 2010, 6, 1472–1482.
54. G. Niggemann, M. Kummrow and W. Helfrich, J. Phys. II, 1995, 5, 413–425.
55. E. C. Cho, L. Au, Q. Zhang and Y. Xia, Small, 2010, 6, 517–522.
56. A. de la Zerda, S. Bodapati, R. Teed, S. Y. May, S. M. Tabakman, Z. Liu, B. T. Khuri-Yakub, X. Chen, H. Dai and S. S. Gambhir, ACS Nano, 2012, 6, 4694–4701.
57. H. Gao, W. Shi and L. B. Freund, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 9469–9474.
58. R. Vácha, F. J. Martinez-Veracoechea and D. Frenkel, Nano Lett., 2011, 11, 5391–5395.
59. C. Knorowski and A. Travesset, Soft Matter, 2012, 8, 12053–12059.
60. Á. M. Cuesta, N. Sainz-Pastor, J. Bonet, B. Oliva and L. Álvarez-Vallina, Trends Biotechnol., 2010, 28, 355–362.
61. N. C. Gauthier, T. A. Masters and M. P. Sheetz, Trends Cell Biol., 2012, 22, 527–535.
62. B. Hong, F. Qiu, H. Zhang and Y. Yang, J. Phys. Chem. B, 2007, 111, 5837–5849.
63. H.-M. Ding and Y.-Q. Ma, Biomaterials, 2014, 35, 8703–8710.
64. Y. Li, X. Zhang and D. Cao, Nanoscale, 2015, 7, 2758–2769.
65. N. Kučerka, Y. Liu, N. Chu, H. I. Petrache, S. Tristram-Nagle and J. F. Nagle, Biophys. J., 2005, 88, 2626–2637.
66. G. Orädd, G. Lindblom and P. W. Westerman, Biophys. J., 2002, 83, 2702–2704.
67. U. Dahal, Z. Wang and E. E. Dormidontova, Macromolecules, 2018, 51, 5950–5961.
68. J. Liang, P. Chen, B. Dong, Z. Huang, K. Zhao and L.-T. Yan, Biomacromolecules, 2016, 17, 1834–1844.
69. M. Dogterom and G. Koenderink, Nat. Mater., 2011, 10, 561.
70. G. T. Charras, J. C. Yarrow, M. A. Horton, L. Mahadevan and T. Mitchison, Nature, 2005, 435, 365.
71. G. Charras and E. Paluch, Nat. Rev. Mol. Cell Biol., 2008, 9, 730.
72. J. Liu, G. E. Weller, B. Zern, P. S. Ayyaswamy, D. M. Eckmann, V. R. Muzykantov and R. Radhakrishnan, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 16530–16535.
73. R. Bradley and R. Radhakrishnan, Polymers, 2013, 5, 890–936.
74. R. Nap, Y.-Y. Won and I. Szleifer, Soft Matter, 2012, 8, 1688–1700.
75. Z. Shen, H. Ye, M. Kröger and Y. Li, Phys. Chem. Chem. Phys., 2017, 19, 13294–13306.
76. Z. Shen, D. T. Loe, J. K. Awino, M. Kröger, J. L. Rouge and Y. Li, Nanoscale, 2016, 8, 14821–14835.
77. J. Y. Wong, T. L. Kuhl, J. N. Israelachvili, N. Mullah and S. Zalipsky, Science, 1997, 275, 820–822.
78. W. Shinoda, T. Nakamura and S. O. Nielsen, Soft Matter, 2011, 7, 9012–9020.
79. C. Huang, D. Quinn, Y. Sadovsky, S. Suresh and K. J. Hsia, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 1910–2915.
80. J. Mao, P. Chen, J. Liang, R. Guo and L.-T. Yan, ACS Nano, 2016, 10, 1493–1502.
81. F. Frey, F. Ziebert and U. S. Schwarz, Phys. Rev. Lett., 2019, 122, 088102.
82. J. Solon, J. Pécréaux, P. Girard, M. C. Fauré, J. Prost and P. Bassereau, Phys. Rev. Lett., 2006, 97, 098103.

---

Footnote

† Electronic supplementary information (ESI) available: Details of the computational models for PEGylated liposomes and lipid membrane, computational methods and additional simulation results. See DOI: 10.1039/C9NR02408E

---