Order parameters

where θ is the angle between the lipid tail bonds and the membrane normal.

The final order parameter for each technique is averaged over all leaflets in the system, and Membrainy will also produce separate order parameters for each lipid type and leaflet. For atomistic and united-atom force fields, Membrainy plots the values of −S _CD for each carbon along the lipid tails. This experiences maximum order at 0.5 and disorder at -1, whereas the Martini force field experiences maximum order at P ₂=1 and disorder at P ₂=−0.5. Membrainy can also produce histograms of the angles measured by each technique. To maximise performance, the order parameter algorithms are multithreaded, where each lipid tail type (e.g. POPE-palmitoyl, POPE-oleoyl, etc.) is assigned its own thread, allowing much of the analysis to be conducted in parallel.

Headgroup orientations

Membrainy calculates lateral and axial headgroup orientations, producing a histogram for each lipid type. The lateral angles are calculated by establishing a headgroup vector from two reference atoms, one being the phosphorous atom and the other being another atom on the headgroup. This vector is then projected onto the membrane normal to produce an angle. The histograms are plotted in the range -90 to 90 degrees, where a value of 0 indicates the headgroup is parallel to the membrane surface and positive angles indicate the headgroup is pointing away from the membrane. Axial angles are calculated by projecting the headgroup vector onto the membrane surface, taken as the xy plane, to produce a radial angle between 0 and 2 π. Each axial angle is plotted for each lipid over time. This algorithm has been multithreaded, where each lipid type is assigned its own thread and run in parallel.

2D surface maps

where A _i,j is the highest atom in cell i,j, \(\phi _{i,j}^{n+1}\) is the resulting scalar value produced by the method, and the final term is the sum of the neighbouring cells’ scalar values. Iterating over this method produces a scalar field of successive displacement, generating a series of Gaussians that can be scaled and mapped to a colour to produce a contour map of the leaflet surface. These maps also behave as density maps, producing more prominent Gaussians in regions of the lattice containing a high density of atoms, such as lipid tails in the gel phase. The scalar field is colour-coded such that blue regions indicate thin or sparsely populated regions of the leaflet, red indicates thick or densely populated regions, with green between the two. Black areas represent a hole or pore in the leaflet, which is identified by unpopulated regions of the lattice. A map for each leaflet is displayed through a graphical interface in real-time and can be saved as an image. Membrainy will also overlay the positions of molecules and ions on the maps. As iterative approaches can be computationally expensive, each leaflet is assigned its own thread allowing the maps to be generated in parallel.

Leaflet/membrane thickness, area per lipid and gel percentage

Membrane thickness is determined by calculating the average height of a user-specified reference atom, typically the phosphorous atom, for each leaflet. The average height of the reference atom for two opposing leaflets can then be subtracted. Leaflet thickness is calculated by subtracting the average height of the reference atom with the geometric centre of the bilayer. A 2D thickness map can also be produced by binning the reference atoms into a 2D lattice and applying the same algorithm used by the 2D surface maps. Membrainy offers a simple area per lipid (APL) calculation by dividing the box area by the number of lipids per leaflet, and will automatically produce multiple APLs for asymmetric bilayers or when lipid flip-flopping is detected. Gel percentages can be approximated by comparing the force field distance between the first and last carbon atoms in the lipid tails with the distance found in the trajectory files. As fluid lipid tails are non-linear, this distance is typically much less than the force field distance. A user-specified tolerance is assigned to the force field distance, and any lipid with a trajectory distance above this tolerance is counted as a ‘gel’ lipid.

Annular shell analysis

Membrainy isolates the annular shell of lipids around molecules by calculating a distance vector between each atom in the bilayer with each atom in the molecule. If the distance between any two atoms is within a user-specified radius, the lipid is counted as being within the shell. These lipids can then be analysed to determine their properties. A control group can also be established by selecting random lipids outside of the shell from the same leaflet, comprising either a fixed number of lipids, an identical number of lipids to those found within the shell or all lipids outside of the shell. An option exists to exclude gel lipids from the control group, as many proteins and peptides are known to show selectivity for inserting into fluid regions [36]. Gel lipids are identified using the same technique described above. If multiple molecules are present, the user may specify one, several or all molecules to construct annular shells for, and Membrainy will assign a thread to each molecule, populating the shells in parallel. The output plots contain an average of all shells in the system. Membrainy is also fitted with an annular shell analysis algorithm to produce detailed records of which lipids occupy the shell at any given time and which lipids spent the longest time in the shell. In mixed bilayer compositions, Membrainy will plot the ratio of lipid types found within the shell over time.

