Em formação

Modelo matemático cinético da glicólise?


Estou procurando artigos que contenham modelos matemáticos cinéticos de glicólise, espero que em células de mamíferos (ou o mais próximo possível, digamos levedura). Os documentos que encontrei fazem análises qualitativas (bifurcações, análises de estabilidade, etc.), mas não ajustam o modelo aos dados, então você não pode ter certeza da precisão das constantes cinéticas envolvidas. O trabalho que preciso fazer depende criticamente de ter estimativas razoavelmente precisas dessas constantes cinéticas.


Um modelo de núcleo cinético da rede de secreção de insulina estimulada por glicose de células β pancreáticas

É descrita a construção e caracterização de um modelo cinético central do sistema de secreção de insulina estimulada por glicose (GSIS) em células β pancreáticas. O modelo consiste em 44 reações enzimáticas, 59 variáveis ​​de estado metabólico e 272 parâmetros. Ele integra cinco subsistemas: glicólise, o ciclo do TCA, a cadeia respiratória, os ônibus NADH e o ciclo do piruvato. Também leva em consideração a compartimentação das reações no citoplasma e na matriz mitocondrial. O modelo mostra o comportamento esperado em suas saídas, incluindo a resposta da produção de ATP à concentração inicial de glicose e a indução de oscilações das concentrações de metabólitos na via glicolítica e nas concentrações de ATP e ADP. A identificação dos pontos de estrangulamento e a análise de sensibilidade dos parâmetros indicam que a via glicolítica e, em menor grau, o ciclo do TCA, são essenciais para o comportamento adequado do sistema, enquanto os parâmetros em outros componentes, como a cadeia respiratória, são menos críticos. Notavelmente, no entanto, a análise de sensibilidade identifica as primeiras reações das vias não glicolíticas como sendo importantes para o comportamento do sistema. O modelo é robusto para a deleção da atividade da enzima málica, que está ausente nas células β pancreáticas de camundongo. O modelo representa um passo em direção à construção de um modelo com parâmetros específicos da espécie que podem ser usados ​​para entender os modelos de diabetes de ratos e a relação desses modelos de ratos com o estado da doença humana.


Modelagem Cinética Estrutural

O comportamento temporal de uma rede metabólica, consistindo em m metabólitos e r reações, podem ser descritas por um conjunto de equações diferenciais (6), onde S denota o mvetor dimensional de reagentes bioquímicos e N a m × r matriz estequiométrica. o rvetor dimensional de taxas de reação ν(S, k) consiste em funções não lineares (e muitas vezes desconhecidas), que dependem das concentrações de substrato S, bem como em um conjunto de parâmetros k.

A seguir, não assumiremos o conhecimento explícito da forma funcional das equações de taxas, mas, em vez disso, objetivaremos uma representação paramétrica do Jacobiano do sistema. Como a única suposição matemática sobre o sistema, exigimos a existência de um estado positivo S 0 que cumpre a condição de estado estacionário (S 0 , k) = 0. Observe que o estado S 0 não precisa ser único ou estável.

Usando as definições (16, 17) com eu = 1, …, m e j = 1, …, r e aplicando a substituição de variável Seu = x eu S eu 0, o sistema pode ser reescrito em termos de novas variáveis x(t)

O Jacobiano correspondente do sistema normalizado no estado estacionário x 0 = 1 é porque as novas variáveis x estão relacionados à S por uma simples constante multiplicativa, Jx pode ser diretamente transformado de volta no jacobiano original.

Qualquer avaliação adicional do Jacobiano agora repousa na interpretação dos termos na Eq. 4 . Começamos com uma análise da matriz Λ: Seus elementos Λeu j têm as unidades de um tempo inverso e consistem nos elementos da matriz estequiométrica N, o vetor de concentrações em estado estacionário S 0, e os fluxos de estado estacionário ν(S 0). Desde que um sistema metabólico seja designado para modelagem matemática, podemos assumir que existe algum conhecimento sobre as concentrações relevantes, ou seja, para cada metabólito, podemos especificar um intervalo S eu − ≤ S eu 0 ≤ S eu +, que define uma faixa fisiologicamente viável da respectiva concentração. Além disso, os fluxos de estado estacionário ν(S 0) estão sujeitos à restrição de balanço de massa (S 0 ) = 0, deixando apenas r - classificação (N) taxas de reação independentes (6). Novamente, um intervalo νeu - ≤ νeu 0 ≤ νeu + pode ser especificado para todas as taxas de reação independentes, definindo um espaço de fluxo fisiologicamente admissível.

A seguir, denotamos S 0 e ν(S 0), geralmente correspondendo a um estado experimentalmente observado do sistema, como o “ponto operacional” em que o Jacobiano deve ser avaliado. Esta informação, juntamente com a matriz estequiométrica N, especifica totalmente a matriz Λ.

A interpretação da matriz θ x µ na Eq. 4 é um pouco mais sutil porque envolve os derivados das funções desconhecidas µ(x) com relação às novas variáveis ​​normalizadas no ponto x 0 = 1. No entanto, uma interpretação desses parâmetros é possível e não depende do conhecimento explícito da forma funcional detalhada das equações de taxas: Cada elemento θxeu µj da matriz θ x µ mede o grau normalizado de saturação da reação νj com respeito a um substrato S eu no ponto operacional S 0 . Em particular, a dependência de quase todas as leis de taxa bioquímica νj(S) em um reagente bioquímico S eu pode ser escrito na forma νj(S, k) = k v S eu n /f m(S, k), Onde n denota um expoente inteiro e f m(S, k) um polinômio de ordem m no S eu com coeficientes positivos k. Todos os outros reagentes foram absorvidos em k (6). Depois de aplicar a transformação da Eq. 2 , nós obtemos

com α ∈ [0, 1] denotando uma variável livre no intervalo unitário. Os casos limites são sempre limitadosSeu 0 → 0 α = 1 e limSeu 0 → ∞ α = 0. Para avaliar a matriz θ x µ , restringimos, portanto, cada parâmetro de saturação a um intervalo bem definido, especificado da seguinte forma: Como para a maioria das leis de taxas bioquímicas n = m = 1, a derivada parcial geralmente assume um valor entre zero e a unidade, determinando o grau de saturação da respectiva reação. No caso de comportamento cooperativo com expoentes n = m ≥ 1, a derivada parcial normalizada encontra-se no intervalo [0, n] e, analogamente, no intervalo [0, -m] para interação inibitória com n = 0 e m ≥ 1. Para exemplos e prova da Eq. 5 , veja, respectivamente, Materiais e métodos e a Apêndice de Apoio e as Figs. 8–11, que são publicados como informações de apoio no site da PNAS.

As matrizes θ x µ e Λ, conforme definido acima, especifique completamente o Jacobiano do sistema. A seguir, ambas as grandezas são tratadas como parâmetros livres, definindo o “espaço de parâmetros” fisiologicamente admissível do sistema. É importante ressaltar que nossa representação do Jacobiano preenche as três condições essenciais a seguir. (eu) O Jacobiano reconstruído representa o Jacobiano exato neste ponto no espaço de parâmetros. Não há aproximação envolvida. (ii) Cada termo no Jacobiano é diretamente acessível experimentalmente, como valores de fluxo ou concentração, ou tem uma interpretação bioquímica bem definida, como um grau normalizado de saturação de uma determinada reação. (iii) O Jacobiano não depende de uma escolha particular de funções de taxa específicas. Em vez disso, abrange todos os modelos cinéticos possíveis do sistema que são consistentes com as considerações acima. Nesse sentido, o Jacobiano reconstruído é exaustivo.


As equações básicas usadas para análise dinâmica

Os balanços dinâmicos de massa

