Бесплатный автореферат и диссертация по биологии на тему

Атомные силовые поля FFSol и QMPFF3

ВАК РФ 03.01.02, Биофизика

Автореферат диссертации по теме "Атомные силовые поля FFSol и QMPFF3"

На правах рукописи

ПЕРЕЯСЛАВЕЦ ЛЕОНИД БОРИСОВИЧ

АТОМНЫЕ СИЛОВЫЕ ПОЛЯ РРвОЬ И ОМРРРЗ: СОЗДАНИЕ, ПАРАМЕТРИЗАЦИЯ И ТЕСТИРОВАНИЕ

03.01.02-биофизика

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

ПУЩИНО - 2010

? 0 "'" -Г.З

Работа выполнена на кафедре прикладной теоретической физики Национального исследовательского университета Московского физико-технического института и в Учреждении Российской академии наук Институте белка РАН

Научные руководители:

доктор физико-математических наук

Зосимов Виктор Васильевич

член-корр. РАН, доктор физико-математических наук

профессор Финкельштейн Алексей Витальевич

Официальные оппоненты:

доктор физико-математических наук Полозов Роберт Валентинович кандидат физико-математических наук Сивожелезов Виктор Семенович

Ведущая организация:

Биологический факультет Московского государственного университета им. М.В. Ломоносова

Защита состоится 26 мая 2010 г. в 15 часов 30 минут на заседании совета Д002.093.01 по защите докторских и кандидатских диссертаций при Учреждении Российской академии наук Институте теоретической и экспериментальной биофизики РАН по адресу: 142290, г. Пущине Московской обл., ул. Институтская, д. 3, ИТЭБ РАН.

С диссертаций можно ознакомиться в Центральной библиотеке Пущинского НЦБИ РАН по адресу: 142290, г. Пущино Московской обл., ул. Институтская, д. 3.

Автореферат разослан 26 апреля 2010 г. Ученый секретарь

диссертационного совета Д002.093.01 кандидат физико-математических наук

//ЛиН.Ф. Ланина

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность проблемы. Силовые поля (СП) представляют собой набор потенциалов и атомных параметров для описания внутри- и межмолекулярных взаимодействий. Силовые поля могуг использоваться для предсказания свойств новых веществ, для моделирования макромолекул (например, белков, РНК и ДНК), для оценки свободной энергии связывания молекул (белок-белок, белок-лиганд), и т.н. Последнее имеет важнейшее применение в медицине, так как предварительный теоретический отбор тысяч лучших лигандов из миллионов позволяет существенно сэкономить на поиске новых лекарств к целевому белку. Оценки констант (или свободной энергии) связывания белок-лиганд делятся по точности на два подхода. Первый подход - более точный, но медленный, - как правило, базируется на методе возмущений в оценке свободной энергии взаимодействий [Kollman, 1993, Chem. Rev., 93:2395]', он требует длительного времени вычислений и точного силового поля, с детальной атомной параметризацией, которая позволяет хорошо описывать различные функциональные группы, встречающиеся в лекарствах. Второй подход - менее точный, но более быстрый, -неявно учитывает свободную энергию гидрофобного и электростатического взаимодействия с растворителем [Roux, Simonson, 1999, Biophys. Chem., 78:1], оценивая ее с помощью параметров, которые, как правило, не имеют серьезного физического базиса в описании взаимодействий, и подогнаны под уже существующие константы связывания различных комплексов. В настоящее время для оценок с помощью первого подхода разработано поляризуемое силовое поле QMPFF3, которое подогнано под квантово-механичсские данные, полученные с помощью метода МР2, учитывающего парные, и только парные корреляции электронов [Donchev et al., 2006, J. Chem. Phys., 125:244107]. Оно нуждается в улучшении параметризации в случаях ароматических систем, где точность «парного» метода МР2 недостаточна. Полученное силовое поле 1акже нуждается в тестировании, чтобы показать, что поле, подогнанное исключительно под «вакуумные», ab initio полученные кваисовые данные, хорошо описывает и конденсированные фазы, - что позволит использовать его для оценок энергии связи молекул в комплексах с помощью первого подхода. Оцененные, с помощью второго подхода, свободные энергии комплексов белок-лиганд обычно представляют собой разницу больших величин (следующих из оценки площади поверхности при контакте белка с лиган-дом и без него), что приводит к большой ошибке определения потенциала взаимодействия. Таким образом, становится очевидным, что необходима разработка новых методов для оценки связывания комплексов и совершенствование существующих силовых полей, - в чем и состоит данное исследование и что определяет его актуальность.

Цель и задачи исследования. Целью данной работы является построение силовых полей для оценки взаимодействия молекул в вакууме и в неявно учитываемой воде для моделирования белков и их взаимодействий с другими молекулами.

Соответственно, были поставлены следующие задачи: I) создать базу кристаллографических и термодинамических данных для определения параметров взаимодействия молекул в воде; 2) построить, на основании сведений из данной базы, силовое поле для оценки взаимодействий в неявно заданном водном окружении (назовем его FFSol); 3) улучшить параметризацию ароматических молекул в существующем поляризуемом силовом поле QMPFF3, которое базируется на квантово-механических данных, и провести тестирование результирующего поля.

Научная новизна работы. Для построения и тестирования силовых полей создана наиболее обширная на настоящее время база данных CRAFT по кристаллическим структурам и термодинамике молекул. Впервые продемонстрирована возможность построения силового поля FFSol для описания взаимодействий в "неявно заданном" водном окружении на основании данных по кристаллическим структурам органических молекул, по сублимации их кристаллов и по их растворимости в воде. Продемонстрировано улучшение описания ароматических систем в результате коррекции дисперсионных параметров ароматического углерода в силовом поле QMPFF3 при использовании наиболее точных на сегодняшний день квантово-механических данных CCSD(T). Проведенное тестирование поляризуемого силового поля QMPFF3 на большом количестве органических молекул в различных агрегатных состояниях показало хорошее согласие с экспериментальными данными.

Практическая значимость работы. Полученное силовое поле FFSol может быть использовано для оценки энергии взаимодействия молекул в воде. Улучшенное нами силовое поле QMPFF3 в настоящий момент уже использовано для вычисления свободной энергии комплексов белок-лиганд с помощью методов термодинамического интегрирования [Khoruzhii et al., 2008, Proc. Natl. Acad. Sei. USA, 105:10378] и для оценки энергии связи в наноструктурах [Donchev, 2006, Phys. Rev. В, 74:235401]. Данное силовое поле также может быть применено для других задач молекулярной механики и динамики.

Апробация работы. Материалы диссертации были представлены на российских и международных конференциях: Ежегодная научная конференция МФТИ (Москва-Долгопрудный, Московская область, 2004, 2007); XVII European Conference on Dynamics of Molecular Systems (St.-Petersburg, Russia, 2008); Ежегодная научная конференция Института белка РАН (Пущино, Россия, 2009); IV Российский симпозиум Белки и пептиды (Казань, 2009); International Moscow Conference on Computational Molecular Biology (Moscow, Russia, 2009); The Second Rosnanotech Nanothechnology international forum (Moscow, Russia, 2009); VII Национальная конференция "Рентгеновское, Синхротронное излучения, Нейтроны и Электроны для исследования наносистем и материалов. Нано-Био-Инфо-Когнитивные технологии" (Москва, 2009).

Публикации. Основные результаты диссертации опубликованы в 11 печатных работах, в том числе в 3 статьях.

Структура диссертации. Диссертация изложена на 150 страницах машинописного текста и состоит из введения, четырех глав (обзор литературы, создание силового поля FFSol, коррекция параметров ароматического углерода в СП QMPFF3, тестирование силового ноля QMPFF3), заключения, выводов, благодарностей и списка цитируемой литературы, включающего 206 ссылок. Работу иллюстрируют 11 рисунков и 23 таблицы. Кроме того, работа содержит два приложения: в первом представлена валентная типизация и параметры взаимодействия, а также предсказания силовых нолей FFSol, FFSubl и сравнение их с экспериментом, всего 6 таблиц; во втором приложении представлено сравнение предсказаний силового поля QMPFF3 с экспериментом и квашиво-механнческими данными, всего 6 таблиц.

СОДЕРЖАНИЕ РАБОТЫ Создание «водного» силового поля FFSol

Для ускорения вычислений конформаций и взаимодействий молекул в водном окружении создаются модели межмолекулярного взаимодействия с неявным учетом воды. Как правило, такие модели просто добавляют к "вакуумным" потенциалам отдельно учитываемое гидрофобное и электростатическое взаимодействие. Такое описание не учитывает ван-дер-Ваальсового взаимодействия с "неявной" водной средой. Взяв стандартный вид взаимодействия Леннард-Джонса и изменив его параметры, мы учитываем данный фактор и, впервые используя определяемый константами Генри термодинамический потенциал растворения кристалла в воде, получаем параметры межмолекулярного взаимодействия с неявным учетом воды.

В данном разделе детально описано создание нами силовых полей FFSubl и FFSol, которые основываются на экспериментальных кристаллографических структурах и экспериментальных термодинамических данных. Первое вспомогательное «вакуумное» СП FFSubl, базирующееся на данных по сублимации крист&тлов, строится для тестирования нашей процедуры получения параметров поля, так как способ получения невалентных параметров СП из кристаллических структур отработан уже давно. Второе «водное» СП FFSol, которое базируется на данных по растворению и сублимации кристаллов, неявно учитывает водное окружение, и построение этого нового СП является одной из целей диссертационной работы.

Функциональная форма силового поля невалентных взаимодействий в воде, FFSol, берется в традиционном для СП виде: В„

^noiiccv % )

(1)

где первое слагаемое - стандартный ван-дер-Ваальсов (в форме Леннард-Джонса), а второе -

электростатический (для среды с диэлектрической проницаемостью е) члены. Суммирование в каждом слагаемом идет по всем парам атомов / и /, разделенных тремя и более валентными связями (или находящимся в разных молекулах); г,у=|гу-Г/|, где г, - координата атома ;; д, -парциальный заряд атома /'.

Величины А у, Д,, зависят от типов атомов <, ] и от окружающей среды. Если обеспечивающие «твердость» атомов члены А,/>0 мало зависят от окружения, то с параметрами Ву все иначе. В вакууме Ии>0, но в другой среде это не обязательно. Притяжение атомов / и_/ в воде можно представить как и = -(В, -В1У)(В). )/г°, где Вцг - эффективный параметр притяжения атома к водной среде (см. задачу 3.4 в [Финкельштейн, Птицын, 2005, Физика белка]).

Таким образом, даже в водной среде величины Ау, Ву можно представить в виде произведений членов, относящихся к отдельным атомам:

А^А'А)- В,*В;В-. (2)

Для ускорения вычислений с применением молекулярного поля взаимодействия атомов в формуле (1) обычно «обрезаются» на расстоянии Ис »10А при помощи функций типа

