# 4.3.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 than 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.

1. Hard Cut: To count as a contact the atoms i and j have to be at least as close as in the reference structure.

2. 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()`.

3. 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.3.2.1. Examples for contact analysis

### 4.3.2.1.1. One-dimensional 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 numpy as np
    import matplotlib.pyplot as plt
    import MDAnalysis as mda
    from MDAnalysis.analysis import contacts
    from MDAnalysis.tests.datafiles import PSF,DCD
    # 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, select=(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.results.timeseries[:, 1])
    print('average contacts = {}'.format(average_contacts))
    # plot time series q(t)
    fig, ax = plt.subplots()
    ax.plot(ca1.results.timeseries[:, 0], ca1.results.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 step-wise fashion (made more clear by the movie AdK_zipper_cartoon.avi).

Notes

Suggested cutoff distances for different simulations

- For all-atom simulations, cutoff = 4.5 Å

- For coarse-grained simulations, cutoff = 6.0 Å

### 4.3.2.1.2. Two-dimensional contact analysis (q1-q2)

Analyze a single DIMS transition of AdK between its closed and open conformation and plot the trajectory projected on q1-q2 [1]

    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.results.timeseries[:, 0], q1q2.results.timeseries[:, 1],
               label='q1')
    ax[0].plot(q1q2.results.timeseries[:, 0], q1q2.results.timeseries[:, 2],
               label='q2')
    ax[0].legend(loc='best')
    ax[1].plot(q1q2.results.timeseries[:, 1],
               q1q2.results.timeseries[:, 2], '.-')
    f.show()

Compare the resulting pathway to the MinActionPath result for AdK [1].

### 4.3.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, select=(sel_acidic, sel_basic),
                           method=is_any_closer,
                           refgroup=(acidic, basic), kwargs={'dist': 2.5})
    nc.run()

    bound = nc.results.timeseries[:, 1]
    frames = nc.results.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.3.2.2. Functions

MDAnalysis.analysis.contacts.hard_cut_q(r, cutoff)

  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:

    float

MDAnalysis.analysis.contacts.soft_cut_q(r, r0, beta=5.0, lambda_constant=1.8)

  Calculate fraction of native contacts Q for a soft cut off

  The native contact function is defined as [2]

  \[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:

    - r (array) – Contact distances at time t

    - r0 (array) – Contact distances at time t=0, reference distances

    - beta (float (default 5.0 Angstrom)) – Softness of the switching function

    - lambda_constant (float (default 1.8, unitless)) – Reference distance tolerance

  Returns:

    Q – fraction of native contacts

  Return type:

    float

MDAnalysis.analysis.contacts.radius_cut_q(r, r0, radius)

  calculate native contacts Q based on the single distance radius.

  Parameters:

    - r (ndarray) – distance array between atoms

    - r0 (ndarray) – unused to fullfill `Contacts` API

    - radius (float) – Distance between atoms at which a contact is formed

  Returns:

    Q – fraction of contacts

  Return type:

    float

MDAnalysis.analysis.contacts.contact_matrix(d, radius, out=None)

  calculate contacts from distance matrix

  Parameters:

    - d (array-like) – distance matrix

    - radius (float) – distance below which a contact is formed.

    - out (array (optional)) – If out is supplied as a pre-allocated array of the correct shape then it is filled instead of allocating a new one in order to increase performance.

  Returns:

    contacts – boolean array of formed contacts

  Return type:

    ndarray

MDAnalysis.analysis.contacts.q1q2(u, select='all', radius=4.5)

  Perform a q1-q2 analysis.

  Compares native contacts between the starting structure and final structure of a trajectory [1].

  Parameters:

    - u (Universe) – Universe with a trajectory

    - select (string, optional) – atoms to do analysis on

    - radius (float, optional) – distance at which contact is formed

  Returns:

    contacts – Contact Analysis that is set up for a q1-q2 analysis

  Return type:

    `Contacts`

  Changed in version 1.0.0: Changed selection keyword to select Support for setting `start`, `stop`, and `step` has been removed. These should now be directly passed to `Contacts.run()`.

## 4.3.2.3. Classes

class MDAnalysis.analysis.contacts.Contacts(u, select, refgroup, method='hard_cut', radius=4.5, pbc=True, kwargs=None, **basekwargs)

  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.

  results.timeseries

    2D array containing Q for all refgroup pairs and analyzed frames

    Type:

      numpy.ndarray

  timeseries

    Alias to the `results.timeseries` attribute.

    Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please use `results.timeseries` instead.

    Type:

      numpy.ndarray

  Changed in version 1.0.0: `save()` method has been removed. Use `np.savetxt()` on `Contacts.results.timeseries` instead.

  Changed in version 1.0.0: added `pbc` attribute to calculate distances using PBC.

  Changed in version 2.0.0: `timeseries` results are now stored in a `MDAnalysis.analysis.base.Results` instance.

  Changed in version 2.2.0: `Contacts` accepts both AtomGroup and string for select

  Changed in version 2.9.0: Introduced `get_supported_backends()` allowing for parallel execution on `multiprocessing` and `dask` backends.

  Parameters:

    - u (Universe) – trajectory

    - select (tuple(AtomGroup, AtomGroup) | 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', 'radius_cut']` or a callable with call signature `func(r, r0, **kwargs)` (the “Contacts API”).

    - pbc (bool (optional)) – Uses periodic boundary conditions to calculate distances if set to `True`; the default is `True`.

    - 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 is `False`.

  n_initial_contacts

    Total number of initial contacts.

    Type:

      int

  r0

    List of distance arrays between reference groups.

    Type:

      list[numpy.ndarray]

  Notes

  Changed in version 1.0.0: Changed selection keyword to select

  classmethod get_supported_backends()

    Tuple with backends supported by the core library for a given class. User can pass either one of these values as `backend=...` to `run()` method, or a custom object that has `apply` method (see documentation for `run()`):

    > - ‘serial’: no parallelization

    > - ‘multiprocessing’: parallelization using multiprocessing.Pool

    > - ‘dask’: parallelization using dask.delayed.compute(). Requires installation of mdanalysis[dask]

    If you want to add your own backend to an existing class, pass a `backends.BackendBase` subclass (see its documentation to learn how to implement it properly), and specify `unsupported_backend=True`.

    Returns:

      names of built-in backends that can be used in `run(backend=...)()`

    Return type:

      tuple

    Added in version 2.8.0.

References

[1] (1,2,3)

Joel Franklin, Patrice Koehl, Sebastian Doniach, and Marc Delarue. MinActionPath: maximum likelihood trajectory for large-scale structural transitions in a coarse-grained locally harmonic energy landscape. Nucleic Acids Research, 35(suppl_2):W477–W482, 2007. doi:10.1093/nar/gkm342.

[2]

Robert B. Best, Gerhard Hummer, and William A. Eaton. Native contacts determine protein folding mechanisms in atomistic simulations. Proceedings of the National Academy of Sciences, 110(44):17874–17879, 2013. doi:10.1073/pnas.1311599110.