4.2.2. Native contacts analysis — MDAnalysis.analysis.contacts
¶
This module contains classes to analyze native contacts Q over a trajectory. Native contacts of a conformation are contacts that exist in a reference structure and in the conformation. Contacts in the reference structure are always defined as being closer then a distance radius. The fraction of native contacts for a conformation can be calculated in different ways. This module supports 3 different metrics listed below, as well as custom metrics.
 Hard Cut: To count as a contact the atoms i and j have to be at least as close as in the reference structure.
 Soft Cut: The atom pair i and j is assigned based on a soft potential
that is 1 if the distance is 0, 1/2 if the distance is the same as in
the reference and 0 for large distances. For the exact definition of the
potential and parameters have a look at function
soft_cut_q()
.  Radius Cut: To count as a contact the atoms i and j cannot be further apart than some distance radius.
The “fraction of native contacts” Q(t) is a number between 0 and 1 and calculated as the total number of native contacts for a given time frame divided by the total number of contacts in the reference structure.
4.2.2.1. Examples for contact analysis¶
4.2.2.1.1. Onedimensional contact analysis¶
As an example we analyze the opening (“unzipping”) of salt bridges when the AdK enzyme opens up; this is one of the example trajectories in MDAnalysis.
import MDAnalysis as mda
from MDAnalysis.analysis import contacts
from MDAnalysis.tests.datafiles import PSF,DCD
import matplotlib.pyplot as plt
# example trajectory (transition of AdK from closed to open)
u = mda.Universe(PSF,DCD)
# crude definition of salt bridges as contacts between NH/NZ in ARG/LYS and
# OE*/OD* in ASP/GLU. You might want to think a little bit harder about the
# problem before using this for real work.
sel_basic = "(resname ARG LYS) and (name NH* NZ)"
sel_acidic = "(resname ASP GLU) and (name OE* OD*)"
# reference groups (first frame of the trajectory, but you could also use a
# separate PDB, eg crystal structure)
acidic = u.select_atoms(sel_acidic)
basic = u.select_atoms(sel_basic)
# set up analysis of native contacts ("salt bridges"); salt bridges have a
# distance <6 A
ca1 = contacts.Contacts(u, selection=(sel_acidic, sel_basic),
refgroup=(acidic, basic), radius=6.0)
# iterate through trajectory and perform analysis of "native contacts" Q
ca1.run()
# print number of averave contacts
average_contacts = np.mean(ca1.timeseries[:, 1])
print('average contacts = {}'.format(average_contacts))
# plot time series q(t)
f, ax = plt.subplots()
ax.plot(ca1.timeseries[:, 0], ca1.timeseries[:, 1])
ax.set(xlabel='frame', ylabel='fraction of native contacts',
title='Native Contacts, average = {:.2f}'.format(average_contacts))
fig.show()
The first graph shows that when AdK opens, about 20% of the salt bridges that existed in the closed state disappear when the enzyme opens. They open in a stepwise fashion (made more clear by the movie AdK_zipper_cartoon.avi).
Notes
Suggested cutoff distances for different simulations
 For allatom simulations, cutoff = 4.5 Å
 For coarsegrained simulations, cutoff = 6.0 Å
4.2.2.1.2. Twodimensional contact analysis (q1q2)¶
Analyze a single DIMS transition of AdK between its closed and open conformation and plot the trajectory projected on q1q2 [Franklin2007]
import MDAnalysis as mda
from MDAnalysis.analysis import contacts
from MDAnalysisTests.datafiles import PSF, DCD
u = mda.Universe(PSF, DCD)
q1q2 = contacts.q1q2(u, 'name CA', radius=8)
q1q2.run()
f, ax = plt.subplots(1, 2, figsize=plt.figaspect(0.5))
ax[0].plot(q1q2.timeseries[:, 0], q1q2.timeseries[:, 1], label='q1')
ax[0].plot(q1q2.timeseries[:, 0], q1q2.timeseries[:, 2], label='q2')
ax[0].legend(loc='best')
ax[1].plot(q1q2.timeseries[:, 1], q1q2.timeseries[:, 2], '.')
f.show()
Compare the resulting pathway to the MinActionPath result for AdK [Franklin2007].
4.2.2.1.3. Writing your own contact analysis¶
The Contacts
class has been designed to be extensible for your own
analysis. As an example we will analyze when the acidic and basic groups of AdK
are in contact which each other; this means that at least one of the contacts
formed in the reference is closer than 2.5 Å.
For this we define a new function to determine if any contact is closer than
2.5 Å; this function must implement the API prescribed by Contacts
:
def is_any_closer(r, r0, dist=2.5):
return np.any(r < dist)
The first two parameters r and r0 are provided by Contacts
when it
calls is_any_closer()
while the others can be passed as keyword args
using the kwargs parameter in Contacts
.
Next we are creating an instance of the Contacts
class and use the
is_any_closer()
function as an argument to method and run the analysis:
# crude definition of salt bridges as contacts between NH/NZ in ARG/LYS and
# OE*/OD* in ASP/GLU. You might want to think a little bit harder about the
# problem before using this for real work.
sel_basic = "(resname ARG LYS) and (name NH* NZ)"
sel_acidic = "(resname ASP GLU) and (name OE* OD*)"
# reference groups (first frame of the trajectory, but you could also use a
# separate PDB, eg crystal structure)
acidic = u.select_atoms(sel_acidic)
basic = u.select_atoms(sel_basic)
nc = contacts.Contacts(u, selection=(sel_acidic, sel_basic),
method=is_any_closer,
refgroup=(acidic, basic), kwargs={'dist': 2.5})
nc.run()
bound = nc.timeseries[:, 1]
frames = nc.timeseries[:, 0]
f, ax = plt.subplots()
ax.plot(frames, bound, '.')
ax.set(xlabel='frame', ylabel='is Bound',
ylim=(0.1, 1.1))
f.show()
4.2.2.2. Functions¶

