2.1. MD Simulation Details

All the simulations are completed using GROMACS 2018.1 and 2024.1 versions. The necessary GROMACS file (.gro) and topology (.top) files are generated with the Insane (INSert membrANE) [55] tool, which uses the forcefield parameter from the Martini force field. A simulation box with dimensions 15 nm x 10 nm x 9 nm, containing lipids and sterol molecules was generated. The concentration of Na⁺ and Cl⁻ ions in the simulation was maintained at 0.15 M to mimic physiological ionic conditions and ensure charge neutrality, thereby providing a stable environment for the lipid bilayer under investigation. The output configuration is written to .gro file in a cubic periodic boundary condition, solvated with water. The same conditions were maintained throughout the experiment. The portable topology files (.itp) (Martini V2.0) which contain the information on all the necessary parameters required for the simulation were taken from the cgmartini website. The topology (.top) file was written using the result generated by the Insane tool and required .itp files were included. In each leaflet of the bilayer, an exact number of lipids is maintained (Table 2). Before the start of the simulation, energy minimization and equilibration steps were performed.

Energy minimization was performed using the steepest descent algorithm, with a maximum of 1000 steps, using the Verlet cutoff scheme for non-bonded interactions and applying periodic boundary conditions in all dimensions. Coulombic interactions were controlled with the reaction-field method, truncating interactions at 1.1 nm and setting the relative dielectric constant to 15. Van der Waals interactions used a cutoff and potential-shift, with a cutoff distance of 1.1 nm. For equilibration, the Martini force field was used to stabilize the system at the desired temperature and pressure. Time steps of 1 fs, 5 fs, 10 fs, and 30 fs were employed. A Berendsen thermostat maintained the temperature at 298 K with a 1 ps time constant, and pressure was set to 1 bar with a 6 ps time constant. Initial velocities were generated at 298 K, and constraints were applied using the LINCS algorithm.

Confirming effective minimization and equilibration run, a newly configured .mdp file was used to initiate the production run. The integrator algorithm chosen is the standard MD integrator with a time step of 0.03 ps and a total of 400 million steps, corresponding to an overall simulation time of 12 µs. Centre of mass motion removal is configured to occur every 100 steps. For energy minimization, a convergence tolerance of 10 and an initial step size of 0.01 are specified, with a maximum of 20 iterations allowed. The steepest descent method is used with corrections applied every 1000 steps. Neighbor searching employs the Verlet cutoff scheme, updating the neighbor list every 20 steps, with a cutoff radius of 1.0 nm. Electrostatics are managed using the reaction-field method with a Coulomb cutoff of 1.1 nm and a relative dielectric constant of 15. Van der Waals interactions are treated with a cutoff of 1.1 nm, without long-range dispersion corrections. Temperature coupling is managed via the v-rescale method with a reference temperature of 300 K and a time constant of 1 ps. Pressure coupling is semi-isotropic using the Parrinello-Rahman method with a reference pressure of 1 bar and a compressibility of 3 × 10⁻⁴/bar, with a time constant of 12 ps. Velocities are generated at the start of the simulation corresponding to a temperature of 300 K. The LINCS algorithm is used for constraint solving with an order of 8 and 2 iterations. No user-defined groups, walls, or external electric fields were applied, ensuring a straightforward and unbiased molecular dynamics simulation setup with the use of MPI (Message Passing Interface), the simulation is performed on multiple processors simultaneously, enhancing performance and reducing simulation time.

To ensure consistent conditions in every simulation, the same mdp files are used for energy minimization, equilibration, and production runs.

2.2. Simulation Analysis

An energy file (.edr), a GROMACS compressed trajectory file (.xtc), and a lossless trajectory file (.trr) were generated as a post-simulation output. The .xtc file contains compressed atomic coordinates at various time steps, making it useful for trajectory visualization and analysis. In contrast, the .trr file stores full-precision, lossless trajectory data, including atomic coordinates, velocities, and forces, and is primarily used for more detailed analyses. This level of data is required for performing precise calculations and detailed analysis. The .edr file records the simulation’s configurations, including integrator settings and time steps, providing a reference for the simulation parameters. Together, these files offer complete post-simulation analysis, enabling tasks like visualizing molecular dynamics to perform accurate computational calculations.

