Моделирование больших биомолекул и биомолекулярных систем с использованием графического процессора

Белковые волокна, например, фибронектин, волокна фибрина, микротрубочки и актиновые филаменты, выполняют важные механические функции при формировании цитоскелета и поддержании работоспособности клетки [1−3], в процессе слияния клеток и формировании внеклеточной матрицы [4−7], а также при свертывании крови [8−10]. Физические свойства капсул растительных и животных вирусов [11—13], ретровирусов [14] и бактериофагов [15, 16], а также переходы между их стабильным и нестабильным состояниями определяют жизненный цикл многих вирусов, в том числе созревание вируса и заражение клеток [17]. Изучение происхождения уникальных упругоэластичных свойств белковых волокон и механизмов перехода от эластичного к пластичному состоянию в капсулах, а также возможность контролировать их динамическое поведение в ответ на механическое воздействие являются важными сферами исследования в биофизике. Современные технологии, предназначенные для изучения одиночных молекул, такие как атомно-силовая микроскопия и оптические пинцеты, широко используются для экспериментального изучения механических свойств белковых волокон [18−21] и капсул вирусов [15, 16, 22, 23]. Однако ввиду сложности строения этих систем 103 — 105 частиц) и их больших размеров.

50 — 200 нм), результаты подобных экспериментов почти невозможно интерпретировать без предварительных знаний о ландшафте свободной энергии [9].

Стандартные вычислительные пакеты Молекулярной Динамики (МД) в полноатомном разрешении, такие как CHARMM [24, 25], NAMD [26−28], и Gromacs [29— 31], широко используются для изучения поведения биомолекул на субмолекулярном уровне. Так как полно-атомное моделирование в настоящее время ограничено размером молекулы в 10−50 нм и длительностью процесса моделирования в 0,1 — 10 мс [32, 33], данный подход хорошо подходит только для моделирования равновесных процессов, а достижение биологически важного временного интервала от микросекунды до секунды практически невозможно даже для малых систем. Что ещё более важно, для подробного изучения ландшафта свободной энергии, лежащего в основе изучаемого биологического процесса, требуется статистически значимое количество траекторий. Одно из возможных решений такой задачи — проведение МД моде л ирования на компьютерных кластерах — требует огромных вычислительных ресурсов и длительного времени выполнения программы. К примеру, для расчёта 20 коротких (1 не) траекторий для южного вируса мозаики бобовых (southern bean mosaic virus), состоящего из более чем 4,5×106 атомов, потребовалось 800 ООО процессоро-часов работы кластера SGI Altix 4700 [34]. Это ограничивает возможность применение вычислительного эксперимента для изучения широкого спектра биологических проблем, таких как деформация волокон белков, формирование биомолекулярных комплексов и агрегатов, механическое повреждение капсул вирусов, для которых экспериментальные данные уже получены, а прямое сопоставление результатов экспериментов и вычислительных расчётов невозможно.

