Реализация модели самоорганизующегося полимера SOP на графическом процессоре

SOP — это приближённая модель структуры полипептидной цепи, описывающая изменение цепи во времени посредством динамики Ланжевена в трехмерном пространстве. В модели используется простое для вычисления силовое поле (force field), при этом оно достаточно точно для получения макроскопических физических свойств биомолекулярной системы. Предыдущие исследования показали, что SOP модель хорошо описывает механические свойства белков, например зелёного флуоресцентного белка (green fluorescient protein)[1], димера тубулина (tubulin dimer)[2] и кинезина (kinesin)[3]. SOP-модель также была применена для изучения кинетики и построения энергетического ландшафта свободной энергии различных биомолекул, таких как рибозим тетрахимена (tetrahymena ribozyme)[4], riboswitch aptamers [5], GroEL [6], дигидрофолат-редуктаза (DHFR)[7], Протеинкиназа А (protein kinase A)[8], Миозин V (myosin V)[9].
В SOP-моделе каждый аминокислотный остаток цепи представляется всего одной Cα-частицей, что сокращает общее число степеней свободы системы. Кроме того, симуляции динамики Ланжевена реализованы на графических процессорах (ГП), что дает вычислительное ускорение в ~10-200 раз, в зависимости от размера изучаемой системы. В комбинации эти свойства SOP модели позволяют наблюдать за динамикой системы на очень большом временном интервале — от миллисекунд до секунды. Этот диапазон соответствует экспериментальной временной шкале атомно-силовой микроскопии (АСМ), проводимой на одиночных белках, белковых тандемах, белок-белковых комплексах и агрегатах, белковом волокне и белковых образованиях. В настоящее время фактически остается невозможным достичь временную шкалу, соответствующую важным биологическим процессам (от миллисекунды до секунды) даже для очень маленькой системы, состоящей всего из нескольких десятков аминокислот, используя традиционный полноатомный метод (явного и неявного растворителя) молекулярной динамики (МД), реализованный на самых мощных компьютерных кластерах. В случае экспериментов in silico, параметры приложенной силы и геометрия кантилевра, используемые полноатомным МД приближением, отличаются от тех, что применяются в экспериментах in vitro. В этом смысле SOP-модель, разработанная в программе SOP-GPU — единственный известный подход, позволяющий исследовать механическо-химические свойства биологических систем

  • на экспериментальной временной шкале (до секунды).
  • используя реальные параметры приложенной силы.

Это достаточно важно ввиду того, что кинетика (траектории, промежуточные состояния), термодинамика (барьеры преходов и рабочие расстояния), и механизмы — всё это меняется в зависимости от скорости rf изменения силы fext = rft и от значения силы fext. Когда рентгеновская структура биологической системы разрешена (из Protein Data Bank или других источников), модель SOP-GPU позволяет исследователям изучать:

  • изменение структуры белков, РНК и ДНК под действием силы>
  • механизм диссоциации биомолекулярных комплексов и аггрегатов;
  • сжатие белковых волокон и капсул вирусов растений и животных;

Такого рода исследования могут быть проведены за приемлемое время с помощью программы SOP-GPU. Поэтому результаты экспериментов на одной молекуле можно непосредственно сравнивать с симуляциями, что помогает исследователям получать содержательную интерпретацию результатов экспериментальных измерений, например динамику изменения силы и графики зависимости силы от смещения кантилевра.

Модель SOP

В моделе SOP [1, 4, 10] каждый аминокислотный остаток описывается одним центром взаимодействия (Cα-атом). Динамика i-ой Cα-частицы получается из решения уравнений Ланжевена в задемпфированном пределе.

Где ξ — коэффициент трения, и f(Ri)=−∂V/∂Ri — молекулярная сила действующая на частицу вследствие изменения потенциальной энергии V. Здесь Gi(t) — случайная сила с гауссовым распределением, средним значением равным нулю и дельта-функцией корелляции (белый шум). Случайная сила учитывает стохастические столкновения аминокислотных остатков белка с молекулами растворителя (воды). Функция потенциальной энергии конформации белка V, выраженная через координаты Cα-частиц {R} = R1, R2, …, Rn, задаётся формулой:

где rri,i+1 — расстояние между двумя взаимодействующими остатками i и i+1, а r0i,i+1 — его значения в нативной (PDB) структуре. Первый член уравнения — действующий на ограниченном расстоянии нелинейный упругий (FENE) потенциал, описывающий связи между Cα в цепи. Здесь R0=2 Å — толерантность ковалентной связи (tolerance in the change of a covalent bond) и k=1.4 N/m — коэффициент упругости. Второй член уравнения — потенциал Леннарда-Джонса (VNBATT), который учитывает нативные взаимодействия, стабилизирующие нативное состояние. Когда в SOP моделе два несвязанных ковалентно остатка i и j (|i — j| > 2) находятся на расстояние менее RC в нативном состоянии, т.е. rij < RC, тогда Δij=1, иначе 0. Значение εh определяет силу взаимодействия несвязанных ковалентно аминокислотных остатков. Обычно значение εh находится в диапозоне 1,0-1,5 ккал/моль, и зависит от вторичной структуры рассматриваемого белка. Все ненативные взаимодействия (четвёртый член уравнения) считаются отталкивающими (VNBREP). На угол формируемый остатками i, i+1 и i+2 накладывается дополнительное ограничение при помощи отталкивающего потенциала с параметрами εl=1 ккал/моль и σ=3.8 Å, которые определяют соответственно силу и диапазон отталкивания. Чтобы удостовериться в отсутствии самопересечений цепи белка выбрано именно такое значение σ.

Методы симуляций