Evolution of the TMV

and is achieved by splitting the simulation box into ‘slices’ along the z-axis and calculating the charge density in each slice [37]. The box is then corrected such that Ψ(0)=0. Membrainy utilises the GROMACS tool g_potential by splitting the full trajectory into smaller trajectories and calculating the electrostatic potential in each trajectory. The TMV can then be extrapolated from each smaller trajectory and recombined to produce a voltage against time measurement over the full trajectory.

Lipid mixing/demixing entropy

where \(\rho _{x_{i}}\) is the density of a lipid of type x _i . A scaled entropy is also produced such that S _max =1.

Order parameters

Order parameters are a measure of the level of order or entropy in the lipid tails and can give insight into the fluidity of the membrane, as gel lipids exhibit a greater degree of order over fluid lipids. Membrainy was used to generate order parameters for various bilayer compositions, which found saturated lipid tails to be comparable to those produced by the GROMACS tool g_order in the CHARMM and GROMOS force fields. For bilayers employing the CHARMM force field, Membrainy yields more accurate order parameters to g_order by utilising each C −H bond vector, whereas g_order ignores these vectors and instead reconstructs them from the C _i−1−C _i+1 vector in a similar approach used by Membrainy with united-atom force fields. The order parameters for unsaturated lipids were also comparable, excluding the region around the double bond in which g_order calculates incorrectly. Membrainy automatically generates the necessary lipid tail information required to calculate order parameters, requiring no user-input. This is a significant improvement over g_order which requires a lengthy setup of user-constructed index files. This also restricts the order parameter analysis to a fixed number of lipids and as such it would be difficult to conduct an annular shell analysis with g_order.

Annular shell analysis

The properties of the lipids in close proximity to other molecules, whether inserted or in close contact with the membrane surface, may be modified by the presence of such molecules and result in local changes to the membrane. Membrainy isolates an annular shell of lipids surrounding a peptide, protein or other molecule, and compares its properties to lipids outside of the shell. As an exemplar, we inserted the MinD membrane targeting sequence (MinD-MTS) into the headgroup region of a POPE/POPG (3:1) double bilayer. This peptide is an 11-residue cationic amphipathic helix located at the C-terminus of the MinD protein that plays an important role in the cell division of Bacillus subtilis [39]. Using a shell radius of 4 Å, the order parameters of the lipids located within the shell were analysed over 50 ns and compared with a control group, comprising an identical number of lipids selected randomly from outside of the shell within the same leaflet. The order parameters reveal an increase in disorder for lipid tails within the shell when compared to those outside of the shell (Figure 1), which suggests the presence of splayed lipid tails. This phenomenon has previously been predicted for amphipathic peptides inserted into the headgroup region of lipid bilayers [26]. The option to ignore all gel lipids from the control group was enabled as our chosen peptide had inserted into a fluid region of the bilayer that contained ∼27% gel at 300 K. Without this option, the control group was observed to sample a more ordered phase of lipids, providing an inaccurate comparison with the fluid lipids found within the shell. The annular shell analysis algorithm revealed that lipids were continuously entering and exiting the shell, and saw approximately 10-14 lipids occupy the shell at any given time. Membrainy can also determine lipid type ratios within the shell for mixed bilayer compositions, which may be useful when studying molecules that give rise to an enrichment of certain lipid types. This measurement revealed an average lipid ratio of 2.6:1 POPE:POPG within the shell, suggesting that MinD-MTS gives rise to an enrichment of POPG lipids, likely mediated by the increased electrostatics between the cationic peptide and anionic POPG headgroups.

Evolution of the TMV

