( PDF ) Análise computacional de redes metabólicas com regulação

Download PDF

ads:

UNIVERSIDADE FEDERAL DE SANTA CATARINA

PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA QUÍMICA

ITAMAR LEITE DE OLIVEIRA

ANÁLISE COMPUTACIONAL DE REDES

METABÓLICAS COM REGULAÇÃO

Florianópolis - SC

2008

ads:

Livros Grátis

http://www.livrosgratis.com.br

Milhares de livros grátis para download.

UNIVERSIDADE FEDERAL DE SANTA CATARINA

PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA QUÍMICA

ITAMAR LEITE DE OLIVEIRA

ANÁLISE COMPUTACIONAL DE REDES

METABÓLICAS COM REGULAÇÃO

Tese submetida à Universidade Federal de Santa Catarina como parte dos requisitos para

obtenção do título de Doutor em Engenharia Química

Prof. Luismar Marques Porto, Ph.D

Orientador

Florianópolis - SC

2008

ads:

ANÁLISE COMPUTACIONAL DE REDES

METABÓLICAS COM REGULAÇÃO

ITAMAR LEITE DE OLIVEIRA

Esta tese foi julgada adequada para a obtenção do Título de Doutor em Engenharia Química e

aprovada em sua forma final pelo Curso de Pós Graduação em Engenharia Química.

Prof. Leonel Teixeira Pinto, Dr. Eng.

Coordenador do Curso

Banca Examinadora:

Prof. Luismar Marques Porto, Ph.D. - Orientador

Universidade Federal de Santa Catarina – EQA/CTC

Presidente

Prof. José Carlos Merino Mombach, D.Sc.

Universidade Federal de Santa Maria – Depto. de Física

Membro Externo

Prof. Eugênio Simão, Dr. Eng.

Universidade do Oeste de Santa Catarina – Centro Tecnológico

Membro Externo

Profa. Sônia Elena Palomino Castro Bean, Dr. Eng.

Universidade Federal de Santa Catarina – MTM/CFM

Membro Externo

Profa. Cíntia Soares, Dr. Eng.

Universidade Federal de Santa Catarina – EQA/CTC

Membro Interno

AGRADECIMENTOS

Ao professor Luismar Marques Porto, pela orientação e incentivo durante o

desenvolvimento deste trabalho. Seu apoio me estimulou a crescer científica, ética,

profissional e pessoalmente.

A minha querida esposa Maristela Corrêa Meller, pelo incansável apoio, dedicação,

carinho e, principalmente, pela sua paciência durante todos estes anos.

Aos meus pais e familiares, pela dedicação e apoio incondicional desde os meus

primeiros momentos de vida.

Aos meus amigos e colegas do grupo de engenharia genômica que contribuíram para

a minha formação científica e pessoal. Um agradecimento especial a Derce O. S. Recouvreux,

Claudimir Antonio Carminatti e Júlia de Vasconcellos Castro.

Ao CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) pelo

apoio financeiro através de financiamento de projetos.

RESUMO

Este trabalho reúne ferramentas computacionais de biologia de sistemas num

aplicativo chamado GEnSys (Genomic Engineering System), desenvolvido usando o

MATLAB

, e que é formado por um conjunto de módulos interdependentes que permite a

análise e simulação de redes de reações bioquímicas. Os módulos implementados realizam as

seguintes tarefas: simulação dinâmica determinística, análise de controle metabólico, análise

de fluxo metabólico, análise de balanço de fluxo e determinação de modos elementares. Uma

característica importante do GEnSys é que não há a necessidade da entrada direta da matriz

estequiométrica que representa o modelo biológico, pois existe um módulo que lê um arquivo

contendo as reações da rede metabólica, escritas em uma sintaxe específica. Este

procedimento gera automaticamente a matriz estequiométrica da rede. Para tanto, foi criada

uma linguagem formal batizada de Linguagem de Reações Bioquímicas (LRB), flexível e

intuitiva. Por exemplo, os nomes dos reagentes e produtos de uma reação podem conter

caracteres especiais, tais como +, -, espaço e apóstrofo. Esses símbolos são comuns nos

nomes dos compostos bioquímicos. Da forma como foi projetado, o GEnSys pode ser

facilmente extendido para outros tipos de métodos de análise envolvendo a matriz

estequiométrica.

A partir de dados da literatura foram modeladas e simuladas diversas redes

metabólicas. Para a Clostridium acetobutylicum foram consideradas todas as reações usadas

para a produção de ácidos, solventes e hidrogênio molecular a partir de dois substratos

distintos: glicose e glicerol. A rede modelada captura importantes aspectos qualitativos, de

acordo com padrão de distribuição de fluxos (velocidades de reação) desta bactéria. Em um

dos casos analisados foi possível obter valores quantitativos próximos aos valores reportados

na literatura. Foram determinados 29 modos elementares de fluxo para a produção de

hidrogênio molecular pela C. acetobutylicum. Destes, há 13 modos distintos envolvidos na

produção de H

: seis a partir da glicose, seis a partir do glicerol e um relativo ao NADH, o

qual foi considerado externo.

Um estudo de caso em que a regulação da expressão gênica é importante foi

conduzido comparando-se as simulações estocástica e determinística da via de biossíntese do

triptofano da Escherichia coli. Todos os mecanismos de regulação conhecidos do operon trp

(transcrição, tradução e inibição enzimática) foram considerados nas simulações estocástica e

determinística. Em ambas, o comportamento da atividade da enzima antranilato sintase foi

avaliado e comparado aos valores experimentais para uma E. coli selvagem.

ABSTRACT

This work considers a set of computational systems biology tools, written as a

MATLAB

toolbox called GEnSys (Genomic Engineering System). GEnSys comprises several

interdependent modules that allow analysis and simulation of biochemical reaction networks.

The developed implemented modules perform the following tasks: deterministic dynamic

simulation, metabolic control analysis, metabolic flux analysis, flux balance analysis and

determination of elementary modes. An important GEnSys feature is that the stoichiometric

matrix representing the biological model does not have to be supplied directly. A module that

reads an external file describing the metabolic network, in a specified, proposed syntax (the

Biochemical Reaction Language – LRB), automatically produces the stoichiometric matrix.

LRB is a flexible but intuitive language, that also allows the use of common biochemical

symbols such as +, -, blank space and primes. Due to its design strategy, GEnSys may be

easily extended for other analytical methods involving the stoichiometric matrix.

Several metabolic networks were modeled from literature data and simulated using

the developed tools. For Clostridium acetobutylicum all important reactions involved in the

production of acids, solvents and molecular hydrogen, using two substrates, glucose and

glycerol, were considered. The modeled network captures essential qualitative behavior,

according to expected flux distribution from that bacterium. In one of the analysed cases,

good quantitative agreement with literature data was found. Twenty-nine flux elementary

modes for the molecular hydrogen production were found, thirteen of them involved in the H

biosynthesis: six from glucose, six from glycerol, and one regarding NADH production.

NADH was considered an external metabolite.

A comparison between stochastic and deterministic simulation was carried out using

a model where genetic regulation is important, the tryptophan biosynthesis by Escherichia

coli. All known regulatory mechanisms for the trp operon (transcription, translation and

enzymatic inhibition) were taken into account. In both cases, the catalytic behavior of the

anthranilate synthase enzyme was evaluated and compared to experimental values reported to

wild type E. coli.

LISTA DE FIGURAS

Figura 1: Representação da matriz A.

a é um elemento situado na linha i e coluna j.............3

Figura 2: Os seis possíveis casos do problema de quadrados mínimos (PQM) em relação ao

posto da matriz A e suas dimensões m e n...........................................................................7

Figura 3: Matriz estequiométrica total. ....................................................................................13

Figura 4: Matriz estequiométrica dos metabólitos internos. ....................................................13

Figura 5: a) Exemplo de rede metabólica, onde os b

representam os metabólitos externos; v

e v

são as reações internas e v

, v

e v

as reações de fronteira. ..................................14

Figura 6: Rede metabólica com 4 metabólitos e 6 reações ......................................................19

Figura 7: Matriz kernel de N correspondente a rede da Figura 6.............................................19

Figura 8: Vetores bases do espaço nulo da matriz estequiométrica representados como uma

distribuição de fluxos na rede............................................................................................20

Figura 9: Uma rede metabólica (A) e seus modos elementares (B) e vias extremas (C)

correspondentes .................................................................................................................21

Figura 10: Matriz estequiométrica vista como uma transformação linear. Os quatros

subespaços vetoriais associados a N também são mostrados............................................27

Figura 11: Coeficiente de elasticidade (ε) de uma reação v

em relação ao substrato x

e um

inibidor I. ...........................................................................................................................30

Figura 12: Ilustração das restrições aplicadas a uma rede metabólica.....................................35

Figura 13: Esquema das diferentes aplicações de análise estequiométrica..............................37

Figura 14: Hierarquia dos módulos que compõem o GEnSys. Os módulos de onde partem as

setas usam (herdam) as informações dos módulos aos quais estão ligados.......................

Figura 15: Esquema gráfico para a classificação dos metabólitos como internos ou externos.

Nos casos (a), (b) e (c) o metabólito A é considerado como externo................................

Figura 16: Via metabólica representando a primeira parte da glicólise da Saccharomyces

cerevisiae

...........................................................................................................................61

Figura 17: Concentrações dos metabólitos internos ao longo do tempo. Estes gráficos foram

produzidos pelo comando

subplotTrans............................................................................70

Figura 18: Gráficos de todas as concentrações dos metabólitos internos num única figura.

Produzido pelo comando plotTransFig. ............................................................................70

Figura 19: Valores dos fluxos de cada reação (enzima) da rede metabólica no tempo............71

Figura 20: Comparação da atividade enzimática obtida pelos modelos determinístico (linha

pontilhada) e estocástico (linha cheia)...............................................................................

Figura 21: Número de moléculas de triptofano ao longo do tempo relacionado ao experimento

de derepression ..................................................................................................................77

Figura 22: Os dois primeiros passos na via de biossíntese de penicilina pela P. chrysogenum.

As linhas tracejadas indicam inibição da reação pelo metabólito .....................................78

vii

Figura 23: Rede bioquímica utilizada pela C. acetobutylicum para a conversão de glicose e

glicerol em hidrogênio, ácidos (acetato e butirato) e solventes (butanol e acetona) em

condições anaeróbicas........................................................................................................82

Figura 24: Esquema gráfico dos modos elementares de produção de H

a partir do glicerol,

considerando a Figura 23...................................................................................................87

Figura 25: Restrição de balanço de fluxo para a rede da C. acetobutylicum mostrada da Figura

23 .......................................................................................................................................93

Figura 26: Distribuição de fluxo da ABF ao se maximizar a biomassa com velocidade de

consumo de glicose 8,72 mmol/gdw

-1

...........................................................................94

Figura 27: Distribuição de fluxo da ABF ao se maximizar o H

com velocidade máxima de

consumo de glicose 8,7 mmol/gdw

-1

.............................................................................95

Figura 28: Cada grupo de barras representa um fluxo..............................................................97

Figura 29: Distribuição de fluxo obtida pela ABF ao se maximizar a biomassa ou a

biomassa+H2 ou a biomassa+butirato+acetato ou o butirato+acetato. .............................98

Figura 30: Dados experimentais versus simulados, ambos na fase solventogênica.................99

Figura 31: Dados experimentais versus simulados, ambos na fase acidogênica....................100

LISTA DE QUADROS

Quadro 1: Algoritmo do método direto....................................................................................17

Quadro 2: Cálculo da matriz T, em linguagem MATLAB

.......................................................45

Quadro 3: Exemplo de equação cinética. .................................................................................46

Quadro 4: Função MATLAB

gerada automaticamente para a equação cinética do Quadro 3. 46

Quadro 5: Exemplos de reações escritas em LRB....................................................................49

Quadro 6: Sintaxe de uma expressão cinética..........................................................................50

Quadro 7: Exemplo de reações da seção 2.5............................................................................51

Quadro 8: Exemplos de equações cinéticas escritas numa sintaxe própria para as equações da

seção 2.5. Os valores dos parâmetros são aleatórios. ........................................................51

Quadro 9: Sintaxe para atribuir valores iniciais aos metabólitos internos. ..............................52

Quadro 10: Exemplo de atribuição de valores iniciais.............................................................52

Quadro 11: Exemplo de reações e outras informações que formam um modelo gravado num

arquivo...............................................................................................................................56

Quadro 12: Exemplo de reações com metabólitos forçados a serem internos ou externos......59

Quadro 13: Reações em LRB da rede mostrada na Figura 16..................................................62

Quadro 14: Sintaxe para a criação de um modelo para a AFM................................................80

Quadro 15: Sintaxe para a construção de um modelo para ABF..............................................91

LISTA DE TABELAS

Tabela 1: Exemplo das informações armazenadas no campo sInf da estrutura modelinf.........57

Tabela 2: Cinéticas pré definidas no GEnSys ..........................................................................67

Tabela 3: Constantes usadas nas simulações determinística e estocástica...............................75

Tabela 4: Valores iniciais para os compostos na simulação estocástica. Estes valores refletem

um experimento de derepression.......................................................................................77

Tabela 5: Coeficiente de controle de fluxo para as quatro enzimas da rede mostrada na Figura

22. ......................................................................................................................................79

Tabela 6: Resultado dos fluxos calculados usando os dados da Tabela 11 como fluxos

medidos..............................................................................................................................81

Tabela 7: Nomes das enzimas (ou fluxos) da rede metabólica da Figura 23. Algumas reações,

sequência linear, foram agrupadas formando uma única reação.......................................

Tabela 8: Significados dos nomes dos metabólitos da rede da Figura 23................................84

Tabela 9: Os 29 modos elementares correspondentes à rede da Figura 23..............................86

Tabela 10: Reações globais relacionadas aos 29 modos elementares. O sufixo _x em alguns

metabólitos significa que eles foram tratados como metabólitos externos........................88

Tabela 11: Dados experimentais relacionados à fase solventogênica......................................96

Tabela 12: Dados experimentais relacionados à fase acidogênica...........................................99

LISTA DE SÍMBOLOS E ABREVIATURAS

Abreviatura Significado

ABF Análise de balanço de fluxos

AcAcetil-CoA Aceto acetil coenzima A

AcetatoExt Acetato externo

ACM Análise de controle metabólico

ACVS ACV sintetase

ADP Adenosina difosfato

AFM Análise de fluxo metabólico

ATP Adenosina trifosfato

bisACV bis-(L-aminoadipil)-L-cysteinil-D-valina

ButiratoExt Butirato externo

CCC Coeficiente de controle de concentração

CCF Coeficiente de controle de fluxo

CoA Coenzima A

EDO Equação diferencial ordinária

Elmo Modos elementares

FdOxi Ferredoxina oxidada

FdRed Ferredoxina reduzida

GDH-3-P Gliceraldeído-3-fosfato

GEnSys Genomic Engineering System

Hidrogênio molecular

i,j Índices das linhas e colunas de uma matriz bidimensional

Matriz identidade de ordem n

IPNS Isopenicilina N sintetase

LactatoExt Lactato externo

LD Linearmente dependente

LI Linearmente independente

LLD-ACV (L-aminoadipil)-L-cisteinil-D-valina

LRB Linguagem das reações bioquímicas

ME Modo elementar

ModelInfo Módulo do GEnSys que contém os dados básicos do modelo

Abreviatura Significado

Matriz estequiométrica dos metabólitos internos

NAD+ Nicotinamida adenina dinucleotídeo

NADH Nicotinamida adenina dinucleotídeo, forma reduzida

PIR Piruvato

PQM Problema de quadrados mínimos

QR Decomposição ortogonal de uma matriz A em outras duas matrizes: Q e R

Posto ou rank de uma matriz

Matriz estequiométrica total

SVD Decomposição por valores singulares (singular value decomposition)

Vetor de fluxos de uma rede metabólica

v1 Velocidade da glicólise

v10 Acetaldeído desidrogenase

v11 CoA transferase 1

v12 CoA transferase 2

v13 Butiril-CoA desidrogenase

v14 Fosfotransbutirilase

v15 Butiraldeído desidrogenase

v16 Velocidade de crescimento celular

v17 ATPase

v18 Transporte de lactato

v19 Transporte de acetato

v2 Glicerol desidrogenase

v20 Transporte de butirato

v3 Velocidade de produção de piruvato

v4 Piruvato-ferredoxina oxidoredutase

v5 NADH-ferredoxina oxidoredutase;

v6 Hidrogenase

v7 Velocidade de produção de lactato

v8 Tiolase

v9 Fosfotransacetilase

VE Vias extremas

xii

Símbolo Significado

ℜ

Espaço vetorial n-dimensional

nm×

ℜ

Espaço vetorial das matrizes nm

<u,v> Produto interno entre os vetores u e v

Norma 2 do vetor u

∏

Produtório

Matriz dos coeficientes de elasticidade de uma rede metabólica

As matrizes e vetores foram escritos com letras maiúsculas e em negrito.

SUMÁRIO

1 INTRODUÇÃO...............................................................................................1

1.1 OBJETIVOS.....................................................................................................................1

1.1.1 Objetivo Geral............................................................................................................1

1.1.2 Objetivos Específicos ................................................................................................2

1.2 ORGANIZAÇÃO DO TEXTO........................................................................................2

2 FUNDAMENTAÇÃO TEÓRICA ..................................................................3

2.1 TÓPICOS DE ÁLGEBRA LINEAR................................................................................3

2.1.1 Matrizes .....................................................................................................................3

2.1.2 Espaços Vetoriais.......................................................................................................4

2.1.3 Subespaços de um Espaço Vetorial ...........................................................................4

2.1.4 Combinação Linear, Conjuntos Geradores, Base e Dimensão..................................4

2.1.5 Ortogonalidade...........................................................................................................5

2.1.6 Decomposição Ortogonal ..........................................................................................6

2.1.7 Problemas de Quadrados Mínimos............................................................................6

2.2 MATRIZ ESTEQUIOMÉTRICA ..................................................................................12

2.3 SIMULAÇÃO DINÂMICA...........................................................................................14

2.3.1 Simulação Determinística........................................................................................14

2.3.2 Simulação Estocástica..............................................................................................15

2.4 INFORMAÇÕES ESTRUTURAIS DA MATRIZ ESTEQUIOMÉTRICA..................18

2.4.1 Espaço Nulo da Matriz Estequiométrica .................................................................18

2.4.2 Modos Elementares de Fluxos e Vias Extremas......................................................20

2.4.3 Relações de Conservação.........................................................................................22

2.5 CINÉTICA DAS REAÇÕES BIOQUÍMICAS..............................................................25

2.6 SUBESPAÇOS RELACIONADOS À MATRIZ ESTEQUIOMÉTRICA....................26

2.6.1 Conteúdo dos Subespaços da Matriz Estequiométrica ............................................28

2.7 ANÁLISE DE CONTROLE METABÓLICO...............................................................28

2.7.1 Coeficiente de Elasticidade......................................................................................29

2.7.2 Coeficientes de Controle..........................................................................................30

2.7.3 Relações Estruturais.................................................................................................31

2.8 ANÁLISE DE FLUXO METABÓLICO.......................................................................33

2.9 ANÁLISE DE BALANÇO DE FLUXO........................................................................34

2.10 VISÃO GERAL DA MODELAGEM ESTEQUIOMÉTRICA .....................................36

2.11 PRODUÇÃO DE HIDROGÊNIO BACTERIANO.......................................................37

2.11.1 Simulação in silico da Produção de Biohidrogênio..............................................38

2.12 SIMULAÇÃO DINÂMICA DO OPERON TRP NA ESCHERICHIA COLI................39

xiv

3 METODOLOGIA..........................................................................................42

3.1 MATLAB

.....................................................................................................................42

3.2 IMPLEMENTAÇÕES EM MATLAB

.........................................................................43

3.2.1 Problemas de Quadrados Mínimos..........................................................................43

3.2.2 Linguagens das Reações Bioquímicas.....................................................................43

3.2.3 Análise de Fluxo Metabólico...................................................................................43

3.2.4 Análise de Balanço de Fluxo ...................................................................................44

3.2.5 Modos Elementares..................................................................................................44

3.2.6 Relação de Conservação..........................................................................................44

3.2.7 Geração da Equação Cinética ..................................................................................45

3.2.8 Simulação Determinística e Análise de Controle Metabólico.................................47

4 RESULTADOS E DISCUSSÃO ..................................................................48

4.1 LINGUAGENS DAS REAÇÕES BIOQUÍMICAS ......................................................48

4.2 EXPRESSÕES CINÉTICAS ENZIMÁTICAS E VALORES INICIAIS......................50

4.3 HIERARQUIA DOS MÓDULOS .................................................................................52

4.3.1 Módulo ModelInfo...................................................................................................54

4.3.2 Pré-processamento da LRB .....................................................................................55

4.3.3 Módulo totalStMat...................................................................................................57

4.3.4 Metabólitos Internos e Externos..............................................................................59

4.3.5 Módulo stMat...........................................................................................................60

4.4 SIMULAÇÃO DINÂMICA COM O GENSYS ............................................................61

4.4.1 Descrição do Exemplo para Simulação Dinâmica...................................................61

4.4.2 Matriz Estequiométrica Total ..................................................................................62

4.4.3 Matriz Estequiométrica Interna................................................................................63

4.4.4 Conjunto de Equações Diferenciais Ordinárias.......................................................64

4.4.5 Cinéticas das Reações..............................................................................................64

4.4.6 Valores Iniciais dos Metabólitos..............................................................................65

4.4.7 Forma Simplificada..................................................................................................66

4.4.8 Simulação Dinâmica do Modelo..............................................................................69

4.5 SIMULAÇÃO ESTOCÁSTICA DO OPERON TRP DA E. COLI ...............................72

4.5.1 Parâmetros ...............................................................................................................75

4.5.2 Resultados................................................................................................................76

4.6 ANÁLISE DE CONTROLE METABÓLICO NA PRODUÇÃO DE PENICILINA....78

4.7 ANÁLISE DE FLUXO METABÓLICO.......................................................................80

4.8 MODELAGEM E SIMULAÇÃO DA PRODUÇÃO DE BIOHIDROGÊNIO.............81

4.8.1 Rede Metabólica da C. acetobutylicum ...................................................................82

4.8.2 Modos Elementares da Rede Metabólica para a Produção de Biohidrogênio.........84

4.8.3 ABF da Produção de Biohidrogênio da C. acetobutylicum .....................................90

5 CONSIDERAÇÕES FINAIS......................................................................101

5.1 CONCLUSÕES............................................................................................................101

5.1.1 GEnSys ..................................................................................................................101

5.1.2 Simulação Estocástica da E. coli ...........................................................................103

5.1.3 Modelagem e Simulação da Rede Metabólica da C. acetobutylicum....................103

5.2 RECOMENDAÇÕES PARA TRABALHOS FUTUROS...........................................104

REFERÊNCIAS................................................................................................106

1 INTRODUÇÃO

Parte dos assuntos abordados neste trabalho situa-se em uma área emergente de

pesquisa conhecida como biologia de sistemas. Os métodos e técnicas usados em biologia de

sistemas requerem algum conhecimento em diversas outras áreas, tais como: biologia

molecular, engenharia, bioquímica, matemática e computação, por exemplo. Isto poderá ser

observado ao longo deste trabalho.

Serão discutidos também alguns métodos estudados em engenharia metabólica,

como análise de fluxo metabólico e análise de controle metabólico. Apesar de haver um forte

componente experimental em ambas as áreas de biologia de sistemas e engenharia metabólica,

serão vistos apenas alguns métodos matemáticos e computacionais para a análise e simulação

de redes metabólicas e de regulação.

Desde o sequenciamento do primeiro genoma completo, em 1995 – depois

impulsionado pelo sequenciamento do genoma humano –, a quantidade de dados públicos

disponíveis cresce rapidamente possibilitando a construção de redes metabólicas em escala

genômica. A obtenção destes dados é possível devido às tecnologias de alto rendimento

aplicado às ciências ômicas – genômica, proteômica, transcriptômica, entre outras.

O papel fundamental de métodos formais computacionais e matemáticos é extrair

informações, e relações entre elas, a partir desta grande quantidade de dados, de forma a

aumentar a compreensão do funcionamento da fisiologia dos organismos vivos. Estes

conhecimentos podem ser aplicados para o desenvolvimento de processos e produtos

biotecnológicos e, principalmente, na saúde humana. Neste contexto, foi desenvolvido o

GEnSys, um programa implementado em MATLAB

que reúne os principais métodos

computacionais para a análise e simulação de redes biológicas.

1.1 OBJETIVOS

1.1.1 Objetivo Geral

Sistematizar os principais métodos matemáticos e computacionais usados na biologia

de sistemas em um único aplicativo formando um conjunto de rotinas computacionais para a

modelagem e simulação de sistemas biológicos.

1.1.2 Objetivos Específicos

Para atingir o objetivo geral proposto, apresentam-se os seguintes objetivos

específicos:

 desenvolver um conjunto de módulos computacionais flexíveis de forma a proporcionar

uma fácil extensão destes módulos para outros tipos de análises que envolvam a matriz

estequiométrica;

 criar um conjunto de rotinas em MATLAB

para simulação dinâmica e no regime

estacionário de vias bioquímicas;

 fazer a análise de balanço de fluxo da rede metabólica para a produção de hidrogênio

bacteriano;

 realizar a análise estrutural, através dos modos elementares de fluxos, da rede metabólica

para a produção de hidrogênio bacteriano;

 desenvolver um modelo da via reguladora do triptofano na Escherichia coli para

simulação estocástica;

 criar uma linguagem para escrever as reações bioquímicas de forma a ser analisada e

processada pelos métodos computacionais.

1.2 ORGANIZAÇÃO DO TEXTO

Este trabalho foi organizado em 5 capítulos, além das referências que se encontram

no último capítulo. No capítulo 2 são abordados os principais conceitos e métodos usados

neste trabalho para as análises de redes metabólicas e reguladoras. No capítulo 3 são descritas

todas as abordagens usadas na implementação do GEnSys para simulação de vias metabólicas

e reguladoras. São descritos, principalmente, os recursos do MATLAB

que possibilitaram a

implementação das técnicas descritas no capítulo anterior. No capítulo 4 apresentam-se alguns

resultados e discussão de estudos de casos usando o GEnSys, bem como o próprio GEnSys.

Um destes estudos é para a produção de hidrogênio através da rede metabólica da C.

acetobutylicum

usando a análise de balanço de fluxo e modos elementares. No capítulo 5 são

apresentadas as conclusões e as recomendações para trabalhos futuros.

2 FUNDAMENTAÇÃO TEÓRICA

Neste capítulo serão abordados os principais elementos e métodos usados, neste

trabalho, para as análises de redes metabólicas e reguladoras.

2.1 TÓPICOS DE ÁLGEBRA LINEAR

Nesta seção será feita uma revisão de alguns tópicos de álgebra linear, uma vez que

muitos dos métodos usados fazem uso intensivo dos conceitos de álgebra linear.

2.1.1 Matrizes

Chama-se de matriz uma tabela de elementos dispostos em linhas e colunas. Dados

dois números m e n naturais, não nulos, chama-se matriz m por n toda tabela constituída por

nm × elementos dispostos em m linhas e n colunas. A Figura 1 apresenta esquematicamente

uma matriz

nm×

A cujos elementos são denotados por

a .

Figura 1: Representação da matriz A.

a é um elemento situado na linha i e coluna j.

Posto de uma Matriz

O posto linha de uma matriz A é o número de linhas linearmente independentes. O

posto linha de A é igual ao posto coluna de A e, como são equivalentes, trabalha-se apenas

com o posto (ou

rank) de A.

Matriz de Permutação

Uma matriz de permutação (P) é a matriz identidade (I) com suas linhas

reordenadas. Se P é uma matriz de permutação, então pré-multiplicar A por uma matriz de

permutação equivale a permutar as linhas da matriz A. Pós-multiplicar a matriz A por uma

matriz de permutação equivale a permutar as colunas de A. As matrizes de permutação são

ortogonais, então

PP =

−

⎥

⎦

⎤

⎢

⎣

⎡

mnmm

aaa

MOMM

22221

11211

m linhas

n colunas

nm×

2.1.2 Espaços Vetoriais

Um espaço vetorial real é um conjunto não vazio V munido de duas operações,

adição e multiplicação por escalar, tais que:

 u∀ e V∈v , então V∈+ vu . A operação de adição é fechada em V e além disso:

V∈∀+=+ vuuvvu ,,

;

()()

∈

∀++

++ wvuwuvwvu ,,, ;

- Existe em V um único vetor denominado vetor nulo e representado por 0, tal que

V∈∀+=+ uu00u ,

;

- Para cada V

∈

u , existe um único vetor denominado vetor oposto e representado por

u− , tal que

()

0=−+ uu .

 Se V∈v e k é um escalar então Vk

∈

v . A operação de multiplicação é fechada em V ,

isto é:

()

vuvu kkk +=+ , onde k é um escalar e

∈

vu,

;

()

uuu baba +=+

, onde

e b são escalares e V

∈

u ;

- para quaisquer escalares

e qualquer vetor

∈

(

)()

uu baab = ;

- para o escalar 1, uu =1 para qualquer vetor V

∈

u .

2.1.3 Subespaços de um Espaço Vetorial

Um subespaço W de um espaço vetorial V é um subconjunto não vazio contido em V

que satisfaz as seguintes propriedades:

1. se

W∈vu, então

W∈+ vu

;

seja a um escalar, e

W∈u

, então

∈

O item (1) acima declara que W é fechado em relação à operação de adição de

vetores. O item (2) estabelece que o conjunto W é fechado em relação à operação de

multiplicação. Como

VW ⊂ , W compartilha as propriedades de V. Todo subespaço W de V

contém o vetor nulo.

2.1.4

Combinação Linear, Conjuntos Geradores, Base e Dimensão

Seja

um espaço vetorial sobre

