NhaA membrane protein in a POPE:POPG lipid bilayer (all atom)
Last executed: May 16th, 2026 with MDAnalysis 2.10.0 and Python 3.13.5
Packages required: MDAnalysis, MDAnalysisData.
Packages for visualization: Matplotlib, SciPy, NGLview.
In this tutorial, we derive 2D curvature plots for a membrane–protein system containing a single copy of the sodium–proton antiporter NhaA embedded in a lipid bilayer with a 4:1 POPE:POPG composition. The system was obtained from an atomistic molecular dynamics (MD) simulation performed using the CHARMM36 force field
import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
from MDAnalysisData import datasets
import matplotlib.pyplot as plt
import nglview as nv
import numpy as np
from scipy import ndimage
%matplotlib inline
This tutorial is divided into six main steps:
1. Download dataset from MDAnalysisData
1. Download dataset from MDAnalysisData
Since the NhaA dataset is available via the MDAnalysisData collection, we need an extra step to fetch our dataset using the datasets.fetch_nhaa_equilibrium() function.
You can find more information on how to access files from MDAnalysisData here.
We retrieve our NhaA dataset with (a progress bar will show up here):
NhaA = datasets.fetch_nhaa_equilibrium()
NhaA_non_water.gro: 0.00B [00:00, ?B/s]
NhaA_non_water.gro: 0%| | 8.19k/4.19M [00:00<04:54, 14.2kB/s]
NhaA_non_water.gro: 2%|▏ | 73.7k/4.19M [00:00<00:32, 127kB/s]
NhaA_non_water.gro: 9%|▉ | 369k/4.19M [00:00<00:05, 671kB/s]
NhaA_non_water.gro: 19%|█▉ | 786k/4.19M [00:00<00:02, 1.38MB/s]
NhaA_non_water.gro: 77%|███████▋ | 3.21M/4.19M [00:01<00:00, 6.35MB/s]
NhaA_non_water.gro: 4.19MB [00:01, 3.70MB/s]
NhaA_non_water.xtc: 0.00B [00:00, ?B/s]
NhaA_non_water.xtc: 0%| | 8.19k/1.15G [00:00<22:37:56, 14.1kB/s]
NhaA_non_water.xtc: 0%| | 73.7k/1.15G [00:00<2:31:08, 126kB/s]
NhaA_non_water.xtc: 0%| | 270k/1.15G [00:00<39:46, 480kB/s]
NhaA_non_water.xtc: 0%| | 721k/1.15G [00:00<14:34, 1.31MB/s]
NhaA_non_water.xtc: 0%| | 2.33M/1.15G [00:01<04:12, 4.54MB/s]
NhaA_non_water.xtc: 1%| | 6.00M/1.15G [00:01<01:36, 11.8MB/s]
NhaA_non_water.xtc: 1%| | 12.7M/1.15G [00:01<00:49, 22.8MB/s]
NhaA_non_water.xtc: 2%|▏ | 18.7M/1.15G [00:01<00:35, 31.5MB/s]
NhaA_non_water.xtc: 2%|▏ | 24.2M/1.15G [00:01<00:30, 36.6MB/s]
NhaA_non_water.xtc: 2%|▏ | 28.6M/1.15G [00:01<00:33, 33.6MB/s]
NhaA_non_water.xtc: 3%|▎ | 36.5M/1.15G [00:01<00:28, 38.7MB/s]
NhaA_non_water.xtc: 4%|▍ | 44.4M/1.15G [00:02<00:26, 42.0MB/s]
NhaA_non_water.xtc: 5%|▍ | 52.3M/1.15G [00:02<00:24, 43.9MB/s]
NhaA_non_water.xtc: 5%|▌ | 59.2M/1.15G [00:02<00:24, 45.2MB/s]
NhaA_non_water.xtc: 6%|▌ | 63.2M/1.15G [00:02<00:24, 44.0MB/s]
NhaA_non_water.xtc: 6%|▌ | 71.1M/1.15G [00:02<00:23, 45.5MB/s]
NhaA_non_water.xtc: 7%|▋ | 79.1M/1.15G [00:02<00:22, 46.4MB/s]
NhaA_non_water.xtc: 8%|▊ | 87.0M/1.15G [00:02<00:22, 47.0MB/s]
NhaA_non_water.xtc: 8%|▊ | 94.9M/1.15G [00:03<00:22, 47.5MB/s]
NhaA_non_water.xtc: 9%|▉ | 103M/1.15G [00:03<00:21, 47.6MB/s]
NhaA_non_water.xtc: 10%|▉ | 111M/1.15G [00:03<00:21, 47.8MB/s]
NhaA_non_water.xtc: 10%|█ | 119M/1.15G [00:03<00:21, 47.9MB/s]
NhaA_non_water.xtc: 11%|█ | 127M/1.15G [00:03<00:21, 48.0MB/s]
NhaA_non_water.xtc: 12%|█▏ | 135M/1.15G [00:03<00:20, 48.2MB/s]
NhaA_non_water.xtc: 12%|█▏ | 140M/1.15G [00:04<00:20, 48.7MB/s]
NhaA_non_water.xtc: 13%|█▎ | 146M/1.15G [00:04<00:21, 46.5MB/s]
NhaA_non_water.xtc: 13%|█▎ | 154M/1.15G [00:04<00:21, 47.1MB/s]
NhaA_non_water.xtc: 14%|█▍ | 162M/1.15G [00:04<00:20, 47.3MB/s]
NhaA_non_water.xtc: 15%|█▍ | 170M/1.15G [00:04<00:20, 47.9MB/s]
NhaA_non_water.xtc: 16%|█▌ | 178M/1.15G [00:04<00:20, 47.9MB/s]
NhaA_non_water.xtc: 16%|█▌ | 186M/1.15G [00:05<00:20, 48.0MB/s]
NhaA_non_water.xtc: 17%|█▋ | 194M/1.15G [00:05<00:19, 48.3MB/s]
NhaA_non_water.xtc: 18%|█▊ | 202M/1.15G [00:05<00:19, 48.5MB/s]
NhaA_non_water.xtc: 18%|█▊ | 210M/1.15G [00:05<00:19, 48.4MB/s]
NhaA_non_water.xtc: 19%|█▉ | 217M/1.15G [00:05<00:19, 48.1MB/s]
NhaA_non_water.xtc: 19%|█▉ | 222M/1.15G [00:05<00:18, 49.5MB/s]
NhaA_non_water.xtc: 20%|█▉ | 229M/1.15G [00:05<00:19, 46.3MB/s]
NhaA_non_water.xtc: 21%|██ | 237M/1.15G [00:06<00:19, 46.9MB/s]
NhaA_non_water.xtc: 21%|██▏ | 244M/1.15G [00:06<00:16, 53.3MB/s]
NhaA_non_water.xtc: 22%|██▏ | 248M/1.15G [00:06<00:19, 46.4MB/s]
NhaA_non_water.xtc: 22%|██▏ | 252M/1.15G [00:06<00:20, 44.2MB/s]
NhaA_non_water.xtc: 23%|██▎ | 260M/1.15G [00:06<00:19, 46.0MB/s]
NhaA_non_water.xtc: 23%|██▎ | 268M/1.15G [00:06<00:18, 46.6MB/s]
NhaA_non_water.xtc: 24%|██▍ | 276M/1.15G [00:06<00:18, 47.4MB/s]
NhaA_non_water.xtc: 25%|██▍ | 284M/1.15G [00:07<00:18, 47.2MB/s]
NhaA_non_water.xtc: 25%|██▌ | 292M/1.15G [00:07<00:17, 47.7MB/s]
NhaA_non_water.xtc: 26%|██▌ | 300M/1.15G [00:07<00:17, 47.8MB/s]
NhaA_non_water.xtc: 27%|██▋ | 307M/1.15G [00:07<00:16, 51.8MB/s]
NhaA_non_water.xtc: 27%|██▋ | 311M/1.15G [00:07<00:17, 47.9MB/s]
NhaA_non_water.xtc: 27%|██▋ | 315M/1.15G [00:07<00:18, 45.6MB/s]
NhaA_non_water.xtc: 28%|██▊ | 323M/1.15G [00:07<00:17, 46.6MB/s]
NhaA_non_water.xtc: 29%|██▉ | 331M/1.15G [00:08<00:17, 47.3MB/s]
NhaA_non_water.xtc: 30%|██▉ | 339M/1.15G [00:08<00:16, 47.6MB/s]
NhaA_non_water.xtc: 30%|███ | 346M/1.15G [00:08<00:16, 47.7MB/s]
NhaA_non_water.xtc: 31%|███ | 354M/1.15G [00:08<00:16, 47.9MB/s]
NhaA_non_water.xtc: 32%|███▏ | 362M/1.15G [00:08<00:16, 48.0MB/s]
NhaA_non_water.xtc: 32%|███▏ | 370M/1.15G [00:08<00:16, 48.1MB/s]
NhaA_non_water.xtc: 33%|███▎ | 378M/1.15G [00:09<00:15, 48.5MB/s]
NhaA_non_water.xtc: 34%|███▎ | 386M/1.15G [00:09<00:15, 48.6MB/s]
NhaA_non_water.xtc: 34%|███▍ | 394M/1.15G [00:09<00:15, 48.9MB/s]
NhaA_non_water.xtc: 35%|███▌ | 402M/1.15G [00:09<00:15, 48.8MB/s]
NhaA_non_water.xtc: 36%|███▌ | 410M/1.15G [00:09<00:15, 48.9MB/s]
NhaA_non_water.xtc: 37%|███▋ | 418M/1.15G [00:09<00:14, 49.0MB/s]
NhaA_non_water.xtc: 37%|███▋ | 426M/1.15G [00:10<00:14, 48.9MB/s]
NhaA_non_water.xtc: 38%|███▊ | 434M/1.15G [00:10<00:14, 48.8MB/s]
NhaA_non_water.xtc: 39%|███▊ | 442M/1.15G [00:10<00:14, 48.9MB/s]
NhaA_non_water.xtc: 39%|███▉ | 450M/1.15G [00:10<00:14, 49.2MB/s]
NhaA_non_water.xtc: 40%|███▉ | 458M/1.15G [00:10<00:14, 48.9MB/s]
NhaA_non_water.xtc: 41%|████ | 466M/1.15G [00:10<00:13, 48.7MB/s]
NhaA_non_water.xtc: 41%|████▏ | 474M/1.15G [00:11<00:13, 48.7MB/s]
NhaA_non_water.xtc: 42%|████▏ | 482M/1.15G [00:11<00:13, 48.6MB/s]
NhaA_non_water.xtc: 43%|████▎ | 489M/1.15G [00:11<00:12, 53.3MB/s]
NhaA_non_water.xtc: 43%|████▎ | 494M/1.15G [00:11<00:13, 47.0MB/s]
NhaA_non_water.xtc: 43%|████▎ | 498M/1.15G [00:11<00:14, 45.2MB/s]
NhaA_non_water.xtc: 44%|████▍ | 506M/1.15G [00:11<00:13, 46.4MB/s]
NhaA_non_water.xtc: 45%|████▍ | 514M/1.15G [00:11<00:13, 46.9MB/s]
NhaA_non_water.xtc: 46%|████▌ | 522M/1.15G [00:12<00:13, 47.4MB/s]
NhaA_non_water.xtc: 46%|████▌ | 530M/1.15G [00:12<00:12, 47.6MB/s]
NhaA_non_water.xtc: 47%|████▋ | 537M/1.15G [00:12<00:12, 47.7MB/s]
NhaA_non_water.xtc: 48%|████▊ | 545M/1.15G [00:12<00:12, 48.1MB/s]
NhaA_non_water.xtc: 48%|████▊ | 553M/1.15G [00:12<00:12, 48.2MB/s]
NhaA_non_water.xtc: 49%|████▉ | 561M/1.15G [00:12<00:12, 48.3MB/s]
NhaA_non_water.xtc: 50%|████▉ | 569M/1.15G [00:12<00:11, 48.6MB/s]
NhaA_non_water.xtc: 50%|█████ | 577M/1.15G [00:13<00:11, 48.2MB/s]
NhaA_non_water.xtc: 51%|█████ | 585M/1.15G [00:13<00:10, 52.3MB/s]
NhaA_non_water.xtc: 51%|█████▏ | 589M/1.15G [00:13<00:11, 49.5MB/s]
NhaA_non_water.xtc: 52%|█████▏ | 593M/1.15G [00:13<00:12, 45.8MB/s]
NhaA_non_water.xtc: 52%|█████▏ | 601M/1.15G [00:13<00:11, 46.5MB/s]
NhaA_non_water.xtc: 53%|█████▎ | 608M/1.15G [00:13<00:11, 48.1MB/s]
NhaA_non_water.xtc: 53%|█████▎ | 613M/1.15G [00:13<00:11, 47.6MB/s]
NhaA_non_water.xtc: 54%|█████▍ | 620M/1.15G [00:14<00:11, 45.6MB/s]
NhaA_non_water.xtc: 55%|█████▍ | 627M/1.15G [00:14<00:11, 46.5MB/s]
NhaA_non_water.xtc: 55%|█████▌ | 635M/1.15G [00:14<00:10, 47.1MB/s]
NhaA_non_water.xtc: 56%|█████▌ | 643M/1.15G [00:14<00:10, 47.5MB/s]
NhaA_non_water.xtc: 57%|█████▋ | 651M/1.15G [00:14<00:10, 47.8MB/s]
NhaA_non_water.xtc: 58%|█████▊ | 659M/1.15G [00:14<00:09, 51.5MB/s]
NhaA_non_water.xtc: 58%|█████▊ | 663M/1.15G [00:14<00:10, 46.6MB/s]
NhaA_non_water.xtc: 59%|█████▊ | 671M/1.15G [00:15<00:09, 49.8MB/s]
NhaA_non_water.xtc: 59%|█████▉ | 675M/1.15G [00:15<00:09, 47.4MB/s]
NhaA_non_water.xtc: 59%|█████▉ | 680M/1.15G [00:15<00:09, 48.9MB/s]
NhaA_non_water.xtc: 60%|█████▉ | 686M/1.15G [00:15<00:08, 51.5MB/s]
NhaA_non_water.xtc: 60%|██████ | 691M/1.15G [00:15<00:10, 45.3MB/s]
NhaA_non_water.xtc: 61%|██████ | 698M/1.15G [00:15<00:09, 45.2MB/s]
NhaA_non_water.xtc: 62%|██████▏ | 706M/1.15G [00:15<00:09, 46.2MB/s]
NhaA_non_water.xtc: 62%|██████▏ | 713M/1.15G [00:15<00:08, 51.2MB/s]
NhaA_non_water.xtc: 63%|██████▎ | 718M/1.15G [00:16<00:09, 46.1MB/s]
NhaA_non_water.xtc: 63%|██████▎ | 725M/1.15G [00:16<00:09, 46.1MB/s]
NhaA_non_water.xtc: 64%|██████▍ | 733M/1.15G [00:16<00:08, 46.7MB/s]
NhaA_non_water.xtc: 65%|██████▍ | 741M/1.15G [00:16<00:08, 47.3MB/s]
NhaA_non_water.xtc: 65%|██████▌ | 749M/1.15G [00:16<00:08, 47.7MB/s]
NhaA_non_water.xtc: 66%|██████▌ | 757M/1.15G [00:16<00:08, 47.4MB/s]
NhaA_non_water.xtc: 67%|██████▋ | 765M/1.15G [00:17<00:07, 47.7MB/s]
NhaA_non_water.xtc: 67%|██████▋ | 773M/1.15G [00:17<00:07, 47.9MB/s]
NhaA_non_water.xtc: 68%|██████▊ | 781M/1.15G [00:17<00:07, 48.2MB/s]
NhaA_non_water.xtc: 69%|██████▊ | 787M/1.15G [00:17<00:07, 50.4MB/s]
NhaA_non_water.xtc: 69%|██████▉ | 791M/1.15G [00:17<00:07, 45.6MB/s]
NhaA_non_water.xtc: 70%|██████▉ | 799M/1.15G [00:17<00:07, 49.1MB/s]
NhaA_non_water.xtc: 70%|███████ | 803M/1.15G [00:17<00:07, 45.7MB/s]
NhaA_non_water.xtc: 71%|███████ | 810M/1.15G [00:18<00:07, 45.2MB/s]
NhaA_non_water.xtc: 71%|███████▏ | 818M/1.15G [00:18<00:07, 46.1MB/s]
NhaA_non_water.xtc: 72%|███████▏ | 826M/1.15G [00:18<00:06, 46.7MB/s]
NhaA_non_water.xtc: 73%|███████▎ | 834M/1.15G [00:18<00:06, 47.3MB/s]
NhaA_non_water.xtc: 73%|███████▎ | 842M/1.15G [00:18<00:06, 47.7MB/s]
NhaA_non_water.xtc: 74%|███████▍ | 850M/1.15G [00:18<00:06, 47.8MB/s]
NhaA_non_water.xtc: 75%|███████▍ | 858M/1.15G [00:19<00:05, 48.1MB/s]
NhaA_non_water.xtc: 76%|███████▌ | 865M/1.15G [00:19<00:05, 48.1MB/s]
NhaA_non_water.xtc: 76%|███████▌ | 872M/1.15G [00:19<00:05, 50.4MB/s]
NhaA_non_water.xtc: 77%|███████▋ | 877M/1.15G [00:19<00:05, 47.3MB/s]
NhaA_non_water.xtc: 77%|███████▋ | 884M/1.15G [00:19<00:05, 48.5MB/s]
NhaA_non_water.xtc: 77%|███████▋ | 888M/1.15G [00:19<00:05, 45.5MB/s]
NhaA_non_water.xtc: 78%|███████▊ | 896M/1.15G [00:19<00:05, 46.4MB/s]
NhaA_non_water.xtc: 79%|███████▉ | 903M/1.15G [00:19<00:05, 47.7MB/s]
NhaA_non_water.xtc: 79%|███████▉ | 908M/1.15G [00:20<00:05, 47.1MB/s]
NhaA_non_water.xtc: 80%|███████▉ | 915M/1.15G [00:20<00:05, 45.9MB/s]
NhaA_non_water.xtc: 80%|████████ | 922M/1.15G [00:20<00:04, 51.7MB/s]
NhaA_non_water.xtc: 81%|████████ | 926M/1.15G [00:20<00:04, 47.7MB/s]
NhaA_non_water.xtc: 81%|████████ | 930M/1.15G [00:20<00:04, 44.9MB/s]
NhaA_non_water.xtc: 82%|████████▏ | 938M/1.15G [00:20<00:04, 46.2MB/s]
NhaA_non_water.xtc: 83%|████████▎ | 946M/1.15G [00:20<00:04, 47.0MB/s]
NhaA_non_water.xtc: 83%|████████▎ | 953M/1.15G [00:21<00:03, 50.2MB/s]
NhaA_non_water.xtc: 84%|████████▎ | 957M/1.15G [00:21<00:04, 45.0MB/s]
NhaA_non_water.xtc: 84%|████████▍ | 964M/1.15G [00:21<00:03, 52.0MB/s]
NhaA_non_water.xtc: 84%|████████▍ | 968M/1.15G [00:21<00:04, 43.8MB/s]
NhaA_non_water.xtc: 85%|████████▍ | 973M/1.15G [00:21<00:03, 45.2MB/s]
NhaA_non_water.xtc: 85%|████████▌ | 979M/1.15G [00:21<00:03, 43.9MB/s]
NhaA_non_water.xtc: 86%|████████▌ | 987M/1.15G [00:21<00:03, 45.2MB/s]
NhaA_non_water.xtc: 87%|████████▋ | 995M/1.15G [00:21<00:03, 46.2MB/s]
NhaA_non_water.xtc: 87%|████████▋ | 1.00G/1.15G [00:22<00:02, 50.4MB/s]
NhaA_non_water.xtc: 88%|████████▊ | 1.01G/1.15G [00:22<00:03, 44.8MB/s]
NhaA_non_water.xtc: 89%|████████▊ | 1.01G/1.15G [00:22<00:02, 46.3MB/s]
NhaA_non_water.xtc: 89%|████████▉ | 1.02G/1.15G [00:22<00:02, 46.8MB/s]
NhaA_non_water.xtc: 90%|████████▉ | 1.03G/1.15G [00:22<00:02, 47.5MB/s]
NhaA_non_water.xtc: 91%|█████████ | 1.04G/1.15G [00:22<00:02, 48.0MB/s]
NhaA_non_water.xtc: 91%|█████████▏| 1.05G/1.15G [00:23<00:02, 48.2MB/s]
NhaA_non_water.xtc: 92%|█████████▏| 1.05G/1.15G [00:23<00:01, 49.2MB/s]
NhaA_non_water.xtc: 92%|█████████▏| 1.06G/1.15G [00:23<00:01, 47.7MB/s]
NhaA_non_water.xtc: 93%|█████████▎| 1.07G/1.15G [00:23<00:01, 46.6MB/s]
NhaA_non_water.xtc: 94%|█████████▎| 1.07G/1.15G [00:23<00:01, 51.2MB/s]
NhaA_non_water.xtc: 94%|█████████▍| 1.08G/1.15G [00:23<00:01, 46.4MB/s]
NhaA_non_water.xtc: 95%|█████████▍| 1.08G/1.15G [00:23<00:01, 46.4MB/s]
NhaA_non_water.xtc: 95%|█████████▍| 1.09G/1.15G [00:23<00:01, 45.8MB/s]
NhaA_non_water.xtc: 96%|█████████▌| 1.10G/1.15G [00:24<00:01, 46.4MB/s]
NhaA_non_water.xtc: 96%|█████████▋| 1.10G/1.15G [00:24<00:00, 47.2MB/s]
NhaA_non_water.xtc: 97%|█████████▋| 1.11G/1.15G [00:24<00:00, 47.7MB/s]
NhaA_non_water.xtc: 98%|█████████▊| 1.12G/1.15G [00:24<00:00, 48.0MB/s]
NhaA_non_water.xtc: 98%|█████████▊| 1.13G/1.15G [00:24<00:00, 49.9MB/s]
NhaA_non_water.xtc: 99%|█████████▉| 1.13G/1.15G [00:24<00:00, 48.0MB/s]
NhaA_non_water.xtc: 99%|█████████▉| 1.14G/1.15G [00:24<00:00, 49.4MB/s]
NhaA_non_water.xtc: 100%|█████████▉| 1.14G/1.15G [00:25<00:00, 48.0MB/s]
NhaA_non_water.xtc: 1.15GB [00:25, 45.7MB/s]
NOTE:
Downloading data from MDAnalysisData may take some time.
Grab a coffee! ☕ 🙃
The resulting Bunch object includes the DESCR attribute with a human-readable description of the dataset:
print(NhaA.DESCR)
.. -*- coding: utf-8 -*-
.. _`nhaa-equilibrium-dataset`:
NhaA equilibrium trajectory dataset
===================================
Molecular dynamics (MD) trajectory of the NhaA membrane protein in a
POPE:POPG 4:1 model membrane. The equilibrium simulation was performed
in the NPT ensemble at T=300K and P=1 bar. The system was simulated
with Gromacs 5.1.4, using the CHARMM36 force field, the TIP3P explicit
water model, and NaCl at approximately 100 mM concentration.
Trajectory frames were saved every 100 ps for a total of 500 ns
simulated time. The topology only contains the protein, membrane and
ions (because the water molecules were stripped from the trajectory to
save space).
The topology is contained in the NhaA_non_water.gro file. The trajectory is contained in the NhaA_non_water.xtc file.
Notes
-----
Data set characteristics:
:size: 1.07 GB
:number of trajectories: 1
:number of frames: 5000
:number of particles: 60702
:creator: Ian Kenney, Shujie Fan
:URL: `10.6084/m9.figshare.7185203.v2 <https://doi.org/10.6084/m9.figshare.7185203.v2>`_
:license: `CC-BY 4.0 <https://creativecommons.org/licenses/by/4.0/legalcode>`_
:reference: [Kenney2018]_
.. [Kenney2018] M. Kenney, Ian; Fan, Shujie; Beckstein, Oliver (2018): Molecular dynamics
trajectory of membrane protein NhaA. figshare. Dataset.
MDAnalysis. figshare. Fileset. doi:
`10.6084/m9.figshare.7185203.v2
<https://doi.org/10.6084/m9.figshare.7185203.v2>`_
Now that we have our NhaA dataset, we can load the corresponding Universe.
2. Load MDA Universe
This tutorial exercise is intended as an example on how to use the tool but does not provide information about the membrane properties around the protein or induced by the protein.
In this particular membrane-protein example, we find that the NhaA antiporter diffuses in the membrane because the simulation was carried out without position restraints on the protein (See Appendix 1). As a consequence, the results obtained from MembraneCurvature should not be used as a means to assess membrane curvature induced by the protein or to study properties of lipids around the protein.
If the interest is to determine membrane curvature effects driven by protein insertion while the protein is allowed to freely diffuse in the membrane, the current version of MembraneCurvature requires further processing of the input trajectory obtained from the MD simulations. For more details see Part I and Part II of the blog posts that discuss how to to calculate membrane curvature in different MD simulation setups.
We can load the data stored in NhaA into an MDAnalysis Universe with:
universe = mda.Universe(NhaA.topology, NhaA.trajectory)
universe
<Universe with 60702 atoms>
A quick check of the Universe gives a better idea of the system used in this tutorial. To check the dimensions of the simulation box, number of residues in the system, and number of frames in the trajectory contained in our Universe, we can run:
print("Universe info:")
print("\nBox dimensions[x,y,z]: {}." \
"".format(universe.dimensions[:3]))
print("\n{} residues and {} atoms." \
"".format(universe.residues.n_residues, universe.residues.n_atoms))
print("\nThe trajectory includes {} frames."\
"".format(universe.trajectory.n_frames))
Universe info:
Box dimensions[x,y,z]: [119.847176 119.847176 99.499016].
1314 residues and 60702 atoms.
The trajectory includes 5001 frames.
To visualize our system, we can make use of the NGL Viewer and the function view_system:
def view_system(system, zoom, surface=False, box=False):
view = nv.show_mdanalysis(system)
view.add_representation('ball+stick', selection='POPE and .P', radius=2, color="salmon")
view.add_representation('ball+stick', selection='POPG and .P', radius=2, color="darkcyan")
view.add_representation('ball+stick', selection='POPG or POPE', radius=0.25, color="grey", opacity=0.05)
view.add_cartoon("protein", color='navy')
if surface:
view.add_surface("protein", color='skyblue', opacity=0.05)
if box:
view.add_unitcell()
view.camera='orthographic'
view.control.zoom(zoom)
view.control.rotate(
mda.lib.transformations.quaternion_from_euler(
-np.pi/2, -np.pi/3, np.pi/12, 'rzyz').tolist())
return view
view_system(universe, 0.35, surface=True, box=False)
In the visualization obtained above, we can identify lipid headgroups from POPE molecules coloured in pink, and lipid headgroups from POPG lipids coloured in cyan. The NhaA antiporter is shown in cartoon and surface representation, blue coloured.
Since our system contains protein and membrane, it is useful to identify the lipids and protein separately. We can make use of the MDAnalysis select_atoms to save the AtomGroup in two different variables:
protein = universe.select_atoms("protein")
print("\nProtein has {} residues." \
"".format(protein.n_residues))
lipids = universe.select_atoms("not protein")
print("\nLipids molecules include {} residues and {} atoms." \
"".format(lipids.n_residues, lipids.n_atoms))
Protein has 752 residues.
Lipids molecules include 562 residues and 49078 atoms.
3. Select Atoms of Reference
The MembraneCurvature algorithm uses an AtomGroup as a reference to derive surfaces. In the NhaA system, we can use the phospholipid headgroups from the POPE and POPG molecules as an AtomGroup of reference. Given that the NhaA system was obtained from an all atom MD simulation, for our AtomGroup we are going to choose the phosphate phosphorus atoms.
Then, we can select the atoms with the name P in our Universe with:
P_headgroups = universe.select_atoms('name P')
The AtomGroup here assigned by the selection name P, however, include the phosphorus atoms of all the lipids in our system instead of the individual leaflets. Since MembraneCurvature calculates curvature from a single surface, we should select an AtomGroup for lower and upper leaflet, separately.
For illustrative purposes, this exercise uses the leaflet module from MDAnalysis to identify the residues belonging to each bilayer leaflet. When using MembraneCurvature, defining accurate atom selections for each leaflet is the responsibility of the user, and the resulting leaflet assignments should always be inspected and validated before performing curvature analysis.
Reference atoms should remain associated with the same leaflet throughout the simulation. Therefore, atoms that undergo flip-flop events are not suitable for curvature analysis. Incorrect or unstable leaflet assignments can introduce artifacts and lead to inaccurate membrane curvature calculations and misleading interpretations of the resulting curvature maps.
To use the MDAnalysis LeafletFinder we do:
from MDAnalysis.analysis.leaflet import LeafletFinder
L_nhaa = LeafletFinder(universe, 'name P', cutoff=20)
nhaa_upper_leaflet = L_nhaa.groups(0) # upper leaflet
nhaa_lower_leaflet = L_nhaa.groups(1) # lower leafet
type(nhaa_upper_leaflet)
MDAnalysis.core.groups.AtomGroup
Here, the variables nhaa_upper_leaflet and nhaa_lower_leaflet are AtomGroups.
Let’s check the residues in each leaflet:
nhaa_lower_leaflet.residues.resids
array([425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437,
438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450,
451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463,
620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632,
633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645,
646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658,
659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 670, 671,
672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684,
685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697,
698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710,
711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723,
724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736,
737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749,
750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762,
763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775])
nhaa_upper_leaflet.residues.resids
array([386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398,
399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411,
412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424,
464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476,
477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489,
490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502,
503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515,
516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528,
529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541,
542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554,
555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567,
568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580,
581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593,
594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606,
607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619])
Then, we can select the phosphorus headgroup atoms from each leaflet:
leaflets = ['Lower', 'Upper']
nhaa_upper_leaflet_P = nhaa_upper_leaflet.select_atoms("name P")
nhaa_lower_leaflet_P = nhaa_lower_leaflet.select_atoms("name P")
for name, new_lf in zip(leaflets,[nhaa_lower_leaflet_P, nhaa_upper_leaflet_P]):
print("{} leaflet includes {} elements.".format(name,
len(new_lf)))
Lower leaflet includes 195 elements.
Upper leaflet includes 195 elements.
Similarly, we can select the selection string from the AtomGroup, and turn them into a human readable string.
sel_upper = " ".join([str(r) for r in nhaa_upper_leaflet.residues.resids])
sel_lower = " ".join([str(r) for r in nhaa_lower_leaflet.residues.resids])
upper_string = "resid {} and name P".format(sel_upper)
lower_string = "resid {} and name P".format(sel_lower)
The selection sel_upper and sel_lower will be useful to visualize the selections later on.
NOTE:
An alternative approach to select residues for each leaflet can be found in Appendix 2.
To double-check we have the right residues in each leaflet, we can visually check the selection of atoms using NGL Viewer. We use the function color_by_leaflet to visualize the AtomGroup for each leaflet according to the selection found above.
def color_by_leaflet(system, sel_up, sel_lw, zoom, box=False):
view = nv.show_mdanalysis(system)
view.add_ball_and_stick(sel_up,
radius=2.00,
color="orange")
view.add_ball_and_stick(sel_lw,
radius=2.00,
color="green")
if box == True:
view.add_unitcell()
view.camera='orthographic'
view.control.zoom(zoom)
view.control.rotate(
mda.lib.transformations.quaternion_from_euler(
-np.pi/2, np.pi/3, np.pi/12, 'rzyz').tolist())
return view
color_by_leaflet(nhaa_upper_leaflet+nhaa_lower_leaflet, sel_upper, sel_lower, 0.05, box=True)
In the widget above we identify the selection of lipid headgroups in the lower (green coloured) and the upper leaflet (orange coloured). Each leaflet is one AtomGroup, which will be used as a reference to derive surfaces and calculate curvature.
NOTE:
Validating that your leaflet selection is correct is a fundamental step to run MembraneCurvature. Before running, double-check the selections do not include mixed residues between leaflets. Keep in mind that incorrect or incomplete selection in the AtomGroup of reference will therefore provide inaccurate calculations of membrane curvature.
Now that we have identified the AtomGroup of reference, we can run MembraneCurvature.
4. Run MembraneCurvature
MembraneCurvature is a Python class that performs multiframe analyses to derive surfaces from the AtomGroup of reference. From the derived surface, MembraneCurvature calculates mean and Gaussian curvature per frame and their respective average over frames.
To derive surfaces and calculate curvature, MembraneCurvature offers two different methods:
4.1 Binning Method
To use MembraneCurvature, we have five parameters to provide:
MembraneCurvature(universe, # universe
select='name P', # selection of reference
surface_method='binning', # use binning as surface method
n_x_bins=10, # number of bins in the x dimension
n_y_bins=10, # number of bins in the y_dimension
wrap=True) # wrap coordinates to keep atoms in the main unit cell
Universe: The Universe that contains our system of interest. In this tutorial, our Universe comprises a protein embedded in POPE POPG lipid bilayer, as described in Section 2.
Atom selection (
select): This is a key parameter to run MembraneCurvature. Based on this selection, surfaces will be derived for every frame in the trajectory. Simultaneously, curvature will be calculated from the derived surface. In this tutorial, we selected the leaflets contained in our Universe as shown in 3 Select Atoms of Reference.
Number of bins (
n_x_bins,n_y_bins): This parameter will determine how many bins are assigned to the grid in each dimension. In Membrane Curvature, the dimensions of the grid are determined by the size of the simulation box contained in the Universe. The grid comprisesn_x_binsxn_y_binsnumber of bins. Choosing the number of bins is also important. Choosing too many bins may introduce undefined values in the grid built by the MembraneCurvature algorithm. On the other hand, a very low number of bins may result in significant loss of information. Since we recommend taking bins of size ~20Å, and the Universe used in this example has dimensions ~120 x 120 Å, we would typically usen_x_bins=n_y_bins=6. However, in this tutorial we are going to usen_x_bins=n_y_bins=10to illustrate key differences bewteen the two surface methods available in Membrane Curvature.Coordinate wrapping (
wrap): Applying coordinate wrapping is useful when we have atoms falling outside the boundaries of the simulation box. Since this is a raw trajectory, and we want to have a high number of lipid headgroups to derive the surface, we usewrap=Trueto put all the atoms in the primary unit cell.
With a clear definition of the AtomGroup associated to each leaflet and the parameters needed to calculate membrane curvature in our system, we are ready to run MembraneCurvature.
Running MembraneCurvature for the NhaA-membrane system is expected to take less than a minute for each leaflet.
NOTE:
For the examples below we are going to use the first 1500 frames. However, when performing your analysis, be aware of the number of frames you provide to the MDAnalysis run method.
We run MembraneCurvature for the upper leaflet with:
curvature_upper_leaflet_binning = MembraneCurvature(universe,
surface_method='binning',
select=upper_string,
n_x_bins=10,
n_y_bins=10,
wrap=True).run(stop=1500) # use the first 1500 frames
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 4 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 1 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 7 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 5 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 6 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 2 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 3 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 9 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 10 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 11 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 8 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/base.py:341: RuntimeWarning: Mean of empty slice
z_average = np.nanmean(self.results.z_surface, axis=0)
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/base.py:350: RuntimeWarning: Mean of empty slice
self.results.average_mean = np.nanmean(self.results.mean, axis=0)
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/base.py:351: RuntimeWarning: Mean of empty slice
self.results.average_gaussian = np.nanmean(self.results.gaussian, axis=0)
As detailed in the Algorithm page, when using surface_method='binning', MembraneCurvature uses the dimensions of the box to set the dimensions of the grid. The cell above shows a warning message indicating that one of the elements in the AtomGroup of reference has negative coordinates, hence falling outside the grid boundaries established by the dimensions of the NhaA-membrane simulation box (~120 Å x 120 Å divided into 100 bins):
Atom coordinates exceed size of grid and element (10,3) can't be assigned. Maximum (x,y) coordinates must be < (119.84717559814453, 119.84717559814453). Skipping atom.
warnings.warn(self.msg.format(*args))
The warning tells the user that the number of atoms used in the upper leaflet is reduced. Instead of deriving a surface from the 195 elements included in the nhaa_upper_leaflet_P selection, the surface uses 194 elements in the AtomGroup. Deriving a surface with sufficient sampling will reduce or eliminate this warning.
Similarly, for the lower leaflet, we run:
curvature_lower_leaflet_binning = MembraneCurvature(universe,
select=lower_string,
surface_method='binning',
n_x_bins=10,
n_y_bins=10,
wrap=True).run(stop=1500)
/home/docs/checkouts/readthedocs.org/user_builds/membrane-curvature/conda/latest/lib/python3.14/site-packages/membrane_curvature/binning_surface.py:108: UserWarning: 12 atoms fall outside the grid boundaries. Skipping atoms.
warnings.warn(f'{atoms_outside_grid} atoms fall outside the grid boundaries. Skipping atoms.')
4.2 Fourier Method
To use the Fourier method, set surface_method='fourier' in MembraneCurvature:
MembraneCurvature(universe, # universe
select='name P', # selection of reference
surface_method='fourier', # run with Fourier
fourier_m=2, # max Fourier mode index in x (default 2)
fourier_n=2, # max Fourier mode index in y (default 2)
)
Unlike binning, the Fourier path does not build the surface by averaging atoms into bins. It fits a truncated periodic series directly to atom heights. The parameters that control the Fourier fit are fourier_m and fourier_n (both default to 2), thus setting n_x_bins or n_y_bins is not required.
These two Fourier-specific parameters are:
fourier_m(default:2): maximum mode index in \(x\). This represents the number of wavelengths along the box length \(L_x\) that may appear in the fit.fourier_n(default:2): maximum mode index magnitude in \(y\). Similarly, this represents the number of wavelengths along the box length \(L_y\) that may appear in the fit.
See the Fourier surface API for the full workflow and mode-count details.
For most leaflet selections, the defaults fourier_m = fourier_n = 2 are sufficient!
Do not use higher mode values without a clear reason. Values above the default can fit noise in the atom height data and make
curvature less stable. Increase fourier_m or fourier_n only if you need shorter wavelengths and your selection has many
more atoms than fitted parameters. Beware when curvature stops improving systematically and starts becoming noisier.
The parameters n_x_bins, n_y_bins are also accepted in the Fourier method, but they do not define binning!
When choosing the Fourier method, the parameters n_x_bins, n_y_bins although accepted, they are not used to group or average atoms. Instead, they are used to set the resolution of the output grid. This is what the determines the shape of the array with the surface and curvature results. Default values are set to n_x_bins=n_y_bins=100. We recommend using these default values.
For the purpose of this tutorial, we are going to use Fourier with n_x_bins, n_y_bins to intentionally match the number of bins in surface_method='binning'.
In general, when using Fourier, we recommend using the default values n_x_bins and n_y_bins.
To get the Fourier surface for the upper leaflet using the default values for fourier_m and fourier_n, and n_x_bins, n_y_bins to intentionally match the same number of bins as in the example shown in Section 4.1 - Binning method, we do:
curvature_upper_leaflet_fourier = MembraneCurvature(universe,
select=upper_string,
surface_method='fourier',
n_x_bins=10,
n_y_bins=10,
).run(stop=1500)
And for the lower leaflet:
curvature_lower_leaflet_fourier = MembraneCurvature(universe,
select=lower_string,
n_x_bins=10,
n_y_bins=10,
surface_method='fourier',
).run(stop=1500)
5. Extract Results
Regardless of the method used to derive surfaces, MembraneCurvature stores the results in the .results attribute.
For surface, mean curvature, and Gaussian curvature, we can find two different types of results:
Multiframe
MembraneCurvature derives the surface from the AtomGroup of reference, as well as mean and Gaussian curvature, for every frame in the trajectory.
curvature_upper_leaflet.results.z_surface
curvature_upper_leaflet.results.mean_curvature
curvature_upper_leaflet.results.gaussian_curvature
Average over frames
The derived surface, and calculated mean and Gaussian curvature are averaged over the n_frames of the trajectory.
curvature_upper_leaflet.results.average_z_surface
curvature_upper_leaflet.results.average_mean
curvature_upper_leaflet.results.average_gaussian
In this tutorial, we are going to extract the average over frames results for surface, mean and Gaussian curvature. You can find an example of MembraneCurvature using the Multiframe results in the Membrane-only system tutorial.
5.1 Surface
In this tutorial, both methods derive the surfaces using the same reference AtomGroup, which corresponds to the phosphorus atoms of the lipids in each leaflet
We can access the results for surface with:
# Surface results with the binning method
nhaa_surface_upper_leaflet_binning = curvature_upper_leaflet_binning.results.average_z_surface
nhaa_surface_lower_leaflet_binning = curvature_lower_leaflet_binning.results.average_z_surface
# Surface results with the Fourier method
nhaa_surface_upper_leaflet_fourier = curvature_upper_leaflet_fourier.results.average_z_surface
nhaa_surface_lower_leaflet_fourier = curvature_lower_leaflet_fourier.results.average_z_surface
Notes on the binning method
For the binning method, the surface of the upper and lower leaflet is derived from the z positions of the elements in the AtomGroup of reference. For more details see get_z_surface() in the API documentation.
It’s important to note that nhaa_surface_upper_leaflet_binning is an array of shape (n_x_bins, n_y_bins). The values in the grid contain the averaged \(Z\) coordinate of the lipid headgroups in each bin.
We can check the shape of the resulting array from the binning method with:
nhaa_surface_upper_leaflet_binning.shape
(10, 10)
A shape of (10, 10) means that we have a total of 100 bins in the grid. The value stored in each element of nhaa_surface_upper_leaflet_binning is the averaged \(Z\) coordinate of the lipid headgroups assigned to that bin.
Since we are retriving the average \(Z\) surface from the upper leaflet, the shape of the array indicates that the grid used to derive the surface was divided into 10 bins in both \(x\) and \(y\) dimensions.
Notes on the Fourier method
For the Fourier method, the surfaces associated with the upper and lower leaflet are derived from the \((x, y, z)\) positions of the AtomGroup of reference by fitting a truncated 2D Fourier series to atom heights (least squares), then evaluating that periodic surface on a grid of bin centres. For more details see fourier_height_from_atoms() in the API documentation.
We can check the shape of the resulting array from the Fourier method with:
nhaa_surface_upper_leaflet_fourier.shape
(10, 10)
A shape of (10, 10) means the continuous surface fit was evaluated on a \(10 \times 10\) grid of bin centres across the simulation box.
Do not interpret those 100 values as 100 independent bin averages. Using n_x_bins and n_y_bins in Fourier does not define spatial binning. Every point is obtained from the same global Fourier fit to the leaflet atoms.
If we had ommited the parameters n_x_bins and n_y_bins in the Fourier method, we would have found that the array shape is (100, 100). Because even when not directly passed by the user, MembraneCurvature in the Fourier mode defaults to n_x_bins=100 and n_y_bins=100.
Hence, a shape of (10, 10) is expected since we set n_x_bins=10, n_y_bins=10 explicitly to match the binning example.
Note that although with surface_method='fourier' the obtained surfaces can be smoother and less sensitive to empty bins than using the surface_method='binning', the level of detail is controlled mainly by the parameter fourier_m and fourier_n, not by how many output bins you choose.
5.2 Curvature
The calculated values of curvature averaged over frames are stored in the attributes .results.average_mean_curvature and
.results.average_gaussian_curvature. These arrays contain the computed values of mean (\(H\)) and Gaussian (\(K\)) curvature averaged over the n_frames of the trajectory. Note that the units of mean curvature is 1/lenght (\(Å^{-1}\)), and the units of Gaussian curvature is 1/[length^2] (\(Å^{-2}\)).
For details of the calculation of mean and Gaussian curvature see the Curvature API documentation.
Mean Curvature (\(H\))
We can save the calculated values (np.ndarrays) of membrane curvature for each leaflet in two new variables called nhaa_mean_upper_leaflet and nhaa_mean_lower_leaflet:
nhaa_mean_upper_leaflet_binning = curvature_upper_leaflet_binning.results.average_mean
nhaa_mean_lower_leaflet_binning = curvature_lower_leaflet_binning.results.average_mean
nhaa_mean_upper_leaflet_fourier = curvature_upper_leaflet_fourier.results.average_mean
nhaa_mean_lower_leaflet_fourier = curvature_lower_leaflet_fourier.results.average_mean
Since we are extracting the results of averaged mean curvature, the shape of these arrays is simply (n_x_bins, n_y_bins):
nhaa_mean_upper_leaflet_binning.shape
(10, 10)
nhaa_mean_upper_leaflet_fourier.shape
(10, 10)
In nhaa_mean_upper_leaflet, each element of the array correspond to the value of mean curvature in the respective bin.
Gaussian Curvature (\(K\))
Since the information of the averaged Gaussian curvature is stored in the results.average_gaussian attribute, we can store the values of \(K\) in another two variables, nhaa_gaussian_*er_leaflet_binning and nhaa_gaussian_*er_leaflet_fourier:
nhaa_gaussian_upper_leaflet_binning = curvature_upper_leaflet_binning.results.average_gaussian
nhaa_gaussian_lower_leaflet_binning = curvature_lower_leaflet_binning.results.average_gaussian
nhaa_gaussian_upper_leaflet_fourier = curvature_upper_leaflet_fourier.results.average_gaussian
nhaa_gaussian_lower_leaflet_fourier = curvature_lower_leaflet_fourier.results.average_gaussian
Now that we have extracted the data containing values of surface, mean and Gaussian curvature, we can proceed to visualize the results.
In Summary:
Surface and curvature results are both 2D arrays of shape (n_x_bins, n_y_bins) and are accessed from .results (.average_z_surface for the height field, .average_mean and .average_gaussian for curvature), the methods differ in how those values are obtained:
Binning: atoms are assigned to bins in a grid. The full grid has shape (
n_x_bins,n_y_bins). Each cell stores the mean \(Z\) of atoms in that bin (empty bins map tonumpy.nan). Mean and Gaussian curvature are computed from numerical derivatives of that grid vianumpy.gradienton the binned surface.Fourier: a truncated periodic Fourier series is fitted by least squares to all atom heights (
fourier_m,fourier_n). The fit uses every atom once.n_x_binsandn_y_binsonly set how many points that surface is sampled on. Each grid point stores the fitted height at a bin centre, not a local average. Curvature is the result of analytic derivatives of the Fourier series.
6. Visualize Results
There are two main approaches to visualize results from MembraneCurvature.
Plotting MembraneCurvature results using contours or imshow requires setting the correct array orientation.
To plot results from the membrane curvature analysis, use the transposed array np.array([...]).T, and place the [0, 0] index of the array in the lower left corner of the axes by setting origin='lower'.
Please check the examples below for more details.
Contours plots
To visualize the results obtained from the MembraneCurvature we can use contourf from Matplotlib. By using contourf to plot MembraneCurvature results, we gather all the points of equal value and the region enclosed by that set of points is coloured according to a colormap of choice.
Binning method vs. Fourier method
We will plot the averaged surface obtained for the lower leaflet using the two methods (binning and Fourier) compared side-to-side:
methods = ["Binning", "Fourier"]
surfaces = [nhaa_surface_lower_leaflet_binning, nhaa_surface_lower_leaflet_fourier]
fig, axes = plt.subplots(1, 2, figsize=(4, 3.5), dpi=200)
plt.suptitle("Lower leaflet", y=0.8, fontsize=9)
for ax, surface, method in zip(axes, surfaces, methods):
ax.set_title(f'{method} method', fontsize=8)
ax.contourf(surface.T, cmap='YlGnBu_r', origin='lower', levels=20)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
Note that since we use phospholipids as atoms of reference, and the system we are using to derive surfaces from is a membrane-protein system, it is expected to find some np.nans when using the binning method. Empty bins (zero samples) are represented as numpy.nan in the returned (n_x_bins, n_y_bins) array and when plotted they appear as empty.
Similarly for the upper leaflet, and for each method:
surfaces = [nhaa_surface_upper_leaflet_binning, nhaa_surface_upper_leaflet_fourier]
fig, axes = plt.subplots(1, 2, figsize=(4, 3.5), dpi=200)
plt.suptitle("Upper leaflet", y=0.8, fontsize=9)
for ax, surface, method in zip(axes, surfaces, methods):
ax.set_title(f'{method} method', fontsize=8)
ax.contourf(surface.T, cmap='YlGnBu_r', origin='lower', levels=20)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
It’s important to note that one of the main advantages of the Fourier approach is that it does not average atoms per cell, so there are no “empty bins.” Each array value comes from the same fitted surface sampled on the full grid, which is why you typically see no np.nan in the surface or curvature maps.
In plot above, the binning method shows an empty region and regions with rough edges. It also lacks a colorbar for reference. A more elaborate way to plot contours uses the zoom function from the scipy.ndimage package for multidimensional processing. The zoom interpolation is used to smooth the contours.
However, a major limitation when working with scipy.ndimage is that most functions do not treat np.nan values. Instead, NaNs propagate through the computations and can therefore distort the output. As mentioned earlier, in membrane–protein systems np.nan values are expected when using surface_method="binning" because the region occupied by the protein does not contribute data to the grid. As a consequence, the resulting surface array contains undefined grid points that must be interpolated before applying scipy.ndimage operations.
For this exercise, we will interpolate the surface and curvature arrays to demonstrate why the Fourier-based approach is recommended for membrane-protein systems.
We strongly discourage interpolating surfaces that contain undefined regions. Interpolation of partially sampled surfaces can introduce noise and artificial features into curvature maps, particularly near sparse, masked, or reconstructed regions. In membrane-protein systems, NaN values in the vicinity of the protein are unavoidable, and any extrapolation across these regions may substantially distort the calculated curvature and lead to misleading interpretations of the resulting plots generated with MembraneCurvature.
An example to interpolate our surface array is in the function interpolate_nans:
from scipy.interpolate import griddata
def interpolate_nans(z):
"""
Fill NaN values in a 2D array using linear interpolation,
falling back to nearest-neighbor interpolation at boundaries.
"""
z = np.asarray(z, dtype=float)
if not np.isnan(z).any():
return z
valid_mask = ~np.isnan(z)
if not valid_mask.any():
return np.zeros_like(z)
y, x = np.indices(z.shape)
points = (x[valid_mask], y[valid_mask])
values = z[valid_mask]
filled = z.copy()
missing_mask = ~valid_mask
filled[missing_mask] = griddata(points,
values,
(x[missing_mask], y[missing_mask]),
method="linear")
remaining_mask = np.isnan(filled)
if remaining_mask.any():
filled[remaining_mask] = griddata(points,
values,
(x[remaining_mask], y[remaining_mask]),
method="nearest")
return filled
Now that we have a function to interpolate the array for the results obtained with the binning method, we can proceed to plot with scipy.ndimage using the plot_contours function:
def plot_contours(title, results, label, levels_, cm):
"""
Function used to plot contours of MembraneCurvature results.
User can determine number of contour lines / regions (levels),
label of the colorbar (label) and colormap (cmap).
Parameters
----------
results: list
List with results by leaflets as elements [lower_leaflet, upper_leaflet]
label: str
Label to add to colorbar.
levels_: int
Determines number of contour lines.
cm: str
Colormap to use in plot.
"""
fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(4, 3.5), dpi=200)
finite = [rs[np.isfinite(rs)] for rs in results]
max_ = max((np.max(np.abs(v)) for v in finite if v.size), default=0.0)
plt.suptitle(title, fontsize=8, y=0.8)
for ax, rs, lf in zip((ax1, ax2), results, leaflets):
rs = interpolate_nans(np.asarray(rs).T)
rs = ndimage.zoom(rs, 3, mode='wrap', order=1)
if np.min(rs) < 0 < np.max(rs):
levs = np.linspace(-max_, max_, levels_)
im = ax.contourf(
rs, cmap=cm, origin='lower', levels=levs, alpha=0.95, vmin=-max_, vmax=max_
)
tcs = [-max_, 0, max_]
else:
vmin, vmax = np.min(rs), np.max(rs)
levs = np.linspace(vmin, vmax, levels_)
im = ax.contourf(
rs, cmap=cm, origin='lower', levels=levs, alpha=0.95, vmin=vmin, vmax=vmax
)
tcs = [vmin, vmax]
ax.set_aspect('equal')
ax.set_title('{} Leaflet'.format(lf), fontsize=6)
ax.axis('off')
cbar = plt.colorbar(
im, ticks=tcs, orientation='horizontal', ax=ax, shrink=0.7, aspect=10, pad=0.05
)
cbar.ax.tick_params(labelsize=4, width=0.5)
cbar.set_label(label, fontsize=6, labelpad=2)
plt.tight_layout()
Note that there are five arguments in the plot_contours function. The first is the method -here set at the title of the plot-, the second is a list that contains the results of curvature by leaflet, the third one is the label to attach to the colorbar, the fourth one determines the number of contours in the plot, and the fifth one is the colormap to use in the plot.
To compare the surfaces derived from each method, we can use plot_contours as in:
plot_contours("Binning method", [nhaa_surface_lower_leaflet_binning, nhaa_surface_upper_leaflet_binning], 'Surface', 35, 'YlGnBu_r')
plot_contours("Fourier method", [nhaa_surface_lower_leaflet_fourier, nhaa_surface_upper_leaflet_fourier], 'Surface', 35, 'YlGnBu_r')
Comparing the surface contour plots produced by the two methods highlights several important differences:
Binning produces a more patchy, cell-like surface because local structure depends directly on the atoms contained within each spatial bin. Empty bins may contain
np.nanvalues, which for visualization, are interpolated before contouring. As a result, portions of the surface may represent inferred heights rather than direct averages from sampled atoms.Fourier produces a smooth, continuous surface across the entire grid using a single global fit to all leaflet atoms. Because the Fourier representation is defined everywhere on the grid, there are generally no empty regions and no need to interpolate missing values for plotting.
Consequently, the differences between the two figures arise not only from the reconstruction methods themselves, but also from how unsampled regions are treated. In the binning approach, interpolation can significantly affect the resulting surface shape, especially in sparsely sampled areas. In contrast, the Fourier method enforces smoothness intrinsically through the truncated Fourier series (fourier_m, fourier_n).
These differences become especially important when calculating membrane curvature.
Curvature estimates rely on spatial derivatives computed using finite-difference operations (via numpy.gradient), and missing or poorly sampled regions can strongly affect the result. In binned surfaces, np.nan values and interpolation artifacts may be amplified when derivatives are taken, since curvature depends on first- and second-order gradients of the surface.
Particularly for membrane-protein systems, the Fourier approach is the recommended surface method. It is generally more stable because the fitted surface is smooth and continuously defined over the full grid before differentiation.
With these considerations in mind, we plot mean curvature using the plot_contours function:
mean_curvs_binning = [nhaa_mean_lower_leaflet_binning, nhaa_mean_upper_leaflet_binning]
mean_curvs_fourier = [nhaa_mean_lower_leaflet_fourier, nhaa_mean_upper_leaflet_fourier]
plot_contours("Binning method", mean_curvs_binning, "$H$ (Å$^{-1}$)", 30, "bwr")
plot_contours("Fourier method", mean_curvs_fourier, "$H$ (Å$^{-1}$)", 30, "bwr")
Taken together, these mean-curvature maps demonstrate that the choice of surface reconstruction method has a major impact on the resulting curvature analysis. In the binning approach, empty or sparsely sampled cells (np.nan) must be treated before visualization and differentiation. Interpolation can smooth over these gaps, introducing features that were not present in the original bin averages.
Because mean curvature depends on first- and second-order spatial derivatives of the surface, computed in MembraneCurvature using numpy.gradient, these interpolation artifacts, sharp bin boundaries, and sampling discontinuities can become amplified in the final curvature map.
This sensitivity helps explain the substantial differences in both spatial pattern and curvature magnitude between the binning and Fourier results.
By contrast, the Fourier approach constructs a continuous global surface without empty bins and evaluates curvature from analytic derivatives of the fitted model itself. The resulting maps are therefore smoother and less sensitive to uneven spatial sampling or missing regions.
For membrane–protein systems in particular, where leaflet atoms are often distributed non-uniformly across the membrane plane and binned surfaces may contain large unsampled regions—we recommend surface_method='fourier' as the more robust default for interpreting local membrane curvature around proteins and individual leaflets.
For the rest of this tutorial, we are going to continue the analysis with the results obtained from surface_method='fourier'.
NOTE:
As mentioned before, in the MD system considered in this tutorial we observe protein diffusion. Therefore, the results of mean and Gaussian here obtained provide information about the membrane while averaging over protein diffusion instead of the properties of the lipids around the protein or membrane curvature driven by protein insertion. To evaluate properties of the local environment around the protein or membrane curvature as an effect of protein insertion, further trajectory processing is required. (See Note in Section 2)
Mean curvature (\(H\)) describes how the fitted leaflet surface bends in the membrane plane. In MembraneCurvature’s Monge convention, \(H>0\) marks bowl-like regions (local minima in \(z\), depressions), \(H<0\) dome-like regions (local maxima in \(z\), protrusions), and \(H\approx0\) shows nearly flat regions. On the colour map selected for our figure, red is \(H>0\), and blue is \(H<0\).
In other words, mean curvature (\(H\)) gives information about the “inverted shape” of the surface. Positive mean curvature indicates valleys or convexe regions in the membrane, while negative mean curvature indicates peaks or concave regions. \(H=0\) means flat curvature.
Hence, from the obtained contour plot of mean curvature we can identify:
In the lower leaflet, from the bottom left to the upper right corner of the simulation box, a region of positive curvature (red coloured).
Overall, positive curvature dominates in the upper leaflet. However, there is a clear domain of negative curvature in the central region of the leaflet.
By comparing the two plots, the curvature induced by NhaA in the lower leaflet is stronger.
Similar to the approach used for mean curvature, to plot the results of averaged Gaussian curvature we first create the list of curvatures:
gaussian_curvs_fourier = [nhaa_gaussian_lower_leaflet_fourier, nhaa_gaussian_upper_leaflet_fourier]
And then call the plot_contours function:
plot_contours("Fourier", gaussian_curvs_fourier, "$K$ (Å$^{-2}$)", 20, "PiYG")
To read plots of Gaussian curvature, a rule of thumb is: negative Gaussian curvature represents saddle points, positive Gaussian curvature represents locally elliptic bending, points of concave regions. Generally speaking, Gaussian curvature is associated to the elasticity of the membrane.
From the obtained contour plot of Gaussian curvature, we can identify:
In the lower leaflet, most of the membrane surface is relatively flat (\(K\approx 0\), shown in pale/white colors), with only a few localized curved regions. A green area in the lower-left indicates positive Gaussian curvature (dome- or bowl-like bending), while a pink region near the upper-right center indicates negative Gaussian curvature (saddle-like bending).
In the upper leaflet, the curvature pattern is more pronounced and spatially varied. A large dark-green region in the upper-right center corresponds to strong positive curvature, while several pink regions—near the top center, along the right edge, and across part of the lower-left area—indicate negative curvature.
The upper leaflet displays stronger curvature and more localized structural features overall, whereas the lower leaflet remains comparatively flatter, with smaller positive and negative deviations from zero curvature.
For a nice graphical review on curvature, including mean and Gaussian curvature, check this link.
The results obtained of surface, mean curvature, and Gaussian curvature can be also be visualized per leaflet using the function plots_by_leaflet defined below.
def plots_by_leaflet(results):
"""
Generate figure with of surface, $H$ and $K$
as subplots.
"""
cms=["YlGnBu_r", "bwr", "PiYG"]
units=['$Z$ $(\AA)$','$H$ (Å$^{-1})$', '$K$ (Å$^{-2})$']
titles = ['Surface', 'Mean Curvature', 'Gaussian Curvature']
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(7,4), dpi=200)
for ax, mc, title, cm, unit in zip((ax1, ax2, ax3), results, titles, cms, units):
mc = ndimage.zoom(mc,3, mode='wrap', order=1)
bound = max(abs(np.min(mc)), abs(np.max(mc)))
if np.min(mc) < 0 < np.max(mc):
im = ax.contourf(mc.T, origin='lower', cmap=cm, levels=40, alpha=0.95, vmin=-bound, vmax=+bound)
tcs = [np.min(mc), 0, np.max(mc)]
else:
im = ax.contourf(mc.T, origin='lower', cmap=cm, levels=40, alpha=0.95)
ax.set_aspect('equal')
ax.set_title(title, fontsize=12)
ax.axis('off')
cbar=plt.colorbar(im, ticks=[np.min(mc), 0, np.max(mc)] if np.min(mc) < 0 < np.max(mc) else [np.min(mc), np.max(mc)], ax=ax, orientation='horizontal', pad=0.05, aspect=15)
cbar.ax.tick_params(labelsize=7, width=0.5)
cbar.set_label(unit, fontsize=9, labelpad=2)
plt.tight_layout()
<>:9: SyntaxWarning: "\A" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\A"? A raw string is also an option.
<>:9: SyntaxWarning: "\A" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\A"? A raw string is also an option.
/tmp/ipykernel_5677/3430547848.py:9: SyntaxWarning: "\A" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\A"? A raw string is also an option.
units=['$Z$ $(\AA)$','$H$ (Å$^{-1})$', '$K$ (Å$^{-2})$']
results_fourier = [nhaa_surface_lower_leaflet_fourier,
nhaa_mean_lower_leaflet_fourier,
nhaa_gaussian_lower_leaflet_fourier]
plots_by_leaflet(results_fourier)
The plot above allow us to directly compare the averaged results obtained for the upper leaflet using MembraneCurvature.
Imshow plots
As an alternative, we can also plot results from MembraneCurvature via imshow. By using imshow, the visualization is generated by plotting each element of the array in a matrix of m x n elements and according to a colormap of reference. In the case of MembraneCurvature, the matrix has the same shape as the arrays stored in the .results attributes. So it will be a matrix of n_x_bins, n_y_bins. The color of each square is determined by the value of the corresponding array element and the color map used.
The most basic plot using imshow to plot the surface in the lower leaflet can be obtained by:
fig, ax = plt.subplots(1,1)
ax.imshow(nhaa_surface_lower_leaflet_fourier.T, origin='lower', cmap='YlGnBu_r')
ax.set_aspect('equal')
ax.set_title('Basic plot: Surface Upper Leaflet')
plt.show()
The basic imshow plot to visualize the surface in the upper leaflet show a different color for each bin in the array. From this plot is very easy to identify we have 12 bins in each dimension.
This plot, however, is not visually pleasing. We can improve an imshow plot by adding an interpolation method. For consistency with the contour plots, we are going to use the 'gaussian' interpolation method. For more inteprolation methods you can read the imshow interpolation Matplotlib docs.
Hence, we can generate better imshow plots by doing something like:
surfaces = [nhaa_surface_lower_leaflet_fourier, # surface in lower leaflet from frame=0
nhaa_surface_upper_leaflet_fourier] # surface fn upper leaflet from frame=0
fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(4,3.5), dpi=200)
for ax, surfs, lf in zip((ax1, ax2), surfaces, leaflets):
im = ax.imshow(surfs.T, origin='lower', interpolation='gaussian', cmap='YlGnBu')
ax.set_aspect('equal')
ax.set_title('{} Leaflet'.format(lf), fontsize=6)
ax.axis('off')
cbar = plt.colorbar(im, ticks=[surfs.min(), surfs.max()], orientation='horizontal', ax=ax, shrink=0.7, aspect=10, pad=0.05)
cbar.set_ticklabels([int(surfs.min()), int(surfs.max())])
cbar.ax.tick_params(labelsize=5, width=0.5)
cbar.set_label("Height lipid headgroups (${\AA}$)", fontsize=6, labelpad=2)
<>:14: SyntaxWarning: "\A" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\A"? A raw string is also an option.
<>:14: SyntaxWarning: "\A" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\A"? A raw string is also an option.
/tmp/ipykernel_5677/1887779788.py:14: SyntaxWarning: "\A" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\A"? A raw string is also an option.
cbar.set_label("Height lipid headgroups (${\AA}$)", fontsize=6, labelpad=2)
To plot curvature using imshow, we can define a function to plot each leaflet based on this plot:
def plot_by_leaflet(results, label, cm):
fig, [ax1, ax2] = plt.subplots(ncols=2, figsize=(4,2), dpi=200)
for ax, rs, lf in zip((ax1, ax2), results, leaflets):
rs = ndimage.zoom(rs, 4, mode='wrap')
im = ax.imshow(rs.T, origin='lower', interpolation='gaussian', cmap=cm, alpha=0.8)
ax.set_aspect('equal')
ax.set_title('{} Leaflet'.format(lf), fontsize=6)
ax.axis('off')
cbar = plt.colorbar(im, ticks=[], orientation='horizontal', ax=ax, shrink=0.7)
cbar.set_label(label, fontsize=6, labelpad=2)
return
plot_by_leaflet(mean_curvs_fourier, "$H$ (Å$^{-1}$)", "bwr")
Similarly, we can plot the values stored in the variables nhaa_gaussian_lower_leaflet and nhaa_gaussian_upper_leaflet using imshow. We can plot the results using the function plot_by_leaflet.
plot_by_leaflet(gaussian_curvs_fourier, "$K$ (Å$^{-2}$)", "PiYG")
Appendix 1
To determine if the NhaA antiporter diffusses along the membrane, we can make use of the MSD MDAnalysis analysis module.
import MDAnalysis.analysis.msd as msd
Given that our universe contains a membrane-protein system, we can calculate the Mean Square Displacement (MSD) of the protein by using the class EinsteinMSD with:
MSD = msd.EinsteinMSD(universe,
select='name CA', # Select the backbone of the protein
msd_type='xy', # select plane xy
fft=True)
MSD.run(stop=1500)
0%| | 0/752 [00:00<?, ?it/s]
44%|████▍ | 329/752 [00:00<00:00, 3286.11it/s]
89%|████████▉ | 673/752 [00:00<00:00, 3371.16it/s]
100%|██████████| 752/752 [00:00<00:00, 3350.61it/s]
<MDAnalysis.analysis.msd.EinsteinMSD at 0x7be13041ee40>
To access the results from the MSD analysis, we check the .results attribute of MSD:
msd = MSD.results.timeseries
Additionaly, we can define the lagtimes based on the number of frames included in the trajectory. This will help to plot the MSD results more conveniently.
n_frames = MSD.n_frames
lagtimes = np.arange(n_frames)
And we define the plot_msd_protein function to plot the MSD of the protein along the n_frames of the simulation:
def plot_msd_protein(lagtime, msd):
fig, ax = plt.subplots()
ax.plot(lagtimes, msd, label='Protein CA')
plt.legend()
plt.xlabel('Lagtime $\Delta t$')
plt.ylabel('MSD')
plt.grid()
plt.show()
<>:5: SyntaxWarning: "\D" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\D"? A raw string is also an option.
<>:5: SyntaxWarning: "\D" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\D"? A raw string is also an option.
/tmp/ipykernel_5677/2003006459.py:5: SyntaxWarning: "\D" is an invalid escape sequence. Such sequences will not work in the future. Did you mean "\\D"? A raw string is also an option.
plt.xlabel('Lagtime $\Delta t$')
plot_msd_protein(lagtimes, msd)
The plot above illustrates the MSD of the protein backbone with respect to the lag-time. It means that the backbone of the protein diffuses in the plane of the membrane (xy plane), and therefore, additional trajectory processing is required to investigate membrane curvature around the protein as well as curvature induced by protein insertion.
Appendix 2
In section 3. Select Atoms of Reference, we used the leaflet MDAnalysis module to automatically identify residues in each leaflet.
We ran:
L_nhaa = LeafletFinder(universe, 'name P', cutoff=20)
nhaa_upper_leaflet = L_nhaa.groups(0) # upper leaflet
nhaa_lower_leaflet = L_nhaa.groups(1) # lower leafet
which allow us to obtain the atom residues in each leaflet.
nhaa_lower_leaflet.residues.resids
array([425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437,
438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450,
451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463,
620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632,
633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645,
646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658,
659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 670, 671,
672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684,
685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697,
698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710,
711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723,
724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736,
737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749,
750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762,
763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775])
nhaa_upper_leaflet.residues.resids
array([386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398,
399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411,
412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424,
464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476,
477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489,
490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502,
503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515,
516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528,
529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541,
542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554,
555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567,
568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580,
581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593,
594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606,
607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619])
However, the list of residues shown above are not consecutive. For example, residues in the lower leafet (stored in nhaa_lower_leafet), go from 425 to 463, and from 620 to 775.
An alternative approach to automatically identify the residues from each leaflet is by using split and ediff1d from the NumPy package. With np.split and np.ediff1d we can store the atom indexes of each leaflet in the membrane_indexes dictionary as in:
leaflets = ["Lower", "Upper"]
lfs = [nhaa_lower_leaflet,
nhaa_upper_leaflet]
membrane_indexes = {key:[] for key in leaflets}
for leaflet, index in zip(leaflets, range(len(leaflets))):
membrane_indexes[leaflet] = np.split(lfs[index].resids, np.where(np.ediff1d(lfs[index].residues.resids) > 1)[0] + 1)
where 'Lower' and 'Upper' are the names assigned to the keys of the membrane_indexes dictionary.
membrane_indexes
{'Lower': [array([425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437,
438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450,
451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463]),
array([620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632,
633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645,
646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658,
659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 670, 671,
672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684,
685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697,
698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710,
711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723,
724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736,
737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749,
750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762,
763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775])],
'Upper': [array([386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398,
399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411,
412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424]),
array([464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476,
477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489,
490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502,
503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515,
516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528,
529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541,
542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554,
555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567,
568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580,
581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593,
594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606,
607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619])]}
Now we can print the range for each leaflet:
for lf in leaflets:
print("\n{} leaflet:".format(lf))
for index in membrane_indexes[lf]:
print("resid {}-{} ".format(index[0],
index[-1]))
Lower leaflet:
resid 425-463
resid 620-775
Upper leaflet:
resid 386-424
resid 464-619