Area per lipid (APL) provides valuable insights into the structural properties, and stability of lipid membranes, also understanding sterols interaction, the packing density of the bilayer, and its behavior under different conditions, which are important for understanding cellular developments and designing membrane-targeted drugs. A direct method of APL calculation can be done using the inbuilt GROMACS tool gmx_energy. By determining the area of the simulation box and dividing it by the total number of lipids in the single leaflets assuming both have an equal number of lipids. FATSLiM v0.2.1 tool uses a different approach to calculate APL by completely identifying lipids and membranes initially. Each lipid has an available area, and the area per lipid may be roughly calculated using the average of the accessible areas. As a result, it calculates the membrane area which is the average of the two leaflet areas, and the leaflet area is just the total of all individual areas [12]. GridMAT-MD 2.0 (Grid-based Membrane Analysis Tool for use with Molecular Dynamics) tool, reads from a coordinate file (.gro) which only has a single frame [13]. A parameter file is fed into the system in which the user has to define the input parameter file, reference point, solvent, ions, number of grid points across each axis, and the type of grid to be generated as an output i.e., column, matrix, or vector. The output is analyzed to get the result of APL.

The Bilayer Thickness (BT) study of a lipid membrane is essential because it influences membrane permeability and the function of embedded sterols. Understanding bilayer thickness helps elucidate the physical properties of the membrane and its role in various biological processes, such as signaling and transport. For example, a thicker bilayer may stabilize certain receptor conformations, enhancing or inhibiting their ability to transmit signals, also may make it harder for small molecules to diffuse. GROMACS tool gmx density is used for the calculation of BT. Provide the proper index file by choosing the phosphate beads as the headgroup. The generated output can be analyzed using XMGrace software (Version: Grace-5.1.25), where the bilayer thickness (BT) is determined by measuring the distance between the headgroup peak positions. Rather than change in the input index file and reference atom, both FATSLiM and GridMAT-MD use similar approaches for the thickness calculation [12,13]. Later results from direct calculation, FATSLiM, and GridMAT-MD are compared.

Studying Partial density is important because it provides insights into the spatial distribution of different components within a system. This information is crucial for understanding the structural and dynamic properties of complex mixtures, such as lipid bilayers or sterol solutions. By analyzing partial densities, one can determine how molecules are arranged and how they interact with each other. Here also the GROMACS gmx density tool is used in the calculation.

Investigating the Order parameter of lipid and sterol simulations is crucial for understanding the degree of lipid tail alignment and membrane fluidity. Defined as

To generate the output, a Python script do-order-gmx5.py is used within the GROMACS environment. The script is called by giving. xtc and .tpr files as input along with the start and end time for analysis (in ps), specific frequency (in steps) at which the order parameters are calculated, type of lipid, and number of lipids in the leaflet. The calculation is performed along the z-axis and averaged over a specified number of simulation steps. The resulting output, including the profile data, can be used to compare with published results.

Radial Distribution Functions (RDF) is an important method for understanding how molecules are arranged around a reference point in a system. It provides a detailed analysis of molecular packing, local structure, and how different molecules interact with each other. For example, in lipid bilayers or sterol mixtures, RDF analysis helps in understanding how lipids interact with each other or with cholesterol. In this study, we used the inbuilt GROMACS tool gmx rdf to analyze these interactions and better understand the structural and organizational characteristics of the system. Later the generated results were visualized using XMGrace software(Version: Grace-5.1.25).

A representative simulation result of the lipid bilayer with various types of lipids i.e., DLPC, POPC, DGPC, and DIPC, with 50% cholesterol, was shown in Figure 1a–d, and bead representation of PC lipid and cholesterol was shown in Figure 1e. This will provide a visual overview of the cholesterol distribution in the lipid bilayer.

3.1. Structural Properties of Phosphatidylcholine Lipids and Cholesterol

