Исследование системы уравнений, описывающих кинетику элементарной ферментативной реакции Текст научной статьи по специальности «Математика»

ISSN 0868-5886 НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2014, том 24, № 2, c. 33-42

-ФИЗИКА И ХИМИЯ ПРИБОРОСТРОЕНИЯ — -

УДК 577.036, 541.128

© А. С. Бердников, А. А. Фёдоров, Д. Г. Сочивко, В. Е. Курочкин

ИССЛЕДОВАНИЕ СИСТЕМЫ УРАВНЕНИЙ, ОПИСЫВАЮЩИХ КИНЕТИКУ ЭЛЕМЕНТАРНОЙ ФЕРМЕНТАТИВНОЙ РЕАКЦИИ

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

Кл. сл.: кинетика ферментативных реакций, математическое моделирование

ВВЕДЕНИЕ

ПОСТАНОВКА ЗАДАЧИ

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

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

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

Рассмотрим элементарную реакцию, протекающую под действием фермента. Пусть S — исходный продукт, Р — конечный продукт, Е — фермент. На первой стадии реакции образуется комплекс ES, на второй стадии образуется конечный продукт Р, а молекула фермента возвращается к исходному состоянию:

k+i к+2 E + S ES ^ E + P. k-i

(1)

Если 5 — концентрация вещества S, е — концентрация свободных молекул фермента Е, с — концентрация связанного комплекса ES, р — концентрация конечного продукта, t — время, k+\ — скорость реакции образования комплекса ES, ^ — скорость реакции распада комплекса ES, k+2 — скорость реакции образования продукта Р, то для рассматриваемого процесса можно записать систему дифференциальных уравнений химической кинетики:

ds , ,

— = -k..se + к ,c, dt +1 4

de

— = -k.se + k c + knc, dt +1 -1 +2

— k+e k—iC k+2c,

dt

= k C dt ~ k+2C

Легко проверить, что у системы уравнений (2) имеются два очевидных первых интеграла: IN = e + c = const — полное число молекул фермента, IC = s + c + p — полное число молекул преобразуемого в результате реакции вещества. Этот факт позволяет упростить систему уравнений (2) за счет подстановки e = IN - c , p = IC - s - c :

dS

— = -k+iS(lN - c) + k_xc, dt

f = k+iS(In - c) - c(k-i + k+2), dt

p = IC - s - c.

(3)

Задача состоит в исследовании поведения решения системы уравнений (3) в зависимости от времени и от параметров уравнения.

ПЕРЕХОД К БЕЗРАЗМЕРНЫМ ВЕЛИЧИНАМ

Для того чтобы в системе уравнений (3) отделить значимые параметры от незначимых, перейдем к масштабированным безразмерным величинам. Пусть с(t) = Sc - х(т), 5 ^) = • у (т) ,

р ^) = Sp - г (т), t = ST - т , где х(т), у(т), z(т) — безразмерные концентрации, т — безразмерное время, £с, £5, £р, "Т — пока неизвестные масштабирующие коэффициенты. Тогда система уравнений (3) преобразуется к виду

^ М!) =

"р dт

= -k+lS5INУ (т) + (т) У (т) +

+ k-1Sсх (т),

^ ах(г) = (4)

"р dт

= ^ау (т)- (т) у (т)-

- ^-1 + k+2)ScX(т) ,

SpZ^) = 1с -^у(т)-ScX(т).

Если теперь выбрать масштабирующие коэффициенты так, чтобы выполнялись условия

" = Sc = Sp , k+1INST = 1 , k+1ScST = 1 (т. е. Sc = 1Ы ,

= 1И , = IN, "Т = 1/k+1IN ), и ввести обозначения / = ^ /(k+lIN ) , V = ^/(k+lIN ) , £ = ^ , то

система уравнений (4) запишется в виде:

= - У (т)[1 - х (т)] + ^ х (т),

^ = У (т)[1 - х (т)]-^)-х т), (5)

г ^ ) = е-у (т)-х (т), у(0) = е, х(0) = 0, г(0) = 0.

Тем самым в системе уравнений осталось всего три значимых параметра, от которых зависит решение: л — безразмерная скорость распада комплекса ES на фермент и конечный продукт, V — безразмерная скорость распада комплекса ES на фермент и исходный продукт, е — безразмерная начальная концентрация исходного вещества.