MDAnalysis.analysis.contacts.
hard_cut_q
(r, cutoff)[source]¶ Calculate fraction of native contacts Q for a hard cut off.
The cutoff can either be a float or a
ndarray
of the same shape as r.Parameters:  r (ndarray) – distance matrix
 cutoff (ndarray  float) – cut off value to count distances. Can either be a float of a ndarray of the same size as distances
Returns: Q – fraction of contacts
Return type:

MDAnalysis.analysis.contacts.
soft_cut_q
(r, r0, beta=5.0, lambda_constant=1.8)[source]¶ Calculate fraction of native contacts Q for a soft cut off
The native contact function is defined as [Best2013]
\[Q(r, r_0) = \frac{1}{1 + e^{\beta (r  \lambda r_0)}}\]Reasonable values for different simulation types are
 All Atom: lambda_constant = 1.8 (unitless)
 Coarse Grained: lambda_constant = 1.5 (unitless)
Parameters: Returns: Q – fraction of native contacts
Return type: References
[Best2013] RB Best, G Hummer, and WA Eaton, “Native contacts determine protein folding mechanisms in atomistic simulations” _PNAS_ 110 (2013), 17874–17879. doi: 10.1073/pnas.1311599110.

MDAnalysis.analysis.contacts.
radius_cut_q
(r, r0, radius)[source]¶ calculate native contacts Q based on the single distance radius.
Parameters: Returns: Q – fraction of contacts
Return type: References
[Franklin2007] (1, 2, 3) Franklin, J., Koehl, P., Doniach, S., & Delarue, M. (2007). MinActionPath: Maximum likelihood trajectory for largescale structural transitions in a coarsegrained locally harmonic energy landscape. Nucleic Acids Research, 35(SUPPL.2), 477–482. doi: 10.1093/nar/gkm342

MDAnalysis.analysis.contacts.
contact_matrix
(d, radius, out=None)[source]¶ calculate contacts from distance matrix
Parameters: Returns: contacts – boolean array of formed contacts
Return type: ndarray

MDAnalysis.analysis.contacts.
q1q2
(u, selection='all', radius=4.5, start=None, stop=None, step=None)[source]¶ Perform a q1q2 analysis.
Compares native contacts between the starting structure and final structure of a trajectory [Franklin2007].
Parameters:  u (Universe) – Universe with a trajectory
 selection (string, optional) – atoms to do analysis on
 radius (float, optional) – distance at which contact is formed
 start (int, optional) – First frame of trajectory to analyse, Default: 0
 stop (int, optional) – Last frame of trajectory to analyse, Default: 1
 step (int, optional) – Step between frames to analyse, Default: 1
Returns: contacts – Contact Analysis that is set up for a q1q2 analysis
Return type:
4.2.2.3. Classes¶

class
MDAnalysis.analysis.contacts.
Contacts
(u, selection, refgroup, method='hard_cut', radius=4.5, kwargs=None, **basekwargs)[source]¶ Calculate contacts based observables.
The standard methods used in this class calculate the fraction of native contacts Q from a trajectory.
Contact API
By defining your own method it is possible to calculate other observables that only depend on the distances and a possible reference distance. The Contact API prescribes that this method must be a function with call signature
func(r, r0, **kwargs)
and must be provided in the keyword argument method.
timeseries
¶ list – list containing Q for all refgroup pairs and analyzed frames
Parameters:  u (Universe) – trajectory
 selection (tuple(string, string)) – two contacting groups that change over time
 refgroup (tuple(AtomGroup, AtomGroup)) – two contacting atomgroups in their reference conformation. This can also be a list of tuples containing different atom groups
 radius (float, optional (4.5 Angstroms)) – radius within which contacts exist in refgroup
 method (string  callable (optional)) – Can either be one of
['hard_cut' , 'soft_cut']
or a callable with call signaturefunc(r, r0, **kwargs)
(the “Contacts API”).  kwargs (dict, optional) – dictionary of additional kwargs passed to method. Check respective functions for reasonable values.
 verbose (bool (optional)) – Show detailed progress of the calculation if set to
True
; the default isFalse
.