Основные подходы, позволяющие напрямую сопоставить экспериментальные данные с результатами численного моделирования биомолекул включают (1) использование неявного растворителя, (2) упрощение математической модели и (3) использование высокопроизводительных вычислительных систем. Методы моделирования на основе неявного растворителя предполагают, что наличие естественной среды (воды) можно описать при помощи специальной эмпирической функции [35, 36]. Эта функция описывает свободную энергию взаимодействия белка с внешней средой, а механическое влияние столкновений молекул воды с биомолекулой описывается при помощи уравнений Ланжевена. В основе упрощённого моделирования лежит предположение, что не все степени свободы биомолекулы одинаково важны. Упрощение производится путём сокращения моделируемых степеней свободы, и за структурную единицу системы берётся не атомы системы, а их группа (например, одна амино-кислота). Такой метод не только позволяет существенно сократить количество производимых вычислений, но и перейти к большему шагу по времени без потери численной стабильности системы. В число высокопроизводительных систем входят вычислительные кластеры, облачные вычисления, а также узкоспециализированные вычислительные устройства и графические процессоры. Вычислительные кластеры уже давно применяются для биомолекулярного моделирования, однако, их покупка и обслуживание зачастую недоступны для средней научно-исследовательской группы. Облачные вычисления требуют разработки и поддержки сложной инфраструктуры для проведения рассчётов и сбора данных [37, 38]. Разработка узкоспециализированных аппаратных средств невозможна без привлечения высококвалифицированных инженеров [39−41]. Графические Процессоры (ГП), которые изначально были спроектированы для ускорения работы с трёхмерной графикой, способны выполнять многие вычислительные задачи, в том числе и те, которые не связанны с обработкой изображения [42−44]. Недавние технологические достижения на аппаратном уровне, поддержка стандарта IEEE для вещественной арифметики позволяет использовать огромные вычислительные возможности ГП в научных приложениях. В отличие от процессоров с привычной архитектурой, большинство логических элементов ГП отведено на выполнение вычислений, а не на кеш-память и управление логикой. Массивная многопоточность, минимальный контекст потоков и высокая пропускная способность памяти делают ГП эффективным массивно-параллельным вычислительным устройством. Программные платформы для современных ГП включают NVIDIA Compute Unified Device Architecture (CUDA) [45, 46] и Open Computing Language (OpenCL) [47]. CUDA, программная среда для параллельных вычислений, является высокоуровневой программной платформой, расширяющей стандартные языки С и С++. Это позволяет разработчику реализовывать процедуры (ядра), которые могут выполняться одновременно во множестве независимых потоков на ГП. Ещё одним преимуществом ГП является их доступность: цена даже самых современных устройств не превышает 500 — 2000 долларов США, а обычный персональный компьютер, оборудованный двумя-четырьмя такими устройствами, способен достигать производительности компьютерного кластера с 200 — 500 вычислительными ядрами.

Из-за фундаментальных различий архитектуры Графических Процессоров (ГП) и Центральных Процессоров (ЦП), методы молекулярного моделирования, разработанные для выполнения на ЦП, не могут быть просто перенесены или адаптированы для работы на ГП. Тем не менее, в молекулярной динамике парные взаимодействия обычно описываются одной и той же эмпирической функцией потенциальной энергии для всех пар взаимодействующих частиц (силовое поле), а динамика системы определяется из численного решения одного и того же уравнения движения для всех частиц. Таким образом, существует прямое соответствие между ОКМД (Одиночный поток Команд, Множественный поток Данныхангл. Single Instruction, Multiple Data, SIMD) архитектурой ГП на аппаратном уровне и вычислительными процедурами молекулярной динамики на программном уровне. Можно выполнить одну и ту же процедуру (вычисление потенциальной энергии или сил, генерация случайных чисел и численное интегрирование уравнений движения) одновременно для различных наборов данных (для всех частиц) за счёт использования множества арифметических логических устройств, работающих параллельно, повторяя вычисления на протяжении множества шагов по времени. Поэтому, молекулярные расчеты — естественный кандидат для реализации на ГП, но, для эффективного выполнения алгоритма на ГП, он должен быть преобразован для запуска множества независимых потоков, выполняющих одинаковый программный код на различных наборах данных одновременно. Действительно, на данный момент, существуют предварительные версии стандартных пакетов для моделирования МД белков, реализованные на ГП, такие как NAMD [48−51], Gromacs [52], и др. [53−56].

