Диссертация (Анализ данных атомно-силовой микроскопии с помощью компьютерного моделирования), страница 3

Для более полного описания поведения системыв конкретных случаях вводятся дополнительные слагаемые [2]. Например,для описания белков следует ввести водородные связи:U hb ij  i , j  rij1210  ij      ,r   ij  (1.9)16авомногихсовременныхпакетахмолекулярнойдинамикиучитывается второй порядок теории возмущений: возмущения параметроввзаимодействийсвязанныхчастиц-углов,длинсвязей;вводятсямежмолекулярные взаимодействия.Каждое слагаемое полной потенциальной энергии задается наборомпараметров, существуют различные наборы для разных типов модели. Ихзначения определяются из учета различных типов экспериментальныхданных(спектральные,калориметрические,кристаллографические)иквантово-механических расчетов. Значения всех параметров для каждоготипа частиц системы, а также набор формул для ее потенциала, потенциала иэнергии всей системы, и некоторые другие параметры – все это сохраняется вконфигурационный файл, который собственно и представляет собой силовоеполе.

Далее созданное силовое поле тестируется на малых молекулах,которые хорошо изучены экспериментально. Примеры силовых полей –OPLS[3], MMFF94[4], AMBER99[5], UFF[6], на их основе созданы пакетыпрограммного обеспечения GROMACS[7], TINKER[8], NAMD[1] и прочие.1.2Стадии компьютерного моделирования биополимеровметодом молекулярной динамики.Как правило, моделирование биомакромолекул, в частности, белков,предполагает следующие стадии:1.Подготовка системы. Помещение объекта моделирования всреду моделирования, соответствующую реальному эксперименту.2.Минимизация энергии системы для достижения наиболееэнергетически выгодной конформации обюъекта.3.Уравновешивание системы с ограничением движенияатомов биополимера в ансамбле с постоянными числом частиц,объемом и температурой (NVT ансамбль)174.Уравновешивание системы с ограничением движенияатомов биополимера в ансамбле с постоянными числом частиц,давлением и температурой (NPT ансамбль)5.Собственно моделирование6.Анализ траекторииОстановимся подробнее на отдельных стадиях процессамоделирования.Подготовка системыДля проведения молекулярной динамики биополимера необходимоиметь координаты и скорости всех его атомов, а также выбрать силовое полеи построить топологию всех элементов системы.Большинство биополимеров имеют трехмерные структуры,полученные либо из экспериментов, либо c помощью вычислений in silico.Существует внушительное разнообразие программ, позволяющихконструировать новые молекулы, основываясь на квантово-механическихрасчетах.

а также справочная литература по длинам связей, зарядам и углампрактически всех возможных комбинаций атомов в органической химии.Для начала, исследователю требуется тем или иным способом получитьтрехмерную структуру объекта изучения с атомистической точностью, ипредставить ее в виде координатного файла, в котором каждому атомуполагается идентификационный номер, номер основания, т.е. группы атомов,к которому он относится, три координаты в нанометрах либо ангстремах, атакже название, соответствующее данному атому в списке атомоввыбранного силового поля. Далее из такого координатного файла с помощьюпараметров, взятых из силового поля, будет создана топология объекта:файл, содержащий данные о зарядах, длинах связей, углах и двугранныхуглах всех атомов системы. Эта информация не будет меняться в процессемоделирования.18Силовое поле в программных пакетах молекулярной динамикипредставляет собой не только набор общих уравнений для потенциальнойэнергии, но также набор коэффициентов для каждого атома, каждой связи,каждого угла, доступного для моделирования в данном силовом поле.Возможности каждого силового поля напрямую зависят от этихкоэффициентов и количества рассмотренных атомов и их комбинаций.Наиболее распространенная проблема, которая может возникнуть на данномэтапе - это невозможность охватить все химические соединения в одномсиловом поле.

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

