Моделирование генных сетей с помощью обобщенного химико-кинетического метода 1

Лихошвай В.А., Матушкин Ю.Г., Ратушный А.В., Ананько Е.А., Игнатьева Е.В., Подколодная О.А.
Институт цитологии и генетики СО РАН, Новосибирск

Аннотация:

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

Графический интерфейс пользователя для представленных моделей доступен по адресу http://wwwmgs.bionet.nsc.ru/systems/MGL/GeneNet/.

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

Введение

Разработка методов расчета и визуализации динамики генных сетей является важным направлением развития методов сбора, хранения и обработки данных, касающихся генных сетей. В данной статье изложен метод математического моделирования биологических систем и приведены результаты применения метода. Методы формализации должны обеспечивать возможность анализа не только исходных генных сетей и их незначительных модификаций (мутации), но и давать возможность расчета динамики поведения практически любых генетических вариантов, которые можно конструировать на их основе. Иными словами, в модели генной сети должны потенциально содержаться модели многих генных сетей, в том числе и таких, которые могут кардинально отличаться от исходной. Этим требованиям вполне соответствует обобщенный химико-кинетический метод моделирования, разработанный нами ранее и применявшийся для моделирования различных генетических систем: онтогенеза фага $\lambda $, одноклеточный цикл развития $\alpha$-вирусов, вируса гриппа, индукции антивирусного действия интерферона. Этот метод обеспечивает точное и эффективное портретное моделирование закономерностей функционирования генных сетей в силу того, что его компьютерная реализация адекватно отображает фундаментальные свойства биосистем.

Краткое описание обобщенного химико-кинетического метода моделирования

Обобщенный химико-кинетический метод моделирования (ОХКММ) ориентирован на формализованное, в первую очередь портретное, описание функционирования произвольных биосистем. Формализация осуществляется на основе блочного принципа, согласно которому моделируемая система расчленяется на элементарные подсистемы, каждая из которых описывается изолированно от других. Описание элементарных подсистем проводится в терминах элементарных процессов. Элементарные процессы описываются на основе формальных блоков. Формальные блоки однозначно характеризуются упорядоченным списком формальных динамических переменных $\overline {X} $, упорядоченным списком формальных параметров $\overline {P} $ и законом преобразования информации $\overline {F} $ (таблица). В таблице приведено десять формальных блоков. Первый блок описывает формальную бимолекулярную реакцию. Список $\overline {X} $ данного блока состоит из трех формальных переменных $х_{1}$, $х_{2}$ и $х_{3}$, имеющих смысл концентраций продуктов реакции. Список $\overline {P} $ включает два формальных параметра $к_{1}$ и $к_{2}$, имеющих смысл констант скоростей прямой и обратной реакции соответственно. Закон переработки информации $\overline {F} $ в данном блоке задается системой обыкновенных автономных дифференциальных уравнений, которые определяют мгновенные скорости изменения значений концентраций $x_{1}, x_{2}, x_{3}$. В блоках 2-7 законы переработки информации также задаются системами дифференциальных уравнений, которые получены исходя из химико-кинетических закономерностей в условиях идеального перемешивания. Однако в общем случае в рамках ОХКММ не вводится каких либо специальных ограничений на язык описания операторов преобразования информации: это могут быть непрерывные, дискретные, арифметические, логические, стохастические и иные преобразования. Конкретный выбор определяется спецификой элементарного процесса, типом решаемой задачи, средствами анализа, выбранными для решения текущих задач и т.д.

Для примера в таблице приведены формальные блоки 8-10. Блок 8 является арифметическим -- здесь закон переработки информации $\overline {F} $ задается операцией суммирования: $x = y_{1} $+...+ y$_{n}$ (он используется как сервисный блок для выдачи суммарной концентрации различных продуктов, например -- суммарной мРНК различных генов, а также для построения автоматных моделей генных сетей), блок 9 -- условный автомат, который ме-

