Комп’ютерне моделювання та числова оптимізація метаболічних процесів: застосування до біосинтезу флавінових коферментів у дріжджів
А. Цимбрик, Д. Санагурський
Львівський національний університет імені Івана Франка
вул. Грушевського, 4, 79005 Львів, Україна
e-mail: biolog@franko.lviv.ua
Рибофлавін, вітамін В2, є попередником флавінмононуклеотиду та флавінаденіндинуклеотиду - кофакторів для широкого класу ферментів метаболізму, має комерційну цінність як домішка в харчовій промисловості, кілька видів флавіногенних мікроорганізмів, зокрема деякі види дріжджів, застосовуються в промисловості для виробництва рибофлавіну шляхом ферментації. З огляду на це актуальні біохімічні [3-5,8,9] та генетичні [10,13,19,25] дослідження флавіногенезу з метою отримання мутантів-надсинтетиків рибофлавіну для його біотехнологічного виробництва. Одним із можливих шляхів подолання обмежувальних факторів експериментального дослідження є аналіз властивостей відповідної моделі метаболізму із застосуванням комп’ютерного моделювання [14,22,31]. Математичний опис метаболічних систем дає змогу розраховувати очікувану відповідь метаболізму на генетичні модифікації та зміни зовнішнього середовища [18,26], ідентифікувати найважливіші об’єкти для метаболічної інженерії [11,12,15,29]. Математичну методологію успішно використовують для визначення перспективних напрямків побудови генетичних мутантів і пояснення неочікуваних експериментальних результатів [6,17,23,26,28,30,32]. Ми досліджували можливості застосування засобів комп’ютерного програмування та методів числової оптимізації для раціональної метаболічної інженерії (зокрема, побудови мутантів дріжджів P.guilliermondii – надсинтетиків рибофлавіну) шляхом аналізу стаціонарної поведінки моделі метаболічної системи біосинтезу флавінових коферментів.
Модельний експеримент проводили в середовищі MathCad 7. Pro, застосовували засоби Gepasi 3.1 [20,21] (програми MS Windows для моделювання стаціонарних станів і часової поведінки біохімічних реакцій). Кінетичні змінні та параметри, необхідні для розрахунків, брали з доступних баз даних (Інтернет) і з експериментальних робіт – досліджень регуляції надсинтезу рибофлавіну [3-5, 8-10, 13, 19, 25].
Математичний опис метаболічної системи. Фундаментальними рівняннями для опису метаболічних систем є рівняння Кірхгофа, відомі як рівняння масового балансу:
(1)
де Si – концентрації метаболітів; vj – швидкості реакцій, що каталізуються ферментами Ej; Ci,j - стехіометрична матриця, утворена малими цілими числами, переважно 1, -1, 0, і фактично містить більше нулів, оскільки значна частина ферментів перетворює лише кілька метаболітів; у системі, що складається з n метаболітів та m ферментів, Ci,j буде n x m матрицею.
Диференціальне рівняння (1), доповнюється кінетичними рівняннями ферментів:
(2)
У загальному випадку vj – це швидкість реакції на одиницю кількості j-го ферменту відносно до одного з його продуктів (Sn). Тоді Ci,j,vj відповідно буде швидкістю реакції j-го ферменту для i-го метаболіту.
Підставивши вираз (2) у рівняння (1), отримуємо:
i=1,2,…,n. (3)
Концентрації метаболітів ( Si) і потоки ( vj ) є системними змінними , які обчислюються за модельними рівняннями (1–3), кінетичні константи (pk) відіграють роль параметрів, вони не визначаються в моделі, а беруться з експериментально чи аналітично отриманих кінетичних даних ізольованих ферментів.
Параметр pk може набувати різного вигляду, наприклад, це може бути концентрація ферментів (Ej); максимальна активність (Vm); константи Міхаеліса-Ментен (Km); константи інгібування (Ki); концентрації зовнішніх ефекторів, активаторів (A) чи інгібіторів (I) тощо.
Рівняння (3) формують систему залежних від параметрів звичайних диференціальних рівнянь першого порядку, що описує поведінку метаболічної системи. Функція f є нелінійною завдяки перебігу бі- та тримолекулярних реакцій у метаболічному ланцюгу перетворень і нелінійному характеру кінетичних рівнянь ферментів. З цього випливає, що у більшості випадків розв’язок такої системи рівнянь не може бути отриманий аналітичним шляхом, а потребує застосування чисельних методів.
Стаціонарний стан системи. Він визначається умовами, за яких часова похідна системи змінних дорівнює нулю, тобто
f ( S; p ) = 0. (4)
Наслідком нелінійності системного рівняння (4) дуже часто є існування кількох мультистаціонарних станів, які можуть бути або стабільними, або нестабільними. Відповідно до визначення стаціонарності динамічних систем стаціонарний стан розглядається як стабільний, якщо після пертурбацій система повертається до вихідного стану ( або близько до нього). Cистема, яка повертається до стабільності при t®
, називається асимптоматично стабільною [6].
Динамічні зміни. Для обчислення динамічних змін необхідно розв'язати систему диференціальних рівнянь (3) для заданих значень параметрів. Динамічні зміни метаболічної системи додатково до впливу інгібіторів чи активаторів можуть бути зумовлені змінами параметрів системи, температури, рН середовища та ін.
Динамічні зміни in vitro дуже часто досягаються за умов, коли система не може досягти стаціонарності внаслідок деградації вихідного субстрату чи накопичення кінцевого продукту, який виступає як інгібітор реакції. Більше того, система може генерувати динамічні зміни автономно, якщо вона не перебуває в стабільному стаціонарному стані. В такому випадку можуть бути виявлені осциляції субстрату, продуктів реакції чи проміжних метаболітів [1,2,7]. Загалом часова залежність системних змінних визначається шляхом чисельного інтегрування (а не аналітичного дослідження) концентрацій метаболітів чи потоків як функцій від часу.
Кінетична модель. Сукупність реакцій, які беруть участь у метаболізмі флавінів у дріжджів, можна зобразити так:
Рис. Схематичне зображення метаболізму рибофлавіну у дріжджів P.guilliermondii. Назви метаболітів і повні назви ферментів наведені у додатку. Пунктирні лінії зображають регуляторні впливи, символи +/- вказують на відповідний ефект (активуючий чи інгібуючий). Швидкості реакцій, що є суттєвими для регуляторних взаємодій, описані нижче.
Пуринові попередники перетворюються у рибофлавін і згодом у флавінові коферменти ФМН і ФАД, піддаються складному регуляторному впливу як зовнішніх, так і внутрішніх ефекторів [20,21]. Структура метаболічної системи та кінетика її ферментів описані відповідно до експериментальних даних [19,22,23]. У моделі прийнято, що:
Фермент першої реакції синтезу рибофлавіну ГТФ-циклогідролаза (Е1) піддається неконкурентному інгібуванню з боку продукту метаболізму Р, тоді кінетичне рівняння швидкості реакції (vE1) матиме такий вигляд:
,
де V – максимальна швидкість реакції; S – концентрація субстрату реакції; I – концентрація інгібітору реакції; Ki – константа інгібування; Km – концентрація субстрату, при якій v=V/2 для I=0 (константа Міхаеліса);
1. Рівняння швидкості двосубстратних реакцій vE2, vE8, vE9, (механізм Ordered BI BI) матимуть такий вигляд:
, де Vf – максимальна швидкість прямої реакції; Vr – максимальна швидкість зворотної реакції; Keq - константа рівноваги; KmA – константа Міхаеліса для першого субстрату; KmВ - константа Міхаеліса для другого субстрату; KmP - константа Міхаеліса для першого продукту; KmQ - константа Міхаеліса для другого продукту; KiА - константа інгібування першим субстратом прямої реакції; KiВ - константа інгібування другим субстратом прямої реакції; KiP – константа інгібування зворотної реакції;
3. Рівняння швидкості реакцій, що каталізуються ферментами Е3–E6, згідно з кінетикою Генрі - Міхаеліс - Ментен, матимуть такий вигляд:
,
де V – максимальна швидкість реакції; Km – константа Міхаеліса; S – концентрації субстратів.
Математична модель. Запропонованій кінетичній моделі (див. рисунок) відповідає така система диференціальних рівнянь, складена відповідно до рівнянь масового балансу (1):
(5.1-5.5)
(5.6-5.11)
Аналізуючи зображений на рисунку метаболічний шлях, ми ставили за мету знайти умови, за яких потік у напрямку утворення продукту реакції Р є найвищим при мінімальних витратах субстратів (S1 і S2). З цією метою ми припустили, ніби технологічно впливаємо на концентрації ферментів (наприклад, клонуванням) і можемо маніпулювати їхніми кінетичними параметрами. Проблема полягає в тому, якими саме параметрами необхідно маніпулювати, щоб досягнути необхідного результату (в нашому випадку - збільшення виходу продукту).
Для цього провели комп’ютерний експеримент: математичний опис біосинтетичного шляху (5.1-5.11) ми задали в середовищі прикладної програми Gepasi 3.1, назначили довільно вибрані числові значення для різних параметрів і промоделювали стаціонарний стан заданої системи. Отримані результати прийняли як вихідні для подальших досліджень (наш “дикий тип”). Тоді вибрали 10 параметрів, які могли бути потенційними мішенями для генетичних маніпуляцій: всі концентрації ферментів, загальна концентрація кофакторів, плюс кінетичні параметри (також ті, які визначають дію Р як інгібітору Е1). Для цих параметрів визначили межі, в яких вони можуть змінювати свої значення в процесі оптимізації:
максимальні швидкості реакцій (V) між 0.001 та 1000;
загальна кількість пулу А ([A]+[AH]) у межах від 1Е-6 до 1Е+6;
регулюючі константи (для Р) у межах від 1Е-7 до 1Е+7.
Під час комп’ютерного експерименту алгоритми числової оптимізації реалізовувались у процесі моделювання стаціонарної поведінки системи, це дало змогу застосувати різні алгоритми оптимізації для досягнення поставленої мети. Для дослідження біосинтетичного шляху (див. рисунок) застосували п’ять оптимізацій-них методів [23,24,27].
На підставі аналізу результатів комп’ютерного експерименту, наведених у таблиці, виявили, що найліпший результат визначення максимуму функції (потоку J8), що на два порядки перевищує значення “дикого типу”, отримано із застосуванням методу випадкового пошуку (Random search). Випадковий пошук –
метод оптимізації, у процесі реалізації якого оптимум знаходиться шляхом перевірки значень досліджуваної функції в серії комбінацій випадкових значень
змінних параметрів. Випадкові значення генеруються у попередньо визначених межах, більше того, будь-які комбінації значень параметрів, що виходять за визначені межі, вилучаються з розрахунку.
Порівнюючи значення кінетичних параметрів “дикого типу” з даними, отриманими у процесі оптимізації, можна зробити висновок, що оптимальні значення потоку J8 одержали за умов, коли концентрації ферментів Е1,Е4, Е8 досягали своїх верхніх граничних значень (від 10 до 100), тоді як концентрації інших ферментів набували нижніх граничних значень (у межах від 0.01 до 0.1).
Розрахункові значення потоку (J8), отримані в результаті чисельного інтегрування з використанням різних оптимізаційних алгоритмів
Застосований метод | Значення потоку (J8) мкM/c, 10-4 | Кількість моделювань |
| |||
| “Дикий тип” | 0.028 | 1 | |||
| Steepest descent | 0.32 | 18 | |||
| L-BFGS-B | 0.30 | 9 | |||
| Hooke and Jeevs | 0.96 | 193 | |||
| Genetic algorithm | 0.78 | 903 | |||
| Random search | 5.55 | 1000 | |||
На підставі цих результатів ми можемо припустити, що для того, щоб досягнути максимуму виходу продукту (Р) за незмінних концентрацій субстратів, необхідно:
забезпечити надекспресію ферментів Е1, Е4, Е8;
максимально збільшити чутливість цих ферментів до їхніх субстратів (тобто, зменшити Km);
знизити інгібуючу дію продукту (Р) на першу реакцію метаболізму.
Додаток
Символ, Xi | Назва метаболіту | |
S1 S2 X1X2 X3 X4 X5 X6 X7 X8 P | Гуанозинтрифосфат (GTP) Рибулозо-5`-фосфат 2,5-диаміно-4-окси-6-рибозиламінопирімидин-5`-фосфат 2,5-диаміно-4-окси-6-рибітиламінопіримідин-5`-фосфат 5-аміно-2,4-диокси-6-рибітиламінопіримідин-5`-фосфат 5-аміно-2,4-диокси-6-рибітиламінопіримідин 3,4-дигідрокси-2-бутанон-4-фосфат 6,7-диметил-8-рибітиллюмазин Рибофлавін Флавінмононуклеотид (FMN) Флавінаденіндинуклеотид (FAD) | |
| Закінчення додатка | |
Символ, Ei | Назва ферменту | ЕС no. |
E1E2E3 E4 E5 E6 E7 E8 E9 | GTP-циклогідролаза Редуктаза Дезаміназа Фосфатаза 3,4-дигідрокси-2-бутанон-4-фосфат синтаза 6,7-диметил-8-рибітиллюмазинсинтаза Рибофлавінсинтаза Рибофлавінкіназа ATФ: ФМН-аденілилтрансфераза | 3.5.4.25 1.1.1.193 3.5.4.26 3.3.1.- ---- 2.5.1.9 2.5.1.9 2.7.1.26 2.7.7.2 |
_______________
1.Гольдштейн Б.Н., Горянин Н.Н. Моделирование колебаний двух активностей 6-фосфофрукто-2-киназы/фрукто-2,6-биофосфатазы в печени крысы // Молек. биол. 1996. Т.30. №4. С.933-940.
2.Иванова А.Н., Гольдштейн Б.Н. Кинетическая модель биофункционального фер-мента /фосфорилирование 6-фосфофрукто-2-киназы/фруктозо-2,6-биофосфатазы/ // Молек. биол. 1986. Т.20. №6. С.1522-1528.
3.Кащенко В.Е., Преображенская Е.Н., Сибирный А.А. Новый метод определения активности рибофлавинкиназы // Биохимия, 1994. Т.59. С.1254-1265.
4.Кащенко В.Е., Преображенская Е.Н., Шавловский Г.М. Получение и свойства иммобилизованной рибофлавинкиназы из дрожжей Pichia quilliermondii // Укр. биохим. журн. 1989. Т.61. №1. С.32-36.
5.Логвиненко Е.М. Шавловский Г.М., Конторовская Н.Ю. О биохимических функциях продуктов генов RIB5 и RIB6, участвующих в биосинтезе рибофлавина у дрожжей Pichia guilliermondii // Генетика. 1987. Т.23..№9. С.1699-1701.
6.Рубин А.Б. Биофизика: у 2 кн.: Учеб.для биол.спец.вузов. Кн.1 Теоретическая биофизика. М.,1997.
7.Сельков Е.Е., Шевелев Е.Л. Колебания и триггерные явления в обмене фруктозо-2,6-биофосфата. Математическая модель // Биофизика. 1987. Т.32. №2. С.242-247.
8.Шавловский Г.М., Логвиненко Е.М. Биогинез флавинов у дрожжей // Укр. биохим. журн. 1985. Т.57. №4. С. 104-112.
9.Шавловский Г.М., Логвиненко Е.М. Сверхсинтез флавинов у микроорганизмов и его молекулярные механизмы // Прикл. биох. и микр. 1988. Т.24. №4. С. 435-447.
10. Шавловский Г.М., Бабяк Л.Я., Сибирный А.А., Логвиненко Е.М. Генетический контроль биосинтеза рибофлавина у дрожжей Pichia quilliermondii // Генетика. 1985. Т.21.№3. С.368-374.
11. Cornish-Bowden A. Metabolic control analysis in biotechnology and medicine // Nature Biotechnology. 1999. Vol. 17. P.641-643.
12. Cornish-Bowden A., Eisenthal R. Prospects for Pharmacological Manipulation of
Metabolism // New Beer in an Old Bottle: Eduard Buchner. and the Growth of Biochemical Knowledge. University of Valencia. 1997. P.215-224.
13. Garcia-Ramirez J.J., Santos M.A., Revuelta. The Saccharomyces cerevisiae RIB4 Gene
Codes for 6,7-Dimethyl-8-ribityllumazine Synthase Involved in Riboflavin Biosynthesis. Molecular characterization of the gene purification of the encoded protein // J. Biol. Chem. 1995. Vol.270. №40. P.23801-23807.
14. Goryanin J., Selkov E. Mathematical simulation of cellular metabolism and regulation / Proc.of the IMACS. Symp. on Mathematical Modelling. Austria, 1994. P.332-336.
15. Hatzimanikatis V. Nonlinear Metabolic Control Analysis // Metabolic engineering. 1999. Vol. 1. P.75-87.
16. Hatzimanikatis V., Emmerling H., Sauer U., Baily J.E. Application of mathematical tools for metabolic design of microbial ethanol prodaction // Biotechnol. Bioeng. 1998. Vol.58. P.154-161.
17. Heinrich R. Mathematical models of metabolic systems: General principles and memb-rane transport in erythrocytes // Biomed. Bioch. Acta. 1985. Vol. 44. №6. P. 913-927.
18. Holzhütfer H., Jacobasch G., Biselorff A. Mathematical modeling of metabolic path-way affected by an enzyme deficiency // Eur. J. Bioch. 1985. Vol.149. №1. P.101-111.
19. Liauta-Teglivets O., Hasslacher M., Boretskii I. et.al. Mole-cular cloning of the GTP-cyclohydrolase structural Gene RIB1 of Pichia guilliermondii involved in riboflavin biosynthesis // Yeast. 1995. Vol. 11. №10. P.945-952.
20. Mendes P. Biochemistry by numbers: simulation of biochemical pathways with Gepasi 3 // Trends Biochem. Sci.1997. Vol.22. P.361-363.
21. Mendes P. GEPASI:A software package for modelling the dynamics, steady states and control of biochemical systems // Comput. Applic. Biosci. 1993. Vol.9. P.563-571.
22. Mendes P., Kell D. Computer Simulation of Biochemical Kinetics // In BioThermo-Kinetics of the living cell. Еds. H.V.Westerhoff, J.L.Snoep, F.E.Sluse, et al. BioTermoKinetics Press. Amsterdam, 1996. P.254-257.
23. Mendes P., Kell D.B. Numerical optimisation and simulation for rational metabolic engeneering // In BioThermoKinetics in the past genomic era Еds. C.Harsson, I.-L.Pahlman and H.Gustaffson, Chalmers Reproservice. Gateborg, 1998. P.348-349.
24. Mendes P., Kell D.B. Nonlinear optimization of biochemical pathways:applications to metabolic engineering and parameter estimation // Bioinform.1998. Vol.14. P.869-883.
25. Mortl S., Fischer M., Richter G., Biosynthesis of Riboflavin Lumazine synthase of Escherichia coli // J. Biol. Chem. 1996. Vol. 271. №52. P.33201-33207.
26. Nieberbeger P., Prasad R., Miozzori G., Kacser H. A strategy to increasing on in vivo flux by genetic manipulation the tryptophan system in yeast // Biochem. J. 1992. Vol.187. P.473-479.
27. Savinell J. M., Palsson B.O. Network analysis of intermediary metabolism using linear optimization. I. Development of mathematical formalism // J.Theor.Biol. 1992. Vol.154. P.421-454.
28. Schauer M., Heinrich R., Rapoport S.M. Mathematical modelling of glicolis and of the adenine nucleotide metabolism of human erythrocytes // Acta biol. Med. Germ. 1981. B.40. S.1659-1697.
29. Schuster S. Use and Limitations of Modular Metabolic Control Analisis in Medicine and Biotechnology // Metabolic Engineering. 1999. Vol. 1. P.232-242.
30. Varner J., Ramkrichna D. Metabolic ehgineering from a cybernetic perspectiv: Aspartate family of amino acids // Metab. Engineering. 1999. Vol.1. P.88-116.
31. Varner J., Ramkrishna D. Mathematical models of metabolic pathways // Current Opinion in Biotechnology. 1999. Vol.10. №2. P.146-150.
32. Vaseghi S., Baumeister A., Rizzi M., Reuss M. In vivo Dynamics of the Pentose Phosphate Pathway in S. cerevisiae // Metab. Engineering. 1999. Vol.1. P.128-140.