ℜ

e V

∈

m21

vvv ,,, K . Qualquer vetor em

forma

m21

vvv

aaa +

+ K

, onde

ℜ

∈

a , é chamado combinação linear de

m21

vvv ,,, K . O conjunto de todas as combinações lineares, denotado por

m21

vvv ,,, K ou

()

m21

vvv ,,, Kger , é chamado espaço gerado por

m21

vvv ,,, K . De modo geral, para

qualquer subconjunto

()

Sger

S consiste de todas as combinações lineares de vetores

de S.

Seja

V um espaço vetorial sobre

ℜ

. Os vetores

V∈

m21

vvv ,,, K

são ditos

linearmente dependentes (LD) sobre

ℜ

se existirem escalares

ℜ

∈

aaa ,,,

K não

simultaneamente nulos tais que 0

m21

vvv

aaa K . Caso contrário, são ditos

linearmente independente (LI) ou, simplesmente, independentes.

Um conjunto

{}

n21

uuu ,,, K=S de vetores é uma base de V se são verdadeiras as

condições:

n21

uuu ,,, K são LI e

n21

uuu ,,, K geram V . Em outras palavras, o conjunto S é

uma base de

V se todo vetor V

∈

v pode escrever-se de maneira única como combinação

linear dos vetores da base

S .

2.1.5

Ortogonalidade

Vetores Ortogonais

Dados dois vetores colunas

x e y do

ℜ quaisquer. Estes dois vetores são ditos

serem ortogonais quando

0, =yx ou 0=yx

, onde yx, é o produto interno dos vetores x

Conjuntos Ortogonais

Seja

{}

n21

uuu ,,, K=B um conjunto de vetores do

ℜ .

é um conjunto ortogonal

0,, ==≠∀

jiji

uuuu

ji . Se em particular i∀= ,1

u , então

é dito ser ortonormal,

ou seja:

⎩

⎨

⎧

≠

se 0

se 1

uu .

Todo conjunto ortogonal é LI e um conjunto de

n vetores ortonormais de um espaço

V n-dimensional é uma base ortonormal para V.

Matriz Ortogonal

Uma matriz

qualquer é ortogonal se

IQQ = e

IQQ = . Para uma matriz

ortogonal, as seguintes afirmações são equivalentes:

 Q possui colunas ortonormais;

 Q possui linhas ortonormais;



QQ =

−1

;



xQx = para todo

x ℜ∈ .

2.1.6 Decomposição Ortogonal

Dada a matriz A nm × com posto k; A pode ser decomposta da seguinte forma:

A = , onde: H é uma matriz ortogonal mm

; R é uma matriz nm × triangular

superior;

K é uma matriz ortogonal nn

. A decomposição de A, como visto acima, é

chamada de decomposição ortogonal.

Decomposição Ortogonal QR

Toda matriz

pode ser fatorada como QRA

, onde

mm×

ℜ∈Q

é ortogonal e

ℜ∈ é uma matriz triangular superior.

Decomposição Ortogonal por Valores Singulares (SVD

)

Qualquer matriz

pode ser escrita como

USVA = , onde U

m×m

e V

n×n

são

matrizes ortogonais e

m×n

é uma matriz diagonal com

(

)

(

)

Sdiag

,,,

, r é o posto

da matriz

A. Os números não negativos

são denominados de valores singulares de A.

Usualmente, adota-se a convenção

≥≥≥ K

. Se a matriz A possui posto r, então A

possui exatamente

r valores singulares estritamente positivos.

2.1.7 Problemas de Quadrados Mínimos

O problema de se achar uma solução, aproximada ou exata, para o sistema bAx

conhecido como o problema de quadrados mínimos (PQM), e é assim definido: dados uma

matriz

de dimensão nm × , com posto

(

)

nmk ,min

≤

, e um vetor m-dimensional b ,

encontrar um vetor

x que minimiza bAx − . Isto é equivalente a achar

x tal que bAx

≅

Dependendo do posto de

, de m e n, o sistema

bAx

pode ter solução exata

ou aproximada. A

Figura 2 exibe os seis possíveis casos para este problema.

SVD: singular value decomposition, em inglês.

Figura 2: Os seis possíveis casos do problema de quadrados mínimos (PQM) em relação ao

posto da matriz A e suas dimensões m e n. Nos casos 1a e 3a, é garantido que o sistema tem

solução exata.

A seguir, será apresentado um método geral para a solução do PQM

bAx =

descrito

por Lawson e Hanson (1995). O uso de transformação ortogonal permite expressar a solução

do PQM da seguinte forma:

Suponha que

A é uma matriz nm

com posto k e que

A = é uma

decomposição ortogonal de

A, onde:

 H é uma matriz ortogonal mm × ;

 R é uma matriz nm × da forma

⎥

⎦

⎤

⎢

⎣

⎡

(2.1)

 R

é uma matriz

kk ×

de posto k;

 K é uma matriz ortogonal nn × ;

m = n

posto(A) = k = m = n

Caso 1a

m = n

posto(A) = k < m = n

≅

Caso 1b

posto(A) = k = n < m

m > n

≅

Caso 2a

posto(A) = k < n < m

m > n

≅

Caso 2b

m < n

posto(A) = k = m < n

Caso 3a

m < n

posto(A) = k < n < m

≅

Caso 3b

Definindo o vetor

}

−

⎥

⎦

⎤

⎢

⎣

⎡

gbH

(2.2)

e introduzindo uma nova variável

}

−

⎥

⎦

⎤

⎢

⎣

⎡

yxK

(2.3)

e fazendo

y a solução única do sistema

1111

gyR =

(2.4)

Então:

 todas as soluções que minimizam bAx − são dadas por

⎥

⎦

⎤

⎢

⎣

⎡

onde y

é arbritário;

(2.5)

 o resíduo r da solução x

⎥

⎦

⎤

⎢

⎣

⎡

=−=

HxAbr

;

(2.6)

 a norma de r satisfaz

gAxbr =−= ;

(2.7)

 a solução de menor distância (da origem) é obtida por

⎥

⎦

⎤

⎢

⎣

⎡

(2.8)

Por questões de desempenho, o método acima será adaptado para os casos descritos

Figura 2. Em cada um dos casos é considerada a decomposição ortogonal QR. Tais

algoritmos foram adaptados de Lawson e Hanson (1995).

Sistema Sobredeterminado e Sistema Exato com Posto Completo

Nos casos 1a e 2a da Figura 2, a matriz nn

triangular superior R

é não singular.

Por isso a solução

para o PQM pode ser obtida computando:

}

−

⎥

⎦

⎤

⎢

⎣

⎡

gQb

(2.9)

e resolvendo:

111

gxR =

(2.10)

para o vetor solução

. O vetor resíduo é xAbr

−

, cuja norma é dada por:

gr ==

(2.11)

Sistema Indeterminado com Posto Completo

Considere o problema PQM,

bAx

, para o caso 3a da Figura 2, m < n e posto(A) =

m. O método de solução é descrito pelos seguintes passos:

[]

⎥

⎦

⎤

⎢

⎣

⎡

QQyx

ybRy

0:RAQ

arbritário com

(2.12)

A matriz Q é ortogonal e R é triangular superior e não singular. O vetor y

, m-

dimensional, é único. Para qualquer vetor

mn−

ℜ∈

y , o vetor x satisfaz bAx = . A solução

com comprimento mínimo é conseguida fazendo-se

Sistema com Posto Deficiente

O algoritmo a seguir é aplicado nos problemas onde A tem posto deficiente (A é

singular), que correspondem aos casos 1b, 2b e 3b da

Figura 2:

{

}

−

⎥

⎦

⎤

⎢

⎣

⎡

≡=

−

RQAP

}

cQb

−

⎥

⎦

⎤

⎢

⎣

⎡

(2.13)

[][]

0:WKR:R

1211

arbritário e

211

ycWy =

PKy

PKx

≡

⎥

⎦

⎤

⎢

⎣

⎡

(

)

0ycyRcAxb

222222

==−=−= se

O algoritmo determina a matriz ortogonal Q e a matriz de permutação P de modo

que R seja triangular superior e R

não singular.

PMQ Usando SVD

Para resolver o PQM

bxA

≅

,com nm ≥ , a decomposição SVD

()

nnm

mmnm

×−

××

⎥

⎦

⎤

⎢

⎣

⎡

(2.14)

pode ser usada substituindo-se o problema

bAx

≅

pelo problema equivalente

()

}

Vpx

bUb

⎥

⎦

⎤

⎢

⎣

⎡

−

⎥

⎦

⎤

⎢

⎣

⎡

≡

≅

⎥

⎦

⎤

⎢

⎣

⎡

(2.15)

A norma residual associada à solução

x é definida como

bAx −=

(2.16)

Pseudo-inversa

Seja

A = uma decomposição ortogonal de A, então

= é a

pseudo-inversa (PENROSE, 1955) de

A e R

é dada por

⎥

⎦

⎤

⎢

⎣

⎡

−

(2.17)

Esta representação explícita da pseudo-inversa é útil para computação da pseudo-

inversa. Geralmente, não é necessário calcular explicitamente a pseudo-inversa de

A. Por

exemplo, para o PMQ, é mais econômico – em termos de tempo e memória – computar a sua

solução usando os três passos a seguir:

ygyR

bHg

11111

⎥

⎦

⎤

⎢

⎣

⎡

para Resolver

(2.18)

PMQ com Restrição de Igualdade (PMQI)

Dados uma matriz

com posto k, um vetor coluna d de dimensão m

, uma

matriz

e um vetor f de dimensão m

. Achar o vetor x que minimiza

fEx − , sujeito a

dCx

(2.19)

Em outras palavras, o problema acima é achar o vetor coluna

x de tamanho m

que

minimiza a norma dada na primeira expressão sujeito a restrição de igualdade da segunda

equação.

PMQ com Restrição de Desigualdade

Dados:

 E uma matriz nm ×

;

 f um vetor coluna de dimensão m

;

 G uma matriz

nm ×

;

 h um vetor de tamanho m.

Os PMQs com restrições lineares de desigualdade podem ser enunciados como

segue:

 PMQ com restrição de desigualdade (PMQD)

Minimizar

fEx − , sujeito a

hGx ≥

;

 PMQ com restrição de não negatividade (NNLS, nonnegative least square)

Minimizar

fEx − , sujeito a 0x ≥ ;

 Problema LDP (least distance programming)

Minimizar

x , sujeito a hGx ≥ ;

Os problemas PMQI, PMQD e NNLS vêm implementados no M

ATLAB

. A seguir

será exibido o algoritmo para o problema LDP:

 Seja a matriz

(

)

m1n

×+

e o vetor coluna f de tamanho n+1

⎥

⎦

⎤

⎢

⎣

⎡

]1,0,0,0[

48476

K=f . Resolver o seguinte problema NNLS para u :

Minimizar

fEu − , sujeito a 0u ≥ ;

 Calcular o vetor fEur −= ;

 Se 0≠r , calcular

−=

x para nj ,,2,1 K

. Caso contrário, o sistema hGx ≥ é

inconsistente, ou seja, problema LDP não tem solução.

2.2 MATRIZ ESTEQUIOMÉTRICA

Nesta seção serão discutidas as propriedades estruturais das redes metabólicas cujos

elementos básicos são:

 compostos ou espécies bioquímicas e suas concentrações;

 reações ou processos de transporte que alteram as concentrações dos compostos

bioquímicos.

Os elementos acima descritos são sempre considerados variáveis em simulações de

redes metabólicas, chamadas de variáveis de estado (KLIPP

et al., 2005). Além destes,

existem ainda os parâmetros do modelo que influenciam a velocidade da reação como, por

exemplo, o nível da enzima, os metabólitos externos (efetores), os parâmetros cinéticos, entre

outros.

As reações bioquímicas e as espécies consideradas, com seus respectivos

coeficientes estequiométricos, compõem os elementos da chamada matriz estequiométrica,

A interseção de uma linha com uma coluna, na matriz estequiométrica, indica se

certa espécie faz parte ou não de uma reação em particular e, de acordo com o sinal do

elemento da matriz (isto é, de seu coeficiente estequiométrico), pode-se saber se este é um

reagente (valor positivo) ou um produto (valor negativo) da reação. E ainda, pela magnitude

do coeficiente sabe-se a quantidade relativa consumida ou produzida em uma determinada

reação. Desta forma, a estequiometria relaciona o número de moles de uma espécie química

envolvida em uma reação, mas não informa nada sobre a sua velocidade. A matriz

estequiométrica é constante e determinada pela constituição genética do organismo

(PALSSON, 2006). A

Figura 3 mostra a estrutura da matriz estequiométrica total.

⎟

⎠

⎞

⎜

⎝

⎛

Figura 3: Matriz estequiométrica total.

Palsson (2006) define a matriz estequiométrica total como mostrado na

Figura 3,

onde as linhas pontilhadas definem partições na matriz. Nesta forma de representação

identificam-se as reações internas (

), as reações de fronteira (v

), os compostos internos (x

) e

externos (

). Considera-se, para algumas análises, a célula viva como um sistema fechado e

trabalha-se apenas com os metabólitos internos, formando a matriz

N(m×n), conforme Figura

4. Esta matriz será usada em todas as análises e métodos descritos neste trabalho. Os

compostos externos (

), se necessários, são considerados constantes em todas as análises.

⎟

⎠

⎞

⎜

⎝

⎛

Figura 4: Matriz estequiométrica dos metabólitos internos.

Os metabólitos da rede formam o vetor coluna

x(m×1); os fluxos, o vetor coluna

v(n×1); e as concentrações das diferentes enzimas, o vetor coluna e(n×1), conforme equações

(2.20).

()

eee

vvv

xxx

,...,,

(2.20)

A título de exemplo, na

Figura 5a apresenta-se uma pequena rede metabólica com

quatro metabólitos internos e três externos. As setas, nesta figura, representam as reações, das

quais, três são de fronteira e as outras três são internas. A matriz estequiométrica total

correspondente é mostrada na

Figura 5b.

⎟

⎠

⎞

⎜

⎝

⎛

−

−−

−

100000

010000

000001

101000

010100

001110

000011

)

Figura 5: a) Exemplo de rede metabólica, onde os b

representam os metabólitos externos; v

e v

são as reações internas e v

, v

e v

as reações de fronteira. A linha pontilhada define o

contorno celular onde x

são os metabólitos internos e b

os produtos (b

e b

) e substrato (b

b) Matriz estequiométrica total da rede metabólica mostrada em (a).

2.3 SIMULAÇÃO DINÂMICA

Serão descritos nesta seção os dois principais métodos de simulação dinâmica de

redes bioquímicas e reguladoras: simulação determinística e simulação estocástica. No

primeiro, são calculados, ao longo de um intervalo de tempo, os níveis de concentração dos

compostos que formam a rede; no segundo, determinam-se o número de moléculas dos

compostos da rede ao longo do tempo.

2.3.1 Simulação Determinística

A variação da concentração dos metabólitos ao longo do tempo é realizada pelo

balanço de massa dinâmico, dado pela equação

(2.21), para os m metabólitos internos de rede

(KLIPP

et al., 2005).

mi para v

jij

,...,1

∑

(2.21)

Na equação

(2.21),

é o coeficiente estequiométrico do metabólito interno i na

reação

j. A equação (2.21) representa um conjunto de m equações diferenciais ordinárias

(EDO). Cada uma das velocidades

deve ter uma expressão matemática que determina a

velocidade desta reação em função das concentrações dos reagentes, produtos (se for uma

reação reversível) e dos parâmetros da reação (ver seção

2.5).

Escrevendo a equação

(2.21) em notação matricial, obtém-se a equação (2.22), onde

N é a matriz estequiométrica correspondente aos metabólitos internos. Esta equação, que

descreve o balanço de massa dos diversos metabólitos, é a base para qualquer simulação de

redes metabólicas.

(2.22)

A velocidade da reação (

) é dada em função dos níveis de concentração dos

metabólitos e dos valores dos parâmetros (parâmetros cinéticos, efetores externos, nível da

enzima, entre outros fatores) que aparecem na equação cinética da reação, equação

(2.23).

Nesta equação, o vetor

p(w×1) contém w parâmetros (ver seção 2.5).

()

px,fv

para

.,,2,1 nj K

(2.23)

2.3.2

Simulação Estocástica

O uso de EDOs para descrever processos bioquímicos faz certas suposições que nem

sempre são justificadas. Uma suposição é que as variáveis possuem valores contínuos. Isto é

uma simplificação, já que as espécies biológicas, moléculas, têm uma natureza discreta

(MCADAMS e ARKIN, 1997; YU

et al., 2006). Se o número de moléculas é suficientemente

grande não há problema com as EDOs, mas se o número de moléculas envolvido na

simulação é da ordem de dezenas ou centenas, essa natureza discreta deve ser levada em

consideração. Outra suposição importante no uso das equações diferenciais é que estas tratam

o processo descrito como sendo determinístico, isto é, as EDOs não levam em consideração as

flutuações aleatórias (MCADAMS e ARKIN, 1997). Novamente, esta suposição não se

sustenta para sistemas com pequenas quantidades de moléculas.

A solução para ambas as limitações é usar simulação estocástica que, explicitamente,

calcula a variação, ao longo do tempo, do número de moléculas das espécies envolvidas em

uma reação química. A seguir é apresentado um algoritmo que descreve a simulação

estocástica para um conjunto de reações químicas.

O algoritmo a seguir está baseado no trabalho de Gillespie (1977). Dado um vetor de

N espécies moleculares S = {S

,...,S

} com uma população inicial X

= {X

,...,X

} e um

vetor de

M reações químicas R = {R

,...,R

}. O algoritmo de simulação estocástica prediz a

evolução no tempo da população

X = {X

,...,X

} num intervalo de tempo de t

(tempo inicial)

end

(tempo final). Cada reação R

é definida pela sua constante de velocidade c

e sua

estequiometria, representada por dois vetores de coeficientes

e p

, onde cada r

ou p

é o

coeficiente do

k-ésimo reagente ou produto da i-ésima reação, Mi ,,2,1 K

. Num dado

intervalo de tempo infinitesimal

dt, a probabilidade de ocorrer a reação R

, dado o estado atual

X, é a

dt. A função a

é chamada de propensidade e é calculada através do produto da

constante

e o número de combinações possíveis entre os reagentes da reação R

. Por

exemplo, sejam as reações:

→

()

se tem 2:

se tem:

22542

21113211

−

caSSR

XXcaSSSR

Na reação R

, o número de combinações entre as moléculas dos reagentes é X

na reação R

(

)

−XX . Uma fórmula genérica para o cálculo da propensidade é dada

pela equação

(2.24), onde

⎟

⎠

⎞

⎜

⎝

⎛

é o coeficiente binomial (ULLAH et al., 2006).

∏

⎟

⎠

⎞

⎜

⎝

⎛

(2.24)

A partir da população inicial X

, o algoritmo de simulação estocástica simula a

evolução de cada estado X respondendo às seguintes perguntas:

 Em que tempo τ a próxima reação irá ocorrer?

 Qual será a próxima reação (índice j) a ocorrer?

A equação

(2.25) responde à primeira pergunta.

⎟

⎠

⎞

⎜

⎝

⎛

(2.25)

O índice j da próxima reação a acorrer é dado pelo menor inteiro que satisfaz a

equação

(2.26).

ara

∑

(2.26)

A população é atualizada pela equação

(2.27).

(

)()

prtXtX

−

(2.27)

Nas equações

(2.25) e (2.26), os símbolos r

e r

são números aleatórios

uniformemente distribuídos entre 0 e 1 e a

é definido pela equação (2.28).

∑

(2.28)

A seguir, no

Quadro 1, é exibido o algoritmo conhecido como método direto

(GILLESPIE, 1977).

Quadro 1: Algoritmo do método direto, proposto por Gillespie (1977).

Este algoritmo, como apresentado acima, é ineficiente em termos de tempo de

execução. Há outras abordagens que visam um melhor desempenho do algoritmo de

simulação estocástica (GIBSON e BRUCK, 2000; RAO e ARKIN, 2003; BURRAGE

et al.,

2004; TURNER

et al., 2004; CAO et al., 2006; MCCOLLUM et al., 2006).

As constantes de velocidades usadas para simulação estocástica e as usadas para

simulação determinística são relacionadas, mas não idênticas. As constantes determinísticas

dependem de concentração e as estocásticas do número de moléculas. A análise dimensional

mostra como as constantes de velocidades clássicas, para as reações de primeira e segunda

ordem, podem ser transformadas apropriadamente para serem usadas em simulações

estocásticas (GILLESPIE, 1977).

Dada uma reação de primeira ordem,

→ . A velocidade da reação, v, é dada por

kXv = , onde k é a constante de velocidade. A dimensão de v em uma descrição

determinística é

mol

sLs

mol 1

= e para a estocástica é

moléculas

moléculas 1

= . A dimensão

da constante de velocidade k é, em ambos os casos, 1/s. Portanto, não há necessidade de

conversão para a simulação entre as duas abordagens.

Agora, imagine uma substância P que é produzida pela reação das moléculas X e Y,

esquematicamente,

PYX

→+

. A velocidade da reação é dada por

kXYv = . A dimensão da

equação determinística é

mol

mols

mol

= e, no caso estocástico, é

1. Inicializar t, X

e o gerador de números aleatórios;

2. Computar a propensidade a

para cada reação R

usando a

equação (2.24) (i=1…M) e a

a partir da equação (2.28);

3. Gerar 2 números aleatórios r

e r

;

4. Determinar τ de acordo com a equação (2.25);

5. Determinar o menor inteiro i (índice da reação ocorrida)

que satisfaz a equação (2.26);

6. Fazer t = t + τ, e atualizar X(t) usando r

e p

, equação

(2.27);

7. Se t>t

end

ou a

=0 então pare, senão vá para o passo 2.

moléculas

moléculasss

moléculas

= (KLIPP et al., 2005). Para converter a constante de

velocidade de

mols

para

moléculass ×

, o valor numérico tem que ser dividido pelo

volume da reação e pela constante de Avogadro (para converter moles em moléculas). Por

exemplo, uma constante de velocidade de

−−

sM , que foi medida em um volume de

−

é convertido para o valor

11-9

101,66

−−

× smoléculas . De forma geral, deve-se usar a fórmula

(2.29)

onde, c é a constante estocástica, k a constante determinística (de segunda ordem), A é a

constante de Avogadro e

V o volume da reação.

2.4 INFORMAÇÕES ESTRUTURAIS DA MATRIZ ESTEQUIOMÉTRICA

A matriz estequiométrica N (Figura 4) contém informações importantes sobre a

estrutura da rede metabólica. Através desta matriz pode-se determinar:

 os fluxos possíveis no regime estacionário (espaço de fluxos admissíveis);

 os trechos não ramificados da via;

 as relações de conservação.

Pode-se determinar, ainda, outras matrizes a partir da matriz estequiométrica

como as matrizes de ligação de metabólitos e de ligação (

link) de fluxos. A seção 2.6 detalha

o uso destas matrizes na análise de controle metabólico.

2.4.1 Espaço Nulo da Matriz Estequiométrica

No regime estacionário não há variação da concentração dos metabólitos internos.

Assim, pode-se escrever a equação

(2.22) da seguinte forma:

0Nv =

(2.30)

Esta equação define um sistema de equações lineares com variáveis

(fluxos) que

tem solução não trivial somente se posto(

N) < n.

Determina-se como matriz

kernel K (espaço nulo à direita, ou simplesmente espaço

nulo) como sendo a matriz que satisfaz a equação

NK = 0 (REDER, 1988). A matriz K

contém todos os

n – posto(N) vetores bases do espaço nulo de N. A escolha da matriz K não é

única, uma vez que qualquer conjunto de

n – posto(N) vetores linearmente independentes com

dimensão

n×1 também satisfaz esta equação. Todas as possíveis distribuições de fluxos, no

regime estacionário, podem ser expressas como uma combinação linear das colunas

de K

(SCHILLING

et al., 1999), como mostra equação (2.31).

∑

