Мутации структуры белковоподобного сополимера Компьютерное моделирование (стр. 8 из 10)

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

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

Гамильтониан системы включает в себя только парные взаимодействия и имеет вид:

N-1 N-2 NN

Н = åH_b(r_ij)+ åå[H_ev(r_ij) + H_a(r_ij)] + 1/2åp_i²(3.1)

i,j = i+1 i=1 j=i+2 i=1

где H_b – потенциал деформации связи, H_evпринимает во внимание исключённый объём, H_a – характеризует энергию притяжения между мономерными звеньями, и i, j изменяются в диапазоне от 1 до N. Расстояние между звеньями определяется как r_ij= çr_i - r_jç, где r_i обозначает расположение радиус-вектора звена в трёхмерном пространстве. Последнее значение в уравнении (3.1) – классическая кинетическая энергия цепи.

Исключённый объём между всеми несвязанными звеньями моделировался при помощи отталкивающего Леннард-Джонсовского потенциала.

4es¹² _ s⁶1 , r_ij£r₀

r_ijr_ij4

H_ev(r_ij) = (3.2)

0, r_ij> r₀

где s = e = 1 как для Н так и для Р мономерных звеньев и r₀ = 2^1/6 s - расстояние «обрезки».

Параметр e, входящий в уравнение контролирует шкалу энергии, тогда как s определяет расстояние на котором взаимодействуют мономерные звенья.

Квазигармонический потенциал связанных звеньев цепи имеет вид:

C_b⁽¹⁾b₀¹²_ 2 b₀ ⁶ _{__} 1 , r_ij£b₀

r_ijr_ij

H_b(r_ij) = (3.3)

C_b^{(2 )}expr_ij²_ 1 _ r_ij² , r_ij> b₀

b₀b₀

где b₀ и r_ij – равновесная и мгновенная длины связей соответственно. Это квазигармоническое уравнение с константами упругости С_b⁽¹⁾ и C_b⁽²⁾ в гамильтониане (3.1) описывает взаимодействие связанных между собой звеньев. Мы использовали комбинированный потенциал вместо обычного квадратичного, так как использование простого гармонического потенциала увеличивает время достижения равновесия. Равновесная длина связи, как и другие величины длин, измерены в еденицах d,типичное значение для реальных полимеров составляет s» 5 А°. Присваиваем b₀ = 1 и C_b⁽¹⁾ = C_b⁽²⁾ = 1. Потенциал, описывающий притяжение между несвязанными мономерными звеньями имеет вид:

_ e_abs 1 _ r_ij ²², r₀ < r_ij£ r_c

r_ij r_c

H_a (r_ij) = (3.4)

0, r_ij> r_c

Параметр e_ab (=e_НН, e_РР, e_НР) устанавливают глубину минимума функции притяжения невалентно связанных звеньев от расстояния и r_c = 2.8 – расстояние на котором «обрезают» потенциал, описывающий притяжение. В уравнении (3.4) параметр описывающий перекрёстные взаимодействия выбирается как: e_НР = (e_НН´e_РР)^1/2 (энергия измерена в единицах k_BT). В нашем эксперименте e_РР – переменный энергетический параметр. При e_НН = e_РР = 2 фактически мы имеем гомополимерную глобулу, при этом образуется плотноупакованная конформация. При e_НН = e_РР = 0 не существует притяжения между звеньями. Для эффективного вычисления использовалась достаточно короткая цепь с N = 128. Предварительные результаты с N = 512 показали, что качественно большие системы сильно не отличаются.

При моделировании системы не включали частицы растворителя. Чтобы сымитировать эффект растворителя и эволюцию системы в контакте с термостатом температуры T, в уравнение движения Ланжевена добавляют некоррелированный член.

m_iF_i = F_i –Гr + R_i, I = 1. 2, …, N(3.5)

где m_i = 1 – масса мономерного звена, F_i = -Ñr_iH (r) - (r) – постоянная сила действующая на звена i, R описывает случайные броуновские силы, действующие на каждое мономерное звено, Г – учитывает вязкость растворителя. Величины R и Г связаны между собой флуктационной-диссипативной теоремой,

áR_ai(0)´R_ai(t)ñ = 2Г_ik_BTd(t), a = x, y, z, и обеспечивает постоянство температуры. Заметим, что если Г включить в уравнение без члена R, то система просто диссипирует. Параметр Г зависит от площади доступной растворителю поверхности (SASA). Чтобы найти значение SASAдля данной конформации, производится аналитическое вычисление площади поверхности A_i для каждого заданного звена. Определив A_i можно найти Г_iкак Г_i = Г₀А_i/A_max , где A_max – максимальное значение площади поверхности доступной растворителю мономера для изучаемой модели и стандартное значение Г₀ принимается равным единице. Весовой коэффициент A_i/A_max показывает степень подверженности растворителю мономера i. Если значение SASA для данного мономера равно нулю, то броуновская сила и сила трения равны нулю и уравнение Ланжевена (3.5) редуцируется в уравнение движения Ньютона. Обычно это происходит, когда звено расположено в ядре глобулы. Наоборот, мономерное звено, находящееся на поверхности глобулы, сильно сольватировано. Это значит, что значение A_i становится близким к A_max, следовательно значение Г_i становиться близким Г₀. Стандартная температура равна T = 1 e/k_B. При интегрировании уравнения движения выбирается шаг интегрирования ∆t = sÖm/e , использовав численный алгоритм Верлета.

3.2.Модель молекулярной эволюции

При моделировании эволюционного процесса используются следующий алгоритм:[22]

1. Чтобы сконструировать начальную конфигурацию (G = 0), нужно сгенерировать цепь (притяжение между мономерами отсутствует) со случайным распределением Н звеньев. Эта цепь является начальной для данного рассчёта.

2. Готовится набухший полимерный клубок, присваивая параметры e_НН и e_РР нулю.

3. Складывание цепи происходит при e_НН = 2 и данном значении e_РР. Эта конформация приходит к равновесию после 4´10⁵ шагов интегрирования.

4. Одна половина звеньев, имеющих наибольшую площадь доступной растворителю (SASA) и одновременно удалённые от центра масс глобулярного ядра, модифицируются в тип Р, остальные звенья, имеющие меньшие значения SASA и находящиеся вблизи центра масс глобулярного ядра превращают в тип Н. Благодаря такому модифицированию последовательность мутирет и переходит в следующую генерацию G®G + 1. Состав последовательности строго определён, поэтому в цепи N/2 гидрофобных и N/2 гидрофильных звеньев. Чтобы вычислить статистические свойства данной последовательности, производят вычисление теоретико – информационные характеристик.

5. Чтобы вычислить термодинамические и структурные свойства, мы производили усреднение по большому числу шагов интегрирования, t₂ = 4´10⁵.

6. Повторяем алгоритм по пунктам 2 – 5 для последней генерации. Это даёт класс глобулярных сополимеров с другой первичной НР структурой.