Overview of analysis tools¶

Here we provide a brief description of the analysis tools currently available in lipyphilic. For more information on each analysis tool, including details of all optional input parameters see the API. To learn more about how to use lipyphilic, check out our interactive tutorials.

Assign leaflets: `lipyphilic.lib.assign_leaflets`¶

This module provides methods for assigning lipids to leaflets in a bilayer. Leaflet assignment is based on the distance in z from a lipid the midpoint of the bilayer. Lipids may be assigned to the upper leaflet (indicated by 1), the lower leaflet (-1) or the bilayer midplane (0).

Below we see how to assign lipids to the upper or lower leaflet of a MARTINI bilayer:

import MDAnalysis as mda
from lipyphilic.lib.assign_leaflets import AssignLeaflets

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

# Find which leaflet each lipid is in at each frame
leaflets = AssignLeaflets(
  universe=u,
  lipid_sel="name PO4 ROH"
)

# Select which frames to use and perform the analysis
leaflets.run(start=None, stop=None, step=None)  # this will use every frame in the trajectory

The results are stored as a NumPy array of shape (n_lipids, n_frames) in the leaflets.leaflets attribute.

If you have used a different force field, you simply need to change the lipid_sel to select the relevant headgroup atoms of your lipids. See the MDAnalysis selection language for more info on how to select atoms.

By default, lipids are only allowed to be in the upper (1) or lower (-1) leaflet. See lipyphilic.lib.assign_leaflets for more information on selecting which molecules are allowed in the midplane.

Note

Assignment of lipids to leaflets is not in itself useful, but it is required in order to calculate, for example, area per lipid, interleaflet correlations, and flip-flop rates.

Flip-flop: `lipyphilic.lib.flip_flop`¶

This module provides methods for detecting the flip-flop of molecules in a lipid bilayer. A flip-flop occurs when a molecule - typically a sterol - moves from one leaflet of a bilayer into the opposing leaflet.

To find all flip-flop events, we first should assign lipids to leaflets as seen in the above example, then:

import MDAnalysis as mda
from lipyphilic.lib.flip_flop import FlipFlop

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

flip_flops = FlipFlop(
  universe=u,
  lipid_sel="name ROH",  # select molecules that may flip-flop
  leaflets=leaflets.filter_leaflets("name ROH")
)

flip_flops.run(start=None, stop=None, step=None)

The results are stored as a NumPy array of shape (n_flip_flops, 4) in the flip_flops.flip_flops attribute. Each row is a single flip-flop event, and the four columns correspond to: the residue index of the flip-flopping molecule; the frame at which the molecule left its original leaflet; the frame at which it entered its new leaflet; the leaflet ID to which it moves.

See lipyphilic.lib.flip_flop for more information on how flip-flop is detected and options such as specifying how long a molecule must residue in the new leaflet for the flip-flop to be considered successful.

Interlealet registration: `lipyphilic.lib.registration`¶

This module provides methods for determining registration of leaflets in a bilayer. Registration is defined by the pearson correlation coefficient of molecular densities in the two leaflets. This is an implementation of the method described by Thallmair et al. (2018).

To calculate the interleaflet correlation of cholesterol, we first need to calculate which leaflet each lipid is in at each frame using lipyphilic.lib.assign_leaflets.AssignLeaflets. Then we pass atom selections for which density correlations will be calculated, along with the relevant leaflet membership data, to Registration:

import MDAnalysis as mda
from lipyphilic.lib.registration import Registration

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

registration = Registration(
  upper_sel="resname CHOL and name ROH",
  lower_sel="resname CHOL and name ROH",
  leaflets=leaflets.filter_leaflets("name ROH")
)

registration.run(start=None, stop=None, step=None)

The results are stored in a NumPy array of shape (n_frames), containing the pearson correlation coefficient of cholesterol densities in the two leaflets. The data are accessible via the registration.registration attribute.

As well as calculating registration of lipid species across the two leaflets, it is also possible to calculate the registration of arbitrary user-defined values across the two leaflets. For example, if you have created a Hidden Markov Model to assign lipids to the Ld or Lo phase, you can calculate the registration of Lo lipids across the two leaflets. See lipyphilic.lib.registration for more details.

Neighbours: `lipyphilic.lib.neighbours`¶

This module provides methods for finding neighbouring lipids in a bilayer. Lipids are neighbours if they are within a user-defined cutoff of one another.

Below we see how to find all neighbours in a MARTINI bilayer based on the ‘GL1’ and ‘GL2’ beads of phospholipids and the ‘ROH’ bead of sterols, using a cutoff of 12 Å:

import MDAnalysis as mda
from lipyphilic.lib.neighbours import Neighbours

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

# Find neighbouring lipids
neighbours = Neighbours(
  universe=u,
  lipid_sel="name GL1 GL2 ROH",
  cutoff=12.0
)

neighbours.run(start=None, stop=None, step=None)

The results are stored in the neighbours.neighbours attribute as a NumPy array of SciPy sparse matrices (of type scipy.sparse.csc_matrix). Each sparse matrix contains the lipid neighbours at a given frame.

