Математическая модель процессов взаимодействия диоксида серы со структурными элементами клеточной мембраны

Официальные оппоненты: доктор физико-математических наук,

профессор

Элькин Михаил Давыдович

· математическая модель и алгоритм математического моделирования процессов взаимодействия диоксида серы со структурными элементами клеточной мембраны;

· база данных и программный модуль BioMolDiagrams для составления молекулярных диаграмм биополимеров мембраны;

· квантово-химические расчеты для установления геометрии молекул и сопоставление результатов с известными данными;

· программный модуль ModelInteractions для моделирования взаимодействий в макромолекулах;

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

Практическая значимость. На основании математической модели и алгоритма разработан программный комплекс Molecular Model (исх. № 17-4 от 13.07.2009), позволяющий:

· структурировать расчеты макромолекул в виде молекулярных диаграмм, которые необходимы химикам-органикам для синтеза новых соединений;

· формировать схемы взаимодействий, на основании которых появляется информация о потенциальных «мишенях» в биологических мембранах, необходимая химикам для установления антидотов к токсичным соединениям, а биоэкологам для изучения экологического последствия.

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

В соавторстве по данному направлению созданы и зарегистрированы базы данных:

· № 2009620009 «Токсическое воздействие на биологические структуры»;

· № 2009620103 «Воздействие блокаторов H-рецепторов и антигистаминных препаратов на биологические структуры»;

· № 2009620395 «Молекулярные диаграммы структурированных поверхностей»;

· исх. № 17-4 от 13.07.2009 «Структурные характеристики биополимеров и двухкомпонентных взаимодействующих систем», которые нашли свое применение при создании новой учебной дисциплины для вузов «Квантовая экологическая химия». Этот предмет изучают студенты химического факультета Астраханского государственного университета. В результате, по программе «У.М.Н.И.К. - 2009» получен грант.

Результаты диссертационной работы используются в учебном процессе по дисциплинам: «квантовая механика и квантовая химия», «экология на стыке математики, физики и химии».

Апробация работы. Результаты исследований доложены на различных Международных и Всероссийских научных конференциях, среди которых: «Средства и методы обеспечения экологической безопасности» (Астрахань, 2005) «Фундаментальные и прикладные проблемы современной химии в исследованиях молодых ученых» (Астрахань, 2006); III школа-семинар «Квантовохимические расчеты: структура и реакционная способность органических и неорганических молекул» (Иваново, 2007); «Фундаментальные и прикладные проблемы современной химии» (Астрахань, 2008); «Инновационные технологии в управлении, образовании, промышленности АСТИНТЕХ-2009» (Астрахань, 2009); IV школа-семинар «Квантово-химические расчеты: структура и реакционная способность органических и неорганических молекул» (Иваново, 2009).

Публикации. По материалам диссертации опубликовано 5 статей в журналах, рекомендованных ВАК и 5 статей в материалах конференций.

Структура и объём диссертации. Диссертация состоит из введения, трёх глав, заключения, приложения и библиографического списка (132 наименования). В приложении имеется 4 акта о внедрении результатов диссертационных исследований в учебный процесс.

СОДЕРЖАНИЕ РАБОТЫ

Во введении обоснована актуальность темы, научная новизна, практическая значимость, сформулированы цели и задачи исследования.

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

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

HШ(q) = EШ(q) (1)

где Ш(q) - волновая функция, удовлетворяющая вероятностному уравнению:

(2)

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

Наиболее ранняя модель, описывающая реакционную способность молекул, представлена в статическом методе [Коулсон и Лонге-Хиггинс, 1947 г.]. Расчеты сводятся к определению параметров р-электронной плотности (р_r) атомов, степени р-электронного связывания (p_rs), индексов свободной валентности (F_r) и значений самополяризуемости (П_rr) атомов:

(3)

(4)

(5)

(6)

где ц_занят - занятая молекулярная орбиталь, n_i - число электронов на i-ой орбитали, с_ir - орбитальный коэффициент, N_max - максимальное число связей атома r, N_r - сумма порядков связей атома r, е_i - энергия электрона i-ой орбитали.

Согласно методу, взаимодействию подвергается тот атом, который обладает максимумом р_r или максимумом П_rr (для электрофильных агентов), минимумом р_r или максимумом П_rr (для нуклеофильных агентов), максимумом F_r (для свободных радикалов).

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

Рассмотренная модель является базой, принципы которой заложены в модель, предложенную нами для описания взаимодействий в макромолекулах. На рисунке 1 приведена структура алгоритма моделирования взаимодействия диоксида серы с компонентами мембраны.

На этапе оптимизации исходных компонентов предлагается разбить процесс на две стадии. На первой стадии, для расчета полной энергии макромолекулы, можно использовать PM3 метод. Большинство расчетов в первом приближении, как правило, имеют разногласия с экспериментальными данными, поэтому задачу можно свести ко второй стадии моделирования. Для дальнейшей обработки z-матрицы и уточнения геометрии предлагается использовать расширенный базис 6-31G^* в приближении самосогласованного поля Хартри-Фока.

