Predicting asymmetric phospholipid microstructures in solutions

Abstract

Asymmetric phospholipid microstructures, such as asymmetric phospholipid membranes, have potential applications in biological and medicinal processes. Here, we used the dissipative particle dynamics simulation method to predict the asymmetric phospholipid microstructures in aqueous solutions. The asymmetric phospholipid membranes, tubes and vesicles are determined and characterized by the chain density distributions and order parameters. The phase diagrams are constructed to evaluate the effects of the chain length on the asymmetric structure formations at equilibrium states, while the average radius of gyration and shape factors are calculated to analyze the asymmetric structure formations in the non-equilibrium processes. Meanwhile, we predicted the mechanical properties of the asymmetric membranes by analyzing the spatial distributions of the interface tensions and osmotic pressures in solutions.

---

Introduction

Asymmetric microstructures, such as asymmetric membranes, possess asymmetric structures on each side.¹ Such asymmetric microstructures result in unusual properties, including automatic cleaning, unidirectional transport, and water–oil separation, and they differ from those in common membranes.^2–5 For example, an asymmetric membrane can be used for unidirectional transmission due to the different degrees of wettability on both sides of the membrane.⁴ Phospholipid membranes with symmetric structures exist extensively in cells and play an important role in various biological processes.² However, asymmetric phospholipid microstructures, either membranes or tubes, could lead to novel properties and potential applications in biological processes. Hence, exploring asymmetric phospholipid microstructures with distinct two surfaces is worthwhile.

Phospholipid molecules have amphipathic properties with one hydrophilic functional group and one or two hydrophobic fatty acid groups in water solutions; these properties can be induced by solution molecules to form rich varieties of structures, such as membranes, tubes, and vesicles.^6–10 For example, a series of lamellar structures, namely, membranes (including liquid crystallized, gel, and fluid lamellae), have been observed for phospholipids in water solutions.⁶ Phase diagrams of these structures have been also constructed using various parameters, where stable structures are separated by phase boundaries distinguished from first-order phase transitions.^8,11–14 Previous experiments have indicated that water concentrations can induce the inverted hexagonal phase transiting to the lamellar phase by using complementary techniques of X-ray scattering and terahertz spectroscopy.¹⁵ Experimental results showed that the phase regions of the lipid membrane are determined by the water concentration in the water solutions.⁸ Other methods, such as electric field and shear flow, have also been used to expand the lamellar phase space in phase diagrams due to the special interest in lamellar structures. For example, under shear flows, lamellar microstructures can orient along the direction of the flows and enlarge the phase spaces for phospholipid membranes in the phase diagrams.¹⁶ These phase diagrams provide information on membrane structure distributions in the phase parameter space under various conditions at equilibrium states.

However, rapid processing usually causes polymers to reach non-equilibrium conformation, which leads to various polymer properties at different molecular levels; thus, the dynamics process plays an important role in deciding the material properties.¹⁷ Many recent experiments and computer simulations have revealed the dynamic processes of phospholipid membranes in solutions.^16–20 For instance, increasing the content of branched fatty acids can increase the fluidity of membranes and change the kinetics of microbial membranes.¹⁸ Another research direction is to explore the mechanical properties, such as local pressure and interfacial permeability, of symmetric phospholipid membranes via experimentations and simulations.^21–25 A previous experiment revealed the pressure response of cholesterol to solid-loaded phospholipid multilayers; it proved that cholesterol at a certain concentration decreases the tendency of separation and presents a new path for the construction of stable model membranes.²¹ Meanwhile, a theoretical approach based on self-consistent field theory has been developed to obtain the local stress across the membranes and interfaces in soft matter; this approach predicts the membrane lateral stress profile in the regions of the head and tail and the stress between hydrophobic and hydrophilic interfaces.²⁴ Investigating asymmetric structures, either asymmetric particles or membranes, is meaningful. Many experiments and simulations have focused on asymmetric particles.^26–29 Asymmetric particles have various potential applications in the biological field. For example, cancer therapeutic materials can be developed by studying Janus particles.³⁰ The research scope should be extended to asymmetric membranes, tubes and vesicles, which exist naturally in biological systems, due to their asymmetric properties.

Several experimental and simulation studies have focused on asymmetric lipid membranes.^31–38 For example, a plasma membrane with asymmetric complex structures was investigated in a previous work through large-scale molecular simulation (MD) simulations, where the membrane model included many lipid species that were asymmetrically distributed across the two leaflets.³⁵ Another atomistic MD simulation revealed that short chains of one bilayer leaflet cannot interdigitate with those of the other leaflet and lead to increased free space in the middle of the bilayer.³⁶ A coarse-grained (CG) MD simulation revealed that asymmetric structures assemble for multiple types of phospholipids in water solutions.³⁴ These contributions concentrated on the structural complexity of individual membranes with multiple types of phospholipid molecules at either the atomistic or CG level. The kinetics and mechanical properties of asymmetric phospholipid membranes, which are relevant to practical applications, still need to be investigated.^39–41 A recent experiment showed that the resulting microstructure can be altered by changing the chain length of phospholipid molecules.⁴² This finding suggests that the novel phenomenon probably appears in the self-assembly of multiple types of phospholipid molecules. In this work, we used the dissipative particle dynamic (DPD) method based on the CG model to predict the structural formation mechanism and performances of several asymmetric phospholipid microstructures in water solutions, which differs from the previous studies that examined the complexity of individual membranes. In the current work, we determined the asymmetric membranes, tubes and vesicles for the phospholipid polymers in the phase diagrams in terms of different chain lengths; we focused on the dynamic processes for the asymmetric phospholipid microstructures; we predicted the mechanical properties of asymmetric phospholipid membranes. Section II describes the model and method, and Section III presents the results and discussion. The summary is presented in Section IV.

