Membrane thickness —
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
Example usage of
An MDAnalysis Universe must first be created before using
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
leaflets = AssignLeaflets( universe=u, lipid_sel="name GL1 GL2 ROH" # assuming we are using the MARTINI forcefield ) leaflets.run()
The leaflets data are stored in the
leaflets.leaflets attribute. We can now create our
MembThickness object by passing the results of
MembThickness along with an atom selection for the lipids:
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" )
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):
memb_thickness.run( start=None, stop=None, step=None, verbose=True )
The results are then available in the
memb_thickness.memb_thickness attribute as a
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( universe=u, leaflets=leaflets.filter_leaflets("resname DOPC and DPPC"), # exclude cholesterol from thickness calculation lipid_sel="resname DPPC DOPC and name PO4", n_bins=10 )
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
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)¶
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