Выделяются следующие критерии оптимизации:

· критерий минимума полной энергии E или равенство градиента функции энергии нулю;

· критерии геометрических параметров (длина связи, валентный угол);

· критерий минимума энергии адсорбции ДE_адс.

Важной характеристикой регистрации стабильной формы является задание значения градиента функции энергии, что соответствует норме изменения i-ой координаты.

(7)

где суммирование производится по всем n атомам; x_i, y_i, z_i - декартовы координаты i - ого атома.

Задача геометрической оптимизации процесса сводится к поиску наиболее устойчивых структур, которым соответствует минимум энергии, поэтому она может быть решена методом следования собственному вектору. Выбор метода основан на рациональном использовании ресурсов ЭВМ и времени расчетов. Критерием окончания поиска точки экстремума на первой стадии является установление «завышенной» величины градиента 10^-2 кДж/(моль•Е), это позволяет выявить локальные стационарные точки, которым соответствуют метастабильные конфигурации. Дальнейшее снижение градиента до 10^-5 кДж/(моль•Е) позволяет скорректировать геометрию найденных конформационных систем с целью обнаружения критической точки. Рассчитанные оптимизированные геометрические показатели (длина связи, валентный угол) должны быть сопоставимы с экспериментальными или справочными данными.

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

На этапе расчёта структурных характеристик оптимизированных моделей следят за условными изменениями в электронной конфигурации по значениям зарядов на атомах (q_i, e), определяющих интегральное значение электронной плотности вблизи каждого атома. Об изменениях электронной плотности, при образовании системы из отдельных подсистем, судят по величине переноса заряда Дq с биополимера на диоксид серы:

(8)

Тепловые эффекты (ДE_адс, кДж/моль) формирования предполагаемых систем рассчитывали как:

(9)

По соотношениям величины ДE_адс и длины связи (r_ij) между атомами, взаимодействующих молекул, судят о стабильности образующихся систем. Перед формированием схем взаимодействий проводится анализ n-количества систем, из которых осуществляется отбор. Рассмотрим критерии оптимального выбора положений и условия их выполнения.

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

Если условие не выполняется то процесс переходит на проверку п. 3.

2. Если Д^/_i+1 = (ДE_min - ДE_i+1) ? 5 кДж/моль, то атом в (i+1)-системе также обозначается ().

Если условие не выполняется, то над атомом фиксируется ().

3. Если r_ij (1,821 - 1,900 Е), осуществляется поиск минимального значения ДE_адс. При выполнении условия, взаимодействующий атом в молекуле биополимера обозначается ().

4. Если Д^/_i+1 = (ДE_min - ДE_i+1) ? 3 кДж/моль, то атом в системе (i+1) фиксируется ().

Если условие не выполняется, то над атомом фиксируется ().

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

Во второй главе представлена реализация модели на ЭВМ, которая позволяет описать структуру и свойства макромолекул. Комплексы программ Mopac и РС Gamess реализуют только численные методы оптимизации геометрии молекул, а разработанная автоматизированная система Molecular Model позволяет генерировать результаты и составлять молекулярные диаграммы биополимеров (модуль BioMolDiagrams). Пользовательский интерфейс программы реализован в среде Microsoft Visual Studio 2008 на языке программирования C#.

После оптимизации макромолекул, из выходного текстового файла РС Gamess проводится выборка данных и внесение их в базу данных. Для построения использовалась система управления базами данных FireBird, версии 2.1. На рисунке 2 приведена схема связей данных.

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

Наименование исследуемых структур хранится в таблице type1, а данные о методе расчета в таблице method_analysis. Таблица structure_one содержит информацию об энергетических параметрах молекулы.

Величина теплоты образования и общая энергия автоматически заполняют базу данных только в том случае, если достигнута сходимость и найдена оптимальная конфигурация молекулы. При удачном окончании процесса оптимизации в программе Gamess, текстовый документ содержит запись “Equilibrium Geometry Located”. Если структура не достигла оптимального положения, то процесс запускается повторно с новыми координатами.

Для характеристики функциональных групп в таблице group_one приведены поля их графического расположения в системе координат xy. Информация о параметрах зарядов на атомах (q_i) и значения электронной плотности (р_i) хранится в таблице atom_one. Длина химических связей (r_ij) между атомами в исследуемой структуре хранится в таблице connection.

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

Входными параметрами модуля BioMolDiagrams являются:

· вычисленные геометрические параметры молекулы в виде z-матрицы для конкретного метода расчета по программе Mopac;

· вычисленные энергетические и зарядовые характеристики молекулы в программе Mopac;

· квантово-химический метод расчета;

· название, формула и тип соединения.

· Выходными данными модуля являются:

· таблица оптимизированных геометрических параметров молекулы;

· минимизированные энергетические параметры молекулы;

· зарядовые характеристики молекулы;

· двухмерная молекулярная модель.

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

В качестве белкового компонента выбрана молекула пентапептида, содержащая различные группировки (рис. 3).

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

Таблица 1

Сравнение результатов расчета геометрии индивидуальных молекул и модели пентапептида методами PM3 и HF/6-31G^*.

Из таблицы видно, что PM3 метод хорошо воспроизводит параметры С?С, С_ар?O, С=O, С?H и С?S связей. Намного сложнее метод описывает связи с гетероатомами, сопряженными с другими функциональными группами, например, длина связи N?C в пептидной группе завышена на 0,1 Е, а валентный угол, включающий атом азота или атом кислорода, как правило, оказывается, занижен на 5-9^є. Расширенный базис в хартри-фоковском приближении намного лучше воспроизводит длину связей и валентные углы с гетероатомами, но напротив не рассчитан на описание серосодержащих и ароматических соединений (табл. 1).

Аналогично проведены расчеты для липидной составляющей - структуры моноглицерида, диглицерида, триглицерида и углеводной составляющей - олигосахарид, лактоза. Сопоставление результатов, позволило установить правомерность использования PM3 метода и расширенного базиса 6-31G^* в хартри-фоковском приближении. Однако полуэмпирический метод, в отличие от расширенного базиса, справляется с поставленной задачей в течение нескольких минут. В таблице 2 приведены данные по затрате времени для установления равновесной конфигурации биополимеров мембраны разными подходами.

Таблица 2

Техническая характеристика применимости методов

Из таблицы 2 следует, что использование результатов расчета (6-31G^*) для формирования схем молекулярных диаграмм может быть ограничено только по времени. В результате, для всех биополимеров в программном модуле BioMolDiagrams составлены молекулярные диаграммы, которые могут быть представлены в виде наглядных иллюстраций (рис. 4).

токсикант молекула математический моделирование

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

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

Энергетические параметры образующейся адсорбционной системы приведены в таблице interaction. Параметры электронной плотности, зарядовые характеристики атомов и информация об их реакционной способности хранится в таблице atom_two. Все параметры получены из текстового файла программы Gamess. Таблица connection_two содержит информацию о длине образующейся связи при взаимодействии молекул.

Принцип подготовки входных данных при работе с адсорбционными системами остается таким же, как и при работе с одной структурой. Программный модуль ModelInteractions реализует расчёт энергии формирования данных систем (ДE_адс) и величины переноса зарядов (Дq) на диоксиде серы, а анализ параметров позволяет оценить активность атомов в макромолекулах и графически сформировать схемы взаимодействия. В результате, статическая модель включает специальные обозначения активности атомов (, ).

Входными параметрами модуля ModelInteractions являются:

· геометрические параметры системы в виде z-матрицы для конкретного метода расчета по программе Mopac;

· вычисленные энергетические и зарядовые характеристики системы в программе Mopac;

· квантово-химический метод расчета;

· название и формула соединения;

· длина связи между атомами взаимодействующих молекул.

Выходными данными модуля являются:

· таблица оптимизированных геометрических параметров системы;

· минимизированные энергетические параметры системы;

· зарядовые характеристики атомов в системе;

· характеристика активности атомов;

· двухмерная молекулярная модель взаимодействия.

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

Рассмотрим модели адсорбционных систем и схемы взаимодействий пентапептида (цистеинил-фенилаланил-аланил-гомоцистеинил-тирозин). Для уменьшения загромождения в системы введены обозначения R, R^', обозначающие остаток углеродного скелета.

В таблице 3 приведена сравнительная характеристика геометрических и энергетических параметров изучаемых адсорбционных систем 1?3, полученных PM3 методом.

Таблица 3

Параметры адсорбционных систем «пентапептид?SO₂»

Полученные результаты позволяют оценить энергию взаимодействия, а так же выяснить наиболее активный центр пентапептида, способный легко подвергаться направленному влиянию со стороны SO₂. В первой системе атака происходит по пептидной группе, устойчивость определяется за счет образования шестичленного цикла. Минимальная величина переноса заряда в системе 1 (-0,210 e) может быть вызвана сильной поляризацией связи >S=O, за счет увеличения на кислороде O¹ электронной плотности q(-0,6474)/q(-0,6979). В системе 2 гидроксогруппа открыта для взаимодействия, поэтому наблюдается перенос заряда на H⁴⁵ q(0,2013)/q(0,2483). Система 3 образована за счет открытого метиленового радикала.

При исследовании 35-ти моделей адсорбционных систем, было отмечено, что в среднем энергия адсорбции диоксида серы на пентапептиде может составить -21,00 кДж/моль. В результате расчетов и выполнения условий выбора оптимальных положений для пентапептида в программном модуле ModelInteractions составлена схема взаимодействий (рис. 6).