Tip

Once the neighbour matrices has been generated, the local lipid compositions or the largest lipids cluster at each frame can be readily.

See lipyphilic.lib.neighbours for more information on this module, including how to calculate local lipid compositions or the lipid enrichment/depletion index, and how to find the largest cluster of a given lipid species over time.

Area per lipid: `lipyphilic.lib.area_per_lipid`¶

This module provides methods for calculating the area per lipid. Areas are calculated via a 2D Voronoi tessellation, using the locality module of Freud to perform the tessellation of atomic positions. See Lukat et al. (2013) a thorough description of calculating the area per lipid via Voronoi tessellations.

Once lipids have been assigned to leaflets, the area per lipid can be calculated as follows:

import MDAnalysis as mda
from lipyphilic.lib.area_per_lipid import AreaPerLipid

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

areas = AreaPerLipid(
  universe=u,
  lipid_sel="name GL1 GL2 ROH",  # assuming we're using the MARTINI forcefield
  leaflets=leaflets.leaflets
)

areas.run(start=None, stop=None, step=None)

The above will use GL1 and GL2 beads to calculate the area of each phospholipid, and the ROH bead to calculate the area of each sterol.

For a more complete description of calculating the area per lipid, and the API of the analysis class, see lipyphilic.lib.area_per_lipid.

Lipid order parameter — `lipyphilic.lib.order_parameter`¶

This module provides methods for calculating the coarse-grained orientational order parameter of acyl tails in a lipid bilayer. The coarse-grained order parameter, \(S_{CC}\), is a measure of the degree of ordering of an acyl tail, based on the extent to which the vector connecting two consecutive tail beads is aligned with the membrane normal.

See Seo et al. (2020) for a definition of \(S_{CC}\) and Piggot et al. (2017) for an excellent discussion on acyl tail order parameters in molecular dynamics simulations.

To calculate \(S_{CC}\), we need to provide an atom selection for the beads in a single tail of lipids in the bilayer — that is, either the sn1 or sn2 tails, not both. If we have performed a MARTINI simulation, we can calculate the \(S_{CC}\) of all sn1 tails of phospholipids as follows:

import MDAnalysis as mda
from lipyphilic.lib.order_parameter import SCC

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

scc = SCC(
  universe=u,
  tail_sel="name ??A"
)

The above makes use of the powerful MDAnalysis selection language. It will select beads such as C1A, C2A, D2A etc. This makes it simple to quickly calculate \(S_{CC}\) for the sn1 tails of all species in a bilayer.

To see how to calculate \(S_{CC}\) using local membrane normals to define the molecular axes, as well as the full API of the class, see lipyphilic.lib.order_parameter.

Lipid \(z\) angles: `lipyphilic.lib.z_angles`¶

This module provides methods for calculating the angle lipids make with the positive \(z\) axis. If we define the orientation of MARTINI cholesterol as the angle between the \(z\)-axis and the vector from the the ‘R5’ bead to the ‘ROH’ bead, we can calculate the orientation of each cholesterol molecule as follows:

import MDAnalysis as mda
from lipyphilic.lib.z_angles import ZAngles

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

z_angles = ZAngles(
  universe=u,
  atom_A_sel="name R5",
  atom_B_sel="name ROH"
)

z_angles.run(start=None, stop=None, step=None)

The results are stored in a numpy.ndarray of shape (n_residues, n_lipids) in the z_angles.z_angles attribute.

For more information on this module, including how to return the angles in radians rather than degrees, see lipyphilic.lib.z_angles.

Lipid \(z\) positions: `lipyphilic.lib.z_positions`¶

This module provides methods for calculating the height in \(z\) of lipids from the bilayer center.

If we have used the MARTINI forcefield to study a phospholipid/cholesterol mixture, we can calculate the height of cholesterol in the bilayer as follows:

import MDAnalysis as mda
from lipyphilic.lib.z_positions import ZPositions

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

z_positions = ZPositions(
  universe=u,
  lipid_sel="name GL1 GL2 ROH",
  height_sel="name ROH",
  n_bins=10
)

z_positions.run(start=None, stop=None, step=None)

lipid_sel is an atom selection that covers all lipids in the bilayer. This is used for calculating the membrane midpoint. height_sel selects which atoms to use for calculating the height of each lipid.

Local membrane midpoints are calculated by creating a grid of membrane patches, with the number of grid points controlled with the n_bins parameter. The distance in \(z\) of each lipid to its local midpoint is then calculated.

Data are returned in a numpy.ndarray of shape (n_residues, n_frames). See lipyphilic.lib.z_positions for more information on this module including the full API of the class.

Lipid \(z\) thickness: `lipyphilic.lib.z_thickness`¶

This module provides methods for calculating the thickness, in \(z\), of lipid tails. This is defined as the maximum distance in \(z\) between to atoms in a tail.

