Note

Go to the end to download the full example code

# Visualization of normal modes from an elastic network model#

The *elastic network model* (ENM) is a fast method to estimate movements
in a protein structure, without the need to run time-consuming MD
simulations.
A protein is modelled as *mass-and-spring* model, with the masses being
the \(C_\alpha\) atoms and the springs being the non-covalent bonds
between adjacent residues.
Via *normal mode analysis* distinct movements/oscillations can be
extracted from the model.

An *anisotropic network model* (ANM), is an ENM that includes
directional information.
Hence, the normal mode analysis yields eigenvectors, where each atom is
represented by three vector components (*x*, *y*, *z*).
Thus these vectors can be used for 3D representation.

In the case of this example a normal mode analysis on an ANM was already conducted. This script merely takes the structure and obtained eigenvectors to add a smooth oscillation of the chosen normal mode to the structure. The newly created structure has multiple models, where each model depicts a different time in the oscillation period. Then the multi-model structure can be used to create a video of the oscillation using a molecular visualization program.

The file containing the eigenvectors can be downloaded via this
`link`

.

```
# Code source: Patrick Kunzmann
# License: BSD 3 clause
from tempfile import NamedTemporaryFile
import numpy as np
from numpy import newaxis
import biotite.database.rcsb as rcsb
import biotite.structure as struc
import biotite.structure.io as strucio
import biotite.structure.io.pdbx as pdbx
# A CSV file containing the eigenvectors for the CA atoms
VECTOR_FILE = "../../../download/glycosylase_anm_vectors.csv"
# The corresponding structure
PDB_ID = "1MUG"
# The normal mode to be visualized
# '-1' is the last (and most significant) one
MODE = -1
# The amount of frames (models) per oscillation
FRAMES = 60
# The maximum oscillation amplitude for an atom
# (The length of the ANM's eigenvectors make only sense when compared
# relative to each other, the absolute values have no significance)
MAX_AMPLITUDE = 5
# Load structure
pdbx_file = pdbx.BinaryCIFFile.read(rcsb.fetch(PDB_ID, "bcif"))
structure = pdbx.get_structure(pdbx_file, model=1)
# Filter first peptide chain
protein_chain = structure[
struc.filter_amino_acids(structure) & (structure.chain_id == structure.chain_id[0])
]
# Filter CA atoms
ca = protein_chain[protein_chain.atom_name == "CA"]
# Load eigenvectors for CA atoms
# The first axis indicates the mode,
# the second axis indicates the vector component
vectors = np.loadtxt(VECTOR_FILE, delimiter=",").transpose()
# Discard the last 6 modes, as these are movements of the entire system:
# A system with N atoms has only 3N - 6 degrees of freedom
# ^^^
vectors = vectors[:-6]
# Extract vectors for given mode and reshape to (n,3) array
mode_vectors = vectors[MODE].reshape((-1, 3))
# Rescale, so that the largest vector has the length 'MAX_AMPLITUDE'
vector_lenghts = np.sqrt(np.sum(mode_vectors**2, axis=-1))
scale = MAX_AMPLITUDE / np.max(vector_lenghts)
mode_vectors *= scale
# Stepwise application of eigenvectors as smooth sine oscillation
time = np.linspace(0, 2 * np.pi, FRAMES, endpoint=False)
deviation = np.sin(time)[:, newaxis, newaxis] * mode_vectors
# Apply oscillation of CA atom to all atoms in the corresponding residue
oscillation = np.zeros((FRAMES, len(protein_chain), 3))
residue_starts = struc.get_residue_starts(
protein_chain,
# The last array element will be the length of the atom array,
# i.e. no valid index
add_exclusive_stop=True,
)
for i in range(len(residue_starts) - 1):
res_start = residue_starts[i]
res_stop = residue_starts[i + 1]
oscillation[:, res_start:res_stop, :] = (
protein_chain.coord[res_start:res_stop, :] + deviation[:, i : i + 1, :]
)
# An atom array stack containing all frames
oscillating_structure = struc.from_template(protein_chain, oscillation)
# Save as PDB for rendering a video with PyMOL
temp = NamedTemporaryFile(suffix=".pdb")
strucio.save_structure(temp.name, oscillating_structure)
temp.close()
```