A reconstrução de uma rede de reação em escala de genoma requer a identificação de todos os seus componentes químicos e as transformações químicas das quais eles participam. Este processo se baseia principalmente em genomas anotados e avaliação bibliômica detalhada (Reed et al, 2006). O resultado deste processo é a matriz estequiométrica, S, que é usado nos balanços dinâmicos de massa

que são a base de todos os modelos cinéticos. Aqui d (.) / Dt denota a derivada do tempo, x é o vetor das concentrações dos compostos na rede e v(x) é o vetor das taxas de reação. Todas as transformações bioquímicas são fundamentalmente uni ou bi-moleculares. Essas reações podem ser representadas pela cinética de ação em massa, ou generalizações das mesmas (Segel, 1975). A taxa de reação líquida para cada reação elementar em uma rede pode ser representada pela diferença entre um fluxo direto e reverso (por exemplo, consulte a Figura 1).

Esta formulação comumente usada é baseada em várias suposições bem conhecidas, como temperatura constante, volume e homogeneidade do meio. Se S, v (x) e as condições iniciais (x0) são conhecidos, então essas equações diferenciais ordinárias podem ser resolvidas numericamente para um conjunto de condições de interesse.

Forma linear

A caracterização dos estados dinâmicos de redes pode ser estudada por meio de simulação numérica ou por meio de análises matemáticas. Uma simulação depende do contexto e representa um estudo de caso. Os métodos matemáticos para a análise das características do modelo normalmente dependem do estudo das propriedades da transformação entre as concentrações e os fluxos. A análise de tais propriedades fundamentais normalmente se baseia na linearização das equações governantes em uma condição definida. A linearização das equações de balanço de massa dinâmico se resume à linearização do vetor de taxa de reação para formar a matriz de gradiente

e então formar a matriz Jacobiana em um estado de referência xref:

Onde x′=xxref e J é a conhecida matriz Jacobiana. A análise das características da matriz Jacobiana é um procedimento padrão na análise matemática da dinâmica do sistema (por exemplo, Strogatz, 1994). A aplicação desses métodos a redes bioquímicas vem sendo realizada há décadas (Heinrich et al, 1977, 1991 Heinrich e Sonntag, 1982 Palsson et al, 1987) e nos últimos anos tem havido um interesse renovado e, recentemente, novos desenvolvimentos na análise dinâmica das propriedades de J apareceram (Teusink et al, 2000 Kauffman et al, 2002 Famili et al, 2005 Bruggeman et al, 2006 Steuer et al, 2006 Grimbs et al, 2007 ).

A matriz Jacobiana para redes de reações bioquímicas é o produto de duas matrizes de dados. Antes de olhar para as propriedades fundamentais de J, consideramos os fluxos de trabalho e propriedades de dados que se relacionam com S e G.


2. Materiais e métodos

2.1. Estrutura do modelo e implementação numérica

