Разработка алгоритмов для оптимизации пространственных структур РНК

Разработка алгоритмов для оптимизации пространственных структур РНК

Введение

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

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

Для определения пространственной структуры молекулы, на практике чаще всего используется рентгеноструктурный анализ. Однако этот физический метод является довольно длительным и затратным, в связи с чем и возникает потребность в менее трудоёмких методах определения пространственной структуры. В работе [7] был предложен способ моделирования сворачивания молекулы РНК на полноатомном уровне. Такой подход позволяет добиться высокой точности результатов, но требует длительного времени для расчёта необходимой пространственной структуры. Одним из менее точных, но зато более быстрых способов является математическое моделирование третичной структуры молекулы в виде набора тонких упругих стержней, который и будет рассмотрен в данной работе.

Из молекулярной биологии нам известно, что молекулярная цепь РНК (ДНК) составляется из нуклеотидов четырёх видов: A - аденин, G - гуанин, C - цитозин, а также U - урацил. Среди поперечных связей, которые образуются при сворачивании молекулярной цепи в пространственную структуру, наиболее сильными являются водородные Уотсон-Криковские связи, возникающие между основаниями некоторых нуклеотидов (между A и U, между C и G и между G и U).

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

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

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

) Шпилечная петля c) Боковая петля

b) Внутренняя петля d) Многозвенная петля

Рис. 1: типы петель

·   Шпилечная петля: соединяет первую и вторую нить в одном стебле. Это простейшая однозвенная петля (она состоит из одного однонитевого участка). Считается, что шпилечная петля всегда содержит не менее трех нуклеотидов.

·   Внутренняя петля: содержит два однонитевых участка. Длины этих участков являются параметрами, определяющими петлю.

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

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

1. Постановка задачи

.1 Свободное состояние молекулы

Молекулярная цепь рассматривается как тонкий упругий однородный стержень, имеющий в свободном состоянии форму винтовой линии (рис. 2). Так как по сравнению с длиной молекулы (до нескольких сантиметров) её толщина (~20 А) незначительна, то поэтому такое представление можно считать вполне допустимым. Длины участков нуклеотидов, которые составляют молекулярную цепь (расстояния между Уотсон-Криковскими связями), предполагаются равными. С шагом, равным длине одного нуклеотида, в стержень вставлены жёсткие перемычки, равные по длине половине Уотсон-Криковской связи.

Рис.2. Тонкий упругий стержень в свободном состоянии

Геометрические параметры винтовой линии и ориентация перемычек (Уотсон-Криковских связей) выбираются так, что два свободных стержня одинаковой длины, которые правильно расположены в пространстве, образуют двойную спираль в A-форме (рис. 3). Пространственная структура молекулы состоит из стеблей и петель в соответствии с заданной вторичной структурой - все элементы вторичной структуры последовательно нумеруются таким образом, что элемент с меньшим номером располагается ближе к начальной точке винтовой линии. В соответствии с краевыми условиями концы стержней, образующих петли, ориентируются так же, как концы нитей в стебле А-формы РНК.

Рис.3: Двойная спираль в A-форме:) шаг витка; b) малая бороздка; c) большая бороздка

1.2 Исходные параметры задачи

В качестве исходных данных для задачи выступают:

§  геометрические параметры (радиус, шаг витка) винтовой линии третичной структуры для свободного состояния молекулярной цепи РНК;

§  длина и ориентация Уотсон-Криковских связей для свободного состояния молекулярной цепи РНК;

§  количество нуклеотидов в каждом элементе вторичной структуры для напряжённого состояния молекулы РНК;

§  номер и тип (стебель либо один из 4-х видов петель) каждого элемента вторичной структуры напряжённого состояния РНК;

§  механические характеристики для элементов вторичной структуры (изгибные и крутильные жёсткости).

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

Так как молекула РНК, в отличии от молекулы ДНК, может принимать практически произвольные формы, то возникает задача определить истинную третичную форму молекулы РНК - форму, которую принимает молекула, находящаяся в равновесии. Из механики твёрдого тела (см. [6]) известно, что равновесному состоянию тела соответствует состояние, при котором значение упругой энергии тела достигает минимума.

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

1.3 Цель работы

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

Все исходные данные поставленной задачи заносятся в текстовые поля ввода программы, которая осуществляет нахождение оптимальной по упругой энергии пространственной структуры РНК, и вычисление координат всех нуклеотидов для этой структуры. Программа написана на языке C++, интерфейс разработан в среде C++ Builder. Получившиеся в результате значения конвертируются в формат программы-визуализатора DS Visualizer 1.7, которая затем строит графическое представление новой пространственной формы молекулы.