−

)(

poston

(2.31)

Os coeficientes

, na equação (2.31), devem ter unidade de fluxo, por exemplo,

mM(min)

-1

. Na Figura 6 é exibida uma rede metabólica com quatro metabólitos e seis

reações. O posto da matriz

N(4×6) correspondente é igual ao seu número de linhas (=4).

Assim, a dimensão do espaço nulo é 2 (=6 – 4); logo, há duas colunas na matriz

kernel K (ver

Figura 7.)

Figura 6: Rede metabólica com 4 metabólitos e 6 reações (PALSSON, 2006).

⎟

⎠

⎞

⎜

⎝

⎛

−

Figura 7: Matriz kernel de N correspondente a rede da Figura 6.

Os dois vetores colunas da Figura 7 formam uma base para o espaço nulo de N. Por

exemplo, fazendo-se

= 2 e α

= 1, pela equação (2.31), determina-se uma distribuição de

fluxos

v como sendo uma combinação linear de α

() ()

⎟

⎠

⎞

⎜

⎝

⎛

⎟

⎠

⎞

⎜

⎝

⎛

−

⎟

⎠

⎞

⎜

⎝

⎛

O problema com a matriz

K é que o fluxo nas reações irreversíveis pode dar um

valor negativo. O exemplo da

Figura 6 mostra que todas as reações são irreversíveis e,

portanto, é desejável que os vetores bases tenham apenas valores positivos (SCHUSTER

al.

, 1999). Este fato é melhor representado na Figura 8, onde o segundo vetor da matriz K tem

sentido oposto ao sentido das reações. Na próxima seção será visto como calcular estes

vetores bases sem violar a condição de irreversibilidade.

Figura 8: Vetores bases do espaço nulo da matriz estequiométrica representados como uma

distribuição de fluxos na rede (PALSSON, 2006).

O conjunto de vetores colunas de K que possui o mesmo valor em cada linha forma

um ramo linear na rede metabólica (KLIPP

et al., 2005). Portanto, analisando a matriz K da

Figura 7, constata-se que os pares de reações (2, 3) e (5, 6) formam um ramo linear da via.

Isto pode ser conferido na

Figura 6 e na Figura 7.

2.4.2 Modos Elementares de Fluxos e Vias Extremas

Como visto na seção anterior, os vetores bases do espaço nulo (base linear) podem

conter valores negativos para as reações irreversíveis, o que não é viável

termodinamicamente. Particionando o vetor

v em reações reversíveis (v

rev

) e irreversíveis

(

irr

) e restringindo v

irr

para que seja positivo obtém-se a equação (2.32) (SCHILLING et al.,

2000; SCHUSTER

et al., 2002; PALSSON, 2006).

irr

(2.32)

O problema agora é encontrar um conjunto de vetores com dimensão

n×1 que

satisfaça as equações

(2.30) e (2.32) simultaneamente. Para resolver este problema é usada a

teoria matemática de análise convexa (SCHUSTER

et al., 1999). Esta abordagem é usada por

diversos algoritmos para determinar as vias possíveis. As análises recentes têm se focado em

dois conceitos: vias extremas (SCHILLING

et al., 2000) e modos elementares (SCHUSTER

et al., 2002; URBANCZIK e WAGNER, 2005). Ambos os algoritmos retornam as arestas de

um cone – a base convexa – como vias da rede metabólica.

Figura 9 são exibidas uma rede metabólica e as suas respectivas vias extremas

(VE) e modos elementares (ME). Nesta figura, cada ME ou VE indica uma via para a

biossíntese de um produto a partir de um substrato e o sentido das reações irreversíveis não é

violado. Vale lembrar que a estequiometria também é calculada para uma via definida como

um modo elementar ou uma via extrema.

Geralmente, existem, para uma determinada via, mais ME do que VE. A diferença

entre os dois tipos de vias está na representação das reações de fronteira (SCHILLING

et al.,

2000). Se estas reações são todas irreversíveis, os ME e VE são equivalentes. Se as reações de

fronteira são reversíveis há mais modos elementares do que vias extremas, como, por

exemplo, na

Figura 9.

Figura 9: Uma rede metabólica (A) e seus modos elementares (B) e vias extremas (C)

correspondentes (PALSSON, 2006).

A determinação das VE e ME para pequenas redes metabólicas é relativamente fácil.

Entretanto, à medida que o tamanho da rede aumenta, o número de VE e ME cresce de forma

exponencial. Portanto, a determinação das VE e ME relacionadas a uma rede é um problema

NP-difícil (KLAMT e

STELLING, 2002).

Os modos elementares podem ser usados para compreender as vias metabólicas em

uma rede, testar um conjunto de enzimas com vista à produção de um determinado metabólito

de interesse, detectar vias não redundantes, reconstruir o metabolismo a partir do genoma

anotado e analisar o efeito de enzimas deficientes (SCHILLING

et al., 1999; SCHUSTER et

al.

, 2000; STELLING et al., 2002; BARRETT et al., 2006).

2.4.3 Relações de Conservação

Existem subgrupos moleculares que são conservados durante a evolução de uma rede

metabólica, isto é, suas concentrações não variam ao longo do tempo (SCHUSTER e

HILGETAG, 1995). Por exemplo, a soma das concentrações de ATP, ADP e AMP pode ser

constante durante a evolução do sistema no tempo. Neste caso, o nucleotídeo adenosina (A) é

conservado (

ADP + AMP + ATP = constante1). Outro exemplo é o sistema formado pelo

NAD+ e o NADH. Neste caso, tem-se

NAD+ + NADH = constante2. Estas somas são

conhecidas como relações de conservação (SCHUSTER e HILGETAG, 1995; CORNISH-

BOWDEN e HOFMEYR, 2002). Cada conjunto de linhas linearmente dependentes (LD) da

matriz estequiométrica forma uma relação de conservação. Em termos computacionais, a

elucidação das relações de conservação em uma rede tem duas implicações práticas (SAURO

e INGALLS, 2004):

 a relação de conservação pode ser usada para reduzir o número de EDOs do sistema. Isto é

feito substituindo uma EDO pela equação algébrica correspondente;

 pode-se calcular diretamente as concentrações dos metabólitos no regime estacionário ao

eliminar as linhas linearmente dependentes da matriz

N(m×n). Sem a eliminação destas

linhas o sistema de equações não lineares pode ter infinitas soluções.

Se há relações de conservação na rede metabólica, então a matriz

N tem linhas LD, o

que implica que

r < m, onde r é o posto da matriz N. As linhas de N podem ser rearranjadas

de forma que as primeiras

r linhas sejam linearmente independentes (LI) e as últimas m – r

linhas LD. Da mesma forma, os metabólitos correspondentes a estas linhas são rearranjados

formando os metabólitos independentes

ind

(r×1) e os dependentes x

dep

(m-r×1). Assim, x

pode ser reescrito de acordo com a equação

(2.33):

⎥

⎦

⎤

⎢

⎣

⎡

dep

ind

(2.33)

Também deve-se particionar

N como descrito acima. Fazendo isso, obtém-se a

equação

(2.34)

⎥

⎦

⎤

⎢

⎣

⎡

(2.34)

Cada linha de

(m-r×n) na equação (2.34) é uma combinação linear das linhas de

(r×n). Por isso pode-se definir a matriz link zero L

(m-r×r) que satisfaz a equação (2.35)

(SAURO e INGALLS, 2004).

R00

NLN =

(2.35)

Pode ser demonstrado que (SAURO e INGALLS, 2004)

[]

onst

dep

ind

⎥

⎦

⎤

⎢

⎣

⎡

−

r-m

(2.36)

onde

m-r

é a matriz identidade de ordem m-r×m-r. Fazendo-se

[]

ILT

−

e usando

equação

(2.33), obtém-se

onst

CTx =

(2.37)

A matriz

T(m-r×m) é chamada de matriz de conservação e

onst

(m-r×1) = x

dep

(0)-

ind

(0) é um vetor determinado pelo valor inicial da concentração dos metabólitos.

Cada linha da matriz

T relaciona um ciclo conservado e, desta forma, o número de

linhas de

T indica o número de relações conservadas (linhas LD da matriz N) na rede. Os

elementos não zero de uma linha de

T indicam quais metabólitos contribuem para um ciclo

em particular.

A equação

(2.37) exibe a relação entre T e x(t), mas não mostra como calcular a

matriz

T. Existe uma variedade de abordagens relacionadas para determinar T (SAURO e

INGALLS, 2004). A seguir duas delas são descritas.

A mais simples delas é através do espaço nulo de

N. Reescrevendo a equação (2.35)

obtém-se a seguinte equação:

[]

⎥

⎦

⎤

⎢

⎣

⎡

−

r-m

(2.38)

Substituindo a equação

(2.34) nesta equação, tem-se

[]

0NIL

−

r-m

(2.39)

Observe que o primeiro termo desta equação é a matriz de conservação

T. Assim,

0TN =

(2.40)

Pela equação

(2.40) vê-se que as linhas de T formam a base do espaço nulo à

esquerda de

N. Como é mais convencional computar o espaço nulo (à direita) de N, faz-se

0TN =

(2.41)

Na equação

(2.41) o sobrescrito T denota a transposta da matriz correspondente.

Com isto pode-se determinar

T simplesmente encontrando uma base para o espaço nulo de

. No MATLAB

, por exemplo, basta fazer T = null(N’)’, onde o apóstrofo (’) representa a

operação transposta.

Outra forma de calcular a matriz de conservação

T é encontrar a forma reduzida por

linha de

N. Pode-se transformar a matriz N em N

se premultiplicar a matriz N por uma série

de matrizes elementares

. Fazendo o produto destas matrizes elementares igual a M(m×m),

tem-se

⎥

⎦

⎤

⎢

⎣

⎡

(2.42)

Pode-se particionar

M como

⎥

⎦

⎤

⎢

⎣

⎡

, substituindo em (2.42), obtém-se

⎥

⎦

⎤

⎢

⎣

⎡

⎥

⎦

⎤

⎢

⎣

⎡

Logo, conclui-se que

0NM

(2.43)

Pela equação

(2.40) e (2.43) tem-se que M

= T. Para encontrar T deve-se aumentar

N com a matriz identidade

[]

IN . Escalonando esta matriz aumentada, por operações

elementares, a matriz identidade é transformada na matriz

M. Feito isto, a matriz de

conservação é obtida tomando-se as últimas

m – r linhas de M. A partir do momento que se

obtém a matriz

T, tem-se também a matriz L

e, consequentemente, a matriz L

, utilizada na

análise de controle metabólico para calcular os coeficientes de controle.

2.5 CINÉTICA DAS REAÇÕES BIOQUÍMICAS

Na seção 2.4 foi visto que a matriz estequiométrica só fornece as relações estruturais

da rede metabólica. A velocidade da reação – isto é, a quantidade de reagente, em unidade de

concentração, transformado em produto por unidade de tempo por uma reação enzimática –

depende de vários fatores intracelulares, como, por exemplo, pH, temperatura e concentração

dos reagentes (MARANGONI, 2003). Todas as reações no interior da célula são catalisadas

por enzimas que, por sua vez, são produtos da expressão dos genes da célula.

Dada a reação irreversível

PS → , catalisada pela enzima E, esta reação pode ser

expandida de forma a representar os eventos que ocorrem durante a catálise da reação:

primeiro o substrato se liga à enzima formando o complexo enzima-substrato (ES) e, depois, a

enzima forma o produto P. Por último, este produto é liberado. Esquematicamente:

PEESSE

cat

+→↔+

É possível – através da lei da ação das massas – elaborar a seguinte equação cinética,

em função do substrato S, para esta reação (LEHNINGER

et al., 2000):

max

(2.44)

A equação

(2.44) é conhecida como cinética de Michaelis-Menten (LEHNINGER et

al.

, 2000). Nesta equação, V

max

representa a velocidade máxima da reação no sentido

substrato-produto,

m,S

é a constante de Michaelis-Menten para o substrato S (afinidade de

ligação entre a enzima e o substrato) e

v é a velocidade da reação.

Da mesma forma, obtém-se a equação

(2.45) para a reação reversível PS ↔ :

⎟

⎠

⎞

⎜

⎝

⎛

−

max

(2.45)

Na equação

(2.45), K

é a constante de equilíbrio para esta reação, K

m,P

a constante

de Michaelis para o produto

P e V

max

como definido. No entanto, a maioria das reações

bioquímicas possui mais de um substrato e/ou mais de um produto.

Quando há mais de um substrato (ou produto), este pode se ligar de forma aleatória

ou ordenada à enzima. Além disso, a reação reversível é classificada segundo o número de

moléculas do substrato (ou produto) utilizadas na reação (uni, bi, ter, ...) (MARANGONI,

2003). Na equação

(2.46) é apresentada a cinética para um mecanismo com um substrato e

dois produtos,

QPS +↔

. Este é conhecido como mecanismo uni-bi ordenado. Nesta

equação,

max1

é a velocidade máxima no sentido substrato-produto e V

max2

a velocidade

máxima no sentido contrário.

m,S

, K

m,P

e K

m,Q

são as constantes de Michaelis para o substrato

S e os produtos P e Q, respectivamente. Além destes, há ainda outros mecanismos. Para mais

detalhes ver Leskovac (2003) e Marangoni (2003).

SVKV

SVV

eqPmeq

1max

2max

,1max

2max,2max

2max1max

+++++

⎟

⎠

⎞

⎜

⎝

⎛

−

(2.46)

Uma relação importante para reações reversíveis, conhecida como relação de

Haldane, é dada pela equação

(2.47)

∏

2max

1max

(2.47)

Esta equação relaciona todos os parâmetros de uma cinética enzimática. Ela permite,

também, eliminar

max2

da equação (2.46), como foi feito na equação (2.45).

2.6 SUBESPAÇOS RELACIONADOS À MATRIZ ESTEQUIOMÉTRICA

Matematicamente, a matriz estequiométrica N pode ser vista como uma

representação de uma transformação linear do vetor de fluxos

v (equação (2.20)) para o vetor

de derivadas, em relação ao tempo, da concentração

x (equação (2.20)) pela equação (2.22),

repetida a seguir (ver

Figura 10).

(2.22)

Existem quatro subespaços vetoriais associados à matriz estequiométrica. Estes

quatro subespaços de

N, mostrados na Figura 10, têm um papel importante na análise de redes

bioquímica (GIANCHANDANI

et al., 2006). Os vetores produzidos pela transformação

linear estão em dois subespaços ortogonais (espaço coluna e espaço nulo à esquerda),

chamado de imagem, e os vetores que estão sendo mapeados estão também em dois espaços

ortogonais (os espaço linha e nulo), chamados domínio da transformação (PALSSON, 2006).

Figura 10: Matriz estequiométrica vista como uma transformação linear. Os quatros

subespaços vetoriais associados a N também são mostrados. Adaptado de Palsson (2006).

As derivadas do espaço coluna de

N, Col(N), pode ser visto pela expansão

vvv

n21

nnn

+++= K

(2.48)

onde

é uma coluna de N. O espaço Col(N) é, portanto, gerado pelos vetores n

. Estes

vetores são as características estruturais da rede e são fixos (LLANERAS e PICÓ, 2008),

enquanto os fluxos

são escalares e variáveis.

Os vetores

do espaço nulo a esquerda, Lnull(N), são ortogonais ao espaço coluna

N, ou seja, 0, =

nl . Os vetores l

representam a conservação da massa.

Os vetores de fluxo podem ser decompostos em 2 componentes: dinâmico e

estacionário:

ssdin

vvv +

(2.49)

Null(N)

Row(N)

Lnull(N)

Coll(N)

din

ℜ

O componente estacionário satisfaz a equação

(2.30), 0Nv

, e assim, v

está no espaço

nulo de

N. O componente dinâmico do vetor de fluxo, v

din

, é ortogonal ao espaço nulo de N e,

consequentemente, ele está no espaço linha de

N. Cada par de subespaços no domínio e

imagem da equação dinâmica do balanço de massa formam conjuntos ortogonais entre si.

2.6.1 Conteúdo dos Subespaços da Matriz Estequiométrica

Os quatros subespaços fundamentais contêm informações sobre uma rede

bioquímica. Seus conteúdos podem ser assim descritos (PALSSON, 2006):

 Espaço nulo. O espaço nulo de N contém todas as distribuições de fluxos no regime

estacionário permitidas por uma rede metabólica. Ver seção

2.4.1;

 Espaço linha. O espaço linha de N contém todas as distribuições dinâmicas dos fluxos de

uma rede e, desta forma, as forças termodinâmicas que alteram a variação das atividades

das reações;

 Espaço nulo à esquerda. O espaço nulo à esquerda de N contém todas as relações de

conservação da rede. Ver seção

2.4.3;

 Espaço coluna. O espaço coluna da matriz N contém todos os possíveis vetores de

derivadas da concentração no tempo e, assim, mostra como as forças termodinâmicas

governam o estado das concentrações dos metabólitos da rede.

2.7 ANÁLISE DE CONTROLE METABÓLICO

Um dos objetivos da análise de controle metabólico (ACM) é medir a sensibilidade

das variáveis – fluxos e concentrações (os controlados) – no regime estacionário em relação

às variações dos parâmetros (controladores) das reações (FELL, 1992). Estas sensibilidades

são propriedades globais, uma vez que dependem de todas as interações dentro do sistema

celular. Além disso, deve-se relacionar tais propriedades globais às propriedades locais de

cada reação.

O termo controle é usado aqui no sentido de ser capaz de influenciar alguma coisa.

Quando se diz que uma enzima (reação) controla o fluxo em uma via metabólica, significa

que uma variação na atividade da enzima resulta numa variação do fluxo. Isso não significa

que a enzima exerce um controle completo, ou nenhum, sobre o fluxo ou sobre a concentração

dos metabólitos na rede (FELL, 1998). Todas as reações exercem um controle sobre a via à

qual pertencem, algumas com maior ou menor grau do que as outras. Os coeficientes de

controle (sensibilidade) representam uma forma de medida deste grau de controle (influência)

sobre uma via metabólica (KACZER e BURNS, 1973).

O relacionamento entre as variáveis (fluxos e concentrações) e os parâmetros que

afetam a atividade das enzimas num sistema é não linear e não pode ser obtido analiticamente

(FELL, 1998).

Em ACM são definidos três coeficientes de controle adimensionais, um local

(coeficiente de elasticidade) e dois globais (coeficientes de controle de fluxo e de

concentração).

2.7.1 Coeficiente de Elasticidade

A atividade de uma enzima é diretamente influenciada por metabólitos específicos

(substratos, produtos, inibidores e ativadores). Estes metabólitos são aqueles que aparecem na

equação da cinética da enzima. O coeficiente de elasticidade (CE) quantifica a sensibilidade

da velocidade local (atividade) de uma enzima a pequenas variações dos metabólitos que

influenciam, de alguma forma, a atividade da enzima (KACZER e BURNS, 1973). Note que

este conceito é aplicado a cada enzima separadamente (local). Formalmente, o coeficiente de

elasticidade

é definido como (STEPHANOPOULOS et al., 1998)

(

)

()

∂

(2.50)

Cada número

determina o coeficiente de elasticidade da j-ésima enzima em

relação ao i-ézimo metabólito. Analisando esta última equação, percebe-se que o coeficiente

de elasticidade é a inclinação da tangente à curva obtida traçando-se o gráfico de ln(v

) contra

ln(x

) (ver Figura 11). O CE descreve uma propriedade isolada de uma unidade funcional

(enzima) do sistema, não informando nada sobre o comportamento global do sistema. Este

coeficiente será usado para determinar os coeficientes de controle.

ln(x

)

ln(v

) ln(v

)

ln(I)

Figura 11: Coeficiente de elasticidade (ε) de uma reação v

em relação ao substrato x

e um

inibidor I.

2.7.2 Coeficientes de Controle

Os coeficientes de controle são uma medida de como a variação do nível de uma

enzima afeta o fluxo (coeficiente de controle de fluxo) ou a concentração de um metabólito

(coeficiente de controle de concentração) (FELL, 1992). Ambos os coeficientes são definidos

em relação a um regime estacionário de referência com níveis de enzima e

, níveis de

metabólitos x

e fluxo global J

constantes. Todos os coeficientes são adimensionais. O

coeficiente de controle de fluxo (CCF) é definido na equação

(2.51) e o coeficiente de

controle de concentração (CCC) na equação

(2.52).

(

)

()

para k,j = 1,...,n

(2.51)

(

)

()

== para j = 1,...,n e i = 1,...,m

(2.52)

As relações entre estes coeficientes são conhecidas como teorema da soma e teorema

da conectividade (KACZER e BURNS, 1973). De acordo com o teorema da soma, a soma dos

CCF de uma via é igual a um e a soma do CCC é igual a zero. Os teoremas de conectividade

relacionam o coeficiente de elasticidade, CE, aos coeficientes de controle CCC e CCF.

Combinados, estes teoremas formam um sistema de equações lineares que permite a

determinação dos coeficientes de controle (propriedade global) a partir dos coeficientes de

elasticidade (propriedade local) (KACZER e BURNS, 1973).

2.7.3 Relações Estruturais

Na seção 2.4.3 (relações de conservação) foi definida a matriz L

. Combinando L

com

(matriz identidade de ordem r) obtém-se a matriz de ligação L

(m×r) dada pela

equação

(2.53).

⎥

⎦

⎤

⎢

⎣

⎡

(2.53)

A matriz

é conhecida como matriz de ligação de metabólitos (REDER, 1988).

Com esta matriz é possível determinar as concentrações dos metabólitos dependentes em

função dos metabólitos independentes, conforme a equação

(2.54):

ind

xLx =

(2.54)

A relação de conservação é derivada como (EHLDE e ZACCHI, 1997):

ind

(2.55)

Integrando a equação

(2.55), obtém-se:

(

)()

AxLx

ind

+= tt

(2.56)

onde

A é um vetor m×1 no qual os primeiros r (posto da matriz N) elementos das linhas são

zeros.

A relação entre os fluxos é derivada de forma similar. Como foi visto na seção

2.4.1,

a equação

(2.30) possui soluções não triviais quando posto(S) < n, isto é, se há dependência

linear entre as colunas de

N. Estas podem ser expressas usando a matriz de ligação de fluxos

(EHLDE e ZACCHI, 1997).

()

ind

JLJ

null

(2.57)

Ehlde e Zacchi (1997) combinaram as relações estruturais com a análise de

sensibilidade geral e derivaram uma fórmula para o cálculo dos coeficientes a partir do qual

eles chamaram de relações cinéticas. A derivação começa a partir da expressão geral para a

velocidade de uma reação, equação

(2.23). Derivando esta equação para v

em relação ao

nível da enzima j, tem-se

∂

(2.58)

Como os parâmetros em p não são acoplados, define-se

⎪

⎩

⎪

⎨

⎧

≠

ep se

(2.59)

Multiplicando a equação

(2.58) por

ve , e levando-se em consideração a equação

(2.59), obtém-se

∑

⎟

⎠

⎞

⎜

⎝

⎛

∂

(2.60)

Se v

for considerado o fluxo no regime estacionário da reação k (v

= J

), pode-

se facilmente reconhecer o coeficiente de elasticidade, equação

(2.50), e os coeficientes de

controle, equações

(2.51) e (2.52). Dessa forma pode-se reescrever a equação (2.60)

()

∑

επ

(2.61)

Em notação matricial:

εCπC +=

(2.62)

Como a velocidade de uma reação enzimática é proporcional à concentração da

enzima, a matriz

π pode ser substituída pela matriz identidade, pois:

⎩

⎨

⎧

≠

kj se

(2.63)

A equação

(2.62) fornece uma relação entre todos os coeficientes definidos até

agora. Nesta equação tem-se

n×n + m×n incógnitas – número de coeficientes de controle – e

apenas

m×n equações. Há relações adicionais que podem ser derivadas das relações

estruturais, equações

(2.56) e (2.57). Pode ser mostrado que (EHLDE e ZACCHI, 1997):

ind

CLC =

(2.64)

ind

CLC =

(2.65)

onde

()

nrn ×−

ind

C e

()

mr ×

ind

C contêm os coeficientes de fluxo e de concentração para as

linhas independentes de

N, respectivamente, e

()

(

)

ind

XLXL diagdiag

1−

(2.66)

()

(

)

ind

JLJL diagdiag

1−

(2.67)

As equações

(2.64) e (2.65) podem ser inseridas na equação (2.62) obtendo-se,

finalmente, a equação

(2.68)

[]

πεLL

ind

1−

−=

⎥

⎦

⎤

⎢

⎣

⎡

(2.68)

onde é possível determinar os CCC (

ind

) e CCF (

ind

) para as linhas independentes

(metabólitos) da matriz estequiométrica. Substituindo estes valores nas equações

(2.64) e

(2.65) determinam-se, portanto, todos os m×n CCC e os n×n CCF de uma via metabólica.

2.8

ANÁLISE DE FLUXO METABÓLICO

A análise de fluxo metabólico (AFM) é um método usado para a determinação

estequiométrica dos fluxos internos de uma rede metabólica no regime estacionário

(STEPHANOPOULOS

et al., 1998). Isto é, na AFM são calculados os fluxos não medidos,

geralmente os fluxos internos, dados alguns fluxos medidos experimentalmente; estes são

geralmente fluxos externos. A equação

(2.30) forma a base para a análise de fluxo metabólico.

Seja v

os fluxos desconhecidos – os fluxos que serão calculados – e v

os fluxos

calculados. Logo, pode-se reescrever a equação

(2.30) segundo o particionamento descrito

acima para se obter a equação

(2.69).

⎥

⎦

⎤

⎢

⎣

⎡

⎥

⎦

⎤

⎢

⎣

⎡

(2.69)

Rearranjando a equação

(2.69) e isolando, na equação resultante, o termo

vN ,

tem-se uma expressão para os fluxos não medidos, equação

(2.70):

ccxx

vNvN −=

(2.70)

A equação

(2.70) forma uma sistema de equações lineares que pode ser (NIELSEN e

VILLADSEN, 1994):

 determinado, se N

for quadrada e não singular, ou seja, N

não pode ter linhas LD.

Assim, a solução é dada por:

()

ccxx

vNNv

1−

−=

(2.71)

 sobredeterminado, se o número de fluxos medidos for maior que o grau de liberdade do

sistema. Neste caso, a matriz N

não é quadrada e não pode ser invertida. Usa-se, então,

um dos métodos do PQM ou a pseudo-inversa (PENROSE, 1955). Com isso, pode-se

reescrever a equação

(2.70):

()( )

ccxx

vNNv

−=

(2.72)

onde

()

N representa a pseudo-inversa de N

;

 indeterminado, quando a quantidade de fluxos medidos é menor que o grau de liberdade

do sistema, ou seja, o número de variáveis é maior que o número de equações do sistema.

Sendo assim, a equação

(2.70) possui infinitas soluções. Uma forma de determinar uma

solução única para este caso é definir uma função objetivo – que deve ser linear – e

otimizar esta função de forma a obter um problema de programação linear, caso similar ao

descrito na seção

2.9. Também, pode-se encontrar

v usando um dos métodos do PQM

ou a pseudo-inversa, equação

(2.72). Nestes casos,

é uma solução exata e com a menor

distância da origem.

Para mais detalhes sobre AFM, bem como para exemplos de uso de AFM,

recomenda-se Nielsen e Villadisen (1994) e Stephanopoulos

et al. (1998).

2.9

ANÁLISE DE BALANÇO DE FLUXO

A análise de balanço de fluxo (ABF), por vezes chamada de análise baseada em

restrições, investiga a capacidade teórica e modos de operação do metabolismo envolvendo

restrições na análise da estequiometria (VARMA e PALSSON, 1994a).

A FBA assume que a rede metabólica atingirá o regime estacionário restrito pela

estequiometria, equação

(2.30). Nesta equação,

nm×

ℜ∈N é a matriz estequiométrica e

1×

ℜ∈

v o vetor de fluxos.

são, respectivamente, o número de metabólito e de reações

(fluxos) da rede. Este sistema, na grande maioria das vezes, é indeterminado e possui várias

soluções, o que implica na formação de um espaço de soluções (KAUFFMAN

et al., 2003).

Pode-se adicionar mais restrições a este espaço especificando-se valores máximos ou mínimos

para os fluxos da rede:

jjj

vvv

maxmin

≤≤ com nj L1

(2.73)

Uma ilustração destas três restrições é apresentada na

Figura 12. Cada uma das

restrições imposta a rede reduz o espaço de distribuição dos fluxos.

Figura 12: Ilustração das restrições aplicadas a uma rede metabólica. Todas as restrições

impostas à rede determinam o espaço de distribuição de fluxo permissível, com isso cada

restrição reduz o tamanho do espaço solução.

Estas restrições confinam os fluxos do regime estacionário a um conjunto de

soluções viáveis – a equação

(2.30) tem infinitas soluções quando n>m ou quando n = m e

posto(N) < m. Mas, usualmente, não se consegue uma única solução para a distribuição dos

fluxos. A determinação de uma distribuição de fluxos em particular é formulada como um

problema de programação linear (VARMA e PALSSON, 1994b), como visto acima.

Para aquelas reações que não possuem valores máximos e/ou mínimos conhecidos

atribuem-se, para as reações irreversíveis,

min

v e

∞

max

v e, para as reversíveis,

−∞=

min

v e

∞

max

v . Dessa forma, a direção das reações (restrição termodinâmica) são

respeitadas. Estas restrições definem a capacidade da rede metabólica (EDWARDS

et al.,

2002). Há várias soluções para os fluxos neste espaço. O objetivo da FBA é encontrar uma

solução, neste espaço, que otimiza uma função objetivo (VARMA e PALSSON, 1994a):

A determinação de uma distribuição de fluxos em particular é formulada como um

problema de programação linear (VARMA

e PALSSON, 1994b). Geralmente, o objetivo é

Nv = 0 v

≥ v

min,

≤ v

max,

maximizar a velocidade de crescimento ou otimizar (max/min) alguns dos

n fluxos da rede

metabólica (LLANERAS e PICÓ, 2008). A idéia é maximizar uma função objetivo

Z que está

sujeita à estequiometria e às restrições de capacidade. A função

Z é definida pela equação

(2.74), onde os w

representam os pesos para as reações individuais. Alguns exemplos usuais

de funções objetivos são: maximização da produção de ATP, minimização do consumo de

nutrientes, maximização do rendimento de um produto desejado, maximização da velocidade

de crescimento, ou uma combinação destas.

máxvwZ

→=

∑

(2.74)

2.10

VISÃO GERAL DA MODELAGEM ESTEQUIOMÉTRICA

A equação geral

(2.30) fornece um conjunto de restrições que ligam os fluxos na rede

metabólica e restringem o espaço de possíveis distribuições de fluxos a um hiperplano,

subespaço

ℜ , onde os eixos correspondem aos fluxos das reações num rede (Figura 13).

Contudo, aos fluxos das reações irreversíveis são impostos à restrição de não negatividade

(por convenção, a direção natural do fluxo é considerada positiva):

≥

irr

(2.75)

Uma vez que estas restrições são de desigualdade, o sistema resultante não pertence

ao escopo da álgebra linear. No entanto, pela análise convexa, mostra-se que o espaço da

distribuição de fluxos formado pelas equações

(2.30) e (2.75) é um cone poliédrico convexo,

conforme a

Figura 13.

Finalmente, pode-se impor os valores de velocidades máximas derivados das

enzimas ou da capacidade de transporte como mais uma restrição:

max

≤

(2.76)

Se essa informação estiver disponível para todos os fluxos da rede, o espaço de

possíveis soluções fica limitado (na verdade, é suficiente que cada aresta do cone tenha seu

fluxo restrito). Em termos matemáticos, o cone poliedro convexo foi convertido em um

politopo (

Figura 13). As equações (2.30), (2.75) e (2.76) representam a restrição invariante no

tempo mais comum. Elas definem um espaço dentro do qual está qualquer distribuição de

fluxo viável. Assim, elas definem um modelo da capacidade do metabolismo que está sendo

estudado.

É possível incorporar mais restrições para restringir mais ainda o espaço de

distribuição de fluxos, inclusive variante no tempo, ou até mesmo prever a própria

distribuição de fluxos. Por exemplo, redes lógicas podem ser incorporadas como restrições

reguladoras (COVERT

et al., 2001).

Figura 13: Esquema das diferentes aplicações de análise estequiométrica. Esta análise pode ser

classificada em duas principais categorias: aquelas que focam nas propriedades do espaço

como um todo das possíveis distribuições de fluxos (ex. ME e VE) e aquelas que determinam

uma distribuição de fluxo particular (ex. ABF e AFM) (LLANERAS e PICÓ, 2008).

2.11

PRODUÇÃO DE HIDROGÊNIO BACTERIANO

O hidrogênio está sendo visto como um componente importante da nova matriz

energética como um portador de energia eficiente e limpa. O hidrogênio molecular, H

, tem o

maior conteúdo de energia por unidade de peso entre todos os combustíveis gasosos

conhecidos e é o único combustível que não possui carbono em sua fórmula (SHAMS e

GONZALEZ, 2007). Estas observações, juntamente com o fato de que os mais importantes

métodos industriais de produção de hidrogênio ainda usam fontes de energia não renovável,

fornece argumentos significativos para a busca de uma forma economicamente viável e limpa

para produção de hidrogênio.

Orientada à

análise do

sistema

Solução da

distribuição de

fluxos

Aplicações da

modelagem

este

uiométrica

Determina o

fenóti

o corrente

Constrói modelos

reditivos

Álgebra linear

(

o nulo

)

Análise Convexa

(

ME e VE

)

Análise de fluxo

metabólico

Espectro de

fluxos

Análise de

balan

o de fluxo

{

0Nv

⎩

⎨

⎧

≥

0Nv

⎩

⎨

⎧

0Nv

⎪

⎩

⎪

⎨

⎧

≤≤

≤

≥

max,min,

max

vvv

0Nv

⎪

⎩

⎪

⎨

⎧

⋅=

≤≤

≤

≥

0Nv

vvv

:max

max,min,

max

A produção de biohidrogênio a partir de microrganismos é favorável pelo fato do

processo de produção ocorrer em temperatura e pressão ambientes. Sendo assim, este

processo requer menos energia que os métodos convencionais (VARDAR-SCHARA e

MAEDA, 2008). Os desafios técnicos para se conseguir tal objetivo incluem baixar os custos

do processo de produção de hidrogênio e melhorar o rendimento destes processos.

O uso de glicerol, gerado como subproduto na produção de biodiesel, representa uma

alternativa promissora para se ter uma viabilidade econômica (SHAMS e GONZALEZ,

2007).

A produção biológica de H

é especialmente atrativa quando se usa o glicerol cru,

isto é, sem qualquer purificação, como fonte de carbono para o microrganismo. O

melhoramento e o desenvolvimento de tecnologias para a produção de biohidrogênio

dependem, principalmente, do uso de microrganismos modelos bem conhecidos. As espécies

do gênero

Clostridium são os mais estudados em processos de fermentação anaeróbicos.

Pelo ponto de vista de bioengenharia, a

C. acetobutylicum é a espécie mais estudada

deste gênero. Os estudos desta bactéria concentram-se no melhoramento da mesma para a

produzir solventes, já que ela é conhecida pela sua habilidade de produzir solventes (JONES e

WOODS, 1986).

A fermentação anaeróbica realizada por esta bactéria tem duas fases: a fase

acidogênica, quando o crescimento é rápido, com maior produção de hidrogênio e os ácidos

acetato e butirato; e a fase solventogênica, quando o crescimento da bactéria é menor e há

diminuição da produção de H

. Nesta fase há a formação de solventes com um consumo

parcial dos ácidos formados durante a primeira fase (RAO e MUTHARASAN, 1987).

2.11.1

Simulação in silico da Produção de Biohidrogênio

Tem havido um esforço considerável de pesquisa na caracterização e manipulação

genética do gênero

Clostridium para elucidar e melhorar a produção de solventes, mas pouco

tem sido feito para a produção de H

. A caracterização e manipulação genética in silico deste

gênero pode ser explorada pelos métodos de engenharia metabólica e biologia de sistema, tal

como ABA e AFM. Estes métodos podem fornecer uma oportunidade para melhorar a

produção de H

e fornecer uma melhor compreensão do metabolismo da C. acetobutylicum.

Através desta abordagem, Senger e Papoutsakis (2008) construíram uma rede

metabólica em escala genômica da

C. acetobutylicum usando ABF para identificar um

conjunto de genes essenciais para o crescimento deste microrganismo, como, por exemplo, o

funcionamento do ciclo do ácido tricarboxílico que é incompleto e parcialmente invertido

nesta bactéria, segundo a anotação de seu genoma (NOLLING

et al., 2001). Manish et al.

(2007) e Oh

et al. (2008) aplicaram AFM e ABF numa rede metabólica da E. coli com o

intuito de melhorar a produção de H

nesta bactéria. Shinto et al. (2007) desenvolveram um

modelo cinético e efetuaram uma análise da produção de acetona-butanol-etanol baseado num

modelo da

C. acetobutylicum, mas não incluíram as reações para a produção de H

. A rede

que eles usaram é semelhante a da

Figura 23, página 82. Em outro trabalho, foi simulado um

modelo estequiométrico para a fermentação da

C. acetobutylicum semelhante a ABF, porém,

as restrições eram formadas por equações não lineares (DESAI; NIELSEN

et al., 1999).

Também foi feito uma AFM da fermentação da

C. acetobuylicum com o objetivo de analisar a

produção de solventes desta bactéria (DESAI; HARRIS

et al., 1999). Nenhuma destas

simulações considerou o glicerol como fonte de carbono.

2.12

SIMULAÇÃO DINÂMICA DO OPERON TRP NA ESCHERICHIA COLI

O triptofano é um dos aminoácidos aromáticos sintetizado por todos os

microrganismos e é utilizado na síntese de proteínas para o crescimento e reprodução celular.

Em algumas bactérias os genes que codificam as enzimas necessárias para a biossíntese do

triptofano estão agrupados formando o

operon trp (YANOFSKY et al., 1981). A via de

biossíntese do triptofano nos microrganismos possui 5 reações (enzimas), mas os genes que

codificam as diferentes subunidades de cada enzima são diferentes. A via do Trp é altamente

regulada naquelas bactérias que os sintetizam. Isto é devido ao fato de que a síntese de uma

molécula de Trp tem um custo energético muito alto (AKASHI e GOJOBORI, 2002).

Escherichia coli os principais mecanismos de regulação da quantidade de

triptofano no meio intracelular ocorrem de três formas distintas: repressão, atenuação

transcricional e inibição enzimática (GOLLNICK e BABITZKE, 2002). Desta forma,

qualquer simulação da via do triptofano deve levar em consideração todos estes mecanismos.

A abordagem conhecida como determinística é realizada através do uso de equações

diferenciais ordinárias. A lei fundamental que governa a velocidade das reações na simulação

determinística é a

lei de ação das massas. Essa lei diz que uma reação em um meio

homogêneo tem velocidade proporcional às concentrações individuais dos reagentes

envolvidos. Usando essa lei as expressões para as velocidades de alteração das concentrações

de cada uma das moléculas podem ser encontradas. Com isso, pode-se expressar qualquer

sistema químico como um conjunto de EDOs acopladas. Geralmente, não se encontra uma

solução analítica para a solução das EDOs; entretanto, pode-se usar um método numérico para

integrar o conjunto de EDOs ao longo do tempo para achar uma solução aproximada da

dinâmica do sistema.

Várias simulações determinísticas foram propostas para a simulação dinâmica do

operon trp da E. coli. Nos modelos de Goodwin (1965), Griffith (1968a) e Griffith (1968b)

havia apenas o mecanismo de repressão. Um outro modelo levou em consideração a inibição

enzimática e o retardo de tempo (BLISS

et al., 1982). As previsões do modelo de Bliss

mostraram uma concordância qualitativa com os resultados experimentais, exibindo

oscilações no nível de concentração do triptofano. Baseado em novos dados experimentais um

outro modelo foi desenvolvido, no qual a repressão ocorre em dois passos (SINHA, 1988).

O modelo de Sinha foi modificado usando-se a função cinética de Michaelis-Menten

para representar a demanda de triptofano para a síntese de proteína (SEM e LIU, 1990). Em

ambos os modelos a inibição enzimática pelo triptofano não foi considerada. Em outros

modelos investigou-se o efeito da velocidade de crescimento da bactéria na síntese de

triptofano e a estabilidade do

operon sob variações de vários parâmetros e em função da

concentração intracelular do próprio triptofano (XIU

et al., 1997; XIU et al., 2002).

Outros modelos recentes e mais completos do

operon trp (SANTILLAN e

MACKEY, 2001b; 2001a; SANTILLAN e ZERON, 2004) levaram em consideração todos os

três mecanismos de regulação e o retardo de tempo.

Nestes modelos, os autores tiveram uma consideração especial na estimação dos

parâmetros e os resultados numéricos do sistema de equações diferenciais ordinárias foram

comparados com dados experimentais. Os modelos previram de forma qualitativa as variações

da atividade enzimática em bactérias cultivadas em meio mínimo com e sem triptofano. Além

disso, as simulações reproduziram, qualitativamente, experimentos em mutantes da

E. coli.

Nas equações

(2.77) a seguir é apresentado um modelo composto por quatro

variáveis dependentes desenvolvido por Santillan e Mackey (2001a) o qual foi simplificado

fazendo-se, em relação ao modelo original,

= 0 e retirando-se o termo exp(-

μτ

) da

expressão da derivada da concentração da enzima antranilato sintase (E) em relação ao tempo,

dE/dt.

(

)

()()

()

()( )

()

[]

()()

()()()

()

tTK

GtE

tTK

tEtMK

tMKbebtOK

tOO

KtTK

KtT

CtT

μγτ

ηη

−

+−−=

+−+−=

⎟

⎠

⎞

⎜

⎝

⎛

−

max

(2.77)

Neste modelo, equações

(2.77), O

, M

, E e T são variáveis dependentes (todas em

função do tempo) e representam a concentração de operador livre, mRNA relativo ao gene

trpE, enzima antranilato sintase (AS) e triptofano, respectivamente.

Não serão mostrados detalhes dos significados de cada termo neste modelo. Para

tanto ver Santillan e Mackey (2001a; 2001b). Este modelo foi resolvido apenas para efeito de

comparação com a simulação estocástica.

3 METODOLOGIA

Neste capítulo serão descritas todas as abordagens usadas na implementação da

aplicação para simulação de vias metabólicas e reguladoras. Serão descritos, principalmente,

os recursos do M

ATLAB

que possibilitaram a implementação das técnicas descritas no

capítulo anterior.

3.1 MATLAB

O software MATLAB

® 2

é uma ferramenta comercial bastante usada em ciências

físicas e engenharias e está disponível em muitos centros de pesquisa. A principal vantagem

do M

ATLAB

sobre outras linguagens de programação – como C++, Java e outras – é que este

possui uma linguagem de script e ambiente de desenvolvimento interativo, o que permite ao

usuário editar, facilmente, as funções escritas nesta linguagem de acordo com as suas

necessidades (ULLAH et al., 2006). O M

ATLAB

tem outras funcionalidades, como: acesso

rápido às ferramentas de análise de dados e otimização; funções para visualização de dados

(amplo recurso para plotagem de gráficos) e recursos para se criar aplicativos contendo uma

interface gráfica. O M

ATLAB

possui, ainda, várias toolboxes

específicas para modelagem e

simulação em várias áreas das engenharias e ciências físicas. Além disso, o M

ATLAB

contém

funções para a representação e manipulação de vetores e matrizes multidimensionais. Há

versões do M

ATLAB

para vários sistemas operacionais, como MS-Windows, Linux, entre

outros. Com estas características é fácil escrever scripts relativamente complexos de forma

rápida e compacta.

Uma alternativa de domínio público e plataforma de software livre é o Octave

(www.octave.org). O Octave foi desenvolvido inicialmente pela comunidade de engenharia

química (University of Wisconsin) e está disponível tanto para MS-Windows quanto para

sistemas operacionais Unix-like. O Octave tem sintaxe muito parecida com a do M

ATLAB

embora seus recursos sejam um pouco mais limitados.

Uma das vantagens de se utilizar o M

ATLAB

, no entanto, é sua cultura amplamente

disseminada e a possibilidade de migração direta para aplicativos cliente-servidor.

MATLAB é marca registrada da MathWorks, http://www.mathworks.com/

Conjunto de rotinas MATLAB projetadas, com ou sem uma interface gráfica, para uma área específica; como

estatística, otimização, processamento de imagens, entre outras.

3.2 IMPLEMENTAÇÕES EM MATLAB

Nesta seção serão descritas as principais implementações, em linguagem MATLAB

dos métodos de análise de rede metabólica discutidos no capítulo

3.2.1 Problemas de Quadrados Mínimos

O problema de quadrados mínimos foi implementado usando a função qr do

ATLAB

. Esta função faz uma decomposição QR da matriz A. Com as matrizes Q e R foram

implementados os algoritmos descritos na seção

2.1.7. Através da matriz R determina-se o

posto de

A e, com isto, a classificação do sistema linear conforme a Figura 2. Comparações

do método usando a decomposição QR com os métodos usando a decomposição SVD e a

pseudo-inversa mostram que o primeiro é por volta de 5 vezes mais rápido que os outros dois.

3.2.2 Linguagens das Reações Bioquímicas

A LRB é uma linguagem formal e, como tal, pode ser processada pelo computador.

Para realizar este processamento foi usado o método conhecido como autômatos finitos

(SIPSER, 1996). Através do uso de autômatos finitos, foi possível analisar e verificar se uma

determinada reação está escrita corretamente e, além disso, produzir a matriz estequiométrica.

Foram criados e implementados 3 autômatos finitos:

 identificador;

 números;

 lado da reação (lado esquerdo ou direito do símbolo de reversibilidade).

Estes autômatos finitos são usados pelo módulo que contém a matriz estequiométrica

total, ou seja, o módulo totalStMatrix. Os autômatos finitos implementados processam a

LRBPP para a análise e geração da matriz

3.2.3 Análise de Fluxo Metabólico

Para a análise de fluxo metabólico foi reaproveitada, com algumas adaptações, uma

toolbox M

ATLAB

já implementada (MENDES, 2006). Nesta toolbox foram implementadas as

equações descritas em Stephanopoulos et al. (1998) para AFM, o que inclui:

 determinação dos fluxos não medidos a partir dos medidos;

 determinação dos novos fluxos medidos;

 análise da sensibilidade;

 reconciliação estatística.

A solução do sistema de equações lineares foi obtida usando-se os métodos de

solução descritos para o problema de quadrados mínimos.

3.2.4 Análise de Balanço de Fluxo

O problema de programação linear da análise de balanço de fluxo é resolvido usando

GLPKMEX, uma interface MEX do pacote glpk (GNU linear programming kit). Além disso,

foram usadas também as rotinas optimizeCbModel e solveCobraLP da toolbox COBRA (em

inglês: constraint-based reconstruction and analysis) (BECKER et al., 2007). A toolbox

COBRA é usada para a análise de fluxo baseada em restrição, como a ABF.

3.2.5 Modos Elementares

Para calcular todos os modos elementares foi usado o programa metatool, também

adaptado para o M

ATLAB

(VON KAMP e SCHUSTER, 2006). Este programa usa uma

linguagem semelhante a LRB, com menos flexibilidade, para a descrição das reações da rede

metabólica e armazena as informações da rede em uma estrutura diferente daquela usada pelo

GEnSys. Sendo assim, foi necessário, então, implementar uma interface para que fosse criada

uma estrutura que pudesse ser usada pelo metatool.

3.2.6 Relação de Conservação

A matriz de conservação (T) foi calculada como descrita na equação (2.43). O

primeiro passo é aumentar a matriz

N com a identidade I

(seção 2.4.3 - Relações de

Conservação) para formar a equação

(3.1):

[]

(3.1)

Reduzindo por linhas esta matriz à forma canônica obtém-se a matriz

T. Usando a

função rref do M

ATLAB

faz-se a redução através de eliminação de Gauss-Jordan. A função

rref tem a seguinte sintaxe:

[

R, jb] = rref(A) (3.2)

onde

R é a matriz reduzida e jb contém os índices das colunas de A referentes às variáveis

pivôs (variáveis independentes) da matriz

A. Logo, o tamanho do vetor jb representa o posto r

A e R(1:r, jb) é a matriz identidade r×r.

No trecho de código usado para calcular

T, que é apresentado no Quadro 2, N é a

matriz estequiométrica dos metabólitos internos.

Quadro 2: Cálculo da matriz T, em linguagem MATLAB

Uma vez calculada a matriz

T obtém-se a matriz L

e, desta última, a matriz L

. Os

metabólitos (ou fluxos) independentes e dependentes são todos determinados a partir dos

índices indepCol e depRow calculados conforme

Quadro 2.

3.2.7 Geração da Equação Cinética

Para cada uma das cinéticas, descritas na seção 4.2, é criada uma função em

ATLAB

que calcula a velocidade da reação, dados os valores dos metabólitos em um

determinado tempo.

Por exemplo, dada a reação v2:S2 <-> P2 (Michaelis-Menten reversível) e a sua

respectiva cinética, como no

Quadro 3.

%determina o número de linhas (m) e colunas (n) de N

[m,n] = size(N);

%calcula o valor de r, posto da matriz N

r = rank(N);

%aumenta a matriz N com uma matriz identidade de ordem m

Na=[N I]

Na = [N eye(m)];

%reduz Na à forma canônica, rb contém os índices das colunas

%dos pivôs

[echelonform, rb] = rref(Na);

%determina os índices das colunas independentes de N

indepCol = rb(1:r);

%determina os índices das linhas dependentes de N

depRow = rb(r+1:end)-n;

%determina, de echelonform, a matriz de conservação

T = echelonform(r+1:m,n+1:end);

Quadro 3: Exemplo de equação cinética.

É automaticamente gerado um arquivo .m (arquivo executável no M

ATLAB

)

contendo uma função equivalente à equação cinética do

Quadro 3.

Quadro 4: Função MATLAB

gerada automaticamente para a equação cinética do Quadro 3.

A primeira parte de função v2 é executada na primeira vez que v2 for chamada;

assim, a inicialização dos parâmetros é feita apenas uma única vez. A instrução F = @v2

atribui à F o handle da segunda função v2 (função aninhada) definida logo abaixo (Fe =

v2(X)); dessa forma, as próximas chamadas a v2 executarão apenas o código dentro da

segunda função v2.

Os nomes dos metabólitos são substituídos pelos valores de suas respectivas

concentrações, que são, por sua vez, passadas à função v2 através de seu argumento de

entrada

X. O vetor X contém as concentrações de todos os metabólitos da rede. Cada

metabólito possui um índice, e o vetor

X é indexado com este valor. No exemplo do Quadro

4, os índices de S2 e P2 são 10 e 4, respectivamente.

v2 % Michaelis-Menten reversível, equação (2.45)

Vmax = 2.0 % mM/min

KmS2 = 0.02 %mM

KmP2 = 1.0E-2 %mM

Keq = 19

v2 = Vmax*(S2 – P2/Keq)/(KmS2*(1 + P2/KmP2) + S2)

%Reaction rate for the v2 reaction

function F = v2(X)

% this block is executed only once

Vmax = 2.0;

KmS2 = 0.02;

KmP2 = 1.0E-2;

Keq = 19;

F = @v2;

%this function contains expressions dependent of metabolite

%concentrations

function Fe = v2 (X)

Fe = Vmax*(X(10) - X(4)/Keq)/(KmS2*(1 + X(4)/KmP2) +

X(10));

end

3.2.8 Simulação Determinística e Análise de Controle Metabólico

A simulação determinística é realizada pela resolução do conjunto de EDOs

definidas pela equação

(2.22). As velocidades v

são calculadas pelas cinéticas descritas vistas

anteriormente. Usou-se a função

ode15s.m, pré-definida no MATLAB

, para integrar

numericamente as EDOs em um intervalo de tempo definido pelo usuário. Esta função é um

solver baseado em fórmulas de diferenciação numérica. Geralmente, usa-se a

ode15s quando

ode45 falha ou é muito ineficiente ou ainda quando supeita-se que o problema é stiff.

Na determinação das elasticidades, equação

(2.50), foi usado o método de diferenças

finitas centrais para avaliar a derivada de v

(j = 1..n) em relação a cada concentração dos

metabólitos x

, conforme a equação (3.3):

()

(

)

xxvxxv

ijij

−

Δ+

(3.3)

Uma vez calculada a matriz ε que contém as elasticidades, os coeficientes de

controle podem ser calculados usando as equações

(2.64) a (2.68), considerando π, na

equação

(2.68), como sendo a matriz identidade.

4 RESULTADOS E DISCUSSÃO

Serão descritos neste capítulo os principais módulos e funções implementadas no

GEnSys. Para esta descrição serão mostrados vários exemplos de uso destes módulos. Alguns

destes exemplos são também resultados de simulações desenvolvidas neste trabalho.

Em vários pontos deste capítulo são usados comandos do M

ATLAB

diretamente no

texto. Os resultados destes comandos são exibidos no texto como saídas do próprio M

ATLAB

4.1 LINGUAGENS DAS REAÇÕES BIOQUÍMICAS

As reações devem ser escritas numa linguagem com sintaxe pré-definida e sem

ambigüidades, proposta aqui com o nome de Linguagem de Reações Bioquímicas (LRB).

Na LRB, cada reação é composta por:

 nome da reação seguida de : (dois pontos);

 lado esquerdo da reação (reagentes ou substratos);

 símbolo de reversibilidade da reação;

 lado direito da reação (produtos).

Esquematicamente, uma reação é escrita da seguinte forma:

[][] [][]

[

]

[

]

[

]

[

][][]

[

]

[

]

444444344444421

4434421

444444344444421

43421

(produtos) direito lado

2211

idadereversibil de símbolo

)(reagentes esquerdo lado

2211

reação da nome

*...**xnSímbolo*...**:

qqpp

yyyRxxxNomeRxn

onde os

são os coeficientes estequiométricos e x

e y

os compostos que participam das

reações. A notação entre colchetes significa que o símbolo ou o caractere é opcional; assim, *

e os coeficientes estequiométricos são opcionais. A palavra RxnSymbol representa um dos

seguintes símbolos: ->, >,=>, <->, = ou <=>.

 Nome da reação: identificador composto por letras, números ou sublinhados (_). Deve-se

começar com uma letra ou sublinhado. Exemplo: v1, r5, rGK e _r3. Após o nome da

reação deve aparecer sempre o símbolo :;

 Coeficiente estequiométrico: representa o número que precede um reagente (substrato)

ou um produto da reação. O coeficiente estequiométrico é um número inteiro ou real

(decimal ou notação científica). Exemplos: 1.0, 2, 1e-3, 2.5e-1 e 0.5;

 Reagentes e produtos da reação: representam os compostos (ou espécies, como às vezes

são chamados) envolvidos numa reação. Eles podem ser formados por letras; números; ‘

(apóstrofo); +; -; espaço; (; ) ou , (vírgula). O nome do composto pode começar com um

dígito, desde que seja sucedido de um ‘ (apóstrofo) ou um dos sinais + ou -. Quando

houver um sinal de + ou – formando um composto, este não pode ser precedido de um

espaço em branco (para não ser confundido com o sinal do coeficiente estequiométrico).

Exemplos: G6P; NAD+; F1,6biP; Rib-5-P; 3’ Gln; H2; H- e Xil 5 P;

 Sinal do coeficiente estequiométrico: os sinais dos coeficientes estequiométricos são os

símbolos + ou -. Estes devem sempre ser precedidos de um espaço em branco. Isto é

necessário pois estes símbolos também podem ser usados para formar os nomes dos

compostos – sem o espaço. Para o primeiro reagente (após o :) ou o primeiro produto

(após o sinal de reversibilidade, setas), estes sinais são opcionais;

 Símbolos de reversibilidade das reações: estes símbolos correspondem, em livros de

bioquímica, às setas que especificam a reversibilidade da reação. Pode-se usar os

seguintes símbolos, ->, > e => para as reações irreversíveis e <->, = e <=> para as reações

reversíveis.

A seguir, no

Quadro 5, apresenta-se um exemplo contendo um conjunto de reações

relacionadas à primeira parte da via glicolítica da Acetobacter xylinum:

Quadro 5: Exemplos de reações escritas em LRB.

r1:Glc(out) -> Glc(in) %transporte de glicose

r2:Glc(in) -> Glc-6-P %glicokinase

r3:Glc-6-P <-> Fru-6-P %fosfoglicoisomerase

r4:Fru1,6,diP -> Fru-6-P %frutose difosfatase

r5:Fru1,6,diP <-> DAcetone-P + GDH-3-P %aldolase

Pode-se usar, no arquivo de reações, os chamados comentários de linha. Isto quer

dizer que, após o comentário, que começa com %, não se pode escrever mais nenhuma reação

nesta linha (ver exemplo acima).

Foram definidos três tipos de reações:

 Reações usuais: são aquelas conforme definidas no Quadro 5, com os lados esquerdo e

direito completos.

 Reações de degradação: são reações do tipo

A -> null. Neste caso, o lado direito é a

palavra reservada null. Isto significa que o reagente será degradado numa determinada

velocidade com estequiometria α.

 Reações de geração: são reações do tipo null ->

A. Neste caso, o composto A é gerado

numa dada velocidade.

Estes dois últimos tipos de reações, geralmente, são usados para representar o

transporte de metabólitos para fora da célula (degradação) e para dentro da célula (geração),

respectivamente.

Uma vez definidas, as reações do modelo devem ser agrupadas em um arquivo texto

– arquivo LRB – de modo que elas possam ser lidas e, consequentemente, gerar a matriz

estequiométrica correspondente. Não há necessidade da utilização de um separador entre as

reações e pode haver várias delas em uma única linha do arquivo.

4.2 EXPRESSÕES CINÉTICAS ENZIMÁTICAS E VALORES INICIAIS

Para fazer uma simulação dinâmica da via metabólica é necessário definir as

expressões cinéticas de cada enzima (ver seção

2.5), bem como os valores iniciais para cada

metabólito interno. Deve-se seguir uma sintaxe própria para escrever a equação da cinética

para enzima da via. A sintaxe geral está descrita no

Quadro 6.

Quadro 6: Sintaxe de uma expressão cinética.

As seguintes definições são aplicadas:

 NomeRxn: corresponde ao nome da reação. Deve ser o mesmo nome dado à reação no

arquivo de reações.

 Inicializações das constantes: este trecho da cinética é reservado para a inicialização das

constantes, como os K

de cada metabólito (constante de Michaelis), as constantes de

inibição K

, a velocidade máxima da reação V

max

, entre outros. Este trecho é opcional, uma

vez que os valores das constantes podem ser usados diretamente na equação da cinética.

 dep: é uma palavra reservada e deve, sempre, estar precedida de : (dois pontos). Esta

palavra, que significa dependente, separa os valores cujo cálculo depende (depois de :

dep)

das concentrações dos metabólitos daqueles que independem, como as constantes

cinéticas.

 Valores dependentes das concentrações: este trecho é usado por cálculos intermediários.

Como, às vezes, a equação da cinética pode ser complexa, usa-se este trecho para dividir a

equação em partes menores.

 Expressão definindo a velocidade: esta é a expressão final que determina a velocidade

da reação em função da concentração atual dos metabólitos.

NomeRxn

Inicializações das constantes, se houverem (opcional).

Valores independentes das concentrações.

:dep

Determinação de valores dependentes das concentrações dos

metabólitos (opcional)

NomeRxn = Expressão definindo a velocidade

As cinéticas devem ser colocadas num arquivo texto que será usado para gerar, para

cada cinética, uma função no M

ATLAB

(arquivo .m). As expressões matemáticas seguem a

sintaxe definida no M

ATLAB

. Neste caso, também usa-se o símbolo % para fazer

comentários no arquivo.

Como exemplos, serão usadas as três equações apresentadas na seção

2.5. Dadas as

reações do

Quadro 7.

Quadro 7: Exemplo de reações da seção 2.5.

Para estas reações pode-se definir as cinéticas mostradas no

Quadro 8.

Quadro 8: Exemplos de equações cinéticas escritas numa sintaxe própria para as equações da

seção

2.5. Os valores dos parâmetros são aleatórios.

Os valores iniciais para cada metabólito interno também devem ser colocados em um

arquivo texto separado. A sintaxe deste arquivo é simples e é mostrada no

Quadro 9.

v1:S1 -> P1 % Michaelis-Menten simples

v2:S2 <-> P2 % Michaelis-Menten reversível

v3:P2 <-> P + Q % Mecanismo uni-bi

% Michaelis-Menten simples; equação

(2.44)

Vmax = 29.0 % mM/min

Km = 0.4 %mM

v1 = Vmax*S1/(S1 + Km)

%ou alternativamente

v1 = 29.0*S1/(S1 + 0.4)

% Michaelis-Menten reversível, equação (2.45)

Vmax = 2.0 % mM/min

KmS2 = 0.02 %mM

KmP2 = 1.0E-2 %mM

Keq = 19

v2 = Vmax*(S2 – P2/Keq)/(KmS2*(1 + P2/KmP2) + S2)

% Mecanismo uni-bi, equação (2.46)

Vmax1 = 58.0 % mM/min

Vmax2 = 21.0 % mM/min

KmP2 = 11.0E-2 %mM

KmP = 1.0 %mM

KmQ = 0.002 %mM

Keq = 9 %mM, neste caso Keq é dimensional

:dep

Num = Vmax1*Vmax2*(P2 – P*Q/Keq) %numerador da expressão

Den1 = Vmax2*KmP2 + Vmax2*P2 + Vmax1*P*KmQ/Keq +

Vmax1*KmP*Q/Keq

Den2 = Vmax2*P2*P/KmP + Vmax1*P*Q/Keq

v3 = Num/(Den1 + Den2)

%Den1 e Den2 fazem parte do denominador.

Quadro 9: Sintaxe para atribuir valores iniciais aos metabólitos internos.

onde Met

representa um dos m metabólitos internos e valor

o seu respectivo valor. Se o valor

de alguns dos metabólitos não for especificado, será atribuído, a estes metabólitos, o valor 0

(zero). No

Quadro 10 são mostradas as atribuições de valores iniciais considerando o exemplo

anterior.

Quadro 10: Exemplo de atribuição de valores iniciais.

Estes metabólitos devem aparecer em pelo menos uma reação. Se não aparecerem,

serão ignorados.

4.3 HIERARQUIA DOS MÓDULOS

Na Figura 14 é apresentada uma hierarquia dos módulos implementados no GEnSys

para a modelagem e simulação de redes metabólicas e reguladoras.

A principal característica desta organização em módulos é que eles podem ser

herdados por outros módulos e estes podem usar todos os dados e métodos do módulo base.

Dessa forma, há reutilização dos códigos do módulo base.

Met1 = valor1

Met2 = valor2

...

Metm = valorm

S1 = 2 %mM

P1 =0.5 %mM

S2 = 89 %mM

P2 = 6.9e-2 %mM

Figura 14: Hierarquia dos módulos que compõem o GEnSys. Os módulos de onde partem as

setas usam (herdam) as informações dos módulos aos quais estão ligados. Os três módulos

superiores, com o título em cinza, têm como função apenas coletar e armazenar as informações

sobre reações, metabólitos e a estequiometria da rede. Modelinfo: informações não

processadas sobre a rede; TotalStMat: estequiometria total; stMatrix: estequiometria interna;

FBA: análise de balanço de fluxo; MCA: análise de controle de fluxo; AFM: análise de fluxo

metabólico; DetSim: simulação dinâmica determinística; elemmode: modos elementares.

Cada módulo nesta hierarquia pode usar os métodos e variáveis do módulo no qual

está ligado. O módulo superior, básico, é o ModelInfo. Este módulo lê o arquivo que contém

as informações para a análise e simulação pretendida. Sendo assim, todos as análises

realizadas pelo GEnSys usam este módulo. A classe totalStMatrix produz a matriz

estequiométrica total (inclui todos os metabólitos da rede metabólica) e a classe stMat a

matriz estequiométrica interna. Os demais são:

 FBA: módulo para a análise de balanço de fluxo;

 MCA: módulo para a análise de controle metabólico;

 MFA: módulo para a análise de fluxo metabólico;

 DetSim: módulo para a simulação dinâmica determinística;

 Elemmode: módulo para análise e determinação dos modos elemetares.

Os outros 3 módulos são ferramentas (toolboxes) auxiliares. O metatool é um

programa que calcula todos os modos elementares dada uma rede metabólica ou reguladora.

O COBRA (constraint-based reconstruction and analysis) é um toolbox que faz análise de

balanço de fluxo estático e dinâmico.

Nas seções a seguir, cada um destes módulos, e seus usos, serão descritos com mais

detalhes.

4.3.1 Módulo ModelInfo

Este é um módulo no topo da hierarquia. Sua principal função é ler todas as

informações do modelo (armazenadas em um arquivo texto ou numa cadeia de caracteres),

fazer um pré-processamento e produzir uma estrutura contendo três campos:

 sRxn: é uma cadeia de caracteres que contém os dados das reações escritas em LRB. Tais

dados já passaram por um pré-processamento.

 sInf: cadeia de caracteres que contém dados não relacionados à equação das reações, isto

é, são outras informações não referentes à estequiometria, nomes de reações ou de

metabólitos. Elas podem ser usadas para alguma análise ou simulação do modelo. Estes

dados estão após as reações (nem todas as reações precisam conter tais dados) e devem

ser, obrigatoriamente, precedidos por um ponto e vírgula (;).

 sBefore: como sInf, é uma cadeia de caracteres que não faz parte das reações. Deve ser

usada para passar informações sobre o modelo. Estas informações devem preceder

qualquer reação do modelo.

As informações contidas em sInf dependem do objetivo do usuário (programador).

Elas podem conter, por exemplo, para uma reação: valor do fluxo (ou atividade) para a análise

de fluxo metabólico; dados cinéticos para simulação dinâmica determinística ou estocástica;

velocidade máxima e mínima para a análise de balanço de fluxo; entre outras. Já em sBefore,

pode ser qualquer informação – também dependente do objetivo da simulação – como por

exemplo, valores iniciais para os metabólitos em uma simulação dinâmica (estocástica ou

determinística).

Este módulo foi definido dessa forma para que possa ser usado na análise e

simulação – implementadas ou não neste trabalho – envolvendo redes metabólicas ou

reguladoras. Qualquer informação extra, não relacionada à estequiometria e nomes das

reações ou dos metabólitos, pode ser passada ao modelo através dos campos sBefore e sInf e

devem ser manipulados de acordo com o objetivo do usuário (programador). Assim definido,

este módulo confere reusabilidade do código, uma vez que ele pode ser facilmente estendido

para vários tipos de simulações in silico que envolvem redes metabólicas ou reguladoras. Nos

módulos que serão descritos a seguir, será visto como estas informações são extraídas e

usadas em vários tipos de simulação in silico de redes metabólicas.

4.3.2 Pré-processamento da LRB

A LRB é uma linguagem na qual o usuário deve escrever, de forma não ambígua,

todas as reações que compõe o modelo. O pré-processamento tem por objetivo preparar as

reações para que possam ser extraídos os nomes dos metabólitos e das reações e,

principalmente, a matriz estequiométrica. No pré-processamento da LRB são realizadas as

seguintes atividades:

 remoção de comentários, linhas e espaços em branco;

 adição de um delimitador das reações (o símbolo @) seguido, entre parênteses, do número

da linha, referente ao arquivo, onde a reação foi lida;

 substituição dos símbolos ->, => e > das reações irreversíveis por um único símbolo: >;

 substituição dos símbolos =, <-> e <=> das reações reversíveis por um único símbolo: =;

 substituição dos símbolos ‘ +’ (espaço seguido de +) e ‘ -’ (espaço seguido de -) por ? e !,

respectivamente. Isto é realizado inclusive onde o sinal ‘ +’ (espaço seguido de +) está

implícito: depois de : ou depois do sinal de reversibilidade da reação. Por exemplo, na

reação DCB A:r1 −=>+ , o símbolo ‘ +’ está implícito antes dos metabólitos A e C.

Na verdade, o pré-processamento traduz a LRB para outra linguagem equivalente

(LRB pré-processada, LRBPP), porém mais simples e mais adequada para a análise sintática

do conjunto de reações – realizada por um conjunto de autômatos finitos – que compõem o

modelo.

Além das atividades descritas acima, outra tarefa do pré-processamento é a

separação dos dados em sRxn (string contendo as reações em LRBPP), sBefore e sInf (estas

duas últimas strings contêm outras informações relacionadas ao modelo).

Por exemplo, dado o conjunto de reações gravado no arquivo exemplo.txt mostrado

Quadro 11:

Quadro 11: Exemplo de reações e outras informações que formam um modelo gravado num

arquivo.

A=0 B=0.3 D = 1e-4

r1:A + B - > C

r2:B + C < -> D; 58

r3:D + E = > E + F; 47

r4:F -> G + D; 789

r5:E + G <=> A

Também, estes dados podem ser atribuídos a uma string strExemplo no M

ATLAB

strExemplo = 'A=0 B=0.3 D = 1e-4 r1:A + B - > C r2:B + C < -> D; 58 r3:D +

E = > E + F; 47 r4:F -> G + D; 789 r5:E + G <=> A'

strExemplo =

A=0 B=0.3 D = 1e-4 r1:A + B - > C r2:B + C < -> D; 58 r3:D + E = > E + F;

47 r4:F -> G + D; 789 r5:E + G <=> A

O número após o ; para algumas das reações representa o fluxo no regime

estacionário, por exemplo. Após o pré-processamento, que descarta os espaços em branco e

comentários, será criada uma estrutura contendo os dois campos sRxn

e sInf. Para criar o

objeto model a partir das reações acima, basta digitar no M

ATLAB

model= modelinfo ('exemplo.txt') ou

model= modelinfo (strExemplo , 's')

model =

sBefore: 'A=0 B=0.3 D = 1e-4 '

sInf: [1x5 struct]

sRxn:

'@(1)r1:?A?B>?C@(2)r2:?B?C=?D@(3)r3:?D?E>?E?F@(4)r4:?F>?G?D@(5)r5:?E?G=?A'

No campo sRxn cada reação é separada pelo delimitador @ seguido, entre

parênteses, pelo número da linha do arquivo onde esta reação foi lida. Este número é usado

para que, caso haja algum erro de sintaxe (falta de um : ou do sinal de reversibilidade, por

exemplo), seja emitido uma mensagem de erro informando o possível erro e o número da

linha onde se encontra a reação no arquivo. No caso de o argumento da função modelinfo ser

uma string de reações, o número da linha corresponde a ordem, na string, da reação. As

demais informações correspondem àquelas substituições descritas anteriormente. Isto é, os

sinais ‘ +’ e ‘ –‘ foram substituídos por ? e !, respectivamente; cada sinal de uma reação

reversível ou irreversível foi substituído por = ou >, respectivamente.

Para o campo sInf é criado um array de estrutura (structure array) que contém, em

cada posição deste array, uma estrutura com os campos: FileLine e RxnData que representam,

respectivamente, a linha do arquivo ou ordem da reação na string (com o mesmo propósito

como descrito anteriormente) e a informação relacionada àquela reação. No exemplo anterior,

sInf: [1x5 struct] é um array com cinco estruturas, uma para cada reação. O conteúdo de sInf

é mostrada na

Tabela 1:

Tabela 1: Exemplo das informações armazenadas no campo sInf da estrutura modelinf.

Índice: 1 2 3 4 5

FileLine: -1 3 4 5 -1

RxnData: '' ' 58' ' 47' '789' ''

Onde não há dados para uma reação, é atribuído -1 e

'' para os campos FileLine e

RxnData, respectivamente.

Os campos sBefore e sInf não podem possuir o caractere @ (meta-símbolo que

delimita as reações) e nem : (dois pontos) precedido de uma letra (ou o sublinhado, _) seguido

de letras ou sublinhados ou dígitos para que não sejam interpretados como nome de reação

seguido de : (dois pontos).

A partir dos dados em sRxn serão produzidos: uma lista com os nomes das reações,

uma lista com os nomes dos metabólitos e a matriz estequiométrica. Ressalta-se que estes

últimos dados são “constantes” em todas as classes onde são implementados os demais

módulos descritos neste trabalho; por exemplo, a matriz estequiométrica é necessária em

todas as classes referentes a um módulo de simulação (AFM, ABF, simulação dinâmica, etc).

Por outro lado, as informações em sBefore e

sInf variam conforme o módulo de

simulação implementado. Por exemplo, na AFM e na atribuição de fluxos, sInf contém os

valores dos fluxos medidos; na simulação dinâmica determinística, estas informações são

nomes de equações cinéticas para as reações; sBefore pode conter, como visto no exemplo, os

valores iniciais na simulação dinâmica e na ABF elas representam o intervalo que um fluxo

pode variar (valores mínimos e máximos). Portanto, em cada uma destas classes deve haver

um método específico para extrair e interpretar tais informações. No entanto, para as

informações em sRxn

– que são úteis em todas as classes e é constante – é necessário apenas

um método que possa extrair, adequadamente, tais informações. Este é o propósito do

próximo módulo a ser descrito.

4.3.3 Módulo totalStMat

Este módulo tem a tarefa de processar o campo sRxn e obtém, como resultado, os

seguintes campos:

 stMatrix

()

qp × : matriz estequiométrica total (todos os metabólitos) contendo

metabólitos eq reações;

 ReacMatrix

()

: matriz estequiométrica contendo apenas os reagentes de todas as

reações. Sendo assim seus coeficientes são todos negativos;

 ProdMatrix

()

: matriz estequiométrica contendo apenas os produtos de todas as

reações. Sendo assim seus coeficientes são todos positivos;

 lComp: listas dos nomes dos

compostos das reações (reagentes e produtos);

 lRxn: listas dos nomes das q reações;

 Irrev_Rxn: irreversibilidades das q reações (1: irreversível, 0: reversível);

 indExt: índices dos metabólitos externos, entre 1 e

, em relação às linhas de stMatrix;

 indInt: índices dos metabólitos internos, entre 1 e

,em relação às linhas de stMatrix;

O construtor desta classe cria uma instância da classe modelinfo, e desta, obtém-se o

campo sRxn de onde são extraídos todos os campos descritos acima. A matriz

stMatrix é

determinada pela soma das matrizes dos regentes (

ReacMatrix) e dos produtos

(

ReacMatrix). Estas duas últimas matrizes serão necessárias na simulação dinâmica

estocástica.

Usando o exemplo anterior, deve-se seguir a seguinte sintaxe no M

ATLAB

model = totalStMat('exemplo.txt') ou

model= totalStMat(strExemplo , 's')

model =

stMatrix: [7x5 double]

ReacMatrix: [7x5 double]

ProdMatrix: [7x5 double]

lComp: {'A' 'B' 'C' 'D' 'E' 'F' 'G'}

lRxn: {'r1' 'r2' 'r3' 'r4' 'r5'}

Irrev_Rxn: [1 0 1 1 0]

indExt: 5

indInt: [1 2 3 4 6 7]

A seguir a matriz stMatrix é exibida:

model.stMatrix

ans =

-1 0 0 0 1

-1 -1 0 0 0

1 -1 0 0 0

0 1 -1 1 0

0 0 0 0 -1

0 0 1 -1 0

0 0 0 1 -1

Observa-se pelo exemplo que há sete metabólitos (7 linhas na matriz

stMatrix),

sendo que um deles foi considerado externo (linha 5, metabólito E) e seis internos (linhas 1, 2,

3, 4, 5 e 7; correspondem aos metabólitos, respectivamente, A, B, C, D, F e G). Há também

cinco reações (r1, r2, r3, r4 e r5) sendo que r1, r3 e r4 são irreversíveis (1) e as demais

reversíveis (0).

O processamento das reações foi implementado usando-se autômatos finitos.

4.3.4 Metabólitos Internos e Externos

É importante identificar os metabólitos considerados internos ou externos numa rede

metabólica. No regime estacionário os metabólitos internos têm concentração constante. Isto

quer dizer que as velocidades de produção e consumo destes metabólitos devem ser iguais.

Nem sempre os metabólitos externos precisam ser extracelulares.

Para ser considerado externo, um metabólito deve satisfazer os seguintes critérios

(

Figura 15):

 aparecer apenas no lado esquerdo das reações irreversíveis da qual participa (entrada);

 aparecer apenas no lado direito das reações irreversíveis da qual participa (saída);

 aparecer apenas uma vez no lado esquerdo ou direito de uma reação reversível.

Nos demais casos os metabólitos são considerado internos.

Figura 15: Esquema gráfico para a classificação dos metabólitos como internos ou externos.

Nos casos (a), (b) e (c) o metabólito A é considerado como externo e no caso (d) como interno.

No entanto, em algumas análises é necessário fazer com que um metabólito torne-se

externo ou interno. No GEnSys adiciona-se o sufixo

_x ou _i ao nome de um metabólito para

que seja considerado externo ou interno, respectivamente. Por exemplo, adicionando os

sufixos

_i ou _x em alguns metabólitos do exemplo anterior:

r1:A_x + B - > C

r2:B + C < -> D

r3:D_x + E_i = > E_i + F

r4:F -> G + D

r5:E_i + G <=> A_x

Quadro 12: Exemplo de reações com metabólitos forçados a serem internos ou externos.

AA A A

(a) (b)

(d)

(c)

tem-se o seguinte resultado:

stMatrix: [8x5 double]

ReacMatrix: [8x5 double]

ProdMatrix: [8x5 double]

lComp: {'A_x' 'B' 'C' 'D' 'D_x' 'E_i' 'F' 'G'}

lRxn: {'r1' 'r2' 'r3' 'r4' 'r5'}

Irrev_Rxn: [1 0 1 1 0]

indExt: [1 5]

indInt: [2 3 4 6 7 8]

Neste exemplo, o metabólito E_i tornou-se externo e os metabólitos A_x e D_x se

tornaram externos. Comparar com o exemplo anterior.

4.3.5 Módulo stMat

Este é o módulo que armazena as informações do modelo relacionadas apenas aos

metabólitos internos. Como foi visto, os metabólitos de modelo são classificados como

externos e internos. Os externos, para a maioria das análises, são removidos da rede

metabólica. Sendo assim, ao se remover tais metabólitos, as respectivas linhas da matriz

estequiométrica total correspondente a estes metabólitos também são removidas, obtendo,

então, uma nova matriz estequiométrica. Além disso, os metabólitos externos são removidos

da lista de metabólitos (lComp) . Esta classe possui os seguintes campos:

 st

()

nm ×

: matriz estequiométrica interna; contém m linhas correspondentes aos

metabólitos internos e

n colunas correspondentes às reações;

 lComp

()

m×1

: listas dos nomes dos m compostos (metabólitos) internos;

 lRxn

()

n×1 : listas dos nomes das reações. Mesma lista da classe totalStMat;

 Irrev_Rxn: vetor de tamanho n (número de reações na rede). Este vetor indica se a j-

ésima reação é irreversível (1) ou reversível (0);

 MoetyCycle

()

mrm ×− )( : é uma matriz onde cada linha representa um ciclo na rede

(relação de conservação) entre os metabólitos. O número de colunas desta matriz é igual

ao número de metabólitos internos (m) e o número de linhas é igual ao número de ciclos

na rede (m – r). Os elementos da i-ésima linha que possuem o valor diferente de 0 fazem

parte do i-ésimo ciclo. Há tantos destes ciclos quanto há linhas linearmente dependentes

na matriz

st, ver seção 2.4.3;

 depRow

()

)(1 rm

−

× : vetor linha que contém os índices, em relação à matriz st, das linhas

linearmente dependentes. r é o posto da matriz

st;

 indepRow

()

r×1 : contém os índices, em relação à matriz st, das linhas linearmente

independente, r é o posto da matriz

st;

 indepCol

()

)(1 rn −× : contém os índices, em relação à matriz st, das colunas linearmente

independentes,

r é o posto da matriz st;

 depCol

()

r×1

: contém os índices, em relação à matriz st, das colunas linearmente

dependentes.

onde:

r é o posto da matriz st, m é o número de metabólitos (número de linhas de st) e n o

número de reações (número de colunas de

st). Será usado como exemplo o modelo da rede da

próxima seção para mostrar os conteúdos dos campos descritos acima.

4.4 SIMULAÇÃO DINÂMICA COM O GENSYS

Para se criar um modelo dinâmico, antes de tudo, é necessário definir todas as

reações que o compõem. Isto é feito estudando-se a via de interesse – objeto de estudo do

modelo – através de consultas à literatura e aos bancos de dados disponíveis na

Internet para

coletar as informações estruturais (estequiometria) e cinéticas (simulação dinâmica)

relacionadas às reações da via a ser simulada.

4.4.1 Descrição do Exemplo para Simulação Dinâmica

Será usada como exemplo a rede metabólica descrita na Figura 16. Esta rede

metabólica representa a parte superior da glicólise da

Saccharomyces cerevisiae (HYNNE et

al.

, 2001).

Figura 16: Via metabólica representando a primeira parte da glicólise da Saccharomyces

cerevisiae.

Glicose

Glic-6-P Fruc-6-P Fruc-1,6-P2

ATP

ADP

ATP

ADP

ATP

ADP

ATP

ADP

ATP + AMP

2ADP

O primeiro passo, uma vez escolhida a rede objeto de estudo, é criar as reações em

LRB (

Quadro 13).

Quadro 13: Reações em LRB da rede mostrada na Figura 16.

v1:Glicose + ATP -> ADP + Glic-6-P

v2:Glic-6-P + ATP -> ADP

v3:Glic-6-P <-> Fruc-6-P

v4:Fruc-6-P + ATP -> Fruc-1,6-P2 + ADP

v5:Fruc-1,6-P2 -> null

v6:ADP -> ATP

v7:ATP -> ADP

v8:ATP + AMP_i <-> 2ADP

Na reação v5, o lado direito,

null, indica que a reação não forma nenhum produto.

Sendo assim, ela é considerada uma reação de saída. Estas reações foram salvas no arquivo

Reactions.m. O seguinte comando, no MATLAB

, produzirá um objeto útil para a simulação

determinística.

model = DetSim('Reactions.m')

model =

totalModel: [1x1 struct]

internalModel: [1x1 struct]

lRxn: {'v1' 'v2' 'v3' 'v4' 'v5' 'v6' 'v7' 'v8'}

hJ: {8x1 cell}

O campo hj, de model, armazena um handle para cada uma das reações em lRxn. Um

handle é um recurso do MATLAB

para guardar o nome de uma função que será executada

pelo próprio M

ATLAB

4.4.2 Matriz Estequiométrica Total

O comando acima produziu a estrutura model, que possui, como campos, outras

estruturas. A estrutura

totalModel representa o modelo total, isto é, informações relacionadas

a matriz estequiométrica total (metabólitos externos e internos):

model.totalModel

ans =

stMatrix: [7x8 double]

lComp: {'Glicose' 'ATP' 'ADP' 'Glic-6-P' 'Fruc-6-P' 'Fruc-1,6-P2'

'AMP_i'}

indExt: 1

indInt: [2 3 4 5 6 7]

O comando a seguir exibe a matriz estequiométrica total (matriz

S na seção 2.2):

model.totalModel.stMatrix

ans =

-1 0 0 0 0 0 0 0

-1 -1 0 -1 0 1 -1 -1

1 1 0 1 0 -1 1 2

1 -1 -1 0 0 0 0 0

0 0 1 -1 0 0 0 0

0 0 0 1 -1 0 0 0

0 0 0 0 0 0 0 -1

Cada coluna j desta matriz representa uma reação; a ordem das reações é a mesma

dada em

model.lRxn. Para cada metabólito i é criado uma linha nesta matriz. Os coeficientes

com sinais negativos significam que o metabólito

i é um reagente da j-esima reação e aqueles

com sinais positivos são produtos.

4.4.3 Matriz Estequiométrica Interna

A matriz estequiométrica interna, ou simplesmente matriz estequiométrica, é obtido

da matriz total excluindo-se desta as linhas correspondentes aos metabólitos externos.

Também, deve-se excluir estes metabólitos da lista de metabólitos. Assim, obtém-se:

model.internalModel

ans =

st: [6x8 double]

lComp: {'ATP' 'ADP' 'Gluc-6-P' 'Fruc-6-P' 'Fruc-1,6-P2'

'AMP_i'}

Irrev_Rxn: [1 1 0 1 1 1 1 0]

MoetyCycle: [1 1 0 0 0 1]

depRow: 1

indepRow: [2 3 4 5 6]

indepCol: [1 2 3 4 8]

depCol: [5 6 7]

Onde: st representa a matriz estequiométrica interna:

model.internalModel.st

ans =

-1 -1 0 -1 0 1 -1 -1

1 1 0 1 0 -1 1 2

1 -1 -1 0 0 0 0 0

0 0 1 -1 0 0 0 0

0 0 0 1 -1 0 0 0

0 0 0 0 0 0 0 -1

Neste caso há um ciclo, isto é, a matriz MoetyCycle contém uma linha. Este ciclo é

formado pelos metabólitos ATP, ADP e AMP_i.

4.4.4 Conjunto de Equações Diferenciais Ordinárias

Uma vez definido o modelo, o seguinte comando exibe as equações diferenciais

ordinárias (ODE, na sigla em inglês) do mesmo:

printODE(model)

ODE system:

d(ADP)/dt = v1 + v2 + v4 - v6 + v7 + 2v8

d(Gluc-6-P)/dt = v1 - v2 - v3

d(Fruc-6-P)/dt = v3 - v4

d(Fruc-1,6-P2)/dt = v4 - v5

d(AMP_i)/dt = -v8

Algebric equations (conservation relationship):

const1 = ATP + ADP + AMP_i

Como se vê, é definida uma EDO para cada metabólito interno da rede. Este

conjunto de EDOs é um exemplo das equações

(2.21) e (2.22). A relação de conservação

(

conservation relationship; que foi originada da matriz MoetyCycle descrita anteriormente) dá

origem a uma equação algébrica. Isto significa que a soma destes metabólitos é constante ao

longo do tempo. Neste caso, o que se conserva, no ciclo de reações, é adenosina (A) em cada

um dos metábolitos envolvidos neste ciclo. O valor de

const1 é determinado a partir dos

valores iniciais dos metabólitos internos ATP, ADP e AMP_i.

4.4.5 Cinéticas das Reações

Para realizar a simulação dinâmica é necessário definir as cinéticas de cada uma das

reações da rede metabólica,

v da equação (2.23). Para este exemplo, será criado o arquivo

Kinetics.m contento as seguintes cinéticas (usando o comando type para exibir o conteúdo do

arquivo):

type 'Kinetics.m'

%kinetic expressions

Vmax = 1398

KmATP = 0.1

KmGlucose = 0.37

v1 = (Vmax*Glucose*ATP)/(1 + ATP/KmATP + Glucose/KmGlucose +

(ATP*Glucose)/(KmATP*KmGlucose))

k2 = 2.26; %mM^-1*min^-1

v2 = k2*ATP*Gluc-6-P

Vf_max = 140.282

Vr_max = 140.282

KmGluc6P = 0.8

KmFruc6P = 0.15

v3 = ((Vf_max*Gluc-6-P/KmGluc6P) - (Vr_max*Fruc-6-P/KmFruc6P)) / (1 + Gluc-

6-P/KmGluc6P + Fruc-6-P/KmFruc6P)

Vmax = 44.7287

KmFruc6P = 0.021; %mM^2

k = 0.15

v4 = Vmax*Fruc-6-P^2/(KmFruc6P*(1 + k*(ATP/AMP_i)^2) + Fruc-6-P^2)

k5 = 6.04662

v5 = k5*Fruc-1,6-P2

k6 = 68.48

v6 = k6*ADP

k7 = 3.21

v7 = k7*ATP

k8f = 432.9

k8r = 133.33

v8 = k8f*ATP*AMP_i - k8r*ADP^2

Neste arquivo, como pode se ver acima, é definido os valores das constantes

cinéticas e a equação da cinética (

reaction law) de cada uma das reações dadas no arquivo

Reactions.m. A sintaxe para a criação do arquivo de cinéticas foi definida na seção 4.2.

Através deste arquivo, o seguinte comando gera as cinéticas – como funções executáveis no

ATLAB

– para cada uma das reações do modelo model:

model = setKinetic(model, 'Kinetics.m')

model =

totalModel: [1x1 struct]

internalModel: [1x1 struct]

lRxn: {'v1' 'v2' 'v3' 'v4' 'v5' 'v6' 'v7' 'v8'}

hJ: {8x1 cell}

Este comando gera, a partir das informações em Kinetics.m, as funções em MATLAB

com os mesmos nomes das reações em

model.lRxn.

4.4.6 Valores Iniciais dos Metabólitos

Por último, é necessário definir os valores iniciais de todos os metabólitos – de todos

os internos e daqueles externos que aparecem numa expressão da cinética –, uma vez que será

resolvido um sistema de equações diferenciais ordinárias. Estes são, também, definidos em

um arquivo separado

InitVal.m:

type 'InitVal.m'

%initial values

Glucose = 12.8174

Gluc-6-P = 1

ATP = 2.01

ADP = 1.4

AMP_i = 0.1

Serão inicializados com zero todos aqueles metabólitos que não estão explicitamente

inicializados no arquivo acima.

Através do comando

initValFromFile todos os valores iniciais são definidos:

model = initValFromFile(model, 'InitVal.m')

model =

totalModel: [1x1 struct]

internalModel: [1x1 struct]

lRxn: {'v1' 'v2' 'v3' 'v4' 'v5' 'v6' 'v7' 'v8'}

hJ: {8x1 cell}

X0: [7x1 double]

LX: [-1.0000 0 0 0 -1]

A: 3.5100

Foram criados dois novos campos em

model: X0 (vetor coluna com m linhas que

contém os valores iniciais de todos os metabólitos) e

LX (matriz ligadora de metabólitos),

definida na seção

2.4.3.

O valor da constante

const1 agora pode ser determinado. Executando o comando

printODE de novo, este valor aparece valendo 3.51:

printODE(model)

ODE system:

d(ADP)/dt = v1 + v2 + v4 - v6 + v7 + 2v8

d(Gluc-6-P)/dt = v1 - v2 - v3

d(Fruc-6-P)/dt = v3 - v4

d(Fruc-1,6-P2)/dt = v4 - v5

d(AMP_i)/dt = -v8

Algebric equations (conservation relationship):

3.51 = ATP + ADP + AMP_i

4.4.7 Forma Simplificada

Existe uma forma mais fácil e simples para a construção de um modelo para

simulação dinâmica determinística. Foi criado no GEnSys algumas funções cinéticas pré

definidas. Estas funções – ou seus nomes – podem ser usadas diretamete no arquivo de

reações. Na

Tabela 2 são mostrados os mecanismos e a sintaxe das cinéticas pré definidas no

GEnSys.

Tabela 2: Cinéticas pré definidas no GEnSys. Quando há dois substratos, estes podem se ligar

à enzima de forma ordenada (um depois o outro) ou aleatória. Por isso, alguns mecanismos

com dois substratos podem ser classificados como ordenado ou aleatório. Além disso, se há

um ou dois substratos, a enzima pode ser classificada como uni ou bi, respectivamente.

Símbolos, Vmax: velocidade máxima da reação; kms1, kms2, kmp1 e kmp2 são as constantes

de Michaelis dos substratos 1 e 2 e dos produtos 1 e 2. Recomenda-se Leskovac (2003) para

mais detalhes.

Sintaxe das cinéticas Mecanismo e observações

hill(Vmax, Kms1, n)

Equação de Hill. Cinética com um único substrato.

Usado onde há cooperatividade.

massaction(kf, kr)

Ação das massas. Cinética reversível. Não há um

número definido de substratos e produtos

mm(Vmax, Kms1)

Michaelis-Menten irreversível. Cinética com um

único substrato.

mwc(Vmax, Kms1, L, n, c)

Monod, Wyman e Changeux. Cinética com um

único substrato. Usado onde há cooperatividade.

orderedbibi(Vmax, Kms1,

Kms2, Kmp1, Kmp2, Keq)

Dois substratos e dois produtos ordenados. Cinética

reversível.

orderedbisub(Vmax, Kms1,

Kms2)

Dois substratos ordenados. Cinética irreversível.

orderedunibi(Vmax, Kms1,

Kmp1, Kmp2, Keq)

Um substrato e dois produtos ordenados. Cinética

reversível.

pingpongbibi(Vmax, Kms1,

Kms2, Kmp1, Kmp2, Keq)

Ping pong bisubstrato. Dois substratos e dois

produtos e reversível

randbisub(Vmax, Kms1,

Kms2)

Bisubstrato aleatório. Irrevesível com dois

substratos.

rmm(Vmax, kms1, kmp1, keq)

Michaelis-Menten reversível. Um substrato e um

produto.

theorellchance(Vmax, Kms1,

Kms2, Kmp1, Kmp2, Keq)

Mecanismo theorellchance. Reversível com 2

substratos e dois produtos

Além das cinéticas, pode-se – neste mesmo arquivo – inicializar todos os

metabólitos. Assim, para o exemplo visto anteriormente, considerando a inicialização e as

cinéticas pré definidas, o arquivo do modelo tem a seguinte forma:

type 'Reactions2.m'

%example: subset of reactions of glycolysis in yeast

%initial values

Glucose = 12.8174

Gluc-6-P = 1

ATP = 2.01

ADP = 1.4

AMP_i = 0.1

%reactions

v1:Glucose + ATP -> ADP + Gluc-6-P

v2:Gluc-6-P + ATP -> ADP + null; massaction(2.26)

v3:Gluc-6-P <-> Fruc-6-P; rmm(140.282, 0.8, 0.15, 0.1875)

v4:Fruc-6-P + ATP -> Fruc-1,6-P2 + ADP

v5:Fruc-1,6-P2 -> null; massaction(6.04662)

v6:ADP -> ATP; massaction(68.48)

v7:ATP -> ADP; massaction(3.21)

v8:ATP + AMP_i <-> 2ADP; massaction(432.9, 133.33)

Observando o arquivo acima, vê-se que para as reações v1 e v4 não foi definida

cinética alguma. Para as demais foram usadas as cinéticas pré-definidas

massaction (lei de

ação das massas) e

rmm (Michaelis-Menten reversível). As cinéticas para as reações v1 e v4

estão definidas no arquivo

Kinetics2.m; este agora menor que o original:

type 'Kinetics2.m'

%kinetic expressions

Vmax = 1398

KmATP = 0.1

KmGlucose = 0.37

v1 = (Vmax*Glucose*ATP)/(1 + ATP/KmATP + Glucose/KmGlucose +

(ATP*Glucose)/(KmATP*KmGlucose))

%/(KmATP*KmGlucose)

Vmax = 44.7287

KmFruc6P = 0.021; %mM^2

k = 0.15

v4 = Vmax*Fruc-6-P^2/(KmFruc6P*(1 + k*(ATP/AMP_i)^2) + Fruc-6-P^2)

É possível, portanto, criar um modelo dinâmico em apenas um arquivo que conterá

as reações, os valores iniciais e as cinéticas. O usuário, se desejar, pode adicionar novas

cinéticas pré-definidas. Para isso, basta criar um arquivo – extensão .m – com a seguinte

forma:

nomePar_1 = par1

nomePar_2 = par2

...

nomePar_N = parN

v = f(par1, par2,…, parN, S1, S2,…, SM, P1, P2,…, PL)

Esta é a sintaxe genérica para a criar uma cinética com

N parâmetros, M substratos e

L produtos. Nesta sintaxe, f(.) significa que a expressão matemática é em função daqueles

parâmetros e variáveis. Por exemplo, a cinética pré-definida de uma reação reversível com

dois substratos e dois produtos usando o mecanismo conhecido como

pingpongbibi é escrita

como:

type pingpongbibi.m

%ping pong bi bi reaction law

%syntax: pingpongbibi(Vmax, Kms1, Kms2, Kmp1, Kmp2, Keq)

Vmax = par1

Kms1 = par2

Kms2 = par3

Kmp1 = par4

Kmp2 = par5

Keq = par6

v = Vmax*(S1*S2 - P1*P2/Keq)/(Kms2*S1 + Kms1*S2 + Kms1*Kms2*P1/Kmp1 +

Kms1*Kms2*P2/Kmp2 + Kms2*S1*P1/Kmp1 + Kms1*Kms2*P1*P2/(Kmp2*Kmp1) +

Kms1*S2*P2/Kmp2 + S1*S2)

4.4.8 Simulação Dinâmica do Modelo

Uma vez definido o modelo – pelos passos descritos acima –, o comando [t, X] =

transient(model, tspan), com tspan = [t

], determina as concentrações dos metabólitos no

intervalo de tempo

tspan integrando o sistema de equações diferenciais de t

a t

com as

condicões iniciais definidas anteriormente. Cada linha na matriz solução

X corresponde aos

tempos retornados no vetor coluna

t e cada coluna de X representa a concentração dos

metabólitos internos. Os comandos abaixo calcularão as concentrações dos metabólitos, na

matriz

X, em cada passo de tempo t, no intervalo tspan. Estes valores não foram exibidos.

tspan = [0 0.5];[t X] = transient(model, tspan);

As concentrações dos metabólitos são melhores visualizadas através de gráficos. O

comando

subplotTrans(model, nGraph, X, t) exibe, em única figura, os gráficos (transientes)

das concentrações de todos os metabólitos internos (

X) ao longo do tempo (t). Cada gráfico

2D relaciona a concentração (eixo y) de cada um dos metabólitos versus o tempo (eixo x).

nGraph representa o número de gráficos na horizontal.

O comando abaixo exibe os gráficos das concentrações dos metabólitos, calculadas

anteriormente, em uma única figura (

Figura 17):

subplotTrans(model,2,X,t);

0 0.1 0.2 0.3 0.4 0.5

1.5

2.5

time (min)

ATP (mM)

0 0.1 0.2 0.3 0.4 0.5

1.5

time (min)

ADP (mM)

0 0.1 0.2 0.3 0.4 0.5

time (min)

Gluc-6-P (mM)

0 0.1 0.2 0.3 0.4 0.5

0.5

time (min)

Fruc-6-P (mM)

0 0.1 0.2 0.3 0.4 0.5

time (min)

Fruc-1,6-P2 (mM)

0 0.1 0.2 0.3 0.4 0.5

0.5

time (min)

AMP

(mM)

Figura 17: Concentrações dos metabólitos internos ao longo do tempo. Estes gráficos foram

produzidos pelo comando subplotTrans.

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

0.5

1.5

2.5

3.5

4.5

time (min)

ATP

ADP

Gluc-6-P

Fruc-6-P

Fruc-1,6-P2

AMP

Figura 18: Gráficos de todas as concentrações dos metabólitos internos num única figura.

Produzido pelo comando plotTransFig.

O comando

plotTrans(model, X, t) exibe os m (número de metabólitos internos)

gráficos 2D, como acima, em figuras separadas (não mostradas).

Já o comando

plotTransFig(model, X, t) cria uma figura contendo todos os gráficos

dos metabólitos sobrepostos (

Figura 18):

plotTransFig(model,X,t)

Também é possível calcular e exibir os valores dos fluxos (atividade enzimática) de

cada reação ao longo do tempo. Os comandos a seguir, como os anteriores, calculam e exibem

um gráfico com os valores dos fluxos ao longo do tempo:

Fluxes = allRxnRates(model, X);subplotFluxes(model, 2, Fluxes, t)

0 0.1 0.2 0.3 0.4 0.5

time (min)

v1 (mM/time)

0 0.1 0.2 0.3 0.4 0.5

time (min)

v2 (mM/time)

0 0.1 0.2 0.3 0.4 0.5

100

time (min)

v3 (mM/time)

0 0.1 0.2 0.3 0.4 0.5

time (min)

v4 (mM/time)

0 0.1 0.2 0.3 0.4 0.5

time (min)

v5 (mM/time)

0 0.1 0.2 0.3 0.4 0.5

100

150

time (min)

v6 (mM/time)

0 0.1 0.2 0.3 0.4 0.5

time (min)

v7 (mM/time)

0 0.1 0.2 0.3 0.4 0.5

-200

200

time (min)

v8 (mM/time)

Figura 19: Valores dos fluxos de cada reação (enzima) da rede metabólica no tempo. Estes

gráficos foram obtidos pelo comando subplotFluxes que exibe os fluxos calculados pela

função AllRxnRates.

Este modelo simula um biorreator funcionando em modo contínuo – a concentração

da glicose (único substrato de entrada) é constante –, por isso faz sentido calcular o valor das

concentrações no regime estacionário:

xss = steadystate(model)

xss =

1.6494

1.4617

4.0252

0.4434

5.3578

0.3990

4.5 SIMULAÇÃO ESTOCÁSTICA DO OPERON TRP DA E. COLI

Como visto anteriormente, há duas abordagens para a simulação dinâmica:

estocástica e determinística. Na abordagem estocástica, um conjunto de reações em um meio

espacialmente homogêneo e com volume constante representando o sistema é simulado pelo

algoritmo de simulação estocástica (SSA) (ver seção

2.3.2 para saber mais detalhes sobre este

algoritmo). Embora a abordagem determinística seja adequada na maioria dos casos, ela não

reflete a natureza estocástica do sistema que é importante em muitos sistemas biológicos,

como naqueles onde há alguma espécie com um número pequeno de cópias numa célula

(DNA ou moléculas reguladoras, por exemplo) (MCADAMS e ARKIN, 1997).

Nesta seção serão apresentadas todas as reações que envolvem a síntese do triptofano

(Trp) e a regulação do

operon trp na E. coli. Será levada em consideração, como no modelo

da equação

(2.77), apenas a subunidade TrpE da enzima antranilato sintase (AS). Isto se

justifica pelo fato desta enzima ser a primeira enzima da via e ainda ser inibida (através da

subunidade TrpE) pelo produto final, o triptofano.

Um dos principais mecanismos de regulação do

operon trp na E. coli é a repressão.

operon trp existem 5 genes estruturais (trpE-trpD-trpC-trpB-trpA, ou trpEDCFA) os quais

codificam as 5 enzimas que catalisam a síntese de Trp na

E. coli a partir do corismato, o

principal precursor da via de biossíntese do Trp. Este

operon também contém um promotor,

um operador (repressão) e uma região líder (atenuação da transcrição). Os genes

trpEDCBA

são transcritos para um mRNA policistrônico que, por sua vez, são traduzidos para suas

respectivas enzimas.

A repressão ocorre da seguinte forma: uma molécula do repressor ativo liga-se ao

operador que bloqueia a ligação da RNA polimerase (RNAP) ao promotor impedindo, assim,

o início da transcrição. O repressor (um dímero formado por duas subunidades codificadas

pelo gene

trpR) deve estar ativado pelo Trp para se ligar ao operador. Esta ativação se dá pela

ligação não cooperativa de 1 molécula de Trp a cada uma das duas subunidades do TrpR. Esta

ativação é representada pelas reações:

RTR

⎯⎯→←+

−+

RTR

⎯⎯→←+

−+

(4.1)

onde,

R, T, R

, R

, K

e K

representam, respectivamente, repressor inativo, triptofano,

complexo repressor_ativo–1Trp, repressor ativo, constante de velocidade pra frente e

constante de velocidade reversa.

A repressão se dá quando um repressor ativo liga-se ao operador numa reação

reversível:

ORO

⎯⎯→←+

−+

(4.2)

Neste caso

é o complexo operador-repressor, O

é o operon livre. A transcrição é

realizada quando uma RNAP se liga ao

numa reação irreversível com velocidade k

MOO

+⎯→⎯

(4.3)

onde M representa o mRNA com tamanho suficiente para conter o sítio de ligação do

ribossomo (RBS). Pelas reações

(4.2) e (4.3) percebe-se que a RNAP concorre com o R

para

se ligar ao

operon livre.

O mRNA

M pode ter sua transcrição interrompida por outro mecanismo de

regulação, a atenuação da transcrição. Após a RNAP transcrever o RBS, o ribossomo se liga

ao mRNA

M e inicia-se a tradução. Se o nível intracelular de Trp estiver alto, a tradução deste

mRNA promove a formação de um grampo que induz a RNAP a interromper a transcrição.

Por outro lado, se houver pouco Trp no meio intracelular, durante a tradução que é realizada

simultaneamente com a transcrição, o ribossomo faz com que outra estrutura seja formada de

forma a não interromper a transcrição. Tem-se a seguinte reação representando a atenuação da

transcrição:

⎯→⎯

(4.4)

Onde

é o mRNA livre que será traduzido para formar a enzima AS. A

, neste

caso, não é a constante de velocidade de transcrição e, sim, a função

()

(

)

ebTA

−

−−= 11,

onde

representa a probabilidade de haver atenuação da transcrição, uma formulação

adequada para simulação estocástica.

O ribossomo ligado ao mRNA livre (

) que não sofreu atenuação dará origem ao

complexo ribossoma–mRNA (

), cuja reação catalítica irreversível é:

FRF

MMM

+⎯⎯→⎯

(4.5)

Nesta reação,

é dividido por dois porque a AS é formada por duas subunidades.

Assim, o nível de AS será a metade do número de subunidade traduzida por molécula de

A partir do complexo

, a enzima AS, representada por E, é alongada pelo

ribossomo através da reação irreversível:

()

EtM

trad

eFR

⎯⎯→⎯−

(4.6)

onde

é o atraso de tempo devido à tradução. A tradução leva

minutos para alongar o

peptídeo

TrpE, por isso o nível de E é proporcional ao nível de M

minutos atrás. Neste

caso, deve-se considerar

em função do tempo t.

A inibição alostérica da AS é o último mecanismo de regulação considerado na

biossíntese do Trp. Duas moléculas de Trp ligam-se de forma cooperativa a cada uma das

subunidades TrpE da AS. As duas reações reversíveis que representam a inibição enzimática

são:

ETE

⎯⎯→←+

−+

ETE

⎯⎯→←+

−+

(4.7)

Finalmente, tem-se a reação da produção de

T a partir da atividade enzimática E. A

produção do triptofano

T é dada pela seguinte reação:

ETE

+⎯→←

−

(4.8)

e o consumo de Trp pela bactéria para a síntese de proteína

pro

⎯⎯→⎯

(4.9)

onde

(

)

gpro

KTPP

é a probabilidade de acontecer esta reação. De novo, foi usado uma

equação em função do nível de triptofano (

T) para denotar a probabilidade de acontecer tal

reação.

Alguns dos compostos presentes nas reações acima são degradados e/ou diluídos. As

reações a seguir representam este fato

Φ⎯→⎯Φ⎯→⎯Φ⎯→⎯Φ⎯→⎯Φ⎯⎯→⎯Φ⎯⎯→⎯

μμμμ

μμ

T TE ET E M M

,2,,,,

(4.10)

onde μ é a velocidade de crescimento específico da bactéria, e k

a velocidade de degradação

do mRNA. O símbolo

Φ, presente no lado direito da reação, representa a degradação do

composto da esquerda com velocidade

. Onde se tem apenas

, como constante

de velocidade da reação, há apenas diluição. Onde há

há degradação com velocidade k

e diluição do composto no meio intracelular.

As reações de

(4.1) a (4.10) acima representam a dinâmica de biossíntese do Trp

levando-se em conta todas as regulações existentes no

operon trp da E. coli.

4.5.1 Parâmetros

Todos os parâmetros para a simulação determinística e estocástica estão na Tabela 3.

A maioria dos parâmetros usados na simulação foi estimada por Santillan e Mackey (2001b) e

estão na unidade de concentração

μM ou μM

-1

min

-1

. Assim, estes valores devem ser

convertidos para número de moléculas (ver

2.3.2). Por isso, para todos os valores iniciais dos

metábolitos e das constantes de velocidades foi usado o seguinte fator de conversão:

−

= para k

em unidades

min

−−

;

VCNN

amolecule

−

= para C em unidades

(4.11)

onde

é o número de Avogrado, k

e C são as constantes a serem convertidas. Usando o

valor para

V = 8,0E-16 litros e N

= 6,02214199E+23, todos os parâmetros foram convertidos

e são mostrados na

Tabela 3.

O parâmetro k

trad

foi estimado considerando que a velocidade de tradução do mRNA

para proteína é de 15 amino ácidos por segundo. Sabendo-se que o TrpE tem 520 amino

ácidos, estimou-se que

trad

é 0,0288 s

-1

(1.73 min

-1

) para cada ribossomo que traduz o mRNA

obtido da transcrição do gene

trpE. Para probabilidade

pro

P da reação (4.9), usou-se P =

1,0E6.

P representa a quantidade de triptofano presente em todas as proteínas da célula no regime

estacionário

é a velocidade específica de crescimento.

Tabela 3: Constantes usadas nas simulações determinística e estocástica.

Nome da

constante

Valor simulação

determinística

Valor

simulação

estocástica

Nome da

constante

Valor

simulação

determinística

Valor simulação

estocástica

460 μM

-1

min

-1

0,96

0,8μM

385

1,2 min

-1

1,2

3,32E-3μM

176 μM

-1

min

-1

0,37

0,2μM

720 min

-1

720

2,0E-2μM

348 μM

-1

min

-1

0,72 OF

1,31E-004μM

2,1E4 min

-1

2,1E4 MF

3,27E-4μM

nu*

10,14min

-1

10,14

0,32μM

nu*

20,01min

-1

20,01

4,1μM

1975

0,6 min

-1

0,6

0,85 0,85

1,0E-2 min

-1

1,0E-2

1,2

nu*

149,44 min

-1

149,44

0,1 min

nu*

25 μM min

-1

1,2E4

0,66 min

0,66

trad

1,73

1,0E6

8E-16 L

*não usado.

determinado por Santillan e Mackey (2001a).

4.5.2 Resultados

O algoritmo usado para a simulação está descrito na seção 2.3.2. Foi usada a vesão

descrita em Gibson e Bruck (2000) e implementada no programa DIZZI (RAMSEY

et al.,

2005).

Com estes parâmetros foram realizadas várias simulações imitando um experimento de

derepression

. Isto é, uma cultura de E. coli de uma linhagem selvagem foi cultivada em meio

mínimo com triptofano adicionado ao meio. Após algum tempo, as células são lavadas e

centrifugadas e depois incubadas no mesmo meio, porém, sem triptofano. A partir deste

momento, é medida a atividade enzimática da AS ao longo do tempo. Para simular este

experimento, resolveu-se o modelo acima com a condição inicial apresentada na

Tabela 4.

Esta condição representa o início da segunda parte do experimento, onde as bactérias – por

terem estado no meio com excesso de triptofano – não possuem mRNA nem enzimas ativas

da via do triptofano. Como em Santillan e Mackey

(2001b), para efeito de comparação, os

valores no regime estacionário de ambos resultados (experimento e simulação) foram

forçados a terem o mesmo valor. A

Figura 20 mostra a comparação dos resultados obtidos

experimentalmente e através da simulação usando os dados da

Tabela 3. Pode observar por

esta figura que o resultado de ambas as simulações prediz qualitativamente a dinâmica da

atividade enzimática.

Tempo (min)

Atividade da AS normalizada

Figura 20: Comparação da atividade enzimática obtida pelos modelos determinístico (linha

pontilhada) e estocástico (linha cheia) com os dados experimentais (círculos e pontos). Estes

valores estão normalizados, pois os dados experimentais estão em unidades arbritárias. Os

valores da simulação estocástica foram obtidos pela média de 22 simulações.

Tabela 4: Valores iniciais para os compostos na simulação estocástica. Estes valores refletem

um experimento de derepression.

= 2

= 383

T =

29625

M = 0

= 0

E =

= 0

Pela simulação estocástica, na média há 25 moléculas de AS ativa na célula no

regime estacionário (dados não mostrados). Isto dá uma atividade final igual a 0.13

μM de

triptofano formado por segundo. Pela

Figura 20 pode-se ver que nas simulações estocástica e

determinísica, a atividade enzimática leva por volta de 6 minutos para ser maior que zero, ou

seja, para começar a sintetizar o Trp. Este período de tempo é porque há bastante trp no meio

e consequentemente não há síntese da enzima AS ou atividade enzimática da AS até que o

triptofano caia até um valor mínimo, neste caso por volta de 1100 molécula na solução

estocástica, veja

Figura 21. Quando o Trp chega a este valor mínimo, os mecanismos de

regulação recrudessem e a síntese de mRNA relativo à AS aumenta e, consequentemente, a

síntese da enzima AS.

Tempo (min)

Moléculas de triptofano

Figura 21: Número de moléculas de triptofano ao longo do tempo relacionado ao experimento

de derepression. Este gráfico mostra os dados da média de 22 simulações. O Trp decresce

rapidamente até um valor mínimo, depois ele aumenta até atingir o regime estacionário.

4.6 ANÁLISE DE CONTROLE METABÓLICO NA PRODUÇÃO DE PENICILINA

Como ilustração, serão calculados os coeficientes da ACM para uma pequena rede

metabólica. Na

Figura 22 é exibida a via metabólica com os primeiros passos na produção de

penicilina pela

Penicillium chrysogenum. Theilgaard e Nielsen (1999) construíram um

modelo cinético para esta via e calcularam uma expressão matemática dos coeficientes de

controle de fluxo para cada reação.

Figura 22: Os dois primeiros passos na via de biossíntese de penicilina pela P. chrysogenum.

As linhas tracejadas indicam inibição da reação pelo metabólito onde se origina a linha. Nesta

via, os metabólitos Isopenicilina N e bisACV inibem a enzima ACVS. Abreviações dos

metabólitos, LLD-ACV: (

L-aminoadipil)-L-cisteinil-D-valina; bisACV: bis-(L-aminoadipil)-L-

cysteinil-

D-valina). Enzimas, J1: ACV Sintetase (ACVS); J2: Isopenicilina N sintetase (IPNS);

J3 (TRS_OX) e J4 (TRS_RED): sistema da tioredoxina redutase.

Nesta via metabólica são incluídos a redução e a formação da bisACV a partir do

primeiro intermediário da via LLD-ACV e a inibição da enzima ACVS (J

) pelos metabólitos

bisACV e isopenicilina N. O modelo cinético é baseado na cinética do tipo Michaelis-Menten

com inibição não competitiva da ACVS por ambos metabólitos.

Para se fazer ACM no GEnSys, isto é, calcular os coeficientes de controle, o

primeiro passo é construir o modelo dinâmico determinístico da via (seção

4.4):

m = DetSim('Rxn.m'); m = setKinetic(m, 'kin.m')

m =

totalModel: [1x1 struct]

internalModel: [1x1 struct]

lRxn: {'rACV_J1' 'rIPN_J2' 'rOX_J3' 'rRED_J4'}

X0: [8x1 double]

LX: [0x2 double]

A: [0x1 double]

hJ: {4x1 cell}

Neste caso, como visto anteriormente, m é uma estrutura que representa o modelo

dinâmico e é criado a partir dos arquivos de reações

Rxn.m e de cinéticas kin.m.

Cisteína

Valina

Ácido Aminoadipídico

3ATP

3AMP

LLD-

Isopenicilina N

NADPH

NADP+

0,5O

bisACV

Uma vez construído o modelo dinâmico, usa-se a função

contrlCoef para calcular o

três coeficientes: elasticidade, controle de fluxo e controle de concentração. Esta função tem a

seguinte sintaxe:

[eps, fcc, ccc] = contrlCoef(m)

onde: eps, fcc e ccc são os coeficientes de elasticidade, de controle de fluxo e de controle de

concentração, respectivamente;

m é o modelo dinâmico. Portanto, para determinar todos os

coeficientes, basta executar o comando no M

ATLAB

[Eps, fcc, ccc] = contrlCoef(m)

Eps =

0.0737 0

0.8228 0

2.0000 0

0 0.9982

fcc =

1.0984 -0.0984 -0.0000 0.0000

2.6697 -2.6697 1.0000 0.0000

2.6684 -2.6684 0.9995 0.0000

ccc =

1.3349 -1.3349 -0.0000 0.0000

2.6733 -2.6733 1.0013 -1.0019

Para efeito de visualização, os FCCs deste exemplo estão reproduzidos na

Tabela 5.

Cada coluna

representa a variação da atividade enzimática da ésimaj

−

enzima e cada linha

ésimok −

fluxo (

J ).

Tabela 5: Coeficiente de controle de fluxo para as quatro enzimas da rede mostrada na Figura

22.

1 2 3 4

1,0984 -0,0984 0,0000 0,0000

2,6697 -2,6697 1,0000 0,0000

2,6684 -2,6684 0,9995 0,0000

Estes coeficientes só podem ser comparados entre si, isto é, não podem ser

comparados com coeficientes de outras vias. Uma pequena pertubação na atividade da enzima

ACVS (j = 1) influencia os fluxos J

e J

com maior intensidade que os outros dois fluxos

(coluna 1 da

Tabela 5). Pela coluna 4, percebe-se que a enzima TRS_RED (j = 4) não

influencia (controla) qualquer fluxo nesta rede.

4.7 ANÁLISE DE FLUXO METABÓLICO

Na análise de fluxo metabólico, os fluxos não medidos são calculados, usando a

equação

(2.70), a partir dos fluxos medidos experimentalmente; estes, geralmente, são os

fluxos dos metabólitos externos. Para criar um módelo para a análise de fluxo metabólico

deve-se seguir sintaxe mostrada no

Quadro 14 (estes dados devem ser gravados num arquivo).

Quadro 14: Sintaxe para a criação de um modelo para a AFM.

r_p

= v_1

r_p

= v_2

...

r_p

= v_p

%n reações envolvendo m metabólitos, em LRB

r_1: reação 1 em LRB

r_2: reação 2 em LRB

...

r_n: reação N em LRB

No arquivo acima,

, com qi ,,2,1 K

, é um subconjunto de

{}

n,,2,1 K onde q é o

número de fluxos medidos e r_

é o nome destes fluxos (reações). Além disso, as n reações

(r_1, r_2, ..., r_n) da rede devem ser fornecidas e escritas em LRB.

Como exemplo será usada a rede mostrada na

Figura 23, página 82, e os dados

experimentais estão na

Tabela 11, página 96; ambos se encontram na seção 4.8. O modelo de

AFM, com os dados experimentais e as reações, é mostrado abaixo (arquivo H2MFA_EX.m):

%fluxos medidos

v1 = 3.94

v2 = 4.65

v15 = 3.86

v19 = 0.42

v20 = 0.66

v6 = 7.24

v10 = 1.3

%reações em LRB

v1: Glicose + 2ATP -> 2ADP_x + 2GDH-3-P

v2: Glicerol + NAD+_x + ATP -> NADH + ADP_x + GDH-3-P

v3: GDH-3-P + 2ADP_x + NAD+_x -> PIR + 2ATP + NADH

v4: PIR + Fdoxi + CoA -> Acetil-CoA + Fdred

v5: Fdoxi + NADH <-> Fdred + NAD+_x

v6: Fdred -> Fdoxi + H2

v7: PIR + NADH <-> NAD+_x + Lactato

v8: 2Acetil-CoA -> CoA + AcAcetil-CoA

v9: Acetil-CoA + ADP_x <-> ATP + CoA + Acetato

v10: Acetil-CoA + 2NADH -> CoA + 2NAD+_x + Etanol

v11: Acetato + AcAcetil-CoA -> Acetil-CoA + Acetone

v12: Butirato + AcAcetil-CoA -> Butiril-CoA + Acetone

v13: AcAcetil-CoA + 2NADH -> Butiril-CoA + 2NAD+_x

v14: Butiril-CoA + ADP_x <-> CoA + ATP + Butirato

v15: Butiril-CoA + 2NADH -> CoA + 2NAD+_x + Butanol

v16: 0.6667Glicose + 0.5833NADH + 9.9ATP -> Biomassa

v17: ATP -> ADP_x

v18: Lactato -> LactatoExt

v19:Acetato -> AcetatoExt

v20:Butirato -> ButiratoExt

O comando a seguir cria uma estrutura representando um módulo de AFM:

modelMFA = mfa('H2MFA_EX.m')

modelMFA =

st: [13x20 double]

lRxn: {1x20 cell}

lComp: {1x13 cell}

vMeasFlux: [7x1 double]

iMeasFlux: [1 2 6 10 15 19 20]

iCalcFlux: [3 4 5 7 8 9 11 12 13 14 16 17 18]

A função mfa retornará a estrutura contendo os seguintes campos:

 st: matriz estequiométrica dos metabólitos internos;

 lRxn e lComp são, como já descritas, as listas dos nomes das reações (fluxos) e dos

nomes dos metabólitos, respectivamente.

 vMeasFlux: vetor coluna contendo os valores dos q fluxos medidos;

 iMeasFlux: vetor linha contendo os índices dos q fluxos medidos;

 iCalcFlux: vetor linha contendo os índices dos n – q fluxos calculados;

No exemplo acima existem 7 fluxos medidos. Como há 20 fluxos, então 13 serão

calculados (não medidos). Para calcular os fluxos não medidos é necessário apenas executar a

função:

solvemfa(modelMFA, pcs). Os fluxos calculados são mostrados na Tabela 6 e a norma

do resíduo é de 3.77E–014.

Tabela 6: Resultado dos fluxos calculados usando os dados da Tabela 11 como fluxos

medidos.

12,53 11,37 -4,13 1,16 4,82 0,63 0,21 0,09 4,52 0,75 1,34 0,61 1,16

Para mais detalhes dos outros valores calculados, recomenda-se Mendes (2006).

Cabe lembrar que, na AFM, todas as reações são consideradas reversíveis; isto significa que,

mesmo escrevendo uma reação como irreversível, o fluxo calculado para esta reação pode dar

um valor negativo.

4.8 MODELAGEM E SIMULAÇÃO DA PRODUÇÃO DE BIOHIDROGÊNIO

Nesta seção será descrita e estudada uma rede metabólica da bactéria Clostridium

acetobutylicum

. Os estudos desta rede serão efetuados através da análise de balanço de fluxo e

por modos elementares, ambos com vista a produção de hidrogênio (H

). Esta rede será

aproveitada como estudo de caso para descrever, no GEnSys, os dois métodos mencionados

anteriormente.

4.8.1 Rede Metabólica da C. acetobutylicum

A rede metabólica para este organismo foi construída a partir de dados da literatura

(SHINTO

et al., 2007; JONES et al., 1986) e possui 20 reações e 25 metabólitos (11 externos

e 14 internos), ver

Figura 23. Foram considerados dois substratos como fonte de carbono:

glicose e glicerol.

Figura 23: Rede bioquímica utilizada pela C. acetobutylicum para a conversão de glicose e

glicerol em hidrogênio, ácidos (acetato e butirato) e solventes (butanol e acetona) em

condições anaeróbicas. O hidrogênio é produzido pela reação v

. São 20 reações e 25

metabólitos; destes, 11 são considerados metabólitos internos e 14 externos. Há duas reações

auxiliares, v

e v

que correspondem, respectivamente, a velocidade de crescimento celular e

a ATPase. O número 2 entre parênteses nas reações v

e v

representa a estequiometria em

relação ao metabólito e a reação próxima a ele. As reações v

e v

são catalisadas pela mesma

enzima: CoA transferase.

Glicose

GlicerolGDH-3-P

NADH NAD+

ATPADP

2ATP

2ADP

(2)

PIR

Acetil-CoA

FdOx

FdRed

NAD+

NADH

CoA

Lactato

NADHNAD+

AcAcetil-CoA

(2)

CoA

2Glicose

1,75NADH

29,7ATP

ADP

ATP

Acetato

ADP

ATP

CoA

Etanol

2NAD+

2NADH

Acetona

Butirato

Butyril-CoA

2NADH

2NAD+

Butanol

2NAD+2NADHATP ADP

CoA

LactatoExt

AcetatoExt

ButiratoExt

CoA

3Biomassa

Muitas reações nesta rede foram agrupadas como uma única reação. Estes grupos

formam um conjunto de reações lineares – sem ramificação – que, no regime estacionário,

possuem os mesmos fluxos; assim, reduz-se o tamanho da rede sem perda de informação.

Como exemplo, a glicólise não tem todas as reações explícitas nesta rede.

Os principais produtos considerados são os solventes (butanol, acetona e etanol),

os ácidos (acetato, butirato e lactato) e o hidrogênio molecular H

. Além destes, a biomassa

também foi considerada como saída e modelada como uma reação. As 20 reações foram

nomeadas por

2021

v,,v,v K e seus significados estão na Tabela 7. As reações v

e v

são,

respectivamente, de entrada da glicose e do glicerol.

Tabela 7: Nomes das enzimas (ou fluxos) da rede metabólica da Figura 23. Algumas reações,

sequência linear, foram agrupadas formando uma única reação. É o caso da glicólise, por

exemplo.

Nome Nome da enzima (fluxo)

Velocidade da glicólise

Glicerol desidrogenase

Velocidade de produçao de piruvato

Piruvato-ferredoxina oxidoredutase

NADH-ferredoxina oxidoredutase;

Hidrogenase

Velocidade de produção de lactato

Tiolase

Fosfotransacetilase

Acetaldeído desidrogenase

CoA transferease 1

CoA transferease 2

Butiril-CoA desidrogenase

Fosfotransbutirilase

Butiraldeído desidrogenase

Velocidade de crescimento celular

ATPase

Transporte de lactato

Transporte de acetato

Transporte de butirato

Os nomes e significados dos 25 metabólitos encontram-se na

Tabela 8. Ali não estão

listados os metabólitos cujos nomes estão completos, como etanol, butirato entre outros. Em

algumas análises realizadas neste trabalho, os cofatores (NAD+, NADH, ATP e ADP) e a

CoA foram considerados externos. Isto significa que eles estão disponíveis em quantidades

necessárias para ocorrer a reação na qual participam como substratos.

Tabela 8: Significados dos nomes dos metabólitos da rede da Figura 23.

Metabólito Significado

ATP Adenosina trifosfato

ADP Adenosina difosfato

GDH-3-P Gliciraldeído-3-fosfato

NAD+ Nicotinamida adenina dinucleotídeo

NADH Nicotinamida adenina dinucleotídeo, forma reduzida

PIR Piruvato

FdOxi Ferredoxina oxidada

CoA Coenzima A

FdRed Ferredoxina reduzida

H2 Hidrogênio molecular

AcAcetil-CoA Aceto acetil coenzima A

LactatoExt Lactato externo

AcetatoExt Acetato externo

ButiratoExt Butirato externo

A velocidade de crescimento celular para a

C. acetobutilycum foi a mesma usada por

Papoutsakis (1984):

3Biomassa 29,7ATP 1,75NADH 2Glicose :v

→

De acordo com o mesmo autor, a glicólise nesta bactéria, durante a fase de

crescimento exponencial e síntese de acetato, produz uma quantidade excessiva de ATP. Isto

implica que um sistema ativo de ATPase deve operar durante o crescimento

(PAPOUTSAKIS, 1984). Sendo assim, a seguinte reação foi adicionada ao modelo:

ADP ATP :v

→ .

4.8.2 Modos Elementares da Rede Metabólica para a Produção de Biohidrogênio

Nesta seção será apresentado o módulo sobre modos elementares. Para tanto, será

usada a rede metabólica da

Figura 23 da Clostridium acetobutylicum como de estudo de caso.

Esta bactéria é muito conhecida e estudada como produtora de solventes e hidrogênio.

Para se criar um objeto deste módulo, após ter definido a rede metabólica, entra-se

com o comando (as reações, escritas em LRB, estão gravadas no arquivo

H2GlucoseGlicerol.m):

model = elemmode('H2GlucoseGlicerol.m')

model =

OverallSt: [14x29 double]

intModel: [1x1 struct]

extModel: [1x1 struct]

metatoolStruct: [1x1 struct]

ems: [20x29 double]

Uma estrutura é criada com os seguintes campos:

 OverallSt: matriz estequiométrica q

das reações globais, onde p é o número de

metabólitos externos e

q é o número de modos elementares. Estas reações envolvem

apenas os metabólitos externos.

 intModel: modelo dos metabólitos internos, visto na seção 4.3.5;

 extModel: modelo dos metabólitos externos, visto na seção 4.3.3;

 metatoolStruct: é a estrutura criada pela toolbox metatool (VON KAMP e SCHUSTER,

2006);

 ems: matriz

qn ×

dos modos elementares, onde n é o número de reações da rede e q o

número de modos elementares. Cada coluna representa um modo elementar e cada

elemento da coluna é a estequiometria da reação que faz parte daquele modo.

Tabela 9 mostra os 29 modos elementares calculados para a rede metabólica da

Figura 23. A coluna da esquerda exibe o número do modo e a coluna da direita todos os

modos elementares com as suas respectivas reações e estequiometria. Todos os modos são

irreversíveis (um modo irreversível é aquele que possui, no mínimo, uma reação irreversível,

caso contrário é considarado reversível). Os modos 8 e 9 são compostos de apenas uma

reação, pois estas reações têm como substratos e produtos apenas compostos externos. Estes

modos são obtidos através do comando

dispElmo(model).

Tabela 9: Os 29 modos elementares correspondentes à rede da Figura 23. Todos os modos são

irreversíveis. Os coeficientes negativos que precedem algumas reações significam que o

sentido do modo é o contrário daquele no qual a reação foi escrita. Em negrito encontram-se os

modos de produção de hidrogênio. ME: modo elementar.

ME Reações de cada modo elementar

{v1} {2*v3} {2*v4} {2*v6} {2*v10}

{v5} {v6}

3 {v2} {v3} {v7} {v18}

4 {v2} {v3} {v4} {-v5} {v10}

5 {2*v2} {2*v3} {2*v4} {-2*v5} {v8} {v9} {v11}

6 {2*v2} {2*v3} {2*v4} {-2*v5} {v8} {v12} {v14}

7 {2*v2} {2*v3} {2*v4} {-2*v5} {v8} {v13} {v15}

8 {v16}

9 {v17}

10 {v2} {v3} {v4} {-v5} {v9} {v19}

11 {2*v2} {2*v3} {2*v4} {-2*v5} {v8} {v13} {v14}

{v20}

12 {v1} {2*v3} {2*v7} {2*v18}

13 {v1} {2*v3} {2*v4} {-2*v5} {2*v10}

14 {v1} {2*v3} {2*v4} {-2*v5} {2*v9} {2*v19}

15 {v1} {2*v3} {2*v4} {-2*v5} {v8} {v9} {v11}

16 {v1} {2*v3} {2*v4} {-2*v5} {v8} {v12} {v14}

17 {v1} {2*v3} {2*v4} {-2*v5} {v8} {v13} {v15}

18 {v1} {2*v3} {2*v4} {-2*v5} {v8} {v13} {v14} {v20}

{v2} {v3} {v4} {v6} {v10}

{v1} {2*v3} {2*v4} {2*v6} {2*v9} {2*v19}

{v2} {v3} {v4} {v6} {v9} {v19}

{v1} {2*v3} {2*v4} {2*v6} {v8} {v9} {v11}

{v1} {2*v3} {2*v4} {2*v6} {v8} {v12} {v14}

{v1} {2*v3} {2*v4} {2*v6} {v8} {v13} {v15}

{2*v2} {2*v3} {2*v4} {2*v6} {v8} {v9} {v11}

{2*v2} {2*v3} {2*v4} {2*v6} {v8} {v12} {v14}

{2*v2} {2*v3} {2*v4} {2*v6} {v8} {v13} {v15}

{v1} {2*v3} {2*v4} {2*v6} {v8} {v13} {v14} {v20}

{2*v2} {2*v3} {2*v4} {2*v6} {v8} {v13} {v14} {v20}

Como visto na

Tabela 9, os modos 1, 2 e de 19 a 29 são as possíveis forma de

produção de H

. O modo 2 é apenas um ciclo que não está ligado, de alguma forma, a um dos

2 substratos externos. Na

Figura 24 são apresentados, em forma gráfica, os modos

elementares de produção de H

a partir de glicerol. No total são 13 modos de produção de H

um a partir do NADH, seis a partir da glicose e seis do glicerol. A diferença entre estes dois

últimos grupos é apenas a troca da reação v

pela v

Figura 24: Esquema gráfico dos modos elementares de produção de H

a partir do glicerol,

considerando a

Figura 23. ME: modo elementar.

O comando

dispOverallRxn(model) exibe todas as reações globais obtidas a partir

dos modos elementares. Os substratos e produtos destas reações são formados apenas pelos

metabólitos externos. A matriz estequiométrica destas reações é obtida multiplicando-se a

matriz dos metabólitos externos pela dos modos elementares. A

Tabela 10 exibe todas as

reações globais obtidas pelo comando

dispOverallRxn(model).

1 2

5 6

ME 1

ME 20

1 2

5 6

ME 22

1 2

5 6

ME 23

ME 24

1 2

5 6

ME 28

Tabela 10: Reações globais relacionadas aos 29 modos elementares. O sufixo _x em alguns

metabólitos significa que eles foram tratados como metabólitos externos.

# Reções globais

Glicose + 2*ADP_x + 2*NADH_x

→ 2*ATP_x + 2*NAD+_x + 2*H2 + 2*Etanol

NADH_x

→ NAD+_x + H2

ADP_x + Glicerol + NAD+_x

→ ATP_x + NADH_x + LactatoExt

ADP_x + Glicerol + NAD+_x

→ ATP_x + NADH_x + Etanol

3*ADP_x + 2*Glicerol + 6*NAD+_x

→ 3*ATP_x + 6*NADH_x + Acetone

3*ADP_x + 2*Glicerol + 6*NAD+_x

→ 3*ATP_x + 6*NADH_x + Acetone

2*ADP_x + 2*Glicerol + 2*NAD+_x

→ 2*ATP_x + 2*NADH_x + Butanol

2*Glicose + 29.7*ATP_x + 1.75*NADH_x

→ 3*Biomassa

9 ATP_x = ADP_x

2*ADP_x + Glicerol + 3*NAD+_x

→ 2*ATP_x + 3*NADH_x + AcetatoEx

3*ADP_x + 2*Glicerol + 4*NAD+_x

→ 3*ATP_x + 4*NADH_x + ButiratoEx

Glicose + 2*ADP_x

→ 2*ATP_x + 2*LactatoExt

Glicose + 2*ADP_x

→ 2*ATP_x + 2*Etanol

Glicose + 4*ADP_x + 4*NAD+_x

→ 4*ATP_x + 4*NADH_x + 2*AcetatoEx

Glicose + 3*ADP_x + 4*NAD+_x

→ 3*ATP_x + 4*NADH_x + Acetone

Glicose + 3*ADP_x + 4*NAD+_x

→ 3*ATP_x + 4*NADH_x + Acetone

Glicose + 2*ADP_x

→ 2*ATP_x + Butanol

Glicose + 3*ADP_x + 2*NAD+_x

→ 3*ATP_x + 2*NADH_x + ButiratoEx

19 ADP_x + Glicerol = ATP_x + H2 + Etanol

Glicose + 4*ADP_x + 2*NAD+_x

→ 4*ATP_x + 2*NADH_x + 2*H2 + 2*AcetatoEx

2*ADP_x + Glicerol + 2*NAD+_x

→ 2*ATP_x + 2*NADH_x + H2 + AcetatoEx

Glicose + 3*ADP_x + 2*NAD+_x

→ 3*ATP_x + 2*NADH_x + 2*H2 + Acetone

Glicose + 3*ADP_x + 2*NAD+_x

→ 3*ATP_x + 2*NADH_x + 2*H2 + Acetone

24 Glicose + 2*ADP_x + 2*NADH_x = 2*ATP_x + 2*NAD+_x + 2*H2 + Butanol

3*ADP_x + 2*Glicerol + 4*NAD+_x

→ 3*ATP_x + 4*NADH_x + 2*H2 + Acetone

3*ADP_x + 2*Glicerol + 4*NAD+_x

→ 3*ATP_x + 4*NADH_x + 2*H2 + Acetone

2*ADP_x + 2*Glicerol

→ 2*ATP_x + 2*H2 + Butanol

Glicose + 3*ADP_x

→ 3*ATP_x + 2*H2 + ButiratoEx

3*ADP_x + 2*Glicerol + 2*NAD+_x

→ 3*ATP_x + 2*NADH_x + 2*H2 + ButiratoEx

É possível associar um rendimento a cada modo elementar calculado para uma via.

Este rendimento, que é associado a um par substrato/produto – ambos externos –, mede a

quantidade em mol de produto produzido por mol de substrato.

Pode-se, no GEnSys, calcular o rendimento – associado a um modo elementar – para

qualquer par substrato/produto com comando

yields(model, substrato, produto)

Rendimento Glicose/H2

Os seis modos vistos na Figura 24 tem um rendimento de dois mol de H

por mol de

glicose.

yields(model,'Glicose','H2')

Yield(Glicose, H2)

Elmo yield

1 2.000

20 2.000

22 2.000

23 2.000

24 2.000

28 2.000

Rendimento Glicerol/H2

Há um rendimento de 1:1 quando se considera o glicerol como fonte de carbono em

cada um dos seis modos de produção do H

yields(model,'Glicerol','H2')

Yield(Glicerol, H2)

Elmo yield

19 1.000

21 1.000

25 1.000

26 1.000

27 1.000

29 1.000

Rendimento Glicose/ATP_x

Os modos 14 e 20 são os modos que têm os maiores rendimentos de ATP por

molécula de glicose. Observar que estes são os modos de produção do acetato (ácido). Os

modos 1 (etanol), 12 (lactato), 13 (etanol), 17 (butanol) e 24 (butano) rendem o menor

número de moléculas de ATP por molécula de glicose e estão associados à produção de

solventes. A produção de ácidos é maior durante o crescimento exponencial da bactéria. Isso

justifica uma maior produção de ATP nesta fase.

yields(model,'Glicose','ATP_x')

Yield(Glicose, ATP_x)

Elmo yield

1 2.000

12 2.000

13 2.000

14 4.000

15 3.000

16 3.000

17 2.000

18 3.000

20 4.000

22 3.000

23 3.000

24 2.000

28 3.000

Rendimento Glicerol/ATP_x

Já o rendimeno de ATP por molécula de glicerol é menor em relação ao da glicose.

Os maiores rendimentos de ATP são através dos modos 10 (acetato) e 21 (acetato e H

). A

diferença entre ambos é que o modo 21 tem a reação v

(reação do H

) e o modo 10 a v

. Se a

reação v

for eliminada, ainda assim, a bactéria poderia crescer normalmente produzindo a

mesma quantidade de ATP, porém produzindo H

yields(model,'Glicerol','ATP_x')

Yield(Glicerol, ATP_x)

Elmo yield

3 1.000

4 1.000

5 1.500

6 1.500

7 1.000

10 2.000

11 1.500

19 1.000

21 2.000

25 1.500

26 1.500

27 1.000

29 1.500

4.8.3 ABF da Produção de Biohidrogênio da C. acetobutylicum

Nesta seção será realizada a análise de balanço de fluxo para a C. acetobutylicum. A

análise de balanço de fluxo foi realizada considerando-se dois substratos: glicose e glicerol.

Esta análise é realizada, como foi visto, através da otimização (maximização ou minimização)

de alguns fluxos e restrição de outros.

A sintaxe geral para a criação de um modelo para ABF é mostrada no

Quadro 15;

onde

{}

p,,p,p

K e

{}

q,,q,q

K são subconjuntos do conjunto

{}

n,,3,2,1 K e não são,

necessariamente, disjuntos.

Quadro 15: Sintaxe para a construção de um modelo para ABF.

maxZ or minZ (minimizar: minZ; maximizar: maxZ)

lista de metabólitos

MIN

r_p

= vmin_1

r_p

= vmin_2

...

r_p

= vmin_k

MAX

r_q

= vmax_1

r_q

= vmax_2

...

r_q

= vmax_s

%n network reactions

r_1: reaction

r_2: reaction

...

r_n: reaction

O bloco maxZ (ou minZ) denota o objetivo da otimização (maximizar ou minimizar

a função objetivo Z) e quais os metabólitos que serão otimizados (lista dos metabólitos). Os

coeficientes para estes metabólitos são formados pela linha da matriz estequiométrica

relacionada aos metabólitos ali listados. O bloco MIN ou MAX contém os valores dos limites

inferiores (LB) ou superiores (UB), respectivamente, para os fluxos especificados nas reações

(velocidades mínimas e máximas). As reações são listadas no último bloco e devem ser

escritas em LRB; é destas reações que será criada a matriz estequiométrica. Os blocos de

dados: maxZ (ou minZ), MIN e MAX devem aparecer nesta ordem. maxZ (ou minZ), um

destes, e as reações são obrigatórias, os demais são opcionais. Caso um fluxo, por exemplo

não aparecer em um dos blocos MIN ou MAX, seus limites, por

default, serão assim

determinados:

 LB(i): ∞− , se r

for uma reação reversível;

 LB(i): 0, se r

for uma reação irreversível;

 UB(i): ∞ , tanto se r

for uma reação reversível ou irreversível.

A seguir, é apresentado um modelo para análise de balanço de fluxo usando a rede

metabólica da

C. acetobutylicum. O modelo, com todas as informações necessárias, foi

gravado no arquivo

H2CAcetobutylicum.m. Digitando o seguinte comando no MATLAB

aparece o conteúdo do arquivo:

type 'H2CAcetobutylicum.m'

maxZ

Biomassa

min

v2 = 0

max

v1 = 8.72

v2 = 0

%vGlyc

v1: Glicose + 2ATP -> 2ADP_x + 2GDH-3-P

%vGlicerol

v2: Glicerol + NAD+_x + ATP -> NADH + ADP_x + GDH-3-P

%vPIR

v3: GDH-3-P + 2ADP_x + NAD+_x -> PIR + 2ATP + NADH

%vPDH

v4: PIR + Fdoxi + CoA -> Acetil-CoA + Fdred

%vNADH

v5: Fdoxi + NADH <-> Fdred + NAD+_x

%vHYD

v6: Fdred -> Fdoxi + H2

%vLAC

v7: PIR + NADH <-> NAD+_x + Lactato

%vTHI

v8: 2Acetil-CoA -> CoA + AcAcetil-CoA

%vPTA

v9: Acetil-CoA + ADP_x <-> ATP + CoA + Acetato

%vADH_Ethanol

v10: Acetil-CoA + 2NADH -> CoA + 2NAD+_x + Etanol

%vCOAT1_Acetone

v11: Acetato + AcAcetil-CoA -> Acetil-CoA + Acetone

%vCOAT2_Acetone

v12: Butirato + AcAcetil-CoA -> Butiril-CoA + Acetone

%vHCOA

v13: AcAcetil-CoA + 2NADH -> Butiril-CoA + 2NAD+_x

%vPTB

v14: Butiril-CoA + ADP_x <-> CoA + ATP + Butirato

%vBADH_Butanol

v15: Butiril-CoA + 2NADH -> CoA + 2NAD+_x + Butanol

%vCrescimento

v16: 0.6667Glicose + 0.5833NADH + 9.9ATP -> Biomassa

%vATPase

v17 : ATP -> ADP_x

%Reações de transporte do interior da bactéria para o ambiente

v18: Lactato -> LactatoExt

v19: Acetato -> AcetatoExt

v20: Butirato -> ButiratoExt

Figura 25 exibe a equação (2.30), 0Nv

, para as reações mostradas acima. Esta

equação é uma restrição de balanço de fluxo para esta análise.

N é a matriz estequiométrica

interna com 13 linhas (uma para cada metabólito interno) e 20 reações, onde as reações

possuem os nomes

v20,v2,v1, K .

Figura 25: Restrição de balanço de fluxo para a rede da C. acetobutylicum mostrada da Figura

23. Esta figura é um exemplo da equação

(2.30): 0Nv

No modelo acima, deseja-se maximizar a biomassa com velocidade máxima de

consumo de glicose (v

) de 8.72 mmol h

-1

gdw

-1

(gdw: grama de massa celular seca).

Maximizar a biomassa corresponde a simular o comportamento da bactéria, supondo que o

objetivo dela é sempre maximizar o seu crescimento. A velocidade de consumo do glicerol é

zero. As velocidades mínimas para as demais reações é zero para as irreversíveis e -

∞ para as

reversíveis. Para as velocidades máximas é atribuído

∞, independentemente da reversibilidade

da reação.

A função

fba deve ser chamada para poder criar um modelo FBA:

FBAModel = fba('H2CAcetobutylicum.m')

FBAModel =

intModel: [1x1 struct]

beq: [13x1 double]

c: [0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 -1 0 0 0 0]

lb: [20x1 double]

ub: [20x1 double]

minMax: -1

A função fba retornará a estrutura FBAmodel contendo os seguintes campos:

 intModel: estrutura que contém o modelo interno (apenas os metabólitos internos). Os

significados dos campos desta estrutura já foram definidos anteriormente. Um dos campos

é a matriz

st, matriz estequiométrica (m mebabólitos por n reações).

 beq: vetor m×1 com todos elementos iguais a zero. Reservado para uso futuro, como

recurso para adicionar fluxos com valores conhecidos;

-2 -1 2 0 0 0 0 0 1 0 0 0 0 1 0 -9,9 -1 0 0 0

2 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 1 1 0 -1 0 -1 0 0 -2 0 0 -2 0 -2 -0,58 0 0 0 0

0 0 1 -1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 -1 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 -1 0 0 0 1 1 1 0 0 0 1 1 0 0 0 0 0

0 0 0 1 0 0 0 -2 -1 -1 1 0 0 0 0 0 0 0 0 0

0 0 0 1 1 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 -1 0 0

0 0 0 0 0 0 0 1 0 0 -1 -1 -1 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 1 0 -1 0 0 0 0 0 0 0 -1 0

0 0 0 0 0 0 0 0 0 0 0 -1 0 1 0 0 0 0 0 -1

0 0 0 0 0 0 0 0 0 0 0 1 1 -1 -1 0 0 0 0 0

v10

v11

v12

v13

v14

v15

v16

v17

v18

v19

v20

 c: vetor n×1 contendo os n coeficientes da função objetivo. É determinado pela soma das

linhas correspondentes aos metabólitos a serem otimizados;

 lb: vetor n×1 onde são armazenados os limites inferiores das velocidades dos fluxos

(velocidade mínima);

 ub: vetor n×1 onde são armazenados os limites superiores das velocidades dos fluxos

(velocidade máxima);

 minMax: tipo do objetivo: -1 para maximizar e 1 para minimizar.

Para determinar a distribuição de fluxos com as restrições impostas ao modelo, usa-

se o comando:

sol = solveFBA(FBAModel)

sol =

f: 3.5232

x: [20x1 double]

y: [14x1 double]

w: [20x1 double]

stat: 1

solver: 'glpk'

time: 0

A solução encontra-se no campo x, vetor coluna de tamanho n (20 neste caso), da

estrutura

sol (solução). O campo f é o valor da função objetivo e vale 3.5232 mmol/gdw

-1

, o

mesmo valor de v

. Para mostrar a solução usa-se o comando (exibe os nomes dos fluxos

seguidos dos valores calculados):

printFluxVector(FBAModel.intModel.lRxn, sol.x). O

resultado encontra-se na parte esquerda da

Figura 26.

Figura 26: Distribuição de fluxo da ABF ao se maximizar a biomassa com velocidade de

consumo de glicose 8,72 mmol/gdw

-1

8,7 0

17,44

15,4

17,44

3,52

v1 = 8.72

v2 = 0

v3 = 17.44

v4 = 17.44

v5 = 15.3849

v6 = 32.8249

v7 = 0

v8 = 0

v9 = 17.44

v10 = 3.2e-006

v11 = 0

v12 = 0

v13 = 0

v14 = 0

v15 = 0

v16 = 3.5232

v17 = 0

v18 = 0

v19 = 17.44

v20 = 0

32,82

O resultado acima mostra que a velocidade de produção de hidrogênio (v

) é de

32,8249 mmol/gdw

-1

, a de consumo de glicose é 8,72 mmol/gdw

-1

). Assim o