Поскольку графические процессоры (ГП) отличаются от центральных процессоров (ЦП) по многим фундаментальным параметрам, методы биомолекулярных симуляций для ЦП не могут быть просто перенесены на ГП. В молекулярных симуляциях взаимодействие между парами частиц описываются одной и той же эмпирической функцией потенциальной энергии (force field), а динамика системы получается из численного решения одинаковых для всех частиц системы уравнений движения. Для модели SOP мы разработали численные методы для генерации списков Верле и случайных сил, чтобы с их помощью для вычислять силы вызываемые ковалентными и нековалентными взаимодействиями, а также численного интегрирования уравнения движения Ланжевена [10]. Детальное описание численных алгоритмов можно найти в этой статье.
На графическом процессоре существуют две основные стратегии оптимизации парных взаимодействий. В первом случае, все силы, действующие на одну частицу вычисляются в одной ните, что требует запуска N нитей, для получения всех сил, действующих на частицу. Мы называем такой подход параллелизацией по частицам. При использование этого подхода требуется дважды вычислять одну и ту же составляющую силы df, действующую на i-ю и j-ю частицы в i-й и j-й нитях. Чтобы описать дальнодействующие силы мы используем общее предположение, что дальнодействующие силы исчезают за пределами некоторого определенного расстояния. Это позволяет нам построить и использовать списки Верле. При параллелизации по частицам списки Верле могут легко обновляться на ГП для ускорения вычисления потенциальной энергии. В стратегии параллелизация по парам взаимодействующих частиц, вычисления сил выполняются для всех пар (i, j), используя суммарно P вычислительных нитей, что соответствует количеству взаимодействующих пар. Далее 2P значений сил сохраняются в различные места глобальной памяти ГП. В этом подходе, генерация списка Верле означает формирование вектора всех пар взаимодействующих аминокислотных остатков для каждого члена потенциальной энергии. Чтобы получить упорядоченный список Верле, список должен быть упорядочен и обновлен на ЦП. Это более эффективно, вычислять расстояния между частицами на ГП, копировать их в память ЦП, и потом строить новый список Верле [10].

Уравнения движения Ланжевена для каждого аминокислотного остатка (Cα-частицы) численно интегрируются с использованием схемы первого порядка аппроксимации (алгоритм Ермака-Маккамона)[10],

где h — шаг интегрирования. Тестирования точности численного интегрирования, оценивающие величину ошибок вычислений, накопленную за многие миллионы шагов, показали, что в длительных симуляциях на ГП арифметика одинарной точности и алгоритм Ермака-Маккамона могут быть использованны для точного описания свойств биомолекул.
Для симуляций динамики Ланжевена требуется хороший генератор (псевдо-)случайх чисел, gi,α (α = x, y, и z, и = 1,2,…,n), для вычисления компонент гауссовой случайной силы:

где T — абсолютная температура (kb — постоянная Больцмана). Используя различные методы, генератор случайных чисел (RNG) выдаёт последовательность случайных чисел ui,α, равномерно расределённых на единичном интервале [0, 1], которая затем преобразуется в последовательность нормально распределённых чисел с нулевым математическим ожиданием и единичной дисперсией (gi,α). В программе SOP-GPU мы реализовали хорошо известное преобразование Бокса-Мюллера. Детальное описание вычислительных алгоритмов для реализации генераторов случайных чисел на ГП можно найти на другой странице.

Результаты тестов производительности пакета sop-gpu, реализованого на ГП можно найти в другом разделе сайта.

Литература

[1] Mickler M, Dima RI, Dietz H, Hyeon C, Thirumalai D, Rief M. Revealing the bifurcation in the unfolding pathways of GFP using single molecule experiments and simulations. Proc Natl Acad Sci USA 2007; 104: 20268–20273.

[2] Dima RI, Joshi H. Probing the origin of tubulin rigidity with molecular simulations. Proc Natl Acad Sci USA 2008.

[3] Hyeon C, Onuchic JN. Internal strain regulates the nucleotide binding site of the kinesin leading head. Proc Natl Acad Sci USA 2007; 104: 2175–2180.

[4] Hyeon C, Dima RI, Thirumalai D. Pathways and kinetic barriers in mechanical unfolding and refolding of RNA and proteins. Structure 2006; 14: 1633–1645.

[5] Lin JC, Thirumalai D. Relative stability of helices determines the folding landscape of adenine riboswitch aptamers. J Am Chem Soc 2008; 130: 14080–14084.

[6] Hyeon C, Lorimer GH, Thirumalai D. Dynamics of allosteric transitions in GroEL. Proc Natl Acad Sci USA 2006; 103: 18939–18944.

[7] Chen J, Dima RI, Thirumalai D. Allosteric communication in dihydrofolate reductase: signaling network and pathways for closed to occluded transition and back. J Mol Biol 2007; 374: 250–266.

[8] Hyeon C, Jennings PA, Adams JA, Onuchic JN. Ligand-induced global transitions in the catalytic domain of protein kinase A. Proc Natl Acad Sci USA 2009; 106: 3023–3028.

[9] Tehver R, Thirumalai D. Rigor to post-rigor transition in myosin V: Link between the dynamics and the supporting architecture. Structure 2010; 18: 471–481.

[10] Zhmurov A, Dima RI, Kholodov Y, Barsegov V. SOP-GPU: Accelerating biomolecular simulations in the centisecond timescale using graphics processors. Proteins 2010; 78: 2984–2999,

[11] A.А. Жмуров, В.А. Барсегов, С.В. Трифонов, Я.А. Холодов, А.С. Холодов. Моделирование микромеханики на графических процессорах с использованием динамики Ланжевена // Математическое моделирование. 2011, т. 23, № 10, с. 133-156.