Динамический молекулярный дизайн био- и наноструктур Текст научной статьи по специальности «Химические науки»

УДК 573.3

Динамический молекулярный дизайн био- и наноструктур

К. В. Шайтан, Е. В. Турлей, Д. Н. Голик, К. Б. Терешкииа, О. В. Левцова, И. В. Федик, А. К. Шайтан, А. Ли, М. П. Кирпичников

КОНСТАНТИН ВОЛ ЬДЕМАРОВИЧ ШАЙТАН — заместитель заведующего кафедрой биоинженерии, доктор физико-математических наук, профессор Биологического факультета МГУ им. М.В.Ломоносова. Область научных интересов: молекулярная биоинженерия, молекулярное моделирование, био- и наноструктуры, теоретическая биофизика. E-mail [email protected]

ЕГОР ВЛАДИМИРОВИЧ ТУРЛЕЙ — аспирант кафедры биоинженерии Биологического факультета МГУ им. М.В.Ломоносова. Область научных интересов: молекулярное моделирование, био- и наноструктуры. E-mail [email protected]

ДМИТРИЙ НИКИТИЧ ГОЛИК — аспирант кафедры биоинженерии Биологического факультета МГУ им. М.В.Ломоносова. Область научных интересов: молекулярное моделирование, био- и наноструктуры. E-mail [email protected]

КСЕНИЯ БОРИСОВНА ТЕРЕШКИНА — ассистент кафедры биоинженерии Биологического факультета МГУ им. М.В.Ломоносова. Область научных интересов: молекулярное моделирование, био- и наноструктуры. E-mail [email protected]

ОЛЬГА ВЛАДИМИРОВНА ЛЕВЦОВА — аспирантка кафедры биоинженерии Биологического факультета МГУ им. М.В.Ломоносова. Область научных интересов: молекулярное моделирование, био- и наноструктуры. E-mail [email protected]

ИГОРЬ ВИКТОРОВИЧ ФЕДИК — аспирант кафедры биоинженерии Биологического факультета МГУ им. М.В.Ломоносова. Область научных интересов: молекулярное моделирование, био- и наноструктуры. E-mail [email protected]

АЛЕКСЕЙ КОНСТАНТИНОВИЧ ШАЙТАН — студент 5-го курса кафедры физики полимеров и кристаллов Физического факультета МГУ им. М.В.Ломоносова. Область научных интересов: молекулярное моделирование, био- и наноструктуры, физика полимеров. E-mail [email protected]

АНЬБАН ЛИ — аспирант кафедры биоинженерии Биологического факультета МГУ им. М.В.Ломоносова. Область научных интересов: молекулярное моделирование, био- и наноструктуры. E-mail [email protected],

МИХАИЛ ПЕТРОВИЧ КИРПИЧНИКОВ — академик РАН, заведующий кафедрой биоинженерии Биологического факультета МГУ им. М.В.Ломоносова, доктор биологических наук. Область научных интересов: белковая инженерия, биоинженерия. E-mail [email protected].

199992 Москва, Ленинские горы, МГУ им. М.В.Ломоносова, Биологический факультет, тел. (495)939-23-74.

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

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

Существуют варианты как равновесной молекулярной динамики, когда исследуется динамика молекул, обусловленная только межатомными взаимодействиями, так и неравновесной (или управляемой) молекулярной динамики (Steered Molecular Dynamics [4]), позволяющей изучать влияние дополнительных внешних силовых воздействий и (или) специальных граничных условий. Применение метода неравновесной молекулярной динамики по отношению к сложным и гетерогенным молекулярным системам более перспективно. Это связано с тем, что появляется принципиальная возможность целенаправленно планировать вычислительный эксперимент и изучать отклик системы на внешнее воздействие.

В данной статье рассмотрены некоторые направления использования управляемой молекулярной динамики для исследования сложных молекулярных систем, развиваемые на кафедре биоинженерии Биологического факультета МГУ.

Динамика биомембран

Биологические мембраны до настоящего времени являются предметом пристального изучения [5—9]. Механизм их функционирования определяется химическим составом, пространственной организацией и молекулярно-кинетическими свойствами. Динамика биомембран изучена относительно мало, данные по кинетическим свойствам фосфолипидного бислоя с учетом его анизотропии и неоднородности практически не доступны.

Ниже обсуждается проблема неоднородности и анизотропности фосфолипидного бислоя по отношению к диффузионным процессам. Эта проблема изучалась в нашей работе на модели мембраны — фосфолипидного бислоя на основе 1-пальмитоил-2-олеоил-от-глицеро-З-фосфатидилхолина (ПОФХ) (64 молекулы), степень сольватации — 44 молекулы воды на молекулу липида. Для молекул воды была выбрана модель TIP3P, в которой длины связи и валентные углы не фиксированы. Площадь в расчете на молекулу липида варьировалась в районе экспериментальных значений (62-68 Á2, [10-13]).

Расчеты молекулярной динамики такой системы проводились с помощью пакета программ PUMA [14], который был специально модифицирован с целью включения управляющих воздействий [6, 7]. Решение системы классических уравнений движения атомов производилось в силовом поле AMBER 99 [15]. Использовались периодические граничные условия, столкновительный термостат [14] (Г= 300 К) и баростат Берендсена [16]. Достижение локального равновесия в системе контролировалось по флуктуациям объема, давления и температуры [6].