The TMV is the electrical potential found across biological cell membranes and plays a crucial role in a wide range of cellular processes, including the transport of nutrients into and out of the cell, biophysical signaling, and cell proliferation [40-42]. Membrainy is capable of measuring the time evolution of the TMV across a double bilayer, which may be of importance in electrophysiology and electroporation simulations. As an exemplar, numerous electroporation simulations were conducted for 30 ns using POPE/POPG (3:1) double bilayers. These systems were initially established with ion imbalances of +20, achieved by moving 10 cations from the inner (anodic) water compartment to the outer (cathodic) water compartment, similar to the approach taken by Sachs et al. [43]. Transient water pores were observed to form after a random time interval, allowing both cations and anions to travel through the pores in opposite directions, resulting in a loss of the initial ion imbalance. Membrainy was used to produce TMV against time measurements for each simulation, one of which is depicted in Figure 2. This measurement reveals an initial TMV of -2.65 V, which lowers to -2.35 V during the first 5 ns due to the lateral expansion of the bilayers undergoing electrostriction. Once a pore had formed, a sharp drop in TMV is observed at a rate of 0.75 V/ns, corresponding to ion transport through the pore at a rate of 3 ions/ns. The resulting TMV is indicative of the remaining ion imbalance (+2). These measurements provide an informative way to monitor changes to the TMV during a trajectory, and can be used in electroporation simulations to determine the time at which a pore is formed and the rate at which the TMV is dissipated.

Leaflet/membrane thickness and area per lipid

Using the same electroporation simulation as above, Membrainy was used to measure the leaflet and membrane thickness over the trajectory (with the membrane thickness depicted in Figure 2). This membrane thickness steadily decreases prior to pore formation, corresponding to the electrostriction effects experienced by the bilayers from the TMV. After pore formation, the bilayer thickness sharply increases as the TMV is dissipated, suggesting that the electrostriction effects had diminished and the bilayers were able to relax towards their initial thickness. Similar results were obtained for the leaflet thickness, and interestingly the anodic leaflets were observed to be consistently thinner than the cathodic leaflets prior to pore formation, which was also observed by Böckmann et al. in a similar electroporation study [44]. Membrainy was used to calculate the APL during the simulation, which revealed an increase in APL prior to pore formation, and a decrease afterwards. This result was identical to the APL produced with the GROMACS tool g_energy, which can output the box dimensions over time to be converted to an APL.

Lipid flip-flopping

Transmembrane lipid translocation, more commonly known as lipid flip-flopping, is the process in which lipids are translocated between the two opposing leaflets of a bilayer [45]. This translocation occurs from both passive and active transport mechanisms and plays a crucial role in the maintenance of asymmetric cell membranes [46]. Lipid flip-flopping has also been observed in simulated DMPC bilayers through electroporation, whereby lipids translocate through transient water pores in both directions [25].

Membrainy was used to detect lipid flip-flopping during a 30 ns simulation of a POPE/POPG (3:1) double bilayer, in which an ion imbalance of +20 was maintained with position restraints. A pore formed within 5 ns and remained open for the duration of the simulation. Figure 3 depicts the TMV and leaflet symmetry measurements during the simulation, where the leaflet symmetry is calculated by subtracting the number of lipids in the cathodic leaflets from the number of lipids in the anodic leaflets, and therefore a value of -2 indicates a single flip-flop to the cathodic leaflet. Upon pore formation, the leaflet symmetry reveals that the toroidal structure of the pore mainly comprised POPE and POPG lipids from the anodic leaflet. This is likely due to the tendency for transient water pores to initiate formation from the anodic water compartments, as observed by Böckmann et al. [44]. After 15 ns, the POPE symmetry returns to zero, indicating that the distribution of POPE lipids across both leaflets has equalised; however, the POPG symmetry steadily decreases, indicating that POPG lipids are translocating through the pore towards the cathodic leaflet. By 30 ns, one POPE lipid in each leaflet had flip-flopped and five POPG lipids had flip-flopped to the cathodic leaflet. This suggests that POPG lipids experience a greater tendency to flip-flop through transient water pores in bilayers subject to a high voltage TMV, in which the POPG lipids are translocated towards the cathodic leaflet, likely due to the additional forces acting on the anionic POPG headgroups from the electric field. Membrainy has therefore detected and interpreted lipid flip-flopping through a transient pore within this system.