rendimento de H

é de 3,7643 (32,8249 mmol/gdw

-1

/8,72 mmol/gdw

-1

) mol de H

por

mol de glicose consumida. Além do hidrogênio, é produzido apenas o acetato (v

= 17,44

mmol/gdw

-1

). Este resultado é condizente com dados experimentais; o hidrogênio é

produzido durante a fase de produção de ácidos que também é a fase de crescimento

exponencial da

C. acetobutylicum (JONES e WOODS, 1986).

A próxima simulação é de maximização do hidrogênio. Para tanto, deve-se mudar a

biomassa para H2 no bloco maxZ do arquivo:

maxZ

min

v2 = 0

max

v1 = 8.72

v2 = 0

Para este novo objetivo, deve-se executar o modelo novamente (são três comandos):

 m = fba('H2CAcetobutylicum.m');

 s = solveFBA(m);

 printFluxVector (m.intModel.lRxn, s.x).

O resultado é apresentado na

Figura 27.

Figura 27: Distribuição de fluxo da ABF ao se maximizar o H

com velocidade máxima de

consumo de glicose 8,7 mmol/gdw

-1

v1 = 8.72

v2 = 0

v3 = 17.44

v4 = 17.44

v5 = 17.44

v6 = 34.88

v7 = 0

v8 = 8.72

