66

66

Государственное образовательное учреждение

высшего профессионального образования

Тверской государственный университет

Химический факультет

Кафедра неорганической и аналитической химии

Направление 510500 - Химия

Новиков Виктор Викторович

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

Магистерская диссертация

Научный руководитель

д.ф.-м.н. проф. Халатур П.Г.

______________________

Работа допущена к защите

« ___ » __________ 2003 г.

Декан химического факультета

д.х.н. профессор Лапина Г.П.

_________________________

Тверь - 2003

План

• 1. Введение
• 2. Литературный обзор
• 2.1. Механическая модель молекулы
• 2.2. Методы компьютерного моделирования полимеров. Метод Монте Карло. Метод Молекулярной динамики
• 2.2.1 Основные подходы к математическому моделированию макромолекул
• 2.2.2. Методы молекулярной динамики (МД)
• 2.2.3. Метод Монте-Карло (МК)
• 2.2.4. Особенности компьютерного эксперимента
• 2.2.5. Решёточные и континуальные модели
• 2.2.6. Трудности машинного эксперимента
• Периодические граничные условия
• 2.2.7. Модернизированные методы компьютерного моделирования
• 2.3. Мотивы укладки цепи в белковых молекулах
• 2.4. Методы анализа белковых структур
• 2.4.1. Определение гомологии первичных структур
• 2.4.2. Нахождение вторичной структуры
• 2.4.3.Метод протягивания
• 2.4.4. Дизайн белковых молекул
• 2.5. Конформационно-зависимый дизайн последовательностей цепи
• 2.5.1. НР - сополимеры, «приспособленные к адсорбции»
• 2.5.2. Молекулярные диспергаторы
• 2.5.3. Моделирование мембранных белков
• 2.7. Белковоподобные сополимеры. Дизайн, структура, свойства
• 2.8. Оптимизация последовательностей белковоподобных сополимеров глобулярных белков
• 2.9. Дальнодействующие корреляции в белковоподобных сополимерах
• 3. Экспериментальная часть
• 3.1. Модель и метод моделирования
• 3.2.Модель молекулярной эволюции
• 3.3. Методы анализа
• 4. Результаты и обсуждение
• 5. Заключение
• 6. Список литературы

1. Введение

Концепция эволюции является одним из краеугольных камней в современных естественных науках: в космологии обсуждается эволюция вселенной, в геологии эволюция Земли, в биологии эволюция живого мира. Эту идею можно использовать и в полимерной науке [1,2]. Соответствующие положения этой проблемы достаточно ясны. В настоящее время биополимеры ( белки, ДНК, РНК) обладают сложной первичной структурой (последовательность мономерных звеньев ) которая определяет их функции и структуру ( в том числе и третичную структуру глобулярных белков). Поэтому эти последовательности ( 20-ти буквенный алфавит в случае белков и 4-х буквенный в случае ДНК и РНК) должны значительно отличаться от случайных и часто проявляют значительные корреляции на различных массштабах. Другими словами, естественно ожидать, что количество информации в таких последовательностях относительно высоко сравнительно высоко по сравнению со случайными (ДНК содержит всю генетическую информацию).

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

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

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

Один из вариантов конформационно-зависимого дизайна сополимеров, который ведёт к статистически достаточно сложным последовательностям рассмотрен в статьях [3,4,5]. Подход заключается в модифицировании поверхности глобулы. Звеньям, находящимся на поверхности присваивается индекс Р (гидрофильные звенья), а находящимся в ядре глобулы - Н (гидрофобные). Конформация полученного сополимера зависит от Н-Н, Н-Р и Р-Р взаимодействий.

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

Были предложены простые различные компьютерные модели, описывающие эволюцию сополимерных последовательностей. Однако, цель большинства этих моделей - решить различные модели физики белка, т.е. проблему конструирования специфической аминокислотной последовательности, которая будет термодинамически стабильна в третичной глобулярной структуре и будет способна быстро сворачиваться в эту конформацию при данной температуре.[6,7] Мы предлагаем другой подход. Именно:

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