Такой метод требуетаккуратности, так как никаких проверок на адекватность системы непроводится.Отдельно следует упомянуть добавление растворителя. Как правило,моделирование проводится в жидкой среде. Растворитель может быть19дискретным, представляющим собой отдельные молекулы, либонепрерывным, который моделируется добавкой к потенциальной энергиисистемы нового члена, отвечающего за энергию сольватации. Как правило, втаком случае используется математическое представление Борна[9],свободная энергия сольватации в таком случае имеет три слагаемых: член,отвечающий за взаимодействие внутри растворителя в заполняемой полости,член, отвечающий за Ван-дер-Ваальсово взаимодействие между раствором ирастворителем, а также член, отвечающий за электростатическуюполяризацию раствора растворителем.Минимизация энергии.Функцияпотенциальнойэнергиимакромолекулярнойсистемыпредставляет собой довольно сложную многомерную гиперповерхность.

Онаимеет одну низшую точку, глобальный минимум, и великое множестволокальных минимумов, в которых все производные потенциальной энергииравны нулю, а все вторые производные не отрицательны[10].Знание всех локальных и глобального минимумов позволило быописать все возможные конформации макромолекулы, а также ее динамику иструктурные превращения. Однако огромное количество координат делаетневозможным описание этой гиперповерхности в сколько-нибудь значимойобласти точек.

На данный момент не существует методов энергетическойминимизации,которыебыгарантировалинахождениеглобальногоминимума. При известной начальной конфигурации системы можно найтитолько ближайший локальный минимум. "Ближайший" не в геометрическомсмысле, а такой, который можно достичь систематическим снижениемлокального градиента.На данный момент существуют три класса методов энергетическойминимизации.201.Методы, требующие анализа только самой энергетическойфункции.

Например, симплексный метод. Каждый последующий шагделается на основе анализа предыдущих. Если требуется информацияо производных энергетической функции, эти методы уступают тем,что используют эту информацию.2.Методы, требующие данные о производных энергетическойфункции. Частные производныепотенциальной энергии, по сутиявляющиеся силами, взятыми с обратным знаком, известны вмолекулярной динамике, поэтому к классу таких методов относятсямодификации алгоритмов молекулярной динамики.3.Методы, использующие также информацию о вторыхпроизводных энергетической функции. Они лучшие по своимсвойствам сходимости близ энергетического минимума: квадратичнаяпотенциальная функция минимизируется в один шаг. Однакопроблема заключается в том, что для молекулы из N частиц требуетсяпосчитать, запомнить и инвертировать матрицу из 3N×3N единиц.Даже не беря в расчет дополнительные программыные ресурсы длявысчитывания вторых производных, такие вычисления доступны длянемногих систем.Ко второй категории относятся наиболее популярные методыэнергетической минимизации, использующиеся, например, в программномпакете Gromacs.

К ним относится, например, метод наискорейшегоспуска[11]. В этом случае следующий шаг минимизации предпринимается внаправлении отрицательного градиента производной (т.е. в направлениисилы), без учета того, что было на предыдущих шагах. Размер шагаподбирается таким образом, чтобы обеспечивать при всей быстротепостоянное уменьшение энергии. Это простой и быстрый, но несколькопримитивный метод, его сходимость довольно низка, особенно вблизилокального минимума. Метод сопряженных градиентов демонстрирует более21быструю сходимость, он использует информацию о прошлых шагах[12].

Этотметод более медленный на ранних стадиях минимизации, чем методнаискорейшего спуска, однако он становится более эффективным по мереприближения к энергетическому минимуму.Алгоритм L-BFGS, или Бройдена — Флетчера — Гольдфарба —Шанно - один из наиболее применимых квазиньютоновских методов[13]. Оноснован на вычислении наилучшего приближения гессиана функции, исходяиз всех предыдущих шагов. Данный метод находит минимум любой дваждынепрерывнодифференциируемойдемонстрируетнаилучшуювыпуклойсходимость,функции.однакоиз-заЭтотметоднеобходимостипроводить коррекцию после каждого шага, не может быть распараллелен ине подходит для больших систем.Уравновешивание системы.Основнаязадачауравновешивания-привестивравновесиерастворитель вокруг объекта, то есть создать нужную температуру иплотность для молекул растворителя.

Приведение системы в состояниеравновесия производится в два шага: на первом, называемом NVTуравновешивание, приводится в равновесие температура системы, затем, спомощьюприложениядавлениявпроцессеNPTуравновешиваниядостигается необходимая плотность.Для того, чтобы объект исследования не участвовал в процессенагревания и перемещения растворителя, на тяжелые атомы объектанакладываются ограничивающие силы.