v9 = 0

v10 = 0

v11 = 0

v12 = 8.72

v13 = 0

v14 = 8.72

v15 = 0

v16 = 0

v17 = 26.16

v18 = 0

v19 = 0

v20 = 0

8,7 0

17,4

17,44

17,44 34,88

8,72

7,72

8,72 0

26,16

O rendimento de H

nesta simulação é 4 mol H

/mol glicose (34,88/8,72). Este é o

rendimento teórico máximo de produção de H

por molécula de glicose. Quando se utiliza

glicerol como único substrato, este rendimento máximo é de 3 mol H

/mol glicerol (dados não

mostrados) quando se maximiza apenas o H

. Mas, como pode ser visto acima, a velocidade

de crescimento, v

, é zero. Isto quer dizer que este rendimento de H

não poderá ser atingido

experimentalmente. Manish

et al. (2007) chegaram a este mesmo resultado ao maximizar a

produção de H

na análise de balanço de fluxo de uma rede metabólica da E. coli.

Foi realizada outra análise de balanço de fluxo usando uma velocidade máxima de

3,94 mmol/gdw

-1

e 4,65 mmol/gdw

-1

para as velocidades de consumo de glicose e

glicerol, respectivamente. Estes valores, v

e v

, foram determinados experimentalmente e