2. Математическое моделирование третичной структуры напряжённого состояния РНК

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

Для этого рассмотрим процесс деформации молекулы:

Рис.4. Свободное состояние молекулы

Под действием Уотсон-Криковских связей часть нуклеотидов спаривается, образуя таким образом стебель в А-форме. При этом каждые две "парные" перемычки объединяются в одну Уотсон-Криковскую связь (для рисунка 4: происходит объединение 1 с 19, 2 с 18 и так далее, вплоть до пары 6-14). Остальные же нуклеотиды образуют шпилечную петлю (в неё входят нуклеотиды с 7 по 13 включительно). В итоге - образуется новая третичная структура для молекулярной цепи РНК, соответствующая напряжённому состоянию. Его вторичная структура выглядит как:

Рис. 5. Вторичная структура молекулы в напряжённом состоянии

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

Рис. 6. Оптимальная третичная структура РНК

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

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

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

При параметризации кривой своей длиной, мы будем отсчитывать её от точки, соответствующей координатам того нуклеотида, которому мы дали первый номер (как на рисунках 4-6).

Для каждой из таких кривых (для петель) в дальнейшем будут задаваться краевые условия - исходя из условий на координаты, соответствующих начальной и конечной точкам петли. То есть, оптимальная по упругой энергии пространственная форма каждой петли вычисляется с таким расчётом, чтобы её начальная и конечная точки совпали с точками другого элемента вторичной структуры, на которые и "крепится" эта петля (под "элементом вторичной структуры" в данной работе будем понимать либо каждый однонитевый участок каждой петли, либо каждый стебель). Для указанной на рисунке 6 шпилечной петли это будут точки 6 и 14, расстояние между которыми равно длине одной Уотсон-Криковской связей.

3. Упругая энергия для одного элемента вторичной структуры напряжённого состояния РНК

Рассмотрим тонкий упругий однородный стержень длины L, который находится в равновесии, а осевая линия которого принимает форму некоторой плоской кривой  в трёхмерном пространстве:

x = x(s); y = y(s); z = z(s);

где  - натуральный параметр на кривой, т.е. длина дуги кривой, измеряемая от начальной точки до текущей точки (x(s), y(s), z(s)). Произведём малые деформации стержня, что позволит вывести формулу расчёта функционала упругой энергии для стержня общего вида из аксиом линейной теории упругости (см [2]).

3.1 Аксиомы линейной теории упругости

·        Упругая энергия W стержня потенциальна - она зависит только от его формы и не зависит от последовательности деформаций, которые привели его к этой форме.

·        Упругая энергия W аддитивна (то есть упругая энергия всего стержня равна сумме упругих энергий его частей), и для каждой формы стержня существует неотрицательная функция  (плотность упругой энергии) такая, что .

·        Для тонкого упругого стержня общего вида выполняются следующие уравнения: существует функция , независимая от формы стержня, такая, что

,

где k(s) - кривизна кривой  в точке (x(s), y(s), z(s)) в напряжённом состоянии, а k₀(s) - кривизна стержня в этой же точке в свободном состоянии.

Тогда полная упругая энергия W вычисляется следующим образом:

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

В проекциях на координатные оси:

где p, q, r и p₀, q₀, r₀ - функции, зависящие от s, а именно:

, , ,

Таким образом:

где x', y', z' относятся к новой винтовой линии, соответствующей напряжённому состоянию молекулы.

.2 Общая формула функционала упругой энергии

Теперь вычислим упругую энергию уже не просто для стержня, а для элемента вторичной структуры молекулы РНК, в поставленной выше математической модели. Так как свободное состояние молекулы в ней моделируется винтовой линией, то для неё (см. [3]):

где a и b - геометрические параметры винтовой линии (радиус и подъём винтовой линии), s - расстояние от начала винтовой линии свободного состояния молекулы до данной точки.

Исходя из того, что:

где L₁ - длина предыдущего участка винтовой линии свободного состояния молекулы (от начала этой винтовой линии до точки, соответствующей началу данного элемента вторичной структуры), а s - уже расстояние от начальной точки элемента вторичной структуры до текущей точки, функционал упругой энергии примет вид:

где  - кривизна винтовой линии свободного состояния молекулы (винтовая линия обладает свойством постоянства значения кривизны в каждой своей точке), L₁ - длина предыдущего участка винтовой линии свободного состояния молекулы.

3.3 Формула функционала упругой энергии для поставленной модели:

Будем вычислять упругую энергию для построенной математической модели пространственной структуры напряжённого состояния для элемента вторичной структуры РНК (кривой в трёхмерном пространстве). Эта задача осложняется тем, что для натуральной параметризации кривой должно выполняться тождество [см. 3]:

(x'(s))² + (y'(s))² + (z'(s))² = 1 (для всех s)

Поэтому поступим так. Введём функции α(s) и β(s):

задаваемые полиномами 3-й степени, такие, что:

Тогда, очевидно, условие натуральной параметризации кривой всегда будет выполнено. То есть:

Поэтому:

Тогда функционал упругой энергии будет иметь вид:

Вводя вспомогательные функции , ,  и  такие, что

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

В итоге, функционал упругой энергии примет окончательный вид:

В ходе минимизации этого функционала, значения содержащихся в нём интегралов будут вычисляться с помощью численного метода Симпсона (см. [5]), так как они не вычисляются аналитически. Таким образом, мы получили функционал от 4+4=8 неизвестных, представляющий собой выражение для упругой энергии выбранного элемента вторичной структуры РНК.

4. Расчёт упругой энергии для всех элементов молекулярной цепи

.1 Общие положения

Так как упругая энергия потенциальна и аддитивна, то упругая энергия всей молекулярной цепи равна сумме упругих энергий всех элементов вторичной структуры молекулы.

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

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

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

4.2 Задание граничных условий

Остаётся только определить значения кривизн в точках перегиба, которые бы соответствовали истинной пространственной форме молекулы. Для этого нам необходимо задать условия на гладкость кривой, соответствующей молекулярной цепи РНК, в этих точках, так как в реальности молекулярной цепи соответствует гладкая и неразрывная кривая. Эта задача осложняется тем, что каждая точка перегиба принадлежит одновременно двум элементам, в каждом из которых уравнение кривой будет иметь свой собственный вид (и соответственно, коэффициенты α₁, α₂, α₃, α₄ и β₁, β₂, β₃, β₄ в формуле функционала упругой энергии будут, скорее всего, различными).

Рассмотрим данную ситуацию на примере соединения двух произвольных кривых:

Рис. 7. Угол между касательными к каждой прямой в точке перегиба

Здесь у нас есть две различные натурально параметризованные кривые, с заданными направлениями отсчёта натурального параметра (то есть направлениями обхода кривой). Точка M на рисунке 7 является точкой соединения этих кривых. Очевидно, что для того, чтобы одна кривая гладко «входила» в другую, необходимо, чтобы угол θ между касательными к каждой кривой был равен нулю. Обозначим φ₁, φ₂ - углы наклона касательных к оси абсцисс. Очевидно, что θ = φ₂ - φ₁. А так как (см. [2]) кривизна k(s) = , то получаем, что k₂(s) - k₁(s) =  = 0. Таким образом, одним из граничных условий, задаваемых при минимизации функционала упругой энергии, будет равенство значений кривизн на концах одних конечных элементов и на началах других (несмотря на то, что эти элементы будут задаваться различными функциями).

5. Минимизация функционала упругой энергии для всей молекулярной цепи

Задача минимизации функционала упругой энергии для одного конечного элемента цепи РНК состоит в нахождении таких неизвестных α_i и β_i, что построенный в пункте 3.3 функционал упругой энергии W достигает своего минимума. Таким образом, используя метод конечных элементов (см. [6]), мы можем провести минимизацию функционала упругой энергии для всей молекулярной цепи.

Будем поочерёдно производить минимизацию каждого конечного элемента (об их задании сказано в пункте 4.1), начиная от элемента, содержащего в себе точку-начало отсчёта натурального параметра s (для всей молекулярной цепи). Для этого при минимизации построенного в пункте 3.3 функционала упругой энергии W будем использовать метод наискорейшего спуска (см. [4]), но также нам необходимо учесть граничные условия, возникающие в «точках перегиба»:

непрерывность молекулярной цепи: координаты X, Y и Z начальной точки (N+1)-го элемента должна совпадать с координатами конечной точки N-го элемента;

гладкость молекулярной цепи: значения кривизн на конце N-го элемента и на начале (N+1)-го элемента должны совпадать.

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

где x^N-1, y^N-1, z^N-1 - полученные на предыдущем (то есть (N-1)-ом) шаге значения координат, s’ соответствует «последней» точке конечного элемента. В качестве стартовой точки для метода Лагранжа возьмём точку {1;1;...;1}. Таким образом, мы будем проводить минимизацию для всех конечных элементов, пока не дойдём до последнего. У каждого из них будут свои параметры α_i и β_i, вычислив которые, мы сможем определить координаты x, y, z для любой точки в молекулярной цепи как:

Эти интегралы, как и интегралы в ограничениях из метода Лагранжа, не вычисляются аналитически, но их значения в заданных точках можно получить с помощью численных методов. В данной работе был использован метод Симпсона (см. [4]). Для вычисления антиградиента в методе наискорейшего спуска был использован метод неопределённых коэффициентов (см. [4]).

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

Рис. 8: Нахождение координат перемычек

Исходя из входных параметров задачи, а также значений, полученных в результате минимизации упругой энергии, для конструкции на рисунке 8 нам будут известны координаты точек A, N и B - они являются координатами соответствующих нуклеотидов, а также значения углов φ и ψ и длина перемычки NM (из входных параметров). Тогда значение точки M очевидным образом находится из формул для скалярного произведения векторов NA и NM, а также NB и NM. Таким образом, мы вычисляем значения координат концов всех перемычек, соответствующих Уотсон-Криковским связям, и подставляем их в программу-конвертер в формат PDB.

6. Расчёт третичной структуры РНК со шпилечными и многозвенными петлями

Определим оптимальную по упругой энергии пространственную структуру молекулы РНК, вторичной структурой которой является шпилечная петля. Зададим количество нуклеотидов в каждом конечном элементе вторичной структуры РНК (23 всего, из них 7 в петле), длины Уотсон-Криковских связей и длину одного нуклеотида в цепи, механические параметры конечных элементов свободного состояния молекулы, а также характеристики свободного состояния молекулы (радиус и шаг витка винтовой линии).

После применения построенного алгоритма минимизации и конвертации полученных значений в формат PDB, полученная в результате ~10 минут вычислений третичная структура выглядит так (вид с разных ракурсов):

Рис. 9 Искомая третичная структура молекулы

Теперь произведём расчёты для молекулы тРНК Phenilalanine Yeast, содержащей 72 нуклеотида, пространственная структура которой была найдена в работе [1] (стр. 29). Она содержит три шпилечные и одну многозвенную петли, её вторичная структура имеет вид:

Рис. 10 : Вторичная структура молекулы (всего 72 нуклеотида)

Рис. 11а. Искомая третичная структура молекулы

Рис. 11б Искомая третичная структура молекулы

Для сравнения, в работе [1] для данной молекулы тРНК была получена такая пространственная структура:

Рис. 12. Найденная ранее третичная структура молекулы

Заключение

В рамках данного проекта была произведена следующая работа:

·        изучены основные принципы строения молекулярной цепи РНК;

·        построена математическая модель свободного и напряжённого состояний молекулы РНК;

·        разработан и реализован алгоритм расчёта упругой энергии напряжённого состояния молекулы РНК;

·        разработан и реализован алгоритм минимизации функционала упругой энергии напряжённого состояния молекулы РНК;

·        реализован конвертер, переводящий полученные данные в формат PDB;

·        построены пространственные структуры для шпилечной петли и тРНК Phenilalanine Yeast.

В результате работы были разработаны программные средства для нахождения координат нуклеотидов, входящих в состав молекулы РНК с заданной вторичной структурой, и для преобразования полученных данных в формат Protein Data Bank.

Литература

молекулярный цепь нуклеотид

1.Н.Н. Козлов, Е.И. Кугушев, Д.И. Сабитов, Т.М. Энеев: «Компьютерный анализ процессов структурообразования нуклеиновых кислот» (Институт прикладной математики имени М.В. Келдыша Российской Академии Наук, Москва, 2002)

.Е.И. Кугушев, Д.И. Сабитов: «О тонких упругих стержнях с быстро меняющимися периодическими характеристиками» (Институт прикладной математики имени М.В. Келдыша Российской Академии Наук, Москва, 2005)

.Ю.А. Аминов: «Дифференциальная геометрия и топология кривых» (изд. «Наука», Москва, 1987)

.И.В Бейко, Б.Н. Бублик: «Методы и алгоритмы решения задач оптимизации» (изд. «Высшая школа», Киев, 1983)

.А.Г. Трифонов: «Постановка задачи оптимизации и численные методы ее решения» (изд. «Дело», Москва, 2003)

.И.Ф. Образцов: «Метод конечных элементов в задачах строительной механики летательных аппаратов» (изд. «Высшая школа», Москва, 1985)

.Rhiju Das, David Baker: «Automated de novo prediction of native-like RNA tertiary structures» (University of California, Berkeley, CA, 2007)