Диффузия модельных сфер в биомембране

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

Á

сравнимы с радиусами атома углерода и небольшой функциональной группы, соответственно), взаимодействующие с атомами системы только посредством сил Ван-дер-Ваальса (константа взаимодействия е = 0,15 ккал/моль). К модельным частицам прикладывалась дополнительная постоянная сила Fext, действующая в направлении нормали и параллельно плоскости мембраны. В первом случае пробный шар предвари-

Á

сти мембраны, во втором — помещался в центр мембраны, и система выдерживалась в течение 2 пс (стадия релаксации). Величина силы FCKÍ составляла от 0,3 до 4 ккалДмоль • Á) [1 ккалДмоль • Á) = 7 "Ю-6 дин].

Á

Á

20 а

10

< 0 sí -10

1 2 К 4 5 , t, пс

-20

30

о с

^ 20 10

ш [ U J

i V 2 3 t,

Рис. 1. Кинетические зависимости прохождения модельной сферы через бислой ПОФХ под действием внешней силы в направлении нормали к бислою.

Положение сферы на оси г (нормали) (а) и скорость продвижения сферы (о).

Центр бислоя находится в точке г=0, границы бислоя — в положениях £=±20 А. Ван-дер-ваальсов радиус шара 2 А, действующая сила 10 ккалДмоль • А). На рис. 1о скорость усреднена по интервалу 0,1 пс. Вверху показана структура бислоя.

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

5 -

4 -

С о

л

ь ,

о 3

О

и

го

К

03

2 -

1 -

2 4 6 8 10

^ехЬ ккал/(моль«А).

С о

л ь о о и го к

со

2,0

1,5

1,0

0,5

Рис. 2. Зависимость вязкости в модельной мембране (ПОФХ) от внешней силы, действующей по нормали (а) и латерально (б) к плоскости мембраны.

Ван-дер-ваальсов радиус пробного шара 2 А. На рис. 2а: 1 — полярные «головки» липида; 2 — алкильные цепи; 3 — водная фаза. Сплошная линия — данные после предварительной релаксации мембраны в течение 500 пс, пунктирная — после окончательной релаксации в течение 1 не

0

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

Метод молекулярной динамики привел к следующим результатам. Проход через мембрану сфер радиусом 2 А при силе 1 — 10 ккалДмоль • А) осуществляется за время менее 2 не. В остальных случаях сферы застревают на поверхности либо успевают погрузиться в мембранный слой только на некоторую глубину. При

АА

А

сравнимо с воздействием приложенной силы.

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

Исследовался также мембранный перенос частиц при условии латерально приложенной силы, равной 1, 2, 4 и

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

Коэффициент вязкого трения у определяется как отношение внешней силы к скорости V дрейфа частицы:

1= ^ехА

Формально коэффициент трения можно пересчитать в терминах коэффициента диффузии, пользуясь известным соотношением Эйнштейна, и в терминах микровязкости среды с использованием формулы Стокса.

Надо сказать, что имеющиеся данные относительно вязкости при движении частицы по нормали к поверхности мембраны и в латеральном направлении в центре бислоя весьма ограничены. Экспериментально установлена усредненная вязкость поверхностного слоя, составляющая от 30 до 190 сПз для различных

липидных мембран [17—19]. Для ПОФХ также известна экспериментальная оценка усредненной вязкости, она составляет 18 сПз [20]. Точную оценку можно получить определением микровязкости для отдельных структурных и динамически неоднородных областей мембраны. В первом приближении можно выделить области липидных «головок» и алкильных цепей. На рис. 2а приведена зависимость микровязкости различных участков мембраны на основе ПОФХ от величины внешней силы, действующей на частицу радиусом А

мембраны в поперечном направлении для частиц такого размера не превышает 6 сПз, вязкость в центральной части бислоя в несколько раз меньше.

На рис. 2а приведены также значения эффективной вязкости воды, рассчитанные с использованием А

0,4 сПз, что в 2—3 раза меньше экспериментального значения. Рассчитанные значения согласуются с известными оценками вязкости воды в модели Т1РЗР [21].

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

А

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

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

А

40 20 0 -20 -40

20 \ 40 60 80 120 пс

Рис. 3. Динамика прохождения остатка аланина через липидную мембрану.

Рассчитывалось положение геометрического центра остатка. Границы бислоя находятся в положениях ±22 А. В точке излома направление внешней силы было изменено на противоположное

С о

л ь о о и го к СО

Рис. 4. Зависимость микровязкости от эффективного радиуса пробных аминокислотных остатков в системе ПОФХ-вода:

1 — полярные «головки» липида; 2 — алкильные цепи;

А

(эффективный радиус 3,1 А, рис. 3) и модельной сферы (ван-дер-ваальсов радиус 2 А). Рассчитанные значения вязкости приведены на рис. 4. В случае многоатомных молекул (триптофан, аланин) сила прикладывалась равномерно ко всем атомам системы.

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

тельна к размеру частиц.

Латеральная самодиффузия липидов в бислое

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

