Usage

This page includes a practical guide to using MembraneCurvature with the surface derivation method of choice (binning or Fourier). It also includes examples showing how curvature can be calculated for three different simulation systems.

1. Surface derivation methods

2. Examples of how to use MembraneCurvature to derive curvature profiles

2.1. Membrane-only systems

2.2 Membrane-protein systems

Warning

Examples included in this page show how to use MembraneCurvature using data files from MDAnalysisTests. To run these examples, MDAnalysisTests must be installed.

1. Surface derivation methods

MembraneCurvature uses an AtomGroup as a reference, user-defined via the select parameter, to derive a surface and calculate mean and Gaussian curvature.

There are two methods available to derive the surface:

  • Fourier (surface_method='fourier', default method) fits a truncated periodic 2D Fourier sum to atom heights by linear least squares at each frame, evaluates the fitted surface, and obtains partial derivatives analytically from that sum (no finite-difference on the grid). Optional arguments fourier_m, fourier_n, tune the truncation for the Fourier fit and the least-squares solve.

  • Binning (surface_method='binning') assigns atoms to a regular n_x_bins x n_y_bins grid, stores the mean \(z\) per cell, and estimates partial derivatives with numpy.gradient() using the physical bin spacing.

Warning

Use fourier_m = fourier_n = 2 (the constructor with default values) unless you need shorter-wavelength structure; increase fourier_m and fourier_n only while curvature improves systematically, rather than becoming noisier.

For copy-paste examples of both methods see:

2. Examples of how to use MembraneCurvature to derive curvature profiles

The following sections offer examples of how to use MembraneCurvature to derive curvature profiles in three types of systems:

2.1. Membrane-only systems

In this example, we show a basic usage of MembraneCurvature in a system that comprises a lipid bilayer of DPPC:CHOL using the Martini force field. Since we have a bilayer, we select atoms of phospholipid head groups in the upper leaflet only using the select parameter and apply coordinate wrapping.

2.1.1 Fourier surface method (default)

We can calculate membrane curvature using the Fourier surface method by either setting surface_method='fourier' with fourier_m=2, fourier_n=2, or omitting surface_method, fourier_m, and fourier_n to rely on the defaults (default fourier_m=2, fourier_n=2):

import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from MDAnalysis.tests.datafiles import Martini_membrane_gro

universe = mda.Universe(Martini_membrane_gro)

curvature_upper_leaflet = MembraneCurvature(universe,
                                            select='resid 1-225 and name PO4',
                                            ).run()

mean_upper_leaflet = curvature_upper_leaflet.results.average_mean

gaussian_upper_leaflet = curvature_upper_leaflet.results.average_gaussian

Note the code to calculate curvature for the upper leaflet with Fourier (default method) is equivalent to:

curvature_upper_leaflet = MembraneCurvature(universe,
                                            select='resid 1-225 and name PO4',
                                            surface_method='fourier',
                                            fourier_m=2,
                                            fourier_n=2,
                                            ).run()

Tip

When using the Fourier method, wrap is not required. Periodic boundary conditions are handled inside the Fourier fit. 🙂 Use wrap=True only with surface_method='binning' to pack atoms into the primary unit cell on raw trajectories.

Advanced: Tuning the Fourier least-squares cutoff (fourier_rcond)

The Fourier surface is fit by solving a linear least-squares system with singular-value truncation SVD. The optional cutoff fourier_rcond controls which singular values are treated as “effectively zero” and therefore removed from the solve. Smaller values keep more directions, and potentially noisier if the system is underdetermined. Larger values regularize more aggressively.

In MembraneCurvature, pass this cutoff as fourier_rcond:

curvature_upper_leaflet = MembraneCurvature(universe,
                                            select='resid 1-225 and name PO4',
                                            surface_method='fourier',
                                            fourier_m=2,
                                            fourier_n=2,
                                            fourier_rcond=1e-12
                                            ).run()

Note that fourier_rcond is a relative singular-value cutoff: singular values \(s \le rcond \cdot s_{\max}\) are treated as zero. If the effective rank of the design matrix is smaller than the number of fitted parameters, a UserWarning is emitted: the solver still returns a well-defined minimum-norm least-squares solution, but the coefficients are not uniquely determined by the data.

