Custom Gamma Functions

This module defines customized gamma functions for elastic network model analysis.

class Gamma[source]

Base class for facilitating use of atom type, residue type, or residue property dependent force constants (γ).

Derived classes:

gamma(dist2, i, j)[source]

Returns force constant.

For efficiency purposes square of the distance between interacting atom/residue (node) pairs is passed to this function. In addition, node indices are passed.

class GammaStructureBased(atoms, gamma=1.0, helix=6.0, sheet=6.0, connected=10.0)[source]

Facilitate setting the spring constant based on the secondary structure and connectivity of the residues.

A systematic study [LT10] of a large set of NMR-structures analyzed using a method based on entropy maximization showed that taking into consideration properties such as sequential separation between contacting residues and the secondary structure types of the interacting residues provides refinement in the ENM description of proteins.

This class determines pairs of connected residues or pairs of proximal residues in a helix or a sheet, and assigns them a larger user defined spring constant value.

DSSP single letter abbreviations are recognized:
  • H: α-helix
  • G: 3-10-helix
  • I: π-helix
  • E: extended part of a sheet
helix:
Applies to residue (or Cα atom) pairs that are in the same helical segment, at most 7 Å apart, and separated by at most 3 (3-10-helix), 4 (α-helix), or 5 (π-helix) residues.
sheet:
Applies to Cα atom pairs that are in different β-strands and at most 6 Å apart.
connected:
Applies to Cα atoms that are at most 4 Å apart.

Note that this class does not take into account insertion codes.

[LT10]Lezon TR, Bahar I. Using entropy maximization to understand the determinants of structural dynamics beyond native contact topology. PLoS Comput Biol 2010 6(6):e1000816.

Example:

Let’s parse coordinates and header data from a PDB file, and then assign secondary structure to the atoms.

In [1]: from prody import *

In [2]: ubi, header = parsePDB('1aar', chain='A', subset='calpha', header=True)

In [3]: assignSecstr(header, ubi)
Out[3]: <AtomGroup: 1aarA_ca (76 atoms)>

In the above we parsed only the atoms needed for this calculation, i.e. Cα atoms from chain A.

We build the Hessian matrix using structure based force constants as follows;

In [4]: gamma = GammaStructureBased(ubi)

In [5]: anm = ANM('')

In [6]: anm.buildHessian(ubi, gamma=gamma)

We can obtain the force constants assigned to residue pairs from the Kirchhoff matrix as follows:

In [7]: k = anm.getKirchhoff()

In [8]: k[0,1] # a pair of connected residues
Out[8]: -10.0

In [9]: k[0,16] # a pair of residues from a sheet
Out[9]: -6.0

Setup the parameters.