estão na

Tabela 11.

Tabela 11: Dados experimentais relacionados à fase solventogênica (VASCONCELOS et al.,

1994).

Fluxos medidos (mmol h

-1

gdw

-1

)

Fonte de carbono Produtos

Glicose

)

Glicerol

)

Etanol

)

Butanol

)

Acetona

e v

)

Acetato

)

Butirato

)

3,94 4,65 1,3 3,86 0 0,42 0,66 7,24

Figura 28 são apresentados os resultados para as três simulações distintas de

maximização de:

 biomassa;

 biomassa e H

;

 biomassa, butirato e acetato;

 butirato e acetato.

Cada grupo de barras na

Figura 28 representa uma velocidade e cada barra de um

grupo uma das simulações descritas acima. O fluxo de produção de acetona é representado

por v11v12, uma vez que ambos os fluxos v11 e v12 produzem este metabólito.

v1 v2 v6 v10 v11v12 v15 v1 6 v1 9 v20

Rea

ões

Fluxo mmol/(gdw)

max. biomassa

max. biomassa + H

max. biomassa + butirato + acetato

max. butirato + acetato

Figura 28: Cada grupo de barras representa um fluxo: v1 e v2 correspondem aos fluxos de

consumo de glicose e glicerol, respectivamente; v6, v10, v11v12, v15, v16, v19 e v20