Коэффициент латеральной самодиффузии липидов определялся в соответствии с известным соотношением

(о - х(0)}2 + О _ ,(0)})= 4/у

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

Усреднение проводится по всем молекулам липида. При вычислении коэффициента латеральной самодиффузии ПОФХ траектория общей длиной 3 не была разделена на 12 участков по 250 пс. Временные зависимости квадратов смещения центров масс молекул ПОФХ (рис. 5) усредняли по всем 12 участкам траектории.

Вычисленное для бислоя ПОФХ значение Охг = = 2,4 • Ю-7 см2/с, что близко к данным по квазиупругому рассеянию нейтронов на бислоях дипальмитоил-фосфатидилхолина (1 • Ю-7 см2/с [22]) и диолеоил-

фосфатидилхолина (2 • Ю-7 см2/с [23]). Результаты измерения методом импульсного ЯМР для бислоев ПОФХ дают значения 2,0 • Ю-7 см2/с при 298 К и 2,5 • Ю-7 см2/с при 303 К [24], что практически совпадает с полученными нами значениями. Отметим, что рассчитанное значение коэффициента латеральной диффузии оказалось гораздо ближе к экспериментальным значениям, чем коэффициенты, получаемые в молекулярно-динамических расчетах в соответствии с другими методиками, использующими, в том числе, другие силовые поля и схемы расчета парных взаимодействий (для дипальмитоилфосфатидилхолиновых бислоев 3±0,6 • Ю-7 см2/с [25]; 4,0-4,5 • Ю^7 см2/с [26] и бислоя ПОФХ 7,3 • 10~5 см2/с [9]).

Оценка гидрофобных свойств аминокислотных остатков на границе вода-мембрана

Оценка параметров гидрофобности используется для построения профилей гидрофобности белков [27, 28], а также широко применяется в фармакокине-тике [29]. Гидрофобность является одним из важней-

Л 2,5

^ 2,0

¡Ъ +

1,5

В 1,0

и

0,5

Рис. 5. Средний квадрат смешения центров масс липидов в плоскости бислоя ПОФХ.

Прямая — линейная аппроксимация

ь о о ас ь о ч с к св ас

J3

ч а> ь

Я

о о

¡E Ь

О

1,0 -■

0,8 -

0,6 -

0,4 -

0,2 -

0 -

-10 0

Ось г, А

Рис. 6. Профили плотности воды (I) и гексана (2) в системе вода/гексан вдоль оси, перпендикулярной разделу фаз.

Профили нормированы к единице

Я

ь о о ас ь к о

О.

<а со

J3

ь о о ас ь о ч С

1,0

0,8

0,6

0,4

0,2

0 5 10 15

Ось г, А

Рис. 7. Вероятность распределения положений центров масс остова (I) и боковой группы (2) аминокислотного остатка в системе вода/гексан

0

ших параметров для реакций, катализируемых ферментами. Сайты связывания ферментов весьма восприимчивы к гидрофобным участкам молекул субстрата [30]. Методами управляемой молекулярной динамики нами исследовалось поведение отдельных аминокислотных остатков на границе модельных систем «вода/неполярный растворитель».

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

В нашей модели аминокислотные остатки принимались в форме бифункционального мономерного звена —NH—CH(R)—СО—, т.е. как незаряженный мономер с формально ненасыщенными связями. Расчет проводили с помощью программного комплекса PUMA. Выбирались периодические граничные условия, силовое поле AMBER 99 и модель воды TIP3P с нефиксированными внутренними степенями свободы. Для системы вода/гексан проводилось баростатирование в направлении, перпендикулярном границе раздела фаз. Система вода/вакуум моделировалась в УУКТ-ансамбле, плотность воды соответствовала нормальным условиям.

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

Рассмотрим поведение фенилаланина на границе раздела фаз вода/гексан. На рис. 6 приведены усредненные по времени профили плотности воды и гексана. Видно, что граница раздела фаз представляет собой переходный слой толщиной около 5 Á. На рис. 7 приведены распределения положений центров масс остова и боковой цепи фенилаланина. Максимумы распределений практически точно соответствуют сере-

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

Рассмотрим далее динамику фенилаланина в системе с виртуальной отталкивающей стенкой (система вода/вакуум). Из рис. 8 видно, что для этой системы крутизна профиля плотности водной фазы больше, чем в системе вода/гексан. Отсутствие неполярной фазы гексана несущественно меняет форму и ширину распределений аминокислотного остатка (рис. 9). Хорошо заметно смещение наиболее вероятного положения центров масс остова и бокового радикала фенилаланина в сторону водной фазы по сравнению с системой вода/гексан. Наиболее вероятное положение молекулы при этом находится в приграничном к разделу фаз слое воды, а не на самой границе раздела фаз. Отметим, что динамика ориентации аминокислотного остатка в системах вода/вакуум и вода/гексан отличается незначительно.

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

1,0

0,8 0,6 0,4 0,2

1,0

я ь о о ас 0,8

о

О.

0> СО

0,6

ь о о ас ь о ч С

0,4

0,2

-40 -20 0 20 40

z, Á

2 А

/ \ 7

> / / \ *

i \ ;

1 \\

/ / \ \

i i \ \

