Энергия связей аденин-урацил в РНК

Постановка задачи

Вот так выглядит данная система

А вот так ее орбитали

Возьмем РНК из двух оснований: w:Аденин - система из 33 атомов и w:Урацил - система из 30 атомов. Они могут образовать различные варианты (модели) связей или просто расположений атомов. Программы молекулярного моделирования дают 3 цифры: энергию аденина, энергию урацила и их общую энергию. Хотелось бы понимать принципы расчета этой энергии. Т.е. имеем задачу: даны все расстояния между атомами, нужно рассчитать энергию систем аденина, урацила и общую.

Подходы к решению

В статье Анализ статьи: Автоматизированное предсказание de novo нативной РНК третичной структуры описывается общая функция оценки энергии для расчета. Для упрощения задачи можно учитывать всего две составляющие - это радиус Ван-дер-Ваальса и радиус Гирации. В более сложных задачах участвуют порядка 10 составляющих.

Оценка радиуса Гирации

Для расчета радиуса Гирации нужна матрица попарных расстояний между центроидами нуклеотидов. В данном случае, соответствующая матрица будет, например, следующей:

$cd={\begin{pmatrix}0&x\\x&0\end{pmatrix}}$

где х - расстояние между центроидами аденина и урацила. Например, х = 43.99 Тогда радиус Гирации:

$rg={\sqrt {\frac {cd_{sum}}{cd_{count}}}}={\sqrt {\frac {43.99}{2}}}=4.69$

$cd_{count}=n^{2}-n$

Оценка радиусов Ван-дер-Ваальса

В расчете радиусов Ван-дер-Ваальса участвуют все атомы системы, но их можно так же представить в виде матрицы пар нуклеотидов, в виде

$vdw={\begin{pmatrix}0&a_{1}\\a_{1}&0\end{pmatrix}}$

Например, если $a_{1}=17.9$, то

$vdw=\sum a*0.8=17.9*0.8=14.3$

Дальше нам будет важна структура каждого нуклеотида. В расчетах компонент vdw (в данном случае только $a_{1}$) участвуют только 9 атомов в зависимости от типа нуклеотида (в таблице эти атомы отмечены желтым цветом).

№	rA - 27	rC - 28	rU - 29	rG - 26
1	P - 1	P - 1	P - 1	P - 1
2	O1P	O1P	O1P	O1P
3	O2P	O2P	O2P	O2P
4	O5*	O5*	O5*	O5*
5	C5* - 2	C5* - 2	C5* - 2	C5* - 2
6	C4*	C4*	C4*	C4*
7	O4*	O4*	O4*	O4*
8	C3* - 4	C3* - 4	C3* - 4	C3* - 4
9	O3*	O3*	O3*	O3*
10	C2*	C2*	C2* - 9	C2*
11	O2* - 9	O2* - 9	O2*	O2* - 9
12	C1* - 3	C1* - 3	C1* - 3	C1* - 3
13	N1 - 8	N1 - 8	N1 - 8	N1 - 8
14	C2	C2	C2	C2
15	N3	O2 - 6	O2 - 6	N2 - 6
16	C4 - 6	N3	N3	N3
17	C5	C4	C4	C4
18	C6	N4 - 5	O4 - 5	C5
19	N6 - 5	C5	C5	C6
20	N7	C6 - 7	C6 - 7	O6 - 5
21	C8 - 7	1H5*	1H5*	N7 - 7
22	N9	2H5*	2H5*	C8
23	1H5*	H4*	H4*	N9
24	2H5*	H3*	H3*	1H5*
25	H4*	1H2*	1H2*	2H5*
26	H3*	2HO*	2HO*	H4*
27	1H2*	H1*	H1*	H3*
28	2HO*	1H4	H3	1H2*
29	H1*	2H4	H5	H1*
30	H2	H5	H6	2HO*
31	1H6	H6		H1
32	2H6			1H2
33	H8			2H2
				H8

В файле Rna vdw parameter.dat содержаться квадраты Ван-дер-Вальсовых радиусов для различных комбинаций атомов и нуклеотидов РНК. В этом файле первых две цифры - это коды соответствующих атомов, следующие две - это коды нуклеотидов, и пятая квадрат радиуса Ван-дер-Вальсса - atomvdw.

Чтобы рассчитать компонент $a_{1}$ нужно проверить все 81 комбинации, и если атомы находятся на расстоянии меньшем, чем Ван-дер-Вальсовый радиус сложить соответствующие оценки. Это можно сделать по следующему алгоритму:

dist2 = atomvdw - ($(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}+(z_{1}-z_{2})^{2}$);
if ( dist2 > 0.0 ) 
{
    vdw_score += ( dist2 * dist2 ) / atomvdw;
}

Дополнительно нужно также учесть возможное перекрытие центроидов:

atomvdw = 3.0 * 3.0;
dist2 = atomvdw - cendist_ij;
if ( dist2 > 0.0 ) vdw_score += ( dist2 * dist2 ) / atomvdw;