The structural configuration of phosphatidylcholine lipids: 1,2-dilauroyl-rac-glycero-3-phosphocholine (DLPC), 1-palmitoyl-2-oleoylphosphatidylcholine (POPC), 1,2-di-(11Z-eicosenoyl)-sn-glycero-3-phosphocholine (DGPC), and 1,2-dilinoleoyl-sn-glycero-3-phosphocholine (DIPC) were analyzed. Since these molecules (DLPC, POPC, DGPC, and DIPC) are all phospholipids, they share a central structure of glycerol, a phosphate group linked to choline, and two fatty acid chains attached to the glycerol. DLPC has two short, saturated lauroyl chains (12 carbons each) (Figure 2a), POPC has one saturated palmitoyl chain (16 carbons) and one unsaturated oleoyl chain (18 carbons with a double bond) (Figure 2b). DGPC has two longer, unsaturated eicosenoyl chains (20 carbons with a double bond at the 11th position) (Figure 2c), and DIPC features two linoleoyl chains (18 carbons with two double bonds each) (Figure 2d). In the .itp files, the atoms are represented in terms of beads. Below is the detailed comparative explanation of molecules as written in the respected .itp files. The positively charged nitrogen atom (NC3) represented by the atom type Q0, bonded to a phosphate group (PO4) with the atom type Qa. The two neutral glycerol atoms (GL1 and GL2) are represented by the atom type Na. In DLPC each fatty acid chain is represented by three beads: C1A, C2A, and C3A for one chain, and C1B, C2B, and C3B for the other, all of type C1, indicating that they are simple, neutral, and do not have any branching or double bonds. POPC has more complex chains, the first chain starts with a carbon atom of type C1 (C1A), followed by a carbon of type C3 (D2A) (which indicates a double bond), and then continues with two more carbon atoms of type C1 (C3A and C4A). The second hydrocarbon chain is simpler, consisting of four consecutive carbon atoms of type C1 (C1B, C2B, C3B, C4B). In DGPC, the backbone connects to two hydrocarbon chains, each chain starts with a carbon atom (C1A and C1B) of type C1, continues with a sequence of carbon atoms (C2A, D3A, C4A, C5A for the first chain, and C2B, D3B, C4B, C5B for the second chain), where the D3A and D3B atoms are of type C3. Among the four lipids, DIPC has the most complex structure. Its first hydrocarbon chain starts with a carbon atom of type C1 (C1A), followed by two carbons of type C4 (D2A and D3A), indicating the presence of double bonds, and ends with a carbon of type C1 (C4A). The second chain has a similar structure (C1B, D2B, D3B, C4B). This structural configuration, with two extended hydrophobic tails and a hydrophilic headgroup, is characteristic of the typical arrangement of a phosphatidylcholine lipid, which is important for forming bilayer structures in membrane simulations.

In cholesterol molecules, the ROH group is represented by the SP1 bead. The subsequent groups are represented by SC1 beads (R1, R3, R4, R5) and SC3 (R2) (Figure 2e). The aliphatic tail segments in cholesterol are represented by SC1 (C1) and C1 (C2) beads. Bond interactions are defined within the molecule to maintain its structure. For example, the bond between C1 and C2 has a length of 0.425 nm and a force constant of 1250. Structural integrity is further confirmed by constraints between key beads (SP1, SC1/SC3, SC1/C1). The correct spatial arrangement and stability of non-linear components in the molecule are defined with virtual sites. Dihedral angles and exclusions focus on maintaining proper orientation and stability. These configurations ensure accurate modeling of the molecular dynamics and interactions of cholesterol within lipid membranes.

3.2. Simulation Result Analysis of Different Concentrations of Cholesterol with Lipids

