This page was generated from examples/pi_stacking.ipynb.

Pi-Stacking Identification

Pi Stacking Analysis on PDB 4BFQ

In this first example, we show pi stacking analysis on the PDB structure 4BFQ. 4BFQ has two pi-stacking interactions between a ligand (resSeq 301) and Y186.

Both aromatic residues on the ligand are placed in such a way that the TYR phenol stacks with them.

Before we show the MDTraj analysis, we will use py3dmol to visualize the structure.

In the visualization below, we will highlight the following:

• In silver: The protein

• In light steel blue: Ligands with name “083”

• In green: The tyrosine residue on the protein.

• In magenta: One ring on an 083 ligand we will be analyzing. The atoms are C16, C17, C18, C19, C20, C21. (pi stacks with TYR)

• In cyan: Another ring on an 083 ligand we will be analyzing. The atoms are N4 C9 C15, C16, C21, and C8. (pi stacks with TYR)

• In yellow: The third ring on an 083 ligand we will be analyzing. The atoms are N5, C10, C11, C12, C13, C14. (does not stack with TYR)

import py3Dmol

# Load the PDB structure from the RCSB PDB database
viewer = py3Dmol.view(query="pdb:4BFQ", width=800, height=600)

# Style the entire protein as silver cartoon
viewer.setStyle({"cartoon": {"color": "silver"}})

# Make residue 186 on chain E visible as bright green sticks
viewer.addStyle({"chain": "E", "resi": "186"}, {"stick": {"color": "#00FF00", "radius": 0.2}})

# Make all of residue 083 visible as sticks
viewer.addStyle({"chain": "E", "resn": "083"}, {"stick": {"color": "lightsteelblue", "radius": 0.15}})

# First highlight group in magenta (pi-stacking)
viewer.addStyle(
    {"chain": "E", "resn": "083", "resi": "301", "atom": ["C16", "C17", "C18", "C19", "C20", "C21"]},
    {"stick": {"color": "magenta", "radius": 0.2}},
)

# Second highlight group in cyan (pi-stacking)
viewer.addStyle(
    {"chain": "E", "resn": "083", "resi": "301", "atom": ["N4", "C9", "C15", "C16", "C21", "C8"]},
    {"stick": {"color": "cyan", "radius": 0.2}},
)

# Third highlight group in yellow (NOT pi-stacking)
viewer.addStyle(
    {"chain": "E", "resn": "083", "resi": "301", "atom": ["N5", "C10", "C11", "C12", "C13", "C14"]},
    {"stick": {"color": "yellow", "radius": 0.2}},
)

# Set the zoom level to focus on our selections
viewer.zoomTo({"chain": "E", "resn": "083", "resi": "301"})

viewer.rotate(90, "y")
viewer.rotate(90, "x")
# Show the viewer
viewer.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

The rings highlighted in magenta and cyan are expected to stack with the TYR residue, while the ring in yellow is not expected to stack with the TYR residue.

import mdtraj as md

Next, we load the structure from the PDB using MDTraj.

t = md.load_pdb("/build/reproducible-path/mdtraj-1.11.0/debian/external/4BFQ.pdb")
print(t)

<mdtraj.Trajectory with 1 frames, 8693 atoms, 1169 residues, and unitcells>

The md.pi_stacking method requires pre-identification of the aromatic residues that one would like to measure pi-stacks between. As such, one must identify aromatic substructures on both, e.g. a ligand of interest and a protein of interest.

We can use the selection language to identify the atom indices of the particular ligand atom indices:

lig_aromatic_atms_to_test = [
    ("C16", "C17", "C18", "C19", "C20", "C21"),  # pi-stacking
    ("N4", "C9", "C15", "C16", "C21", "C8"),  # pi-stacking
    ("N5", "C10", "C11", "C12", "C13", "C14"),  # NOT pi-stacking
]
lig_grps = []
for lig_grp in lig_aromatic_atms_to_test:
    lig_grps.append(
        tuple(
            int(idx)
            for idx in t.top.select(
                f"chainid 5 and resSeq 301 and (name {lig_grp[0]} or name {' or name '.join(lig_grp[1:])})"
            )
        )
    )
print(f"{lig_grps=}")

lig_grps=[(8192, 8193, 8194, 8195, 8196, 8197), (8191, 8192, 8197, 8198, 8199, 8200), (8201, 8202, 8203, 8204, 8205, 8206)]

Let’s also get the atom indices of the aromatic atoms in the sidechain of TYR 186

protein_grps = []
# TYR186
protein_stacking_atms = [("CE2", "CD2", "CG", "CD1", "CE1", "CZ")]
for protein_grp in protein_stacking_atms:
    protein_grps.append(
        tuple(
            int(idx)
            for idx in t.top.select(
                f"chainid 0 and resSeq 186 and (name {protein_grp[0]} or name {' or name '.join(protein_grp[1:])})"
            )
        )
    )
print(f"{protein_grps=}")

protein_grps=[(1472, 1473, 1474, 1475, 1476, 1477)]

Now we can calculate if there is a stacking interaction with each of our protein/ligand groups:

stacking_interactions = md.pi_stacking(t, lig_grps, protein_grps)
print(f"{stacking_interactions=}")
for frame_num, frame in enumerate(stacking_interactions):
    print(f"{frame_num=}")
    for lig_grp, protein_grp in frame:
        print(f"{[t.top.atom(atm) for atm in lig_grp]} <-> {[t.top.atom(atm) for atm in protein_grp]}")

stacking_interactions=[[((8192, 8193, 8194, 8195, 8196, 8197), (1472, 1473, 1474, 1475, 1476, 1477)), ((8191, 8192, 8197, 8198, 8199, 8200), (1472, 1473, 1474, 1475, 1476, 1477))]]
frame_num=0
[083301-C21, 083301-C20, 083301-C19, 083301-C18, 083301-C17, 083301-C16] <-> [TYR186-CG, TYR186-CD1, TYR186-CD2, TYR186-CE1, TYR186-CE2, TYR186-CZ]
[083301-C8, 083301-C21, 083301-C16, 083301-C15, 083301-C9, 083301-N4] <-> [TYR186-CG, TYR186-CD1, TYR186-CD2, TYR186-CE1, TYR186-CE2, TYR186-CZ]

Note how the ligand group ("N5", "C10", "C11", "C12", "C13", "C14") was not found, even though we supplied it as a potential stacking aromatic group, since it didn’t meet the geometric criteria for pi-stacking.

Analyzing T-Stacking Interactions

4BFQ has clear face-to-face pi-stacking interactions. T-stacking (edge-to-face) interactions are also supported in md.pi_stacking.

Before showing this analysis, we will visualize the structure again using py3Dmol.

import py3Dmol

# Load the PDB structure from the RCSB PDB database
viewer = py3Dmol.view(query="pdb:6A22", width=800, height=600)

# Style the entire protein as silver cartoon
viewer.setStyle({"cartoon": {"color": "silver"}})

# Make the whole ligand visible as sticks
viewer.addStyle({"chain": "C", "resi": "9000"}, {"stick": {"color": "lightsteelblue", "radius": 0.2}})

# First ligand group (Ring 1) in magenta (t-stacking, edge-to-face with PHE)
viewer.addStyle(
    {"chain": "C", "resi": "9000", "atom": ["S10", "C9", "C8", "N7", "C6"]},
    {"stick": {"color": "magenta", "radius": 0.2}},
)

# Second ligand group (Ring 2) in yellow (NOT pi-stacking)
viewer.addStyle(
    {"chain": "C", "resi": "9000", "atom": ["C14", "C15", "C16", "C17", "C18", "C19"]},
    {"stick": {"color": "yellow", "radius": 0.2}},
)

# Highlight the PHE residue 378 in green
viewer.addStyle({"chain": "C", "resi": "378"}, {"stick": {"color": "#00FF00", "radius": 0.2}})

# Zoom specifically to the green PHE residue
viewer.zoomTo({"chain": "C", "resi": "9000"})

# Apply rotation after zooming
viewer.rotate(95, "y")
viewer.rotate(150, "x")

# Show the viewer
viewer.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

t = md.load_pdb("/build/reproducible-path/mdtraj-1.11.0/debian/external/6A22.pdb")
print(t)

<mdtraj.Trajectory with 1 frames, 8066 atoms, 1064 residues, and unitcells>

lig_aromatic_atms_to_test = [
    ("S10", "C9", "C8", "N7", "C6"),  # t-stacking, edge-to-face with PHR, shown in magenta in the visualization
    ("C14", "C15", "C16", "C17", "C18", "C19"),  # NOT pi-stacking shown in yellow in the visualization
]
lig_grps = []
for lig_grp in lig_aromatic_atms_to_test:
    lig_grps.append(
        tuple(
            int(idx)
            for idx in t.top.select(
                f"chainid 9 and resSeq 9000 and (name {lig_grp[0]} or name {' or name '.join(lig_grp[1:])})"
            )
        )
    )
protein_grps = []
# PHR118/378
protein_stacking_atms = [("CE2", "CD2", "CG", "CD1", "CE1", "CZ")]
for protein_grp in protein_stacking_atms:
    protein_grps.append(
        tuple(
            int(idx)
            for idx in t.top.select(
                f"chainid 2 and resSeq 378 and (name {protein_grp[0]} or name {' or name '.join(protein_grp[1:])})"
            )
        )
    )
print(f"{lig_grps=}")
print(f"{protein_grps=}")

lig_grps=[(7879, 7882, 7883, 7900, 7901), (7887, 7888, 7903, 7904, 7905, 7908)]
protein_grps=[(2922, 2923, 2924, 2925, 2926, 2927)]

stacking_interactions = md.pi_stacking(t, lig_grps, protein_grps)
print(f"{stacking_interactions=}")
for frame_num, frame in enumerate(stacking_interactions):
    print(f"{frame_num=}")
    for lig_grp, protein_grp in frame:
        print(f"{[t.top.atom(atm) for atm in lig_grp]} <-> {[t.top.atom(atm) for atm in protein_grp]}")

stacking_interactions=[[((7879, 7882, 7883, 7900, 7901), (2922, 2923, 2924, 2925, 2926, 2927))]]
frame_num=0
[9P69000-C6, 9P69000-C8, 9P69000-C9, 9P69000-N7, 9P69000-S10] <-> [PHE378-CG, PHE378-CD1, PHE378-CD2, PHE378-CE1, PHE378-CE2, PHE378-CZ]