2D surface maps and gel percentage

Surface maps were generated for POPE/POPG (3:1) bilayers at 297 K, 300 K and 320 K, and DPPC and POPC bilayers at 297 K using the CHARMM force field (Figure 4). The 300 K POPE/POPG bilayer contained the MinD-MTS inserted into the headgroup region, and the 297 K POPE/POPG bilayer was imaged before and during electroporation. As 297 K is approximately the transition temperature for POPE/POPG bilayers [47], Membrainy detected ∼53% gel within the bilayer. This percentage is accurately portrayed in the surface map (Figure 4a) where approximately half of the map appears as gel, represented by hexagonally packed red dots (the hexagonal packing of lipid tails occurs naturally in gel domains). The same bilayer was analysed at 320 K and showed ∼14% gel, which is also portrayed in the surface map (Figure 4b) by showing fewer gel clusters. The POPC bilayer at 297 K shows a highly fluid bilayer with ∼16% gel (Figure 4c), whereas the DPPC bilayer at 297 K shows ∼85% gel (Figure 4d). These measurements are in agreement with their corresponding transition temperatures of 271 K for POPC bilayers and 314 K for DPPC bilayers [48]. Interestingly, the DPPC bilayer exhibits gel lipids in the tilted \(\phantom {\dot {i}\!}L_{\beta ^{\prime }}\) phase which is portrayed in the surface map by the smudged appearance of the gel clusters. Finally, surface maps were generated for the bilayer containing an inserted MinD-MTS peptide (Figure 4e) and the bilayer from the previous electroporation simulation containing a transient water pore (Figure 4f). These 2D surface maps provide an alternative representation of the bilayer, capable of producing both a contour and density map that portrays gel clusters, pores, surface undulations and defects that may not be easily seen in 3D visualisation software such as VMD [49] and Pymol [50].

Headgroup orientation

The lipid headgroup is the polar interface between the membrane core and the intracellular/extracellular spaces and has been observed to exhibit sensitivity to electric charges, dipole fields, and temperature effects [51,52]. Experimental techniques have shown the lipid headgroup to sit roughly perpendicular to the lipid tails with a variation of around 30 degrees to the membrane surface [19,53]. In MD simulations, measurements of lipid headgroup orientations can provide an effective means to compare bilayers undergoing environmental perturbations, such as those under the influence of a TMV. As an exemplar, a POPC double bilayer was equilibrated for 100 ns without a TMV, after which a 30 ns simulation was conducted using an ion imbalance of +28, achieving a TMV of -1.95 V. Membrainy was used to measure the headgroup orientations before and after applying a TMV, which yielded a mean angle of 23 degrees in both leaflets without a TMV, and mean angles of 25.5 and 21.5 degrees in the anodic and cathodic leaflets respectively after applying a TMV. This reveals a shift of +2.5 degrees in the anodic leaflets and -1.5 degrees in the cathodic leaflets, suggesting that the headgroups in both leaflets are tending to align with the electric field. These measurements are comparable to those found by Böckmann et al. [44].

Lipid mixing/demixing entropy

The entropy of lipid mixing/demixing provides a measure of the two-dimensional spatial heterogeneity of any lipid bilayer system, and a means to study changes following an environmental perturbation. A bilayer was constructed containing 1512 POPE and 504 POPC Martini lipids, where the POPC lipids were initially clustered together in a quadrant of the bilayer creating a perfectly demixed system. This bilayer was simulated for 200 ns and its trajectory was analysed by Membrainy to produce a scaled mixing entropy over time (Figure 5). These measurements reveal an initial entropy of 0.3, which immediately increases as the system began to mix. By ∼150 ns the entropy settles just below the theoretical maximum entropy, indicating the bilayer was fully mixed. Membrainy has therefore quantified the level of mixing/demixing in this system.