Method and model

Method

The DPD method based on the CG model is suitable for simulating the hydrodynamic behavior of complex systems and soft matter.^11,43–46 In DPD simulations, each DPD bead represents a cluster of molecules and obeys Newton's law. Forces, including conservation, dissipative, and random, appear in a pairwise manner between the i-th and j-th pairs of particles. Specifically, conservation force F^C_ij is derived from a potential, dissipative force F^D_ij is introduced to decrease the radial velocity difference between particles, and random force F^R_ij is a stochastic force along the line connecting the particle centers. The total force on the i-th bead in the DPD system can be expressed aswhere the subscript index ij refers to the relative quantity between the i-th and j-th beads. a_ij is a particle interaction constant representing the maximum repulsive force between the i-th and j-th beads. r_ij is the relative distance, r_ij = r_i − r_j is the relative position, and v_ij = v_i − v_j is the relative velocity, where between the i-th and j-th beads. γ is the friction coefficient governing the dissipative force, and σ is the noise amplitude determining the strength of random forces. The relationship between the two parameters is σ² = 2γk_BT, where k_B is the Boltzmann constant and T is the absolute temperature. ζ_ij represents a random uniformly distributed value between 0 and 1 with a Gaussian distribution and unit variance. Usually, standard values of σ = 3.0 and γ = 4.5 are used in simulations.⁴⁷ The weight function w(r_ij) can be expressed aswhere r_c is a cut-off distance.

Model

An asymmetric membrane (Fig. 1) was developed based on the CG model by taking several carbon atoms together and grouping them into one bead. Two types of phospholipids (type one on the left and type two on the right in Fig. 1) were considered. The hydrophilic beads (HB) in types one and two are shown in red and pink, respectively, and the hydrophobic tails (TB) in types one and two are shown in blue and yellow, respectively. The two types of phospholipids are set to have one linear head chain and two linear tail chains, which have been extensively used in previous simulations.^15,16,48,49 The bonded beads in phospholipids are connected by an additional elastic harmonic forcewhere k_s and r_s are the spring constant and equilibrium bond length between the i-th and j-th beads, respectively. In the present simulation, we used k_s = 100.0 and r_s = 0.7r_c, which are similar to the values in previous studies.^48,49 To achieve the bending force of the phospholipid molecule imposed on two consecutive bonds, we considered the following relationship:

F^θ = −∇[k_θ(θ − θ₀)²], (4)

where k_θ is the bending constant and θ₀ is the equilibrium angle. In the present simulation, we used k_θ = 6.0 and θ₀ = π for three consecutive HBs or TBs in each chain, k_θ = 3.0 and between the head and tail chains are set for the last two beads in head chain and the first bead in tail chain, and k_θ = 4.5 and between the two tail chains are set for the last one bead in head chain and the two first beads in each tail chain. We draws the schematic diagram for the phospholipid CG model in two-dimensional plane, as shown in Fig. 1. This model with multiple types of phospholipid membrane is similar to the symmetric membrane models used in previous studies.⁵⁰

Fig. 1 The schematic diagram for the phospholipid molecules with head group and two tail groups, as well as the asymmetric structure. The first type of phospholipid is demonstrated in the left, where the reds denote the head beads and blues are the tail beads. The second type of phospholipid is in the right, where the pinks represent the head beads and the yellows are the tail beads. The asymmetric membrane is taken as an example for the asymmetric structure and modeled by these types of phospholipid molecules is shown in the center.

Parameters

The simulations were performed in a cubic box with volume V = L × L × L under periodic boundary conditions. For simplification, we assumed the phospholipid and solvent water beads to have the same mass, and set this mass as the mass scale. We set the interaction cut-off radius for the phospholipid molecule and solvent beads r_c as the length scale, which is suitable for the phospholipid model. The cut-off radius can be estimated according to the expresses, r_c = (ρV_b)^1/3, where ρ is bead density and V_b is the volume of one DPD bead. We set ρ = 3 in our simulation, similar to previous simulations.⁵¹ The unit of time τ is defined aswhere k_BT is the energy scale. In the current simulations, we set one DPD time step, Δt = 0.01τ.

We took the repulsive interaction parameters a_ii = 25 for the same type of beads and a_ii = 100 for different types of beads; this value has been used extensively in previous simulations.^11,48,52 Then, the Flory–Huggins parameter χ, was estimated by the relationship between a_ii and a_ij, i.e.,χ = 0.286(a_ij − a_ii), which has been adopted in previous simulations.^53–55 To avoid the finite-size effect, we optimized the simulation box sizes by varying the box size L from 25r_c to 35r_c. We selected the one with the lowest energy as the most optimized box size (please see an example in Fig. S1 and S2 of the ESI†). All simulations were performed on NVT ensembles, and the numbers of chains for the first and second types of phospholipids were set to n₁ = n₂ = 600. Usually, the system can reach the equilibrium state when the process is about 200 000 DPD time steps, as shown in Fig. 2. At this state, the energy decreases to the lowest value. Moreover, pressure and temperature were uniformly distributed in the simulation boxes, similar to a previous simulation.¹⁶ An exmaple data is shown in Fig. S3 of the ESI,† where the pressure and temperature distribute uniformly in the space. For the same system with fixed parameters, we inputted different initial distributions of phospholipid chains, which led to different output energies in the system. Then, we selected the one with the lowest energy as the final stable structure for this parameter point. For example, we inputted three initial chain distributions, i.e., lamellar, cylindrical, and spherical, and obtained three output energies in one phase parameter point, as shown in Fig. 2. After comparing the output energies, we selected the one with the lowest energy as the stable structure.