ПРЕДЕЛЬНЫЕ СЛУЧАИ

В предельных случаях ¡л ^ 0 и ¡л ^ да для системы (5) можно найти аналитическое решение.

Рассмотрим случай ¡л ^ да. В этом случае целесообразно искать решение в виде разложения по обратным степеням параметра ¡ :

1

1

У(т) = Уо (т) + -У (т) +—У2 (т)+

И

И

x (т) = x0 (т) +—x1 (т) + ^у x2 (т) + •••,

И

И

z (т) = z0 (т)+ — Z1 (T) + ^ Z2 (т) + -И И2

0.8

0.6

ч:

и

я н

0

V

я я

1 0.4

Я

«

я я

0.2

И

10

2 4 6 8

Врет, отн.ед. Рис. 1. Предельные решения у(т) и г(т) при

/и ^ +<» .

На рисунке показан случай е = 1, однако в общем случае предельное решение просто масштабируется пропорционально значению е

Решением системы (5) является

Уо (т) = еехР(-т) , х0 (т) = ^ ^ (т) = е[1 - ехР(-т)] , у1 (т) = е(е + ут)ехр(-т)-е2 ехр(-2т) , X (т) = еехр(-т),

г1 (т) = -е(1 + е + ут)ехр(-т) + е2 ехр (-2т) , ...

(цепочка разложения легко может быть продолжена и дальше). Предельное решение при р = да равно х(т) = 0, у(т) = ее- , z(т) = е(1 - е-) (граф показан на рис. 1).

ик

2 4 6

Время, отн.ед.

ч

V

X

н

о «

н я

СЗ &

X

V

я к

о

и

0.8

0.6

0.4

0.2

4 6

Время, отн.ед.

ю

4 6 8

Время, отн.ед.

4 6 £

Время, отн.ед.

ю

Кг)

*(т)

4 6 8

Время, отн.ед.

Рис. 2. Предельные решения х(т) и у(т) при ^ ^ 0. а — для случая е = 0.5; б — е = 1.0; в — е = 1.75; г — е = 1.0; е = 2.0; д — е = 2.5

б

а

в

г

д

Рассмотрим случай ¡л = 0. Решением, как легко проверить прямой подстановкой, является

у(t) = р + Я, еХР ч, К> 1 - уехр (^)

х(t) = е-р-Яу еХР( Я) ч, г(t) = 0, Н Г1 - уехр

где Я = , р = ((Я + е)-(1 + v))|2,

у = (е - р)/(е + Я - р) (техника решения соответствующей системы дифференциальных уравнений не принципиальна).

Равновесным значением при t ^ да очевидным образом является ум = р , хм = е- р , где • при е <1

е 2 е

у » V--V -- +----,

1 „ /1 \3 ' 1 - е (1 - е)

е 2 е

х » е - V--+ V •

1 - е

(1 - е)3

при е > 1

1

ум*(е- ^гЧ+ •

е-1 (е -1)

1 -V-

- + V

е-1 (е -1)3 при е = 1

у™=(хя4+7)-v)/2, хм = 1 -(^4+7)-V)/2.

Характеристический показатель X при е Ф 1 ве-

2е

/- Т II I 1 + е 2 дет себя как Я» 11 - е +V--г -V

1 - е

-+•

1 - е

а при

е = 1: .

При V = 0 решением являются функции

