Cholesterol affected dynamics of lipids in tailor-made vesicles by ArcVes software during multi micro second coarse grained molecular dynamics simulations

1.   Introduction

Vesicles and liposomes represent spherical objects comprised of a lipid membrane. In biotechnology, they are important structural features in both living cells and artificial systems [1]–[3]. In the former case, vesicle-like systems are constantly generated and depleted by budding and fusion, respectively, and are part of the endosomal [4] and lysosomal pathways [5]. Within the cellular life cycle, vesicular systems harbor membrane proteins. Within the interior of the vesicle, essential metabolic biochemical reactions occur, such as post processing folding, metabolic degradation, and virus assembly [6]. These cellular tasks have triggered the development of manmade analogs [7] to be used for cargo delivery [8]–[10], including the delivery of either membrane proteins or soluble products to target cells, as well as biocontainment mechanisms for (bio)chemical reactions [11],[12].

One can mimics the dynamics of such vesicles by applying structure based molecular models of lipids, proteins and the surrounding water molecules. Some of the aforementioned cellular features of vesicles have already been modeled [13], such as virus encapsulation [14], lipid dynamics in virus-like liposomes [15], membrane disruption [16], fusion [17],[18], membrane protein diffusion [19] and their function as bio containers [20]. On the experimental side, a series of precise mixtures and technical properties [7],[21]–[25] need to be screened and assessed for technical applications regarding stability, as well as lipid and compound molecular dynamics [26]. Thus, there is a high demand for a computational method, which allows for the formulation of tailor-made lipid vesicles.

In many of the existing computational methods, the lipid molecules are set by operators. Several methods have been used to generate vesicular systems using coarse-grained (CG) lipid models and applying molecular dynamic (MD) simulations for further analysis. One of the methods is to place lipid molecules into an aqueous solution and let them self-assemble due to their amphipathic character [27]–[29]. Another method is to place lipid molecules along the inner and outer radii [30]–[32] and allowing for flip-flopping using transient pores to avoid frustration of the vesicle [33],[34]. Especially in the first case, it is difficult to create lipid mixtures that also align across the leaflets. More recently, a step-wise build-up method, which uses a similar strategy to generate vesicles, was developed independently from this study [35].

In the presented software, Arc-corrected Vesicle construction (ArcVes), parameters of lipids known from planar systems, such as the area-per-lipid (apl) and membrane thickness, are mapped onto the surface of a sphere to produce vesicles with customized and well-defined lipid numbers, composition and eventually embedded proteins. Albeit shown here for CG models, the method can also be applied for united atom model systems and be applied to any force field (ff).

The set-up can be used to explore the possibility of ff methods to quantify dynamic properties. In the present study, the effect of cholesterol on lipid membranes in a vesicle-like geometry are explored based on key values such as radius, membrane thickness, lipid flipping and diffusivity.

2.   Materials and methods

2.1.   Generation and simulation of planar lipid patches for CG models

The planar lipid patches were generated by using the CHARMM-GUI web server (Chemistry at HARvard Molecular Mechanics – Graphical User Interface) [36]–[38]. Each planar patch contains 125 lipid molecules per leaflet and around 5,200 to 5,500 water molecules, hence the name “wet” patches. Five different patches were generated, including two single-lipid patches and three mixed-lipid patches. The two single-lipid patches were a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) patch and a 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) patch. The mixed-lipid patches were made of both POPC and DOPC (1: 1) (POPC/DOPC), POPC and cholesterol (CHOL) (7:3, POPC/CHOL), and DOPC-CHOL (7:3, DOPC/CHOL). The lipid composition was chosen to satisfy (i) the abundance of POPC in lipids and plasma membranes [39], (ii) the presentation of a simple dual mixture and (iii) the occurrence of 30% CHOL in both leaflets of epithelia cells [40]. All patches were duplicated, and water molecules and additional ion atoms were removed in the duplicated patches to be used as the dry-MARTINI (nickname of the city of Groningen, NL) system. Since they did not contain any water molecules, these duplicated dry patches were energy minimized for 10,000 steps by the steepest descent with the dry MARTINI v2.1 ff [41]. The original “wet” patches were operated through the MARTINI v2.2 ff [42] for 10,000 steps steepest descent minimization and additional 5 steps equilibration runs. The equilibration runs were operated by releasing the lipid restraints force from 200 to 10 kJ mol⁻¹ nm⁻². The restraints force was completely released to 0 kJ mol⁻¹ nm⁻² in the production run.

The production runs were 10 µs, with time steps of 30 fs for both dry and wet patches. The cut off distances for the Coulomb and van der Waals force were 1.2 nm for dry patches and 1.1 nm for wet patches. The dry patches were run under the NVT condition (canonical ensemble with fixed number of atoms, N, volume, V, and temperature, T) with system pressure equals to 1 bar. The wet patches were run under the NPT condition (isobaric-isothermal ensemble with fixed number of atoms, N, pressure, P, and temperature, T) with a system pressure of 1 bar, which was coupled by the Parrinello-Rahman semi-isotropic coupling scheme (x-y directions). The system temperature was 323 K, coupled by Berendsen coupling for dry patches, and 320 K coupled by v-rescale coupling for wet patches. The equilibration runs and production runs were performed on the Gromacs 2021 platform [43].