Fig. 2 An example for obtaining the equilibrium and stable state in the dynamic processes, where the parameters are set to be N_HB1 = 3, N_TB1 = 3, N_HB2 = 3, and N_TB2 = 4. Three initial inputting distributions, i.e., the lamella, cylinder and sphere inputting, are considered, respectively.

Results and discussions

For a system with two types of phospholipid chains in water solutions, the parameter space is too large to explore. Hence, we reduced the globular space into a simplified parameter space by fixing several parameters. In particular, we fixed the interaction parameter (a_ij), numbers of phospholipid chains (n₁ and n₂), and head bead numbers (N_HB1 = N_HB2 = 3). Meanwhile, we varied the tail bead numbers for the two types of phospholipid chains (N_TB1 and N_TB2) from 2 to 10. We concentrated on the observed asymmetric structures (Fig. 3–5). The corresponding dynamic processes (Fig. 6–8) and mechanical properties (Fig. 9–11) were derived by varying the tail bead numbers.

Fig. 3 Typical asymmetric membrane with N_HB1 = 3, N_TB1 = 4, N_HB2 = 3, and N_TB2 = 9. Snapshot for the phospholipid bilayer with (a) the top view and (b) side view are demonstrated, and (c) the density profiles of head and tail beads and (d) the order parameters are shown along the z axes. The dotted line indicates the fitting line.

Fig. 4 Typical asymmetric tube and vesicle. (a1–a3) tube with N_HB1 = 3, N_TB1 = 3, N_HB2 = 3, and N_TB2 = 3, and (b1–b3) vesicle with N_HB1 = 3, N_TB1 = 2, N_HB2 = 3, and N_TB2 = 2. The top view and side view are list while the density profiles of head and tail beads are list in the right.

Fig. 5 Characteristic for the asymmetric structures depending on the chain length. A phase diagram with N_HB1 = 3 and N_HB2 = 3, in terms of the tail chain length, N_TB1and N_TB2. The symbols , , represent the asymmetric vesicle, asymmetric tube and asymmetric membrane, respectively.

Fig. 6 The dynamic processes for the typical (a) asymmetric membrane with N_TB1 = 4 and N_TB2 = 9, (b) asymmetric tube with N_TB1 = 3 and N_TB2 = 3, and (c) asymmetric vesicle N_TB1 = 2 and N_TB2 = 2. The typical structures are shown in three stages with the average energies and duration times list below.

Fig. 7 The average radius of gyration of the asymmetric membrane with N_HB1 = 3, N_TB1 = 4, N_HB2 = 3, and N_TB2 = 9. Three components for (a) the first type of phospholipid chain and (b) the second type of phospholipid chain as functions of time steps.

Fig. 8 The shape factor of the asymmetric vesicle with N_HB1 = 3, N_TB1 = 2, N_HB2 = 3, and N_TB2 = 2. The data circled by the red dotted box is enlarged and shown in the inserted.

Fig. 9 The interface tension 〈σ_z〉 as function of distance along x-direction for the asymmetric membrane, with N_HB1 = 3, N_TB1 = 3, N_HB2 = 3, and N_TB2 = 10. The asymmetric membranes are also inserted.

Fig. 10 The interface tension 〈σ_z〉 as function of distance along z-direction for the asymmetric membrane, with N_HB1 = 3, N_TB1 = 3, N_HB2 = 3, and N_TB2 = 7. (a) The first type of phospholipid chain, (b) the second type of phospholipid chain. The asymmetric membranes are also inserted.

Fig. 11 Osmotic pressure as function of distance along z-direction for the asymmetric membrane with N_HB1 = 3, N_TB1 = 3, N_HB2 = 3, and N_TB2 = 7. (a) The first type of phospholipid chain, (b) the second type of phospholipid chain.

Asymmetric structures and phase diagrams