х(т) = е 1 -ехр(-(' - Ф) ,

У ' 1 -еехр(-(1 -е)т)

'« = еС-е). г(т)-0

с очевидными предельными значениями: ут = 0, хт = е при е <1 и ум = е-1, хм = 1 при е >1 (на рис. 2, а-д, показаны графики соответствующих решений).

ч щ

Я 1.5 Н

те к

я 1.0

я

Ц. 0.5 Я

к

о

* о

Рис. 3. Численные решения уравнений системы (5) при /и = 1.0 и различных е ~ 1. а — х(т), б — у(т), в — г(т)

(е =0.1, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0; кривые расположены снизу вверх)

1

Б

х

б

в

Интерес представляет поиск решения в предположении, что л и V малы, в виде ряда по малым параметрам л и V. К сожалению, найти решение для соответствующей цепочки систем дифференциальных уравнений в аналитическом виде не удалось.

Во всех остальных случаях система уравнений (5) может быть проинтегрирована только численно.

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ БЕЗРАЗМЕРНОЙ СИСТЕМЫ УРАВНЕНИЙ

Рассмотрим частный случай V = 0 (в этом случае остается только два безразмерных параметра).

Прежде всего следует заметить, что решения системы (5) существенным образом отличаются для случаев е <1 (число молекул реагирующего вещества не превышает числа молекул фермента) и е >> 1 (число молекул реагирующего вещества много больше числа молекул фермента). В первом случае следует ожидать, что всегда имеются несвязанные молекулы фермента и что лимитирующим фактором будет являться соотношение между скоростью связывания исходного вещества и молекулы фермента в комплекс и скоростью распада

а

комплекса на конечный продукт и свободную молекулу фермента. Во втором случае следует ожидать, что лимитирующим фактором будет недостаток свободных молекул фермента.

Рассмотрим случай е <1. На рис. 3, а-в, показаны результаты численного решения системы уравнений (5) при л = 1 и различных значениях параметра е . Наглядно виден рост и последующий спад функции х (т) (причем даже при е = 2 решение не выходит на уровень насыщения х(т) ~ 1). Функция у(т) экспоненциально спадает до нуля, а функция z(т) монотонно растет, выходя на уровень насыщения, равный начальному значению у(т). На рис. 4, а-в, показаны результаты численного решения системы уравнений (5) при фиксированном е = 1 и различных значениях параметра ¡. Видно, что с ростом л не только растет скорость спада функции х(т), но и стремится к нулю ее максимум (а само положение максимума все больше сдвигается к точке т = 0). Соответственно с ростом л увеличивается скорость спада к нулю функции у(т) и скорость выхода на уровень насыщения функции z(т).

б

Рис. 4. Численные решения уравнений системы (5) при е = 1.0 и различных ц . а — х(т), б — у(т), в — z(т)

(¡л= 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0; кривые расположены сверху вниз на а, б и снизу вверх на в)

£=80

4 6 8

Время, отн. ед.

£■=80

40 60

Время, отн. ед.

100

е= 80

40 60

Время, отн. ед.

Рис. 5. Численные решения уравнений системы (5) при ¡и = 1.0 и различных е > 1. а — х(т), б — у(т), в — z(т)

(е =1, 10, 20, 30, 40, 50, 60, 70, 80; кривые расположены слева направо на а, б и снизу вверх на в)

Рассмотрим случай е >> 1. На рис. 5, а-в, показаны результаты численного решения системы уравнений (5) при л = 1 и различных значениях параметра е. Наглядно виден режим насыщения х(т) ~ 1, а также начальный и заключительный переходные процессы. На рис. 6, а-в, показаны результаты численного решения системы уравнений (5) при фиксированном е = 50 и различных значениях параметра л. Из графиков видно, что при увеличении параметра л отрезок времени, в течение которого система пребывает в режиме насыщения х(т) ~ 1, быстро сокращается.

ПРИБЛИЖЕННЫЕ АНАЛИТИЧЕСКИЕ РЕШЕНИЯ

Практический интерес представляет случай, когда концентрация исходного вещества много меньше, чем концентрация фермента. Хотя точное аналитическое решение для системы уравнений (5), по всей видимости, отсутствует, для случая е <1 можно найти приближенное аналитическое выражение, которое будет достаточно точным для практических целей.

Из физических соображений понятно, что если у(0) = е, то все функции, входящие в систему (5), тоже будут иметь величину не больше е при любом т. Поэтому будем искать решение системы (5) в виде ряда по параметру е:

х(т) = ех1 (т) + е2х2 (т) + е3х3 (т) +—, У (т) = еу (т) + е2у2 (т) + е3у3 (т) + - (6)

z(т) = еz1 (т) + е2z2 (т) + е3z3 (т) +—

После подстановки (6) в систему (5) и группировки членов, соответствующих одному и тому же показателю степени ек, получим цепочку линейных уравнений для функций х^т), Ук(т), zk(т). Так, для членов порядка е полученные уравнения имеют вид

х1 = -(и+у) х1 + У^

У1=ух1 - Уl, z1' = ¡х1,

(7)

б

а

в

ц = 0.1

40 60 80

Время, отн. ед.

40 60 80

Время, отн. ед.

40 60 80

Время, отн. ед.

Рис. 6. Численные решения уравнений системы (5) при е = 50 и различных /и . а — х(т), б — у(т), в — г(т)

(/и = 0.1, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0; кривые расположены справа налево)

а

в

откуда с учетом начальных условий х1 (0) = 0, у1 (0) = 1, г1 ( 0) = 0 получается решение

х (т)= Я (еХР (-Ят/)- еХР (Я)),

у1 (т) = ЯЯ((Яр -1)ехР(Я) +

+ (1 - Ят) ехР (-Яpt)),

(8)

1 (т) = 1 - Я (ЯР ехР (-Я) - Ят ехР (Я)),

где введены обозначения Я = -^(1 + u + v)2 -4/ , Яр =(1 + /и + v + Я)/2, Ят =(1 + /и+v -Я)/2. Аналогичным образом уравнения для степени е2 имеют вид

с2 = -(u+v) х2 + у 2 +

2

+ ехР (-(Яр + Ят) t)

t )-