Warning

fourier_rcond controls how aggressively we ignore poorly constrained combinations of Fourier coefficients. We strongly recommend using fourier_rcond with its default value None. Larger values keep fewer singular-value directions (more stable / more regularized). Smaller values keep more directions (closer fit but potentially noisier).

Rough intuition:

  • fourier_rcond=None: sensible default; uses NumPy’s heuristic cutoff.

  • fourier_rcond=0: truncate only exactly zero singular values.

  • fourier_rcond=1e-12 or 1e-10: more aggressive truncation; can reduce noise when the fit is underdetermined.

Note

All Fourier least-squares steps (design matrix, SVD, coefficients) use 64-bit floating point numpy.float64(). The cutoff defined by fourier_rcond is a relative threshold: singular values with \(s \le rcond \cdot s_{\max}\) are dropped. The meaningful scale is therefore relative to the largest singular value \(s_{\max}\), not absolute coordinates or heights.

With fourier_rcond=None, the cutoff scales with the size of the least-squares problem, typically the number of atoms in the selection, or the number of fitted coefficients if that is larger. Passing a value much smaller than \(\sim 10^{-16}\), or much smaller than that automatic cutoff, usually has no visible effect. To smooth an underdetermined fit on purpose, use larger values such as 1e-12 or 1e-10 (see warning above).

2.1.2 Binning surface method

Alternatively, set surface_method='binning' and provide the values for the binning grid n_x_bins and n_y_bins. Note that you need to apply coordinate wrapping with wrap=True:

import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from MDAnalysis.tests.datafiles import Martini_membrane_gro

universe = mda.Universe(Martini_membrane_gro)

curvature_upper_leaflet = MembraneCurvature(universe,
                                            select='resid 1-225 and name PO4',
                                            surface_method='binning',
                                            n_x_bins=8,
                                            n_y_bins=8,
                                            wrap=True
                                            ).run()

mean_upper_leaflet = curvature_upper_leaflet.results.average_mean

gaussian_upper_leaflet = curvature_upper_leaflet.results.average_gaussian

By default, binning uses fft_filter='auto', so the average_z_surface map above is computed from the FFT-filtered time-averaged height (see 2.1.2.1 Binning with FFT filtering).

You can find more detailed examples in the notebooks available in the Tutorials page.

2.1.2.1 Binning with FFT filtering

Warning

FFT filtering is available only with surface_method='binning'. It does not filter each frame.

The order of operations is:

  1. Per frame: bin atoms, store height field in results.z_surface.

  2. After the trajectory is processed, time-average z_surface over frames, optionally apply one brick-wall filter in reciprocal space to that average, then compute results.average_mean and results.average_gaussian from the (possibly filtered) average height.

Binning analyses use fft_filter='auto' by default (low-pass (0, 0.5 * q_Nyq) from dx and dy). Pass fft_filter=None to disable filtering on the average map.

curvature_upper_leaflet = MembraneCurvature(universe,
                                            select='resid 1-225 and name PO4',
                                            surface_method='binning',
                                            n_x_bins=8,
                                            n_y_bins=8,
                                            wrap=True,
                                            fft_filter={'q': (0, 0.5 * q_Nyq)}
                                            ).run()

Warning

As shown above, MembraneCurvature allows custom values for the FFT filtering by passing a tuple of (q_low, q_high) in rad/Å to fft_filter={'q': (q_low, q_high)}. However, this is not recommended! Custom values should be used with caution.

The recommended way is to use the automatic mode fft_filter='auto' with the default low-pass (0, 0.5 * q_Nyq) from dx and dy.

2.2 Membrane-protein systems

Tip

To improve sampling when passing raw trajectories:

  • In systems of membrane-only or membrane-protein with position restraints, set wrap=True to translate the atoms of the AtomGroup back in the unit cell.

  • In membrane-protein systems with no position restraints, set wrap=False and preprocess the trajectory with rotational/translational fit.

