Моделирование проводимости ионных каналов на основе методов молекулярной и броуновской динамики (1103919), страница 3
Текст из файла (страница 3)
Наши исследования [2]показали возможное локальное изменение концентрации ионов в примембранной области до 10 раз. Очевидно, что в таком случае эффективные коэффициенты диффузии данных ионов будут меняться. Для учета этого эффекта нашадиффузионная модель была дополнена эмпирическими соотношениями Эйринга и Джонса-Дола для описания концентрационных зависимостей вязкости идиэлектрической проницаемости среды от растворенного электролита: = 01 + (4)где — мольная доля электролита в растворе, и — некоторые константы,подбираемые эмпирически.=1+(︀√ )︀1+ ∑︀(5)Здесь — мольная доля -ого иона в растворе, — ионная сила, выраженнаячерез мольную долю, а = 1 + 2 = 1 + 2 где 1 , 2 определяются из экспериментальных данных, а 1 и 2 для широкогоспектра электролитов принимают постоянные значения 1.44 · 106 и −1389 K−1 .13Следующим важным аспектом является тот факт, что для большинстваразностных схем интегрирования уравнения Ланжевена, применяемых на данный момент (это методы Эйлера (Euler), Хейна (Heun) и неявной средней точки(implicit midpoint) была обнаружена зависимость среднего квадрата скорости⟨︀ 2 ⟩︀⟨︀ ⟩︀ и перемещения 2 частицы [3] от размера временного шага интегрирования.
Подробнее эти вопросы обсуждаются в Главе 2 диссертации. Это является причиной отклонения от равновесной температуры термостата и изменениякоэффициентов диффузии описываемых частиц. Для устранения данного недостатка мы использовали полученную нами разностную схему с учетом скоррелированности стохастических приращений:(︂)︂−Δ2−2Δ1−(1−) = −1 −Δ + + 0,2[︂]︂−Δ1−∆1 − −Δ(6) = −1 + −1+ −+ 2[︂]︂)︂(︂32−2Δ−Δ+ 0, 3 ∆ − + 2−22где нормировочная константа находится из предположения того, чтосистема находится в термодинамическом (тепловом) равновесии:√︂2=а коэффициент корреляции ():2 (1 − − )() = (7)Однако неявный учет растворителя возможен не во всех случаях. Известно, что молекулы воды играют ключевую роль в процессе ионного транспортадля некоторых каналов.
Так, для калиевого канала KcsA участие молекул воды приводит к специфическому механизму проводимости, разделяя в областиселективного фильтра ионы + , что уменьшает характерные величины амплитуды токов. В большей степени это характерно для ионных каналов, обладающих ярко выраженной селективностью к определенным ионам и достаточноузкой проводящей порой, что приводит к явно не диффузионному механизмупроводимости. Поэтому в данной работе мы используем упрощенные методы молекулярной динамики (МД) для моделирования движения заряженных частицв области ионной поры (компартмент II).14Ионная пора, в соответствии с описывающей геометрией, заполнена растворителем в явном виде (Рис.
1). Число молекул растворителя высчитываетсяна основе объема ионной поры и плотности воды в нормальных условиях. После начальной случайной генерации координат молекул растворителя в областипоры проводится оптимизация системы до равновесного состояния. Это происходит путем численного моделирования данной структуры до выравниванияравновесной температуры системы. После этого данная конфигурация сохраняется и используется для дальнейшего моделирования движения иона в областипоры. Необходимо отметить, что молекулы растворителя, покидающие областьII компартмента в ходе моделирования заново генерируются в данной областисо стороны, противоположной стороне выхода молекулы воды.В качестве модели растворителя была выбрана универсальная 3-х частичная модель воды SPC, движение описывалось классической схемой Эйлера,а взаимодействие между частицами определялось потенциалом, состоящим изэлектростатического и парного потенциала взаимодействия типа Леннард-Джонса (2).Помимо описания движения заряженных частиц в области ионной поры методы молекулярной динамики применялись для разрешения следующей проблемы.
Как отмечалось выше, использование методов броуновской динамики длязадач биологического моделирования параметризуется выбором усредненногодействия растворителя на частицу. В классической модели Ланжевена такимпараметром раствора является коэффициент трения , вид которого зависитот используемого приближения.
Однако, если в случае простых частиц (таких,как ионы + , + , − ) подобные приближения могут дать приемлемый результат, то в случае многоатомных молекул ( 4+ , 42− ) ситуация заметноусложняется. Это связано с неоднозначностью выбора аппроксимирующей геометрии, распределения заряда для данных теорий. Поэтому единственной возможностью определить коэффициент диффузии D и трения таких сложныхчастиц—это использование прямого МД моделирования.Нахождение коэффициента диффузии исследуемого вещества производилось двумя способами. В первом случае использовался средний квадрат перемещения частицы:{︃ }︃2∑︁∑︁1⟨[⃗(0) − ⃗()] ⟩ = lim= lim|⃗ ( + 0 ) − ⃗ (0 )|2(8)→∞→∞66 =0 =00где N — число исследуемых частиц в системе, n — число шагов моделирования,15причем = · ∆.Другой метод предполагает построение автокорреляционной функции скорости (VACF) и использование соотношения Грина-Кубо :}︃∞∞ZZ {︃∑︁ ∑︁1 ⃗1=⟨ (0) · ⃗ ()⟩ =⃗ (0 ) · ⃗ ( + 0 ) (9)33 =0 =0000В свою очередь, для определения коэффициента трения броуновской используется автокорреляционная функция стохастической силы (FACF) :∞Z1=⟨⃗ (0 )⃗ ()⟩(10)3 0Использование данных методов позволяет расширить применимость броуновской динамики на молекулы произвольной сложности, устраняя существующие недостатки приближенных теорий, связанные с геометрической аппроксимацией и сложностью расчета коэффициентов трения.
В ходе работы былиполучены характерные макроскопические характеристики основных сложныхионов и нейромедиаторов, участвующих в процессе синаптической передачинервного импульса [6], что позволяет использовать данные частицы в нашейработе. Полученные результаты обсуждаются в Главе 4 диссертации.Разработанная методология комбинированного моделирования применялась для исследования ионной проводимости пуринергических рецепторов P2X2 ,P2X4 и P2X7 рецепторов. Выбор данных ионных каналов в качестве объектаисследования не случаен: среди семейства P2X каналов на сегодняшний деньтолько для P2X4 получена пространственная структура (Hattori, 2012) .
ДляP2X2 ионного канала только недавно была получена 3D реконструкция данныхэлектронной микроскопии (Mio, 2009) , а пространственную организацию P2X7рецептора подобрали по гомологии (Jiang, 2013). Поэтому проведенные в ходеданной работы численные эксперименты призваны не только подтвердить работоспособность методологии, но и внести ясность в пространственную организацию селективных фильтров P2X2 и P2X7 каналов. Ко всему прочему, с данными ионными каналами было проведено достаточно большое количество «patchclamp» измерений (Evans, 1996; Wong, 2000), что представляет собой идеальную экспериментальную базу для имитационного моделирования.
Обсуждениюданных вопросов и полученных результатов посвящена Глава 5 диссертации.Для начала были построены упрощенные геометрии данных каналов. Вкачестве входных данных мы использовали PDB структуру P2X4 в открытом16состоянии (PDB: 4DW1), идентичную геометрию для P2X7 (с отличием тольков аминокислотных остатках в селективном фильтре), а геометрия P2X2 рецептора была получена на основе пространственной реконструкции электронноймикроскопии.Для проверки чувствительности модели к конфигурации селективного фильтра, для каждого ионного канала были определены заряженные или полярныеаминокислотные остатки, которые могут быть повернуты в область ионной поры. Поскольку для P2X4 рецептора известна пространственная организация, томоделирование «нативной» и измененной геометрии канала позволит ответитьна данный вопрос.
Согласно PDB структуре данного канала, в области селективного фильтра отсутствуют какие-либо заряженные или полярные аминокислотные остатки. Однако поворот второго трансмембранного домена (ТМ2) почасовой стрелке на 35∘ развернет остаток аспарагиновой кислоты Asp357 в область поры, что приведет к изменению потенциала вдоль оси канала и к возможному увеличению уровня проводимости. Поскольку для P2X7 конфигурация селективного фильтра получена моделированием по гомологии, то вопрос оего пространственной организации остается открытым, в связи с чем были рассмотрены все 7 различных конфигураций фильтра, образованного остаткамиSer339 , Ser342 и Asp352 .
Анализ аминокислотной последовательности ТМ2 P2X2рецептора выявил возможное участие остатков Thr339 , Ser345 и Asp349 в формировании 7 различных конфигураций.Согласно предлагаемой методологии, аминокислотные остатки и фосфолипиды формируют систему фиксированных зарядов II компартмента. В отдельных работах, посвященных изучению влияния распределения заряда вдоль осиканала на величины проводимости (Boda at al.) заряд аминокислотных остатковаппроксимируют величиной в ±0.5 или ±1 заряда электрона, однако предпосылки, заложенные в таком предположении, остаются далеко не очевидными.Поэтому в данной работе для вычисления распределения заряда по данныммолекулам мы использовали один из наиболее популярных методов квантовойхимии — метод функционала плотности DFT. Данные аминокислоты рассматривались в депротонированном состоянии.
Сначала была проведена предварительная геометрическая оптимизация молекулы в вакууме. Затем электроннаяплотность основного состояния для данной геометрии была спроецирована наположения атомов методом Хиршфельда. Использовался локальный функционал в форме PWC, неограниченные по спину волновые функции, двойнойчисленный базис с учетом валентных орбиталей для атомов водорода. Ради17ус обрезания 3.7 Å.
Проводился полноэлектронный расчет без использованияпсевдопотенциалов. Порог сходимости для геометрической оптимизации составил 10−5 Хартри, 0.005 Å по расстоянию. Относительная точность сходимостина этапе расчета самосогласованного поля 10−6 . Мультипольное разложение —до октуполя. Для улучшения сходимости применялся DIIS. Вычисления проводились с использованием пакета Accelrys Materials Studio, модуль DMol.
Полученный по нашим расчетам заряд рассматриваемых в работе аминокислотныхостатков составил: Asp ( = −0.73e, = 1.7 Å), Ser ( = −0.43e, = 1.2 Å)и Thr ( = −0.50e, = 2.1 Å). Для каждой молекулы определялся размерbounding-сферы, аппроксимирующий полярный/ заряженный участок молекулы, используемый в нашей модели (Рис. 1).Аналогичные вычисления были проведены для молекул основных мембранообразующих фосфолипидов в модели растворителя COSMO. Результаты представлены в Таблице 1. Учитывая, что «patch clamp» измерения проводились на реальных биологических тканях, состав фосфолипидного бислоя,используемый в моделировании соответствовал фосфолипидному составу теланейрона крысы: PC 28 %, PE 20%, SPH 4%, PI 6% и PS 4 %, что соответствует поверхностной плотности заряда на внутренней стороне мембраны ≈ 40мКл/м2 .Состав компартментов I и III соответствовал «patch clamp» протоколам измерений одиночных каналов данного типа.
Для P2X4 : [ ] = 147, [ ] =5, [ ] = 140 (мМ/л), внешний потенциал −150 мВ. P2X2 : [ ] =145, [ ] = 5, внешний потенциал −100 мВ. P2X7 : [6 11 7 ] = 150,[6 11 7 ] = 150, внешний потенциал −110 мВ. Параметры моделирования системы были: размер ячейки 300 Å, температура 298∘ K, шаг броуновскойдинамики ∆ = 0.1 нс, шаг молекулярной динамики ∆ = 1 пс, времямоделирования 1 мкс.Зная число ионов, прошедших через ионный канал за время моделирования системы, легко определить величины проводимости и амплитуды тока .Для чистоты модельного эксперимента концентрации растворенных веществ вкомпартментах I и III были фиксированы: после прохождения через канал ионзаново помещался в свой стартовый компартмент со случайными координатами.В Таблице 2 представлены результаты численного моделирования проводимости ионных каналов P2X2 , P2X4 и P2X7 типа при различных конфигурацияхселективных фильтров.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