Construímos um modelo cinético do metabolismo das células do câncer pancreático usando modelos publicados anteriormente de metabolismo de vários tipos de células (Mulukutla et al., 2010 Mar & # x000EDn-Hern & # x000E1ndez et al., 2011 Mulukutla et al., 2012 Mar & # x000EDn-Hern & # x000E1ndez et al., 2014 Shestov et al., 2014 Mulukutla et al., 2015). Nosso modelo é composto por um total de 46 metabólitos e 53 reações enzimáticas, incluindo glicólise, glutaminólise, o ciclo de TCA, o PPP e malato-aspartato-cetoglutarato-glutamato entre os compartimentos citosólico e mitocondrial (Figura 1). Cada etapa da via metabólica é modelada de acordo com reações enzimáticas conhecidas, que incluem mecanismos de reação que variam de Michaelis-Menten simples a uma cinética bi-bi aleatória complicada, expressa como diferentes formulações matemáticas. As leis de taxa para cada mecanismo de reação são incorporadas a um sistema de 46 equações diferenciais ordinárias não lineares (ODEs) que descrevem como as concentrações de metabólitos evoluem ao longo do tempo. Existe uma única ODE para cada metabólito, representando a taxa de variação da concentração da espécie, que depende das taxas nas quais a espécie é produzida e consumida na rede de reação. Usamos um solucionador de equações diferenciais implícitas no MATLAB (Guide, 1998) para integrar numericamente as equações e resolver as concentrações de metabólitos. Trata-se de um modelo determinístico, que simula as concentrações em um conjunto homogêneo de células que vivenciam, em média, condições ambientais intra e extracelulares semelhantes. Ao integrar os EDOs, calculamos a dinâmica média da população de células.

Figura 1. Esquema do modelo. A rede metabólica é composta por 46 metabólitos interagindo por meio de 53 reações enzimáticas. As principais vias envolvem a glicólise, a glutaminólise, o ciclo do TCA, o PPP e as reações de transporte entre os compartimentos mitocondrial (retângulo sombreado) e citoplasmático. As abreviaturas dos metabólitos, cofatores e enzimas são fornecidas no Arquivo Suplementar S3. Os nódulos coloridos representam os metabólitos para os quais a alteração de dobra foi medida experimentalmente durante o knockdown da enzima GOT1 (mostrado em vermelho). As setas representam a direção dos fluxos de reação no modelo de linha de base no início da simulação.

Resumimos brevemente as equações do modelo abaixo e o conjunto completo de ODEs é fornecido no Material Suplementar (arquivos de modelo: & # x0201CS1.m & # x0201D e & # x0201CS2.xml & # x0201D). Abreviações para os metabólitos e nomes de reações são fornecidas no Arquivo Suplementar S3 e os valores dos parâmetros fixos estão listados no Arquivo Suplementar S4. As equações de taxa detalhadas para glicólise e suas constantes cinéticas correspondentes são baseadas principalmente no modelo de glicólise para células HeLa (Mar & # x000EDn-Hern & # x000E1ndez et al., 2011, 2014). Esta rede de reação da glicólise foi estendida para incluir reações do ciclo de TCA e PPP usando leis de taxa cinética e parâmetros de Mulukutla e colegas de trabalho (Wu et al., 2007 Mulukutla et al., 2010, 2012, 2015). As reações que envolvem o consumo e a produção de ATP no citoplasma foram definidas no modelo de Shestov et al. (2014), e as concentrações de ATP e ADP no compartimento mitocondrial foram mantidas constantes como em Mulukutla et al. Além disso, os parâmetros de transporte de glutamina foram obtidos de Pingitore et al. (2013).

AKT é um forte promotor da tumorigenicidade do câncer pancreático mediado por KRAS (Asano et al., 2004) devido à sua influência nas taxas de reações metabólicas na glicólise. Sabe-se que as células PDAC têm aumento da captação de glicose (Ying et al., 2012), que é mediada pela suprarregulação de enzimas glicolíticas específicas, incluindo o transportador de glicose-1 (GLUT1), hexoquinase (HK) e lactato desidrogenase A (LDHA ) Além disso, AKT promove aumento da captação de glicose ao ativar GLUT1, HK e a enzima fosfofrutocinase (PFK) (Rathmell et al., 2003 Elstrom et al., 2004). Incorporamos o efeito de AKT em nosso modelo metabólico, simulando a atividade glicolítica aumentada mediada por AKT. Especificamente, as atividades das enzimas GLUT1, HK e PFK (representadas por seus Vmax valores) são modelados para ter 20% de atividade basal, enquanto 80% de sua atividade é devido à ativação por AKT (Mosca et al., 2012 Mulukutla et al., 2012).

A fim de prever como as vias metabólicas intracelulares influenciam o crescimento celular, incorporamos o número de células às reações catalisadas por enzimas. Especificamente, o modelo é aumentado para incluir uma 47ª ODE que descreve a evolução temporal do número de células cancerosas, CN. O crescimento celular é implementado como uma equação logística (Enderling e Chaplain, 2014) que é responsável pela capacidade máxima de carga do tumor, KCC (Equação 1).

O número de células cancerosas está diretamente ligado ao metabolismo, onde a taxa de crescimento depende das concentrações intracelulares de três metabólitos primários conhecidos por influenciar a proliferação celular: glicose, glutamina e ATP (Venkatasubramanian et al., 2008 Zhu et al., 2012) . A dependência desses três metabólitos é modelada assumindo funções do tipo Monod (Higuera et al., 2009) (Equação 2).

Os parâmetros de crescimento e morte & # x003B1atp, & # x003B1glc, & # x003B1gln, e & # x003B1d estão em unidades de min & # x022121. Os parâmetros de concentração kgc, kgn, e kap têm unidades de mM.

2.2. Condições iniciais

Simulamos o modelo com vários conjuntos de concentrações de metabólitos iniciais para identificar a faixa apropriada de condições iniciais. As informações sobre as concentrações iniciais de metabólitos intracelulares em células cancerosas pancreáticas são limitadas. Portanto, permitimos que a concentração inicial de cada metabólito variasse dentro de uma faixa especificada. Especificamos os intervalos de concentração com base em medições publicadas obtidas a partir de várias linhas celulares, incluindo câncer cervical humano (Mar & # x000EDn-Hern & # x000E1ndez et al., 2011, 2014), amostras de tecido normal doente e circundante de pacientes com câncer de estômago e cólon (Hirayama et al., 2009), extratos de células de câncer de mama (Le Guennec et al., 2012), amostras de pacientes com câncer PDAC (Fontana et al., 2016) e mieloma de camundongo e linhas celulares CHO (Mulukutla et al., 2012, 2015) . Incerteza adicional para células cancerosas pancreáticas foi considerada aumentando e diminuindo os limites superior e inferior, respectivamente, em 20%. Devido à falta de medições que distinguem os níveis de metabólitos em diferentes compartimentos celulares, as concentrações iniciais de metabólitos que estavam presentes nos compartimentos mitocondrial e citosólico foram consideradas as mesmas. Os intervalos de concentrações de metabólitos apresentados na Tabela 1 são responsáveis ​​pela variabilidade nas medições da literatura, bem como uma incerteza adicional para a concentração intracelular desconhecida de linhas de células de câncer pancreático em particular.

Tabela 1. Limites para as condições iniciais usadas nas simulações do modelo.

O Latin Hypercube Sampling (McKay et al., 2000 Oguz et al., 2013) foi aplicado à amostra dentro dos intervalos selecionados para cada metabólito. O LHS separa a faixa de concentrações dos metabólitos em vários intervalos e amostras de cada intervalo exatamente uma vez, explorando assim com eficiência toda a gama possível de condições iniciais para cada metabólito. Selecionamos para obter 100 conjuntos de condições iniciais para cada metabólito para análise de identificabilidade de parâmetro (Seção 3.1.1) e, em seguida, selecionamos aleatoriamente 50 desses conjuntos para serem usados ​​na estimativa de parâmetro (Seção 3.1.3). Este procedimento explora adequadamente os possíveis intervalos de condições iniciais enquanto equilibra os recursos computacionais necessários para a otimização de parâmetros globais.

2.3. Estimativa de Parâmetro

O modelo de linha de base, adaptado da literatura, tem um total de 372 parâmetros, que inclui 71 velocidades de reação (as taxas de avanço e reverso, Vf e Vr, respectivamente). As velocidades de reação refletem a quantidade de enzima presente e a atividade enzimática correspondente. Convencionalmente, as velocidades de reação são pensadas para distinguir o metabolismo em diferentes tipos de células. Portanto, dos muitos parâmetros cinéticos incluídos nas equações de taxa de reação, apenas as velocidades de reação foram ajustadas aos dados de treinamento, e as outras constantes de taxa foram mantidas em seus valores de literatura. Também ajustamos os parâmetros de crescimento celular mostrados nas Equações (1) e (2). A seguir, descrevemos os dados experimentais usados ​​para treinar o modelo e o método usado para realizar a estimativa dos parâmetros.

O modelo é adequado para medições experimentais de Son et al. (2013), que mediu as concentrações de 14 metabólitos intracelulares usando análise metabolômica direcionada. Filho e colegas de trabalho procuraram entender o metabolismo não canônico da glutamina nas células do câncer pancreático após o knockdown do GOT1, uma enzima importante no metabolismo da glutamina. As concentrações de metabólitos foram medidas quando a enzima GOT1 foi desativada, em relação à condição sem desativação. Assim, eles quantificaram a alteração nas concentrações de metabólitos.

O protocolo experimental usado por Filho e colegas de trabalho é ilustrado na Figura S1. Simulamos a mesma sequência de etapas para prever a alteração nas concentrações dos 14 metabólitos. Uma vez que o nível de expressão relativa da enzima pode ser correlacionado com a regulação dos níveis de atividade da enzima, simulamos o knockdown da enzima multiplicando o Vf pelo fator (1 - & # x003B1) (Nolan e Lee, 2012). Consideramos o valor de & # x003B1 como 0,85, com base no nível de expressão GOT1 médio de dois experimentos de knockdown relatados em Son et al. (2013). O modelo é simulado para GOT1 knockdown para prever a alteração na concentração dos 14 metabólitos em relação ao caso sem knockdown. Procuramos minimizar a soma ponderada do erro quadrático (WSSR) entre os dados experimentais e as previsões do modelo.

Além disso, Son e colegas usam em vitro cultura de células para investigar como o metabolismo intracelular influencia a proliferação celular. Eles medem o número de células com e sem knockdown de GOT1 e na presença de concentrações de nutrientes extracelulares variáveis. Também simulamos seus protocolos experimentais e comparamos as previsões do modelo com suas medidas experimentais.

A otimização do enxame de partículas (PSO) foi usada para identificar os valores dos parâmetros necessários para permitir que as previsões do modelo se ajustem melhor aos dados e minimizem o WSSR. PSO (Iadevaia et al., 2010 Kennedy, 2010 Tashkova et al., 2011) é uma técnica de otimização global estocástica de inspiração biológica desenvolvida por Kennedy e Eberhart (1995). Baseia-se no conceito de comportamento social observado na natureza. Em PSO, muitos partículas, conjuntos de parâmetros são constantemente atualizados a partir de seus valores iniciais aleatórios para identificar os valores dos parâmetros que melhor se ajustam aos dados experimentais. Cada partícula tem uma posição, representando a localização no espaço de parâmetros multidimensional, e uma velocidade com a qual ela se move em direção a um mínimo local no WSSR. As partículas se comunicam entre si para atualizar sua posição e velocidade, em última análise, movendo-se em direção ao mínimo global no WSSR. Usamos PSO para estimar as velocidades de reação para o modelo de linha de base. Cada partícula representa um vetor de todas as velocidades de reação a serem otimizadas, onde os valores dos parâmetros iniciais são retirados de um espaço bem amostrado com os limites dados. Para especificar os limites, as velocidades de reação foram permitidas variar 100 vezes para cima e para baixo de seus valores iniciais (retirados da literatura, consulte Materiais e Métodos). Cada execução do algoritmo PSO executa 2.500 iterações, um valor definido pelo usuário para equilibrar a despesa computacional da pesquisa de parâmetro. Realizamos a estimativa do parâmetro duas vezes para cada conjunto de condições iniciais (ou seja, um total de 5.000 iterações por condição inicial) e, para cada caso, selecionamos o conjunto de parâmetros que gerou o menor erro. Isso deu um total de 50 conjuntos de parâmetros de melhor ajuste, um conjunto para cada condição inicial.

Estimar as velocidades de reação para cada condição inicial foi a primeira etapa do ajuste do modelo. Na segunda etapa do ajuste do modelo, buscou-se estimar os parâmetros de crescimento minimizando o WSSR. Como há menos parâmetros para ajustar em comparação com a primeira etapa de ajuste, usamos a otimização de mínimos quadrados não linear. Realizamos o ajuste 100 vezes para cada condição inicial para se aproximar do mínimo global no erro do modelo. Dado o conhecimento prévio limitado da faixa de valores de base para os parâmetros de crescimento (Higuera et al., 2009), pesquisamos em um espaço de parâmetro abrangendo sete ordens de magnitude para cada parâmetro. As simulações do modelo para otimizar o crescimento celular foram conduzidas de modo que o mesmo conjunto de sete parâmetros de crescimento pudesse se ajustar à curva de crescimento experimental para ambas as condições sem knockdown e GOT1 knockdown.

2.4. Extração de dados

Dados experimentais para treinamento e validação de modelo foram extraídos de Son et al. (2013) usando o programa MATLAB GRABIT (Guia, 1998). Os dados de treinamento incluem a variação nas concentrações de metabólitos e o número de células sob o knockdown de GOT1. Os dados de validação incluem o número de células em privação de nutrientes.

2,5. Análise de Identificabilidade de Parâmetro

Usamos a análise de identificabilidade dos parâmetros estruturais (Maly e Petzold, 1996 Ascher e Petzold, 1998 Shampine et al., 1999 Finley et al., 2011 Berthoumieux et al., 2013) para reduzir o número de parâmetros do modelo adequados aos dados de treinamento. A identificabilidade do parâmetro determina dependências implícitas entre os parâmetros. Se dois parâmetros forem encontrados correlacionados, podemos especificar uma relação matemática entre os parâmetros e apenas ajustar um no procedimento de estimativa de parâmetro. Aqui, nós apenas especificamos a relação entre as velocidades de reação direta e reversa correlacionadas, onde a velocidade de reação reversa, Vr, é expresso como uma função da velocidade de reação direta, Vf, com a constante de equilíbrio, Veq: Vr = Vf/Veq. Nestes casos, apenas a velocidade de reação direta é adequada aos dados experimentais, reduzindo assim o número de parâmetros ajustados. o Veq é calculado usando os trabalhos publicados dos quais nosso modelo é derivado (Wu et al., 2007 Mulukutla et al., 2010, 2012, 2015 Mar & # x000EDn-Hern & # x000E1ndez et al., 2011, 2014).

2.6. Análise sensitiva

Aplicamos a análise de sensibilidade global (Saltelli et al., 2008) para determinar quais dos parâmetros do modelo influenciam mais significativamente as concentrações de metabólitos previstas. Especificamente, usamos o método estendido de Fourier Amplitude Sensitivity Test (eFAST) (Marino et al., 2008), uma abordagem baseada em variância, para entender a robustez dos resultados do modelo (concentrações de metabólitos) dada a variância nas entradas do modelo (a reação velocidades) (Zi, 2011). Permitimos que as entradas do modelo variassem duas ordens de magnitude para cima e para baixo em relação aos valores da literatura. O método eFAST calcula dois índices que fornecem uma estimativa da sensibilidade dos resultados do modelo em relação aos parâmetros do modelo. O índice de primeira ordem, Seu, quantifica a variação da saída do modelo em relação às variações de cada entrada individual e o índice FAST total, Sti, quantifica a variação da saída do modelo em relação às variações de cada entrada e covariâncias entre combinações de entradas. o Seu, então, é uma medida da sensibilidade local da saída do modelo para cada entrada individual, enquanto Sti é uma medida da sensibilidade global, contabilizando as interações ou correlações entre múltiplas entradas.


Discussão

Considerações gerais

Neste trabalho, usamos um modelo cinético detalhado do metabolismo da glicose hepática para investigar a contribuição do fígado para a homeostase da glicose sanguínea em condições fisiológicas e patológicas. Validamos o modelo comparando simulações de modelo com uma variedade de dados experimentais sobre fluxos de troca de glicose, concentrações de metabólitos e enchimento / esvaziamento do estoque de glicogênio intra-hepático. Nosso objetivo central foi dissecar a importância relativa dos mecanismos reguladores de enzimas individuais para a resposta adequada do fígado a níveis variáveis ​​de glicose no plasma. O achado geral desta análise é que as mudanças no estado de fosforilação induzida por hormônio de enzimas regulatórias importantes, bem como as mudanças na concentração de efetores alostéricos, são pelo menos da mesma importância que as mudanças na abundância da enzima para o ajuste da produção metabólica à demanda metabólica definida pelas condições externas. A forte influência dos reguladores além das mudanças na abundância de proteínas é presumivelmente a principal razão para a fraca correlação geralmente observada entre os fluxos metabólicos e a abundância das enzimas catalisadoras associadas.

A necessidade de uma ação coordenada de diferentes modos de regulação enzimática pode ser fundamentada considerando que mesmo desafios extremamente elevados para a adaptação do metabolismo celular (por exemplo, durante períodos de fome ou supernutrição, exposição a agentes tóxicos, inflamação ou proliferação) nunca são constantes com o tempo, mas flutuando e, portanto, não pode ser satisfeita apenas aumentando ou diminuindo a abundância de enzimas. Portanto, é justo afirmar que qualquer conceito teórico que vise a um melhor entendimento da regulação das redes metabólicas deve levar em consideração a regulação das atividades enzimáticas além do nível de expressão gênica.

Consequências funcionais das mudanças na abundância de proteínas durante a transição entre os estados nutricionais em jejum e alimentado e no diabetes

Em primeiro lugar, analisamos como as mudanças na abundância das principais enzimas metabólicas relatadas para o fígado do rato em jejum e alimentação em condições nutricionais influenciam a produção metabólica em várias condições fisiológicas. Alterações na abundância das principais enzimas regulatórias das vias glicolíticas e gliconeogênicas foram encontradas para acarretar diferenças significativas na relação entre os níveis de glicose plasmática e os fluxos de troca de glicose hepática - o fígado "em jejum" torna-se um produtor de glicose mais forte, enquanto o "alimentado" o fígado se torna um usuário mais forte da glicose. Esta adaptação é alcançada por mudanças na abundância enzimática e é vantajosa para a homeostase da glicose plasmática enquanto a situação fisiológica prevista persistir. No entanto, pode ser uma desvantagem se ocorrer uma mudança repentina (inesperada). Um fígado adaptado ao jejum de 1-2 dias está menos preparado para responder a um forte aumento súbito da glicose plasmática do que um fígado que apresenta níveis de glicose plasmática continuamente elevados (Fig. 15). Uma complicação clínica bem conhecida que ocorre durante a realimentação em pacientes fortemente desnutridos é a intolerância à glicose [12], um tipo de desregulação metabólica que normalmente ocorre nos estágios iniciais do diabetes tipo 2. De acordo com esta observação, a relação entre os níveis de glicose plasmática e hepática as taxas de troca de glicose previstas pelo modelo são muito semelhantes em condições de jejum e diabetes (Fig. 10), ou seja, a gliconeogênese é aumentada e a capacidade glicolítica reduzida.

Produção e utilização diurna de glicose pelo fígado

As taxas de HGP e HGU dependem do nível plasmático de nutrientes e hormônios que mudam permanentemente ao longo do dia. Tomando perfis de plasma diurnos medidos como entrada, o modelo cinético permite a simulação de mudanças diurnas de HGP e HGU e o estado de enchimento do glicogênio (Figs. 12, 13 e 14). Essas simulações mostram que, no estado nutricional alimentado, o fígado é capaz de alternar entre HGP e HGU. Durante o dia, o fígado funciona predominantemente como um utilizador de glicose e durante a noite é um produtor de glicose. Dependendo do momento da ingestão de alimentos e da duração e intensidade do exercício físico, os perfis metabólicos individuais podem diferir significativamente do perfil genérico usado como entrada para nossas simulações. Em contraste com o estado alimentado, durante a inanição de longo prazo, prevê-se que o fígado funcione constantemente como produtor de glicose para garantir que os níveis de glicose no plasma permaneçam suficientemente altos para abastecer os consumidores de glicose, como o cérebro e os eritrócitos. No caso do diabético, as alterações na abundância das enzimas glicolíticas e gliconeogenéticas, juntamente com as respostas prejudicadas do hormônio da glicose, levam a uma mudança patológica em direção à gliconeogênese (Fig. 11). A boa concordância entre a mudança do ponto de ajuste e os níveis de glicose plasmática observados reforça a importância do fígado na determinação dos níveis de glicose plasmática. Essas diferenças fundamentais nos estados metabólicos basais do fígado também se refletem nos estados de enchimento do glicogênio. No estado nutricional alimentado, o glicogênio é degradado durante o HGP e reabastecido durante o HGU, o enchimento do estoque de glicogênio varia entre 50% e 80% da capacidade total de armazenamento. No estado de jejum, o estado de enchimento geral do estoque de glicogênio é muito menor e varia apenas entre 10% e 25%.

Avaliação da importância relativa da abundância enzimática variável e regulação cinética das atividades enzimáticas para a regulação das taxas de troca de glicose hepática

To investigate the relative importance of the different regulation modes of enzyme activities we simulated the glucose exchange flux of the liver for different nutritional states (Fig. 16 and Table 2). In the fed state, the strongest regulatory influence is exerted by the changes of enzyme abundances, whereas in the fasted state, reversible phosphorylation has the largest impact. An important finding is that the short-term metabolic adaptation of the liver can be largely attributed to hormonal regulation as the glucose exchange fluxes become almost constant when we fix the phosphorylation state of the interconvertible enzymes. Nevertheless, the impact of allosteric regulation is substantial both in the fasted and the fed state, accounting for approximately 50 % of flux changes brought about by reversible phosphorylation in the fasted and fed state, respectively. It has to be noted that the role of allosteric regulation is certainly underestimated in our model as the concentration values of important cofactors (adenosine tri-/di-/monophosphate, nicotinamide adenine dinucleotide, and its reduced form NADH) and of allosterically important metabolites of the citric acid cycle (e.g. citrate inhibition of the phosphofructokinase) were not taken into account.

Metabolic control analysis (MCA)

To dissect the importance of individual enzymes for hepatic glucose exchange rates under different conditions we used the MCA concept. Calculation of flux control coefficients for the ‘fed’ and ‘fasted’ states revealed that enzymes carrying significant control are those showing significant changes of their abundance under different physiological conditions. Furthermore, the flux control is shared between different groups of enzymes in different conditions – enzymes being important in the glycolytic phase of liver metabolism are different from the ones central during gluconeogenesis (Table 3). Importantly, the control coefficients for the glucose exchange flux exhibit significant fluctuation over one day and diverge (by definition) when the glucose exchange flux is zero (Fig. 17).

We also calculated π-elasticity coefficients to quantify the relative share of reversible phosphorylation and concentration changes of reactants and allosteric effectors in the regulation of individual enzymes of hepatic glucose metabolism (Fig. 18). This analysis revealed a large variability in the relative contribution of the three fast regulatory modes to the control of regulatory enzymes and hence the control of glucose exchange flux.

Direct experimental validation of the computed elasticities in vivo is unfeasible because this would require monitoring of the glucose exchange flux of the liver at clamped plasma levels of glucose and hormones in response to the gradual variation of an effector specifically influencing one kinetic parameter of the target enzyme under study. As a surrogate, we checked whether the predicted elasticities are concordant with observed changes in plasma glucose levels induced by targeting a single key regulatory enzyme either by drugs or genetic interventions.

The maximal control that can be exerted by GK is low in the fasted-hypoglycemic state but becomes large in the fed-hyperglycemic state. Experimentally, glucosamine-induced inhibition of GK caused only marginal reduction of glucose uptake in euglycemia, whereas in hyperglycemia a significant reduction of the net hepatic glucose uptake of about 40 % was observed [13], confirming our simulation results. Clinically, the regulatory importance of GK in hyperglycemia is used to target this enzyme in diabetes type 2 [14, 15].

Glycogen phosphorylase (GP)

Torres et al. [16] investigated the effects of a GP inhibitor (GPeu) and metformin on hepatic glucose in presence of basal and four-fold increased levels of plasma glucagon in 18-h fasted conscious dogs. In euglycemic conditions, no change in the net hepatic glucose balance and plasma glucose was observed in the presence of GPeu. However, after glucagon stimulation, the presence of GPeu significantly diminished the glucose output. Both findings confirm the predicted relatively high control of the GP in hypoglycemic conditions as well as the large share of allosteric regulation of this enzyme.

Phosphofructokinase-1 (PFK1)

To our knowledge, a rate-limiting role of this enzyme in the liver is not reported. In several non-hepatic tissues, PFK1 exerts only insignificant control of glycolytic flux [17], which agrees with the predicted very small values of control coefficients in both the hypo- and hyperglycemic cases.

PFK2/FBP2

The enzyme PFK2/FBP2 exerts control of the glucose exchange flux mainly by changes in its phosphorylation state, whereas changes in the abundance of this enzyme have no impact of the glucose exchange flux. The reason for the latter is the bifuncionality of this enzyme – the phosphorylated enzyme (PFK2) acts as a kinase catalyzing the formation of fructose 2,6-bisphosphate (Fru26P2), an efficient allosteric activator of PFK1, whereas the non-phosphorylated enzyme (FBP2) acts as a phosphatase catalyzing the degradation of Fru26P2 to fructose 6-phosphate (Fru6P). These two opposite reactions create a futile cycle Fru6P → Fru26P2 → Fru6P that consumes one molecule of ATP. Obviously, unequal modulation of these opposite activities cannot be achieved by changes of protein abundance because any change in enzyme amount influences both activities to the same extent and thus leaves the net flux unchanged. However, reversible phosphorylation enhances FBP2-activity and in parallel diminishes the PFK2-activity of this enzyme, resulting in a high sensitivity of the glucose exchange flux to changes in the phosphorylation state of the PFK2/FBP2.

Phosphoenolpyruvate carboxykinase (PEPCK)

Metabolic control of liver gluconeogenesis was quantified in groups of mice with varying PEPCK protein content. Surprisingly, livers with a 90 % reduction in PEPCK content showed only a 40 % reduction in gluconeogenic flux, indicating a lower than expected capacity for PEPCK protein content to control gluconeogenesis (estimated control coefficient of about 0.18) [18]. This is in good agreement with our theoretical predictions. She et al. [18] concluded that the liver PEPCK functions more as an integrator of hepatic energy metabolism than as a determinant of gluconeogenesis.

Knowledge of the flux control exerted by a specific enzyme and of the regulatory mechanisms that contribute to its control is valuable information for the design of new drugs. For example, our analysis revealed that changing the protein abundance of the bifunctional enzyme PFK2/FBP2 should have no influence on the stationary glucose exchange flux of hepatocytes. Hence, drugs targeting this enzyme as non-competitive inhibitors can be expected to have little impact on the modulation of the hepatic glucose exchange flux. However, our analysis suggests that drugs specifically targeting only the phosphorylated enzyme (phosphatase) or non-phosphorylated enzyme (kinase) have a strong impact on the glucose exchange flux. These theoretical findings are supported by the fact that cancer cells express a specific phosphatase (TIGAR) that catalyzes the degradation of the glycolytic activator Fru26P2 in order to suppress glycolysis and to redirect the glucose flux through the oxidative pentose phosphate pathway. Flux control by PFK2/FBP2 may serve as a good example of why up or down regulation of the abundance of an enzyme does not necessarily imply corresponding flux changes as is assumed in numerous publications dealing with the potential metabolic consequences of varying protein abundances.


Fundo

A kinetic model of a biochemical reaction network consists of a set of ordinary differential equations describing the dynamics of the concentrations of all metabolites in the reaction network. Most biochemical reaction networks are complex and involve many enzyme-catalyzed processes with non-linear kinetics and intricate stoichiometric and regulatory interactions between the enzymes. Consequently, the mathematical models of such networks contain high-dimensional sets of coupled rational differential equations, which sometimes require huge computational effort to analyze. The current state-of-the-art numerical tools for stability analysis, for bifurcation study and for other types of dynamical analysis are known to suffer from a so-called curse-of-dimensionality. For example, the largest biological model that has numerically been analyzed for bifurcation in[1] consists of 25 metabolites and 37 parameters and the one in[2] has 22 metabolites. Moreover since complex models of biochemical reaction networks involve a huge number of parameters, the task of identifying these parameters (in addition to those parameters that have been identified in the literature) is enormous and requires large datasets. The complexity of this task is further compounded by the fact that often not all the metabolite concentrations can be measured. Thus, there is a need for techniques that can reduce a given kinetic model of a biochemical reaction network to a simplified version that mimics the behaviour of the original model satisfactorily, but contains less differential equations and parameters.

For biochemical reaction networks, a number of model reduction techniques are known. See[3] for a detailed review of some of the well known methods of model reduction. Here we list only a few of the known methods in the literature. The singular perturbation method, the time-scale separation technique[4–8], the rapid-equilibrium approximation, also known as the quasi-equilibrium approximation (see[9]) and the quasi steady-state approximation (QSSA) (see for e.g.,[10]) are the most commonly used techniques. The reduced state vector obtained by any of these techniques contains a subset of the metabolite concentrations (state vector components) of the full model. Härdin[11] extends the QSSA approach by considering the higher-order approximation in the computation of the quasi steady-state. In[12], some approaches for reduction of complexity in biochemical reaction networks are considered for use in SYCAMORE, which is a computational research environment in systems biology. One of the approaches considered in[12] is the intrinsic low-dimensional manifold (ILDM) approach, which is an improved time-scale separation technique. More recently[13], a computational singular perturbation (CSP) algorithm was developed to analyze the time-scales of the NF- κ B signaling network which could be used in order to reduce its model. The application of singular perturbation, time-scale separation, quasi equilibrium and quasi steady state approaches to general enzyme-kinetic rate laws, such as Michaelis-Menten and ping-pong bi-bi is difficult and leads to complicated rate law expressions in the reduced models. Some of the time-scale separation techniques are either based on a priori experimental information or eigenvalues of the Jacobian corresponding to the model and hence are only locally effective. Another recent approach for model reduction uses tropical geometry (see e.g.,[14]), wherein the polynomial occuring in every rate equation is replaced by a monomial which is equal to the largest, in absolute value among the monomials that constitute the polynomial.

One of the ways of reducing the complexity of a biochemical model is to reduce the number of parameters in it. This can be done by carrying out a parameter sensitivity analysis (see for e.g.,[15]) and eliminating those parameters whose variations have least effect on the dynamics of a given network. In[8], a time-scale analysis is first done using experimental data to identify the fast and the slow metabolites and this information is then used to carry out an a priori parameter sensitivity analysis to obtain a reduced kinetic model of a biochemical reaction network. In[16], an implicit multiparametric variability analysis (MPVA) method is used to search the entire parameter space in order to determine the set of parameters that can be eliminated. In[17], the authors go one step further by identifying a region in the parameter space where some of the parameters are zero-valued, others have readjusted values and where nevertheless the outputs of the original model match those of the reduced model within a certain tolerance. They further show that parameter sensitivity analysis approach for model reduction may not be always successful. Another way of reducing the complexity of a biochemical reaction network model is to reduce the number of reactions. In[18–20], optimization techniques are used to determine which reactions can be eliminated from the original model without a substantial alteration of the model behaviour.

The method proposed in[21] simplifies a given chemical reaction network by linearly combining reactions in such a way that the resulting network has lesser number of species. However, the kinetics for the reduced set of reactions are determined by the rate limiting steps in the original network and this requires prior biological knowledge of the network. Danø et al.[22] propose reduction by identification and elimination of variables in such a way that the basic dynamic properties of the original model are preserved. In[23], model reduction is achieved by simplifying the rate equations of individual enzymes in the network. A limitation of this method is that it yields a reduced model that predicts measurement data only under specific experimental conditions.

In this paper, we describe a new model reduction method that reduces the number of reactions, metabolites and parameters, such that the dynamics of the metabolite concentrations of the reduced model are close to those of the original model. This method proceeds by a simple stepwise reduction in the number of ‘complexes’, which are defined as the left and right-hand sides of the reactions in the network. The effect of this stepwise reduction is monitored by an error integral, which quantifies how much the behaviour of the reduced model deviates from the original. The method is based on the reduction of the underlying weighted Laplacian (see[24] for a definition) describing the graph structure of the complexes and reactions in the chemical reaction network. A similar technique is also employed in the Kron reduction method for reduction of resistive electrical network models[25].

Our model reduction technique is easy to implement and can be used to reduce reaction networks governed by a vast majority of enzyme kinetic rate laws including Hill kinetics and reversible Michaelis-Menten kinetics. It does not rely on prior knowledge about the dynamic behaviour or biological function of the network. Consequently, it can be automated. Furthermore, the reduced model largely retains the kinetics and structure of the original model. This enables a direct biochemical interpretation and yields insight into which parts of the network have the highest influence on its behaviour. It also accelerates computations and facilitates parameter fitting, especially when we deal with models of huge biochemical reaction networks.

We show the application of our model reduction technique to a yeast glycolysis model[26] and a model of beta-oxidation in rat liver[27]. We have simulated the transient behaviour of the metabolites that are not eliminated during the model reduction procedure. In both the cases, a 30% reduction of the number of variables still retained about 92.5–96.5% agreement between the outputs of the full and the reduced networks.


Organism-level constraints

Organism-level constraints are based on properties that are unique to a specific organism while being consistent for all environmental or experimental conditions. The metabolic network of an organism (determined by DNA sequence) itself can be seen as an organism-level constraint. Organism-level constraints are mostly based on knowledge about physiological limitations and peculiarities of a particular organism or the assumption that the modified organism design is feasible if resources and/or parameters of the existing organism will not be exceeded. No information about experimental conditions is needed to determine values of constraints, as they are applicable for all experimental conditions.

To take into account limited enzyme-building resources, the total enzyme activity constraint [19,20] can be implemented by setting limits for the sum of enzyme concentrations without detailed experiment-specific analysis. This type of constraint is based on assumption that the modified organism should be able to produce as much (or a small fraction more) protein as the initial one. The total enzyme activity constraint is used in several kinetic model studies [15,21–24]. The total enzyme activity constraint is implemented also in stoichiometric models [25,26].

The application of steady-state constraints in both kinetic and stoichiometric models enables the synergy between two model types. Most kinetic models are relatively small but include kinetic parameters and take into account metabolite concentrations. The steady-state fluxes found in kinetic models can be put into stoichiometric models as constraints to test the feasibility of the kinetic model steady-state flux distribution at genome scale where mass and energy balance can be taken into account at a larger scale [27]. The range of metabolite concentrations in kinetic models can be also used to calculate lower and upper constraints of reaction fluxes in stoichiometric models where concentrations cannot be directly applied.

The ability of kinetic models to calculate metabolite concentrations enable the application of metabolite concentration-related constraints. In most cases, these are organism-specific and not related to particular experimental conditions. Each organism may have specific metabolites that are cytotoxic above some concentration that can be used as an upper limit constraint for metabolite concentrations [28]. Similarly, there can be unrealistic or unfeasible upper or lower levels of metabolite concentration prohibited as new steady-state concentrations or even as concentrations during the transition process [29].

Another metabolite concentration-related constraint is the homeostatic one [30,31] used to limit the impact of large changes of internal metabolite concentration in pathway-scale kinetic models to other reactions outside the model's scope via gene expression, reaction flux, reaction directionality, and other potential effects [15,28]. The homeostatic constraint limits optimized steady-state concentrations of metabolites within some range around the steady-state concentrations of the initial model. It can be applied for (1) a pool of internal metabolites [15,21,24,32], (2) each metabolite separately [23,29,33], or (3) a combination of both [22]. Heavy reduction in objective function values by the homeostatic constraint can take place [23]. It might be useful to choose an acceptable concentration range for each metabolite separately taking into account their degree of involvement in other cellular processes outside of the model's scope [23]. A variety of homeostatic constraints is used by Magnus et al. [15] applying metabolite concentration-dependent Gibbs-free energy calculations to assure feasibility of reactions.

A minimal set of adjustable parameters is used as a design constraint, assuming that it reduces the number of unpredictable potential side effects not included in the model [32–34]. It can be applied both for kinetic and stoichiometric models. In case of kinetic models, the comparison of optimization results for a different number of adjustable parameters may be required [35]. The optimization runs for many adjustable parameter combinations and can be done also in an automated way [36].

Application example of organism-level constraints

The impact of the total enzyme activity constraint and the homeostatic constraint looking for minimal set of adjustable parameters (total optimization potential (TOP) approach) is demonstrated in detail [23] (Figure 2), optimizing the model to increase sucrose accumulation in sugarcane culm [37]. The aim is to maximize objective function: proportion of sucrose accumulation in the vacuole relative to sucrose hydrolysis by invertase [23]. The best objective function value for a combination of all five adjustable parameters without constraints was 2.6 × 10 6 (Figure 2B), requesting an unrealistic 1500-fold increase in glucose concentration and a 5-fold increase in enzyme concentrations that are used as adjustable parameters. Introduction of the total enzyme activity constraint did not allow increase in total concentration of enzymes and reduced the objective function value 10-fold to 0.16 × 10 6 (Figure 2B), still relying on an unrealistic 118-fold increase in fructose concentration. Introducing the homeostatic constraint allowing changes of metabolite concentrations just by ±20% reduced the objective function to 4.7 (Figure 1C), demonstrating a dramatic decrease in the objective function value compared with previous cases. Despite that it brings an 34% increase in the objective function value of the original model. Implementation of both the total enzyme activity and homeostatic constraints (Figure 2D) did not reduce further the objective function value of the full set of five adjustable parameters (concentrations of five enzymes) [23]. This illustrates how highly promising but unrealistic designs are made less promising by the value of objective function but are improved in terms of biological viability by implementing additional constraints. The same study demonstrated also application of the TOP approach, revealing that homeostatic and total enzyme activity constraints heavily influence the rank of best adjustable parameter combinations when a minimal set of adjustable parameters is sought (Figure 2). Applying the homeostatic constraint the objective function increase is smaller, but almost all optimization potential can be reached by manipulating values of just two adjustable parameters (Figure 2C and 2D).

Objective function values for the increase in sucrose accumulation in sugarcane culm tissue (maximization of the proportion of sucrose accumulation in the vacuole relative to sucrose hydrolysis by invertase) using up to five enzyme concentrations as adjustable parameters [23]: no constraints applied (A), total enzyme activity constraint applied (B), homeostatic constraint applied (C), total enzyme activity constraint and homeostatic constraint applied (D).

Objective function values for the increase in sucrose accumulation in sugarcane culm tissue (maximization of the proportion of sucrose accumulation in the vacuole relative to sucrose hydrolysis by invertase) using up to five enzyme concentrations as adjustable parameters [23]: no constraints applied (A), total enzyme activity constraint applied (B), homeostatic constraint applied (C), total enzyme activity constraint and homeostatic constraint applied (D).

Metabolite concentration related constraints directly can be applied only in kinetic models while total enzyme activity constraint can be implemented also in stoichiometric models. Despite the fact that the above described organism level constraints are proposed and applied several decades ago, they are also frequently ignored. The calculation of constraints for stoichiometric models using kinetic models has not been applied practically, as far as the authors are aware.


Towards a kinetic model of the Entner-Doudoroff pathway in Zymomonas mobilis

ENGLISH ABSTRACT: Metabolic networks of cellular systems are complex, in that there are numerous components with multiple non-linear interactions. To understand how these networks work they are often split into manageable pieces and studied individually. However, an individual part is unable to account for the complex properties of systems. In order to study these interactions the eld of systems biology has developed. Systems biology makes use of computers to construct models as a method to describe aspects of living systems. Using cellular pathways, kinetic models of metabolic pathways can be constructed and used as a tool to study the biological systems and provide a quantitative description. This thesis describes the quantitative analysis of a bacterium using a systems biology approach. Zymomonas mobilis is a rod shaped, Gram-negative, non-mobile facultative anaerobe and has one of the fastest observed fermentations, yet least energy e cient extractions found in nature. Furthermore it is the only known micro-organism to use the Entner-Doudoro (2-keto-3-deoxy-6- phosphogluconate) pathway anaerobically. The low energy yield of fermentation in Z. mobilis is a result of the usage of the Entner-Doudoro glycolytic pathway, which has half the energy yield per mol substrate compared to the well known Embden-Meyerhof-Parnas glycolytic pathway. The work presented in this thesis forms part of a larger project to compare glycolytic regulation in di erent micro-organisms Z. mobilis, Escherichia coli, Saccharomyces cerevisiae and Lactococcus lactis. These organisms were chosen based on their usage of di erent glycolytic mechanisms. Kinetic models are suitable tools to draw a comparison between these organisms. The emphasis here is on the construction of a kinetic model of the Entner-Doudoroff glycolytic pathway as it occurs in Z. mobilis. The aim of this thesis was to characterise as many of the Entner-Doudoro pathway enzymes as possible, under standard conditions. This was done using enzyme assays, to obtain the kinetic parameters of each of the enzymes. Microtitre plate assays were used to characterise most of the enzymes of the Entner-Doudoro pathway. However, not all characterisations could be done using plate assay methods, as some intermediates were not commercially available to perform coupled assays. Nuclear magnetic resonance (NMR) spectroscopy was used to characterise these enzymes. These experimentally obtained parameters were then incorporated in a mathematical framework. Time simulations on the initial model were unable to reach a steady-state, with a build up of metabolic intermediates. A secondary model was constructed (using calculated maximal activities) which allowed us to identify discrepancies in the initial model. This showed that the experimentally determined maximal activities of three enzymes in lower glycolysis were unrealistically low, which might be due to protein denaturation by sonication. A nal model was constructed which incorporated a correction factor for these three enzymes. The models' predicted output (steady-state concentrations and ux) was compared to that of either literature or experimentally determined values, as a method to validate the model. The model output compared well to literature values. The constructed and partially validated kinetic model was then used as an analytical tool to identify points of control and regulation of glycolysis in Z. mobilis. The model presented in this work was also compared to published models. Our model relies much less on literature obtained values, and uses kinetic parameters experimentally determined under the same conditions. The parameters of the published models were obtained from the literature and in many instances, the assay conditions for these parameters were set-up to yield the maximum activity under non-physiological conditions. Furthermore, the number of excluded or assumed parameters is much less in our model. However, introduction of a milder, more predictable extraction technique for preparing cell lysates, should be considered for future work, to obtain the parameters that was not determined during this study. The published models do include reactions not included in our model (e.g ATP metabolism), which should be considered for inclusion, as we strive to construct a detailed kinetic model of glycolysis in Z. mobilis in the future.

AFRIKAANSE OPSOMMING: Sellul^ere metaboliese netwerke is komplekse stelsels, omdat hulle bestaan uit talle komponente met verskeie nie-lineêre interaksies. Om die funksionering van hierdie netwerke te verstaan, word hulle dikwels in hanteerbare stukke verdeel en individueel bestudeer. 'n Enkele komponent is egter nie in staat om die komplekse eienskappe van sulke stelsels te verklaar nie. Die veld van sisteembiologie het ontwikkel met die doel om sulke stelsels te bestudeer. Sisteembiologie maak gebruik van rekenaarmodelle as 'n metode om aspekte van lewende sisteme te beskryf. Kinetiese modelle van metaboliese paaie word gebou en gebruik as gereedskap om die biologiese stelsels te bestudeer en 'n kwantitatiewe beskrywing te bekom. Hierdie tesis beskryf die kwantitatiewe ontleding van 'n bakterie deur middel van 'n sisteembiologiese benadering. Zymomonas Mobilis is 'n staafvormige, Gram-negatiewe, nie-mobiele fakultatiewe ana erobe, en het een van die vinnigste waargenome fermentasies, maar met die minste energie-doeltre ende ekstraksie wat in die natuur aangetref word. Verder is dit die enigste bekende mikro-organisme wat die Entner-Doudoro (2-keto-3-dioksi-6-fosfoglukonaat) pad ana erobies gebruik. Die lae-energieopbrengs van fermentasie in Z. mobilis is 'n gevolg van die gebruik van die Entner-Doudoro metaboliese pad, wat die helfte van die energie-opbrengs per mol substraat lewer, in vergelyking met die bekende Embden-Meyerhof-Parnas pad. Die werk wat in hierdie tesis aangebied word, vorm deel van 'n groter projek om glikolitiese regulering in verskillende mikro-organismes te vergelyk, naamlik Z. mobilis, Escherichia coli, Sac- charomyces en Lactococcus lactis. Hierdie organismes is gekies op grond van hul gebruik van verskillende glikolitiese meganismes. Kinetiese modellering is 'n handige metode om 'n vergelyking tussen hierdie organismes te trek. Hierdie werk fokus op die bou van 'n kinetiese model van die Entner-Doudoro glikolitiese metaboliese pad soos dit in Z. mobilis voorkom. Die doel van hierdie tesis was om so veel moontlik van die Entner-Doudoro ensieme onder standaard-toestande te karakteriseer. Die kinetiese parameters van elk van die ensieme is met behulp van ensimatiese essai's bepaal. Vir die meeste essai's is 96-put mikrotiterplate gebruik, maar nie al die karakteriserings kon met behulp van hierdie metode gedoen word nie, omdat sommige intermediate nie kommersieel beskikbaar was om gekoppelde essai's mee uit te voer nie. Kernmagnetiese resonansie (KMR) spektroskopie is gebruik om hierdie ensieme te karakteriseer. Die eksperimenteel bepaalde parameters is opgeneem in 'n wiskundige raamwerk. Tydsimulasies op die aanvanklike model was nie in staat om 'n bestendige toestand te bereik nie, omdat metaboliete opgebou het. 'n Sekond^ere model is gebou (met behulp van berekende maksimale aktiwiteite) wat ons toegelaat om teenstrydighede in die aanvanklike model te identi seer. Dit het getoon dat die eksperimenteel bepaalde maksimale aktiwiteite van drie ensieme in die laer gedeelte van glikolise te laag was, waarskynlik as gevolg van prote en denaturering tydens die ultrasoniese disintegrasieproses. 'n Finale model is gebou waarin 'n korreksiefaktor vir hierdie drie ensieme opgeneem is. Die modelle se voorspelde uitset (bestendige toestand konsentrasies en uksie) is vergelyk met waardes uit die literatuur of wat ons self bepaal het, as 'n metode om die model te valideer. Die model uitset was in goeie ooreenstemming met hierdie waardes. Die gedeeltelik gevalideerde kinetiese model is voorts gebruik as 'n analitiese instrument om beheer en regulering van glikolise in Z. mobilis te ondersoek. Die model wat in hierdie werk ontwikkel is, is ook vergelyk met die vorige gepubliseerde modelle. Ons model berus baie minder op waardes uit die wetenskaplike literatuur, en maak gebruik van parameters wat eksperimenteel bepaal is, onder identiese toestande. Die parameters van die gepubliseerde modelle is meesal verkry uit die literatuur, en in baie gevalle was die eksperimentele kondisies vir hierdie analises opgestel om die maksimale aktiwiteit te lewer onder nie- siologiese toestande. Verder bevat ons model minder parameters wat of uitgesluit is of wie se waardes aangeneem moes word. In toekomstige werk sal daar egter klem gel^e moet word op 'n minder wisselvallige ekstraksietegniek vir die verkryging van selekstrakte, om sodoende parameters te identi seer wat nie in hierdie werk bepaal kon word nie. Die gepubliseerde modelle sluit ook reaksies in wat nie ingesluit is in ons model nie (bv. ATP metabolisme). Hierdie sou in ag geneem moet word vir insluiting in 'n toekomstige uitgebreide model, om daarna te streef om 'n gedetailleerde kinetiese model van glikolise in Z. mobilis te bou.


Resumo

Kinetic modeling of metabolic pathways has become a major field of systems biology. It combines structural information about metabolic pathways with quantitative enzymatic rate laws. Some of the kinetic constants needed for a model could be collected from ever-growing literature and public web resources, but they are often incomplete, incompatible, or simply not available. We address this lack of information by parameter balancing, a method to complete given sets of kinetic constants. Based on Bayesian parameter estimation, it exploits the thermodynamic dependencies among different biochemical quantities to guess realistic model parameters from available kinetic data. Our algorithm accounts for varying measurement conditions in the input data (pH value and temperature). It can process kinetic constants and state-dependent quantities such as metabolite concentrations or chemical potentials, and uses prior distributions and data augmentation to keep the estimated quantities within plausible ranges. An online service and free software for parameter balancing with models provided in SBML format (Systems Biology Markup Language) is accessible at www.semanticsbml.org. We demonstrate its practical use with a small model of the phosphofructokinase reaction and discuss its possible applications and limitations. In the future, parameter balancing could become an important routine step in the kinetic modeling of large metabolic networks.

SPECIAL ISSUE

This article is part of the B: Robert (Bob) A. Alberty Festschrift special issue.


A comparative analysis of kinetic models of erythrocyte glycolysis.

  • APA
  • Autor
  • BIBTEX
  • Harvard
  • Padrão
  • RIS
  • Vancouver

Research output : Contribution to Journal › Article › Academic › peer-review

T1 - A comparative analysis of kinetic models of erythrocyte glycolysis.

N2 - Since the 1970s, with Heinrich as a pioneer in the field, numerous kinetic models of erythrocyte glycolysis have been constructed. A functional comparison of eight of these models indicates that the production of ATP and GSH in the red blood cell is largely controlled by the demand reactions. The rate characteristics for the supply and demand blocks indicate a good homeostatic control of ATP and GSH concentrations at different work loads for the pathway, while the production rates of ATP and GSH can be adjusted as needed by the demand reactions. © 2007 Elsevier Ltd. All rights reserved.

AB - Since the 1970s, with Heinrich as a pioneer in the field, numerous kinetic models of erythrocyte glycolysis have been constructed. A functional comparison of eight of these models indicates that the production of ATP and GSH in the red blood cell is largely controlled by the demand reactions. The rate characteristics for the supply and demand blocks indicate a good homeostatic control of ATP and GSH concentrations at different work loads for the pathway, while the production rates of ATP and GSH can be adjusted as needed by the demand reactions. © 2007 Elsevier Ltd. All rights reserved.


Assista o vídeo: Matematyczny pokaz naukowy cz 1 (Dezembro 2021).