The. xtc,. tpr,. edr, and .gro files were used to analyze the results such as bilayer thickness, area per lipid, and order parameter. Sterols are a rigid, planar ring structure that inserts between lipid molecules, promoting the extension and alignment of lipid chains, which results in lipid tails being more stretched. The data demonstrate that sterols significantly influence the structure of lipid membranes, though the extent of their effects varies across different lipid types. Sterol addition generally results in thicker bilayers, reduced lipid area, and better alignment of lipid tails, enhancing membrane stability (Table S1). However, the degree of these changes depends on the specific lipid composition, with some membranes responding more dramatically to sterol incorporation than others. The influence of sterol concentration on lipid membranes was explored across four different types: DLPC, POPC, DGPC, and DIPC.

The effect of sterol concentration on bilayer thickness showed distinct trends across the different lipids. For DLPC, the thickness increased from 3.445 nm at 0% sterol to 4.013 nm at 50% (Figure 3a), reflecting the membrane’s expansion due to sterol incorporation. In POPC, the change in bilayer thickness with increasing cholesterol concentration is relatively minimal. The lowest thickness observed is 3.897 nm at 30% cholesterol, while the highest is 4.281 nm at 10%. Unlike DLPC, POPC does not exhibit a consistent increasing trend. This may be because of the presence of an unsaturated oleoyl chain in POPC (Figure 3b). DGPC and DIPC showed a more pronounced change, with bilayer thickness increasing from 4.370 nm to 4.606 nm (Figure 3c) and 3.461 nm to 3.732 nm (Figure 3d), peaking at 40% sterol concentration (4.667 nm and 3.795 correspondingly), suggesting that the effect of sterols on bilayer expansion may level off around 40-50%. The cholesterol aligns itself between the acyl chain tails more tightly due to its rigid structure, which results in a stiffer bilayer and reduces kinks, leading to a thicker bilayer [50,57]. The deformity observed at 50% cholesterol in lipids other than DLPC is attributed to phase behavior and the presence of a double bond in the lipid. The lipid tails can only be extended and aligned to a certain degree, beyond which additional sterols may not contribute to further structural changes.

Sterol addition led to a decrease in the area per lipid across all lipid types, indicating tighter molecular packing within the membranes. For DLPC, the area per lipid dropped from 0.591 nm² to 0.399 nm² as sterol concentration increased (Figure 4a). POPC also saw a substantial reduction, from 0.642 nm² to 0.425 nm² (Figure 4b). DGPC’s area per lipid decreased from 0.666 nm² to 0.428 nm² (Figure 4c), while DIPC showed a greater decline from 0.741 nm² to 0.474 nm² (Figure 4d), reflecting less pronounced membrane condensation. This decrease in area per lipid significantly alters the membrane surface properties. The increase in cholesterol helps increase the stability of the membrane, for which the membrane shows resistance to changes in temperature and mechanical stress, which in turn contributes to the membrane’s structural integrity and functional efficiency.

The order parameter, which reflects the alignment of lipid acyl chains, increased with rising sterol concentration in all cases; however, the extent of this increase varied among different lipid types. In DLPC, the order parameter rose significantly from 0.479 to 0.746, indicating enhanced tail alignment (Figure 5a). POPC exhibited a smaller increase, from 0.375 to 0.585 (Figure 5b). For DGPC, the order parameter climbed from 0.331 to 0.445 (Figure 5c), showing improved alignment, while DIPC experienced a more modest rise from 0.227 to 0.262 (Figure 5d), but the highest being at 40% concentration i.e., 0.272, suggesting a limited impact of sterols on tail alignment compared to the other lipids.

The results are further confirmed by comparing the values obtained from different computational tools such as FATSLiM v0.2.1 and GridMAT-MD 2.0. (Table S1).

3.3. Partial Density and Radial Distribution Function (RDF) Analysis of Phosphatidylcholine Lipids Bilayer

The simulation analysis showed a significant observation at the 50% cholesterol and Lipid composition. The increase in cholesterol leads to an increase in order parameters which results in the reaching the saturation point. However, in DLPC the bilayer thickness follows the increasing trend, and the order parameter is higher compared to other lipids. This proves that DLPC may still have some space in between bilayer lipid tails, till it reaches its saturation point. Figure 6a depicts the partial density of the headgroup of lipid PO4, indicating their well-defined position in the bilayer, close to the water-lipid interface. From the previously observed phenomenon of DGPC and sterol interaction i.e., a decrease in the bilayer thickness at the 50% sterol concentration, the same pattern is observed in the case of DIPC (3.732 nm) (Figure 7a), which is very low compared to other lipids.

