CC BY
12Математические заметки СВФУ Июль—сентябрь, 2014. Том 21, № 3
УДК 517.926; 517.27
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ЗАДАЧИ МИНИМИЗАЦИИ РАСХОДА РЕСУРСОВ ДЛЯ ЛИНЕЙНЫХ СИСТЕМ С ПОСТОЯННЫМИ ЗАПАЗДЫВАНИЯМИ В СОСТОЯНИИ И УПРАВЛЕНИИ Г. В. Шевченко
Аннотация. Предлагается численный метод решения задачи минимизации расхода ресурсов для линейных систем с постоянным временным запаздыванием как в фазовых состояниях системы, так и в управлении. Этот метод является обобщением метода решения указанной задачи с запаздыванием только в фазовом состоянии системы. Доказана глобальная сходимость предлагаемого метода к е-оптимальному решению. Под е-оптимальным решением понимается допустимое управление и(Ь), Ь £ [0, Т], переводящее систему в ^-окрестность начала координат и доставляющее функционалу задачи значение, отличное от оптимального не более чем на е.
Ключевые слова: дифференциальное уравнение с запаздыванием, функционал, оптимальное управление, численный метод.
Пусть управляемый объект описывается системой линейных обыкновенных дифференциальных уравнений с постоянными запаздываниями Н1 > 0 и > 0 в фазовом состоянии и управлении соответственно вида
х(г) = Л(г)х(г) + с(г)х(г - + в(г)и(г) + Б(г)и(г — н2), г е [0, т],
х(г) = <^(г), г е [—0], (1)
где х — фазовый вектор состояния объекта, <^(г) — заданная непрерывная функция, Л(г), С (г) и В (г) — непрерывные матрицы размеров п х п, п х п и п х в соответственно, и(г) — измеримое управление, стесненное ограничением
Ы*)1<1, з=~з, *е[о,г], (2)
и ч(Ь) = 0 на интервале [—¡г2, 0].
Предполагается, что система (1) переводима при заданной на [—Л-1, 0] функции ^ и нулевом управлении на [—¡г2, 0] в начало координат допустимыми управлениями.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 13—01—00329) и Сибирского отделения Российской академии наук (проект №80).
© 2014 Шевченко Г. В.
Задача 1. Требуется найти допустимое управление и0(Ь), Ь € [0,Т], переводящее систему (1) за фиксированное время Т в нулевое конечное состояние ж(Т) = 0 и минимизирующее функционал
Т 8
^(«)=/]Т а,- К (Ь)| <И, (3)
0 ¿=1
где аз = 0, а, > 0. ¿=1
Очевидно, что задача 1
имеет решение только при Т > Топт, где Топт время оптимального по быстродействию перевода системы (1) в начало координат. Естественно предположить, что Т > Топт.
Введем согласно принципу максимума Понтрягина в задаче с запаздыванием [1, гл. 4, § 27; 2] сопряженную систему
. = Г -А*(Ь)^(Ь) - С * (Ь + ^(Ь + Ь € [0,Т - ы (4)
= 1 -А*(Ь)^, 4 € [Т - ЛЬТ], ()
и выпишем функцию Понтрягина вышепоставленной задачи:
Я (^(Ь),ж(Ь),ж(Ь - - = аз К (Ь)1
¿=1
+ (^>(Ь), А(Ь)ж(Ь) + С(Ь)ж(Ь - ^)> + (^(Ь), В(Ь)и(Ь) + - ^)>, (5)
где (•, •) — скалярное произведение векторов.
Тогда для оптимальности управления и0(Ь) и траектории ж0(Ь), Ь € [0,Т], необходимо существование ненулевой вектор-функции ^0(Ь), являющейся решением сопряженной системы (4) при некотором вполне определенном граничном условии ^(Т) = с0, такой, что
Я(^0(Ь),ж0(Ь),ж0(Ь - ^1),и0(Ь),и0(Ь - Л2))
тах[Я(Ь), ж0(Ь), ж0(Ь - и, и0(Ь - + Я+ иеи (6) ж0(Ь + Л2), ж0(Ь - 1ц + Л2), и0(Ь + Л2), и)], Ь € [0, Т - Л2],
тах Я(^0(Ь), ж0(Ь), ж0(Ь - ^), и, и0(Ь - Ь € [Т - Т],
иеи
где и = {и | < 1, ] = 1, в}. Отметим, что ж°(£) = <р(Ь) при I £ [—/11, 0],
и0(Ь) = 0 при Ь € [-^2, 0].
Из (6) и (5) следует, что
и0(Ь)
агг тах
иеи
агг тах
иеи
-2 £ а и | + (-0°(Ь), В(Ь)и> + + Л2), + ^)и>
¿=1
Ь € [0,Т - Л2], Ь € [Т - ^2,Т].
-£ а,-кi + (^0(Ь),В(Ь)и> ¿=1
Последнее соотношение можно записать в виде
-1, ВД < -2а,,
-Чоч < < з=Т^, [0,Т-/г2], (7) 1, (Ь) > 2а,,
-1, _
щ{1) = { 0, -<х,- < {ф°{1),Вз{1)) < а^ з = 1, I € [Т-к2,Т], (8)
1, (¿)) > а.з,
гдеФ^) = + В ¿{г) и£>^+/г2) -^'-е столбцы
матриц -В(£) и + /12) соответственно, j = 1, в, — решение сопряженной системы (4) с граничным условием
^(Т )= с0. (9)
Здесь с0 — некоторый ненулевой вектор из М".
В силу однородности системы (4) в дальнейшем можно ограничиться рассмотрением только граничных условий (9) с единичными нормами ||с|| = 1, заменив (7) и (8) выражениями
Г -1, ВД < -2уа„
| -2ус^- < < 2 ус^-, 3 = 1, I € [0,Т- к2], (Ю)
[ 1, (¿) > 2уа^-,
Г -1, (Ф0(1),В^)) _
= I о, -ус^- < {ф°{1),Вз{1)) < уо^-, 3 = 1,ге [т-н2,т], (и)
[ 1, (¿)) >уа^-,
где у — неотрицательное действительное число.
Пусть Б с М" — единичная сфера с центром в начале координат. Обозначим через с, у) управление, компоненты которого удовлетворяют (10) и (11) при некоторых векторе с € Б и неотрицательном числе у. Из (10) и (11) следует, что если а2 > 0, то при
имеет место тождество щ (4, с, у) = 0, 4 € [0, Т], а при
у > у(с) = тах
все компоненты вектор-функции с, у) такие, что соответствующие а2- > 0, тождественно равны нулю. Рассмотрим функцию
Г s
у) = у)) = / ^^aj|uj(t, с, у)| dt
0 j=i
при фиксированном c из S, где u(c, у) = u(t, с, у), t g [0,T]. Функция у) на интервале [0,у(с)] непрерывна и монотонно убывает по у в силу (10) и (11). Более того, существует единственное число у(с), 0 < у(с) < у (с), такое, что на интервале [0, у(с)] функция у) постоянна и равна
у -"Г а-т
^ max — / J 5
j=i
а на интервале [//(с),/(с)] строго убывает по /. Следовательно, при любых фиксированном с из Б и положительном ^ < /шах существует единственное число (с) € [//(с),/(с)] такое, что 5^(с, (с)) =
Обозначим через О(^) множество правых концов траекторий системы (1) при всевозможных допустимых управлениях, доставляющих функционалу (3) значение, меньшее или равное , т. е.
О(^) = {ж = ж(Т,и) | и(Ь) € и, Ь € [0,Т], ^(и) <
где ж(Т, и) — решение задачи (1) при допустимом управлении и в момент времени Ь = Т. Границей </") является множество
= {ж = ж(Т,и^(с)) 1 (с) = и^(С)(Ь) = и(Ь,с,/^г(с)) Ь € [0,T], с € 5}.
Множество О(^тах) совпадает с областью достижимости ^(Т) системы (1), поскольку их границы совпадают и они телесны.
Так как функция (с, /) строго убывает по / на отрезке [//(с), /(с)] при любом с € Б, справедливы включения
О(^) с О(^) при любых ^ > Л > 0. (12)
Из полученных включений и телесности ^(Т) следует, что О(^) — тело при любом 0 < ^ < ^тах. Кроме того, в силу (12) существует наименьшее Д^п такое, что начало координат лежит на границе множества О(^"тт), а при любом ^ € (^тт, </"тах] — внутри О(^). Это значит, что существуют с* € Б и //* такие, что ж(Т, и(с*,//*)) = 0.
Поиск таких с* € Б и //* предлагается проводить с использованием симплексов из симплексных покрытий множеств О(^) [3]. Описание и обоснование предлагаемого метода требует введения некоторых понятий.
Пусть г1,..., 2П+1 € К" — различные точки такие, что их линейная выпуклая оболочка а = [г1,..., гп+1] является телом в К". Будем называть множество а п-мерным симплексом с вершинами г1,.. ., Два п-мерных симплекса а1 и а2 называются смежными, если у них п общих вершин и их пересечение есть (п - 1)-мерный симплекс. Из определения следует, что пересечение является общей гранью максимальной размерности.
Пусть В с К" — компактное тело, а0 = [г0, ..., 2"+1] — п-мерный симплекс с вершинами, лежащими на границе В. По каждой его грани максимальной размерности ст® = ..., 1, + ] = 1,71 + 1, строим смежный
1 ¿—1 ~7 7+1 "+1
ему симплекс с вершинами -0,..., -0 , г7, -0 ,..., ^ , у которого «новая»
вершина г7 является граничной точкой В, максимально удалена от гиперплоскости, проходящей через остальные вершины, и расположена по разные стороны с точкой -0 относительно этой гиперплоскости.
Назовем построенные симплексы симплексами первого слоя, а симплекс а0 будем считать симплексом нулевого слоя. Затем для каждого симплекса первого слоя строим по его (п - 1)-мерным граням, которые не являются общими с (п - 1)-мерными гранями симплекса а0, по той же схеме смежные симплексы. Построенные симплексы составят второй слой. Ясно, что во втором слое содержится ровно п(п + 1) симплексов. Аналогично для каждого симплекса к-го слоя (к > 2) строятся смежные ему симплексы (к + 1)-го слоя.
Обозначим через объединение всех симплексов k-го слоя. Ясно, что в
оо
k-м слое будет nk-1(n + 1) симплексов. Назовем U выпуклым покрытием
k=0
множества B и обозначим его через П®.
По построению в силу компактности множества 58 имеем со 58 = П®, где со «В — выпуклая оболочка множества 58, D — замыкание множества D.
Таким образом, внутренность множества B содержится в выпуклом покрытии П®, a его граница содержится в границе замыкания покрытия П®. (В дальнейшем, говоря о покрытии n-мерными симплексами, будем подразумевать построенное покрытие.) Тем самым имеет место
Теорема 1 (о покрытии). Выпуклая оболочка любого компактного тела совпадает с замыканием выпуклого покрытия тела.
Следствие 1. Пусть B — компактное тело в Rn и z0 g int B. Тогда в любом покрытии П® тела B имеются конечное ко > 0 и n-мерный симплекс о g &k0 такие, что z0 g о.
Легко видеть, что симплексное покрытие любого компактного тела B из Rn определяется выбором (единственного) симплекса нулевого слоя. Не ограничивая общности, будем считать, что этот симплекс о0 = [z1,..., zn+1] построен по следующей схеме.
Пусть с1 — произвольная единичная точка из Rn. В качестве первой вершины симплекса о0 выбирается точка z1 из B такая, что (c^z1) = max(с1, ж).
же®
Очевидно, что z1 является граничной точкой множества B.
Вершины z3, j = 2, n + 1, выбираются так. Вначале находится с3 — решение системы линейных алгебраических уравнений
{c,zl) = -1, i=TJ^l, (13)
которое нормируется. Затем в качестве j-й вершины выбирается точка zj такая, что
{с3, z3) = max(cJ,x), j = 2,n+l. (14)
же®
Очевидно, что z3, j = 2,n + 1, являются граничными точками множества 58.
При таком построении симплекса о0 остается некоторый произвол в выборе векторов с3, j = 2,п, поскольку системы уравнений (13) при j < n + 1 недо-определены и имеют не единственное решение. Исключить произвол в этом выборе можно следующим образом. Пусть Z — квадратная матрица размера j — 1, составленная из первых j — 1 линейно независимых столбцов матрицы системы (13), а ж = (ж1,... ,xj-1) — решение системы линейных алгебраических уравнений ZxT = bT, где b = (—1,..., —1). Тогда берем cj = (cj,..., с^) такое, что
{Xi, если ¿-й столбец матрицы системы (13)
входит в матрицу Z, i = 1 ,п, 0 в противном случае,
и нормируем его.
При отсутствии произвола в выборе решения системы (13) не исключается произвол в выборе вершин симплекса нулевого слоя, поскольку в общем случае в силу невыпуклости множества B максимум в (14) может достигаться в нескольких точках. Если это так, то среди них выбираем в качестве j-й вершины точку с минимальной нормой.
Область достижимости R(T) (или, что одно и то же, 0(^тах)) является компактным телом в Rn. Следовательно, к ней применима теорема о покрытии и ее следствие. Кроме того, в силу линейности правой части системы (1) по фазовым состояниям и по управлению R(T) — строго выпуклое компактное тело.
В разд. 2 описаны метод и вычислительный алгоритм решения поставленной в разд. 1 задачи управления.
2. Метод и вычислительный алгоритм решения
Предлагаемый метод решения основан на построении последовательности симплексов {<rk} с вершинами, лежащими на границах множеств О(^). Перед полным описанием вычислительной схемы метода дадим его краткое формальное описание.
Первый этап. Строится такая последовательность смежных симплексов {<rk}, вершины которых являются граничными точками R(T), что р(<тк) > p(CTk+i) для любого k > 0, где р(<т) — расстояние от начала координат до симплекса <г. Отметим, что в силу следствия эта последовательность конечна.
Симплекс <г0 строится по вышеприведенной схеме, где в качестве тела B взято множество достижимости R(T). Построение заканчивается, когда очередной симплекс ak° поглотит начало координат. Симплекс ak, k > 1, строится следующим образом. Пусть точка z* е такова, что ||z*|| = min ||z||. Так
zEo-k-1
как 0 е <rk-1 (в противном случае построение закончилось бы при ko = k — 1), точка z* лежит на границе симплекса <rk-1 в некоторой его грани максимальной размерности. Пусть эта грань содержит вершины z1,..., zi0-1, zi0+1,..., zn+1. Эти вершины будут n первыми вершинами симплекса ak. Находим решение с системы линейных алгебраических уравнений
(c, zl) = —1, i = 1,. .., i0 — 1, i0 + 1,.. . ,n +1,
и нормируем его. В качестве (n + 1)-й вершины берется конец траектории Х(Т) системы (1) при управлении
u(t) = argmax(c,x(t)), t е [0,T]. (15)
и£ U
При таком управлении справедливо равенство (с, X(T)) = max (с, ж).
x RR (Т)
Пусть zl — вершина симплекса ak°, U и сг — управление и сопряженные параметры, соответствующие этой вершине, г = 1, п + 1, и пусть ..., + 1 — решение следующей системы линейных алгебраических уравнений:
n+1 n+1
]ГА^г = 0, ]ГАг = 1. (16)
i=1 i=1
max •
Так как 0 £ ак°, имеем > О для всех г = 1, п + 1. Положим = ¿Р^
Второй этап. Если Л = {г | Л0 = 0} = 0, то переходим к шагу 2.
Шаг 1. Находим 7, 0 < 7 < 1, такое, что 0 € [х(Т, 7м1),..., х(Т, 7и"+1)] (сг, х(Т, 7иг)} > 0, и полагаем
п+1
£ А°сг * ¿=1
n+1
£ A0ci
i=1
zl = x(T,jul), i = l,n + l.
Затем находим решение A°,..., АП+1 системы уравнений (13) и переходим к шагу 2.
Шаг 2. Полагаем ¿0 := min i, находим решение с системы линейных ал-
ieA
гебраических уравнений (с, Ж + 0,5(zi — Ж)} = —1, i = ¿0,1 < i < n +1, где Ж — правый конец траектории свободного движения системы (1) (при управлении и = u(t) = 0, t g [0, T]), и нормируем его. Если (с, zi0} < 0, то меняем знак у вектора с. Полагаем с* := си переходим к следующему шагу.
Шаг 3. Находим число «0(с*) > 0 такое, что ^(и(с*, «0(с*))) = J.
Пусть z* = ж(Т,и(с*,«0(с*))). Если ||z*|| < е, где е > 0 — необходимая точность попадания в начало координат, то процесс вычислений заканчивается. Полученное управление и(с*,«0(с*)) = u(t, с*, «0(с*)), t g [0,T], приближенно оптимально, а ^"(и(с*,«0(с*))) — приближенно оптимальное значение функционала (3).
Если ||.z*|| > £, то среди точек zl, г = 1, п + 1, выбираем п точек z11,..., zln таких, что симплексу о* = [zi1,..., zin, z*] принадлежит начало координат. Точки zi1,..., zin, z* и соответствующие им параметры сп ,..., с^, с^+1 = с* перенумеровываются по порядку. Находим решение (A°,..., АП+1) системы (16). Если |Л| = 1, то переходим к шагу 1, иначе — к шагу 2.
Перед рассмотрением описания вычислительного алгоритма решения задачи, поставленной в разд. 1, введем следующие обозначения:
с(к) — k-е приближение оптимальных значений граничных условий для сопряженной системы (4);
— решение сопряженной системы (4) с граничным условием ^(T) = с(к);
k — k-е приближение оптимального значения J^m функционала (3);
1 2
у«1, у2 — нижняя и верхняя границы локализации решения у уравнения
^■(и(с* ,«)) = J»k; (17)
к — номер итерации.
Вычислительный алгоритм решения указанной выше задачи с точностью до обозначений совпадает с описанием вычислительного алгоритма решения задачи минимизации расхода ресурсов с запаздыванием только в фазовом состоянии системы [4]. Доказательство сходимости аналогично доказательству сходимости в [4].
и
ЛИТЕРАТУРА
1. Понтрягин Л. С., Болтянский В. Г., Гамкрелидзе Р. В., Мищенко Е. Ф. Математическая теория оптимальных процессов. 4-е изд. М.: Наука, 1983.
2. Боков Г. В. Принцип максимума Понтрягина в задаче с временным запаздыванием // Фундамент. и прикл. математика. 2009. Т. 15, № 5. С. 3—19.
3. Шевченко Г. В. Метод нахождения оптимального по минимуму расхода ресурсов управления для нелинейных стационарных систем // Автоматика и телемеханика. 2009. Т. 70, № 4. С. 119-130.
4. Шевченко Г. В. Метод нахождения оптимального по минимуму расхода ресурсов управления для линейных стационарных систем с постоянным запаздыванием // Автоматика и телемеханика. 2014. Т. 75, № 10. С. 25-38.
Статья поступила 15 ноября 2014 г. Шевченко Геннадий Васильевич
Институт математики им. С. Л. Соболева СО РАН, пр. Академика Коптюга, 4, Новосибирск 630090 зЬеусЬФша^. пбс . ги
Математические заметки СВФУ Июль—сентябрь, 2014. Том 21, №3
УДК 51.76
ОЦЕНКА ПРОФИЛЕЙ ПОТЕНЦИАЛА СРЕДНЕЙ СИЛЫ МЕТОДОМ UMBRELLA SAMPLING ДЛЯ ТРАНСМЕМБРАННОГО ПЕРЕНОСА МОЛЕКУЛЫ ВОДЫ
М. Ю. Антонов, Т. В. Науменкова, А. В. Попинако, И. Н. Николаев, К. В. Шайтан
Аннотация. Моделирование пассивного трансмембранного транспорта малых молекул методами молекулярного моделирования сопряжено с рядом сложностей по той причине, что за доступное время моделирования практически невозможно наблюдать процесс спонтанной диффузии и измерять усредненные макроскопические характеристики процесса. В связи с этим интерес вызывают подходы и методы, используя которые, можно за достижимое время получать результат, позволяющий проводить сравнительное изучение или сравнение с натурными экспериментальными данными. В данной работе с помощью методов управляемой молекулярной динамики и метода Umbrella sampling изучен процесс трансмембранного транспорта молекул воды через липидный бислой. Использован бислой 1,2-димиристоил^п-глицеро-3-фосфатидилхолина (ДМФХ) в силовом поле OPLS-AA. Профили средней силы построены из данных расчетов с длинами траекторий 6—10 нс. Произведено сравнение получаемых профилей как при различном пространственном положении пробной (проникающей через бислой) молекулы воды, так и при разном шаге выбора стартовых конфигураций, а также при различных начальных значениях скоростей атомов в системе. Общее время моделирования составило более 10 мс. Показана значительная зависимость результатов расчетов от начальных данных и протокола вычислительного эксперимента. В рамках подхода произведены оценка кинетических параметров проникновения и оценка профиля потенциала средней силы, а также изучена зависимость вычисляемых значений от начальных данных эксперимента. Оценена погрешность вычислений на основе анализа данных ряда численных экспериментов.
Ключевые слова: молекулярное моделирование, метод зонтичной выборки, биомембраны.
Одной из основных функций клеточной мембраны является выполнение барьерной функции — управляемого обмена веществом между клеткой и межклеточной средой. Структурную основу биологических мембран составляют бислои липидов, при этом свойство полупроницаемости липидного бислоя для
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 12-04-31934\13), Научно-образовательного фонда поддержки молодых ученых Республики Саха (Якутия) (грант 2014-01-0008) и Российского научного фонда (код проекта 14-24-00172).
© 2014 Антонов М. Ю., Науменкова Т. В., Попинако А. В., Николаев И. Н., Шайтан К. В.
некоторых молекул, таких как молекулы воды, является весьма важным для выполнения мембраной барьерной функции. Методы молекулярного моделирования достаточно давно доказали свою эффективность при изучении систем, имитирующих биологические, в связи с чем представляет интерес их использование для изучения процессов диффузии в липидных бислоях [1].
При экспериментальном изучении проницаемости липидных бислоев измеряется коэффициент проницаемости P (см/сек) [2,3]. Однако компьютерное моделирование пассивного трансмембранного транспорта малых молекул методами молекулярного моделирования сопряжено с рядом сложностей по той причине, что за доступное время моделирования практически невозможно наблюдать процесс спонтанной диффузии и измерять усредненные макроскопические характеристики процесса. В связи с этим интерес вызывают подходы и методы, используя которые, можно за достижимое время получать результат, позволяющий проводить сравнительное изучение или сравнение с натурными экспериментальными данными.
Достаточно широко распространенной методикой является изучение трансмембранного транспорта с использованием методов управляемой молекулярной динамики (УМД). В рамках данного подхода к системе прикладывается дополнительный потенциал, позволяющий стимулировать систему по выбранным степеням свободы. В частности, метод позволяет получать значения «эффективной» микровязкости бислоя, а также оценку коэффициента диффузии, которые можно использовать для сравнительного анализа. Недостатком этого подхода является наблюдаемая зависимость результатов от протокола вычислительного эксперимента, а также отсутствие оценки коэффициента разделения растворитель-мембрана [4,5].
Другим подходом является использование метода зонтичной выборки (Umbrella Sampling, US). Использование этого метода позволяет проводить статистический анализ высокоэнергетических областей конфигурационного пространства, которые практически не заполнены в равновесных МД экспериментах. В данной работе с помощью методов управляемой молекулярной динамики и метода US изучен процесс трансмембранного транспорта молекул воды через липидный бислой. Использован бислой 1,2-димиристоил-8п-глицеро-3-фосфатидилхолина (ДМФХ) в силовом поле OPLS-AA. Произведено сравнение получаемых профилей как при различном пространственном положении пробной (проникающей через бислой) молекулы воды, так и при разном шаге выбора стартовых конфигураций, а также при различных начальных значениях скоростей атомов в системе.
1. Материалы и методы
Суть метода молекулярной динамики сводится к представлению изучаемой системы как множества материальных точек, взаимодействие которых описывается классическими законами механики. Каждая точка при этом представляет собой атом или группу атомов изучаемой системы. На рис. 1 показан вид изучаемой системы, а также вид молекулы ДМФХ и ее химическое строение.
Взаимодействие в системе описывается в виде суммы частичных потенци-
Рис. 1. Вид изучаемой системы вид и строение молекулы ДМФХ (Ь).
алов:
и (Г) = 53 Ц3 (Ь,,, ) + £ )
+ £ 3 (ф<>з->к>,) + £ (гг>з) + 3 (Г4,3-
г,з,к,1 г=з
(1)
где первые три слагаемых представляют взаимодействия между химически связанными атомами (энергии валентных связей, валентных углов и торсионных углов), а последние два представляют существующие между парами валентно не связанных атомов силы — электростатические и Ван-дер-Ваальсовы.
Энергии валентных связей и валентных углов, как правило, описываются в виде гармонического потенциала с заданной константой жесткости:
(2)
величина константы жесткости выбирается таким образом, чтобы воспроизводить геометрию и колебательные частоты модельных соединений.
Торсионный угол четырех последовательно связанных атомов представляет собой двугранный угол между плоскостями, определяемыми крайними тройками атомов. Вид данного потенциала описывается соотношением
3 (Фг,з,м) = £ 4Ып(1 + осз(пф - ф0)).
(3)
В методе МД электроны не учитываются явно, каждый атом обладает парциальным зарядом, а кулоновская энергия взаимодействия между атомами опи сывается стандартным выражением:
№з
Цс (Гг,з ) =
4пеоГгз'
(4)
парциальные заряды атомов при этом подбираются таким образом, чтобы лучше описывать распределение электронной плотности в системе.
Ван-дер-Ваальсовы взаимодействия представляют собой силы отталкивания и притяжения электродинамической природы, описываемые потенциалом Леннард — Джонса следующего вида:
ЧШ -ШУ <«
В целом общий вид функции потенциальной энергии практически не изменяется в разных вариантах силовых полей, однако параметризация сил может различаться. В данной работе использовалось силовое поле ОРЬ8-ЛЛ [6, 7].
Таким образом, движение атомов в потенциальном поле описывается системой обыкновенных дифференциальных уравнений второго порядка:
= — (6) \дхдуг дгг)'
Данная система уравнений решается численно. Могут применяться различные численные методы решения уравнений, различающиеся точностью и вычислительной сложностью. Компромиссом между точностью процедуры и скоростью вычислений является широко используемый в МД метод Верле [1] и его варианты, при этом значения координат атомов на новом временном слое вычисляются по следующим формулам:
п{1, + Аг) = 2п(г) - п{1, - М) +
тг (7)
Применение схем более высокой точности в методе МД нецелесообразно в связи с высокой вычислительной сложностью, которая не может быть компенсирована использованием большей величины шага интегрирования. Помимо этого схема Верле обладает свойством сохранения некоторых интегралов движения и, в целом, лучше сохраняет энергию системы.
Исследование систем в ансамблях с постоянной энергией, как правило, не представляет значительного интереса, поскольку большинство экспериментальных данных получено в условиях постоянства температуры и/или давления. Поэтому для поддержания условий постоянства температуры и давления в МД экспериментах используются специальные алгоритмы: термостаты и баростаты. В данной работе использовался термостат стохастической динамики. Его суть заключается в том, что к атомам системы прилагается дополнительная случайная (стохастическая) сила, описывающая взаимодействие системы с внешней средой заданной температуры [8]. Для контроля условий постоянства давления использовался баростат Берендсена [9]. В рамках данного подхода объем системы может изменяться: к основному уравнению движения добавляются уравнения движения границ периодической ячейки.
В экспериментальных работах при исследовании пассивного трансмембранного транспорта малых молекул удается измерить коэффициент проницаемости мембраны Р, а также коэффициент распределения вода-гидрофобная среда Кр [2, 3,10,11]. Эти величины связаны соотношением
6
где D — коэффициент диффузии в мембране, Kp — коэффициент межфазного разделения вещества, Дж — толщина мембраны. В МД удается различными способами прямо измерить коэффициент диффузии D [4], однако измерить прямыми способами коэффициент межфазного распределения Kp не удается.
Для описания перехода системы между состояниями удобно ввести координату реакции позволяющую описать непрерывно изменение системы. В этом случае для свободной энергии Гиббса G справедливо соотношение
G(£) = -kT • ln(p(0), (9)
где p(£) — плотность вероятности нахождения системы в состоянии c заданным значением При моделировании трансмембранного транспорта в качестве как правило, выбирается координата z центра массы малой молекулы, где ось Oz перпендикулярна латеральной плоскости мембраны (рис. 2a).
Таким образом, из равновесной МД-траектории исследуемой системы можно получить вид профиля свободной энергии G, оценив величину плотности вероятности при помощи статистических методов. Проблемой является то, что за доступное время моделирования практически невозможно наблюдать спонтанную диффузию малых молекул через липидный бислой и невозможно набрать адекватную статистику по вероятности нахождения в той или иной фазе, поскольку регионы конформационного пространства с высокими значениями энергии G не будут представлены в выборке или будут представлены недостаточно.
Выходом здесь является метод зонтичной выборки, позволяющий за достижимое время моделирования оценить значение коэффициента распределения вода-мембрана. В рамках данного подхода к системе прилагается дополнительный потенциал, способный удержать исследуемую молекулу внутри потенциально невыгодной области [12-14] (рис. 2b). Модификация исходной функции потенциальной энергии заключается в добавлении внешнего потенциала — как правило, используется гармонический потенциал, приложенный к координате реакции:
U(r) = U(r) - W(r), где W(r) = kw(£(r) - £o)2. (10)
Поскольку вид функции W (r) известен, значение свободной энергии G для невозмущенной системы может быть получено по формуле
G(£) = -kT • ln(p(£)) - W(£) + const, (11)
где p(£) — плотность вероятности нахождения системы в состоянии c заданным значением £ в возмущенной системе.
Поскольку в рамках такого подхода возможно достоверно восстановить вид профиля свободной энергии только в небольшой окрестности £o, выбирается множество начальных точек вдоль координаты реакции таким образом, чтобы функции распределения системы вокруг каждой из начальных точек перекрывались функциями распределения в соседних точках.
Исходный вид профиля свободной энергии может быть восстановлен из множества полученных распределений вокруг каждой начальной точки одним из методов. В данной работе использовался метод анализа взвешенных гистограмм (WHAM, [15]).
(а) (Ь)
Рис. 2. Ось координаты реакции (а), вид гармонических потенциалов для удержания молекулы воды внутри мембраны (области высоких значений свободной энергии) (Ь).
2. Протокол молекулярной динамики
В работе использовалась модель воды Т1Р4Р/2005, которая достаточно хорошо описывает жидкое и кристаллическое состояние, диффузные и вязкостные свойства, а также эффекты поверхностного натяжения [16-18]. Парциальные заряды атомов в липидах оценивались с использованием неограниченного метода Хартри — Фока, базиса 6-3Ю* и метода Малликена. Радиус обрезания невалентных взаимодействий составлял 1.8 нм. Вычисления проводились в NPT-ансамбле, для поддержания условий постоянства температуры и давления использовался термостат стохастической динамики и баростат Берендсена [9]. Для учета эффектов поверхностного натяжения бислоя величина давления в латеральной плоскости выбиралась отрицательной, достаточной для поддержания величины удельной площади на липид в известных из эксперимента пределах [19]. Давление баростата вдоль нормали мембраны составляло 1 бар, перпендикулярно нормали мембраны — 50 бар. В качестве модельной мембраны эукариотической клетки использовался бислой из 80 молекул ДМФХ. Для построения бислоя использовались по 5 конформаций молекул каждого типа. Молекулы липидов помещались в узлы гексагональной решетки и поворачивались на случайный угол вокруг своей оси. Площадь поверхности на одну молекулу липида составляла 0.6 нм2. Для сборки мембранных систем использовалось оригинальное программное обеспечение Б1Ьауег.
Вычислительные эксперименты проводились в программном пакете Сгошася 4.6.7 [20]. Перед проведением процедуры ИЯ производилась релаксация модельной мембраны в течение 20 нс в NPT ансамбле при температуре 300К. Начальные скорости атомов определялись с помощью генератора случайных
чисел по распределению Максвелла. Шаг интегрирования составлял 1 фс. Для процедуры зонтичной выборки стартовые конфигурации для молекулы воды выбирались с шагом 0.05 и 0.1 нм вдоль оси О,г (координаты реакции). Использовался гармонический потенциал с силовой константой 1000 кДж/(моль-нм2). Длины траекторий для набора статистики составляли 6-10 нс. Проводился анализ зависимости вычисляемых профилей энергии от начальных данных, координат, длин траекторий, а также зависимость от начальных скоростей атомов в системе.
Обработка и анализ траекторий осуществлялись с помощью встроенных модулей Сгошаев [20], для визуализации использовался пакет УМБ [21].
3. Результаты
В табл. 1 приведены основные структурные характеристики бислоев после релаксации в течение 20 нс. Усреднение значений проводилось по последним 2 нс траекторий. Для сравнения приведены данные, полученные экспериментально (Э) [22].
Таблица 1. Основные структурные характеристики моделируемых бислоев
Параметр Вислой ДМФХ
Э мд
Толщина бислоя, Е 33.8 37.1±0.2
Площадь поверхности на липид, Е2 65.4 55.2±0.2
Параметр порядка 0.184 0.218
Как видно, использованный протокол МД позволяет моделировать мембранные системы, структурные характеристики которых близки к экспериментально наблюдаемым. Полученные мембраны и разработанный протокол молекулярной динамики использовались для дальнейших вычислений.
Был проведен анализ зависимости рассчитываемых значений профиля свободной энергии от начальных данных. Рассматривалась зависимость от длин траекторий моделирования возмущенных систем, шага выборки стартовых конфигураций вдоль координаты реакции, начальных скоростей атомов в системе, а также выбора самой траектории движения малых молекул через бислой ли-пидов.
Была исследована зависимость вычисляемых результатов от положений начальных координат молекул воды в толще липидного бислоя. Для этого было выбрано 4 возможных траектории прохождения молекулы воды через липид-ный бислой, и начальные точки были выбраны на этих траекториях с шагом 0.1 нм, длины траекторий составляли 6-10 нс. Было проведено 9 расчетов с разными начальными значениями скоростей атомов в системе (рис. 3). Рассчитанные высоты потенциального барьера, приведенные в табл. 2, составили для разных расчетов от 20 до 36 кДж/моль, среднее рассчитанное значение — 27 кДж/моль, среднеквадратическое отклонение — 5.5 кДж/моль. Известные экспериментальные оценки коэффициента межфазного распределения
Кр вода-гидрофобный растворитель для воды лежат в пределах от 4.2•Ю-5 (вода-гексадекан [2]) до 1.4•Ю-3 (вода-оливковое масло [2]). Полученная нами средняя величина высоты барьера 27 кДж/моль по 9 измерениям соответству-
т-^ша-Ь/тетЬгапе п 1Л — л
ет расчетному значению Кр ' =1.6 • 10 4 и не противоречит экспе-
риментальным оценкам. В то же время довольно большая величина средне-квадратического отклонения может говорить как о стохастическом характере исследуемого процесса, так и о возможных погрешностях выбранной процедуры измерений.
Таблица 2. Рассчитанная высота барьера в экспериментах при шаге 0.1 нм и К^агт = 1000 кДж/(моль-нм2)
Номер вычислительного эксперимента 1 2 3 4 5 6 7 8 9
Величина барьера, кДж/моль 23.5 19.6 25.0 22.8 32.0 27.7 35.5 25.0 36.0
На рис. 4 приведен вид гистограмм, характеризующих статистическое покрытие конформационного пространства вдоль координаты реакции в вычислительном эксперименте. Можно отметить, что при выборе шага 0.1 нм и величине
Рис. 4. Вид гистограммы для одного из экспериментов с шагом 0.1 нм.
Заметны области низкой заселенности (отмечены на оси Ог)
силовой константы Кьагт гармонического потенциала 1000 кДж/(моль-нм2) в гистограмме присутствуют области с низкой заселенностью.
Для исследования влияния качества выборки конформационного пространства проведены дополнительные вычислительные эксперименты с шагом 0.05 нм, величина силовой константы Кьагт не менялась. Было проведено 4 измерения с различными начальными скоростями атомов. Как видно из табл. 3, уменьшение величины шага до 0.05 нм не приводит к значительному изменению величины среднеквадратического отклонения (рассчитанная величина 4.64 кДж/моль). Исследование гистограмм, характеризующих статистическое покрытие конформационного пространства вдоль координаты реакции (не показано), не указывает на наличие областей, слабо представленных в итоговой выборке.
Таблица 3. Рассчитанная высота барьера в экспериментах при шаге 0.05 нм и К^агт = 1000 кДж/(моль-нм2)
Номер вычислительного эксперимента 1 2 3 4
Величина барьера, кДж/моль 32.0 19.5 28.0 24.0
4. Выводы
Таким образом, в рамках подхода произведены оценка кинетических параметров проникновения и оценка профиля потенциала средней силы, а также изучена зависимость вычисляемых значений от начальных данных вычислительного эксперимента.
Показано, что метод Umbrella Sampling позволяет проводить исследование трансмембранного транспорта молекул воды через бислой ДМФХ и произво-
дить оценку профилей средней силы. При этом вычисленная средняя величина энергетического барьера порядка 27 кДж/моль лежит в пределах допустимых значений, известных из натурных экспериментов по измерению Кр^а1 межфазного распределения вода-гидрофобный растворитель.
Сэмплирование с шагом 0.1 нм и длинами траекторий 6-10 нс может приводить к появлению областей конформационного пространства с низкой заселенностью вдоль координаты реакции, при шаге 0.05 нм такого не наблюдалось. Однако изменение шага незначительно повлияло на измеряемое значение среднеквадратического отклонения, что может указывать в данном случае на стохастичность процесса трансмембранного транспорта.
При этом точность по 9 экспериментам с шагом 0.1 нм можно охарактеризовать как условно достаточную для приближенных оценок величины потенциального барьера, но недостаточную для оценки коэффициентов межфазового разделения.
Вычислительные расчеты выполнены с использованием ресурсов суперкомпьютерного комплекса МГУ им. М. В. Ломоносова [23] и суперкомпьютерного комплекса СВФУ «Ариан Кузьмин».
ЛИТЕРАТУРА
1. Verlet L. Computer experiments on classical fluids. I. Thermodynamical properties of Len-nard-Jones molecules // Phys. Rev. 1967. V. 159, N 1. P. 98-103.
2. Walter A., Gutknecht J. Permeability of small nonelectrolytes through lipid bilayer membranes // J. Membr. Biol. 1986. V. 90, N 3. P. 207-217.
3. Diamond J. M., Katz Y. Interpretation of nonelectrolyte partition coefficients between dimyri-stoyl lecithin and water // J. Membr. Biol. 1974. V. 17, N 2. P. 121-154.
4. Shaitan K. V., Antonov M. Y., Tourleigh Y. V., Levtsova O. V., Tereshkina K. B., Nikolaev I. N., Kirpichnikov M. P. Comparative study of molecular dynamics, diffusion, and permeability for ligands in biomembranes of different lipid composition // Biochem. Suppl. Ser. A Membr. Cell Biol. 2008. N 2. P. 73-81.
5. Антонов М. Ю., Науменкова Т. В., Левцова О. В., Шайтан К. В. Исследование динамики и трансмембранной диффузии методами компьютерного моделирования // Суперкомпьютерные технологии математического моделирования: тр. II Междунар. конф. 2014. С. 16-29.
6. Jorgensen W. L., Chandrasekhar J., Madura J. D., Impey R. W., Klein M. L. Comparison of simple potential functions for simulating liquid water //J. Chem. Phys. 1983. V. 79, N 2. P. 926.
7. Jorgensen W. L., Maxwell D. S., Tirado-Rives J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids //J. Amer. Chem. Soc. 1996. V. 118, N 45. P. 11225-11236.
8. Holm C., Kremer K. Advanced computer simulation. V. 173. Berlin; Heidelberg: Springer-Verl., 2005.
9. Berendsen H. J. C., Postma J. P. M., van Gunsteren W. F., DiNola A., Haak J. R. Molecular dynamics with coupling to an external bath // J. Chem. Phys. 1984. V. 81, N 8. P. 3684.
10. Olbrich K., Rawicz W., Needham D., Evans E. Water permeability and mechanical strength of polyunsaturated lipid bilayers // Biophys. J. 2000. V. 79, N 1. P. 321-327.
11. Bean R. C., Shepherd W. C., Chan H. Permeability of lipid bilayer membranes to organic solutes // J. Gen. Physiol. 1968. V. 52, N 3. P. 495-508.
12. Buch I., Sadiq S. K., de Fabritiis G. Optimized potential of mean force calculations for standard binding free energies // J. Chem. Theory Comput. 2011. V. 7, N 6. P. 1765-1772.
13. Kastner J. Umbrella sampling // Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011. V. 1, N 6. P. 932-942.
14. Torrie G. M., Valleau J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling // J. Comput. Phys. 1977. V. 23, N 2. P. 187-199.