В данной работе представлена комбинация двух методов ускорения биомолекулярных рассчётов: упрощения вычислительной модели и использования графических процессоров в качестве вычислительной платформы. Так как молекулярные преобразования в белках, подверженных внешнему механическому воздействию, определяются главным образом топологией и общей структурой системы, для описания белков была использована упрощённая [57—59] Модель Самоорганизующегося Полимера (англ. Self Organized Polymer, SOP) [60, 61]. В методологической части работы описаны численные методы, использованные для расчёта потенциалов взаимодействия модели SOP, динамики Ланжевена и генерации случайных чисел полностью реализованные на ГП. Методы были тщательно исследованы на наличие вычислительных и аппаратных ошибок, проверена возможность использования чисел с плавающей точкой одинарной точности, проведено тестирование на случайность для генераторов псевдо-случайных чисел. Далее разработанная методология применена для моделирования силовой денатурации фибриногена и его фрагментов, а также капсулы вируса H К 97. Результаты были сравнены с доступными экспериментальными данными и не только позволили их по-новому интерпретировать, но также предоставили объяснение происходящих процессов на субмолекулярном уровне.

Целью данной работы является разработка программной реализации численных методов упрощённого моделирования белковых молекул с использованием графических процессоров. Такая комбинация позволяет моделировать большие белковые системы на экспериментальных временных интервалах без использования дорогих в покупке и обслуживании вычислительных кластеров. Так как основной целью являлось моделирование экспериментов на единичной молекуле, таких как Атомная Силовая Микроскопия, для реализации была выбрана модель SOP. Эта модель обладает простой потенциальной функцией, а её реализация на ЦП уже позволяла моделировать небольшие системы 50−300 амино-кислот) на экспериментальных временных интервалах. Таким образом, использование модели SOP полностью реа-лизованноё на ГП дало возможность наблюдать микромолекулярную динамику белковых систем размером до 106 аминокислот.

Предметом исследований являются мономер и димер фибриногена, его фрагменты, а также капсид вируса НК97. Экспериментальные данные, полученные при помощи Атомной Силовой Микроскопии, показали, что молекулы фибриногена обладают сложными механическими свойствами. Однако, полученные данные не позволили точно определить источник динамических характеристик этой молекулы. Именно поэтому, стало необходимым произвести молекулярное моделирование силовой денатурации фибриногена. Однако, ни один из существующих вычислительных методов не позволял произвести молекулярное моделирование в условиях воздействия на молекулу, идентичных экспериментальным. Поэтому было решено адаптировать вычислительные методы для эффективной работы на ГП. Полученная реализация также была применена для объяснения динамических характеристик капсида вируса НК97.

Таким образом, данная тема является актуальной с точки зрения развития методологии молекулярного моделирования, поскольку представленная программная реализация является единственной существующей, позволяющей моделировать большие белковые системы на экспериментальных временных интервалах. Результаты численных исследований, полученные для мономера, димера фибриногена и его фрагментов, позволили по-новому интерпретировать экспериментальные данные и лучше понять микромеханику данного белка. Моделирование капсида НК97 показало, что динамические характеристики данной белковой системы зависят от скорости и геометрии воздействия. Полученная программная реализация была организована как программный пакет БОР-СРи, который распространяется в виде открытого кода, сопровождается детальным руководством пользователя и не требует навыков программирования для своего использования. Таким образом, он может быть использован другими научными группами для исследования интересующих их биомолекулярных систем.

6.3 Выводы.

Реализация упрощённой модели SOP на ГП позволила-детально изучить микромеханику капсида вируса НК97. В процессе моделирования был использован как протокол, достаточно близкий к АСМэкспериментам, так и тот, в котором внешнее воздействие прилагалось быстрее. Полученные результаты показали, что ответ вируса на внешнее воздействие сильно зависит и от скорости движения зонда ^ и от геометрии внешнего воздействия (радиуса зонда R).