/ . \ ч

15

20 Ось z, Á

25

Рис. 8. Профиль плотности воды в системе вода/вакуум вдоль оси, перпендикулярной отталкивающей стенке

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

характе-

Динамика функционирования ионных каналов

Ионные каналы, формируемые полипептидными структурами, являются новыми и сложными объектами для молекулярного моделирования, и здесь методы управляемой (неравновесной) динамики оказываются наиболее эффективными.

Ниже рассмотрены примеры определения ионной подвижности в катион(Ыа+)-проводящем грамициди-новом канале и канале ацетилхолинового рецептора, а также в анион(СР)-проводящем канале рецептора глицина.

Динамика грамицидинового канала

Грамицидин А — природный антибиотик, полипептид, активная форма которого является димером. При встраивании в биологическую мембрану он образует канал, проводящий протоны и одновалентные катионы по градиенту их концентрации, что приводит к понижению мембранного потенциала покоя и сопротивления мембраны [31—33].

Известны две конформации грамицидинового канала в мембране: спиральный димер и двойная спираль [34], различающиеся по стабильности и функциональной активности [35]. В нашей работе рассматривалась полноатомная структура комплекса фосфолипидного бисяоя ПОФХ с грамицидино-вым каналом, состоящим из двух молекул грамицидина А в конформации спирального ди-мера. Расчеты подвижности ионов в канале проводились с помощью пакета программ молекулярной динамики Gro-macs 3.1.4 (потенциальное поле Gromos-96) [36]. Использовалась стохастическая модель динамики с временным параметром затухания стохастических кор-

Рис. 9. Вероятность распределения положений центров масс остова (I) и боковой группы (2) аминокислотного остатка в системе вода/вакуум

