Модель и метод моделирования

Модель и метод моделирования. Глобулярные белки можно грубо считать как разновидность НР-сополимеров.

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

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

Прежде всего, определим модель и алгоритм, используемый для моделирования эволюционного процесса. Мы рассмотрим модель сополимерной цепи, которая включает только два типа звеньев, Н гидрофобные и Р гидрофильные или полярные. Для простоты зафиксируют НР состав последовательности как 1 1. Рассмотрим континуальную модель, противоположную широко используемой решёточной модели, так как динамика последней дискретна и релаксация плотной глобулы происходит достаточно медленно.

Время эволюции системы определялись уравнениями Ньютона, которые решались при помощи метода молекулярной динамики МД . Мономерные звенья, связанные упругими связями, образуют линейный НР сополимер длины N. Гамильтониан системы включает в себя только парные взаимодействия и имеет вид N-1 N-2 N N Н еHb rij е е Hev rij Ha rij 1 2еpi2 3.1 i, j i 1 i 1 j i 2 i 1 где Hb потенциал деформации связи, Hev принимает во внимание исключённый объём, Ha характеризует энергию притяжения между мономерными звеньями, и i, j изменяются в диапазоне от 1 до N. Расстояние между звеньями определяется как rij зri - rjз, где ri обозначает расположение радиус- вектора звена в трёхмерном пространстве.

Последнее значение в уравнении 3.1 классическая кинетическая энергия цепи. Исключённый объём между всеми несвязанными звеньями моделировался при помощи отталкивающего Леннард-Джонсовского потенциала. 4e s 12 s 6 1 , rij Ј r0 rij rij 4 Hev rij 3.2 0, rij r0 где s e 1 как для Н так и для Р мономерных звеньев и r0 21 6 s - расстояние обрезки. Параметр e, входящий в уравнение контролирует шкалу энергии, тогда как s определяет расстояние на котором взаимодействуют мономерные звенья.

Квазигармонический потенциал связанных звеньев цепи имеет вид Cb 1 b0 12 2 b 0 6 1 , rij Ј b0 rij rij Hb rij 3.3 Cb 2 exp rij 2 1 rij 2 , rij b0 b0 b0 где b0 и rij равновесная и мгновенная длины связей соответственно. Это квазигармоническое уравнение с константами упругости Сb 1 и Cb 2 в гамильтониане 3.1 описывает взаимодействие связанных между собой звеньев.

Мы использовали комбинированный потенциал вместо обычного квадратичного, так как использование простого гармонического потенциала увеличивает время достижения равновесия. Равновесная длина связи, как и другие величины длин, измерены в еденицах d, типичное значение для реальных полимеров составляет s 5 А . Присваиваем b0 1 и Cb 1 Cb 2 1. Потенциал, описывающий притяжение между несвязанными мономерными звеньями имеет вид eabs 1 rij 2 2 , r0 rij Ј rc rij rc Ha rij 3.4 0, rij rc Параметр eab eНН, eРР, eНР устанавливают глубину минимума функции притяжения невалентно связанных звеньев от расстояния и rc 2.8 расстояние на котором обрезают потенциал, описывающий притяжение.

В уравнении 3.4 параметр описывающий перекрёстные взаимодействия выбирается как eНР eННґeРР 1 2 энергия измерена в единицах kBT . В нашем эксперименте eРР переменный энергетический параметр. При eНН eРР 2 фактически мы имеем гомополимерную глобулу, при этом образуется плотноупакованная конформация.

При eНН eРР 0 не существует притяжения между звеньями. Для эффективного вычисления использовалась достаточно короткая цепь с N 128. Предварительные результаты с N 512 показали, что качественно большие системы сильно не отличаются. При моделировании системы не включали частицы растворителя. Чтобы сымитировать эффект растворителя и эволюцию системы в контакте с термостатом температуры T, в уравнение движения Ланжевена добавляют некоррелированный член. miFi Fi Гr Ri , I 1. 2 N 3.5 где mi 1 масса мономерного звена, Fi -Сri H r - r постоянная сила действующая на звена i, R описывает случайные броуновские силы, действующие на каждое мономерное звено, Г учитывает вязкость растворителя.

Величины R и Г связаны между собой флуктационной-диссипативной теоремой, ?Rai 0 ґRai t с 2ГikBTd t, a x, y, z, и обеспечивает постоянство температуры. Заметим, что если Г включить в уравнение без члена R, то система просто диссипирует. Параметр Г зависит от площади доступной растворителю поверхности SASA . Чтобы найти значение SASA для данной конформации, производится аналитическое вычисление площади поверхности Ai для каждого заданного звена.

Определив Ai можно найти Гi как Гi Г0Аi Amax, где Amax максимальное значение площади поверхности доступной растворителю мономера для изучаемой модели и стандартное значение Г0 принимается равным единице. Весовой коэффициент Ai Amax показывает степень подверженности растворителю мономера i. Если значение SASA для данного мономера равно нулю, то броуновская сила и сила трения равны нулю и уравнение Ланжевена 3.5 редуцируется в уравнение движения Ньютона.