1 Яр ехР(-2Яmt)-ЯехР(-2Яр^,

Я2

у 2 = '^х2 - у 2 -

Я2

Я2 ехР (

2 - ЯР - ^"("(Яр + Ят )t)

1 - Я

+

1 - Я

Я2 -ехР(-2Я„/) + яехР(-2Яpt),

г2 — /их,2,

(9)

t)+

а их решение с учетом начальных условий х2 (0) = 0 , у2 (0) = 0, г2 (0) = 0 имеет вид

Х2 (т)= Хтт ехР (-2^) +

+ Xтр ехр (-(лт + Лр) t) + Хрр ехр (-2 Лpt) +

+ хт ехр (Л) + хр ехр (-Л), У2 (т)= Ущт ехр (-2Л„/) + + ?тР ехр (-(Л + Лр) t) + Урр ехр (-2Лpt) + (10) + Ут ехр (-Л) + Гр ехр (-Л),

Z2 (т) = Ятт ехр (-2Л) +

+ ^тр ехр (-(Лт + Лр) t) + 2рр ехр (-2Л^) ■ + Ят ехр (-Л) + Яр ехр (-Л),

+

где

Х_ =

2 (и-Лт )

Хрр =

X =

тр

х_ =

хр =-

У_ =

Л2 (л2 - Л(Лр + Лт) + и)'

2 (и-Лр) л2 (Л2 + л(Лр + Лт) + и), (1 + и+у)(1 - ¡и-у)

¡Л2 '

Лт (1 - 2и+у - Л(1 - 2 и - 2у)) ¡Л(2Л2- и) Лр (1 - 2и + у + Л(1 - 2и- 2у))

¡Л(2Л2- и)

2Лт + иЛр - 3и

Урр =

Л2(л2 -л(Лр + Лт) + и)'

2Лр + иЛт - 3и л2 (Л2 + л(Лр + Лт) + и),

У™ =-

У_ =

ур =

= -

Я- =

(1 +у)(1 - ¡-у)

¡Л2 '

(и-Лт )(1 - 2и + у - Л(1 - 2 и - 2у))

иЛ(2Л2- и) (Лр -и)(1 -2и+у + Л(1 -2и-2у)) иЛ(2Л2 - и)

и(Лр - и)

Л2 (л2 - л(Лр + Лт) + и)

и(и-Лт )

Л (Л2 + л(Лр + Лт) + и)' 1 - ¡-У

Л2

ур =

1 -2и+у - Л(1 -2и-2у) Л(2Л2 - и) 1 - 2и + у + Л(1 - 2и-2у) Л(2Л2- и) .

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

К сожалению, эти ряды не слишком хорошо сходятся, как показывает сравнение приближенных аналитических выражений и результатов численного решения системы дифференциальных уравнений: качественно поведение решений совпадает, но количественное совпадение оставляет желать лучшего. Во многих случаях разложение по обратным степеням параметра л, использованное в разделе "ПРЕДЕЛЬНЫЕ СЛУЧАИ", дает лучшие результаты. Является целесообразным использовать дробно-линейное представление для приближенного решения (см. результаты упомянутого раздела). В этом случае решение ищется в виде