This subsection discusses three types of asymmetric structures, i.e., asymmetric membrane, asymmetric tube, and asymmetric vesicle. Fig. 3 shows a typical asymmetrical membrane structure with parameters of N_TB1 = 4 and N_TB2 = 9. This asymmetric bilayer membrane constitutes two types of monolayer membranes (top and side views shown in Fig. 3a and b, respectively). Naturally, the hydrophilic head beads are arranged near the outer water regions, and the hydrophobic tail beads are packed into the center regions inside the membranes due to the amphiphilicity of phospholipid molecules in water solutions. To intuitively describe the asymmetric membrane structure, we plotted the density distributions of the hydrophilic and hydrophobic beads along the z-direction, as shown in Fig. 3c. The density profiles for the head and tail bead distributions have four peaks. The peaks in the head bead profiles are separated for the first and second types of phospholipids, which clearly demonstrates the bilayer structures in the membranes. Meanwhile, the peaks overlap in the tail bead profiles, which indicates that the tail chains are interdigitated for the two types of phospholipids. Such a bimodal structure has also been confirmed in previous phospholipid membrane simulations.^56–58 In addition, we investigated the orientation ordering of the head beads for asymmetric phospholipid membranes as follows:^59,60where θ is the angle between the chain direction and the z-axis direction and the bracket represents an ensemble average. In general, when the direction of the chain is completely parallel to the z-axis, the value of the degree of order is equal to 1. When the direction of the chain is completely perpendicular to the z-axis, the value of the degree of order approaches −0.5.¹¹ The orientation distributions are shown in Fig. 3d, where the red dot represents the simulation results of the first type of head beads, the pink dot denotes the simulation results of the second type of head beads, and the dotted lines are the fitting results. The fitting curve is close to the Gaussian distribution, and the simulation value is between 0 and 0.5, but not 0. These facts indicate that the angles between the chain length and z-axis are between 0 and π, thus showing that the asymmetrical phospholipid membranes have well-ordered structures. However, the order degree of head bead is bigger for the second type of phospholipid than those of the first type phospholipid. Although the numbers of head beads are the same for these two type of phospholipids, they have different tail chains. The longer tail chains for the second type phospholipid, i.e., N_TB2 = 9, have more obvious effects on the ordering degree of the head chain than those with the short tail chain of N_TB1 = 4.

The two other asymmetric structures, namely, asymmetric tube and asymmetric vesicle, are shown in Fig. 4. The top and side views and bead density distributions of the two typical asymmetric structures are presented. The side and spatial cross-sectional views of the asymmetric tube structure, with parameters of N_TB1 = 3 and N_TB2 = 3, are shown in Fig. 4a1 and 4a2, respectively. From these cross sections, we observed that the cylindrical tube has an asymmetric structure, with the first type of beads distributed in the inner layer. To analyze the spatial structure of the asymmetric tube, we plotted the density profile of the beads along the z-axis cross-section in Fig. 4a3. For the inner layer, the double peaks appear in the profiles of head and tail bead distributions, i.e., ϕ_HB1 and ϕ_TB1, which indicates that the asymmetric column is hollow. The peak of tail bead profile ϕ_TB1 for the inner layer is close to the peak of the tail bead profile in outer layer ϕ_TB2. This fact indicates that the chains in the two layers overlap. For the asymmetric vesicle, we plotted front and spatial cross-sectional views, as shown in Fig. 4b1 and 4b2, respectively. Similarly, the spherical vesicle has a bilayer structure, with the first type of chains in the inner regions and the second type of chains in the outer regions forming a hollow spherical-like structure. This distribution is due to the hydrophilic heads and hydrophobic tails of the phospholipid molecules in water solutions. Then, we analyzed the asymmetric spherical structure by plotting the bead density distributions in the y-axis cross-sectional direction, as shown in Fig. 4b3. Although the double peaks in these bead profiles are not obvious, we could still distinguish the bilayer structures from the obvious peaks near the outer regions, and the phospholipid molecules are distributed in a certain order. Our observations on the asymmetric vesicles and membranes are similar to those observed for planar and vesicle membranes in multiple types of phospholipids, among which POPC, POPE, POPS, PPCS, CHOL, and fatty acid molecules with fixed chain lengths participate in assembly.³⁴ However, we observed the cylindrical membranes, i.e., the tubes, for only two types of phospholipids in water solutions. We focused on the structure distributions that depended on chain lengths for the two types of phospholipids.

To observe such a chain length effect, we changed the tail chains for the two types of phospholipids and constructed the phase diagrams, as shown in Fig. 5. Here, we assume only three candidate phases, the asymmetric membrane, asymmetric tube and asymmetric vesicle, appearing in the phase diagram. In particular, we changed the number of tail beads in the simulation to range from N_TB1 = 2–10 and N_TB2 = 2–10. Given that the number of beads in the two tail chains is symmetrical, only half of the phase diagrams are plotted. In the phase diagrams, three phase regions are divided with a dotted line for visual guidance, and the asymmetric membranes, asymmetric tubes, and asymmetric vesicles are separated. The phase distributions in the phase diagram have several characteristics. The first one is that the asymmetric vesicles are located in the bottom of the phase diagram where the tail chains are shorter for one type of phospholipid molecules, mostly even with N_TB1 = 1, as shown in Fig. 5. Usually, spherical structures have high curvatures, which need high asymmetric components in the inner and outer regions. This condition requires small numbers of tail beads in the outer regions and large numbers of tail beads in the inner regions. The second characteristic is that the asymmetric tube is distributed near the diagonal regions where the tail chains are nearly equal to the two types of phospholipids. This distribution indicates that the two types of tail chains have similar features and are inputted into the center regions, whereas the head chains are attracted to the outermost and innermost regions by the water molecules. The third characteristic is that the asymmetric membranes appear in regions where one tail chain is long, whereas the other chains are short. This condition indicates that the relatively long tail chains can easily form such structures. This fact guides us in obtaining asymmetric membranes for multiple types of phospholipid molecules in water solutions. In additional, we have studied the change in the average radius of gyration of the structure with N_TB1 = 3 along the dotted line in Fig. 5 (see Fig. S4 in the ESI†). Here, the average radius of gyration will be defined in eqn (7) and (8) in the next subsection. The results indicated that the average radius of gyration initially decreases as N_TB2 increases, where the structure changes from an asymmetric tube to an asymmetric vesicle. Then, the average radius of gyration smoothly increases as N_TB2 increases when the structures are always asymmetric membranes.