Таблица 1:
1. Обратимая биомолекулярная реакция:
$x_{1} + \;x_{2} \;{\mathop { \Leftrightarrow} \limits_{k_{2}} ^{k_{1}}} \;x_{3}$		 $ \begin{array}{l} \bar {X} = (x_{1} ,\;x_{2} ,\;x_{3} ),\;\;\;\bar {P}\; = \;(... ...;{\frac{{dx_{2}}} {{dt}}} = \; - \;{\frac{{dx_{3}}} {{dt}}}.\, \\ \end{array}$
2. Необратимая мономолекулярная реакция:
$x\;{\buildrel {k} \over \longrightarrow} \;y_{1} \; + \;y_{2} \; + ... + \;y_{n}\;\;$		 $%%\[ \begin{array}{l} \bar {X} = (x,y_{1} ,\;y_{2} ,\;...\;,\;y_{n} ),\;\;\bar... .....\; = - {\frac{{dy_{n}}} {{dt}}} = \; - k \cdot x,\,\,n \ge 0.\, \end{array}$
3. Конститутивный синтез: $\;{\buildrel {k} \over \longrightarrow} \;x_{1} \; + ... + \;x_{n} \;\; $ $%%\[ \bar {X} = \,(x_{1} ,...,x_{n} ),\;\;\bar {P} = \,(k),\;\;F:\,\;{\frac{{dx_{i}}} {{dt}}} = k,\,\,i = 1,...,n,\,\,n \ge 1. $
4. Обобщенная схема 1 Михаэлиса - Ментен $ \begin{array}{l} \bar {X} = \;(x_{1} ,\;x_{2} ,\;...\;,\;x_{m} ,\;y_{1} ,\;y_... ...k_{d} + \;x_{m} ) - \;x_{1} ...x_{m}}} }\;,\,\,m \ge 2,\,\,n \ge 0. \end{array}$
   