Some points to keep in mind when calculating membrane curvature in membrane-only and membrane-protein systems with position restraints are addressed in this blog post.

2.2.1 Membrane-protein systems, protein with position restraints

In this example, we have a simulation box comprising a copy of the Yiip transporter, embedded in a lipid bilayer of POPE:POPG. Similar to the example for membrane-only, we select the atoms for the upper leaflet to run the analysis.

Fourier surface method (default)

We can calculate membrane curvature using the default values with:

import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from MDAnalysis.tests.datafiles import XTC_MEMPROT, GRO_MEMPROT

universe = mda.Universe(GRO_MEMPROT, XTC_MEMPROT)
curvature_upper_leaflet = MembraneCurvature(universe,
                                            select='resid 297-517 and name P'
                                            ).run()

avg_mean_curvature_upper_leaflet = curvature_upper_leaflet.results.average_mean

avg_gaussian_curvature_upper_leaflet = curvature_upper_leaflet.results.average_gaussian

Binning surface method

The same trajectory and selection can use the binning method by setting surface_method='binning' with the values for n_x_bins and n_y_bins and apply coordinate wrapping with wrap=True.

import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from MDAnalysis.tests.datafiles import XTC_MEMPROT, GRO_MEMPROT

universe = mda.Universe(GRO_MEMPROT, XTC_MEMPROT)
curvature_upper_leaflet = MembraneCurvature(universe,
                                            select='resid 297-517 and name P',
                                            surface_method='binning',
                                            n_x_bins=2,
                                            n_y_bins=2,
                                            wrap=True
                                            ).run()

avg_mean_curvature_upper_leaflet = curvature_upper_leaflet.results.average_mean

avg_gaussian_curvature_upper_leaflet = curvature_upper_leaflet.results.average_gaussian

2.2.2. Membrane-protein systems, protein with no position restraints

For membrane-protein systems where the simulation setup has no position restraints on the protein, a trajectory preprocessing by the user is required. If the goal is to assess membrane curvature induced by the protein, the preprocessed trajectory should have the protein centered in the simulation box with translational and rotational fit.

for example, in Gromacs, the trajectory would be preprocessed with:

gmx trjconv -pbc whole -ur compact -c
gmx trjconv -fit rot+trans
gmx trjconv -fit transxy

After you have preprocessed the trajectory, use wrap=False (the trajectory is already fitted and centered).

Fourier surface method (default)

With the default surface_method='fourier' and fourier_m=2, fourier_n=2. Omit surface_method, fourier_m, and fourier_n unless you need shorter wavelengths:

import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from membrane_curvature.tests.datafiles import XTC_MEMBPROT_FIT, GRO_MEMBPROT_FIT

universe = mda.Universe(GRO_MEMBPROT_FIT, XTC_MEMBPROT_FIT)

curvature_lower_leaflet = MembraneCurvature(universe,
                                            select='resid 2583-3042'
                                            ).run()

avg_mean_curvature = curvature_lower_leaflet.results.average_mean

avg_gaussian_curvature = curvature_lower_leaflet.results.average_gaussian

Binning surface method

Set surface_method='binning' with the values for n_x_bins and n_y_bins, and in this case set wrap=False to avoid the warning message since the trajectory is already fitted and centered:

import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from membrane_curvature.tests.datafiles import XTC_MEMBPROT_FIT, GRO_MEMBPROT_FIT

universe = mda.Universe(GRO_MEMBPROT_FIT, XTC_MEMBPROT_FIT)

curvature_lower_leaflet = MembraneCurvature(universe,
                                            select='resid 2583-3042',
                                            surface_method='binning',
                                            n_x_bins=10,
                                            n_y_bins=10,
                                            wrap=False
                                            ).run()

avg_mean_curvature = curvature_lower_leaflet.results.average_mean

avg_gaussian_curvature = curvature_lower_leaflet.results.average_gaussian

Note

Since you are providing a preprocessed trajectory with translation/rotational fit you can ignore the warning message: WARNING   `wrap == False` may result in inaccurate calculation of membrane curvature.

More information on how to visualize the results of the MDAnalysis Membrane Curvature tool can be found in the Visualization page.