2. Мы исследовали применимость теоретико-информационных характеристик для описания молекулярной эволюции.

2. Литературный обзор

2.1. Механическая модель молекулы

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

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

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

В приближении МО ЛКАО ССП энергия молекулярной системы выражается как сумма одно-центровых взаимодействий (), двуцентровых взаимодействий (), трёхцентровых взаимодействий ( ) и т.д.

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

Многочастичные взаимодействия вносят малый вклад в энергию системы, и ими часто пренебрегают.

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

Предполагается, что энергию взаимодействия можно разделить на ряд составляющих.

Энергия деформации связи.

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

U(R) = K(R-R_i)² (2.3)

Где К - силовая постоянная,

R - равновесная длина связи

R_i - мгновенная длина связи

Рассмотрим вопрос о способе описания валентных связей. Здесь существует альтернатива[9]: 1) моделировать связи каким-либо гладким потенциалом (например гармоническим или потенциалом Морзе), 2) учитывать их как абсолютно жёсткие (геометрические) ограничения. Оба способа имеют свои преимущества и недостатки. В первом случае справедливы простые уравнения Ньютона. Однако из-за высокой "жёсткости" обычных валентных потенциалов связанные частицы могут совершать быстрые колебательные движения, учесть которые можно при очень малом шаге интегрирования. Поэтому даже после весьма большого числа шагов эволюции системы будет на самом деле прослежена лишь в течении короткого интервала времени. Проблема решается путём удаления соответствующих степеней свободы при рассмотрении химических связей как постоянных геометрических ограничений. При этом в системе действует единственный невалентный потенциал, который и задаёт минимальный временной масштаб. Однако в подобной постановке движение частиц описывается уже не уравнениями Ньютона, а уравнениями Лагранжа, что усложняет вычислительную схему и делает её более трудоёмкой. Кроме того, абсолютно жёсткие связи не вполне правильно отвечают реальной ситуации.

Оптимальный вариант состоит в подборе таких валентных потенциалов, жёсткость которых мало отличается от невалентных. Один из них, известный как потенциал FENE, имеет вид

-kr₀² ln[1-(r/r₀)] при r < r₀

U_ij(r) = (2.4)

при r > r₀

Здесь r - расстояние между парой связанных частиц i и j, а параметры k и r₀ принимаются равными 10/² и 1.95. Другой удачной аппроксимацией является конструкция из двух потенциалов Леннарда-Джонса, действующих навстречу друг к другу и описывающих только отталкивание.