The ROH groups show broader and less defined peaks, indicating a more dispersed distribution within the bilayer (Figure 6b). The graph depicts the position of cholesterol presence in the bilayer wherein 0 represents the center of the bilayer and −2 and 2 are the two leaflets. The more stretched the ROH is towards the bilayer, the more it shows the order parameter (Figure 7b). In Figure 6b, the distribution of cholesterol within the bilayer shows an increasing order which can also be visualized in Figure 1 and can be compared with the graphs, resulting in the decrease of order parameter from DLPC to DIPC (0.746 to 0.246) (Figure 5). The cholesterol moving in-between the bilayer rather than aligning itself with the lipid headgroup is due to the presence of double bonds in the lipid tail (Figure 1). ‘DLPC allows more cholesterol to align toward the headgroup than DIPC. The double bond in the hydrophobic tail creates a kink, which increases the area per lipid but disrupts the packing of neighboring lipids. This disruption affects the bilayer’s structure, potentially hindering cholesterol’s ability to pass through. The impact of the double bond on bilayer movement can be visualized in the density map of ROH along with the double bond and terminal carbon atom of the lipid tail (Figure 7). POPC has only one double bond in a single chain (D2A) (Figure 7b), it has a lesser distribution of cholesterol towards the headgroup compared to DLPC (Figure 7b). DIPC has the highest number of double bonds (D2A, D3A, D2B, and D3B) (Figure 7d) making it have a more flexible tail, which results in the highest distribution of cholesterol between the bilayer than the other four simulated lipids. The flexibility of the acyl chain strongly influences cholesterol alignment and distribution within the lipid bilayer. Lipids like DLPC, with fully saturated chains, tend to exhibit more stretched acyl chains compared to lipids such as DIPC, which contain double bonds. An increasing number of double bonds in lipid tails reduces chain stretch and increases flexibility, particularly in the presence of cholesterol. This increased flexibility disrupts the alignment of cholesterol between the lipid’s hydrophilic headgroups, leading to lower bilayer order parameters and faster saturation of order parameter values at lower cholesterol concentrations. Figure 8 visualizes this phenomenon confirming the movement and distribution of molecules inside the bilayer depending on its structural properties.

The Radial Distribution Function (RDF) graphs (Figure 9) provide a clear picture of how cholesterol interacts with different lipids—DLPC, POPC, DGPC, and DIPC—revealing the role of lipid structure in shaping the bilayer’s molecular organization. DLPC shows the strongest interaction with cholesterol (Figure 9a), as seen in the sharp, high RDF peak around 0.5 nm and consistent values beyond 1 nm. This is due to its fully saturated lipid tails, which promote tighter packing and better cholesterol alignment. On the other hand, DIPC exhibits the weakest interaction (Figure 9d), with a smaller peak and a rapid drop-off, reflecting the flexibility of its tails with multiple double bonds. These double bonds create kinks that disrupt packing and make it harder for cholesterol to align, resulting in reduced bilayer order. POPC and DGPC fall in between (Figure 9b,c), with POPC performing slightly better due to having only one double bond, causing less disruption. These findings align with observations of bilayer thickness and order, where DLPC forms a more structured bilayer, while DIPC’s flexible tails lead to a looser, more disorganized arrangement. Compared to cholesterol-lipid RDFs, lipid-lipid interactions are generally stronger and more consistent, especially for DLPC and POPC. DIPC remains the most disordered lipid in both cases, highlighting how double bonds in the lipid tails disrupt packing and reduce order. Saturated lipids like DLPC promote tight packing, while unsaturated lipids like DIPC introduce more flexibility and disorder (Figure S1). Overall, this analysis highlights how lipid tail saturation and flexibility directly affect cholesterol distribution, lipid packing, and bilayer stability.