Parameters:
  • atoms (Atomic) – A set of atoms with chain identifiers, residue numbers, and secondary structure assignments are set.
  • gamma (float) – Force constant in arbitrary units. Default is 1.0.
  • helix (float) – Force constant factor for residues hydrogen bonded in α-helices, 3,10-helices, and π-helices. Default is 6.0, i.e. 6.0`*gamma.
  • sheet (float) – Force constant factor for residue pairs forming a hydrogen bond in a β-sheet. Default is 6.0, i.e. 6.0`*gamma.
  • connected (float) – Force constant factor for residue pairs that are connected. Default is 10.0, i.e. 10.0`*gamma.
gamma(dist2, i, j)[source]

Returns force constant.

getChids()[source]

Returns a copy of chain identifiers.

getResnums()[source]

Returns a copy of residue numbers.

getSecstrs()[source]

Returns a copy of secondary structure assignments.

class GammaVariableCutoff(identifiers, gamma=1.0, default_radius=7.5, **kwargs)[source]

Facilitate setting the cutoff distance based on user defined atom/residue (node) radii.

Half of the cutoff distance can be thought of as the radius of a node. This class enables setting different radii for different node types.

Example:

Let’s think of a protein-DNA complex for which we want to use different radius for different residue types. Let’s say, for protein Cα atoms we want to set the radius to 7.5 Å, and for nucleic acid phosphate atoms to 10 Å. We use the HhaI-DNA complex structure 1mht.

In [1]: hhai = parsePDB('1mht')

In [2]: ca_p = hhai.select('(protein and name CA) or (nucleic and name P)')

In [3]: ca_p.getNames()
Out[3]: 
array(['P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P',
       'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'P', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA', 'CA',
       'CA', 'CA', 'CA', 'CA', 'CA'], dtype='|S6')

We set the radii of atoms:

In [4]: varcutoff = GammaVariableCutoff(ca_p.getNames(), gamma=1,
   ...:     default_radius=7.5, debug=False, P=10)
   ...: 

In [5]: varcutoff.getRadii()
Out[5]: 
array([10. , 10. , 10. , 10. , 10. , 10. , 10. , 10. , 10. , 10. , 10. ,
       10. , 10. , 10. , 10. , 10. , 10. , 10. , 10. , 10. , 10. , 10. ,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,
        7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5,  7.5])

The above shows that for phosphate atoms radii is set to 10 Å, because we passed the P=10 argument. As for Cα atoms, the default 7.5 Å is set as the radius (default_radius=7.5). You can also try this with debug=True argument to print debugging information on the screen.

We build ANM Hessian matrix as follows:

In [6]: anm = ANM('HhaI-DNA')

In [7]: anm.buildHessian(ca_p, gamma=varcutoff, cutoff=20)

Note that we passed cutoff=20.0 to the ANM.buildHessian() method. This is equal to the largest possible cutoff distance (between two phosphate atoms) for this system, and ensures that all of the potential interactions are evaluated.

For pairs of atoms for which the actual distance is larger than the effective cutoff, the GammaVariableCutoff.gamma() method returns 0. This annuls the interaction between those atom pairs.

Set the radii of atoms.

Parameters:
  • identifiers (list or numpy.ndarray) – List of atom names or types, or residue names.
  • gamma (float) – Uniform force constant value. Default is 1.0.
  • default_radius (float) – Default radius for atoms whose radii is not set as a keyword argument. Default is 7.5

Keywords in keyword arguments must match those in atom_identifiers. Values of keyword arguments must be float.

gamma(dist2, i, j)[source]

Returns force constant.

getGamma()[source]

Returns the uniform force constant value.

getRadii()[source]

Returns a copy of radii array.

class GammaED(atoms, Ccart=6.0, Ex=6, Cseq=60.0, Slim=3)[source]

Facilitate setting the spring constant based on sequence distance and spatial distance as in ed-ENM [OL10]. This ENM is refined based on comparison to essential dynamics from MD simulations and can reproduce flexibility in NMR ensembles.

It has also been implemented in FlexServ [CJ09] and used in MDdMD [SP12] and GOdMD [SP13].

The sequence distance-dependent term is Cseq/(S**2) for S=abs(i,j) <= Slim

The structure distance-dependent term is (Ccart/dist)**Ex

[OL10]Orellana L, Rueda M, Ferrer-Costa C, Lopez-Blanco JR, Chacón P, Orozco M. Approaching Elastic Network Models to Molecular Dynamics Flexibility. J Chem Theory Comput 2010 6(9):2910-23.
[CJ09]Camps J, Carrillo O, Emperador A, Orellana L, Hospital A, Rueda M, Cicin-Sain D, D’Abramo M, Gelpí JL, Orozco M. FlexServ: an integrated tool for the analysis of protein flexibility. Bioinformatics 2009 25(13):1709-10.
[SP12]Sfriso P, Emperador A, Orellana L, Hospital A, Gelpí JL, Orozco M. Finding Conformational Transition Pathways from Discrete Molecular Dynamics Simulations. J Chem Theory Comput 2012 8(11):4707-18.
[SP13]Sfriso P, Hospital A, Emperador A, Orozco M. Exploration of conformational transition pathways from coarse-grained simulations. Bioinformatics 2013 29(16):1980-6.

Example:

Let’s parse coordinates from a PDB file.

In [1]: from prody import *

In [2]: ubi = parsePDB('1aar', chain='A', subset='calpha')

In the above we parsed only the atoms needed for this calculation, i.e. Cα atoms from chain A.

We build the Hessian matrix using GOdMD spring constants as follows;

In [3]: gamma = GammaGOdMD(ubi)

In [4]: anm = ANM('')

In [5]: anm.buildHessian(ubi, gamma=gamma)

Setup the parameters.

Parameters:
  • atoms (Atomic) – A set of atoms
  • Ccart (float) – Multiplication constant inside the exponential in the spatial distance-dependent term Default is 6.
  • Ex – Exponent in spatial distance-dependent term Default is 6
  • Cseq (float) – Multiplication constant in sequence distance-dependent term Default is 60.
  • Slim (float) – Sequence distance limit for sequence distance-dependence This limit is used with a less-or-equal operator Default is 3
gamma(dist2, i, j)[source]

Returns force constant.

GammaGOdMD

alias of GammaED