Membrane thickness — lipyphilic.lib.memb_thickness


Paul Smith




GNU Public License v2

This module provides methods for calculating the membrane thickness over time.

The membrane thickness is defined as the mean distance in \(z\) between lipids in the upper and lower leaflets. A discrete intrinsic surface is constructed for each leaflet based on user-defined lipid headgroup atoms, and the mean distance in \(z\) between the two surfaces defines the membrane thickness.


  • universe : an MDAnalysis Universe object

  • leaflets : a NumPy array containing the leaflet membership of each lipid at each frame

  • lipid_sel : atom selection for lipids in the upper leaflet to used in the thickness calculation

Options: - n_bins : a discrete intrinsic surface of each leaflet is created with n_bins * n_bins patches


  • thickness : the mean membrane thickness at each frame

Thickness data are returned in a numpy.ndarray, containing the mean membrane thickness at each frame.

Example usage of MembThickness

An MDAnalysis Universe must first be created before using MembThickness:

import MDAnalysis as mda
from lipyphilic.lib.memb_thickness import MembThickness

u = mda.Universe(tpr, trajectory)

Then we need to know which leaflet each lipid is in at each frame. This may be done using lipyphilic.lib.assign_leaflets.AssignLeaflets:

leaflets = AssignLeaflets(
  lipid_sel="name GL1 GL2 ROH" # assuming we are using the MARTINI forcefield

The leaflets data are stored in the leaflets.leaflets attribute. We can now create our MembThickness object by passing the results of lipyphilic.lib.assign_leaflets.AssignLeaflets MembThickness along with an atom selection for the lipids:

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

To calculate the membrane thickness, based on interleaflet PO4 to PO4 distances, we need to use the run() method. We select which frames of the trajectory to analyse (None will use every frame) and choose to display a progress bar (verbose=True):

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

Changing the resolution of the 2D grid

By default, the lipid positions of each leaflet are binned into a two-dimensional histogram using 1 bins in each dimension. This is equivalent to calculating the mean height of all headgroup atoms in the bilayer, without discretising the surface.

It is also possible to specify the bins to use for binning the data:

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

This will use 10 bins in each dimension for creating the two-dimensional histogram.

Interpolate missing values in a grid with many bins

This is useful only if you would like a very high resolution grid. Having a higher resolution grid may be useful if you would like to later calculate, for example, the correlation between local membrane thicknesses and the local membrane area per lipid. The area per lipid can be projected onto the membrane plane using the class ProjectionPlot, and the height of the bilayer as a function of \(xy\) can be obtained from MembThickness by setting the return_surface keyword to True.

A grid with a small bin size (large n_bins) will lead to bins with no atom, and thus no height value. In this instance, the interpolate keyword should be set to True. However, interpolation substantially decreases performance and should be left as False unless it is strictly necessary.

The class and its methods

class lipyphilic.lib.memb_thickness.MembThickness(universe, lipid_sel, leaflets, n_bins=1, interpolate=False, return_surface=False)[source]

Calculate the bilayer thickness.

Set up parameters for calculating membrane thickness.

  • universe (Universe) – MDAnalysis Universe object

  • leaflets (numpy.ndarray (n_lipids,)) – An array of leaflet membership in which: -1 corresponds to the lower leaflet; 1 corresponds to the upper leaflet; and 0 corresponds to the midplane. If the array is 1D and of shape (n_lipids), each lipid is taken to remain in the same leaflet over the trajectory. If the array is 2D and of shape (n_lipids, n_frames), the leaflet to which each lipid is assisgned at each frame will be taken into account when calculating the area per lipid.

  • lipid_sel (str, optional) – Selection string for lipid atoms to be used in the thickness calculation. The default is None, in which case all atoms of the lipids will be used.

  • n_bins (int, optional) – Number of bins in x and y to use to create a grid of membrane patches. The intrinsic surface of a leaflet is constructed via the height in z of each patch. The default is 1, which is equivalent to computing a single global leaflet height.

  • interpolate (bool, optional) – If True, interpolate the two intrinsic surfaces to fill missing values. This substantially decreases performance but allows for the construction of higher-resolution grids. The default is False.

  • return_surface (bool, optional) – If True, the height of the bilayer at grid point at each frame is returned as numpy ndarray. The default is False.


Leaflet membership can be determined using lipyphilic.lib.assign_leaflets.AssignLeaflets.