Обычно это происходит, когда звено расположено в ядре глобулы. Наоборот, мономерное звено, находящееся на поверхности глобулы, сильно сольватировано. Это значит, что значение Ai становится близким к Amax, следовательно значение Гi становиться близким Г0. Стандартная температура равна T 1 e kB. При интегрировании уравнения движения выбирается шаг интегрирования ?t sЦm e, использовав численный алгоритм Верлета. 3.2.Модель молекулярной эволюции При моделировании эволюционного процесса используются следующий алгоритм 22 1. Чтобы сконструировать начальную конфигурацию G 0 , нужно сгенерировать цепь притяжение между мономерами отсутствует со случайным распределением Н звеньев.

Эта цепь является начальной для данного рассчёта. 2. Готовится набухший полимерный клубок, присваивая параметры eНН и eРР нулю. 3. Складывание цепи происходит при eНН 2 и данном значении eРР. Эта конформация приходит к равновесию после 4ґ105 шагов интегрирования. 4. Одна половина звеньев, имеющих наибольшую площадь доступной растворителю SASA и одновременно удалённые от центра масс глобулярного ядра, модифицируются в тип Р, остальные звенья, имеющие меньшие значения SASA и находящиеся вблизи центра масс глобулярного ядра превращают в тип Н. Благодаря такому модифицированию последовательность мутирет и переходит в следующую генерацию G G 1. Состав последовательности строго определён, поэтому в цепи N 2 гидрофобных и N 2 гидрофильных звеньев.

Чтобы вычислить статистические свойства данной последовательности, производят вычисление теоретико информационные характеристик. 5. Чтобы вычислить термодинамические и структурные свойства, мы производили усреднение по большому числу шагов интегрирования, t2 4ґ105. 6. Повторяем алгоритм по пунктам 2 5 для последней генерации. Это даёт класс глобулярных сополимеров с другой первичной НР структурой.

В нашем исследовании шаги 2 5 независимо повторяем 20 раз, начиная от различных случайных конформаций и тогда все результаты усредняются по этим расчетам, чтобы статистика была лучше.

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

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

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

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

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

Также возможно образование длинной петли или двух ІхвостовІ.Можно предположить, что на начальном этапе эволюции глобулы образуются длинные петли. Затем одна из длинных петель вырождается в длинный хвост. 26 На рис. 3.1. можно видеть типичную морфологии глобулы белково-подобного сополимера и структуру типа головастик. Для количественного описания этого перехода было предложено ряд характеристик. а в Рис. 3.1. а Мицелоподобная структура, в структура типа ІголовастикІ. 3.3. Методы анализа Рассмотрим поподробнее характеристики, которые были предложены для описания перехода глобула головастик . 1. Среднеквадратичный радиус инерции.

Его можно рассчитать по формуле R2g N-1 ri r0 2 3.6 где N число частиц, r0 радиус вектор центра масс, ri радиус вектор i-й частицы Этот параметр характеризует размер макромолекулы. Косвенным образом может характеризовать форму молекулы. Часто используют для изучению перехода клубок - глобула. 2.Длины хвостов. Под хвостом понимают непрерывный участок Н или Р звеньев, который берёт начало с конца полимерной цепи. Так как в головастике длина хвоста достаточна велика, то вероятно это удачная характеристика для описания этого перехода.

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

Так в первичной структуре можно выделить достаточное число петель, то усреднение будет лучше и характеристика меньше флуктуирует во времени. 4. Размер заархивированного файла характеристику можно объяснить следующим образом. Обозначим Н звенья как 0, а Р звенья - как 1. В результате первичную структуру сополимера можно представить как последовательность единиц и нулей. Такой цифровой код записывается в файл и подвергается архивированию. Размер заархивированного файла Lв, выраженным в байтах и характеризует первичную структуру, в частности, распределение в ней единиц и нулей.

Мы использовали стандартный архиватор GZIP. Определение размера заархивированного файла показано на схеме в приложении 6. Эта характеристика удобна тем, что при вырождении глобулы в головастик значение её резко уменьшается. Это обусловлено увеличением длин петель и хвостов . 5. Индекс Шеннона. I Индекс Шеннона Shannon s index вычисляется по формуле I Nlog2N - Nilog2Ni 3.7 где Ni - количество элементов сорта i, n - количество сортов элемента, N - общее количество элементов, Например в цепочке имеется блоков длиной 1 - 2 штуки N1 2 блоков длиной 2 - 1 штука N2 1 блоков длиной 3 - 2 штуки N3 2 блоков длиной 4 - 1 штука N4 1 блоков длиной 5 - 1 штука N5 1 блоков длиной 6 - 1 штука N6 1 блоков длиной 7 - 1 штука N7 1 всего блоков различной длины - 9 штук N 2 1 2 1 1 1 1 9 , сортов блоков - 7 сортов n 7 Индекс Шеннона равен I 9 log29 - 2 2 log22 7 1 log21 9 3.1699 - 2 2 1 7 1 0 24.5 Индекс Шеннона характеризует информационную энтропию.

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

На рис. 3.2. и 3.3 можно видеть, что в билогарифмических координатах эти зависимости можно апроксимировать прямой. Рис. 3.2. Зависимость размера заархивированного файла от длин Н- и Р- петель. Из рисунков можно видеть, что LB Loop m cons t, где m может меняться от -0.45до 0.7. Рис. 3.3. Зависимость размера заархивированного файла от индекса Шеннона.

Можно видеть, что эта зависимость имеет вид LB In cons t. Параметр n меняется в интервале от 0.63 до 0.7. 4.