Моделирование перехода А-формы ДНК в В-форму в воде

1. Параметры моделирования

  1. Силовое поле используемое при построении топологии -- amber99sb.
  2. Заряд системы -- -10 (из-за заряженных фосфатов остова).
  3. Размер и форма ячейки -- "кубическая", 5.09200 х 4.90600 х 5.46200 нм.
  4. Минимизация энергии:
    1. Алгоритм минимизации энергии -- l-bfgs (quasi-Newtonian algorithm for energy minimization according to the low-memory Broyden-Fletcher-Goldfarb-Shanno approach).
    2. Алгоритм расчёта электростатики -- Cut-off.
    3. Алгоритм расчёта Ван-дер-Ваальсовых взаимодействий -- Cut-off.
    4. Модель, которой описывался растворитель -- tip3p.
  5. Утряска растворителя:
    1. Параметр который обуславливает неподвижность биополимера -- define = -DPOSRES.
    2. Число шагов -- 1000.
    3. Длина шага -- 0.0002 пс.
    4. Алгоритм расчёта электростатики -- PME.
    5. Ван-дер-Ваальсовых взаимодействий -- Cut-off.
    6. Алгоритм термостата -- V-rescale.
    7. Алгоритм баростата -- нету баростата.
  6. Основной расчёт МД:
    1. Время моделирования -- 6h13:53.
    2. Количество процессоров -- 8 OpenMP Threads.
    3. Эффективность масштабирования -- Wait + Comm. 1%.
    4. Длина траектории -- 10000 пс.
    5. Число шагов -- 5000000.
    6. Длина шага -- 0.002 пс.
    7. Алгоритм интегратора -- md (leap-frog algorithm for integrating Newton’s equations of motion).
    8. Алгоритм расчёта электростатики -- PME.
    9. Алгоритм расчёта Ван-дер-Ваальсовых взаимодействий -- Cut-off.
    10. Алгоритм термостата -- V-rescale.
    11. Алгоритм баростата -- Berendsen.

2. Моделирование

1) Минимизация энергии

При минимизации силы и потенциальная энергия изменились следующим образом:

Step 0, Epot=-2.693403e+03, Fnorm=1.938e+03, Fmax=4.684e+03 (atom 12)
Step 93, Epot=-6.034576e+03, Fnorm=1.004e+02, Fmax=7.389e+02 (atom 176)
In [1]:
from IPython.display import Image
display(Image('em.png'))

Рис. 1. Сравнение структуры ДНК после минимизации энергии (цианом покрашена оригинальная структура, магентой -- после минимизации).

Как видно из Рис. 1, ДНК почти не подвинулась, что не удивительно -- обычно она сворачивается в A-форму при пониженных концентрациях воды.

2) Утряска воды

In [7]:
display(Image('eq.png'))

Рис. 2. Сравнение структуры ДНК после утряски воды (цианом покрашена структура после минимизации энергии, магентой -- после утряски воды).

Из Рис. 2 сложно что-либо понять, но вполне очевидно, что вода больше всего подвинулась вблизи ДНК, что повлекло за собой небольшое движение даже на окраинах ячейки.

3) Молекулярная динамика

Делаем гифку:

In [4]:
import imageio
with imageio.get_writer(f'movie.gif', mode='I', fps=15) as writer:
        for x in range(0, 50):
            filename = f'movie/{x+1:04}.png'
            image = imageio.imread(filename)
            writer.append_data(image)
In [5]:
display(Image('movie.gif', format='png'))

Упс... fit прописали, а -pbc nojump -- забыли. Из-за этого очень сложно визуально определить точку перехода из одной формы в другую. Но ничего, анализ количественных признаков должен помочь нам решить эту проблему. А пока посмотрим на первый и последний фреймы:

In [6]:
display(Image('md1.png'))
display(Image('md2.png'))

Первый фрейм -- безусловно А-форма, -- а последний очень напоминает B-форму. Вроде все замечательно. Посмотрим теперь на количественные признаки.

3. Анализ результатов моделирования

0) Импорты

In [20]:
import pandas as pd
import seaborn as sns
from  jupyterthemes import jtplot
jtplot.style(figsize=(14, 9))

1) RMSD

По сравнению с изначальной структурой:

In [28]:
df = pd.read_csv('rms_1.tsv', sep='\t')
sns.lineplot(x='Time (ps)', y='RMSD (nm)', data=df)
_ = 0

Хммммм... Опять факап -- все тот же pbc нас губит. Тем не менее, видно, что RMSD подскакивает до 3 ангстрем практически сразу и далее держится примерно на этом уровне. Теперь посмотрим на RMSD по сравнению со структурой 400 кадров (4000 пс) назад:

In [29]:
df = pd.read_csv('rms_2.tsv', sep='\t')
sns.lineplot(x='Time (ps)', y='RMSD (nm)', data=df)
_ = 0

Мда, получилось не очень информативно, поскольку 4000 пс -- это почти половина траектории, а у нас ДНК часто скакала между ячейками.

2) Гидрофобная и гидрофильная поверхность

In [32]:
df = pd.read_csv('sas_dna.tsv', sep='\t')
sns.lineplot(x='Time (ps)', y=r'Area (nm\S2\N)', hue='type', data=df)
_ = 0

Фух, хотя бы тут все норм. Можно отметить тренд увеличения гидрофильной поверхности и уменьшения гидрофобной. Также бросаются в глаза гидрофильный пик при 8000 пс и гидрофобная яма при 9500 пс.

3) Водородные связи

Для начала, посмотрим на водородные связи внутри ДНК:

In [35]:
df = pd.read_csv('hbond_dna.tsv', sep='\t')
sns.lineplot(x='Time (ps)', y='Number of bonds', data=df)
_ = 0

Тут хорошо видно, что при 7000 пс начинаются какие то сильные изменения в структуре, что приводит к уменьшению количества водородных связей внутри ДНК. Можно заметить и второй скачок где-то около 9500 пс. Теперь посмотрим на водородные связи с водой:

In [36]:
df = pd.read_csv('hbond_dna_sol.tsv', sep='\t')
sns.lineplot(x='Time (ps)', y='Number of bonds', data=df)
_ = 0

Тут совсем ничего интересного.

4. Вывод

Много данных было потеряно в результате форс мажора с pbc, но основываясь на графиках поверхностей и водородных связей в ДНК, можно утвержать, что конформационный переход начался около 7000-8000 пс, но видимо он так и не закончился.