SOP model

Overview

In SOP model, each residue is described using a single interaction center (C-alpha atom). The potential energy function of a protein conformation V, specified in terms of the coordinates <math>\{\bf r\} = {\bf r}_1, {\bf r}_2, \ldots {\bf r}_N</math> , is given by

<math>V = V_{FENE}+V^{ATT}_{NB}+V^{REP}_{NB}</math>, where

<math>V_{FENE} = -\sum_{covalent}{\frac{k}{2}R_0^2\log{\left(1-\frac{\left(r_{ij}-r^0_{ij}\right)^2}{R_0^2}\right)}}</math>

<math>V^{ATT}_{NB} = \sum_{native}\varepsilon_h\left[\left(\frac{r^{0}_{ij}}{r_{ij}}\right)^{12}-2\left(\frac{r^{0}_{ij}}{r_{ij}}\right)^{6}\right]</math>

<math>V^{REP}_{NB} = \sum_{repulsive}{\varepsilon_r\left(\frac{\sigma}{r_{ij}}\right)^{6}}</math>.

The distance between two residues, i and j, is <math>r_{ij}</math>, while <math>r_{ij}^0</math> is its value in native (PDB) structure. The finite extensible nonlinear elastic (FENE) potential <math>V_{FENE}</math> describes covalent bonds (both backbone and disulfide bonds) with tolerance value in the change of covalent bond <math>R_0 = 2</math>Å.

We used Lennard-Jones potential (<math>V^{ATT}_{NB}</math>) to account for the noncovalent interactions that stabilize the native state (i.e. native interactions). Two non-covalently linked residues form native contact if, in native state (PDB structure), the distance between their C-alpha atoms is less then cut-off distance (simple Go definition) or/and centers of heavy atoms of amino-acids side chains are within some cut-off distance (full Go definition). Most common definition is a simple Go with a cut-off distance of 8Å. The value of <math>\varepsilon_h</math> quantifies the strength of non-bonded interactions.

All other pairs of atoms (i.e. that neither connected covalently nor form a native contact) are considered to have repulsive Lennard-Jones potential <math>V^{REP}_{NB}</math> (<math>\sim 1/r^6</math>). The strength or repulsive potential is described by parameter <math>\varepsilon_r = 1</math> kcal/mol, repulsive distance is equal to the length of amino-acid <math>\sigma</math> = 3.8Å.

Units

Characteristic time for underdamped motion of spherical particle of mass m and radius a with energy scale <math>\varepsilon_h</math> is <math>\tau_L = \sqrt{m a^2/\varepsilon_h}</math>. Then, <math>\alpha = \zeta m/\tau_L</math> is the friction coefficient for the spherical particle (ζ = dimensionless friction coeff). <math>\alpha = 6\pi\eta a</math> is the Stokes-Einstein friction for spherical particle of radius a (in Å) in liquid with viscosity η. Therefore <math>\zeta = (6\pi\eta a^2)/\sqrt{m\varepsilon_h}</math>. Bulk water viscosity <math>\eta = 0.01</math> gs^-1cm^-1. In the program ζ = 50.

This was obtained for <math>a \sim 5</math> Å and <math>m \sim 3 \cdot 10^{-22}</math> g (mass of a residue) and using bulk water viscosity. In general, a varies between 3.8Å to 5Å, while m varies between 3⋅10^-22 g to 5⋅10^-22 g. In simulations a = 3.8Å. Because of the fact that ζ depends on <math>\varepsilon_h</math>, every time <math>\varepsilon_h</math> changed, valid m value shold be found. This value should give ζ = 50.

Example: for <math>\varepsilon_h</math> = 1 kcal/mol from the above equation for ζ we find that <math>m = 4.3 \cdot 10^{-22}</math> g which is a valid value. For <math>\varepsilon_h</math> = 1.5 kcal/mol, we get <math>m = 3 \cdot 10^{-22}</math> g which is still a valid value. After finding the mass m, we can go back to the expression for <math>\tau_L</math> and get its value. For example, for <math>\varepsilon_h</math> = 1 kcal/mol we get <math>\tau_L = 3</math> ps while for <math>\varepsilon_h</math> = 1.5 kcal/mol, we get <math>\tau_L = 2</math> ps.

<math>\tau_H = \zeta\varepsilon_h\tau_L/kT</math> is the characteristic time for the overdamped Langevin dynamics. To get it in ps, both <math>\varepsilon_h</math> and kT need to be in the same units. Because <math>\varepsilon_h</math> is in kcal/mol, kT is also in kcal/mol (T = 300K corresponds to 0.6 kcal/mol). Therefore an elementary time step for the simulation is <math>h\times\tau_H</math> in ps.

After simplification, we get <math>\tau_H = \frac{6 \pi \eta a^3}{k T}</math> = 248 ps (assuming <math>\eta = 0.01</math> gs^-1cm^-1, T = 300 K and a = 3.8Å)

The pulling speed is <math>v_f = \Delta x/(n_{av}\times h\times\tau_H)</math> where <math>\Delta x</math> is a displacement during n_av timesteps in Å. The force is in kcal/(mol Å), to get a force in pN, value in kcal/(mol Å) should be multiplied by 70.

Indentation potential

The cantilever tip is presented by Lennard-Jones sphere:

<math>\varepsilon_l\left(A\times\left(\frac{\sigma}{r}\right)^{12}+B\times\left(\frac{\sigma}{r}\right)^{6}\right)</math>