реляций 0,02 пс, периодические граничные условия, баростат Берендсена (вдоль оси нормали к плоскости мембраны поддерживалось давление 1 бар, перпендикулярно нормали — 260 бар, коэффициент баростатирова-ния 10 пс), модель воды Т1РЗР, Г = 300 К, диэлектрическая проницаемость бралась равной 2. Катион натрия помещался в водную среду на расстоянии 5 А над порой канала. Вокруг катиона достаточно быстро формировалась гидратная оболочка из шести молекул воды. К иону натрия прикладывалось ускорение 25 нм/пс2 вдоль оси по нормали к мембране (что эквивалентно внешней А

При вхождении в полость поры катион натрия раздвигает атомы полярных «головок» липидов и молекул грамицидина. На этой стадии происходит замена шести молекул гидратной оболочки иона натрия в водной среде на две молекулы воды в поре ионного канала. Уменьшение скорости движения иона соответствует локальным энергетическим минимумам (рис. 10), связанным с притяжением катиона к отрицательно заряженным

7,0

2,0

10

15

/, пс

20

25

30

Рис. 10. Динамика прохождения иона через пору грамицидинового канала.

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

0

0

0

5

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

Подвижность, или эффективный коэффициент диффузии единичного катиона в ионном канале можно оценить, исходя из модели вязкого трения:

D(z)= — , J=-

Y v

где D(z) — коэффициент диффузии по оси z, Y — коэффициент трения; F — приложенная внешняя сила; v — средняя скорость прохождения иона через канал.

Оценка по этой формуле дает для коэффициента диффузии в грамицидиновом канале значение 0,6 • Ю-5 см2/с, что более чем в два раза меньше значения коэффициента диффузии иона в воде TIP3P, полученного тем же методом (1,5 • Ю-5 см2/с). Отметим, что расчеты, проведенные для катиона натрия в канале ацетилхолинового рецептора и хлорид-аниона в канале глицинового рецептора, дают близкие значения коэффициентов диффузии. Таким образом, применение методов управляемой динамики позволяет определять как сайты связывания ионов в каналах, так и вычислять характеризующие проводимость величины за доступные времена.

Динамика ионов в каналах глицинового и ацетилхолинового рецепторов

Глициновые и ацетилхолиновые рецепторы относятся к семейству лиганд-зависимых ионных каналов, обеспечивающих быструю передачу сигнала через нейроны в различные отделы центральной нервной системы [37]. Глициновые рецепторы наряду с рецепторами Y~aMMHOMacnHHott кислоты ГАМКд и ГАМ К-формируют тормозящие синапсы, а ацетилхолиновые рецепторы — возбуждающие синапсы.

Все вышеперечисленные лиганд-зависимые рецепторы имеют пентамерную белковую структуру, как и многие другие лиганд-зависимые рецепторы [38]. Каждая субъединица рецепторного белка (рис. 11а) состоит из четырех a-спиралей и содержит четыре трансмембранных домена ТМ1—ТМ4 [39]. Активация глициновых рецепторов, когда происходит открытие хлорных каналов, обеспечивается связыванием рецептора с глицином. Для ацетилхолиновых рецепторов лигандом является ацетилхолин. Считается, что основную роль в миграции ионов играют домены ТМ2, непосредственно формирующие пору канала (рис. 11 б) и таким образом являющиеся наиболее важным компонентом рецептора с точки зрения функциональной активности [40, 41]. С мутациями в субъединицах глицинового рецептора связаны некоторые патологии, в частности, рефлекс гиперстраха в ответ на неожиданные раздражения [42]. Точечные мутации в субъединицах ацетилхолинового рецептора способны приводить к развитию ночной эпилепсии лобной доли мозга и синдрома Дауна.

Для изучения диффузии ионов и аквакомплексов через ионные каналы глицинового и никотинового ацетилхолинового рецепторов (nAChR), а также для изучения механизма открытия канала nAChR была построена модель ионного канала, состоящая из ТМ2-доменов а,в,5,а,7-субъединиц. Модель канала nAChR строилась на основе данных рентгеноструктурного анализа, взятых из банка белковых структур (PDB: ÍOED) [43]. В этой структуре внутрь канала обращены два незаряженных кольца, образованных боковыми

aa Согласно результатам исследований мутагенных форм рецептора, именно эти остатки выполняют в канале роль «ворот». Канал в модели IOED находится в закрытом состоянии.

Модель канала глицинового рецептора была реконструирована по гомологии с каналом ацетилхолиново-

А

NHn

Рис. 11. а-Субъедииица глицинового рецептора. Мембранная топология (А); более темный на рисунке — трансмембранный домен ТМ2, формирующий пору (по [39]). Схематическое представление двух трансмембранных доменов ТМ2 и соответствующих сечений пор гомомерных рецепторов (слева направо): глицинового рецептора дикого типа; глицинового рецептора с тройной мутацией из а1-субъедиииц; никотинового ацетилхолинового рецептора из а7-субъединиц (по [41]) (Б). Показаны три заряженных кольца, оказывающие влияние на миграцию ионов: внеклеточное (внеш.), внутреннее (внутр.) и цитоплазматическое (цит.)

Рис. 12. Пентамерная структура глицинового канала, стабилизированная углеводородным кольцом

го рецептора. При этом использовались данные рент-геноструктурного анализа по расшифровке структуры а-спирзди ТМ2, подтипа al (PDB:lMOT). Для стабилизации пентамерной структуры канала использовалось неразветвленное углеводородное кольцо C^Hjoq, выполнявшее роль скрепляющего обруча (рис. 12). Для предотвращения «схлопывания» канала использовалась специальная процедура усиления энергии невалентных взаимодействий между некоторыми атомами

a

кими к ним атомами углерода кольца. Все расчеты проводились методом управляемой молекулярной динамики в виртуальной среде с вязкостью 1 сПз при температуре 300 К с использованием потенциального поля AMBER 99.

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

женных, так и незаряженных атомных групп, которые имитировали атомы и ионы натрия и хлора, а также их гидратированные комплексы. Последние моделировались сферами, масса и радиус которых соответствовали гидратированным ионам. Перенос атомных групп осуществлялся под действием постоянной силы 1 ккалДмоль • А), направленной вдоль оси канала. Оказалось, что незаряженные атомы С1 и Na свободно проходят через канал в конформации lOED. Однако для заряженных частиц этот канал закрыт. Ион Na+ притягивается к остаткам B(P)Gly272, C(8)Gly280. Комплекс Na(H20)6+ остается в канале между цепями A(a) и В(Р), притягиваясь к остаткам A(a)Leu263, A(a)Thr267, В(Р)^р268, B(P)Pro271, B(P)Gly272. СГ и C1(H20)~6 покидают область канала между спиралями. Наиболее узкое место в канале играет роль ворот, которые образованы кольцом из следующих остатков: A(a)Val255, B(P)Val261, C(S)Val269, D(a)Val255 и E(y)^u264. Диаметр «ворот» около 6 А. На рис. 13 показаны динамика прохождения незаряженного гид-ратного комплекса под действием постоянной силы А

«ворот». Здесь обращает на себя внимание то обстоятельство, что стерические затруднения не являются единственным фактором, определяющим проницаемость канала. Очень важная роль отводится электростатическим эффектам.

На рис. 14 представлена динамика прохождения иона СГ через пору канала глицинового рецептора под

А

ной вдоль оси канала. Как оказалось, пора глицинового рецептора имеет три основных минимума энергии для отрицательных ионов. Два минимума соответствуют положительным кольцам внутри канала, образованным аргининовыми остатками (см. рис. 116). Выявляется также минимум энергии в поре канала, соответствующий положительно заряженными атомами бокового радикала остатка метионина Met. Значение эффективного коэффициента диффузии иона СГ составляет Ю-5— 10~6 см2/с в различных частях канала.

t, пс

Рис. 13. Динамика прохождения незаряженной гидратированной атомной группы через пору канала ацетилхолинового рецептора.

Канал находится в закрытой конформации. Средней стрелкой показана область «ворот»

0 50 100 150 200 250

/, пс

Рис. 14. Динамика прохождения иона С1 через пору канала глицинового рецептора

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

Взаимодействие углеродной нанотрубки с фосфолипидным бислоем

Наноструктуры и их комплексы с биологическими макромолекулярными системами — новые объекты для применения методов неравновесной (управляемой) молекулярной динамики. В работе [44] рассматривается, в частности, процесс спонтанного встраивания модельной нанотрубки в фосфолипидный бислой

Рис. 15. Прохождение углеродной трубки через фосфолипидный бислой.