Dynamic process

We investigated several typical dynamic processes of asymmetric structures, namely, asymmetric membrane, asymmetric tube, and asymmetric vesicle, with a time step. Three examples are shown in Fig. 6, where the main stages are listed with evolutionary states and time and the corresponding average energies. First, all three dynamic processes of asymmetric structures can be divided into three stages: random generation, mutual adaptation, and formation. Specifically, for asymmetric membrane the initial random stage lasts for 26τ at a relatively high energy of 7.95k_BT and subsequently evolves into the mutual adaption stage, which sustains a relatively long time of 374τ at a low average energy of 6.53k_BT, as shown in Fig. 6a. In this stage, the structure is continuously adjusted and evolved. Afterward, the system develops into the structure-formatting stage, where the relatively regular membranes are completely formatted with the lowest energy of 6.09k_BT. The energy is gradually reduced and tends to stabilize during the formation of the asymmetric membrane, which indicates that the asymmetric membrane structure in the formation stage has reached a stable structure. This evolution process is similar to that in previous simulations.^11,56 In the dynamic processes for the asymmetric tube and vesicle, as shown in Fig. 6b and c, the durations and energies differ from those observed in the asymmetric membranes. The energy changes in the three stages are 6.58k_BT, 5.96k_BT, and 5.94k_BT, which tend to stabilize the structure as the time step increases, as shown in Fig. 6b. The random generation phase of the asymmetric tube structure is longer than that of the asymmetric membrane structure, but the asymmetric tube structure can enter the formation phase more quickly. This result suggests that the process can rapidly drown in the polymer fiber where the polymer chains are highly stretched.¹⁷ Our observation on the asymmetric tube agrees with this prediction. For the dynamic process of the asymmetric vesicle structure, as shown in Fig. 6c, the structure occludes the upper and lower parts, gradually joins the connection over time, and finally forms a hollow asymmetric spherical structure with energy changes of 6.16k_BT, 5.90k_BT, and 5.89k_BT. The asymmetric spherical structure requires a longer simulation time to achieve a stable structure compared with the asymmetric membrane and asymmetric cylinder.

To describe the dynamic process in detail, we considered the average radius of gyration 〈R_g〉. The radius of gyration tensor R²_g can be expressed as⁶¹

The elements R²_gαβ can be written as

with α,β∈{x,y,z};N and r_i,x are the number of chains and x-coordinate of i-th bead, respectively; and r_c is the mass center. In the current simulations, we focused on the average radius of gyration of the asymmetric membrane, as shown in Fig. 7, where the parameters are N_TB1 = 4 and N_TB2 = 9. The average radius of gyration for the two types of phospholipid molecules are shown in Fig. 7a and b. The average value of 〈R_gzz〉 is smaller than those of 〈R_gxx〉 and 〈R_gyy〉, and the average values of 〈R_gxx〉 and 〈R_gyy〉 are basically the same, indicating that the polymer chains are compact in the z-direction and stretched in the x- and y-directions, as shown in Fig. 7a. During the self-assembly period, the average values of 〈R_gxx〉 and 〈R_gyy〉 increase initially then decrease to stable values, as shown in Fig. 6a. This result is consistent with the microstructural dynamics observed in Fig. 6a, which is a layered distribution on the x–y plane. For the second type of phospholipids, which is unlike the first type, the average radius of gyration in the z direction 〈R_gzz〉 is almost the same as those in the x- and y-directions, 〈R_gxx〉 and 〈R_gyy〉. The values of 〈R_gxx〉 and 〈R_gyy〉 in the x- and y-directions are basically the same, and they smoothly decrease to stable values in the structural formation stage. As observed in Fig. 6a, the second type of phospholipids is distributed in the above layers, most parts of which complete the chain arrangement in the last two stages. A small adjustment was observed in the distances between the chain beads, where the distance between the beads decreases and becomes compact. This result is similar to the formation process we observed previously in phospholipid symmetric membranes.¹⁶

Furthermore, we extracted the shape factor of the structure by using the average radius of gyration and investigated its evolution in the dynamic processes. Shape factor 〈δ〉 can be expressed as follows:^62,63

where L₁², L₂², and L₃² are the three eigenvalues of the square radius of gyration tensor R_g². For convenience, we sorted them as L₁² ≤ L₂² ≤ L₃². We investigated the shape factor for the chains in the asymmetric vesicle. An example is shown in Fig. 8, where the parameters are set to N_TB1 = 2 and N_TB2 = 2. The shape factor increases then stabilizes at 0.5, as shown in Fig. 8, where the chain conformations are also plotted. In the beginning, the chain exhibits an elliptical conformation then a circular conformation at the stable stage. This is evidence that shape factor 〈δ〉 begins at 0.47 and then approaches 0.5 at the stable stage, a result that is similar to previous simulation results.⁶⁴ To observe the shape factor in detail, we drew the variation of the shape factor between 0τ and 500τ as the inserted part. By investigating the variation of the shape factor, we found that the asymmetric vesicle structure easily forms when the structure of the single chain exhibits spherical conformation.

Mechanical properties

