Membrainy has been written in Java, which provides maximum compatibility across a range of operating systems, requires no compilation and enables the safe and efficient execution of multithreaded code. Membrainy contains various multithreaded algorithms to optimise efficiency and processor use across a range of architectures. These include algorithms for using multiple threads to load larger trajectory files, for preloading the next frame in the trajectory while the current frame is being analysed, and for running each analytical technique in parallel. Membrainy has been primarily designed for use with the GROMACS MD package [13], and contains a user interface that should be intuitive to GROMACS users. Membrainy is capable of reading GROMACS xtc, trr, tpr, cpt and gro trajectory file types, along with the standard pdb trajectory file type used by other MD packages (e.g. AMBER [27], CHARMM [28], NAMD [29], etc.). Membrainy has been implemented with the CHARMM36 [30], Berger/GROMOS87 [31] and Martini v2.0 [32] force fields, and is expandable to include other force fields and trajectory formats. Asymmetric bilayers and lipid flip-flops are detected by assigning each lipid to a corresponding leaflet depending on the height of its phosphorous atom relative to the geometric centre of the bilayer. All output graphs are readable by the Grace plotting software [33] and are preprogrammed with appropriate axis labels and titles. Double bilayer systems are automatically detected and incur additional output plots which contain averages of the inner and outer leaflets for certain analytical techniques.
Order parameters
Order parameters for saturated and unsaturated lipid tails in atomistic force fields are calculated from the equation
$$ S_{CD} = \left \langle \frac{3cos^{2}\theta - 1}{2} \right \rangle $$
((1))
where θ is the angle the C −H bond vectors along the lipid tails make with the membrane normal [34], taken as the z-axis for planar bilayers. This approach utilises each individual C −H bond in the lipid tails. As united-atom force fields lack non-polar hydrogen atoms, the above equation is modified to include the relation
$$ S_{CD} = \frac{2}{3}S_{xx} + \frac{1}{3}S_{yy} $$
((2))
which is derived from the order parameter tensor [35], and achieved by defining molecular axes where the z-axis encompasses the C
i−1−C
i+1 vector, the y-axis lies on the plane containing C
i−1−C
i
−C
i+1, and the x-axis is orthogonal to the y and z axes. The angles that the x and y axes make with the membrane normal is then used to determine S
xx
and S
yy
from Equation 1. Martini order parameters are calculated from the equation
$$ P_{2} = \frac{1}{2} \left (3\:cos^{2} \left \langle \theta \right \rangle -1 \right) $$
((3))
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
2=1 and disorder at P
2=−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
The membrane surface can be represented in a 2D map by binning the heights of each atom in each leaflet into a 2D lattice and applying the Gauss-Seidel method
$$ \phi_{i,j}^{n+1} = -\frac{1}{4}\left [ A_{i,j} - \left (\phi_{i-1,j}^{n} + \phi_{i+1,j}^{n} + \phi_{i,j-1}^{n} + \phi_{i,j+1}^{n} \right)\right ] $$
((4))
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
In double bilayer systems, the TMV can be extrapolated from the average electrostatic potential between the two bilayers, which is calculated from a double integral of Poisson’s equation
$$ \Psi (z) = - \frac{1}{\varepsilon_{0}}{\int_{0}^{z}}dz^{\prime} \int_{0}^{z^{\prime}} \rho \left(z^{\prime\prime}\right)dz^{\prime\prime} $$
((5))
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
Membranes containing two or more lipid types can have their lipid mixing/demixing quantified as an entropy with the equation
$$ S(x_{1},\!..,x_{N}) = N \sum\limits_{x_{i},nb_{i}} p(x_{i},nb_{i})\: log\: p(x_{i}\mid nb_{i}) $$
((6))
as described by Brandani et al. [38], where p(x
i
,n
b
i
) is the probability to find a lipid of type x
i
neighboured to a lipid of type n
b
i
, and p(x
i
∣n
b
i
) indicates the conditional probability that a lipid is of type x
i
given that its neighbour is of type n
b
i
. To calculate the entropy, a distance vector is established between the phosphorous atoms on each lipid in a leaflet to determine the nearest neighbouring lipid and its type. This information is then binned into a probability matrix and normalised such that the total probability is always 1, and then used with Equation 6 to produce an entropy. A theoretical maximum entropy can be calculated from
$$ S_{max} = -\sum \rho_{x_{i}} \: log \: \rho_{x_{i}} $$
((7))
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.