Молекулы фосфолипида, выталкиваемые из слоя, изображены более подробно

с использованием крупномодульного (соагее-§гат) приближения. Однако до сих пор практически нет публикаций по динамике взаимодействия биомембран с углеродными нанотрубками в полноатомном приближении.

В нашей работе изучалась динамика проникновения углеродной нанотрубки (диаметр 13,5 А и длина 32 А) через модельную мембрану — гидратированный бислой ПОФХ (см. выше) под действием внешней силы. Нанотрубка ориентировалась по нормали к плоскости мембраны. К атомам углеродной нанотрубки равномерно прилагалась сила в направлении нормали. Результирующее давление составляло порядка 7 • 104 бар (для сравнения отметим, что давление детонации тротила существенно выше, порядка 2 • 105 бар). При таком давлении нанотрубка проходит через фосфолипидный бислой за время порядка 200 пс (рис. 15), выталкивая из бислоя две молекулы фосфолипида. Проведенные расчеты позволяют оценить соответствие временных масштабов силовым воздействиям, которое важно учитывать, например, при постановке аналогичных экспериментов атомно-силовой микроскопии. В данном случае внешняя сила соответствует напору кантилевера с присоединенной на конце нанотрубкой.

Взаимодействие полипептида с углеродной нанотрубкой

Изучение взаимодействия нанотрубок с биологическими молекулами имеет большое практическое значение в связи с проблемой селективной доставки лекарств в клетки. Модельные исследования в этом важнейшем направлении уже начинают развиваться. Так, в работе [45] приводится расчет процесса проникновения однонитевого ДНК-олигонуклеотида в углеродную нанотрубку в водной среде. В работе [46] изучалось прохождение РНК под действием приложенной силы через отверстия в коротких нанотрубках, составляющих монослой. В [47] изучалась адсорбция амилозы на нанотрубке в водной среде и ее проникновение внутрь нанотрубки.

О пс 100 пс 200 пс

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

Опс 50 пс 100 пс 150 пс

Рис. 17. Последовательные стадии самосборки комплекса пептида и нанотрубки при 1000 К

В наших численных экспериментах было обнаружено явление самосборки полиаланина* и углеродной нанотрубки с образованием структуры в виде спирали полиаланина внутри нанотрубки. Показано, что при 300 К за время порядка 200 пс происходит адсорбция (рис. 16) полипептида в а-спиральной конформации на поверхности нанотрубки (исходное положение пептида на расстоянии 30 А от нанотрубки).

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

На рис. 17 показана детальная картина акта самосборки структуры полипептид-нанотрубка. Находясь на внешней стенке нанотрубки, молекула пептида достаточно продолжительное время остается около ее центральной части. Изредка один из концов пептида оказывается около отверстия. Постепенно вследствие флуктуаций молекула пептида начинает перемещаться вдоль нанотрубки. Затем конец пептида приближается

* NH2CH(CH3)CO-[ NHCH(CH3)C01„- NHCH(CH3)COOH.

к отверстию, после чего вся макромолекула быстро проникает в нанотрубку. Эта стадия при 1000 К длится 130 пс. При 2000 К самосборка идет по такому же механизму, однако при повышении температуры процесс становится более обратимым, и за счет этого время от начала до завершения акта встраивания увеличивается до 300 пс.

Продолжительность формирования активной для самосборки конфигурации при температуре 1000 К составила 4,64 не, при 2000 К — 0,655 не. Это дает оценку величины энергии активации порядка 7,8 ккал/моль. Ожидаемое время самосборки при 300 К составляет 43 мкс. Отметим, что рассматриваемый процесс моделировался в вакууме. Для процесса, протекающего в растворителях, энергия активации самосборки должна быть, по-видимому, ниже вследствие влияния энергии сольватации полипептида и нанотрубки.

Динамика функциональных наноструктур. Наношнриц

Рассмотренный выше комплекс полипептида и нанотрубки с закрытым концом может быть в принципе использован для доставки пептида (или иной молекулы) через биологическую мембрану в клетку или отдельный компартмент. Зачастую требуется селективная доставка низкомолекулярных синтетических веществ, имитирующих действие природных биомакромолекул и имеющих терапевтический потенциал. Эти конструкции могут быть также использованы и для изучения механизмов молекулярного распознавания [48].

Рис. 18. Последовательные стадии выталкивания пептида в мембрану (А) и в воду (Б)

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

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

А

рения 26 и 13 пс соответственно) до значений порядка радиуса нанотрубки. При таких условиях создавался практически «нановзрыв», и система срабатывала как «нанопушка».

На рис. 18 приведен сценарий выброса пептида при таких экстремальных параметрах «выстрела» с максимумом давления в нанотрубке порядка 105 бар. В момент «выстрела» нанотрубка несколько деформируется, но эти деформации не выходят за пределы ее прочности. По окончании процесса выталкивания пептида нанотрубка полностью восстанавливает первоначальную конформацию за время порядка 3 пс.

В ходе рассматриваемого процесса молекула поли-аланина испытывает конформационные изменения.

Начальная спиральная конформация наиболее сильно деформируется при выбросе полипептида в вакуум и менее всего — в мембрану. По-видимому, среда играет в этом процессе демпфирующую и структурирующую роль. При уменьшении скорости «нановзрыва» конформационные изменения молекулы полипептида обычно уменьшаются.