We focused on the mechanical properties of asymmetric membrane structures in solutions by analyzing the interface tension. The interface tension of the phospholipid symmetric membrane has elicited much interest in recent years.^25,65,66 According to the Irving–Kirkwood definition, the formulation of tension σ_z along the z direction is as follows:⁶⁷where the element p_xx is a component of the pressure tensor and can be achieved aswhere N_p is the number of DPD particles and F_ijx and x_ij are the force and relative position between particles i and j along the x-direction. Elements p_yy and p_zz have the same definition, except for the replacement of the corresponding subscripts. Fig. 9 shows the surface tension along the x-direction of an asymmetric membrane structure with parameters N_TB1 = 3 and N_TB2 = 10, where surface tensions σ_z are averaged over the y- and z-directions. The average value of σ_z is nearly equal to zero, indicating that the asymmetric membrane is a free and planar membrane in the solutions despite of x changes. This result is similar to those of previous theoretical and simulation observations of symmetric membranes for phospholipids in solutions.^11,68

Given that the interface tensions are equal to zero by taking the average over y- and z-directions, the detailed distribution in the z-direction is still necessary. We investigated the mechanical properties for the asymmetric membrane by calculating the interface tension distributed along the z-direction for the two types of phospholipid molecules, as shown in Fig. 10. The example parameters were set to N_HB1 = 3 and N_TB1 = 3 for the first type of phospholipids and to N_HB2 = 3 and N_TB2 = 7 for the second type. As shown in Fig. 10a, first-type phospholipid molecules are predominantly distributed in a range of 10.7–13.9r_c, leading to the distribution of interface tension in the same region. For the head beads, the interface tension decreases in the beginning then increases, namely, the interface tension is distributed in the interfacial regions between the mixing molecules. However, the interface tension only exists near the regions between the water and phospholipid molecules for the tail beads. This condition means that the tensions in the interface regions between the head and tail chains are mainly contributed by the parts of the head chains because they are more rigid than those of the tail chains. The second type of phospholipids, which is mainly distributed in the regions from 13.5r_c to 16.8r_c, makes up the other asymmetric layer, as shown in Fig. 10b. Similarly, we plotted the interface tensions at the tail and head chains for the second type of phospholipids, as shown in Fig. 10b. The head chains have similar distributions as those in the first type due to the rigidity of the head chains, but the tension of the tail chains only fluctuates within a small value. This result was obtained probably because the second type has more flexibility than the first type; the tail chains are longer for the second type than for the first type. Comparison of the interface tensions of the two types of phospholipids shows that the interior of the asymmetric membrane structure is essentially in a tension-free state, but interfacial tension exists near the boundary.

In addition to surface tension, the permeability of the asymmetric membrane structure is also important. The formula for calculating osmotic pressure is as follows:

Here, we plotted the osmotic pressures as functions of position z for the asymmetric membrane, as shown in Fig. 11, where the parameters are set to N_TB1 = 3 and N_TB2 = 7, the same as those in Fig. 10. For the first type of phospholipids, the two peaks appear in the two interfacial regions between the head beads and water beads, head beads and tail beads, respectively, as shown in Fig. 11a. However, the tail beads mainly contribute one maximum value near the tail beads from the second type of phospholipid. Actually, the osmotic pressure results from the gradient distribution of interfacial tension, it obviously has the maximum value where the interfacial tension varies most intensively, according to the distributions in Fig. 10a. Similar cases were observed for the second type of phospholipids, as shown in Fig. 11b. The difference is that the osmotic pressure in the center regions results from the head chains not the tail chains. Namely, the tail beads of the second type phospholipid contribute a very small percentage for the osmotic pressure, which is probably due to its short length. This result might help us understand the non-directional transmission through asymmetric membranes.

Summary

In this work, we used DPD simulation to predict an asymmetric membrane, asymmetric tube, and asymmetric vesicle for phospholipids in aqueous solutions. The asymmetric structures assembled from two types of phospholipid molecules were investigated by analyzing their bead density distributions and order parameters. The results showed double layers in these asymmetric structures and good ordering near the boundaries. The asymmetric structures distributed in the phase diagram were also analyzed by changing the tail chain lengths of the two types of phospholipids. The phase diagram indicated that the asymmetric tubes were distributed along the diagonal regions, and the asymmetric membranes were located in the regions with long tail chains due to the chain asymmetry and amphiphilicities in the solutions. However, the phase diagram in the current work is limited in the variances of tail chain lengths. We expect a richer phase behavior when the parameters of the system expand beyond those considered here. For example, we assumed comparable interactions between the beads, i.e., the Flory–Huggins parameter. When the Flory–Huggins parameter become stronger, we expect that the membrane phase can be extended in the phase diagram. In particular, our observation on the subset of the full parameter space is a starting point to explore the phase behavior of asymmetric structures in the full parameter space.

Then, we investigated the dynamics of the asymmetric structures and observed the formation of the asymmetric membrane, asymmetric tube, and asymmetric vesicle in the solutions. The structural evolution was divided into three stages for the three asymmetric structures. The stages were initial random, adaptation, and formation, which were determined by the energy of the structure. Furthermore, the average radius of gyration of the asymmetric membrane structure was investigated. The dynamic formation process of the asymmetric membrane structure was analyzed, and the formation process of the membrane was evaluated intuitively. The shape factor of the asymmetric vesicle structure proved that the self-assembly of the conformation phospholipid molecules formed a spherical conformation.