х (т) =

У (т) =

[еХ1т + е2Х2т + -) + еЛ (еХ,г + е2Х2р + -) + +Лр > (е2Х2рт + -) + -1 + (еРт + е2Р2т + -) + е^ (еРрр + е2Р2г + -)

е Л Цт + е2У2т + -) + е^ (еУ1 р + е2У2 р + -) + е^+Лр > (е2У2 рт + -) + -

1 + (едш + е2^2т + -) + е Л (е^р + е2^2р + -)

z (т) = е-х (т)-У (т),

(11)

где Хр, Хт — изначально неизвестные характеристические показатели экспонент (вообще говоря, являющиеся функциями параметра е), а Х1т, Х2т, Y'm, Y2m, ... — неопределенные коэффициенты. Если ограничиваться квадратичными членами в числителе и линейными членами в знаменателе, то анализ показывает, что при такой параметризации решения образуются четыре дополнительных (свободных) параметра, с помощью которых можно управлять точностью решения.

ЗАКЛЮЧЕНИЕ

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

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

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

Работа выполнена при финансовой поддержке программ фундаментальных исследований Президиума РАН № 24 "Фундаментальные основы технологий наноструктур и наноматериалов", проект "Исследования кинетики и аналитических характеристик полимераз-ной цепной реакции в реальном времени ".

СПИСОК ЛИТЕРАТУРЫ

1. Эмануэль Н.М., Кнорре Д.Г. Курс химической кинетики. М.: Высш. шк., 1984. 463 с.

2. Корниш-Боуден Э. Основы ферментативной кинетики М.: Мир, 1979. 282 с.

3. Gevertz J.L., Dunn S.M., Roth C.M. Mathematical model of real-time PCR kinetics // Biotechnol Bioeng. 2005. Vol. 92. P. 346-355.

4. Booth C.S., Pienaar E., Termaat J.R., Whitney S.E., Louw T.M., Viljoen H.J. Efficiency of the polymerase chain reaction // Chem. Eng. Sci. 2010. Vol. 65. P. 4996-5006.

5. Cobbs G. Stepwise kinetic equilibrium models of quantitative polymerase chain reaction // BMC Bioinformat-ics. 2012. Vol. 13. P. 203-216.

Институт аналитического приборостроения РАН, г. Санкт-Петербург (Бердников А.С., Фёдоров А.А., Курочкин В.Е.)

ЗАО "Синтол", г. Москва (Сочивко Д.Г.)

Контакты: Фёдоров Алексей Александрович, [email protected]

Материал поступил в редакцию: 24.01.2014

ANALYSIS OF ELEMENTARY ENZYMATIC REACTION

EQUATION SYSTEM

A. S. Berdnikov1, A. A. Fedorov1, D. G. Sochivko2, V. E. Kurochkin1

1 Institute for Analytical Instrumentation of RAS, Saint-Petersburg, RF 2Joint Stock Company "Syntol", Moscow, RF

Mathematical equations for enzymatic reaction kinetics usually cannot be solved analytically and require numeric methods of analysis. However, successful modeling some problems critically depends on our ability to obtain at least approximate analytical solution. This applies, for instance, to the mathematical model of real-time

polymerase chain reaction which could be used to develop an efficient algorithm for experimental data processing. In the present paper, a novel approach to the analysis of enzymatic kinetics equation system is described which allows obtaining approximate analytic solutions expressing time course of reactant concentrations, with some limitations on reaction parameters.

Keywords: enzymatic kinetics, mathematical modeling

REFERENСES

1. Gevertz J.L., Dunn S.M., Roth C.M. Mathematical model of real-time PCR kinetics. Biotechnol Bioeng, 2005, vol. 92, pp. 346-355.

2. Booth C.S., Pienaar E., Termaat J.R., Whitney S.E., Louw T.M., Viljoen H.J. Efficiency of the polymerase

chain reaction. Chem. Eng. Sci., 2010, vol. 65, pp. 4996-5006.

3. Cobbs G. Stepwise kinetic equilibrium models of quantitative polymerase chain reaction. BMC Bioinforma-tics, 2012, vol. 13, pp. 203-216.

Contacts: Fedorov Alexey Alexandrovich, [email protected]

Article arrived in edition: 24.01.2014