Фолдинг модельных полимерных структур

Закономерности формирования пространственной структуры полимеров и биополимеров остаются малоизученными, несмотря на их важнейшую роль при функционировании биологических систем [49, 50]. Метод молекулярной динамики позволяет исследовать сворачивание (фолдинг) относительно простой модели полимерной цепи. Отметим, что специфическое сворачивание характерно даже для обычной свободно сочлененной молекулярной цепи с ван-дер-вальсовыми взаимодействиями между звеньями [51—53].

Укажем работы в этом направлении. Проведены вычислительные эксперименты по сворачиванию го-мополимерной цепи [54]. Взаимодействия соседних атомов определялись энергией химической связи. Взаимодействия атомов (звеньев цепи), разделенных двумя и более связями, описывались потенциалом Лен-нард-Джонса. Наблюдалось возникновение разнообразных структур: всевозможных слоев, правых и левых спи-

/=Опс

( = 25 :

Рис. 19. Рефолдинг полимерной цепи при взаимодействии с модельной

нанотрубкой

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

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

от начальной конфигурации полимера. Следует отметить, что в биологических системах молекулярные комплексы, способствующие специфическому сворачиванию, часто имеют форму трубки (шапероны группы ШрбО, мембранные каналы) [55—57].

Рассмотрен рефолдинг полимерной глобулы в спираль при специфическом взаимодействии полимера с модельной нанотрубкой. Полимерная цепь, проходя через нанотрубку, разворачивается, а затем формирует новую упорядоченную структуру (рис. 19). Модельное устройство функционирует при определенном строении поверхности потенциальной энергии для взаимодействия нанотрубки с полимерной цепью. Устройство этой поверхности иллюстрируется сечениями, изображенными на рис. 20.

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

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

35 х

30

X..

25 х

20 х.

15 х

10 х

5 0

-5 х

-2 "

0 2.Х -5 у 4 -10 ^

-800 -600 -400

200 0

-200 -400 -600 -800

800

400

и и

к я

и

о.

<а ас Г)

-400

-800 -10"

10

15 20

С, 30 25

35

10 -5

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

0

г

0

х

0

Заключение

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

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

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

* * *

Авторы признательны Российскому фонду фундаментальных исследований (грант 04-04-49645), Мин-обрнауки РФ, Роснауке и Правительству Москвы за финансовую поддержку проводимых исследований.

ЛИТЕРАТУРА

1. Gao H. In: European white book on fundamental research in materials science. Eds. M.H. Van der Woorde e. a. Max Planck Gesellschaft, 2001, p. 144-148.

2. Frenkel D., Smit B. Understanding molecular simulation: from algorithms to applications, 2nd ed. San Diego: Academic Press, 2002.

3. Шайтан К.В., Терешкина К.Б. Молекулярная динамика белков и пептидов (методическое пособие). М.: Ойкос,

2004, 103 с.

4. Park S., Schulten К. J. Chem. Phys., 2004, v. 120, p. 5946-5961.

5. Chizmadzhev Y.A. Bioelectrochem., 2004, v. 63, p. 129—136.

,

2005, т. 50, № 6, c. 491-502.

7. Турлей E.B., Шайтан К.В., Балабаев H.K. Биофизика, 2005, т. 22, № 6, с. 1042-1047.

8. Rabinovich A.L., Balabaev N.K., Alinchenko M. G., Voloshin V.P., Medvedev N.N., Jedlovszky P. J. Chem. Phys., 2005, v. 122, p. 84906.

9. Heller H., Schaefer M., Schulten К. J. Phys. Chem., 1993, v. 97, p. 8343-8360.

10. Pabst G., Rappolt M., Amenitsch H., Laggner P. Phys. Rev. E.,

2000, v. 62, p. 4000-4009.

11. Isralewitz В., Gao M., Schulten К. Curr. Opin. Struct. Biol.,

2001, v. 11, p. 224-230.

12. Smaby J.M., Mornsen M.M., Brockman H.L., Brown R.E. Biophys. J., 1997, v. 73, p. 1492-1505.

13. Hyslop P.A., Morel В., Sauerheber R.D. Biochem., 1990, v. 29, p. 1025-1038.

14. LemakA.S., Balabaev N.K Mol. Simul., 1995, v. 15, p. 223-231.

15. Wang J., Cieplak P., Kollman P.A. J. Comp. Chem., 2000, v. 21, p. 1049-1074.

16. Berendsen H.J.C., Postma J.P.M., van Gunsteren W.F., DiNola A., Haak J.R. J. Chem. Phys., 1984, v. 81, p. 3684-3690.

17. Nagle J.F, Zhang R., Tristram-Nagle S., Sun W.J., Petrache H.I., SuterRM. Biophys. J., 1996, v. 70, p. 1419-1431.

18. Wiener M.C., White S.H. Ibid., 1992, v. 61, p. 434-447.

19. Reiter G., Siam M., Falkenhagen D., Gollnentsch W., Baurecht D., Fringeli U.P. Langmuir., 2002, v. 18, p. 5761-5771.