5.Обобщенная схема 2 Михаэлиса - Ментен $%%\[ \begin{array}{l} \overline {X} = (e,x_{1} ,...,x_{m} ,y_{1} ,...,y_{n} ),... ...m}}} }}}{{{\frac{{1}}{{x_{1}}} } + ... + {\frac{{1}}{{x_{m}}} }}}} \end{array}$
6. Обобщенная схема 3 Михаэлиса - Ментен $%%\[ \begin{array}{l} \overline {X} = (e_{1} ,...,e_{l} ,x_{1} ,...,x_{m} ,y_{... ......e_{l} x_{1} ...x_{m}}} {{(s_{1} + x_{1} )...(s_{m} + x_{m} )}}}, \end{array}$
7. "Реакция"
$ x_{1} + ... + x_{m} < = > y_{1} + ... + y_{n} $		 $%%\[ \begin{array}{l} \overline {X} = (x_{1} ,...,x_{m} ,y_{1} ,...,y_{n} ),\,... ... r_{yj} k_{2} y_{1}^{b_{1}} ...y_{n}^{b_{n}} ,\,\,j = 1,...,n. \\ \end{array}$
8. Суммирование:
$ x = x_{1} + ... + x_{n}$		 $%%\[ \overline {X} = (x,x_{1} ,...,x_{n} ),\,\,\overline {F} :\,\,x = x_{1} + ... + x_{n} .\, $
9. Пороговый автомат:
$%%\[ \begin{array}{l} \overline {X} = (x,x_{1} ,x_{2} ),\,\,\overline {F} :\,\,x = x_{1}, \mbox{если}\; x_{1} \ge x_{2}.\,\ \end{array}$
10. Стохастический автомат:
$%%\[ \begin{array}{l} \overline {X} = (x,x_{1},y_{1},y_{2} ),\,\,\overline {F}... ...зке}\; [0,1], \; \par y_2\; \mbox{ --- число из интервала}\; [0,1]. \end{array}$


Основные формальные блоки ОХКММ

няет значение переменной $x$ на $y_{1}$, если значение $y_{1}$ превзошло значение $y_{2}$ (используется для построения автоматных моделей генных сетей) и, наконец, блок 10 является стохастическим автоматом: переменной x присваивается значение х$_{{\rm 1}}$, если случайное значение $y_{1}$ превзошло пороговое значение $y_{2}$ (используется для построения стохастических автоматных моделей генных сетей). Отметим, что приведенный список формальных блоков не исчерпывает всего списка блоков (и не может исчерпать по принципиальным соображениям), которые могут быть использованы при моделировании реальных биосистем.

Процесс моделирования биосистемы включает два этапа: этап описания элементарных процессов и этап конструирования из элементарных процессов моделей конкретных модификаций биосистемы. Описание элементарного процесса состоит в том, что выбирается формальный блок и его формальным переменным и параметрам придается биологический смысл, т.е. они становятся конкретными переменными и константами. Один и тот же формальный блок может быть использован при моделировании многократно. Аналогичным образом при написании программы на языке программирования достаточно высокого уровня можно один раз формально описать процедуру-функцию а потом использовать ее многократно, подставляя новые фактические параметры вместо формальных. Например, блок "бимолекулярная реакция" можно применять при описании любых бимолекулярных реакций. Конкретизация формальных переменных и параметров достигается путем задания бинарных наименований. Бинарное наименование представляет пару "Name{Reduction}". Между "Name" и "Reduction" имеется существенное смысловое различие, аналогичное различию в родовом и видовом названии при классификации организма. "Name" определяет тип (класс, разновидность) биомолекулы: "ген", "мРНК", "белок" и т.д. Через "Reduction" осуществляется выбор конкретной биомолекулы в классе "Name". Например, "мРНК" будет сопровождаться конкретным значением "Reduction" в зависимости от того, мРНК какого гена рассматривается: "мРНК{ген1}", "мРНК{ген2}", "мРНК{ген3}" и т.д.

Еще на стадии описания элементарных процессов каждая формальная переменная (параметр) обязательно наделяется определенным "Name", а "Reduction" может быть либо задано сразу, либо остаться неопределенным. Во втором случае "eduction" является специальным служебным словом (комбинацией символов), необходимым для определения корректного "Reduction" уже на втором этапе - конструирования моделей. Служебными являются те и только те слова (комбинации символов), которые включены в специальный список. Список служебных слов может состоять только из пустого слова {}. И хотя только пара "Name{Reduction}" полностью характеризует переменную или константу в модели, тем не менее пара "Name{}" также несет вполне конкретный биологический смысл. Например, продукт "gene{}" обозначает "ген вообще", "mRNA{}" обозначает "мРНК вообще". В тоже время переменная "gene{а1}" обозначает конкретный "ген" а1, переменная "RNA{pol}" обозначает конкретную РНК-полимеразу и т.д.

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

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

Алгоритмы расчета

Модели, которые строятся из элементарных процессов, описанных в терминах формальных кинетических блоков, (например, в терминах блоков 1-7, представленных в таблице), относятся к классу систем автономных дифференциальных уравнений. В настоящей версии ОХКММ для численного интегрирования моделей применяются два метода: метод Гира и метод дробных шагов Яненко.

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

Математическая модель регуляции биосинтеза холестерина в клетке

Холестерин является незаменимым структурным компонентом клеточных мембран и наружного слоя липопротеинов плазмы крови. Одновременно холестерин является предшественником целого ряда других стероидов, а именно кортикостероидов, половых гормонов, желчных кислот и витамина D. Холестерин синтезируется во многих тканях из ацетил-CoA и большая его часть в плазме крови содержится в составе липопротеинов низкой плотности (ЛНП). Свободный холестерин удаляется из тканей при участии липопротеинов высокой плотности (ЛВП) и транспортируется в печень, где превращается в желчные кислоты. При коронарном атеросклерозе наблюдается высокий уровень соотношения (холестерин ЛНП)/(холестерин ЛВП) в плазме.

Генная сеть регуляции внутриклеточного биосинтеза холестерина изучена на настоящий момент достаточно подробно. Сведения о закономерностях ее функционирования накапливаются в базе данных GeneNet (http://wwwmgs.bionet.nsc.ru/systems/mgl/genenet/). Нами разработана модель динамики функционирования данной генной сети. В модели описаны стадии биосинтеза холестерина и механизмы обмена внутриклеточного холестерина с холестерином крови. Учтены отрицательные обратные связи, благодаря которым на уровне транскрипции происходит контроль скорости синтеза холестерина и ЛНП рецепторов в зависимости от внутриклеточного содержания холестерина. Всего модель содержит 65 элементарных процессов. Модель динамики функционирования генной сети, представленной в базе данных GeneNet (http://wwwmgs.bionet.nsc.ru/systems/mgl/genenet/) содержит 40 продуктов (динамических переменных) и 93 константы. Значения констант ряда элементарных процессов были оценены на основании доступных литературных данных. Значения остальных параметров модели определены в результате численных экспериментов, в которых критерием адекватности служили количественные и качественные характеристики поведения данной системы, известные из литературных источников. Модель позволяет рассчитывать равновесное состояние биосистемы. При отсутствии нового вмешательства, система, под действием обратных отрицательных связей, возвращается в исходное состояние: примерно через 3 часа восстанавливается концентрация свободного холестерина в клетке, а полное восстановление исходного состояния наступает через 10-15 часов.

Недостаток рецепторов ЛНП приводит к значительному повышению уровня этих частиц в крови. Исключительно высокий уровень холестерина ЛНП в крови наблюдается при гомозиготной форме гиперхолестеринемии, когда у больного отсутствуют рецепторы ЛНП. Высокая концентрация ЛНП, циркулирующих в крови на протяжении 6 суток у гомозигот (вместо 2,5 суток в норме), облегчает процессы их модификации (химических изменений -- перекисного окисления, десиалирования, гликозилирования и других реакций). Такие изменения имеют важное значение для развития атеросклероза, т.к. именно измененные, а не нативные ЛНП, приобретают атерогенные свойства. На рис. 1 приведено сравнение результатов теоретического расчета нормально функционирующей системы регуляции биосинтеза холестерина и потребления его клеткой с мутантной системой. В последнем случае промоделирована мутация, приводящая к уменьшению скорости наработки рецепторов ЛНП на 30% с вытекающими отсюда последствиями. Ответная реакция клетки на повышение содержания в плазме крови частиц ЛНП на $\sim $300 мг/дл при мутации оказывается менее выраженной, чем в нормальной клетке (рис.1, а, а', в, в'). Поступление в клетки экзогенного холестерина заметно понижается из-за недостатка рецепторов ЛНП. Количество рецепторов ЛНП на поверхности мутантной клетки уменьшается примерно в 2 раза (рис.1, г, г', д, д'). Равновесная концентрация частиц ЛНП в плазме крови увеличивается тоже примерно в 2 раза (рис.1, б, б'), что вполне согласуется с ситуацией, описанной выше.

Рис. 1: Сравнение кинетики изменения количества основных компонентов (а, в-д) нормальной и мутантной (уменьшение скорости наработки рецепторов ЛНП на 30% ) систем биосинтеза холестерина при реакции клетки в ответ на повышение на 30-ой минуте содержания в плазме крови частиц ЛНП (б, б') (расчет по модели). По оси Х - время в часах, по оси Y - число молекул на клетку.

Математическая модель регуляции созревания эритроцита

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

Взаимодействие эритропоэтина с клеточным рецептором приводит к активации транскрипционного фактора GATA-1, ключевого регулятора дифференцировки эритроцитов. GATA-1 стимулирует синтез $\alpha$- и $\beta$-глобинов, а также ферментов, участвующих в биосинтезе гема. Кроме того, GATA-1 активирует свой собственный ген и ген эритропоэтинового рецептора (положительная обратная связь). Гемоглобин, основной компонент зрелого эритроцита, образуется из гема, $\alpha$- и $\beta$-глобинов.

Принципы функционирования данной биосистемы совершенно иные по сравнению с рассмотренной в предыдущем разделе. Биосистема находится в неактивном состоянии до тех пор, пока не будет запущена под воздействием эритропоэтина. После этого программа созревания эритроцита развивается как необратимый процесс, в котором ведущую роль играют положительные связи. В модели функционирование генной сети описано как последовательность событий, которые протекают в клетках созревающего эритроидного ряда. В начале ряда находится эритропоэтин-чувствительная клетка-предшественник, а в конце -- зрелый эритроцит. Для определения значений различных констант модели использовали прямые экспериментальные данные, описанные в литературе. Значения остальных параметров модели определены в результате численных экспериментов, в которых критерием адекватности служили количественные и качественные характеристики поведения данной системы, известные из литературных источников. Из поведения модели прогнозируется, что ряд компонентов системы в ходе созревания эритроцитов имеют осциллирующую динамику наработки, которая объясняется наличием в генной сети интегрального баланса положительных и отрицательных обратных связей. Проведен численный анализ влияния некоторых мутаций на систему. На рис. 2 приведено сравнение динамики наработки гемоглобина при нормальном (рис. 2, ${\it 1}$) и мутантных (рис. 2, ${\it 2,3}$) типах системы дифференцировки эритроидной клетки. Кривая 2 отражает динамику концентрации гемоглобина при мутации гена EKLF при которой полностью прекращается синтез EKLF (фактора, стимулирующего транскрипцию гена $\beta$-глобина). Кривая 3 - то же при мутации гена GATA-1 - в данном случае происходит снижение активности фактора GATA-1 на 30% . Как и следовало ожидать, содержание гемоглобина в зрелом эритроците мутантных систем понижено, так как при данных мутациях снижается скорость продуцирования тех или иных компонентов системы, которые являются предшественниками гемоглобина в эритроидной клетке - предшественнике. Сравнительный анализ двух мутаций показывает, что роль гена GATA-1 важнее роли гена EKLF.

Рис. 2: Сравнение динамики наработки гемоглобина при нормальном (1) и мутантных (2, 3) типах системы дифференцировки эритроидной клетки (расчет по модели). $2$ - Мутация гена EKLF (прекращение синтеза транскрипционного фактора EKLF); $3$ - мутация гена GATA-1 (снижение активности транскрипционного фактора GATA-1 на 30% ). Hb - гемоглобин. По оси Х - время в часах, по оси Y - число молекул на клетку (1 = 10$^{{\rm 8}} $молекул).

Заключение

Таким образом, разработанный нами метод ОХКММ позволяет описывать и исследовать динамику генных сетей, используя естественное разделение системы (структурируемость, иерархичность) как в физическом смысле (на отдельные элементы), так и в информационном (на отдельные уровни организации). С помощью метода построены адекватные динамические модели самых разных систем: онтогенеза бактериофага $\lambda $, системы индукции и противовирусного действия интерферонов, онтогенеза вируса гриппа, созревания эритроцита, регуляции биосинтеза холестерина в клетке, а также ряда других биосистем.

Примечание

... метода 1
Работа получила финансовую поддержку Государственной научно-технической программы "Геном человека" (N 106), Междисциплинарного интеграционного проекта N 65 Сибирского отделения Российской академии наук "Моделирование фундаментальных генетических процессов и систем", Российского фонда фундаментальных исследований (98-04-49479, 00-07-90337).
Ваши комментарии
[SBRAS]
[Головная страница]
[Конференции]
[СО РАН]

© 2001, Сибирское отделение Российской академии наук, Новосибирск
© 2001, Объединенный институт информатики СО РАН, Новосибирск
© 2001, Институт вычислительных технологий СО РАН, Новосибирск
© 2001, Институт систем информатики СО РАН, Новосибирск
© 2001, Институт математики СО РАН, Новосибирск
© 2001, Институт цитологии и генетики СО РАН, Новосибирск
© 2001, Институт вычислительной математики и математической геофизики СО РАН, Новосибирск
© 2001, Новосибирский государственный университет
Дата последней модификации Sunday, 09-Sep-2001 15:40:03 NOVST