₀[(b/r)¹² - 2(b/r)⁶ +1 при r < b

U_ij(r) = (2.5)

₀ b¹² _ 2b⁶ +1 при r > b

(2b-r)¹² (2b-r)⁶

Здесь b - заданная длина связи между частицами i и j. Данная функция имеет нулевой минимум при r = b и становится бесконечной при r = 0 и r = 2b.

Энергия деформации углов.

Для каждого атома в молекуле существуют некоторые идеальные углы, отклонение от которых требуют затрат энергии угловых деформаций (Е_угл). Предполагается, что эта энергия аддитивна, причём при малых отклонениях справедлив закон Гука. Аналогично энергии деформации связей можно записать.

U() = C(_i - )² (2.6)

Где С - силовая константа,

- идеальное значение угла связи,

_i - мгновенное значение угла связи

Деформация углов происходит значительно легче, чем деформация связей ( константа К на порядок выше, чем С).

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

Невалентные взаимодействия атомов (Е_с).

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

E_c = f_ij(r) (2.7)

i < j

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

Отталкивание между атомами апроксимируется обычно функциями вида B/rⁿ или Cexp(-Dr), а притяжение - функцией A/r^m (A, B, C, D, m, n - постоянные).

Отсутствие строгого выражения для f(r) вынуждает искать приближённые аналитические формы. Наибольшее распространение получили потенциалы типа "6-12" и "6 - exp".

Потенциал Леннарда - Джонса ( m - n)

f(r) = n(n-m)^-1(n/m)^m/(n-m)[(/r)ⁿ - (/r)^m] =

n(n-m) -¹[(m/n)(r₀/r)ⁿ - (r₀/r)^m],

где r₀ = 2^1/6, m и n - численные коэффициенты. Коэффициент m выбирается равным 6. Величина n лежит в интервале от 10 до 14 и чаще всего принимается равной 12.

В этом случае данный потенциал называют потенциалом "6-12". При m = 6 и n = 12 имеем

f(r) = 4[(/r)¹² - (/r)⁶] (2.8)

Для малых r функция f(r) имеет большое положительное значение, а с ростом r - становится отрицательной, проходя при r = r₀ через минимум ( где f(r₀) = -) и асимптотически приближается к нулю. Значение r при котором функция f(r) пересекает ось r , обозначено . Наличие дальнодействующей части создаёт расчётные трудности, ибо в плотной системе для каждой частицы становится необходимым учёт её взаимодействий с большим числом окружающих частиц. Поэтому функцию 2.8. обрывают на некотором заранее выбранном расстоянии r_c, полагая, что при r > r_c f(r) = 0.

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

В качестве примера представлена подобная функция, первая и вторая части которой сшиваются в точке r = 2^1/6, где f(r) = -, а производная df(r)/dr обращается в нуль при r = 2^1/6 и r = r_c :

f(r) = 4[(/r)¹² - (/r)⁶] при r < 2^1/6

f(r) = 3(r-2^1/6)² _ 2(r-2^1/6)³ - 1 при 2^1/6 < r < r_c (r_c-2^1/6)² (r_c-2^1/6)³ (2.9)

0 при r > r_c

Ван-дер-ваальсовые взаимодействия имеют тот же порядок, что и тепловое движение атомов (кТ). Если ещё учесть, что число этих взаимодействий равно N(N-1)/2, то можно сделать вывод, что при согласованном действии этих сил могут происходить значительные крупномасштабные конформационные перестройки макромолекулы.

Торсионная энергия.

Потенциалы невалентных взаимодействий дают слишком малые значения барьеров внутреннего вращения.

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

U() = U₀/2(1 - cos(n)) (2.10)

Где U₀ - барьер внутреннего вращения. Для карбоцепных полимеров этот барьер (между транс- и гош- состоянием) составляет примерно 3 ккал/моль

n - порядок вращения. Для карбоцепных полимеров n =3.

-значение двухгранного угла.

Электростатические взаимодействия (Е_эл).

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

E_эл = Kq_iq_j/r (2.11)

i < j

Где q_i и q_j - парциальные заряды ( в долях заряда электрона) на атомах i и j, разделённых расстоянием r; - эффективная диэлектрическая постоянная среды; K = 332 - переводной множитель, позволяющий выразить E_эл в ккал/моль, если r - в A.

Парциальные заряды получаются из дипольных моментов связей

q = 0,208/l (2.12)

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

Парциальные заряды на связанных атомах A и B могут быть также получены по приближённой формуле Смита

q = 0,16(X_A - X_B) + 0.035(X_A - X_B)² (2.13)

Здесь X - относительные электроотрицательности атомов ( по Полингу).

Если расстояние, разделяющее два точечных заряда, превышает толщину мономолекулярного слоя растворителя, в котором находятся рассматриваемая молекула, то коэффициент в 2.11 близок к диэлектрической постоянной этого растворителя. В противном случае 1- 4. Чаще всего принимают = .

Электростатические взаимодействия могут проявляются в полиэлектролитах. За счёт ионизации ионогенных групп на полимерной цепочке возникают заряды, а следовательно электростатические взаимодействия. В этом случае эти взаимодействия являются дальнодействующими ( убывают пропорционально 1/r) и определяют значительные крупномасштабные конформационные перестройки.

Другие составляющие энергии молекулы.

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

Взаимодействие двух неподелённых электронных пар (принадлежащим атомам O, S, N и др.) описываются эмпирическим потенциалом, аналогичным (2.11). Неподелённая пара представляется точечным зарядом (q = 0.1 e), расположенным на расстоянии l от ядра. Для атомов O, S,.. расстояние l отсчитывается в направлении продолжения биссектрисы угла с вершиной при данном атоме. Для атома N направление неподелённой пары совпадает с продолжением нормали к плоскости, отсекающей равные отрезки связей при атоме N. Величина l определяется сортом атома.. Если расстояние r, разделяющее две неподелённые пары, превышает 3 А, то параметр в (2.11) равен единице; в противном случае = (3.5r² - 21r + 32.5) -¹

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

2.2. Методы компьютерного моделирования полимеров. Метод Монте Карло. Метод Молекулярной динамики

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

2.2.1 Основные подходы к математическому моделированию макромолекул

Основная задача статистической теории - вычисление средних значений различных величин, которые характеризуют поведение системы в состоянии равновесия. Существуют два подхода к решению этой общей задачи. В первом случае среднее значение <А> некоторого свойства A( r, v ), которое предполагается зависящим от совокупности координат- {r} и скоростей {v} частиц, определяют путем усреднения множества "мгновенных" значений A[r(t),v(t)], наблюдаемых в последовательные моменты времени t на достаточно протяженном интервале :[11]

А = (2.14)

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

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

Такой подход называют усреднением по ансамблю.

Вероятность (или частота) возникновения того или иного состояния пропорциональна его статистическому весу =e^-U/kT, где U - потенциальная энергия данной конфигурации, k - константа Больцмана, Т - абсолютная температура. В этом случае наблюдаемые средние значения даются общим выражением

А =?…? A(r)(r)dr / ?…?(r)dr (2.15)

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

2.2.2. Методы молекулярной динамики (МД)

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

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

m_i a_i(t) = f_i(t) (масса ускорение = сила). (2.16)

Здесь f_i(t) - сила, действующая в данный момент времени t на частицу i со стороны всех остальных частиц системы (эта сила связана с потенциальной энергией известным соотношением f_i(t) = -U/r_i); ускорение определяется как а_i(t) = dv_i(t)/dt или a_i(t) = d²r_i/dt². Если эти производные заменить их конечно-разностными аналогами, то систему уравнений Ньютона, записанных для всех частиц, можно решить на компьютере. То есть, зная координаты частиц r(t) и отвечающие им силы f(t) в некоторый момент времени t, можно через небольшой промежуток времени t найти новые координаты r(t+t) и силы f(t+t) в следующий момент времени t+t и т.д., шаг за шагом. Очевидно, что скорости оцениваются как v [r(t+t)-r(t)]/t. Вычисляя на каждом шаге интересующий нас параметр А можно проследить его эволюцию во времени, а усреднив по достаточно большому числу сделанных шагов s получаем искомые равновесные свойства. Такую схему расчета принято называть численным экспериментом динамического типа или просто методом молекулярной динамики (МД). Используются также различные вариации метода МД, в которых наряду с "внутренними" силами, обусловленными взаимодействием атомов друг с другом, включаются те или иные внешние силы. Подобные схемы моделирования составляют группу методов неравновесной молекулярной динамики.

2.2.3. Метод Монте-Карло (МК)

Вычислительную схему, в основе которой лежит альтернативный (вероятностный) принцип определения средних значения, называют методом статистических испытаний или методом Монте-Карло (МК). В этом методе переходы между состояниями системы осуществляются следующим образом. На каждом шаге случайным образом выбирается частица (или группа частиц) и перемещается на небольшое расстояние в случайном направлении. Это приводит к изменению потенциальной энергии системы на некоторую величину U, которая и определяет вероятность перехода р ~ е^-U/kT из "старого" в "новое" состояние системы. Интересующие характеристики вычисляются на каждом шаге и усредняются по большому числу сделанных шагов.

2.2.4. Особенности компьютерного эксперимента

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

Компьютерная имитация методами молекулярной динамики или Монте-Карло модели физической системы с целью изучения ее характеристик в зависимости от заданных параметров представляет собой численный (компьютерный) эксперимент с этой моделью

Как уже упоминалось, что число частиц при моделировании методами Монте-Карло и молекулярной динамики с помощью современных суперкомпьютеров может достигать колоссальных величин. Даже без суперкомпьютеров достаточно типичны численные эксперименты для значений N порядка десятков и сотен тысяч. Примеры успешного применения методов Монте-Карло и молекулярной динамики для моделирования равновесных составов смесей при постоянном давлении, фазовых равновесий, адсорбции на поверхности твердых тел, свойств жидкостей в микропорах и т.д. достаточно многочисленны. Этими же методами решаются задачи поиска устойчивых конформаций (поворотных изомеров) полимерных молекул, чрезвычайно важные для биохимических приложений.[12,13]

N,V, U(r₁, r₂, ...,r_n)

Монте-Карло Молекулярная динамика

N, V, T N, V, E

Генератор случайных Решение уравнений динамики

движений F = ma

Отбор с вероятностями P = e^-U/kT Траектории r(t), v(t)

Усреднение Усреднение

Равновесные свойства Равновесные и неравновесные

свойства

РИС. 2.2. Схема расчетов методами Монте-Карло и молекулярной динамики.

2.2.5. Решёточные и континуальные модели

На основе решёточной модели выполнено множество теоретических построений, в частности связанных с решением классической и, в каком то смысле, основной задачи физикохимии полимеров о влиянии объемных взаимодействий на конформацию и, соответственно, на свойства гибкой полимерной цепи. Под объемными взаимодействиями обычно подразумевают короткодействующие силы отталкивания, которые возникают между удаленными вдоль по цепи звеньями, когда они сближаются в пространстве за счет случайных изгибов макромолекулы [14,15]. В решеточной модели реальную цепь рассматривают как ломаную траекторию, которая проходит через узлы правильной решетки заданного типа: кубической, тетраэдрической и др. Занятые узлы решетки соответствуют полимерным звеньям (мономерам), а соединяющие их отрезки - химическим связям в скелете макромолекулы. Запрет самопересечений траектории (или, иными словами, невозможность одновременного попадания двух и более мономеров в один решеточный узел) моделирует объемные взаимодействия (Рис. 2.3.). В методе МК при смещении случайно выбранного звена оно попадает в уже занятый узел, то такая новая конформация отбрасывается и уже не учитывается в вычислении интересующих параметров системы. Различные расположения цепи на решетке соответствуют конформациям полимерной цепи. По ним и проводится усреднение требуемых характеристик, например расстояния между концами цепи R.

Исследование такой модели позволяет понять, как объемные взаимодействия влияют на зависимость среднеквадратичной величины <R²> от числа звеньев в цепи N. Конечно величина <R²>, определяющая средние размеры полимерного клубка, играет основную роль в разных теоретических построениях и может быть измерена на опыте; однако до сих пор не существует точной аналитической формулы для расчета зависимости <R²> от N при наличии объемных взаимодействий. Можно также ввести дополнительно энергию притяжения между теми парами звеньев, которые попали в соседствующие узлы решетки. Варьируя эту энергию в компьютерном эксперименте, удается, в частности, исследовать интересное явление, называемое переходом "клубок -- глобула", когда за счет сил внутримолекулярного притяжения развернутый полимерный клубок сжимается и превращается в компактную структуру - глобулу, напоминающую жидкую микроскопическую каплю. Понимание деталей такого перехода важно для развития наиболее общих представлений о ходе биологической эволюции, приведшей к возникновению глобулярных белков [16].

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

66

66

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

2.2.6. Трудности машинного эксперимента. Периодические граничные условия

Любой современный компьютер способен оперировать таким количеством частиц N, которое обычно неизмеримо меньше, чем в реальных макроскопических системах, где N имеет порядок числа Авогадро (~10²³). Пределом технических возможностей наиболее мощных ЭВМ являются совокупности из N~10⁶-10⁷ частиц. Поэтому если речь идет не о замкнутых микрообъемах вещества (например, микрокаплях, то есть кластерах) или об отдельных полимерных молекулах, то необходимо решение вопроса о том, как исходя из моделирования совокупности сравнительно малого числа частиц, интерпретировать свойства макросистемы. Такая проблема решается с помощью специального технического приема, суть которого состоит в том, что из макроскопического объема вещества мысленно вырезается небольшой объем, называемый расчетной (или базовой) ячейкой, а затем отслеживается поведение частиц только внутри этого выделенного объема.

Базовую ячейку определяют как прямоугольный параллелепипед с ребрами длиной L_x, L_y, и L_z, которые ориентированы по осям X, У и Z лабораторной системы координат. Объем ячейки V = L_x·Ly·L_z выбирается таким, чтобы среднечисленная плотность частиц p=N/V равнялась заданной макроскопической плотности. Основной вопрос при конструировании базовой ячейки связан с описанием поведения частиц вблизи ее границ. Нетрудно понять, что если сделать грани ячейки проницаемыми, то в ходе эволюции системы, движущиеся частицы со временем покинут первоначально занимаемый объем. Вместе с тем, в случае непроницаемых границ система, по сути, является малой и конечной, причем значительная доля частиц будет взаимодействовать со стенками (например, при N=1000 и с=1 вблизи стенок кубической ячейки размещено около 60% всех частиц). Поэтому для устранения поверхностных эффектов чаще всего используют так называемые тороидальные или периодические граничные условия (ПГУ). При таком описании противоположные грани ячейки объявляются тождественными. Иначе говоря, производится воображаемое попарное "склеивание" противоположных граней, в результате чего ячейка замыкается сама на себя и преобразуется в некую торообразную фигуру, у которой вовсе нет границ. Поясним суть этого приема на примере двумерной системы.

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

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

2.2.7. Модернизированные методы компьютерного моделирования

Разработаны также имитационные методы, сочетающие принципы динамического и вероятностного (стохастического) моделирования. Они предназначены для изучения влияния растворителя на конформацию и динамические свойства растворенной макромолекулы без явного учета малых молекул растворителя. Такой подход, основанный на численном решении диффузионных уравнений, позволяет расширить вычислительные возможности за счет резкого уменьшения общего числа атомов в системе. В этом случае сила, входящая в правую часть уравнений (2.16.), складывается из обычных потенциальных сил взаимодействия атомов рассматриваемой цепи, сил вязкого трения, пропорциональных скорости движения атомов в среде с заданной вязкостью, а также случайных сил, отражающих хаотичные толчки со стороны молекул среды. Соответствующие расчеты называют моделированием методом стохастической динамики (СД).[17] Поскольку плотная среда изменяет взаимодействие атомов по сравнению с их взаимодействием в вакууме, возникает вопрос о вычислении эффективных потенциалов, которые описывают такое влияние. В настоящее время эта непростая задача чаще всего решается с привлечением специальных представлений статистической теории жидкостей - метода интегральных уравнений. В таком комбинированном ("гибридном") подходе представляющая главный интерес подсистема - макромолекула - моделируется обычными методами, а влияние окружающих её малых молекул учитывается путем численного решения интегральных уравнений. На пути построения комбинированных вычислительных схем возможны и разного рода "логические инверсии": например, для выделенной (не слишком большой) молекулы можно проводить квантово-химические расчеты. а окружающую среду моделировать стандартными методами МД или МК. Таким образом, удается, в частности, исследовать химические реакции, протекающие в конденсированных средах.

2.3. Мотивы укладки цепи в белковых молекулах

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

Рассмотрим некоторые особенности строения глобулярных белков. [18]

1. Глобулярные белки образуют достаточно компактную глобулу. Структура глобулы очень плотная и строго организована. В этом смысле белковая молекула очень похожа на кристалл. Так говоря словами Шрёдингера, белковая структура - это трёхмерный «апериодический кристалл».

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

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

4. Во - вторых возможностью образованию связей с молекулами воды.

5. Это особенность связана с первичной структурой белковой цепи. Плотной и стабильной глобуле соответствует достаточно большое количество первичных структур. Конечно можно придумать первичные структуры, которым будут соответствовать очень стабильная структура ( стабильность выше, чем у других глобулярных структур), но их будет немного. То есть в природе заложена определённая вариативность первичной структуры. Важно, чтобы глобула была не только стабильной, но и ей соответствовало достаточное количество последовательностей цепи.

Отметим ещё некоторые особенности денатурации белковой молекулы.

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

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

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

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

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

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

Понимание общих черт строения белков необходимо для создания адекватных моделей. Основные мотивы строения закладываются в математические модели и служат базой для построения структур белков.

2.4. Методы анализа белковых структур

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

2.4.1. Определение гомологии первичных структур

Этот метод применим не только к белкам, но и к нуклеиновым кислотам. Разработано множество программ, ищущих гомологии. Все они строят выравнивание (alignment) последовательностей, добиваясь наибольшего сходства между ними. При этом за повышение сходства часто приходится платить "разрывом" последовательностей.

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

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

2.4.2. Нахождение вторичной структуры

Зная вклады отдельных взаимодействий в стабильность -спирали, мы можем рассчитать свободную энергию спирализации любого участка цепи, а следовательно -- и Больцмановскую вероятность образования спирали в любом месте полипептидной цепи, еще не свернувшейся в глобулу. Суммируя и усредняя эти вероятности, мы можем рассчитать и среднюю спиральность такого "несвернутого" полипептида. Потом результат можно сравнить с опытными данными -- например, с КД спектрами.

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

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

2.4.3.Метод протягивания

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

Такой подход называют "методом протягивания" (threading method). Он был предложен Б.А.Ревой в 1990 г. и -- независимо, в более простом и более удобном варианте -- Д. Айзенбергом и его группой в 1991 г. Сейчас метод протягивания стал весьма популярным методом опознавания структур "новых" белков по их аналогии со "старыми".

В общем, работа по протягиванию напоминает поиск гомологии, -- только на этот раз "выравниваются" не две первичные структуры, "новая" и "старая", а "новая" первичная структура со "старым" белковым скелетом.

Здесь возникают аналогичные проблемы, как в любом предсказательном методе.

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

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

2.4.4. Дизайн белковых молекул

Задача дизайна -- обратная по отношению к задаче предсказания структуры. Если при предсказании мы должны найти пространственную структуру, наиболее пригодную для рассматриваемой последовательности -- то при дизайне мы должны найти, сконструировать последовательность, годную для создания желаемой пространственной структуры.

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

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

2.5. Конформационно-зависимый дизайн последовательностей цепи

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

2.5.1. НР - сополимеры, «приспособленные к адсорбции»

Такие сополимеры получают по следующей схеме. Гомополимерную цепь, состоящую из звеньев Н, адсорбируем на плоской поверхности. Для этого вводится сильное притяжение между звеньями цепи и поверхностью. Затем звенья, лежащие ближе к поверхности, обозначаем как Н, другие звенья как Р. В компьютерном эксперименте это значит, что _{Н Р}. Таким образом, звенья Н считаются лиофобными, то есть избегают взаимодействия с растворителем, а звенья Р лиофильными. [19] На рис. 2.4 проиллюстрирована схема получения таких сополимеров.

Рис. 2.4. Получение первичной структуры НР -сополимера, «приспособленного к адсорбции».[19]

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

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