В процессе моделирования-наблюдался спектр всевозможных реакций капсида на внешнее воздействие — от быстрого вминания и равномерного продавливания при малых силах до механического разрыва при больших. Теоретически полученные динамические параметры могут характеризовать физические свойства, общие для целого ряда оболочек вирусов. Было обнаружено, что динамический ответ вируса на внешнее механическое воздействие является стохастической динамической характеристикой, которая также зависит от внешнего воздействия. Полученное при близкой к эксперименту скорости движения зонда v ~ 2,5мкм/с значение коэффициента упругости К ~ 0.01 — 0.02Н/м оболочки вируса НК97 хорошо согласуется с экспериментальными значениями для пустых оболочек вирусов [222]. Коэффициент упругости К является локальной характеристикой, так как он зависит от размера зонда. Большие зонды воздействуют на большее количество структурных единиц оболочки вируса, которые взаимодействуют между собой, сопротивляясь приложеной силе. Поэтому использование больших зондов, размеры которых сопоставимы с размерами капсида, позволяет получать усреднённые характеристики. В дополнение, при малых размерах зонда, важно знать точку его соприкосновения с капсидом. Исследование этого вопроса осталось за рамками данной работы.

Было обнаружено, что потеря сопротивления на внешнее воздействие со стороны капсида (К) может быть обусловлена как проминанием последнего, так и его механическим разрывом. В случае проминания капсид довольно быстро восстанавливает свои упругие свойства, что характеризуется значительным возрастанием К. Быстрое падение значения коэффициента упругости К скорее всего является показателем разрыва капсида, когда нативные контакты частично разрываются, а небольшие части оболочки 3 — 5нм) капсида денатурируют. Это обратимый процесс, так как часть нативных контактов восстанавливается почти сразу после прохождения зондом его оболочки (Рис. 22). Это хорошо согласуется с экспериментальными результатами, полученными на других капсидах [15, 222]. Также нами было обнаружен ожидаемый переход от упругого поведения при небольших силах внешнего воздействия к вязкому режиму при больших силах [34], что также отражается в уменьшении К с увеличением силы. Этот эффект не так хорошо заметен, что может быть объяснено присутствием дополнительных ковалентных связей, стабилизирующих оболочку [224]. Эта структурная особенность капсида НК97 усиливает её упругие свойства. Даже сильное падение количества нативных контактов (<3) при механическом разрыве капсида обуславливается в основном разрывами контактов внутри белковых структурных единиц, а не между ними.

Было показано, что модуль Юнга У зависит от скорости движения зонда и от его размеров (Таблица 7). При скорости, максимально близкой к экспериментальной (г- = 2,5мкм/с), модуль Юнга для капсида НК97 равен У — 60 — 160МПа, что хорошо соотносится с экспериментальным значением, полученным для пустого капсида ССМУ, где У = 140МПа [222], но несколько больше значения, полученного для капсида 029, где У = 1,8ГПа [15]. Пустая оболочка вируса НК97 способна выдерживать давление в 60 — 140МПа, что хорошо согласуется со значением давления упакованного ДНК на стенки капсида 029. Эти результаты показывают, что численное моделирование индентации под действием силы может быть использовано для того, чтобы предсказать пределы упругости капсид вирусов, а также чтобы оценить максимальное внутреннее давление, которое они способны выдерживать.

Заключение

.

Механические свойства больших белков, белковых комплексов и образований представляет большой интерес в современной биофизике. Размер и составная структура этих биомолекулярных систем усложняют интерпретацию доступных результатов экспериментальных исследований, таких как эксперименты на единичных молекулах с применением Атомной Силовой Микроскопии (АСМ) и оптических пинцетов. Методы молекулярного моделирования, совместно с экспериментальными методами способны предоставить детальную картину происходящих микромолекулярных процессов. В силу ограничений АСМ оборудования, экспериментальные результаты легче получать для больших систем, а скорость применяемого механического воздействия ограничена сверху и колеблется в районе 0,1−10, Омкм/с. В молекулярном моделировании, наоборот, размер системы негативно влияет на производительность программы, а скорость движения зонда ограничена снизу. В полноатомном разрешении, где шаг интегрирования по времени составляет 1 — 2фс, получение результатов на экспериментальном или биологическом временном интервале невозможно даже с применением новейшего вычислительного оборудования. Упрощённые модели используют шаг интегрирования в 1−100пс, поэтому на данный момент это единственный способ получить результаты молекулярного-моделирования в условиях, ин-дентичных экспериментальным. Но для больших белковых систем, таких как нити фибрина и капсиды вируса, даже упрощённые модели требуют серьёзных вычислительных затрат.