correspondem aos fluxos de produção de H

, etanol, acetona, butanol, biomassa, acetato e

butirato. Em todas as maximizações, o valor máximo de consumo de glicose (v1) e de glicerol

(v2) é de 3,94 mmol/gdw

-1

e de 4,65 mmol/gdw

-1

, respectivamente.

Em todos os casos, sem exceção, as distribuições de fluxos calculadas para a rede

metabólica coincidiram (

Figura 28). São produzidos apenas H

= 28,23 mmol/gdw

-1

) e

acetato (v

= 12,53 mmol/gdw

-1

). A velocidade de crescimento (v

) é de 2,53

mmol/gdw

-1

(Figura 29).

Como já foi dito, durante a velocidade de crescimento exponencial é que há a maior

produção de ácidos e de H

. Nestas simulações foi obtido este comportamento; ao se

maximizar a biomassa produz-se muito ácidos e H

; ao se maximizar H

e biomassa, produz-

se muito ácido; e, por último, ao se maximizar ácido produz-se muito H

e biomassa.

Portanto, este comportamento

in silico está condizente com o comportamento padrão da C.

acetobutylicum

Figura 29: Distribuição de fluxo obtida pela ABF ao se maximizar a biomassa ou a

biomassa+H2 ou a biomassa+butirato+acetato ou o butirato+acetato.

A distribuição de fluxo apresentada na

Tabela 11 corresponde a dados experimentais

obtidos na fase solventogênica (VASCONCELOS

et al., 1994). Para simular esta fase, foi

realizada a análise de balanço de fluxo na qual a velocidade de produção do solvente butanol

foi maximizada. As restrições de velocidades máximas de consumo de glicose e gliceral

permaneceram as mesmas. Foi adicionada a restrição para a velocidade mínima de etanal

(v10) de 1,3 mmol/gdw

-1

A distribuição de fluxo para esta simulação é mostrada na

Figura 30. Pode-se ver

que, qualitativamente, os valores obtidos pela simulação se parecem com os obtidos

experimentalmente. Os valores determinados pela análise de balanço de fluxo são v6, v11v12,

v15, v19 e v20. Com a exceção do fluxo v11v12 (acetona), cujos valores coincidem – ambos

são zero –, os demais fluxos formam um padrão semelhante aos experimentais.

v1 = 3.94

v2 = 4.65

v3 = 12.53

v4 = 12.53

v5 = 15.70

v6 = 28.23

v7 = 0

v8 = 0

v9 = 12.53

v10 = 0

v11 = 0

v12 = 0

v13 = 0

v14 = 0

v15 = 0

v16 = 2.53

v17 = 0

v18 = 0

v19 = 12.53

v20 = 0

3,94 4,65

12,53

15,7 28,23

12,53

2,53

12,53

v1 v2 v6 v10 v11v12 v1 5 v19 v2 0

Reações

Fluxo mmol/(gdw)

dados experimentais

max. butanol