20. SeeligA., See tig J. Biochem., 1977, v. 16, p. 45-50.

21. Seelig J., SeeligA. Q. Rev. Biophys., 1980, v. 13, p. 19-61.

22. Sackmann E. In: Handbook of Biological Physics. V. 1. Structure and Dynamics of Membranes. Eids. R. Lipowsky, E. Sackmann. Amsterdam: Elsevier Sei. Pubs. B.V., 1995, p. 213-304.

23. Pfeiffer W., Henkel Т., Sackmann E., Knoll W., Richter D. Europhys. Lett., 1989, v. 8, p. 201-206.

24. Filippov A., Oradd G., Lindblom G. Biophys. J., 2003, v. 84, p. 3079-3086.

25. Essmann U., Berkowitz. M.L. Ibid., 1999, v. 76, p. 2081-2089.

26. Ryg Т., Murzyn K., Pasenkiewicz-Gierula M. Acta Biochim. Polon., 2003, v. 53, p. 789-798.

27. KyteJ., Doolittle R.F. J. Mol. Biol., 1982, v. 157, p. 105-142.

28.Hopp T.P., Woods KR. Proc. Nat. Acad. Sei. USA, 1981, v. 78, p. 3824-3828.

29. Sotomatsu-Niwa Т., Ogino A. J. Mol. Struct. (Theochem.), 1997, v. 392, p. 43-54.

30. Stryer L. Biochemistry. 4th. Ed. W.H. Freeman & Co. N.Y., 1995. Part I. Chapter 1.

31 .DoeblerJA. Cell Biol. Toxicol., 1999, v. 15, p. 279-289.

32. Rokitskaya T.I., Kotova E.A., Antonenko Y.N. Biophys. J., 2002, v. 82, p. 865-873.

33. Williams R.J.P. J. Theor. Biol., 2002, v. 219, p. 389-396.

34. Pascal S.M., Cross T.A. J. Mol. Biol., 1994, v. 241, p .431-439.

35. Groot B.L., Tieleman D.P., Pohl P., Grubmuller H. Biophys. J., 2002, v. 82, p. 2934-2942.

36. Lindahl E., Hess В., van der Spoel D. J. Mol. Mod., 2001, v. 7, p. 306-317.

37. Corringer P.J., Le N.N., Changeux J.P. Annu. Rev. Pharmacol. Toxicol., 2000, v. 40, p. 431-458.

38. Langosch D., Thomas L., and Betz H. Proc. Natl. Acad. Sei. USA, 1988, v. 85, p. 7394-7398.

39. Legendre P. Cell. Mol. Life Sciences, 2001, v. 58, p. 760-793.

40. Yushmanov V. E., Mandal P. K, Liu /.. Tang P., Xu Y. Biochem., 2003, v. 42, p. 3989-3995.

41. Keramidas A., Moorhouse A.J., French C.R., Schofleld P.R., Barry P.H. Biophys. J., 2000, v. 79, p. 247-259.

42. Jentsch T.J., Stein V., Weinreich F., Zdebik A.A. Physiol. Rev.,

2002, v. 82, p. 503-568.

43. Miyazawa A., Fujiyoshi Y., Unwin N. Nature, 2003, v. 423, p. 949.

44. Lopez. C.F., Nielsen S.O., Moore P.В., Klein M.L. Proc. Natl. Acad. Sei. USA, 2004, v. 101, p. 4431-4434.

45. Gao //., Kong Y, Cui D., Ozkan C.S. Nano Letters, 2003, v. 3, p. 471-473.

46. Yeh I.-C., Hummer G. Proc. Natl. Acad. Sei. USA, 2004, v. 101, p. 12177-12182.

47. Xie Y.H., Soh A.K. Materials Letters, 2005, v. 59, p. 971-975.

48. Li II, Dowd V., Stewart DJ., Burton S.J., Lowe C.R. Nat. Biotech., 1998, v. 16, p. 190-195.

49. Fersht A. Enzyme structure and mechanism. 2nd. Ed. W.H. Freeman & Co. N.Y., 1985.

50. Финкельштейн A.B., Птицьш О.Б. Физика белка. М.: Книжный дом «Университет», 2002, 376 с.

51. Khalatur P.G., Novikov V.V., Khokhlov A.R. Phys. Rev. E.,

2003, v. 67, p. 051901.

52. Ivanov V.A., Chertovich A. V., Lazutin A.A., Shusharina N.P., Khalatur P.G., Khokhlov A.R Macromol. Symp., 1999, v. 146, p. 259.

53. Chertovich A.V., Ivanov VA., Zavin B.G., Khokhlov A.R. Macromol. Theory Simul., 2002, v. 11, p. 751.

54. Шайтан К.В., Турлей E.B, Голик Д.Н., Терешкина К.Б., Левцова О.В., Федик И.В., Шайтан А. К., Кирпичников М.П. Хим. физика, 2006, т. 25, X» 7.

55. Koronakis V. FEBS Lett., 2003, v. 555, p. 66-71.

56. Levitt M. Ann. Rev. Biochem., 1997, v. 66, p. 549—579.

57. Miranker A.D., Dobson C.M. Curr. Opin. Struc. Biol., 1996, v. 6, p. 31-42.