How to calculate residue‑membrane contact profile (residue index vs contact count) in GROMACS‑2018 +

User discussions
1

GROMACS version:
GROMACS modification: Yes/No
Here post your question
Hello everyone,

I am working on coarse‑grained MD simulations (Martini3 force field, GROMACS‑2018).
I want to compute a residue‑membrane contact profile plot:

  • X‑axis: residue index of protein (1‑185 for K‑Ras4B)
  • Y‑axis: contact count / contact probability between each protein residue and inner‑leaflet lipid head‑groups.

Although gmx mindist can output results on a per‑residue basis, it only outputs minimum distance values instead of contact counts.

2

Hi,

gmx mindist does have an option to output contact counts (-on). If you tried that and, for some reason, it did not work, then I recommend checking out the PyLipID library for protein-membrane interaction analysis; it also does exactly what you want: 0-Application walk thourgh - PyLipID documentation

Cheers!

3

I ended up having to use MDAnalysis for this task because there’s no Gromacs tool which provides information about each contact. Here’s some code which might help:

import MDAnalysis as mda
import numpy as np

# Define universe, selection groups, and a list of headgroup atom ids
u = mda.Universe(top, traj)        
protein_atoms = u.select_atoms('<some select string>')
head_atoms = u.select_atoms('<some select string>'
dist_to_head_aid  = np.array([a.id for a in head_atoms], dtype=int)


# Compute contacts
cutoff = 3.5
print("Calculating distances...")
all_distances = np.empty((len(protein_atoms), len(head_atoms)))
for i, frame in enumerate(u.trajectory):
    # Get all distances between protein and lipid heads
    # This will parallelize over however many processes are available
    distances.distance_array(protein_atoms, head_atoms, box=frame.dimensions, result=all_distances, backend="OpenMP")

    # Get atom IDs for all atoms where the distance array is less than the cutoff
    indices = np.nonzero(all_distances<cutoff)
    for j in range(len(indices[0])):
        protein_id = indices[0][j]
        lipid = dist_to_head_aid[indices[1][j]]
        dist = all_distances[(protein_id, indices[1][j])]
        # I used this to build a database of contacts, but you could just accumlate into an array of len(protein_atoms)