Figura 30: Dados experimentais versus simulados, ambos na fase solventogênica. v1 e v2

correspondem aos fluxos de consumo de glicose e glicerol, respectivamente; v6, v10, v11v12,

v15, v19 e v20 correspondem aos fluxos de produção de H

, etanol, acetona, butanol, acetato e

butirato. O valor máximo de consumo de glicose (v1) e de glicerol (v2) é de 3,94 mmol/gdw

-1

e de 4,65 mmol/gdw

-1

, respectivamente. O valor mínimo da velocidade de produção de

etanol (v10) é de 1,3 mmol/gdw

-1

Foi feita também uma simulação para a fase acidogênica e os resultados foram

comparados aos dados experimentais desta mesma fase (

Tabela 12). Como foi visto acima,

independentemente do que for maximizado – se biomassa ou ácidos – o resultado é o mesmo.

Assim, optou-se por, primeiro, maximizar a biomassa juntamente com o H

e, depois, apenas

a biomassa. Em ambas as simulações, a única fonte de carbono é a glicose, com uma

velocidade máxima de consumo (v1) igual a 8,72 mmol/gdw

-1

. Na segunda simulação –

maximização de biomassa – impôs-se uma velocidade máxima de 3,67 mmol/gdw

-1

para o

acetato.

Tabela 12: Dados experimentais relacionados à fase acidogênica (VASCONCELOS et al.,

1994).

Fluxos medidos (mmol h

-1

gdw

-1

)

Fonte de carbono Produtos

Glicose

)

Glicerol

)

Etanol

)

Butanol

)

Acetona

e v

)

Acetato

)

Butirato

)

8,72 0 0,26 0,01 0,02 3,67 6,09 19,1

100

Os resultados estão mostrados na

Figura 31. Ao se maximizar biomassa e H

não

houve produção de butirato (v20) e a velocidade de produção de H

e acetato foi maior que as

suas velocidades experimentais.

Ao se maximizar apenas a biomassa e limitando a velocidade de consumo da glicose

(v1) e do acetato (v19), os resultados da simulação são muitos próximos dos experimentais;

comparar as barras cinzas com as pretas na

Figura 31.

v1 v6 v10 v11v12 v15 v19 v20

Reações

Fluxo mmol/(gdw)

max. biomassa + H

max. biomassa

dados experimentais

Figura 31: Dados experimentais versus simulados, ambos na fase acidogênica. v1 corresponde

ao fluxo de consumo de glicose; v6, v10, v11v12, v15, v19 e v20 correspondem aos fluxos de

produção de H

, etanol, acetona, butanol, acetato e butirato. O valor máximo de consumo de

glicose (v1) é de 8,72 mmol/gdw

-1

. A velocidade máxima de acetato (v19) é de 3,67

mmol/gdw

-1

Estes resultados mostram que o modelo proposto poderá ser útil para guiar pesquisas

e análises de propriedades biológicas da

C. acetobutylicum, o que pode levar na direção de

uma metodologia para melhorar a produção de H

e solventes por esta bactéria.

5 CONSIDERAÇÕES FINAIS

5.1 CONCLUSÕES

Como neste trabalho foi descrito um programa para análise e simulação de redes

metabólicas e alguns estudos de caso de modelagem e simulação, as conclusões são separadas

em seções distintas.

5.1.1 GEnSys

Neste trabalho foi descrito o GEnSys (Genomic Engineering System) que é formado

por um conjunto de módulos interdependentes que permite a análise e simulação de redes

bioquímicas. Os módulos implementados realizam as seguintes tarefas:

 processamento da LRB para a obtenção da matriz estequiométrica, da lista dos nomes dos

metabólitos e da lista de reações da rede;

 determinação dos metabólitos internos e externos da rede e suas matrizes estequiométricas

correspondentes;

 determinação das linhas e colunas, em relação à matriz estequiométrica interna,

linearmente independentes e dependentes. Isto é particularmente útil para a determinação

dos

pools de metabólitos conservados, isto é, aqueles metabólitos cuja soma de suas

concentrações é invariante no tempo. A determinação destes metabólitos reduz o número

de EDOs na simulação determinística;

 modelagem e simulação determinística de uma forma compacta, onde, num mesmo

arquivo, pode-se definir o modelo contendo as concentrações iniciais, as cinéticas pré-

definidas e as reações da rede metabólica;

 análise de controle metabólico para uma rede a partir de um modelo determinístico. Todos

os coeficientes de controle são calculados invocando uma única função;

 análise de fluxo metabólico: a partir dos fluxos medidos é possível calcular todos os

demais fluxos da rede. O modelo de AFM é facilmente definido num arquivo onde se

colocam tais fluxos medidos e todas as reações da via. A solução do sistema da AFM é

por quadrados mínimos; isto significa que sempre há uma solução e o vetor solução tem

norma mínima;

 análise de balanço de fluxo: para esta análise é necessário especificar, além da equação

geral de conservação da massa, alguns valores de fluxos mínimos e máximos e o

metabólito/fluxo que se deseja otimizar. Todos estes componentes são colocados no

102

arquivo do modelo para serem processados e a distribuição de fluxo determinada de forma

a otimizar a função objetivo;

 análise estrutural através de modos elementares: um modo elementar é definido como o

menor conjunto de reações operando conjuntamente para manter o regime estacionário.

Isto significa que um modo elementar é uma via que liga um ou mais metabólitos de

entrada a um ou mais metabólitos de saída. Para determinar os modos elementares de uma

rede usando o GEnSys basta especificar as reações e os metabólitos externos da rede.

Para efetuar qualquer uma das análises descritas anteriormente, é necessário,

priori

, criar uma instância (definição dos valores dos campos internos) do módulo de análise.

Somente depois de criada tal instância é que, de fato, realiza-se a análise, usando, para tanto,

as funções definidas especificamente para aquele módulo.

O GEnSys faz uso intensivo da LRB (Linguagem das Reações Bioquímicas) que é

uma linguagem de reações de propósito geral e flexível, desenvolvida neste trabalho. Até

onde se sabe, não há uma linguagem tão flexível e com tantos recursos para a escrita dos

nomes dos metabólitos. As linguagens de reações existentes são bastante limitadas quanto aos

símbolos usados nos nomes dos compostos e nos sinais usados nas reações.

Da forma que foi projetado, o GEnSys pode ser facilmente extendido para outros

tipos de métodos de análise envolvendo a matriz estequiométrica. Todos os módulos foram

cuidadosamente definidos para conterem um conjunto de dados útil de forma que,

dependendo do tipo de análise, sua reutilização seja simples. Por exemplo, o módulo

totalStMat contém, separadamente, as duas matrizes estequiométricas: uma relativa a todos os

substratos e outra a todos os produtos; estas duas matrizes são necessárias para a simulação

dinâmica estocástica.

Finalmente, o GEnSys pode ser muito útil, tanto para usuário que deseja desenvolver

algum programa/módulo, quanto para aquele que deseja fazer apenas uma modelagem e

simulação que envolve a matriz estequiométrica. O primeiro, geralmente, é da área de

ciências exatas e o segundo da área de biologia ou bioquímica. Estes últimos podem

facilmente usar o GEnSys, pois para se criar um modelo neste programa basta editar um texto

no formato ASCII e conhecer as sintaxes usadas pelo GEnSys. Além disso, necessita conhecer

um mínimo do ambiente M

ATLAB

para invocar uma função neste ambiente.

103

5.1.2 Simulação Estocástica da E. coli

Foi apresentado e construído um conjunto de reações bioquímicas representando a

via de biossíntese do Trp em

E. coli. Usou-se o algoritmo de simulação estocástica

desenvolvido por Gibson e Bruck (2000) para fazer as simulações estocásticas do

operon trp

desta bactéria. A maioria das constantes de velocidades e concentrações iniciais foi obtida de

Santillan e Zeron (2001b). Em ambas as simulações, o resultado, especificamente a atividade

da enzima AS, previu qualitativamente, para uma

E. coli selvagem, valores próximos aos

dados experimentais.

Todos os mecanismos de regulação conhecidos do

operon trp foram levados em

consideração nesta simulação. Mas foram feitas algumas simplificações, como considerar

apenas a enzima AS, bem como um número constante de moléculas para o repressor

trpR, que

é auto regulado. Outra simplificação foi não considerar a importação ou exportação do Trp

do, ou para, o meio externo à célula.

A simulação estocástica é uma excelente ferramenta de simulação para um conjunto

de reações químicas, principalmente para redes de sinalização, onde o número de moléculas é

pequeno, como o operador

trp ou o TrpR (repressor), comparado ao volume da célula e o

número de outros metabólitos que podem estar presentes em grandes quantidades.

O principal problema a ser superado pela simulação estocástica é o tempo de

execução do algoritmo. A solução do modelo determinístico – solução das EDOs por métodos

numéricos – é muito mais rápida do que a simulação estocástica. Na simulação estocástica,

cada reação é executada separadamente em cada passo do tempo, o que torna sua execução

mais lenta. Sendo assim, pode-se dizer que a simulação estocástica é mais adequada para

pequenas redes metabólicas e que os compostos envolvidos não devem ter altas

concentrações.

5.1.3 Modelagem e Simulação da Rede Metabólica da C. acetobutylicum

Foi construída, a partir de dados da literatura, uma rede metabólica da C.

acetobutylicum

envolvendo todas as reações usadas para a produção de ácidos, solventes e

hidrogênio molecular a partir de dois substratos distintos: glicose e glicerol. No total, a rede

possui 25 metabólitos e 20 reações.

Foram efetuadas duas abordagens diferentes para o estudo da rede metabólica da

bactéria

C. acetobutylicum; a saber: análise de balanço de fluxo e a determinação dos modos

104

elementares desta rede. Até onde se sabe, nenhuma destas análises já foi realizada levando-se

em conisderação o glicerol e a glicose como fontes de carbono.

Foram determinados 29 modos elementares de fluxo para esta via. Destes, há 13

modos distintos envolvidos na produção de H

; seis a partir da glicose e seis a partir do

glicerol, o décimo terceiro é a partir do NADH. Foi possível também calcular o rendimento de

por molécula de substrato. Todos os seis modos de produçao de H

a partir de uma

molécula de glicose rendem duas de H

. Já para o glicerol este rendimento é de 1:1. Isto

acontece porque para cada molécula de glicose são produzidas duas de piruvato e para o

glicerol é produzida apenas uma.

A determinação dos modos elementares foi útil para elucidar as possíveis vias de

produção de H

por esta rede. Dependendo do tamanho da rede, não são possíveis de serem

determinadas.

Através da ABF foram realizadas nove simulações para esta rede metabólica. Estas

análises, com exceção de uma – que deu velocidade de crescimento zero quando se maximiza

o H

–, obtiveram, qualitativamente, um resultado que está de acordo com padrão de

comportamento, em termos de distribuição de fluxos, desta bactéria. Na última simulação, os

valores da distribuição de fluxos obtidos foram muito próximos aos valores dos dados

experimentais.

Estes resultados mostram que este modelo pode ser usado para estudar

in silico

possíveis formas de melhorar a produção de H

, ou até mesmo de solventes, na C.

acetobutylicum

. O modelo poderá ser usado para simular nocautes de genes (remoção de

reações) e/ou adição de novas reações.

5.2 RECOMENDAÇÕES PARA TRABALHOS FUTUROS

Conforme mostrado na primeira parte deste capítulo, o GEnSys foi projetado para

facilitar a incorporação de outros tipos de análises. Uma delas é a simulação dinâmica

estocástica que seria facilmente implementada usando os módulos do GEnSys.

Uma outra abordagem seria a atribuição de fluxos a modos elementares. Por

exemplo, foi visto na ABF que o rendimento máximo teórico é de 4 moléculas de H

por

molécula de glicose. Também foi visto que há vários modos elementares de produção de H

No entanto, não se sabe quais dos modos dão este rendimento de 4:1. Na verdade, é uma

combinação de mais de um dos modos elementares. A atribuição de fluxos aos modos

elementares pode ser uma solução para este problema, pois se determinaria um fluxo (com

105

unidade) para cada um dos módulos da rede. Este seria um módulo facilmente incorporado ao

GEnSys. Esta atribuição pode ser resolvida como um problema de quadrados mínimos com

restrição de desigualdade. Também pode ser visto como um problema de otimização e ser

resolvido por uma técnica de inteligência artificial conhecida como algoritmos genéticos, por

exemplo.

Conforme demonstrado, a AFM considera todas as reações reversíveis. Seria

possível também implementar um módulo com AFM levando em consideração as restrições

de irreversibilidades de alguma reação. Este, também, pode ser visto como problema dos

quadrados mínimos com restrição de desigualdade.

Com relação à produção do H

, existe a possibilidade de realizar análises levando em

consideração a adição ou supressão de enzimas na rede, o que equivale a uma adição ou

nocaute de genes.

REFERÊNCIAS

AKASHI, H.; GOJOBORI, T. Metabolic efficiency and amino acid composition in the

proteomes of

Escherichia coli and Bacillus subtilis. Proceedings of the National Academy

of Sciences, USA

, v.99, n.6, p.3695-3700. 2002.

BARRETT, C. L.; PRICE, N. D.; PALSSON, B. O. Network-level analysis of metabolic

regulation in the human red blood cell using random sampling and singular value

decomposition.

BMC Bioinformatics, v.7, p.132. 2006.

BECKER, S. A.; FEIST, A. M.; MO, M. L.; HANNUM, G.; PALSSON, B. O.;

HERRGARD, M. J. Quantitative prediction of cellular metabolism with constraint-based

models: the cobra toolbox.

Nature Protocols, v.2, n.3, p.727-738. 2007.

BLISS, R. D.; PAINTER, P. R.; MARR, A. G. Role of feedback inhibition in stabilizing the

classical operon.

Journal of Theoretical Biology, v.97, n.2, p.177-193. 1982.

BURRAGE, K.; TIAN, T. H.; BURRAGE, P. A multi-scaled approach for simulating

chemical reaction systems.

Progress in Biophysics & Molecular Biology, v.85, n.2-3, p.217-

234. 2004.

CAO, Y.; GILLESPIE, D. T.; PETZOLD, L. R. Efficient step size selection for the tau-

leaping simulation method.

Journal of Chemical Physics, v.124, n.4, p.044109. 2006.

CORNISH-BOWDEN, A.; HOFMEYR, J. H. The role of stoichiometric analysis in studies of

metabolism: an example.

Journal of Theoretical Biology, v.216, n.2, p.179-191. 2002.

COVERT, M. W.; SCHILLING, C. H.; PALSSON, B. Regulation of gene expression in flux

balance models of metabolism.

Journal of Theoretical Biology, v.213, n.1, p.73-88. 2001.

DESAI, R. P.; HARRIS, L. M.; WELKER, N. E.; PAPOUTSAKIS, E. T. Metabolic flux

analysis elucidates the importance of the acid-formation pathways in regulating solvent

production by

Clostridium acetobutylicum. Metabolic Engineering, v.1, n.3, p.206-213.

1999.

DESAI, R. P.; NIELSEN, L. K.; PAPOUTSAKIS, E. T. Stoichiometric modeling of

Clostridium acetobutylicum

fermentations with non-linear constraints. Journal of

Biotechnology

, v.71, n.1-3, p.191-205. 1999.

EDWARDS, J. S.; COVERT, M.; PALSSON, B. Metabolic modelling of microbes: the flux-

balance approach.

Environmental Microbiology, v.4, n.3, p.133-140. 2002.

EHLDE, M.; ZACCHI, G. A general formalism for metabolic control analysis.

Chemical

Engineering Science

, v.52, n.15, p.2599-2606. 1997.

FELL, D. A. Metabolic control analysis: a survey of its theoretical and experimental

development.

Biochemical Journal, v.286, n.31, p.313–330. 1992.

______. Increasing the flux in metabolic pathways: a metabolic control analysis perspective.

Biotechnology and Bioengineering, v.58, n.2-3, p.121-124. 1998.

107

GIANCHANDANI, E. P.; PAPIN, J. A.; PRICE, N. D.; JOYCE, A. R.; PALSSON, B. O.

Matrix formalism to describe functional states of transcriptional regulatory systems.

PLoS

Compuational Biology

, v.2, n.8, p.e101. 2006.

GIBSON, M. A.; BRUCK, J. Efficient exact stochastic simulation of chemical systems with

many species and many channels.

Journal of Chemical Physics A, v.104, n.9, p.1876-1889.

2000.

GILLESPIE, D. T. Exact stochastic simulation of coupled chemical-reactions.

Journal of

Physical Chemistry

, v.81, n.25, p.2340-2361. 1977.

GOLLNICK, P.; BABITZKE, P. Transcription attenuation.

Biochimica et Biophysica Acta-

Gene Structure and Expression

, v.1577, n.2, p.240-250. 2002.

GOODWIN, B. C. Oscillatory behavior in enzymatic control processes.

Advances in

Enzyme Regulation

, v.3, p.425-438. 1965.

GRIFFITH, J. S. Mathematics of cellular control processes. I. Negative feedback to one gene.

Journal of Theoretical Biology, v.20, n.2, p.202-208. 1968a.

______. Mathematics of cellular control processes. II. Positive feedback to one gene.

Journal

of Theoretical Biology

, v.20, n.2, p.209-216. 1968b.

HYNNE, F.; DANO, S.; SORENSEN, P. G. Full-scale model of glycolysis in

Saccharomyces

cerevisiae

. Biophysical Chemistry, v.94, n.1-2, p.121-163. 2001.

JONES, D.; WOODS, D. Acetone-butanol fermentation revisited.

Microbiology and

Molecular Biology Reviews

, v.50, n.4, p.484-524. 1986.

KACZER, H.; BURNS, J. A. The control of flux.

Journal of Society of Experimental

Botany

, v.27, n.256, p.65-104. 1973.

KAUFFMAN, K. J.; PRAKASH, P.; EDWARDS, J. S. Advances in flux balance analysis.

Current Opinion of Biotechnology, v.14, n.5, p.491 - 496. 2003.

KLAMT, S.; STELLING, J. Combinatorial complexity of pathway analysis in metabolic

networks.

Molecular Biology Reports, v.29, n.1-2, p.233-236. 2002.

KLIPP, E.; HERWIG, R.; KOWALD, A.; WIERLING, C.; LEHRACH, H.

Systems biology

in practice: concepts, implementation and application

. Weinheim: Wiley-VCH. 2005. xix,

465 p.

LAWSON, C. L.; HANSON, R. J.

Solving least squares problems: Englewood Cliffs (N.J.):

Prentice-Hall. 1995

LEHNINGER, A. L.; NELSON, D. L.; COX, M. M.

Principles of biochemistry. New York:

W.H. Freeman & Company. 2000

LESKOVAC, V.

Comprehensive enzyme kinetics. New York: Kluwer Academic Press.

2003. xi, 438 p.

108

LLANERAS, F.; PICÓ, J. Stoichiometric modelling of cell metabolism.

Journal of

Bioscience and Bioengineering

, v.105, n.1, p.1-11. 2008.

MANISH, S.; VENKATESH, K. V.; BANERJEE, R. Metabolic flux analysis of biological

hydrogen production by

Escherichia coli. International Journal of Hydrogen Energy, v.32,

n.16, p.3820-3830. 2007.

MARANGONI, A. G.

Enzyme kinetics: a modern approach. New Jersey: John Wiley &

Sons. 2003. 241 p.

MCADAMS, H. H.; ARKIN, A. Stochastic mechanisms in gene expression.

Proceedings of

National Academy of Science, USA

, v.94, n.3, p.814-819. 1997.

MCCOLLUM, J. M.; PETERSON, G. D.; COX, C. D.; SIMPSON, M. L.; SAMATOVA, N.

F. The sorting direct method for stochastic simulation of biochemical systems with varying

reaction execution behavior.

Computational Biology and Chemistry, v.30, n.1, p.39-49.

2006.

MENDES, E.

Uso racional de ferramentas de engenharia metabólica: produção

bacteriana de hidrogênio e de violaceína

. 2006. 119 f. Dissertação (Mestrado em

Engenharia Química) - Programa de Pós-Graduação em Engenharia Química, Universidade

Federal de Santa Catarina, Florianópolis.

NIELSEN, J. H.; VILLADSEN, J.

Bioreaction engineering principles. New York: Plenum

Press. 1994. xxiv, 456 p.

NOLLING, J.; BRETON, G.; OMELCHENKO, M. V.; MAKAROVA, K. S.; ZENG, Q.;

GIBSON, R.; LEE, H. M.; DUBOIS, J.; QIU, D.; HITTI, J.; WOLF, Y. I.; TATUSOV, R. L.;

SABATHE, F.; DOUCETTE-STAMM, L.; SOUCAILLE, P.; DALY, M. J.; BENNETT, G.

N.; KOONIN, E. V.; SMITH, D. R. Genome sequence and comparative analysis of the

solvent-producing bacterium

Clostridium acetobutylicum. Journal of Bacteriology, v.183,

n.16, p.4823-4838. 2001.

OH, Y.-K.; KIM, H.-J.; PARK, S.; KIM, M.-S.; RYU, D. D. Y. Metabolic-flux analysis of

hydrogen production pathway in

Citrobacter amalonaticus Y19. International Journal of

Hydrogen Energy

, v.33, n.5, p.1471-1482. 2008.

PALSSON, B.

Systems biology: properties of reconstructed networks. Cambridge; New

York: Cambridge University Press. 2006. xii, 322 p.

PAPOUTSAKIS, E. T. Equations and calculations for fermentations of butyric acid bacteria.

Biotechnology and Bioengineering, v.26, n.2, p.174-187. 1984.

PENROSE, R. A generalized inverse for matrices.

Proceedings of Cambridge Philosophy

Society

, v.51, n.406-413. 1955.

RAMSEY, S.; ORRELL, D.; BOLOURI, H. Dizzy: stochastic simulation of large-scale

genetic regulatory networks.

Journal of Bioinformatics and Computational Biology, v.3,

n.2, p.437-454. 2005.

109

RAO, C. V.; ARKIN, A. P. Stochastic chemical kinetics and the quasi-steady-state

assumption: application to the Gillespie algorithm.

Journal of Chemical Physics, v.118,

n.11, p.4999-5010. 2003.

RAO, G.; MUTHARASAN, R. Altered electron flow in continuous cultures of

Clostridium

acetobutylicum

induced by viologen dyes. Applied Environmental Microbiology, v.53, n.6,

p.1232-1235. 1987.

REDER, C. Metabolic control theory: a structural approach.

Journal of Theoretical Biology,

v.135, n.2, p.175-201. 1988.

SANTILLAN, M.; MACKEY, M. C. Dynamic behavior in mathematical models of the

tryptophan operon.

Chaos, v.11, n.1, p.261-268. 2001a.

______. Dynamic regulation of the tryptophan operon: a modeling study and comparison with

experimental data.

Proceedings of the National Academy of Sciences, USA, v.98, n.4,

p.1364-1369. 2001b.

SANTILLAN, M.; ZERON, E. S. Dynamic influence of feedback enzyme inhibition and

transcription attenuation on the tryptophan operon response to nutritional shifts.

Journal of

Theoretical Biology

, v.231, n.2, p.287-298. 2004.

SAURO, H. M.; INGALLS, B. Conservation analysis in biochemical networks:

computational issues for software writers.

Biophysical Chemistry, v.109, n.1, p.1-15. 2004.

SCHILLING, C. H.; LETSCHER, D.; PALSSON, B. O. Theory for the systemic definition of

metabolic pathways and their use in interpreting metabolic function from a pathway-oriented

perspective.

Journal of Theoretical Biology, v.203, n.3, p.229-248. 2000.

SCHILLING, C. H.; SCHUSTER, S.; PALSSON, B. O.; HEINRICH, R. Metabolic pathway

analysis: basic concepts and scientific applications in the post-genomic era.

Biotechnology

Progress

, v.15, n.3, p.296-303. 1999.

SCHUSTER, S.; DANDEKAR, T.; FELL, D. A. Detection of elementary flux modes in

biochemical networks: a promising tool for pathway analysis and metabolic engineering.

Trends Biotechnology, v.17, n.2, p.53-60. 1999.

SCHUSTER, S.; FELL, D. A.; DANDEKAR, T. A general definition of metabolic pathways

useful for systematic organization and analysis of complex metabolic networks.

Nature

Biotechnology

, v.18, n.3, p.326-332. 2000.

SCHUSTER, S.; HILGETAG, C. What information about the conserved-moiety structure of

chemical reaction systems can be derived from their stoichiometry?

The Journal of Physical

Chemistry

, v.99, n.20, p.8017-8023. 1995.

SCHUSTER, S.; HILGETAG, C.; WOODS, J. H.; FELL, D. A. Reaction routes in

biochemical reaction systems: algebraic properties, validated calculation procedure and

example from nucleotide metabolism.

Journal of Mathematical Biology, v.45, n.2, p.153-

181. 2002.

110

SEN, A. K.; LIU, W. M. Dynamic analysis of genetic, control and regulation of amino-acid

synthesis - the tryptophan operon in

Escherichia coli. Biotechnology and Bioengineering,

v.35, n.2, p.185-194. 1990.

SENGER, R. S.; PAPOUTSAKIS, E. T. Genome-scale model for

Clostridium

acetobutylicum

. Part 1: metabolic network resolution and analysis. Biotechnology and

Bioengineering

. Publicado on line em 05/junho. 2008.

SHAMS, S. Y.; GONZALE, R. Anaerobic fermentation of glycerol: a path to economic

viability for the biofuels industry.

Current Opinion in Biothechnology, v.18, p.213-219.

2007.

SHINTO, H.; TASHIRO, Y.; YAMASHITA, M.; KOBAYASHI, G.; SEKIGUCHI, T.;

HANAI, T.; KURIYA, Y.; OKAMOTO, M.; SONOMOTO, K. Kinetic modeling and

sensitivity analysis of acetone-butanol-ethanol production.

Journal of Biotechnology, v.131,

n.1, p.45-56. 2007.

SINHA, S. Theoretical study of tryptophan operon - application in microbial technology.

Biotechnology and Bioengineering, v.31, n.2, p.117-124. 1988.

SIPSER, M.

Introduction to the theory of computation. Boston: PWS Publishing. 1996. xi,

239 p.

STELLING, J.; KLAMT, S.; BETTENBROCK, K.; SCHUSTER, S.; GILLES, E. D.

Metabolic network structure determines key aspects of functionality and regulation.

Nature,

v.420, n.6912, p.190-193. 2002.

STEPHANOPOULOS, G. N.; ARISTIDOU, A. A.; NIELSEN, J.

Metabolic engineering:

principles and methodologies

. San Diego: Academic Press. 1998. xvii, 707 p.

THEILGAARD, H. A.; NIELSEN, J. Metabolic control analysis of the penicillin biosynthetic

pathway: the influence of the LLD-ACV:bisACV ratio on the flux control.

Antonie van

Leeuwenhoek

, v.75, n.1, p.145-154. 1999.

TURNER, T. E.; SCHNELL, S.; BURRAGE, K. Stochastic approaches for modelling in vivo

reactions.

Computational Biology and Chemistry, v.28, n.3, p.165-178. 2004.

ULLAH, M.; SCHMIDT, H.; CHO, K. H.; WOLKENHAUER, O. Deterministic modelling

and stochastic simulation of biochemical pathways using matlab.

Systems Biology, IEE

Proceedings

, v.153, n.2, p.53-60. 2006.

URBANCZIK, R.; WAGNER, C. An improved algorithm for stoichiometric network

analysis: theory and applications.

Bioinformatics, v.21, n.7, p.1203-1210. 2005.

VARDAR-SCHARA, G.; MAEDA, T. Metabolically engineered bacteria for producing

hydrogen via fermentation.

Microbial Biotechnology, v.1, n.2, p.107-125. 2008.

VARMA, A.; PALSSON, B. O. Metabolic flux balancing: basic concepts, scientific and

practical use.

Nature Biotechnology, v.12, n.10, p.994-998. 1994a.

111

______. Stoichiometric flux balance models quantitatively predict growth and metabolic by-

product secretion in wild-type

Escherichia coli W3110. Applied Environment

Microbiology

, v.60, n.10, p.3724-3731. 1994b.

VASCONCELOS, I.; GIRBAL, L.; SOUCAILLE, P. Regulation of carbon and electron flow

Clostridium acetobutylicum grown in chemostat culture at neutral PH on mixtures of

glucose and glycerol.

Journal of Bacteriology, v.176, n.5, p.1443-1450. 1994.

VON KAMP, A.; SCHUSTER, S. Metatool 5.0: fast and flexible elementary modes analysis.

Bioinformatics, v.22, n.15, p.1930-1931. 2006.

XIU, Z. L.; CHANG, Z. Y.; ZENG, A. P. Nonlinear dynamics of regulation of bacterial

trp

operon: model analysis of integrated effects of repression, feedback inhibition, and

attenuation.

Biotechnology Progress, v.18, n.4, p.686-693. 2002.

XIU, Z. L.; ZENG, A. P.; DECKWER, W. D. Model analysis concerning the effects of

growth rate and intracellular tryptophan level on the stability and dynamics of tryptophan

biosynthesis in bacteria.

Journal of Biotechnology, v.58, n.2, p.125-140. 1997.

YANOFSKY, C.; PLATT, T.; CRAWFORD, I. P.; NICHOLS, B. P.; CHRISTIE, G. E.;

HOROWITZ, H.; VANCLEEMPUT, M.; WU, A. M. The complete nucleotide sequence of

the tryptophan operon of

Escherichia coli. Nucleic Acids Research, v.9, n.24, p.6647-6668.

1981.

YU, J.; XIAO, J.; REN, X.; LAO, K.; XIE, X. S. Probing gene expression in live cells, one

protein molecule at a time.

Science, v.311, n.5767, p.1600-1603. 2006.

Livros Grátis
( http://www.livrosgratis.com.br )
 
Milhares de Livros para Download:
 
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas

Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
 
 