We also investigated the mechanical properties of the asymmetric membrane in solutions by analyzing the interface tension and osmotic pressure. Evaluation of the mechanical properties of the asymmetric membrane structure revealed that the average tension for the membrane structure was zero, but interface tension and osmotic pressure were distributed along the z-direction. The simulation suggested that tension and pressure existed near the boundaries between the phospholipid and water domains due to the rigidities of the head chains. This work provides a method of designing asymmetric membranes in biological processes and understanding the mechanism of directional transport via membranes.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was supported by the program of National Natural Science Foundation of China (Grant No. 21973070, 21474076, 21674082, and 11875205).

Notes and references

1. H.-C. Yang, J. Hou, V. Chen and Z.-K. Xu, Angew. Chem., Int. Ed., 2016, 55, 13398–13407.
2. H. I. Okur, O. B. Tarun and S. Roke, J. Am. Chem. Soc., 2019, 141, 12168–12181.
3. M. Aron, O. Vince, M. Gray, C. Mannaris and E. Stride, Langmuir, 2019, 35, 13205–13215.
4. H. Zhou and Z. Guo, J. Mater. Chem. A, 2019, 7, 12921–12950.
5. Y. N. Liu, R. Qu, W. Zhang, X. Li, Y. Wei and L. Feng, ACS Appl. Mater. Interfaces, 2019, 11, 20545–20556.
6. C. V. Kulkarni, Nanoscale, 2012, 4, 5779–5791.
7. H. Qiu and M. Caffrey, Biomaterials, 2000, 21, 223–234.
8. C. V. Kulkarni, T.-Y. Tang, A. M. Seddon, J. M. Seddon, O. Ces and R. H. Templer, Soft Matter, 2010, 6, 3191–3194.
9. Z. Wang, X. Wang, Y. Ji, X. Qiang, L. He and S. Li, Polymer, 2018, 140, 304–314.
10. A. Saitov, S. A. Akimov, T. R. Galimzyanov, T. Glasnov and P. Pohl, Phys. Rev. Lett., 2020, 124, 108102.
11. X. Qiang, X. Wang, Y. Ji, S. Li and L. He, Polymer, 2017, 115, 1–11.
12. S. Ma, Y. Hu and R. Wang, Macromolecules, 2016, 49, 3535–3541.
13. H. Tan, C. Yu, Z. Lu, Y. Zhou and D. Yan, Soft Matter, 2017, 13, 6178–6188.
14. M. Podewitz, Y. Wang, P. Gkeka, S. Grafenstein, K. R. Liedl and Z. Cournia, J. Phys. Chem. B, 2018, 122, 10505–10521.
15. M. Hishida, K. Tanaka, Y. Yamamura and K. Saito, J. Phys. Soc. Jpn., 2014, 83, 044801.
16. Y. Shan, X. Wang, Y. Ji, L. He and S. Li, J. Chem. Phys., 2018, 149, 244901.
17. S. Chandran, J. Baschnagel, D. Cangialosi, K. Fukao, E. Glynos, L. M. C. Janssen, M. Müller, M. Muthukumar, U. Steiner, J. Xu, S. Napolitano and G. Reiter, Macromolecules, 2019, 52, 7146–7156.
18. B. Mostofian, T. Zhuang, X. Cheng and J. D. Nickels, J. Phys. Chem. B, 2019, 123, 5814–5821.
19. R.-H. Teresa, G. F. Thomas and L. Mahadevan, Phys. Rev. Lett., 2019, 123, 038102.
20. S. Gao, Z. Kang, R. Yuan, N. Liu, P. Zhu and B. Wang, J. Mol. Struct., 2019, 1192, 35–41.
21. G. Surmeier, M. Paulus, P. Salmen, S. Dogan, C. Sternemann and J. Nase, Biophys. Chem., 2019, 252, 106210.
22. A. Mescola, N. M. Medina, G. Ragazzini, M. Accolla and A. Alessandrini, J. Colloid Interface Sci., 2019, 553, 247–258.
23. M. Manisegaran, S. Bornemann, I. Kiesel and R. Winter, Phys. Chem. Chem. Phys., 2019, 21, 18533–18540.
24. C. L. Ting and M. Müller, J. Chem. Phys., 2017, 146, 104901.
25. L. Guillaume, A. Antoine, P. Juan, L. Sid, L. Martin and C. Clément, Phys. Rev. X, 2020, 10, 011031.
26. S. Jiang, Q. Chen, M. Tripathy, E. Luijten, K. S. Schweizer and S. Granick, Adv. Mater., 2010, 22, 1060–1071.
27. X.-C. Luu, J. Yu and A. Striolo, J. Phys. Chem. B, 2013, 117, 13922–13929.
28. A. Perro, S. Reculusa, S. Ravaine, E. Bourgeat-Lami and E. Duguet, J. Mater. Chem., 2005, 15, 3745–3760.
29. A. Walther and A. H. E. Müller, Soft Matter, 2008, 4, 663–668.
30. L. Zhang, Y. Chen, Z. Li, L. Li, P. Saint-Cricq, C. Li, J. Lin, C. Wang, Z. Su and J. I. Zink, Angew. Chem., Int. Ed., 2016, 55, 2118–2121.
31. M. Nau, N. Herzog, J. Schmidt, T. Meckel, A. Andrieu-Brunsen and M. Biesalski, Adv. Mater. Interfaces, 2019, 6, 1900892.
32. X. Yang, L. Yan, F. Ran, A. Pal, J. Long and L. Shao, J. Membr. Sci., 2019, 576, 9–16.
33. Y.-L. Yang, Y.-J. Sheng and H.-K. Tsao, Phys. Chem. Chem. Phys., 2018, 20, 27305–27313.
34. S. Sharma, B. N. Kim, P. J. Stansfeld, M. S. P. Sansom and M. Lindau, PLoS One, 2015, 10, 1–21.
35. H. I. Ingólfsson, M. N. Melo, F. J. Eerden, C. Arnarez, C. A. Lopez, T. A. Wassenaar, X. Periole, A. H. Vries, D. P. Tieleman and S. J. Marrink, J. Am. Chem. Soc., 2014, 136, 14554–14559.
36. R. Gupta, B. S. Dwadasi and B. Rai, J. Phys. Chem. B, 2016, 120, 12536–12546.
37. M. Paloncyova, K. Vavrova, Z. Sovova, R. H. DeVane, M. Otyepka and K. Berka, J. Phys. Chem. B, 2015, 119, 9811–9819.
38. R.-X. Gu, S. Baoukina and D. P. Tieleman, J. Am. Chem. Soc., 2020, 142, 2844–2856.
39. X. Yue, T. Zhang, D. Yang, F. Qiu, G. Wei and H. Zhou, Nano Energy, 2019, 63, 103808.
40. H.-C. Yang, Y. Xie, J. Hou, A. K. Cheetham, V. Chen and S. B. Darling, Adv. Mater., 2018, 30, 1801495.
41. T.-D. Pan, Z.-J. Li, D.-H. Shou, W. Shou, J.-T. Fan, X. Liu and Y. Liu, Adv. Mater. Interfaces, 2019, 1901130.
42. N. Tran, A. M. Hawley, J. Zhai, B. W. Muir, C. Fong, C. J. Drummond and X. Mulet, Langmuir, 2016, 32, 4509–4520.
43. B. D. Todd and P. J. Daivis, Mol. Simul., 2007, 33, 189–229.
44. X. Song, H. Guo, J. Tao, S. Zhao, X. Han and H. Liu, Chem. Eng. Sci., 2018, 189, 75–83.
45. P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett., 1992, 19, 155–160.
46. P. Español and P. B. Warren, J. Chem. Phys., 2017, 146, 150901.
47. L. Zhang, M. Becton and X. Wang, J. Phys. Chem. B, 2015, 119, 3786–3794.
48. H.-M. Ding and Y.-Q. Ma, Small, 2015, 11, 1055–1071.
49. M. Venturoli, B. Smit and M. M. Sperotto, Biophys. J., 2005, 88, 1778–1798.
50. M. Cebecauer, M. Amaro, P. Jurkiewicz, M. J. Sarmento, R. Šachl, L. Cwiklik and M. Hof, Chem. Rev., 2018, 118, 11259–11297.
51. X. Li, Y. Liu, L. Wang, M. Deng and H. Liang, Phys. Chem. Chem. Phys., 2009, 11, 4051–4059.
52. H.-M. Ding and Y.-Q. Ma, Biomaterials, 2013, 34, 8401–8407.
53. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
54. A. Maiti and S. McGrother, J. Chem. Phys., 2004, 120, 1594–1601.
55. A. S. Özen, U. Sen and C. Atilgan, J. Chem. Phys., 2006, 124, 064905.
56. L. He, L. Zhang, Y. Ye and H. Liang, J. Phys. Chem. B, 2010, 114, 7189–7200.
57. M.-J. Huang, R. Kapral, A. S. Mikhailov and H.-Y. Chen, J. Chem. Phys., 2012, 137, 055101.
58. G. J. A. Sevink and J. G. E. M. Fraaije, Soft Matter, 2014, 10, 5129–5146.
59. Å. A. Skjevik, B. D. Madej, R. C. Walker and K. Teigen, J. Phys. Chem. B, 2012, 116, 11124–11136.
60. K. R. Hadley and C. MCabe, Soft Matter, 2012, 8, 4802–4814.
61. W. Huang and V. Zaburdaev, Soft Matter, 2018, 15, 1785–1792.
62. H. W. Diehl and E. Eisenriegler, J. Phys. A: Math. Gen., 1989, 22, L87–L91.
63. G. Zifferer and W. Preusser, Macromol. Theory Simul., 2001, 10, 397–407.
64. Y. Shan, X. Qiang, J. Ye, X. Wang, L. He and S. Li, Sci. Rep., 2019, 9, 15393.
65. A. D. Petelska, Cent. Eur. J. Chem., 2012, 10, 16–26.
66. T. Takei, T. Yaguchi, T. Fujii, T. Nomoto, T. Toyota and M. Fujinami, Soft Matter, 2015, 11, 8641–8647.
67. J. H. Irving and J. G. Kirkwood, J. Chem. Phys., 1950, 18, 817–829.
68. T. Yamamoto and S. A. Safran, Soft Matter, 2011, 7, 7021–7033.

---

Footnote

† Electronic supplementary information (ESI) available: Examples for the finite box size effect, distribution of pressure and temperature, and average radius of gyration with varying tail beads. See DOI: 10.1039/d0ra03732j

---