This page provides examples of how to visualize membrane curvature via Matplotlib. Two different approaches are suggested.

1. imshow

2. contourf

1. imshow

A simple plot using imshow() can be obtained like so:

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.imshow(avg_mean_curvature, cmap='bwr', interpolation='gaussian', origin='lower')
ax.set_title('Mean Curvature')

With imshow(), each element of the array is plotted as a square in a matrix of m x n elements. The color of each square is determined by the value of the corresponding array element and the colormap of preference.

For example, to visualize the results obtained in 1. Membrane-only systems, we can run:

In [1]: import MDAnalysis as mda
   ...: from membrane_curvature.base import MembraneCurvature
   ...: from MDAnalysis.tests.datafiles import Martini_membrane_gro
   ...: import matplotlib.pyplot as plt

In [2]: universe = mda.Universe(Martini_membrane_gro)

In [3]: curvature_upper_leaflet = MembraneCurvature(universe,
   ...:                                             select='resid 1-225 and name PO4',
   ...:                                             n_x_bins=8,
   ...:                                             n_y_bins=8,
   ...:                                             wrap=True).run()

In [4]: mean_upper_leaflet = curvature_upper_leaflet.results.average_mean

In [5]: curvature_lower_leaflet = MembraneCurvature(universe,
   ...:                                             select='resid 226-450 and name PO4',
   ...:                                             n_x_bins=8,
   ...:                                             n_y_bins=8,
   ...:                                             wrap=True).run()

In [6]: mean_lower_leaflet = curvature_lower_leaflet.results.average_mean

In [7]: leaflets = ['Lower', 'Upper']

In [8]: curvatures = [mean_lower_leaflet, mean_upper_leaflet]

In [9]: fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(6,3), dpi=200)
   ...: for ax, mc, lf in zip((ax1, ax2), curvatures, leaflets):
   ...:     ax.imshow(mc, interpolation='gaussian', cmap='seismic', origin='lower')
   ...:     ax.set_aspect('equal')
   ...:     ax.set_title('{} Leaflet'.format(lf))
   ...:     ax.axis('off')

2. contourf

You can use contour plots using contourf(). With this approach, contour lines and filled contours of the obtained two-dimensional data are plotted. A contour line connects points with the same curvature values.

When plotting using contourf(), an extra step is required to perform an interpolation. We suggest using scipy.ndimage.gaussian_filter() as in:

In [10]: from scipy import ndimage

In [11]: leaflets = ['Lower', 'Upper']

In [12]: fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(5,3))
   ....: for ax, mc, lf in zip((ax1, ax2), curvatures, leaflets):
   ....:     ax.contourf(ndimage.gaussian_filter(mc, sigma=1, order=0, mode='reflect'),
   ....:                 cmap='bwr',
   ....:                 levels=30)
   ....:     ax.set_aspect('equal')
   ....:     ax.set_title('{} Leaflet'.format(lf))
   ....:     ax.axis('off')