(Г' [О (при г > Дс)"

значения (и производные) которых падают до нуля при г-»Кс. В результате, рассматриваемое нами СП представляется в виде

U

noncov L ' J---V-' -/ J- - \ -j -/j p.

('•;) ь С,;)

В отличие от других работ по СП, где обрезание потенциалов делается уже после оптимизации параметров поля, заданного формулой (1), мы, для повышения точности поля с «обрезанными» взаимодействиями, будем оптимизировать параметры именно такого, уже «обрезанного» (при Rc = 10Â) поля, заданного формулой (4).

«Тип атома», определяющий его свойства в различных взаимодействиях, зависит не только от химического сорта этого атома, но и от его ковалентного окружения. Мы использовали типизацию атомов для невалентных взаимодействий, принятую в СП ENCAD [Levitt et al., 1995, Comput. Phys. Commun, 91:215] (при этом мы ограничились только атомами, встречающимися в белках). В качестве парциальных зарядов всех атомов всех молекул были использованы принятые в силовых полях, рассчитанные на основе квантовой механики «RESP заряды» [Cornell et al, 1993, J. Am. Chem. Soc., 115:9620], дающие оптимальную картину пространственного распределения электростатического потенциала молекулы.

Энергия ковалентных взаимодействий Ucm взята - и для «водного» случая, и для «вакуумного» - в форме и параметризации, присущей СП ENCAD. В тех случаях, когда валентных параметров ENCAD не хватало для описания рассматриваемых в данной работе моле-

кул, недостающие параметры добавлялись по подобию с имеющимися в ENCAD, с учетом равновесных значений длин связей и углов между ними в кристаллах Кембриджской базы данных структур (Cambridge Structure Database, CSD) [Allen, 2002, Acta Cryst., B58:380],

Свободная энергия растворения молекул кристалла в воде складывается из энтальпии сублимации (перехода молекулы из кристалла в пар) и изменения действующего на молекулу потенциала средней силы в воде но сравнению с таковым в парс (причем последний потенциал, в достаточно разреженном паре, практически равен нулю).

Энтальпия сублимации одного моля молекул Л//5иы = //vap - #ciyst, где #vap = Е*ар + РКар и //cryst = Ecr)st + PVcnst, - соотвегствепно, энтальпия моля молекул в насыщенном паре и в кристалле (при этом £vap, Еаyst и Fvap, Fcryst - соответствующие энерг ии и объемы, а Р - давление насыщенного пара). Так как при невысоком давлении PVv,р = RT(где R - газовая постоянная, а Т- абсолютная температура пара), a Kayst « Kvap, то PV,cr>st близко к нулю, то разность энергий молекул в кристалле и паре составляет (в расчете на моль)

E^-E,^-älUabl + RT. (5)

Если Сцц, концентрация растворенных в воде молекул, находится в равновесии с концентрацией cvap = 1/f'vap тех же молекул в паре, то потенциал средней силы, дополнительно (по сравнению с ситуацией в парс) действующий на моль молекул в воде, есть A<Paq-vap = -RT ln(caq/cV4,). Отметим, что величина c,4/cv,р есть, при с->0, константа Генри Таким образом, перенормированная взаимодействием с водой разность «энергий» молекул в кристалле и водном растворе составляет

<^<PCiysl-aq = ¿'cryst ■ Еуар - ДФаЧлаР ~ -ЛЯ5цы + RT+ RT 1п(£н,сс). (6)

Для приведения потенциала (5) к разнице потенциальных энергий основных состояний в кристалле и в газе необходим учет освобожденных степеней свободы в газовой фазе. Пусть у нелинейной молекулы в газе есть «,„, свободных внутренних вращений, кроме того, у нее как целого в газе есть ло=6 свободных трансляционных и вращательных степеней свободы. Значит, дополнительная потенциальная энергия возбужденных степеней свободы молекулы в кристалле на (no+nrai)RT/2 больше, чем в газе, то ссть UCF! st - Uvzp ~ /'cryil ~ ^vap - (nr,^nM)RT/2.

В результате, разность потенциальных энергий основных состояний моля молекул в кристалле и в газе есть

Ai/cmt •vap — ^cryst " ^Лар -&HsM + RT(l -(«о+"го.)/2), (7)

в то время как перенормированная водой разность «потенциальных энергий» основных состояний моля молекул в кристалле и в воде есть

AC/cryst-aq* -AHsM + RT(l -(«o+"rot)/2) + Ärin(*H.a) (8)

Созданне базы данных "Crystals for Force Field Development and Testing" (CRAFT). В данной работе рассматриваются только простые, состоящие из одного типа молекул кри-

сталды, для которых известны или могут быть вычислены (по табличным термодинамическим данным) величины A#Subi и

Для получения «вакуумного» силового поля необходима база данных, состоящая из кристаллических структур, полученных при невысоком давлении и температуре (~1 атм, -300 К), которые можно найти в базе данных CSD, и энтальпий сублимации Д#5иы(7) приведенных к «нормальной» температуре Г= 298 К (большинство табличных данных относятся именно к 25° С, но для веществ, плавящихся при более низкой температуре 7'тыь A#subi(Г) можно вычислить с помощью литературных данных путем экстраполяции от Mfsubi(^mcii); подробнее о разных способах получения ЛЯ5иы(Т) см. [Chickos, Aeree, 2002, J. Phys. Chem. Ref. Data, 31:537]). При этом отобранные молекулы не должны быть очень гибкими (чтобы избежать сильных внутримолекулярных невалентных взаимодействий в газе и в растворе) и не должны ассоциировать или диссоциировать ни в растворе, ни в газе.

Созданная нами база данных является самым большим на сегодняшний день пересечением баз по кристаллическим структурам и соответствующим им термодинамическим данным (в том числе и вычисленным при стандартной температуре 25°С) и может быть использована не только для получения «вакуумного» силового поля FFSubl, но также для тестирования качества других силовых полей, как это сделано для СП QMPFF3.

Для получения «водного» силового поля FFSol, помимо кристаллических структур и энтальпии сублимации, необходимы еще константы Генри кцсс(Г) [Sander, 1999, http://mvw. henrvs-laiv. ore 1 также приведенные к «нормачьной» температуре Т =298 К. Так как множество экспериментальных данных для построения FFSol меньше, чем для FFSubl, а последнее служит тестом процедуры получения парамезров, то в данной работе количество молекул для построения FFSubl было сужено до множества данных, используемых при построении FFSol.

В результате, для получения FFSol и FFSubl было отобрано 58 молекул. Детальное описание молекул с рисунками, термодинамическими параметрами и деталями вычисления этих параметров находится в файле htro://phvs.protres.ru/resources/FFS/A2.pdf

Процедура получения параметров. Параметры искомого силового поля (е, А,, B¡; / = 1, 2,..., М) получаются в результате минимизации различий («невязок») между экспериментальными и вычисленными (с помощью этих параметров поля) величинами. При этом для каждого кристалла к в расчет принимаются:

а) среднее квадратичное различие кристаллографических Зсг, Ьсг, сСг и вычисленных 3eq, b<¡4, ceq равновесных векторов ячеек решетки кристалла к, нормированное на средний диаметр атома S~ ЗА,

A« = [[(.£> -а£)2 +(b¡» -Ь£)г + «>-c«)2]/3]/¿2; (9)

б) среднее квадратичное различие кристаллографических и вычисленных равновесных

6

координат всех /V'',} атомов ячейки кристалла, нормированное на средний диаметр атома ¿КЗА,

/б1- (10)

в) нормированное квадратичное отличие экспериментальной оценки А0"^ «энергии» образования кристалла к из моля отдельных молекул от ее вычисленной величины Д