If we have used the MARTINI forcefield to study a DPPC/DOPC/cholesterol mixture, we can calculate the thickness of DPPC and DOPC sn1 tails, as well as the thickness of cholesterol, as follows:

import MDAnalysis as mda
from lipyphilic.lib.z_positions import ZThickness

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

z_thickness = ZThickness(
  universe=u,
  lipid_sel="(name ??1 ??A) or (resname CHOL and not name ROH)"
)

z_thickness.run()

The above makes use of the powerful MDAnalysis atom selection language to select the DPPC and DOPC sn1 tails along with cholesterol.

The thickness data are stored in a numpy.ndarray of shape (n_residues, n_frames) in the z_thickness.z_thickness attribute. See lipyphilic.lib.z_thickness for the full API of the class.

Membrane \(z\) thickness: `lipyphilic.lib.memb_thickness`¶

This module provides methods for calculating the bilayer thickness. It is defined as the peak-to-peak distance of lipid headgroup density in \(z\).

Lipids must first be assigned to the upper and lower leaflets. This can be done with the class lipyphilic.lib.assign_leaflets.AssignLeaflets. Then, to calculate the membrane thickness we need to define which atoms to treat as headgroup atoms and pass the leaflet membership information to MembThickness. If we have studied a DPPC/DOPC/cholesterol mixture with MARTINI, we could calculate the membrane thickness as follows:

import MDAnalysis as mda
from lipyphilic.lib.z_positions import ZThickness

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

memb_thickness = MembThickness(
  universe=u,
  leaflets=leaflets.filter_leaflets("resname DOPC and DPPC"),  # exclude cholesterol from thickness calculation
  lipid_sel="resname DPPC DOPC and name PO4"
)

memb_thickness.run()

The results are then available in the memb_thickness.memb_thickness attribute as a numpy.ndarray.

For more information on calculating membrane thickness, including options to calculating local membrane thicknesses rather than a single global thickness, see lipyphilic.lib.memb_thickness.

Lateral diffusion `lipyphilic.lib.lateral_diffusion`¶

This module contains methods for calculating the mean squared displacement (MSD) and lateral diffusion coefficient, \(D_{xy}\),of lipids in a bilayer.

The MSD of all lipids in a DPPC/DOPC/cholesterol MARTINI bilayer can be calculated using lipyphilic.lib.lateral_diffusion.MSD:

import MDAnalysis as mda
from lipyphilic.lib.lateral_diffusion import MSD

# Load an MDAnalysis Universe
u = mda.Universe('production.tpr','production.xtc')

msd = MSD(
  universe=u,
  lipid_sel="name PO4 ROH"
)

 msd.run()

The MSD of each lipid is then available in the msd.msd attribute as a numpy.ndarray, and the lagtimes are stored in the msd.lagtimes attribute.

For more information on this module, including how to calculate the lateral diffusion coefficient, see lipyphilic.lib.lateral_diffusion.

Plotting utilities: `lipyphilic.lib.plotting`¶

lipyphilic can produce joint probability density plots (or PMFs if a temperature is provided), as well as density maps of membrane properties projected onto the membrane plane. The former may be used to plot, for example, the PMF of cholesterol orientation and height in a bilayer. The latter may be used to generate plots of, for example, the area per lipid as a function of \(xy\) in the membrane plane.

See lipyphilic.lib.plotting for the full API of lipyphilic.lib.plotting.JointDensity and lipyphilic.lib.plotting.ProjectionPlot.

On-the-fly transformations `lipyphilic.transformations`¶

lipyphilic contains a module for applying on-the-fly transformation to atomic coordinates while iterating over a trajectory. These are available in the module lipyphilic.transformations.

There are three transformations available in lipyphilic:

lipyphilic.transformations.nojump, which prevents atoms from jumping across periodic

boundaries. This is useful when calculating the lateral diffusion of lipids.
lipyphilic.transformations.center_membrane, which can take a membrane that is split

across periodic boundaries, make it whole and center it in the box.
lipyphilic.transformations.triclinic_to_orthorhombic, which transforms triclinic coordinates

into their orthorhombic representation.

See lipyphilic.transformations for full details on these transformations including how to apply them to your trajectory.

Overview of analysis tools¶

Assign leaflets: lipyphilic.lib.assign_leaflets¶

Flip-flop: lipyphilic.lib.flip_flop¶

Interlealet registration: lipyphilic.lib.registration¶

Neighbours: lipyphilic.lib.neighbours¶

Area per lipid: lipyphilic.lib.area_per_lipid¶

Lipid order parameter — lipyphilic.lib.order_parameter¶

Lipid \(z\) angles: lipyphilic.lib.z_angles¶

Lipid \(z\) positions: lipyphilic.lib.z_positions¶

Lipid \(z\) thickness: lipyphilic.lib.z_thickness¶

Membrane \(z\) thickness: lipyphilic.lib.memb_thickness¶

Lateral diffusion lipyphilic.lib.lateral_diffusion¶

Plotting utilities: lipyphilic.lib.plotting¶

On-the-fly transformations lipyphilic.transformations¶