2.2.   Generation of the vesicles using newly developed Arc-corrected Vesicle construction (ArcVes) software

A single coarse-grained representative for each of the lipid molecules, POPC, DOPC, and CHOL [44] was obtained from the Martini General Purpose Coarse-Grained Force Field (http://cgmartini.nl/index.php, last used June 2021) (see Suppl. Figure 1 for naming of the spheres and their individual parameters). The hydrophobic tails of the downloaded lipid molecules, POPC and DOPC, were straight structures and the tails parallel to each other. These lipid molecules were used as the building blocks in the ArcVes software (Suppl. Figure 1, last used November 2022) to generate vesicles.

The software allows the following parameters to be entered in the following sequence: lipid ratio, radius, density of each leaflet providing the apl values, thickness of the bilayer, and the position of the lipids other than the default value (Figure 1). In addition, the protein can be added (more than 3 protein molecules). Finally, there is the option to open another window and customize the number of lipid molecules in each leaflet by individually changing the numbers for each leaflet. The software can be downloaded for free including the Windows-, Mac- and Linux-versions, https://wfischer.web.nycu.edu.tw/downloads-2).

Figure 1.  Arc-corrected Vesicle construction (ArcVes): The two entry windows to generate tailor-made vesicle by placing the lipid molecule and proteins as well as ‘customizing’ the distributing of the lipid molecules in each leaflet.

DownLoad: Full-Size Img PowerPoint

The vesicles were generated by defining the radius, r, of the vesicles. The geometric center of each molecule was calculated (g_editconf, Gromacs). The geometric center of the respective lipid molecule was then used to place the molecule. Each of the lipid molecules were placed by either plus (r+, here 2.0 nm) or minus (r-, here −2.0 nm) half of the calculated distance between the geometric centers of the lipid molecules in each leaflet of the lipid bilayer (here 4 nm). This procedure generated the outer (r+) and inner (r-) leaflet. The geometric center of each lipid molecule was chosen as the center for translation and rotation for the vesicle construction.

The number of lipids that needs to be placed was calculated by dividing the surface area, which was defined from the two radii (r+ and r-), by the apl of the respective lipid molecules (here 0.66 nm²). The obtained real numbers were converted into integers. For decimals < 0.5, the real number was directly taken as the integer. For decimals ≥ 0.5, the number was rounded up and consequently an additional lipid molecule count was added. When placing the lipid molecules, the distance (square root of the apl) between the lipid molecules was adjusted in the arc length. Due to the different radii for the lipids being positioned in the inner and outer layers, there are different numbers of lipid molecules in each layer. One lipid molecule from the most abundant lipid molecules of the lipid vesicle was placed at the r+ (r- for building the inner layer) position to start with the placement for the vesicle construction. The other lipid molecules were placed in a ‘layer’ (e.g., in x/y plane) around this central lipid molecule at an arc-corrected distance. Then, the next layer was built by placing the lipid at an arc corrected position in the z-direction. This leads to a layer-by-layer build-up of the two leaflets of the vesicle. Next, mixed vesicles were generated by randomly (Python module) replacing a respective number from the total number of lipid molecules calculated above. The number of lipids to be replaced is based on the total number of the most abundant lipid molecules in each leaflet. The protocol described above was merged into the ArcVes software.

2.3.   Simulation of vesicles for CG models

The simulations were performed using Gromacs 2021 with the dry MARTINI v2.1 ff. Simulation were run as a canonical (NVT) ensemble. The integration time step was Δt = 30 fs. Periodic boundary conditions were applied. The cut-off distance of the non-bonded interaction was 1.2 nm. The lipid patches were energy minimized with the steepest decent (5,000 steps). The duration of the production runs was 10 µs. The parameter setup for the production run is the same as described in the dry MARTINI patch.

The initial radii chosen for the vesicles were 10 nm for POPC and DOPC and 5 nm for POPC only (Supplemental Figure 1). The initial radius for the vesicles of the lipid mixtures was 10 nm.

2.4.   Data analysis: apl, radii, membrane thickness (d), number of lipids, estimated diffusivity (D), and lipid density

The apl-value for the membrane patches was generated by taking the number of lipid molecules for each leaflet divided by the area of simulation boxes. For vesicles, the apl values were calculated through the following equation: apl = (4πr²)/(lipid count). Values from all time points were averaged to obtain the averaged apl values for both inner and outer leaflets of the membrane patches and vesicles.

The thickness, d, of the planar lipid bilayer patches was obtained by subtracting the averaged numbers of the z-coordinate from the PO4 sphere of the POPC and/or DOPC molecules for each leaflet in a single time frame. The average was obtained from all time frames. The thickness of the leaflets of the vesicles were calculated based on the respective r+ and r- coordinates for all POPC, DOPC, and CHOL molecules. The center of the beads representing the hydrophilic residues of the head group was taken for the calculations.

The Δz-coordinate was used to identify whether the lipid molecules were in either the inner or the outer leaflet of the patches and for calculating the number of lipids in the planar patches. The CG spheres representing the PO4 and either D2A or D2B were used to calculate the Δz-coordinates of the POPC and DOPC molecules. For the POPC and DOPC lipid molecules at 0 ns in the outer leaflet the Δ-values, between PO4 and either D2A or D2B (PO4 − D2A or D2B), were positive. Flipping to the inner leaflet turned these values negative. For the POPC and DOPC lipid molecules at 0 ns in the inner leaflet, the Δz coordinate is negative and turns positive upon flipping to the outer leaflet. The CG spheres for calculating the Δ-values for CHOL were ROH and C2. The outer and inner locations of each lipid molecule were analyzed at each time point, and the flipping through the simulation was calculated based on the distribution of negative and positive values between successive time steps.

The CG spheres used for monitoring the flipping were PO4 and C3B of POPC and DOPC in the vesicles from ArcVes. CG spheres used for CHOL were ROH and C2. The lipids were monitored by extracting the Cartesian coordinates from the production run. The center point of the vesicle was re-calibrated by calculating the center of mass of all lipids, including those in the inner and outer layers for each frame. The position vector was re-generated for each lipid molecule based on the center point of the vesicle. The internal vector of the lipid was directed from the hydrocarbon tail bead to the phosphate head bead. The inner product of the lipid positional vector and the lipid tail-to-head vector was calculated to distinguish between the inner and outer layers. A positive value for the inner product indicated a lipid molecule within the outer leaflet, while a negative value indicated its position in the inner leaflet.

The translational estimated diffusivity (D) of the lipids in the planar lipid patches was calculated from trajectory-files with a frame resolution of 1.5, 3, 4.5, 6, and 7.5 ns using g_msd (Gromacs). The center-of-mass was calculated for each lipid molecule and was used to derive the position of the molecules in space. The D-values obtained from each resolution were used to calculate the averaged D-values.

In the vesicles, the mean squared displacement (MSD), MSD = < (θ − θ₀)²> × r² [45], of each lipid (using center-of-mass for position in space) was calculated for the Δt = 30 ns time frames over the entire simulation and averaged. Therefore, θ is the angle for which the vector between center and center-of-mass of the lipid has moved during the time step; in reference to the radial coordinates defined for the vesicle, θ₀ is the reference angle with respect to the radial coordinates at the position when the simulation has started and r is the distance of the center-of-mass of the lipids from the center of the vesicle. The great arc was multiplied with the radius and radian span at two time points. The radius was the average of two time points. In the protein data bank (PDB) file, the unit in length is angstroms (Å). The time for the GROMACS (GROningen Machine for Chemical Simulations) output is ‘ps’. The dimension was converted into cm² s⁻¹. The same calculation was performed for time frames of Δt = 60 ns, 90 ns, 120 ns, and 150 ns. D was derived from the slope of the least square fitted line over the four time-steps. The deviation from linearity of the fitted data was less than 1 %. Regarding the number of data points, the five timeframes were equally weighted. The averaged D value was generated by averaging the D values from the five time-frames. The protocol allowed for obtaining enough data for statistic averaging [46],[47].

A grid (100 × 100 bins) was mapped onto the membrane for calculating the frequency-of-occurrence, f, of lipids in the planar system for generating the lipid density map. The center of mass based on the CG spheres was calculated and assigned to the respective bin for each lipid. Accumulation in the bins was performed by using all time frames. The units were f / nm² and referred to as ‘density-of-occurrence’, ρ_f, (count/nm²). Bins were mapped onto the vesicle based on the longitudinal and lateral separation for calculating ρ_f in the vesicles. The non-uniform bin-area was adjusted by correcting each area by multiplying with sin⁻² θ, with θ being the azimuthal angle of the spherical coordinates. Calculations were performed over all time frames and averaged. With an in-house developed software, the color coding in the figures was achieved, with blue indicating low occurrence and red indicating high occurrence.

A spherical area (radius 3 nm) was defined onto the two-dimensional data of the density maps and the ρ_f of each lipid molecule was calculated within each of the respective bins. The values for ρ_f in the bins were normalized to the maximum number of occurrences found in each of the maps. The ρ_f values were averaged over three one-microsecond periods ranging from 7 to 10 µs.

2.5.   Hardware

Data preparation and the analysis were performed on a DELL Precision 5820 Tower workstation with a Nvida GeForce RTX3080 (GPU).

3.   Results

3.1.   Generation of mixed lipid vesicles with cholesterol

Figure 2.  Coarse grained (CG) models of mixed vesicles with a setup-radius of 10 nm. (A) Upper panel: vesicles of POPC/cholesterol (POPC/CHOL) generated by ArcVes (left) and at the end of the simulation (right, 10 µs) with POPC represented in grey (outer leaflet) and light blue spheres (inner leaflet). Lower panel: vesicles of DOPC/cholesterol (DOPC/CHOL) with DOPC represented in green and yellow for outer and inner leaflet, respectively. Cholesterol is represented by blue (outer leaflet) and red (inner leaflet) spheres with the ROH spheres highlighted in magenta for outer and black for inner leaflet to represent the head group of cholesterol. The vesicles present insight into the interior by omitting half of the molecules of the outer leaflet and one quarter of the molecules of the inner leaflets. (B) The same mixed vesicles showing the lipids separately for respective POPC and DOPC (‘PO4’-spheres in blue for outer and red for inner leaflet) and the respective cholesterol molecules (‘ROH’-spheres in blue for outer and red for inner leaflet).

DownLoad: Full-Size Img PowerPoint

Vesicles with various lipid compositions of POPC, DOPC and cholesterol, with radii of 10 nm, were generated using ArcVes (Figure 2A, left side). Vesicles consisted of a binary mixture (e.g., POPC/CHOL and DOPC/CHOL) in a 70: 30 (lipid: CHOL) ratio, and were generated by applying the ratio to both leaflets, thereby leading to a random distribution of all the lipids (Figure 2B, left side).

After an initial reduction within the first 2 to 3 µs of the vesicles with radii of10 nm, the vesicles remained intact during the remaining MD simulation time, adopting mean radii of (8.9 ± 0.1) nm and (9.0 ± 0.1) nm for POPC/CHOL and DOPC/CHOL, respectively (Figure 1A, right side and Table 1). The CHOL molecules in both POPC/CHOL and DOPC/CHOL vesicles seemed to cluster (Figure 2B, right side).

Table 1.  Averaged values for radii [nm] of the vesicles consisting of the mixtures POPC/CHOL, DOPC/CHOL, and POPC/DOPC, the thickness [nm] of the membranes, the number of lipids [count], as well as their diffusivity [10⁻⁷ cm² s⁻¹]. Values are given for the equilibrated initial conditions (initial) before the production run as well as averaged values from the MD simulations (simulation) and presented with standard deviation. The difference between the first and the last 1.5 µs of the simulation is calculated and presented as Δ-values for the radii (Δr) and thickness (Δnm). “/” = different types of lipids or molecules, “i\o” = separation for inner (i) and outer leaflet (o).

	POPC/CHOL	DOPC/CHOL	POPC/DOPC
Radius [nm]
Initial	10.1 ± 0.9	10.1 ± 0.9	9.8 ± 0.9
Simulation	8.9 ± 0.5	9.0 ± 0.5	10.2 ± 0.6
Δr	-0.1 ± 0.0	-0.0 ± 0.0	0.0 ± 0.0
Thickness [nm]
Initial	3.9 ± 0.1 / 3.7 ± 0.1	3.8 ± 0.1 / 3.6 ± 0.1	3.9 ± 0.2/ 3.9 ± 0.2
Simulation	4.2 ± 0.0 / 3.3 ± 0.1	4.2 ± 0.0 / 2.9 ± 0.0	3.9 ± 0.0 / 3.9 ± 0.0
Δnm	0.0 ± 0.0 / 0.1 ± 0.0	0.0 ± 0.0 / 0.0 ± 0.0	0.0 ± 0.0 / 0.0 ± 0.0
Number of lipids [count]
Initial (i\o)	1049/449 \ 1581/678 70 : 30	1049/449 \ 1581/678 70 : 30	749/749 \ 1130/1129 50 : 50
Simulation (i\o)	992 ± 31/551 ± 15 \ 1638 ± 31/576 ± 15 64 ± 1 : 36 ± 1 \ 74 ± 1 : 26 ± 1	1031 ± 12/524 ± 13 \ 1599 ± 12/603 ± 13 66 ± 1 : 34 ± 1 \ 73 ± 1 : 27 ± 1	750 ± 1/750 ± 1 \ 1129 ± 1/1128 ± 1 50 ± 0 : 50 ± 1 \ 50 ± 0 : 50 ± 0
Diffusivity [10−7 cm2 s−1]
i	5.3 ± 0.4 / 5.7 ± 0.4	4.4 ± 0.3 / 5.9 ± 0.3	12.0 ± 0.2 / 12.2 ± 0.3
o	9.8 ± 0.7 / 6.3 ± 0.4	8.9 ± 0.5 / 7.1 ± 0.4	16.7 ± 0.3/ 17.1 ± 0.3

---

 | Show Table

DownLoad: CSV

The time dependence of the radii indicated a very small decrease of the value for both POPC/CHOL and DOPC/CHOL mixture vesicles (Figure 3A, Table 1). The thickness of the bilayer increased during the simulation from the equilibrated initial values of about 3.9 nm (POPC/CHOL) and 3.8 nm (DOPC/CHOL) to 4.2 nm for both POPC and DOPC (Figure 3B, Table 1). The initial membrane thickness of the vesicles was set-up as 4 nm by ArcVes.

Figure 3.  Time dependent calculate properties from the simulations. (A) Radius [nm], (B) thickness [nm], and (C) number of lipids [count × 10²]. The left graphs represent the mixed vesicle POPC/CHOL, and the right graphs present the DOPC/CHOL. Thickness values for POPC and DOPC are shown in green line and cholesterol in red line. Standard deviations (SD [nm]) are shown in grey area representation in the radius plot and for POPC and DOPC in thickness plot. The SD for cholesterol is shown by a black line in the thickness plots. The value for the ArcVes setup-radius is marked as a light grey dashed line. In the plots for the number of lipids, values for POPC, DOPC and cholesterol are all shown in blue and red for outer (o) and inner leaflets (i), respectively.

DownLoad: Full-Size Img PowerPoint

In the POPC/CHOL vesicle, the cholesterol number decreased by approximately 15% in the outer leaflet (initial 678, averaged 576) (Table 1). The POPC number decreased in the inner leaflet by 5% (initial 1049, averaged 992) and increased in the outer leaflet by 4% (initial 1581, averaged 1638). A similar trend is observed for the DOPC/CHOL vesicles, with a move of 11% of the cholesterol molecules from the outer leaflet (initial 678, averaged 603) to the inner leaflet by the average number of lipids and the reverse flipping of the DOPC (2%) from the inner (initial 1049, averaged 1031) to the outer leaflet.

The diffusivity of the POPC and DOPC lipid molecules, is approximately two times larger in the outer leaflet than in the inner leaflet (5.3 × 10⁻⁷ cm² s⁻¹ (i) \ 9.8 × 10⁻⁷ cm² s⁻¹ (o) POPC, 4.4 × 10⁻⁷ cm² s⁻¹ (i) \ 8.9 × 10⁻⁷ cm² s⁻¹ (o) DOPC) (Table 1). The difference for cholesterol between the inner and outer leaflet is less pronounced (5.7 × 10⁻⁷ cm²/s (i)\6.3 × 10⁻⁷ cm² s⁻¹ (o) CHOL in POPC, 5.9 × 10⁻⁷ cm² s⁻¹ (i) \ 7.1 × 10⁻⁷ cm² s⁻¹ (o) CHOL in DOPC). The apl values during the MD simulation are smaller than the values obtained after the equilibration run for both lipids (Suppl. Table 1). Although the apl is about 0.52 nm²/lipid (POPC/CHOL) and 0.49 nm²/lipid (DOPC/CHOL) in the inner leaflets, it is 0.49 nm²/lipid and 0.56 nm²/lipid for POPC/CHOL and DOPC/CHOL in the outer leaflets, respectively.

Density plots of the mixed vesicles identified a mild negative correlation for POPC and DOPC with cholesterol in both leaflets (Figure 4). The negative correlation is more pronounced in the POPC/CHOL vesicle than in the DOPC/CHOL vesicle. The cholesterol molecules show clustering in both leaflets, with a slightly larger density value for the outer leaflets compared to the inner leaflets. There is evidence that the spatial extension of the cluster is larger in the inner leaflets with larger patches for the POPC/CHOL vesicle compared to the DOPC/CHOL vesicle in the outer leaflets.

Figure 4.  2D projection of the density [count/nm²] of the lipids in the mixed vesicles POPC/CHOL (left) and DOPC/CHOL (right) calculated from a 10 µs MD simulation. (A) Density of both leaflets. (B) Density of the outer (upper row, o) ad inner leaflet (lower row, i). The calculation is done by flattening the entire vesicles resulting in axis ranging from 0 to 180 ° and to 360 °. The color coding goes from low to high density as black < dark blue < light blue < green < yellow < red < white.

DownLoad: Full-Size Img PowerPoint

3.2.   Other mixed and pure vesicles as well as membrane patches

The key values of radii, membrane thickness, and number of lipids calculated for mixed vesicle of POPC/DOPC (Suppl. Figure 2A, B) do not show time dependent changes during the 10 µs MD simulation (Suppl. Figure 3A-C). There is almost no change of the radii between the equilibrated initial radius (9.8 ± 0.9) nm and the average radius during the simulation of (10.2 ± 0.1) nm (Table 1). Additionally, the thickness calculated separately for the POPC and DOPC molecules in the POPC/DOPC vesicles do not largely deviate from the initial values of approximately 3.9 nm (appr. 3.9 nm and 3.9 nm, respectively, during the simulation). The number of lipids remains constant for both the lipids and the leaflets (~749/749\1129/1128) (Table 1).

The diffusivity of the lipids in the outer leaflet of the POPC/DOPC mixed vesicle is 16.7/17.1 × 10⁻⁷ cm² s⁻¹, which is higher than the diffusivity of 12.0/12.2 × 10⁻⁷ cm² s⁻¹ in the inner layer (Table 1). The values are higher than those values for the mixed vesicles containing cholesterol, which is not largely different from the values for the pure vesicles. The apl values during the MD simulation are similar to the equilibrated apl values for the lipids (Suppl. Table 1). There is almost no difference between the inner and the outer leaflet.

In respect to the key values, the pure vesicles of POPC and DOPC do show similar behaviors as the POPC/DOPC vesicles. Neither the radii and thickness nor the number of lipids show time dependent changes during the simulations (Suppl. Figure 3 D-F, Suppl. Table 2). The flipping remains therefore very low (~0.1 % for both lipids). The diffusivity remains at 10.4 × 10⁻⁷ cm² s⁻¹ (POPC) and 13.2 × 10⁻⁷ cm² s⁻¹ (DOPC) for the inner leaflets, and shows larger values, 17.5 × 10⁻⁷ cm² s⁻¹ (POPC) and 18.4 × 10⁻⁷ cm² s⁻¹ (DOPC), for the outer leaflet.

Reducing the size of the vesicles consisting of pure POPC molecules does not change the key values during the simulations as compared to the equilibrated initial values (Suppl. Table 2). The flipping activity in both leaflets is more frequent than the flipping in the larger POPC vesicle. The diffusivity of the lipids in the outer layer, which is approximately 15.8 × 10⁻⁷ cm² s⁻¹, remains in the same range as those in the larger vesicles, though it is very much reduced in the inner leaflet to about 8.9 × 10⁻⁷ cm² s⁻¹; this trend is the same for the POPC/DOPC vesicle, which show lipid clustering (Suppl. Figure 4). The apl values during the MD simulation are almost unchanged compared to the equilibrated value for the lipids (Suppl. Table 1). Some exceptions are values that are calculated from the outer leaflet of POPC vesicles with a 10 nm diameter (0.74 nm²/lipid) and in the inner leaflet of the POPC-5 vesicle (0.79 nm²/lipid).

Comparing the calculated membrane thicknesses of the vesicles with lipid patches of the same type of lipid and mixtures (Suppl. Figure 2A–D) shows that the thickness values for the dry patches are in the same range (e.g., ~3.7 nm for pure and mixed POPC/DOPC patch (Suppl. Table 3)), though the values are lower (3.3 nm (POPC) and 3.2 nm (DOPC)) in the presence of CHOL (Suppl. Table 4). When conducting the simulation with the fully hydrated system, the thickness is slightly larger. No change in the number of lipids, which is due to flipping of the molecules, is observed independent of the simulated patch system. Moreover, clustering of the lipids is not observed (Suppl. Figure 4).

In the presence of cholesterol, similar to the behavior in the vesicles, diffusivity of the POPC and DOPC molecules is reduced (e.g. 2.7\2.7 × 10⁻⁷ cm² s⁻¹, dry, POPC in POPC/CHOL patch and 2.7\2.3 × 10⁻⁷ cm² s⁻¹, dry DOPC in DOPC/CHOL patch) compared to the values of the respective pure patches (e.g. 4.1\4.0 × 10⁻⁷ cm² s⁻¹, dry, POPC patch and 5.1\5.7 × 10⁻⁷ cm² s⁻¹, dry DOPC patch) (Suppl. Tables 3 and 4).

In summary, cholesterol induces a thickening of the membrane in the vesicles and membrane patches. In consequence, the presence of cholesterol in lipid vesicles and patches leads to (i) the decreased diffusivity of lipids compared to the diffusivity calculated for systems without cholesterol, and (ii) a larger diffusivity in the vesicles of the lipids in the outer leaflets compared to the inner leaflet. The key values such as the radii and thickness and the number of lipids change during the simulation, thereby indicating the adaptation of the POPC/CHOL and DOPC/CHOL systems to the thermodynamically stable systems. Similar adaptations are not observed for lipid mixtures such as POPC/DOPC, the pure vesicles and patches.

4.   Discussion

4.1.   Assessment of the system

There is a variety of software and methodologies available to generate curved membranes up to the point when the curvature turns into a vesicle (reviewed in [48]). The methodologies cover the formation via self-assembly [27]–[29],[49], using triangular surfaces as building units [50], and Charmm-Gui Martini Maker [30]–[32]. With each of these methods, it is difficult to control the exact numbers of lipids in each leaflet during the generation of the vesicles prior to the production run. With the ArcVes software, it is possible to generate the vesicles on demand by distributing lipids individually and stepwise across the leaflet. Consequently, the stability and dynamics of the lipids can be analyzed after multi microsecond MD simulations. Various types of mixtures can be generated, and proteins can be inserted. Any type of ff can be applied after the generation of the vesicles, and whether they are built with either CG or united-atom models of the molecules. With this option, ArcVes fills a niche in the ensemble of methodologies, thus allowing for a thorough analysis of mechanical properties of vesicles in silico.

At this stage, dry Martini ff in combination with ArcVes is chosen for the MD simulations, since the ff allows for fast screening of vesicle compositions, thereby fulfilling requirements of simplicity and speed for applications in bioengineering. The software allows for the application of newer ffs, including Martini 3.0, since ArcVes is a software that places molecules according to a macroscopic physical parameter, such as the apl.

The lipid dynamics of the vesicles are based on the CG lipid models operated under the MARTINI ff. The ff is optimized to reproduce partition free energies of the molecules [51]. It is especially used to simulate planar lipid membranes as well as lipid-protein systems. It is anticipated that the molecular models of the lipid molecules are able to provide estimated data presenting dynamical properties such as diffusivity, as well as when it is used in spatially different topologies [34],[52]. The spheres in the different lipid tails of POPC and DOPC are labeled with different interaction strengths: three spheres of the palmitoyl tail with 4.5 kJ/mol in POPC versus the same number of spheres of the oleoyl tail with 4.0 kJ/mol in DOPC under dry Martini ff, as reported in the literature [53]. The different interaction parameters explain the higher diffusivity of the mixed vesicle POPC/DOPC compared to the cholesterol rich mixed vesicles, as well as the pure DOPC vesicle compared to the POPC vesicle.

The positioning of each of the lipid molecules is based on their centers-of-mass onto the surface of a sphere at distances calculated from the apl values, which results in the stabilization of the vesicles during multi microsecond MD simulations. A first approximation is the use of the surface-apl to define the distance between the lipids, assuming that the lipids will adjust their positions during the MD simulations accordingly. The apl values are taken from experiments, in which a combination of small-angle neutron and X-ray scattering techniques on unilamellar vesicles (ULVs) with approximate diameters of 60 nm are used at various temperatures [54]. The apl value of 0.66 nm² is adopted from the data of POPC at the chosen temperature for the simulation. Moreover, this value leans on values identified by experiments for other lipid molecules (e.g. 0.725 nm² for DOPC and 0.653 nm² for DOPS [55], both in the fluid phase). Values of 0.633 nm² and 0.660 nm² for DPPC (measured at 50 °C) and DSPC (measured at 65 °C), respectively, are reported from NMR spectroscopic measurements on large multilamellar vesicles, with approximated diameters in the range of 6 to 10 mm [56].

Additionally, similar apl values are calculated in MD simulations on planar lipid bilayers and vesicles, applying the same lipid models and the same ff [57] or other ffs solely used in planar systems [58]–[60]. The apl values decrease in the presence of physiological ions, which has been shown by a set of simulations of planar lipid bilayers at similar temperature reported in this study [61]. In the absence of a substantial number of additional ions other than those to compensate for the negative charges of the DOPS molecules, the apl value of 0.66 nm² applied in this study is considered as a compromise.

The size of the vesicles generated in this study (about 20 nm diameter, 10 nm radii) is still smaller than experimentally derived diameters for large unilamellar vesicles (LUVs) of POPC (about 30–100 nm diameter) [62] and cellular vesicles (about 100 nm diameter [63]) and about 2,000 times smaller compared to artificially derived vesicles of e.g. DOPC vesicles (about 20 µm diameter [64]). The vesicles remain stable during the simulation, and no rupture is observed since the vesicles are not inflated by water beads, which would lead to a thinning of the membrane otherwise [65]. A large self-diffusion coefficient is reported for much larger vesicles with a diameter of 200 nm and are reported to be slower in reversed vesicles with a diameter of 20–30 nm [66].

Experimentally, however, it has been shown that the diffusivity of lipids is higher when they are within vesicle compared to the situation where they are a part of supported planar lipid membranes [67],[68]. In the context of this study, the planar system can also be seen as a ‘free hanging’ system thus representing an enormously large unilamellar vesicle. Thus, it is anticipated that the values of the diffusivity elevate with an increasd radius of the simulated vesicles.

The diffusivity of lipid molecules calculated here for planar patches are in the range of reported values from CGMD simulations of very large POPC patches [52]. However, the values seem to be dependent on the modeling conditions, especially on the thickness of the water layer in which they are embedded into [52]. In smaller patches, they are also dependent on the cut-off conditions (e.g., POPC 4.0–4.5 × 10⁻⁷ cm² s⁻¹ [69] Amber), lipid types (e.g., united atom MD simulations on POPC, 0.62–0.64 × 10⁻⁷ cm² s⁻¹, DMPC 0.64–0.55 × 10⁻⁷ cm² s⁻¹ [70], (accMD and Charmm) and DOPC (0.84 ± 0.04) × 10⁻⁷ cm² s⁻¹ [71] (Charmm)) and temperature (e.g., experimental study [72]). The presence of proteins reduces the diffusivity of the lipid molecules from ‘bulk-like’ values in a lipid-type dependent manner. For example, the diffusivity for DPPC (1,2-dipalmitoyl-sn-glycero-3-phosphocholine (16:0 PC)) decreases from approximately 0.3 down to slightly less than 0.3 × 10⁻⁷ cm² s⁻¹. Additionally, the diffusivity decreases from approximately 1.4 down to approximately 1.1 × 10⁻⁷ cm² s⁻¹ for DLiPC (1,2-dilinoleoylphosphatidylcholine (di-18:2 (cis) PC) [19] using CGMD, or POPC and POPS from the range of 10–12 × 10⁻⁷ cm² s⁻¹ for a POPC/POPS/POPE/CHOL mixture down to 7–9 × 10⁻⁷ cm² s⁻¹ at a distance below 2 nm from the membrane protein [73]).

A low diffusivity of lipid molecules such as POPC (0.9–1.0 × 10⁻⁷ cm² s⁻¹) is experimentally calculated for giant unilamellar vesicles (GUVs) [74] and virus-like vesicle systems [75]. Other reports show that the POPC diffusivity decreases from 3.10 × 10⁻⁸ cm² s⁻¹ (1 mM) to around 2.3 × 10⁻⁸ cm² s⁻¹ (25 mM) depending on lipid concentration [62], as well as for planar systems when immobilized on a solid (mica) support (for DOPC the diffusion coefficient changes from (0.78 ± 0.08) × 10⁻⁷ cm² s⁻¹ to (0.31 ± 0.03) × 10⁻⁷ cm² s⁻¹) [67].

Overall, the vesicles generated in this study mimic in vitro conditions of unilamellar vesicles or liposomes.

4.2.   Cholesterol in vesicular and planar systems

Cholesterol affects many characteristics of lipid bilayers. It induces the thickening of the lipid membrane and induces intrinsic lipid curvatures, which consequently affects other characteristics such as fluidity and phase separations [76]. In relation with the calculated key values, phase separation leads to the formation of domains enriched in cholesterol [77],[78]. This feature is also observed in the vesicles studied here. Similarly, the flip flop of cholesterol is larger than other phospholipids. This feature can be used to create asymmetric lipid vesicles [22] and induces larger biological events such as budding and fusion [79]. Flipping is also proposed by MD simulations in planar lipid patches [80] and the clustering behavior of lipid molecules [81], which is not observed in the patches investigated in this study.

The concentration of cholesterol in either leaflet of the plasma membrane is still controversially discussed. Depending on the stimulations of the cell, some investigations identify more cholesterol in the inner leaflet [82], and other investigations show the opposite results, with more cholesterol in the outer leaflet [83],[84]. MD simulations on larger mixed vesicles (150 nm in diameter), such as mimicking the HIV-1 envelope, report more cholesterol in the outer leaflet [15].

The lateral diffusivity of cholesterol within lipid bilayer patches is proposed to be dependent on the location of the cholesterol within the lipid bilayer. Cholesterol positioned within the leaflets is slower than when positioned in the center of the bilayer [85].

Calculations of the mean square displacement (MSD) of lipids in large unilamelar complex mixed vesicles (150 nm in diameter), such as mimicking the HIV-1 envelope using CG MD simulations, reveals that the diffusivity values of the lipids are in the range of 10⁻⁶ cm² s⁻¹. The diffusivity shows larger values for the inner leaflet than for the outer leaflet [15]. This is attributed to the lipid rafts identified in the outer leaflet of the vesicles and the higher cholesterol content in the outer leaflet. These vesicles were generated by using a packmol library, which uses a spatial based packing optimization [86],[87], in combination with an alchemical phasing procedure [88]. Wet Martini is used for all atom models of the lipids.

In this study, the mixed vesicles consist of binary mixtures of lipids. In these mixed vesicles, the cholesterol tends to flip to the inner leaflet more than to the outer leaflet, and the lipids POPC or DOPC flip to the outer leaflet to a lesser extent.

The lower number of lipids in the inner leaflets supports the formation of cholesterol patches. The slight increase in membrane thickness shows domain and raft-like formations [77],[78]. Due to its chemical moieties of a polar headgroup, rigid steroid ring structure and nonpolar hydrocarbon tail, cholesterol is suggested to be attracted by saturated lipids and repelled by unsaturated lipids leading to its clustering and in consequence to raft formation [78]. Furthermore, lipid molecules with low-melting characteristics due to their unsaturated chains (here POPC one chain, DOPC two chains) ‘push’ cholesterol away, in constrast to high-melting lipids exerting the opposite effect of ‘pulling’ the cholesterol towards them (‘push-pull mechanism) [89]. Experimental studies validate the combination of lipid and cholesterol characteristics in patches and rafts and for the formation of those structures, as well as simulation methods ([90] and references therein).

The data obtained from the simulations in this study support the explanation that the lower diffusivity of the lipid molecules in the inner leaflet and the repelling (‘push’) character of the unsaturated and low-melting lipids, POPC and DOPC, support the patch formation of cholesterol molecules. The cholesterol patches consequently enhance the flipping of more cholesterol molecules to the inner leaflet.

Overall, the dynamics of the lipid molecules decreases, thereby resulting in decreased diffusivity values of the lipid molecules in the presence of cholesterol, and within the inner leaflet.

4.3.   Using the area-per-lipid value

The apl value, which is an intrinsic property of the lipid molecules, was used to identify the number of each lipid in either the inner or the outer layer. Thereby, a single apl value is used for all lipid types (i.e., POPC, DOPC and CHOL). This value is used to calculate the numbers of lipids from the available surface area of each leaflet by taking care of the arc correction when applying the distance between the lipids to be placed. This approach is seen to be an optimized mimicry of in vivo conditions. During the simulations of the mixed vesicles with CHOL, the apl adopts lower values compared to the initial value, though with no tendency of larger or lower apl values for the inner or outer leaflet. With the suggestion that a higher tension increases the diffusivity and apl [91], it is reasonable that the higher diffusivity for the outer leaflet observed in this study represents the higher tension in the outer leaflet of the vesicles. Additionally, the calculation of the diffusivity uses the radius dependent angular displacement. Despite using almost similar molecular conditions compared to the large liposome simulations [15], the diffusivity in this study is larger in the outer layer compared to the inner layer.

5.   Conclusions

Using the presented software ArcVes, in-vivo-like vesicles can be generated and used in MD simulations to explore the mechanics and characteristics of lipids and their mixtures in vesicles. The simulations generate key parameters that are similar to those reported in in-vitro/vivo vesicular systems, and also indicated by MD simulations in the literature.

The surface area of the lipid molecules is not altered when simulating lipid molecules as CG models on strongly curved surfaces as those used in this study compared to simulations of the same molecules as CG models in flat lipid membranes. Consequently, the value of apl can be used as a parameter to map lipid molecules onto a sphere and generate vesicles that remain stable during MD simulation. The protocol presented in this study can be used for screening a series of other lipid compositions as well as when proteins or peptides are inserted.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.