Диссертация (1103763), страница 8
Текст из файла (страница 8)
17. Энергетические потенциалы взаимодействия в стенке микротрубочки. А. Потенциалпоперечного взаимодействие между димерами. Б. Потенциал продольного взаимодействия.2.1.3Изгибные взаимодействия.Изгиб протофиламентов в плоскости их движения описывается с помощью квадратичногопотенциала:g kbending,nB( k ,n T ( D ) ) 22(16)45Здесь В – параметр жесткости связи на изгиб, k ,n i i 1 - угол между двумя соседнимимономерами в одном протофиламенте (Рис.16 Г,Д). Параметр T ( D ) соответствуетравновесному углу между двумя мономерами и равен 0,2 рад для ГДФ состояния и 0,0рад для ГТФ [23,74]. Значение параметра B было выбрано так, чтобы в выпрямленномсостоянии в одном ГДФ-димере была запасена энергия гидролиза одной молекулы ГТФ: В= 300 kT/rad2 [52].2.1.4Нахождение полной энергии системыПолная энергия системы описывается как сумма энергий поперечного, продольного,изгибного взаимодействий:13 K nU total vklat,n vkint,ner vkint,nra g kbending,n(17)n 1 i 1Суммирование происходит по всем 13 протофиламентам (индекс n) и по каждомумономеру в протофиламенте (индекс i).2.1.5Вычисление градиента энергииПроизводная полной энергии системы по независимой переменнойследующим образом:вычислялась(18)Эта формула включает слагаемое, ответственное за поперечное взаимодействие,продольное (между димерами и внутри димеров), а также энергию изгиба протофиламента(уравнения (13)-(17)).
Производная вычислялась следующим образом:46.(19),Производная продольной энергии взаимодействия между димерами вычисляласьаналогично производной поперечной энергии.2.2Алгоритм вычисленийИспользованный в программе подход состоит из двух составляющих – динамической икинетической. Динамическая составляющая основана на броуновской динамике.Использование алгоритма броуновской динамики [112] в моделировании микротрубочкистало возможным с развитием технологий высокопроизводительных вычислений ипоявлением мощных суперкомпьютеров, таких как суперкомпьютерный комплексМосковского Университета.В динамической части алгоритма для нахождения положения частиц системы взаданные моменты времени решается уравнения движения Ланжевена в броуновскомприближении.2.2.1Метод броуновской динамикиДля молекулярно-динамической модели микротрубочки нами был выбран методброуновской динамики. Для взаимодействующих частиц размера белковых молекулпомимо метода Монте-Карло часто используется также метод броуновской динамики[113,114].
Этот альтернативный методу Метрополиса Монте-Карло подход имеет одноважное преимущество. А именно, в методе Метрополиса[108] переходам междусостояниями в явном виде не сопоставлено время, поэтому для калибровки временной47шкалы нужно использовать какую-либо известную из эксперимента зависимость.Например, для микротрубочки такой зависимостью может быть скорость роста илиразборки.
При использовании метода броуновской динамики такая дополнительнаякалибровка не требуется, что в нашем случае является существенным преимуществом, таккак скорости роста и разборки микротрубочки могут быть использованы для калибровкипараметров потенциалов взаимодействия в решетке микротрубочки.Метод броуновской динамики основывается на решении уравнения Ланжевена:md 2qdq F (t ) Frand (t )2dtdt(20)где m – масса частицы, γ – коэффициент вязкости броуновской частицы, q – ее координата,Frand (t ) - случайная сила, F (t ) dU- систематическая сила, действующая на частицу вdqэнергетическом поле U.Свойства случайной силы таковы, что [115]1) средняя по ансамблю частиц случайная сила равна нулю (сбалансированностьфлуктуаций):2) ДляFrand (t ) 0 .временныхинтерваловпродолжительнееслучайная сила не обладает корреляцией:длительностистолкновенияFrand (t1 ) Frand (t 2 ) (t1 t 2 ) , где (t1 t2 ) - функция ДиракаВ формуле (20) пренебрегаем членом, содержащим вторую производную координаты повремени ввиду быстрого затухания ускорения на рассматриваемых временных интервалах[112], и оставляем в левой части только дифференциал координаты:dq dt dU dt Frand dq1 dq2 , dq (21)48dq1 гдеdt dU dq(22)- слагаемое, имеющее смысл приращения координаты, вызванного систематическимисилами,dq2 dtFrand(23)- слагаемое, имеющее смысл стохастического смещения броуновской частицы за малоевремяdt .В случае же если рассматривается достаточно большой интервал времени, которыйсущественно дольше характерного времени одного столкновения частиц раствора сброуновской частицей, мы можем найти среднеквадратичную амплитуду этого смещения,используя формулу Эйнштейна: dq 2 2 Ddt , и выражая коэффициент диффузии по2формуле D k BT( k B -константа Больцмана, T – температура в единицах Кельвина),получим среднеквадратичную амплитуду стохастического члена в уравнении: dq2 22k BTt(24)Здесь вместо бесконечно малого дифференциала dt уже используется приращение t >>τ, где τ – время взаимодействия частиц среды.Сам стохастический член уравнения выражается через свою амплитуду умножением наслучайное число из нормального распределения [112]:dq2 2k BTtN (0,1) ,где N(0,1) – случайное число из нормального распределения.(25)49q qТогдаii 1dq i 1dtt ,(26)где i – номер шага итерации.
И разностная схема для нахождения координат на каждомследующем интервале по времени:q i q i 1 гдеt U totalt 2k BT N (0,1) .iq(27)Δt – шаг по времени; U total - полная энергия системы); r - радиус мономера, -коэффициент вязкости среды. В зависимости от того, рассматривается вращательная илипоступательная составляющая, в формулу подставляется q 6r - поступательныйкоэффициент вязкости частицы, либо 8r 3 - вращательный коэффициентвязкости частицы [116],U total- производная полной энергии системы по переменнойqki ,nqki ,n .Кинетическая часть алгоритма включает в себя проверку события прикрепления игидролиза нового димера тубулина к концу любого протофиламента и происходит синтервалом 1 мс.
Вероятность события прикрепленияp = Kon × ctub × 1 мс,(28)где Kon – константа присоединения, ctub – концентрация тубулина в растворе (см. Таблицу1). Для осуществления данной кинетической схемы каждую миллисекунду времени вмодели с помощью генератора случайных чисел [117] выбирается случайное число изравномерного распределения [0, 1), а затем это число сравнивается с вероятностьюсобытия. Если сгенерированное случайное число меньше заданной вероятности, тособытие считается наступившим.
Координаты вновь присоединенного димера задаютсятаким образом, чтобы его продольная и изгибные связи были ненапряженны: угол междувновь присоединенным и нижележащим димерами равен равновесному (0 рад) ирасстояние между продольными сайтами связывания между этими двумя димерами равнынулю.Вероятность события гидролиза для каждого ГТФ-димера не зависит от концентрации.50Было выбрано самое простое правило ГТФ-гидролиза: гидролиз в ГТФ-димерепроисходил с фиксированной вероятностью вне зависимости от положения и окруженияприсоединенного ГТФ-димера.
При гидролизе равновесное значение угла 0 в формуле(16) меняется от 0,0 рад на 0,2 рад.Начальная конфигурация, с которой начинаются вычисления – это микротрубочкадлиной в 10 слоев димеров, каждый из которых находится в ГТФ-состоянии и никогда неиспытываетпереходвГДФ-состояние.Этаконфигурацияописываеткороткиемикротрубочки-затравки, использованные в экспериментах in vitro [44] для наблюдениядинамической нестабильности микротрубочек.Каждые 0,2 нс в модели обновляются координаты мономеров.
Однако для ускорениярасчетов и экономии вычислительных ресурсов рассчитываются и обновляютсякоординаты только тех мономеров, которые находятся не более чем на 10-15 слоевмономеров от начала микротрубочки, если считать относительно самого близкого кначалу димера без поперечных связей. Такое упрощение было сделано на основании того,что колебания димеров в теле трубке дальше чем на 3-4 слоя от ее конца не оказываютвлияния на связи у конца микротрубочки [52].Координаты и нуклеотидные состояния димеров сохранялись с шагом повремени 14 мсдля анализа скоростей роста и частот катастроф, 2 мс для вычисления скорости разборкимикротрубочки.Вычислительный алгоритм реализован на языке С++ с применением технологиипараллельного программирования OpenMP.
Для визуализации результатов вычислений иобработки выходных данных нами был создан программный пакет на языке MATLAB2012. Для вычислений был задействован главным образом суперкомпьютерный комплексМосковского Государственного Университета им. М.В. Ломоносова [118]. Расчет 1секундымодельноговременизанимаетоколо25суперкомпьютерном комплексе МГУ им. М.В. Ломоносова.часоввычисленийна512.3Определение параметров моделиСписок параметров приводится в Таблице 1.