лГ =[д^-д^>((г},{г''",})]!/[дс/^]г. (И)

Для «вакуумных» потенциалов А(/'"? - ДСм^мар (ДЛЯ кристалла к). а для «водных» потенциалов = ДЦ^^-ач (для того же кристалла); расчет ДУ^ц-ар и Дищ Л-аа СМ. ВЫШе.

Величина

ЛОМДГ^"» = ^{г}-^^'} (12)

вычисляется для данных («водных» или «вакуумных») потенциалов при рассматриваемых координатах г атомов молекул кристалла к и координатах г(/",е) атомов освобожденных из этого кристалла молекул

Таким образом, общая функция невязки определяется суммой невязок, усредненной по всем Л" рассматриваемым кристаллам:

л = д1+д!+дз = 1£д;'ч1£лг+±£д<'>. (13)

Л- 1-.1 Л- ¿.^ Л А.г1

Оптимальные параметры СП (с, А В,', где / = 1, ..., М) получаются в результате процесса минимизации функции невязки Д. Оптимизация включает в себя следующие этапы:

1) задание «затравочных» параметров СП;

2) задание для каждого кристалла кристаллографических координат его ячейки {г^}»» и векторов асг, его решетки;

3) локальную минимизацию величины С/С1у5[ {г} для каждого кристалла путем описанной выше локальной оптимизации набора векторов (г}000 и а,Ь,с при заданных параметрах СП, и получение равновесных для данного кристалла значений ае1,, Ьч, ссч и {г,,,,}^,;

4) локальную минимизацию величины и6а{г</'™>} для одной молекулы, извлеченной из того же кристалла путем аналогичной локальной оптимизации, описывающий ее набор векторов {г(/'"°} при тех же параметрах СП (если в асимметричной части кристаллической ячейки находится несколько молекул, энергия каждой из них минимизируется отдельно, и минимальная из этих энергий берется в качестве С'6ее);

5) расчет по формулам (9) - (13) величины ДсЧ = Д(ач,Ьч,счЛге,}<х>о>{г^'}) (теперь,

7

когда для всех кристаллов и извлеченных из них молекул определены равновесные ач.Ь«,,^, {гч}000 и {rif"°}, Да, является функцией только параметров СП);

6) модификацию параметров СП и повторение (снова стартующее для каждого кристалла с соответствующих [га}т, ает,Ь„,ссг) вышеописанных шагов (3), (4) (5) - с целью поиска параметров СП, отвечающих минимуму величины Д.

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

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

Параметры А',В',е получаются в результате последовательного (по типам атомов) расширения набора молекул, входящих в оптимизацию заданной уравнением (13) функции Д. На первом шаге рассматриваются только молекулы, содержащие атомные типы H и С. Для значений А^,А'Н,В'С,В'Н строится четырехмерная логарифмическая сетка в пределах охватывающих область от половины минимального до удвоенного максимального значения каждого из таких параметров в СП ENCAD [Levitt et al., 1995, Comput. Phys. Commun, 91:215], AMBER99 [Wang et al., 2000, J. Comput. Chem., 21:1049], CHARMM27 [Foloppe, MacKerell, 2000, J. Comput. Chem. 21:86]. Количество точек логарифмической шкалы для каждого параметра берется равным 10. По анализу величин Дсч, полученных для равновесных конфигураций кристаллов и отдельных молекул при параметрах СП, заданных точкой четырехмерной решетки, в этой решетке можно локализовать точку с минимальным Aeq Из этой точки производится локальная минимизация Acq в пространстве параметров поля методом Нелдера-Мида; после этого найденные оптимальные параметры сохраняются. Далее, при добавлении молекул, содержащих помимо типов II и С, также атомный тип А (углерод в sp2 гибридизации), производится, с учетом оптимальных параметров А'С,А'Н,В'С,В"Н, полученных выше, вычисление величин Acq в узлах логарифмической сетки параметров А\,В'Л (полученной на основании исследования таких параметров в других силовых полях, как описано выше), где также берется по 10 точек шкалы для каждого параметра. Потом выбирается точка с минимальным ДеЧ и производится локальная оптимизация Деч в пространстве всех параметров Aç,А'н,А'а,В'., В'п,В'л. Далее, последовательным добавлением молекул, которые содержат атомные типы D, V; О; N; L; S; W получаем все искомые «вакуумные» параметры.

Вся данная процедура проводится для разных значений параметра s [б = 0.25, 0.5, 1, 2, 4, 8, 16,32, 64, 128], и после получения для всех этих s всех искомых «атомных» параметров

8

A'(s),B'(e) все эти параметры оптимизируются методом Пелдера-Мида вместе с параметром s (причем нами было показано, что получающийся сдвиг от начального значения параметра е невелик). После этого выбирается окончательный, с минимальной невязкой Л«,, набор всех параметров А', В',е для «вакуумного» СП FFSubl.

Получение «водных» параметров проводится почти так же, как и «вакуумных». Отметим, что «водные» параметры В] могуг сильно отличаться от вакуумных. В принципе, они могут принимать и отрицательные значения. Поэтому в качестве начальной сетки для каждого В' берется и логарифмическая шкала в положительной области, как описано для «вакуумного» случая, и симметричная ей шкала в отрицательной области, а также нулевое значение параметра В'. В результате получаются параметры А', В',с для «водного» СП FFSol.

Параметры FFSubl и тестирование их качества. Расчет «вакуумных» потенциалов и их тестирование является отправной точкой работы и тестированием метода получения параметров. Мы получили параметры с сохранением и без сохранения кристаллической симметрии в кристаллах. Таблица 1 показывает, что полученные в обоих случаях параметры СП примерно одинаковы, однако общая функция невязки Д (между вычисленными и экспериментальными свойствами кристаллов, см. формулы (9-13)) существенно ниже в случае минимизации с соблюдением кристаллографической симметрии.

Таблица 1. Набор атомных параметров «вакуумного» СП FFSubl, соответствующий минимальной общей функции невязки Л при минимизации с учетом симметрий в кристаллах и без такого учета

Тип незаряженного атома

«Вакуумная» оптимизация FFSubl е учетом _симметрий_

в;

ккал д моль

«Вакуумная» оптимизация FFSubl без учета симметрий

д

Iккал д I моль

В,

ккал д моль

II (Водород неполярный) С (Углерод в «р3 гибридизации) А (Углерод в ер2 гибридизации) О (Водород полярный) V (Кислород в -ОН, -О- группах) О (Кислород в >С=0 группе) N (Азот неароматический) Ь (Азот ароматический) в (Сера в -в-, -БН группах) (Кислород в воде)

106.04 1017.86 1243.03 0.28 406.23 196.96 1254.41 820.44 1759.38 684.14

6.69 21.78 26.58 0.89 17.68 14.81 34.09 35.31 52.64 33.98

103.57 1165.64 1167.43 0.27 425.17 205.36 1273.73 406.94 2267.30 748.08

6.62 22.44 25.65 0.89 18.35 13.69 34.03 26.69 58.76 27.44

е=1.958

е=1.997

Общая функция иевязки

Д =0.021

Д =0.046

Частные функции невязки

0.0052

0.0092

Д.

0.0065

0.021

0.014

А,

0.0109

При этом основной выигрыш достигается за счет лучшего описания геометрии ячейки кристалла. По-видимому, это объяснятся тем, что оптимизация с учетом симметрий работает в пространстве с меньшим числом параметров и ее более «усредненный» путь не застревает в малых «случайных» локальных минимумах энергетической поверхности. Предсказанные СП FFSubl величины &Ucryst.y,р как с учетом симметрий, так и без него, достаточно хорошо согласуются с экспериментом: коэффициент корреляции г = 0.954 в первом случае и /•=0.947 во втором. Полезно также сравнить результаты работы лучшего (полученного с учетом симметрий кристаллов) из наших силовых полей, СП FFSubl, с результатами работы других СП, для которых в научной литературе можно найти соответствующие данные. Сравнение (таблица 2) проводилось с СП CFF II [Ewig et al, 1999, J. Phys. Chem. B, 103:6998] (подогнанным под кристаллические структуры без учета симметрий в них) и с СП QMPFF3 [Donchev et al., 2008, J. Comput. Chem., 29:1242] (подогнанным под квантово-механические расчеты взаимодействий простых молекул, также не учитывающим симметрии кристаллов).

Таблица 2. Относительные точности воспроизведения энергий взаимодействий и объемов кристаллических ячеек, общих для нар силовых полей CFF II - QMPFF3, CFF II - FFSubl и QMPFF3 - FFSubl

Силовые поля К, количество общих кристаллов в двух СП АЕ/Е("Л) Л V/V{"/a)

СП- 11 6.66 нет данных

ГТвиЫ (без учета симметрий) 7 2.80 1.73

РРвиЫ (с учетом симметрий) 3.54 2.04

О.МРГРЗ 7.49 4.09

РРвиЫ (без учета симметрий) 30 6.28 2.52

ГКБиЫ (с учетом симметрии) 6.11 1.97

СТГII (ЗМРРРЗ 21 5.20 5.98 нет данных 3.79

Примечание. ДЕ1Е = £ р^-АС/^)2/(л^)2 /К АУ/У= £ ^-¿ОЧС)''*-

1-1. ,;: ¡.-к к

Параметры РГво! и тестирование их качества. Получение «водных» потенциалов СП РР8о1 является одной из целей данной работы. В таблице 3 даны атомные параметры невалентных взаимодействий для «водного» силового поля РРЭо!, соответствующего минимизации общей функции невязки Д на полном наборе 58 используемых молекул при растворении их кристаллов в воде. Для контроля там же даны аналогичные параметры, полученные при минимизации Д на двух половинках этого набора. "Полунаборы" выделялись таким образом, чтобы каждый представлял приблизительно половину молекул для атомов каждого типа. Первый из наборов состоит из 28 органических молекул, а также молекулы воды, а второй -из оставшихся 29 органических молекул и молекулы воды (которая входит в оба набора, так как только она включает атом типа \У). Во всех случаях минимизация проводилась с учетом симметрий кристаллов, в соответствии с опытом, полученным нами при получении «вакуумных» параметров силового поля РРЭиЫ (таблица 1).

10

Сравнение таблиц 1 и 3 показывает, что параметры «водного» СП ГР5о1 довольно сильно отличаются от параметров «вакуумного» СП ГР5иЫ, причем основная разница связана с ожидаемым резким (по сравнению с РР5иЫ) ростом е (и, соответственно, ослаблением электростатических взаимодействий) в «водном» СП РРБо1 и с уменьшением в нем параметров В , отвечающим притяжению полярных водородных (П) атомов.

Ишересно, что в "вакуумном" СП РРЯиЫ оптимизированная величина г лежит между типичной величиной диэлектрической проницаемости кристаллов органических молекул (~3) и диэлектрической проницаемостью вакуума (1.0), а в "водном" СП РТБо! е лежит между диэлектрической проницаемостью кристаллов органических молекул и воды (=¡80).

Таблица 3. Набор атомных параметров «водного» силового поля РКЯо!, соответствующий минимизации обшей функции невязки А па полном наборе используемых молекул н на двух половинках этого набора.

Тип незаря- Оптимизация «водного» СП 1^01 по всем 58 молекулам Оптимизация «водного» СП РРво! по первой половине набора молекул Оптимизация «водного» СП РРБо! по второй половине набора молекул

женного атома 4' в' К в' 4' Д'

|жал X" V моль мекал д6 V моль 1 ккал д12 Умоль ккал д6 V моль 1 о \2 гасап А V моль V ккал д6 моль

II С Л О V о N ь в \У 102.66 1281.69 1154.92 0.23 722.92 162.13 1309.67 390.02 2014.41 397.50 7.27 22.27 23.82 0.33 27.88 9.85 36.26 17.29 46.81 24.47 107.72 797.50 1134.30 0.50 441.32 59.59 1043.00 355.07 1853.42 416.35 6.95 20.03 24.30 0.65 20.43 5.38 30.39 19.76 46.37 25.60 100.33 1747.57 1077.40 0.25 705.24 166.21 1500.45 1307.84 5778.23 95.83 8.25 22.52 21.32 0.28 26.86 11.52 36.32 29.03 81.32 25.13

£=32.06 £ =32.00 £-32.00

Общая функция невязки А Д=0.049 Д=0.055 Д=0.060

Частные функ- д, | д, 1 Д3 А, Аг д3 А, А; А,

ции невязки 0.0110 10.01721 0.0207 0.0135 0.0204 0.0207 0.0199 |0. 0190 ¡0.0214

Таблица 3 показывает также, что получаемые параметры СП РР5о1 обычно мало зависят от того, на каком наборе проводилась их оптимизация. Эта таблица показывает также, что величина функции невязки Д, полученной на всем наборе молекул, мало отличается от тех, что получены на двух «половинных» наборах. При этом сравнение таблиц 1 и 3 показывает, что функция невязки для «водного» случая вдвое больше, чем для вакуумного, что отвечает увеличению пофешностей приблизительно на 40% (см. ДЕ/Е в таблице 4). Детальное предсказание для всех трех наборов представлено на рисунке 1. Рост невязки в СП РР5о1 естественен, так как наша модель, неявно описывающая водное окружение, заметно менее строга, чем модель, описывающая взаимодействие в вакууме. Гораздо более любопытно то,

11

что этот рост относительно невелик, - это показывает, что используемая нами модель влияния водного окружения далеко не плоха.

Таблица 4 показывает, что и у использующихся «вакуумных» СП СРР II, ОМРРРЗ и у нового «вакуумного» СП ГРЭиЫ и нового «водного» СП РР5о1 относительная погрешность расчета энергии составляет 6.3 - 9.7%, а относительная погрешность расчета объема кристаллической ячейки - 2.8 - 4.1%.

Таблица 4. Точность воспроизведения энергий взаимодействий и объемов кристаллических ячеек (в расчете на одну молекулу) силовыми полями СРР И, ОМРРРЗ, РРБиЫ и КГво! для всех кристаллов, использованных в каждой работе

Силовое поле Количество кристаллов 5£ (ккал/моль) АЕ/Е (%) SV(k3) АК/К(%)

«Вакуумное поле» РРБиЫ (без учета симметрии) 58 1.78 7.91 6.78 2.95

«Вакуумное поле» РРвиЫ (с учетом симметрии) 58 1.48 6.27 5.93 2.82

«Водное ноле» (Т8о1 (с учетом снмметрий) 58 1.71 9.71 5.78 3.22

«Вакуумное поле» СРР II 34 2.44 7.13 нет данных нет данных

«Вакуумное поле» ОМРРРЗ 78 2.14 9.13 8.04 4.14

Примечание. Ж= Г ^(д^-А^^/ЛГ; <Г= I -/К, Z - число моле-

кул в кристаллической ячейке; АЕ/Е, AVIV- см. примечание к таблице 2.

Д^сryst-aq, эксперимент ~

7 1г / 1»м

/

и/U /

У X

33 17

33 г-

/

Параметры были получены по:

■ по полному набору молекул X по первой половине набора □ по второй половине набора

AUcryst-aq, эксперимент

Рисунок 1. Полученная силовым полем РР'йо! величина ДЦ^м-ач в сравнении с ее значением, полученным из опыта по формуле (8). Все единицы в ккал/моль. Слева: параметры СП РР'Бо! рассчитаны по всем 58 молекулам, приведенным в 1Щр://рЬу5.рпЛге5.ги/ге50цгсе5/РРЗ /A2.pdf. Справа: параметры СП РРЭо! рассчитаны: 1) по всем 58 молекулам (точки, коэффициент корреляции г - 0.920); 2) по первой половине полного набора (крестики, г = 0.912); 3) по второй половине полного набора (квадратики, г = 0.920).

Уточнение параметров ароматического углерода CII QMPFF'3

Построение поляризуемого силового поля общего назначения QMPFF3 (Quantum mechanical polarizable force field, version 3), разработанного в компании Algodign, представляет собой весьма объемную задачу. Я принимал участие в двух из ее этапов - в получении параметров взаимодействий ароматических углеводородов (несущих на себе обобществленные многоэлсктронные системы) и в тестировании результатов работы полученного силового поля при расчете различных свойств молекул, в том числе белков, и свойств различных фазовых состояний вещества.

Прежде чем переходить к описанию проделанной нами работы, нужно пояснить принципиально новые особенности силового ноля QMPFF3. В нем есть три отличные от остальных СП особенности:

1) параметризация силового ноля только под квантово-механические данные (КМД), что позволяет избежать зависимости ог ограниченною объема экспериментальных данных;

2) достаточно физически обоснованное описание 4 компонент невалентного взаимодействия (ES - электростатическое, ЕХ - обменно-оггалкивающее, DS - дисперсионное, PL -поляризуемое), в том числе моделирование электронной плотности распределенным экспоненциально спадающим (диффузным) облаком;

3) отдельная параметризация под каждую из четырех компонент взаимодействия, с помощью специальной процедуры декомпозиции энергии ab initio квантово-мехаиического метода МР2 (Moller-Plcsset) [Meiler, Plesset, 1934, Phys. Rev., 46:618].

Дисперсионная компонента энергии взаимодействия двух атомов а и Ъ на расстоянии Кь - jRj > Б силовом поле QMPFF3 определена в следующем виде:

и? (Ks)=-(c.rCV6 (14)

где CDS6 и CDS8 - силовые константы, которые характеризуются атомным типом (причем параметры CDS6 в QMPFF3 аналогичны по смыслу параметрам В в FFSubl), а

2л

Vi„ [r> va' vb) = г'2"(I -g"2W(v-,+v,)^(2г/(v, + vB))'/к!) - функции Танга-Тоениса (определен-

А=1

ные для всех целых п> 3) [Tang, Toennies, 1984, J. Chem. Phys., 80:3726], которые были выведены как члены разложения дисперсионного взаимодействия, учитывающего поправку на перекрывание электронных облаков. Параметры vA и vB являются параметрами ширины электронной плотности для дисперсионного взаимодействия DS. При устремлении параметров ширины к нулю нетрудно видеть, что выражение (14) стремится к классическому дисперсионному выражению [Stone, 2008, Science. 321:787]: C/f (ИЛ)=-(Сff'Cf6/г6+Cf8Cfs/г8)

В компании Algodign [Donchev et al, 2006, J. Chem. Phys., 125:244107] была разработана специальная процедура декомпозиции МР2 энергии на 4 основные части: UEs, Uex, Uds.

13

Upl. Подг онка под различные части энергии, выделяемые из КМД, упрощает процедуру параметризации СП QMPFF3 и предполагает более точное описание каждой из частей. Известно, однако, что в некоторых случаях метод МР2, который учитывает кулоновскую корреляцию только двух электронов, дает некорректное описание обобществленных, многоэлектронных систем, где дисперсионная DS-часть взаимодействий имеет существенное значение. Хорошо известным примером таких систем являются ароматические молекулярные комплексы, в которых из-за обобществленного я-резонансного облака необходимо учитывать корреляцию, как минимум, 6 электронов. Вычисление энергии димера бензола в оптимальной конфигурации в приближении МР2 дает 4.8 ккал/моль, что почти в два раза больше, чем в соответствующих результатах, полученных методом CCSD(T) [Bartlett, ¡989, 93:1697] (самых точных на сегодняшний день квантово-механических данных для димера бензола).

Поэтому нами была проведена дополнительная двухшаговая процедура параметризации СП QMPFF3 для ароматических углероводородов.

Предварительная параметризация на основе МР2 метода. На первом шаге мы произвели стандартную параметризацию QMPFF3 с использованием тренировочного набора, основанного на МР2 расчетах, чтобы найти первичные параметры (далее QMPFF3Mra) путем подгонки ES, EX и DS компонент энергии к соответствующим квантово-механическим (КМ) компонентам энергии в димерах. PL компонента была подогнана под молекулярные поляризуемости и другие КМД (молекулярные свойства, энергии и неаддитивные эффекты в муль-тимерах).

Любая последовательная декомпозиция КМ энергии является достаточно условной процедурой, поэтому интересно сравнить МР2 подход, первоначально использовавшийся в QMPFF3, с альтернативной декомпозицией, основанной на методе SAPT2 (Symmetry Adapted Perturbation Theory) [Moszynski et al, 1996, Molecular Physics, 88:741], которая была применена в работе [Sinnokrot, Sherrill, 2004, J. Phys. Chem. A, 108:10200] к димеру бензола. Таблица 5 отображает дополнительную информацию для Т и PD конфигураций при МР2 оптимизированных геометриях димера бензола (центры масс мономеров разделены на 4.9, 3.8 Â, соответственно).

Во-первых, как следует из таблицы S декомпозиции' энергии SAPT2 и QMPFF3MP2 находятся в разумном согласии друг с другом. Электростатическая (ES), поляризуемая (PL) и обменно-отгалкивающая (ЕХ) компоненты различаются максимум на 0.3 ккал/моль, а дисперсионная (DS) - значительно больше, на 0.7 - 1.6 ккал/моль. Также было показано, что КМ квадрупольные моменты и средние электронные поляризуемости, которые косвенно отвечают за качество электростатической (ES) и поляризуемой (PL) компонент, тоже хорошо согласуются с экспериментом (п пределах 5%). Данный результат можно рассматривать как один из вариантов подтверждения достоверности схемы декомпозиции энергии QMPFF3 для

ЕЭ, ЕХ и РЬ компонент, в то время как дисперсионная компонента нуждается в значительном уточнении.

Таблица 5. Сравнение МР2 и 8АРТ2 схем декомпозиции энергии для Т и РГ) конфигураций димера бензола. Энергии дана в ккал/моль.

Метод декомпозиции ES | EX | DS I PL

T-конфигурация

SAPT2 -2.2 4.9 4.4 -0.7

МР2 -2.2 4.7 1-5.1 -0.7

PD-конфигурация

SAPT2 -2.8 I 8.7 -7.9 -0.9

МР2 -2.9 | 9.0 -9.5 -1.1

Коррекция дисперсионного параметра с помощью данных CCSD(T). Далее, на втором шаге, параметры дисперсионной компоненты DS мы скорректировали так, чтобы описание оптимального (PD) димера бензола совпадало с наиболее точными КМД, которые были найдены в литературе для полной энергии димера бензола в устойчивых конфигурациях. Обнаружено, что хорошее согласие с CCSD(T) вычислениями может быть получено путем варьирования всего лишь одной силовой константы CDS6 для углеродного атомного типа при члене R"6 в уравнении (14), тогда как другие параметры мы оставили фиксированными. Скорее всего похожие результаты могли бы быть получены, если бы изначально был использован метод CCSD(T) для получения КМД для начальной параметризации вместо переопределения на втором шаге. Результирующее СП с улучшенными параметрами, описывающими ароматический углерод, обозначено далее как QMPFF3ccst><T).

Таблица 6. Энергии и межмономерные расстояния стационарных PD конфигураций димера бензола для МР2 и CCSD(T) квантово-механическнх методов и подогнанных под эти КМД соответствующих вариантов СП QMPFF3

Метод R(A) E (ккал/моль) CDS6( V E6 - ккал / моль )

КМД: MP2 3.68 -4.58

СП QMPFF3MFÏ 3.69 -4.51 23.58

КМД: CCSD(T) 3.94 -2.0,-2.48,-2.78 '

СП QMPFF3llm,(1) 3.96 -2.46 18.35

"Данные из работ: [Hobza et ai, ¡996, J. Phys. Chem, 100'l8790\, (Tsuzuki et al., 2001, J. Am. Chetn. Soc., ¡24:104], [Sirmokrot et al., 2002, J. Am. Chem. Soc., 124:10887]

Скорректированная величина CDS6 для ароматического углерода стала равна 18.35 VE6-ккал/моль, что на 28% меньше, чем значение, полученное из МР2 данных (см. таблицу 6). Полное подтверждение улучшения ситуации при использовании QMPFF3CCSD(T) показано на рисунке 2, где полная энергия димера, посчитанная с помощью QMPFF3CCSD(T), сравнивается с CCSD(T) литературными данными для PD, Т и S конфигурации. Можно увидеть, что СП QMPFFЗCCSD(T, достаточно хорошо воспроизводит энергии для наиболее выгод-

ной РО конфигурации, тогда как в двух других случаях есть систематические отклонения до нескольких десятых ккап/моль.

К(Л) R(A) R(A)

Рисунок 2. Зависимость энергии от расстояния R между центрами мономеров в димере бензола для PD, Т, S конфигураций. Полые круги отображают МР2 данные. Пунктирные линии описывают предсказания СП QMPFF3MP2. Сплошные квадраты, ромбы и круги отображают CCSD(T) данные из работ [Hobza et al, 1996, J. Phys. Chem, 100:18790], [Sinnokrot et al., 2002, J. Am. Chem. Soc., 124:10887], [Tsuzuki et al, 2001, J. Am. Chem. Soc., 124:104] соответственно, a сплошные кривые отображают предсказания перепараметризованного СП QMPFF3CCSDm.

Последующее сравнение СП QMPFF3MP2 (полученного только на основе МР2 данных) и СП QMPFF3CCSDm, демонстрирует качество коррекции для ПАУ веществ в трех фазах: газе, жидкой фазе и кристаллах.

Подтверждение качества коррекции параметров ароматического углерода в газовой, жидкой и твердых фазах. На этом шаге представлены результаты вычислений с помощью СП QMPFF3 второго вириального коэффициента fij, который характеризует первую поправку, связанную с дисперсионным притяжением молекул, к уравнению идеального газа с концентрацией п и давлением р при температуре Т: р/ (к„Т) = и + В2 (Т)п2 (кв - постоянная Больцмана). Вычисления выполнялись с помощью метода Монте-Карло, рассматривая молекулы как твердые тела. Среди всех ароматических веществ наше внимание было ограничено только газом бензола, так как экспериментальные данные доступны только для него. На рисунке 3 предсказания QMPFF3 и QMPFF3CCSD(T) сравнены с собранными экспериментальными данными в температурных пределах 300-800К. Как можно видеть, QMPFF3MP2 переоценивает второй вириальный коэффициент в 2 раза, a QMPFF3CCSD(T) обеспечивает лучшее описание, с систематической переоценкой экспериментальных данных всего на -15% для всех исследуемых температур.

Т(К)

Рисунок 3. Второй вириалышй коэффициент для пара бензола вычисленный с помощью QMPFF3MP2 (пунктир) и QMPFF3CCSI^1) (сплошная кривая), крестики отображают экспериментальные данные [Lide, 2001, CRC Handbook of Chemistry and Physics, 82" ed.].

Нами были вычислены также термодинамические свойства жидкого бензола, для которого есть экспериментальные данные. Моделирование с помощью молекулярной динамики было проведено в изотермическом, изобарическом ансамбле для 30Ä кубического бокса с соблюдением периодических условий. Потенциалы QMPFF были обрезаны на расстоянии 9А с использованием стандартной даиьнодейстаующей коррекции дисперсионных взаимодействий [Donchev et al, 2006, Proc. Natl. Acad. Sei. USA., 103:8613]. Уравнения движения интегрировались в течение 300 пикосскунд с фемтосекундым шагом, используя схему скоростей Верлета [Swope et al., 1982, J. Chem. Phys., 76:637] с термостатом Нозе-Гувера [Martyna et al., 1992, J. Chem. Phys., 97:2635] и баростатом Берендсена [Berendsen et al., 19S4, J. Chem. Phys., 81:3684]. Моделируемая область содержала 187 молекул бензола. Энтальпия испарения была посчитана с помощью стандартного соотношения AHvap = Egas-E|iq+RT, где Egas и Ецч представляют собой энергии газового и жидкостного состояния в расчете на одну молекулу. Газовая энергия EgaS была вычислена на основании короткой (100 пс) МД симуляции одиночной молекулы в газе. Таблица 7 сравнивает вычисленные плотности р и энтальпии испарения АНуар для жидкого бензола с экспериментом для температур Т=293 и Т=353К.

Таблица 7. Свойства жидкого бензола для температур Т=293К и 1=353К

Плотность (г/см3) Энтальпия испарения Hvap (ккал/моль) Коэффициент самодиффузии D(10"Vc"')

Т=293К Т=353К Т=293К Т=353К Т=293К Т=353К

QMPFF3m1>2 1.070 1.023 13.4 12.5 0.07 0.4

Q\ll'FK3'lM,(l) 0.881 0.804 8.3 7.52 1.3 3.1

Эксперимент* 0.879 0.81 8.2 7.3 2.0 5.5

* См. работу [Donchev et al., 2006, J. Chem. Phys., 125:244107

Для тестирования твердой фазы энергии когезии Есо-„ (энергии невалентных взаимодействий молекулы со всеми ее соседями в кристалле) и специфические объемы молекул в кристаллических ячейках были отобраны для 15 полиароматических углеводородных кристаллов из созданной нами базы данных CRAFT. QMPFF3MP2 и QMPFF3CCSD(T) предсказания

сравниваются с экспериментальными данными в таблице 8, где показаны (^МРРРЗ результаты для всех рассмотренных веществ, и их сравнение с экспериментальными величинами, как объемами так и энергиями когезии.

Далее энергии когезии были вычислены с помощью добавления к усредненным энтальпиям общепринятой аппроксимирующей величиной 2КТ ЕС0/1-=Н1„ь^2РТ без учета внутренних вращений, так как в полиароматических соединениях нет внутренних свободных степеней свободы. Дополнительные неточности относятся к энергии нулевых колебаний ячейки, оцениваемой для ПАУ как 0.5 ккал/моль, что не учитывается в наших вычислениях.

Таблица 8. Энергии когезии и объем ячейки на одну молекулу (объем молекулы) в полиароматических кристаллах

Формула Вещество Энергия когезии Е„Дккал/молъ) Объем молекулы V(Ä3)

QMPFF3 Эксп.а QMPFF3 Эксп.

МР2 CCSD(T) МР2 CCSD(T)

С6Н6 Бензол 18.6 12.4 11.3 108 117 123

C\oHg Нафталин 28.7 18.5 17.8 160 173 178

СнНю Фепаптрен 36.5 23.4 21.8 219 236 242

С14Н10 Антрацен 38.9 24.9 23.9 212 230 238

С1&И10 Флуорантен 39.4 24.9 22.6 239 259 268

С|бНш Пирен 40.1 24.9 24.9 234 253 263

С18Н,2 Трифенилен 45.8 28.4 30.8 265 286 290

С|8Н,2 Кризен 47.9 30.6 32.0 268 288 294

С,8Н,2 Тетрацен 50.1 32.0 34.9 264 285 291

С20Н12 Перилен 52.0 32.3 32.9 279 300 309

С22Н12 Бензоперилен 54.1 33.2 31.8 303 327 337

С22Н14 Пицен 57.7 36.8 34.5 322 347 353

С22Н14 Пснтацен 60.8 38.6 38.4 315 341 343

С24Н12 Коронен 58.1 33.9 35.0 310 348 358

С32Н14 Овален 76.7 45.0 51.8 404 436 448

8 Усредненные экспериментальные энергии когезии, т.е. (для молекул без внутренних степеней свободы) усредиегшые энтальпии сублимации из [Chickos, Aeree, 2002, J. Phys. Chem. Ref. Dala, 31:537] плюс 2кТ

Если проанализировать таблицу 8, можно заметить, что там находятся 4 изомерных группы: СнНю, С1бНю, С^Нц и С22Н14 - в каждой два или три вещества. Можно видеть, что QMPFF3ccsd(t) прекрасно воспроизводит порядок изомеров одновременно и по их энергиям когезии, и по объемам молекул.

Качество нового дисперсионного параметра для ароматического углерода также подтверждается последующим успешным моделированием с помощью QMPFF3 кристаллов графита и фуллерена и их взаимодействий с полиароматическими веществами [Donchev, 2006, Phys. Rev. В, 74:235401], а также точностью предсказывания свободных энергий связывания белков с лигандами, содержащими ароматические кольца [Khoruzhii et al., 2008, Proc. Natl. Acad. Sei. USA, 105:10378]. Таким образом, использование наилучших на сегодняшний день КМ данных может служить надежным базисом для разработки реалистичных высоко-

переносимых СП, которые могут быть применены в широком спектре научных н технологических проблем.

Тестирование СП ОМРРТЗ

Нами было проделано всестороннее тестирование СП (}МРРРЗ (в варианте дМ1'РРЗссадт)) молекулярных свойств в различных агрегатных состояниях на самом широком спектре органических молекул, чтобы подтвердить хорошую переносимость СП ОМРРР'З на все фазы.

Тестирование ОМРРКЗ: Предсказание свойств молекул с помощью СП ОМРРГЗ.

Сравнение предсказания С>МРРРЗ с КМД и с экспериментальными данными но средней поляризуемости и величине днпольных моментов представлено в таблице 9, которая отображает статистические характеристики отклонений ОМРРРЗ от эксперимента и от КМД. Можно видеть, что относительное среднее отклонение для рМРРРЗ от эксперимента является следствием неточности КМД.

Таблица 9. Статистика относительных отклонений (%) средних электронных ио-ляризуемостен и днпольных моментов для всех молекул, для которых есть экспериментальные данные

QMPFF3 от Клад | QMPFF3 от эксп. | КМД ог эксп.

Средняя поляризуемость (57 молекул)

S -0.2 % -4.0% -3.8 %

RRMSD 0.9 % 6.0 % 5.7 %

Дипольныи момент (78 молекул)

5 0.4 % 3.4% 2.9 %

RRMSD И.5% 15.5% 8.4 %

8 = (х!х0-\). ЯКККО=^(х7ха-1)2 , где .г- данные дМРРРЗ либо КМД, ха - данные КМД (МР2) либо эксперимент

Тестирование ОМРРРЗ: Свободная энергия Гиббса димеризации молекул в газах представляется в виде: О - -КТ1п(-Вг/У111), где Вг - второй вириальный коэффициент (его вычисление описано выше), а Ум - молярный объем идеального газа при температуре Т. Результаты представлены на рисунке 4.

Рисунок 4. Сравнение экспериментальных «-> 1 Гиббсовых свободных энергий димеризации ь, молекул [Lide, 2001, CRC Handbook ojChemis- g try and Physics, S2"d ed с вычисленными для 52 О веществ. Все величины даны в ккал/моль. Корреляция между эксперимепг-альными и вычисленными нами с помощью СП QMPFP3 Q g величинами равна 0.89. Систематический сдвиг 0.124 ккал/моль

0.4

19 0.4 0.S 0.8 1.0 1.2 1.4

Эксперимент

Свободная энергия димеризации

Тестирование ОМРЬТЗ: Моделирование однородных жидкостей; сравнение с результатами других СП. Моделирование МД было проведено в изотермическом, изобарическом ансамбле с такими же условиями, как и для бензола (см. выше). Сравнение с экспериментом представлено на рисунке 5.

Энталышя испарения Плотность жидкостей

4 7 10 13 16 19 22 0.4 0.6 0.8 1.0 1.2 1.4 Эксперимент (ккал/моль) Эксперимент (г /см^)

Рисунок 5. Сравнение вычисленных с помощью СП (ЗМРБРЗ свойств однородных жидкостей с экспериментальными значениями (см. данные в \Donchev е/ а1, 2008, 3. Сотрм. СЬет., 29:1242]). Слева: сравнение 50 экспериментальных плотностей (г/см3) с вычисленными (корреляция 0.97, систематический сдвиг 0.013 г/см3). Справа: сравнение 102 экспериментальных энтальпий испарения жидкостей (ккал/моль) с вычисленными (корреляция 0.97, систематический сдвиг -0.22 ккал/моль).

Таблица 10. Сравнение экспериментальных плотностей и энтальпий испарении для жидкостей, вычисленных с помощью СП СМРГГЗ, ГО-СНАКММ, 1ТГ-ОР1$

Энтальпия испарения Н„„ (ккал/моль) Плотность р (г/см3)

Вещество Т(К) ОМРККЗ РГР-ОРГ.У % ОХ и. < в и Эксперимент ОМРП'З РРР-ОРЬв ' РО- | СНАИММ Эксперимент

Этан 184 3.52 3.32 3.60 3.51 0.562 0.529 0.547 0.546

Пропан 231 4.42 4.79 4.02 4.55 0.594 0.591 0.585 0.581

Бутан 271 5.22 5.62 5.21 5.36 0.611 0.614 0.606 0.602

Метанол 298 8.56 8.84 8.91 8.95 0.826 0.794 0.767] 0.791

Ацетон 298 6.53 7.92 7.91 7.41 0.758 0.756 0.780 0.790

Формамид 293 16.60 15.50 16.35 15.70 1.262 1.173 1.067 1.133

н-Мегилацетамид 373 12.97 13.90 13.91 13.80 0.866 0.946 0.868 0.894

Метантиол 279 6.19 5.96 7.20 6.15 0.892 0.864 0.904 0.888

Этантиол 298 6.44 6.73 7.80 6.52 0.827 0.846 0.871 0.832

Бензол 298 8.30 8.20 7.50 8.09 0.881 0.873 0.884 0.877

1Ш80 - 0.50 0.24 0.65 - 0.045 0.026 0.028 -

КМЬО - корень среднеквадратичного отклонения от эксперимента (см. [Вопскеу е/ а!., 2008, 1 Сотрм. СИет., 29:1242])

В таблице 10 представлены свойства жидкостей, вычисленные с помощью СП QMPFF3, СП FQ-CHARMM [I'atel, Brooks, 2004, J. Comput. Chem., 25:1} и СП PFF-OPLS [Kaminski et al., 2004, J. Phys. Chem A, 108:621], а также сравнение с экспериментальными данными для 10 жидкостей, рассмотренных в этих статьях.

Тестирование QMPFF3: Оптимизация энергии кристаллов.

Экспериментальные энергии когезии оценены стандартным для твердых молекул образом как Ecoh - Ifsub 2RT, где Я,„г, -энергии сублимации, взятые из нашей базы данных CRAFT.

Энергия когезии Объем ячейки на одну молекулу

Рисунок 6. Сравнение предсказаний СП QMPFF3 с экспериментальными значениями из нашей базы CRAFT для 78 кристаллических структур. Слева: сравнение энергии когезии кристалла в ккал/моль (корреляция 0.975, систематический сдвиг 1.74 ккал/моль). Справа: сравнение объема ячейки на одну молекулу в А3 (корреляция 0.998, систематический сдвиг 1.58

А3)

Моделирование кристаллов было выполнено путем минимизации энергии кристаллической ячейки без сохранения симметрии кристалла с помощью метода L-BFGS |Liu, Nocedal, 1989, Math. Program., -/J.-503], начиная с экспериментальных параметров ячейки и внутренних атомных координат. На рисунке 6 дано сравнение вычисленных предсказаний QMPFF3 для двух основных свойств кристалла: энергии когезии и объема одной молекулы в кристалле.

Тестирование QMPFF3: Минимизация энергии белков. Можно ожидать, что белковые структуры, оптимизированные даже без учета окружающей белок воды, могут быть близки к нативному белку; поэтому среднеквадратичное отклонение (RMSD) между экспериментально определенной и минимизированной по энергии геометрией может быть рассмотрено как один из критериев качества СП. Для проведения этого теста мы взяли из Protein Dala Bank (PDB) 22 белковых структуры, оптимизированные с помощью СП PFF-OPLS

{Kaminskietal., 2002, J. Comput. Chem., 23:1515] и СП FQ-CHARMM [Pate! et al.,2004, ./.

21

Comput. Chem., 25:1504]. Молекулярно-механическая оптимизация с помощью CII QMPFF3 была проведена методом L-BFGS; геометрии PDB были взяты как начальные. Атомы водо-родов были добавлены с нейтрализацией всех аминокислотных остатков, как это сделано в работах для других СП. Полученные отклонения оптимизированных от экспериментальных структур представлены в таблице 11.

Таблица 11. Отклонение RMSD (А) для белковых структур, оптимизированных в вакууме с помощью СП FQ-CHARMM, СП PFF-OPLS, СП QMPFF3. ВВ (BackBone) -отклонение вычисленное по пептидному скелету, Heavy - по всем тяжелым атомам.

QMPFF3 PFF-OPLS FQ-CHARMM

PDB BB Heavy BB Heavy BB Heavy

1сс5 1.40 1.56 1.91 2.22 1.86 1.96

lcrn 0.90 0.99 1.63 1.83 0.98 1.12

lctf 0.64 0.88 1.65 2.18 1.18 1.84

ldur 0.88 1.04 1.94 2.43 1.52 1.81

lpaz 0.80 1.04 1.04 1.71 1.25 1.63

lpcy 0.69 0.81 1.05 1.77 1.22 1.53

lpax 2.22 2.46 2.15 4.02 1.41 1.81

Ippt 0.65 1.07 1.70 1.90 1.58 1.79

lr69 0.70 0.90 1.13 1.70 1.21 1.60

lrnt 0.95 1.15 1.30 2.30 1.25 1.45

lsn3 0.90 1.08 1.17 2.28 0.89 1.32

lubq 1.10 1.21 1.16 1.97 2.15 2.29

2cdv 2.00 2.20 2.37 3.53 2.37 2.76

2gn5 2.15 2.86 2.06 2.53 2.34 3.00

2lzm 1.26 1.54 1.37 1.81 1.75 2.09

2ovo 0.99 1.12 1.55 1.76 1.61 2.06

2rn2 1.34 1.60 1.09 1.54 1.36 1.93

3\cb 0.98 1.11 1.47 1.82 1.70 2.13

4fdl 0.89 1.22 1.76 2.48 2.04 2.52

4pti 0.83 1.06 1.58 2.30 1.32 1.74

5c pv 1.07 1.19 1.30 1.97 1.41 1.70

7rxn 0.82 1.16 1.31 1.74 1.30 1.74

Среднее 1.10 1.33 1.53 2.17 1.53 1.90

выводы

1. Созданная нами база данных CRAFT является самой обширной на сегодняшний день базой, содержащей экспериментальные данные по органическим молекулам, для которых есть сведения и по кристаллографическим структурам, и по термодинамике растворения и сублимации. База данных CRAFT может быть использована для создания и тестирования различных силовых полей.

2. На основании данных из базы CRAFT построено силовое поле FFSol для оценки потенциалов взаимодействия молекул в неявно заданном водном окружении, а также силовое поле FFSubl для оценки энергии взаимодействия молекул в вакуумном окружении.

3. С помощью точных квантово-механических расчетов энергии димеров бензола проведена коррекция дисперсионного параметра взаимодействий ароматического углерода в кван-тово-механическом поляризуемом силовом поле QMPFF3. При этом показано улучшение описания свойств полиароматических углеводородов во всех фазах.

4. Проведено тестирование СП QMPFF3 в различных агрегатных состояниях, для наибольшего на сей день числа органических вен;еств> путем сравнения предсказаний этого СП с экспериментальными данными, а также с предсказаниями других СП. По всем критериям СП QMPFF3 оказалось сравнимо с другими современными силовыми полями, а по ряду характеристик превосходит их.

Список публикаций по теме диссертации

1. Donchev, A.G., Gallun, N.G., Pereyaslavets, LB., Tarasov, V.l., Quantum mechanical polarizable force field (QMPFF3): Refinement and validation of the dispersion interaction for aromatic carbon, i. Chem. Phys., 2006. 125. P. 244107.

2. Donchev, A.G., Galkin, N.G., Illarionov, A.A., Khoruzhii, O.V., Olevanov, М.Л., Ozrin, V.D., Pereyaslavets, L.B., Tarasov, V.l., Assessment of performance of the general purpose polarizable force field OMPFF3 in condensed phase. 1. Comput. Chem., 2008. 29. P. 1242-1249.

3. Переяславец, Л.Б., Финкельштейн, A.B., Разработка нового, основанного на растворимости молекулярных кристаллов атомного силового поля FFSol для расчета взаимодействий молекул в водном окружении. Молекулярная биология, 2010. 44(2). Стр. 340-354.

4. Переяславец, Л.Б., Финкельштейн, A.B., Разработка модели межмолекулярного взаимодействия на основе водной растворимости молекулярных кристаллов. XLVII Научная конференция МФТИ. 2004. Стр. 53. Россия, Долгопрудный.

5. Переяславец, Л.Б., Дончев, А.Г., Галкин, Н.Г., Тарасов, В.И., Процедура коррекции параметров дисперсионного взаимодействия для ароматического yziepoda в квантово-

механических поляризуемых потенциалах. 50-я научная конференг/ия МФТИ. 2007. Стр. 185186. Россия, Москва-Долгопрудный

6. Donchev, A.G., Galkin, N.G., Illarionov, А.А., Khoruzhii, O.V., Olevanov, M.A., Ozrin, V.D., Pereyaslavets, L.B., Tarasov, V.I., Assessment of Performance of the Ceneral Purpose Polarizable Force Field QMPFF3 in Condensed Phase. MOLEC-2008 XVII European Conference on Dynamics of Molecular Systems. 2008. P. 112. Russia, St. Petersburg.

7. Финкелыптейн, A.B., Переяславец, Л.Б., Разработка нового силового поля "FFS" для оценки взаимодействия белков в воде, рассчитанного на основе растворимостей молекулярных кристаллов. IVРоссийский симпозиум Белки и пептиды. 2009. Стр. 37. Казань.

8. Переяславец, Л.Б., Финкелыптейн, А.В., Оценка свободной энергии связывания молекул при помощи нового силового поля "FFS", полученного на основе растворимостей молекулярных кристаллов в воде. IV Российский симпозиум Белки и пептиды. 2009. Стр. 390. Казань.

9. Finkelstein, A., Pereyaslavets, L., A new atomic force field "FFS" for protein interactions, computed from solubility of molecular crystals in water. Moscow Conference on Computational Molecular Biology. 2009. P. 101. Russia, Moscow.

10. Pereyaslavets, L.B., Prokudin, M.S., Development of a new force field for modeling (-CFr )„-like perfluoropolymers on the basis of a precisely defined conformational and energetic properties of perfluoroalkanes calculated by quantum-mechanical methods. Rosnanotech Nanothechnology international forum. Abstracts. The Second International Competition of Scientific Papers in Nanotechnologyfor Young Researchers. 2009. P. 202. Moscow.

11. Переяславец, Л.Б., Финкелыптейн, А.В., База данных кристаллических структур и термодинамических величин для создания и тестирования "вакуумных" и "водных" силовых полей. VII Национальная конференция "Рентгеновское, Синхротронное излучения, Нейтроны и Электроны для исследования наносистем и материалов Нано-Био-Инфо-Когнитивные технологии.'" 2009. Стр. 508. Москва.

ООО «Фотон-век», ИНН 5039008988 г.Пущино, Московская область АБ-21/145 тел. (4967) 73-94-32 beornot@rambler.ru Мр://р!ю(оп-чек. пагоб. ги

Подписано в печать 25.04.2010 Формат 60x90/16 Тираж 100 экз.

Содержание диссертации, кандидата физико-математических наук, Переяславец, Леонид Борисович

Список сокращений.

Введение.

Глава 1. Обзор литературы.

1.1. Классические силовые поля.

1.1.1. Стандартный вид невалентных взаимодействий.

1.1.2. Валентные взаимодействия.

1.1.3. Явные и неявные модели воды.

1.2. Поляризуемые силовые поля.

1.2.1. Модели с точечными диполями.

1.2.2. Модель флуктуирующих зарядов.

1.2.3. Модель на основе осцилляторов Друде.

1.2.4. Поляризуемость на основе модели Толе.

1.2.5. Возможные улучшения поляризуемой модели.

1.3. Методы получения параметров СП из экспериментальных данкых. 19 1.3.1. Методы моделирования различных агрегатных состояний вещества.

1.4. Квантовые вычисления и их использование для построения СП.

1.4.1. Уравнение Шредингера.

1.4.2. Метод Хартри-Фока.

1.4.3. Наборы базисных функций.

1.4.4. Метод МР2 и его модификации.

1.4.5. Другие квантовые методы и их применение.

1.5. Функциональная форма СП QMPFF3.

1.5.1. Зарядовая плотность QMPFF3.

1.5.2. Электростатическое взаимодействие.

1.5.3. Обменно-отталкивающее взаимодействие.

1.5.4. Дисперсионное взаимодействие.

1.5.5. Поляризационное взаимодействие.

1.6. Параметризация QMPFF3 при помощи квантово-механических данных.

1.6.1. Квантовые вычисления для построения QMPFF3.

1.6.2. Декомпозиция МР2 энергии для выделения четырех компонент энергии.

1.7. Типизация и процедура параметризации QMPFF3.

Глава 2. Создание силового поля FFSol, неявно учитывающего водное окружение молекул.

2.1 Поле FFSol и необходимые для его параметризации экспериментальные данные.

2.1.1 Функциональная форма силового поля FFSol.

2.1.2 Свободная энергия растворения молекул кристалла в воде.

2.1.3 Вклад заторможенных в кристалле степеней свободы.

2.1.4 Необходимые экспериментальные данные.

2.1.5 Создание базы данных CRAFT.

2.2 Оптимизация энергии кристалла.

2.2.1 Типизация и заряды СП.

2.2.2 Энергия кристаллической ячейки и ее минимизация.

2.3 Получение невалентных параметров СП.

2.3.1 Функция невязки и ее минимизация.

2.3.2 Процедура получения параметров.

2.4 Тестирование силового поля.

2.4.1 Параметры вспомогательного силового поля FFSubl и тестирование их качества.

2.4.2 Параметры «водного» силового поля FFSol и тестирование их качества.

Глава 3. Уточнение параметров взаимодействий ароматического углерода в СП

QMPFF3.

3.1 Коррекция параметров ароматического углерода в СП QMPFF3.

3.1.1 Предварительная параметризация на основе МР2 метода.

3.1.2 Коррекция дисперсионного параметра с помощью данных CCSD(T).

3.2 Подтверждение качества коррекции.

3.2.1 Газовая фаза.

3.2.2 Жидкая фаза.

3.2.3 Кристаллическая фаза.

3.3 Качество силового поля определяется качеством лежащих в его основе квантово-механических данных.

Глава 4. Тестирование силового поля QMPFF3.

4.1 Предсказание QMPFF3 свойств молекул.

4.1.1 Поляризуемость.1.

4.1.2 Дипольный момент.

4.2 Энергия Гиббса димеризации в газах.

4.3 Моделирование однородных жидкостей и сравнение с другими СП.

4.4 Оптимизация кристаллов и сравнение с другими СП.

4.5 Вакуумная минимизация белков как тест силового поля.

Введение Диссертация по биологии, на тему "Атомные силовые поля FFSol и QMPFF3"

В настоящее время сфера биотехнологий и создания новых материалов приближается к новому технологическому прорыву. Для предсказания новых свойств материалов, для моделирования молекулярных, и, в частности, биомак-ромолекулярных систем все чаще обращаются к компьютерному эксперименту, использующему силовые поля. Силовым полем (СП) называют набор потенциальных функций и параметров взаимодействий, описывающих валентное и невалентное взаимодействие атомов в молекулах и между молекулами [1]. С помощью молекулярной динамики (МД) можно моделировать как простые жидкости [2-7], так и сложные системы, включающие несколько макромолекул и растворитель [8]. В последнем случае, реалистичность описания требует учета очень большого числа молекул растворителя, что - если молекулы рассматриваются в явном виде - резко, на порядки увеличивает время расчетов.

Таким образом, разработка силовых полей наталкивается на две принципиальные трудности. С одной стороны, СП должно работать достаточно быстро. С другой стороны, оно должно достаточно точно описывать рассматриваемые взаимодействия. Точность описания разнообразных систем с помощью МД зависит от качества СП. Высокая точность описания систем необходима, в особенности, в том случае, когда нужно выделить одну лучшую структуру из большого множества: сотню лучших лигандов из миллиона, потенциально способных связаться с целевым белком [9-12]; стабильную нативную структуру белка среди множества других [13]; наиболее подходящие, с точки зрения протонной проводимости, сульфат-ионные (-R-SO3") боковые группы в перфторполимерной мембране [14, 15] и т.п.

Точность современных СП, даже модифицированных специально под белки, до сих пор не позволяет описать моделируемые белки в воде при помощи МД с расхождением (по Са атомам) от первоначальной, экспериментально-определенной структуры менее, чем на 1.5 А [16]. Последний раунд соревнования CASP по предсказанию белковых структур по их аминокислотным последовательностям не показал существенного прогресса за последние годы, следовательно, такие предсказания все еще нуждаются в более точных, чем существующие, силовых полях [17]. Последнее десятилетие было посвящено разработке поляризуемых СП, в которых учитывается влияние электрического поля на поляризуемые атомные облака, так как только в таких системах может возникнуть корректная диэлектрическая проницаемость [18], а также другие эффекты, не являющиеся попарно-аддитивными [19-21], что является необходимым условием для моделирования таких сложных гетерогенных объектов, как белки.

Недавно было создано поляризуемое СП QMPFF3 общего назначения (т.е. для широкого круга задач), которое базируется исключительно на квантовых данных, что позволяет ему не зависеть от все еще недостаточного количества экспериментальных данных. QMPFF3 имеет сложную, физически-обоснованную функциональную форму. Чем сложнее и физически-обоснованнее функциональная форма силового поля, тем, потенциально, точнее может оно описать энергию взаимодействия в различных системах, и тем больше требований предъявляется к точности подбора численных параметров такого СП.

В частности, силовое поле QMPFF3 не позволяет достоверно описывать ароматические структуры, а также полиароматические углеводороды (ПАУ) [22]. Правильное описание взаимодействия ПАУ является существенной частью разработки СП общего назначения, так как ПАУ имеют большое значение для решения широкого спектра научных и технологических задач (например, разработка лекарств [12], самообразование графитных нанопроволок [23] и т.п.) и являются частями белков, ДНК, РНК, многих лекарств, наноматериалов и т.п. Большую важность ПАУ представляют и для медицинских и экологических наук из-за своей существенной канцерогенности и токсичности [24].

При исправлении недостатка, связанного с плохим описанием ароматических структур, и после качественного тестирования, которое сможет подтвердить хорошую предсказательную силу СП QMPFF3, в том числе и для белков, его можно было бы использовать для задач моделирования макромолекул, а также для предсказания свободных энергий связывания комплексов белок-лиганд [12].

Детальные и сложные СП обычно требуют существенно больше времени для вычислений, чем простые. Поэтому для МД моделирования сворачивания белков (притом, что, характерное время их сворачивания - секунды-миллисекунды, а моделирование таких процессов занимает многие тысячи часов) необходимы «быстрые» модели, желательно учитывающие воду неявно, и при этом правильно описывающие все основные компоненты взаимодействий. Несколько моделей взаимодействий с неявным учетом воды уже было представлено до сегодняшнего дня [25-29], однако до сих пор они не столь точны, чтобы использовать их для выделения нативной структуры белка. Поэтому очень интересной и важной задачей представляется создание такого силового поля для невалентных взаимодействий, которое систематически учитывало бы взаимодействие с водой в неявном виде (так, как учитывается диэлектрическая проницаемость в макроскопических моделях веществ), т.е. было бы основано на принципиально новом для атомных силовых полей подходе.

Целью данной работы является построение силовых полей для оценки взаимодействия молекул в вакууме и в неявно учитываемой воде для моделирования белков и их взаимодействий с другими молекулами.

Соответственно, были поставлены следующие задачи: 1) создать базу кристаллографических и термодинамических данных для определения параметров взаимодействия молекул в воде; 2) построить, на основании сведений из данной базы, силовое поле для оценки взаимодействий в неявно заданном водном окружении (назовем его FFSol); 3) улучшить параметризацию ароматических молекул в существующем поляризуемом силовом поле QMPFF3, которое базируется на квантово-механических данных, и провести тестирование результирующего поля.

Заключение Диссертация по теме "Биофизика", Переяславец, Леонид Борисович

Выводы

1. Созданная нами база данных CRAFT является самой обширной на сегодняшний день базой, содержащей экспериментальные данные по органическим молекулам, для которых есть сведения и по кристаллографическим структурам, и по термодинамике растворения и сублимации. База данных CRAFT может быть использована для создания и тестирования различных силовых полей.

2. На основании данных из базы CRAFT построено силовое поле FFSol для оценки потенциалов взаимодействия молекул в неявно заданном водном окружении, а также силовое поле FFSubl для оценки энергии взаимодействия молекул в вакуумном окружении.

3. С помощью точных квантово-механических расчетов энергии димеров бензола проведена коррекция дисперсионного параметра взаимодействий ароматического углерода в квантово-механическом поляризуемом силовом поле QMPFF3. При этом показано улучшение описания свойств полиароматических углеводородов во всех фазах.

4. Проведено тестирование СП QMPFF3 в различных агрегатных состояниях, для наибольшего на сей день числа органических веществ, путем сравнения предсказаний этого СП с экспериментальными данными, а также с предсказаниями других СП. По всем критериям СП QMPFF3 оказалось сравнимо с другими современными силовыми полями, а по ряду характеристик превосходит их.

Благодарности

Прежде всего, хочу выразить мою искреннюю благодарность моим научным руководителям Виктору Васильевичу Зосимову и, особенно, Алексею Витальевичу Финкелыптейну за поддержку, советы и опыт, который они передали мне в течение своего руководства.

Я благодарю моего научного руководителя Алексея Витальевича Фин-кельштейна, моего друга Маргариту Геннадиевну Никитину и мою жену Елену за внимательное, неоднократное прочтение данной диссертации, высказанные пожелания и советы для ее улучшения. Большое спасибо Дмитрию Петровичу Харакозу за интерес к моей работе, ценные замечания и советы по ее редактированию.

Различные части данной работы были выполнены в компании Алгодайн и в Лаборатории физики белка Института белка РАН. Считаю своим приятным долгом поблагодарить замечательных сотрудников компании Алгодайн с которыми я работал над разными проектами, в первую очередь над созданием СП QMPFF3: Виктора Васильевича Зосимова, Александра Георгиевича Дончева, Владимира Драгановича Озрина, Владимира Ивановича Тарасова, Олега Владимировича Хоружего, Николая Георгиевича Галкина, Михаила Александровича Олеванова, Алексея Аркадьевича Илларионова, Олега Владимировича Бутина, а также Сергея Константиновича Нечаева, принявшего участие в первом этапе разработки СП FFSol и СП FFSubl.

Также хочу поблагодарить теоретическую часть Лаборатории физики белка: Алексея Витальевича Финкельштейна, Оксану Валериановну Галзитскую, Дмитрия Николаевича Иванкова, Наталью Сергеевну Богатыреву, Сергея Александровича Гарбузинского, Михаила Юрьевича Лобанова, Никиту Владимировича Довидченко.

Каждый из них помогал советом, давал прекрасные комментарии, ставил интересные вопросы. Конечно, все перечисленные - не единственные, кто поддерживал меня в трудные минуты; остальные члены коллективов компании Алгодайн и Лаборатории физики белка создали ту незабываемую атмосферу, в которой было так интересно и приятно работать.

Особенно хотелось бы еще раз поблагодарить Александра Георгиевича Дончева, Алексея Витальевича Финкелыптейна и Олега Владимировича Хору-жего, каждый из которых по своему стал для меня «идеалом» ученого и у каждого я научился многому из того, что я умею, и без влияния, которых данная работа никогда не была бы закончена. С благодарностью и скорбью хочу упомянуть

Александра Михайловича Дыхне, который очень сильно повлиял на, то чем мне пришлось заниматься в науке.

У каждого дела есть начало, которое, если посмотреть с конца практически незаметно. Но еще в школе мои учителя, а именно Николай Васильевич Дзумедзей, Валентина Николаевна Дончик), Светлана Сергеевна Жугина привили мне такую сильную любовь к естествознанию, что я выбрал ту дорогу, по которой иду.

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

Заключение

Создание силовых полей очень важно для моделирования структуры и поведения различных макромолекул и предсказания свойств разнообразных материалов, в первую очередь для предсказания констант связывания белок-лиганд, что в свою очередь имеет большое значение для медицины, так как предварительный теоретический отбор тысяч лучших лигандов из миллионов позволяет существенно сэкономить на поиске новых лекарств к целевому белку. Данная работа посвящена созданию, параметризации и тестированию двух силовых полей: FFSol и QMPFF3, которые могут быть применены для оценок связывания комплексов белок-лиганд.

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

Методика получения силового поля на основе экспериментальных кристаллических структур была оттестирована с помощью построения вспомогательного СП «классического» типа FFSubl. В процессе получения параметров мы впервые использовали полную поатомную минимизацию энергии кристаллов с учетом их симметрий, что, помимо ускорения вычислений, приводит к сглаживанию функции невязки, что, в свою очередь, приводит к более точному нахождению глобального минимума. Данное «вакуумное» силовое поле описывает имеющиеся энтальпии сублимации с точностью ~8%.

Методика получения параметров для «водного» силового поля FFSol была протестирована путем деления тренировочного набора кристалла на две равных части. Силовые поля, полученные на половинках наборов, дают почти такое же, по качеству, описание кристаллических структур. Созданное нами силовое поле FFSol описывает экспериментальный термодинамический потенциал, учитывающий и внутреннюю энергию, и сольватацию молекул с суммарной точностью -10%, и в дальнейшем может быть использовано для быстрой оценки свободной энергии связывания молекулярных комплексов.

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

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

СП QMPFF3 было протестировано путем сравнения предсказаний энергии Гиббса для взаимодействия молекул в газе с экспериментом. Мы провели моделирование также большого количества жидкостей с помощью молекулярной динамики и сравнили энтальпии испарения молекул и плотности жидкостей с экспериментом. Минимизация энергии множества кристаллических структур позволила нам оценить оценили качество описания экспериментальной энергии когезии и объема ячеек кристаллических решеток. Для всех полученных величин отклонение от экспериментальных значений было не более 10%, что является одним из самых лучших показателей для силового поля общего назначения, не подогнанного под конкретное вещество и конкретное агрегатное состояние.

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

Библиография Диссертация по биологии, кандидата физико-математических наук, Переяславец, Леонид Борисович, Пущино

1. Полозов, Р.В., Метод полуэмпирического силового поля в конформационном анализе биополимеров. 1981, Москва: "Наука". 123.

2. Rigby, D., Fluid density predictions using the COMPASS force field. Fluid Phase Equilibria, 2004. 217: p. 77-87.

3. Kaminski, G.,Jorgensen, W.L., Performance of the AMBER94, MMFF94, and OPLS-AA Force Fields for Modeling Organic Liqitids. J. Phys. Chem., 1996. 100: p. 1801018013.

4. Martin, M.G.,Thompson, A.P., Industrial property prediction using Towhee and LAMMPS. Fluid Phase Equilibria, 2004. 217(1): p. 105-110.

5. Martin, M.G., Comparison of the AMBER, CHARMM, COMPASS, GROMOS, OPLS, TraPPE and UFF force fields for prediction of vapor-liquid coexistence curves and liquid densities. Fluid Phase Equilibria, 2006. 248(1): p. 50-55.

6. Freddolino, P.L., Arkhipov, A.S., Larson, S.B., McPherson, A.,Schulten, K., Molecular dynamics simulations of the complete satellite tobacco mosaic virus. Structure, 2006.14: p. 437-449.

7. Gilson, M.K., Given, J.A., Bush, B.L.,McCammon, J.A., The statistical-thermodynamic basis for computation of binding affinities: A critical review. Biophys. J., 1997. 72: p. 1047-1069.

8. Shoichet, B.K., Virtual screening of chemical libraries. Nature, 2004. 432: p. 862-865.

9. Reddy, M.R., Erion, M.D.,Agarwal, A., Free energy calculations: Use and limitations in predicting ligand binding affinities. Rev. Comput. Chem., 2000. 16: p. 217-304.

10. Фиикельштейн, А.В.,Птицын, О.Б., Физика белка, 3-е издание. 3 ed. 2005, М.: Книжный дом "Университет". 376.

11. Krieger, E., Darden, Т., Nabuurs, S.B., Finkelstein, A.,Vriend, G., Making Optimal Use of Empirical Energy Functions: Force-Field Parameterization in Crystal Space. Proteins: Structure, Function and Bioinformatics, 2004. 57: p. 678-683.

12. Moult, J., Fidelis, K., Kryshtafovych, A., Rost, B.,Tramontano, A., Critical assessment of methods of protein structure prediction Round VIII. Proteins: Structure, Function and Bioinformatics, 2009. 77(S9): p. 1-4.

13. Vorobyov, I.V., Anisimov, V.M.,MacKerell, A.D., Jr., Polarizable Empirical Force Field for, Alkanes Based on the Classical Drude Oscillator Model. J. Phys. Chem. B, 2005. 109(40): p. 18988-18999.

14. Axilrod, B.M.,Teller, E., Interaction of the van der Waals Type Between Three Atoms. J. Chem. Phys., 1943. 11: p. 299.

15. Finkelstein, A.V., Average and extreme multi-atom Van der Waals interactions: Strong coupling of multi-atom Van der Waals interactions with covalent bonding. Chemistry Central Journal, 2007. 1: p. 21.

16. Finkelstein, A.V., Lobanov, M.Y., Dovidchenko, N.V.,Bogatyreva, N.S., Many-atom van der Waals interactions lead to direction-sensitive interactions of covalent bonds. Journal of Bioinformatics and Computational Biology, 2008. 6(4): p. 693-707.

17. Perera, F.P., Environment and Cancer: Who Are Susceptible? Science, 1997. 278: p. 1068-1073.

18. Chothia, C., Hydrophobic bonding and accessible surface area in proteins. Nature, 1974. 248: p. 338-339.

19. Still, W.C., Tempczyk, A., Hawley, R.C.JIendrickson, Т., Semianalytical Treatment of Solvation for Molecular Mechanics and Dynamics. J. Am. Chem. Soc., 1990. 112: p. 6127-6129.

20. Gallicchio, E., Zhang, L.Y.,Levy, R.M., The SGB/NP Hydration Free Energy Model Based on the Surface Generalized Born Solvent Reaction Field and Novel Nonpolar Hydration Free Energy Estimators. J. Comput. Chem, 2002. 23: p. 517-529.

21. Ponder, J.W.,Case, D.A., Force Fields for Protein Simulations. Adv. Prot. Chem., 2003. 66: p. 27-85.

22. Mitomo, D., Watanabe, Y.S., Kamiya, N.,Higo, J., Explicit and GB/SA solvents: Each with two different force fields in multicanonical conformational sampling of a 25-residuepolypeptide. Chem. Phys. Lett., 2006. 427: p. 399-403.

23. Levitt, M., Hirshberg, M., Sharon, R.,Daggett, V., Potential energy function and parameters for simulations of the molecular dynamics ofproteins and nucleic acids in solution. Comput. Phys. Commun., 1995. 91: p. 215-231.

24. Lee, F.S., Chu, Z.T.,Warshel, A., Microscopic and semimicroscopic calcidations of electrostatic energies in proteins by the POLARIS and ENZYMIXprograms. J. Сотр. Chem., 1993. 14: p. 161-185.

25. Vardi-Kilshtain, A., Roca, M.,Warshel, A., The Empirical Valence Bond as an Effective Strategy for Computer-Aided Enzyme Design. Biotechnol J., 2009. 4(4): p. 495-500.

26. Allen, M.P.,Tildesley, D.J., Computer Simulation of Liquids. 1991, Oxford: Oxford University Press.

27. McCammon, J., Gelin, G.B.,Karplus, M., Dynamics of folded proteins. Nature, 1977. 267: p. 585-590.

28. Jayachandran, G., Vishal, V.,Pande, V.S., Folding Simulations of the Villin Headpiece in All-Atom Detail. J. Chem. Phys., 2006. 124: p. 164902.

29. Elsen, E., Houston, M., Vishal, V., Darve, E., Hanrahan, P.,Pande, V. N-Body simulation on GPUs. Conference on High Performance Networking and Computing Proceedings of the 2006 ACMI/IEEE conference of Supercomputing. 2006. P. Tampa, Florida, USA.

30. Jones, J.E., On the Determinant of Molecidar Fields. II. From the Equation of State of a Gas. Proc. R. Soc. London, Ser. A, 1924. 106: p. 463-477.39.42.