Графические процессоры, которые изначально были спроектированы для ускорения работы с трёхмерной графикой, обладают беспрецендентными вычислительными возможностями. Пиковая вычислительная производительность современных ГП пре.

I • восходит ЦП аналогичного уровня технологий на два порядка. С появлением высокоуровневых языков программирования на ГП, таких как CUDA и OpenCL, разработка программного обеспечения, полностью адаптированого под работу на ГП, больше не требует детального знания графического интерфейса разработки приложений. Но, в силу новизны данных устройств в качества вычислительной платформы, необходима адаптация программного обеспечения для работы на ГП.

В данной работе представлена реализация упрощённой модели БОР, полностью работающая на ГП (БОР-СРи). Рассмотрена реализация расчёта сил межчастичного взаимодействия, составления списков Верле, интегрирования уравнений движения. Детально описаны способы генерации случайных чисел, необходимых для численного решения уравнений Ланжевена. Также описан подход «много запусков на ГП», позволяющий получать статистически важную выборку результатов за короткое время. Программная реализация была тщательно проверена на предмет ошибок, как аппаратных, так и программных. Была отсеяна возможность возникновения корреляций в силу огромного количества случайных чисел требуемых для таких высокопроизводительных вычислений. Скорость работы реализации была сравнена с работой аналогичной программы на ЦП, а результирующее ускорение расчётов достигало 200раз для систем большого размера. Приведена также и скорость вычислений для систем разного размера — от домена IVIV (34 аминокислоты) до капсулы вируса НК97 (115 140 амино-кислот).

Полученная реализация была использована для моделирования механических свойств фибриногена и капсида вируса НК97. Результаты моделирования на мономере и димере фибриногена позволили по-новому интерпретировать результаты АСМ-экспериментов. Было показано, что типичная «пилообразная» форма графика зависимости силы от растяжения, полученная в экспериментах по механической денатурации мономеров и олигомеров фибриногена, возникает из-за последовательной денатурации 7-модулей белка. При этом, каждый 7-модуль денатурирует в три этапа, а растяжение и сжатие свёрнутых «-спиралей упруго, без существенных падений силы. Численные исследования различных фрагментов фибриногена показали важность белок-белковых взаимодействий между различными доменами молекуля. В частности, была охарактеризована важность дополнительной стабилизации 7-модуля его нековалентными взаимодействиями с /3-модулем.

Численные исследования капсида вируса НК97 показали все три возможных сценария механической индентации — быстрый прогиб, упругий режим и механический разрыв оболочки. Оказалось, что вероятность возникновения того или иного режима зависит от скорости внешнего воздействия (скорости движения зонда) и геометрии радиуса) зонда. Использование большого зонда (Я = 25нм) приводило к быстрому прогибу или упругой реакции капсида, маленький зонд (Я = 5нм) прорывал оболочку капсида. Увеличение скорости движения зонда также приводило к большей вероятности механического разрушения оболочки вируса. В целом, полученные численные оценки оказались схожими к оценкам, полученным на других капсидах экспериментально.

Реализация программы БОР-СРи оформлена в виде готового программного продукта и размещена в открытом доступе в сети интернет. Для взаимодействия с программой используется текстовый файл параметров симуляции, а входные (выходные) данные принимаются (сохраняются) в общепринятых форматах, используемых другими программами для молекулярного моделирования [24, 26, 29]. Таким образом, для использования пакета ЭОР-СРи не требуется навыков программирования, а результаты моделирования могут быть показаны и обработаны при помощи популярных средств (например, программы У1УГО[225]).