Em formação

Derivação da cinética de Michaelis-Menten para simulações estocásticas discretas

Derivação da cinética de Michaelis-Menten para simulações estocásticas discretas


We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

De acordo com este artigo,

[…] A função de propensão para a reação de conversão S → P no caso estocástico discreto bem misturado pode ser escrita $ a (S) = frac {V_ {max} cdot S} {K_m + S / Omega} $ onde $ Omega $ é o volume do sistema.

Eu não entendo muito bem como esta fórmula é derivada do não discreto Cinética de Michaelis-Menten $ v = frac {V_ {max} cdot [S]} {K_M + [S]} $ (ver Wikipedia). De acordo com meu entendimento, $ [S] = S / Omega $. Se aplicarmos isso à fórmula da Wikipedia, obteremos $$ frac {V_ {max} cdot [S]} {K_M + [S]} = frac {V_ {max} cdot S / Omega} {K_M + S / Omega} = frac {V_ {max} cdot S} { Omega cdot K_M + S} $$ que não é o mesmo que $ frac {V_ {max} cdot S} {K_m + S / Omega} $. Então, como se pode derivar a fórmula do artigo citado (se estiver correta)? Se não, como podemos obter corretamente uma função de propensão discreta a partir da cinética de Michaelis-Menten?


Acho que descobri sozinho. Como $ a (S) $ está em $ frac {mol} {s} $ e $ v $ está em $ frac { frac {mol} {l}} {s} $, temos que $ a (S ) = v * Omega $. Portanto, $$ a (S) = frac {V_ {max} cdot [S]} {K_M + [S]} cdot Omega = frac {V_ {max} cdot S / Omega cdot Omega } {K_M + S / Omega} = frac {V_ {max} cdot S} {K_M + S / Omega}. $$


Cinética enzimática estocástica dependente do tempo: solução de forma fechada e bimodalidade transitória

Derivamos uma solução aproximada de forma fechada para a equação química mestre que descreve o mecanismo da reação de Michaelis-Menten de ação enzimática. Em particular, assumindo que a probabilidade de um complexo se dissociar em enzima e substrato é significativamente maior do que a probabilidade de um evento de formação de produto, obtemos expressões para as distribuições de probabilidade marginal dependentes do tempo do número de substrato e moléculas de enzima. Para as condições iniciais da função delta, mostramos que a distribuição do substrato é unimodal em todos os momentos ou torna-se bimodal em momentos intermediários. Essa bimodalidade transitória, que não tem contrapartida determinística, se manifesta quando o número inicial de moléculas de substrato é muito maior do que o número total de moléculas de enzima e se a frequência de eventos de ligação enzima-substrato é grande o suficiente. Além disso, mostramos que nossa solução de forma fechada é diferente da solução da equação mestre química reduzida por meio da aproximação estocástica discreta de Michaelis-Menten amplamente usada, onde a propensão para o decaimento do substrato tem uma dependência hiperbólica do número de moléculas do substrato . As diferenças surgem porque o último não leva em consideração as flutuações do número de enzimas, embora nossa abordagem as inclua. Confirmamos por meio da simulação estocástica de todas as etapas elementares da reação no mecanismo de Michaelis-Menten que nossa solução de forma fechada é precisa sobre uma região maior do espaço de parâmetros do que a obtida usando a aproximação estocástica discreta de Michaelis-Menten.


Resumo

É apresentada uma análise das suposições de aproximação que fundamentam uma derivação recentemente proposta da equação de taxa de reação determinística tradicional a partir de uma formulação estocástica discreta de cinética química. É mostrado que se o sistema está próximo o suficiente do limite termodinâmico, no qual as populações moleculares e o volume de conteúdo se aproximam do infinito de tal forma que as concentrações moleculares permanecem finitas, então as suposições aproximadas necessárias serão justificadas para praticamente todas espacialmente sistemas homogêneos que provavelmente serão encontrados.


Resumo

Os sistemas naturais são, quase por definição, heterogêneos: isso pode ser uma dádiva ou um obstáculo a ser superado, dependendo da situação. Tradicionalmente, ao construir modelos matemáticos desses sistemas, a heterogeneidade tem sido normalmente ignorada, apesar de seu papel crítico. No entanto, nos últimos anos, os métodos computacionais estocásticos tornaram-se comuns na ciência. Eles são capazes de explicar apropriadamente a heterogeneidade, de fato, eles são baseados na premissa de que os sistemas contêm inerentemente pelo menos uma fonte de heterogeneidade (a saber, heterogeneidade intrínseca).

Nesta mini-revisão, damos uma breve introdução à modelagem teórica e simulação em biologia de sistemas e discutimos as três diferentes fontes de heterogeneidade em sistemas naturais. Nosso tópico principal é uma visão geral dos métodos de simulação estocástica em biologia de sistemas.

Existem muitos tipos diferentes de métodos estocásticos. Nosso foco é um grupo que se tornou especialmente popular em biologia de sistemas, bioquímica, química e física. Esses métodos estocásticos de estado discreto não seguem os indivíduos ao longo do tempo, ao invés disso, eles rastreiam apenas as populações totais. Eles também assumem que o volume de interesse é espacialmente homogêneo. Fornecemos uma visão geral desses métodos, com uma discussão sobre as vantagens e desvantagens de cada um, e sugerimos quando cada um é mais apropriado para uso. Também incluímos referências a implementações de software deles, para que iniciantes possam começar a usar métodos estocásticos para problemas práticos de interesse.


Fundo

Os sistemas bioquímicos com um pequeno número de componentes em interação têm sido cada vez mais estudados nos últimos anos, visto que são alguns dos sistemas mais básicos da biologia celular [1–3]. Os efeitos estocásticos podem influenciar fortemente a dinâmica de tais sistemas. A aplicação de modelos de equações diferenciais ordinárias determinísticas (EDO) a eles, que aproximam o número de partículas como concentrações contínuas, pode levar a resultados confusos [4, 5]. Em alguns casos, mesmo os sistemas com grandes populações não podem ser modelados com precisão por EDOs. Por exemplo, quando próximo a um regime de bifurcação, as aproximações ODE não podem reproduzir o comportamento do sistema para alguns valores de parâmetros [6]. Os sistemas estocásticos podem ser modelados usando processos de Markov discretos. A densidade de estados de um sistema de reação química estocástica bem agitada em cada ponto no tempo é dada pela equação mestre química (CME) [7, 8]. O algoritmo de simulação estocástica (SSA) [9] é um método exato para simular trajetórias do CME conforme o sistema evolui no tempo.

O SSA pode ser computacionalmente intensivo para ser executado em problemas realistas e métodos alternativos, como o τ-leap foram desenvolvidos para melhorar o desempenho [10]. Em vez de simular cada reação, o τ-leap executa várias reações ao mesmo tempo, 'saltando' assim ao longo do eixo histórico do sistema. Isso significa que, ao contrário do SSA, o τ- salto não é exato A precisão é mantida, não permitindo que muitas reações ocorram por etapa. O tamanho de cada passo de tempo, τ, determina o número de reações que ocorrem durante essa etapa, dado por um número aleatório de Poisson.

Esse ganho de velocidade deve ser equilibrado com perda de precisão: etapas maiores significam menos cálculos, mas precisão reduzida. Muitos comuns τ-leap implementações empregam um tamanho de passo variável, usando o tamanho de passo ideal τ em cada ponto é crucial para a precisão do método [10-12]. No entanto, uma implementação de etapa fixa pode ser útil em alguns casos. Embora possa ser menos eficiente, é muito mais fácil de implementar do que equivalentes de etapa variável. Mais importante ainda, a estrutura de extrapolação que descrevemos neste artigo requer um método de etapas fixas.

O original τ-leap como descrito por Gillespie [10] é conhecido como o Euler τ-leap, pois pode ser comparado ao método de Euler para resolver EDOs. Foi demonstrado que existe uma ordem de convergência fraca sob a escala τ → 0 (escala tradicional) [13, 14] e V → ∞ (escala de grande volume) [15], onde V é o volume do sistema. No mesmo jornal, Gillespie também propôs o ponto médio τ-leap method [10], que possui convergência de ordem superior em alguns casos [15, 16]. Tian e Burrage [17] propuseram uma variante conhecida como Binomial τ- método de salto que evita problemas com espécies químicas se tornando negativas. Só recentemente mais trabalho foi feito na construção de métodos estocásticos de ordem superior. Um desses métodos é o corrigido aleatoriamente τ-leap [18], onde a cada passo de tempo uma correção aleatória é adicionada ao número aleatório de Poisson que determina o número de reações naquele passo. Dada uma correção aleatória adequada, os erros de ordem mais baixos nos momentos podem ser cancelados. Desta forma, métodos com convergência de até ordem fraca dois para média e covariância foram construídos. Mais recentemente, Anderson e Koyama [19] e Hu et al.[20] propôs outro método fraco de segunda ordem, o θ-trapezoidal τ-leap, que é uma adaptação do solver de equação diferencial estocástica (SDE) de Anderson e Mattingly [21] para o caso estocástico discreto.

Neste artigo, apresentamos uma estrutura para melhorar a ordem de precisão dos métodos estocásticos de passo fixo existentes, usando-os em conjunto com a extrapolação de Richardson, uma técnica bem conhecida para melhorar a ordem de precisão de um solucionador numérico combinando conjuntos de simulações com tamanhos de passos diferentes [22]. A extrapolação foi originalmente desenvolvida para solucionadores de ODE, mas também foi aplicada com sucesso aos métodos SDE [23]. Nossa abordagem tem três vantagens principais:

Aumenta a ordem de precisão dos métodos fornecidos a ele. Isso é desejável pela razão óbvia de que as soluções resultantes são mais precisas, bem como por que passos de tempo maiores podem ser usados ​​para atingir um determinado nível de precisão, reduzindo o custo computacional. Isso é discutido mais detalhadamente em nossas Conclusões.

Ele pode ser aplicado a qualquer solucionador de etapas fixas, por exemplo, métodos inerentemente de ordem superior, como o θ-trapezoidal τ-leap ou métodos com uma região de estabilidade estendida, como os métodos estocásticos de Runge-Kutta [24].

As soluções de ordem superior resultantes podem ser extrapoladas novamente para fornecer soluções com uma ordem ainda maior, pois não há limite (teórico) para o número de vezes que um método pode ser extrapolado (embora erros estatísticos possam obscurecer os resultados se o método for muito preciso - veja a seção Erro de Monte Carlo).

Nossos métodos extrapolados podem ser úteis para pesquisadores em biologia e bioquímica, pois são fáceis de implementar e podem simular com precisão e rapidez sistemas estocásticos discretos que, de outra forma, poderiam ser muito intensivos em computação.

Mostramos como a estrutura de extrapolação pode ser aplicada a algoritmos estocásticos de passo fixo usando os exemplos de Euler de passo fixo τ-leap, midpoint τ-leap [10] e θ-trapezoidal τmétodos -leap [20]. O procedimento de extrapolação depende fortemente da existência de uma expansão de erro global apropriada para o erro fraco do método numérico. Uma vez que isso é conhecido, a extrapolação consiste em aritmética simples. Calculamos essa expansão para um método de primeira ordem fraco arbitrário, o que nos permite usar a extrapolação para obter soluções de ordem superior. A ordem fraca de todos os momentos de tais métodos pode ser melhorada por extrapolação. Para reforçar isso, realizamos uma análise de erro simples, comparando as equações para a média verdadeira e numérica de Euler τMétodo -leap, vemos que seu erro global é de ordem um, e extrapolando-o aumenta a ordem para dois para o caso de reações de ordem zero e de primeira ordem. Usando simulações numéricas, demonstramos que isso é verdade para dois sistemas de teste de primeira ordem e três de ordem superior com o Euler, ponto médio e θ-trapezoidal τmétodos -leap. Além disso, os métodos extrapolados apresentam erros consistentemente menores e, em muitos casos, convergência de ordem visivelmente superior nos primeiros dois momentos (a falta de convergência em algumas das simulações é discutida na Seção Erro de Monte Carlo). Finalmente, demonstramos que a estrutura de extrapolação pode ser usada para fornecer soluções numéricas de ordem ainda mais elevada, aplicando uma segunda extrapolação ao método de Euler τmétodo -leap.

O restante deste trabalho está organizado da seguinte forma. Começamos com uma visão geral do SSA e do τMétodos -leap que usaremos mais tarde. Em seguida, discutimos a extrapolação de Richardson para EDOs e SDEs e apresentamos a estrutura estocástica discreta extrapolada. Fornecemos resultados numéricos para apoiar nossas afirmações de que a extrapolação reduz o erro dos métodos de etapas fixas. Finalmente, discutimos o erro de Monte Carlo e damos nossas conclusões. As derivações das expansões de erro globais para SDEs e métodos estocásticos discretos e material relacionado são apresentadas no Apêndice.

Visão geral dos métodos de simulação estocástica

O SSA de Gillespie [9] é um método estatisticamente exato para simular caminhos de um processo de salto de Markov. As duas suposições básicas do SSA são (i) que as moléculas individuais não são explicitamente rastreadas e (ii) há colisões não reativas frequentes. Assim, assumimos que o sistema é bem misturado e homogêneo.

O SSA simula um sistema de reações bioquímicas com N espécies e M reações, interagindo dentro de um volume fixo V a temperatura constante. As populações de espécies químicas (como números moleculares, não concentrações) no tempo t são representados como um vetor de estado x ≡ X (t) ≡ (X 1,…, X N) T. As reações são representadas por uma matriz estequiométrica ν j ≡ (ν 1 j, ..., ν Nj) T, onde j = 1,…,M, composto de M vetores estequiométricos individuais. Cada vetor estequiométrico representa uma reação j ocorrendo e o sistema mudando de estado x para x + ν j. Cada reação ocorre em um intervalo t t + τ) com probabilidade relativa uma j(x)dt, Onde uma j é a função de propensão do j-ª reação. As funções de propensão são dadas pela cinética de ação em massa das espécies químicas reagentes. Para mais detalhes, o leitor deve consultar a Ref. [9]. As variáveis X,ν j e uma j(X) caracterizam completamente o sistema em cada momento.

Algoritmo 1. Método direto SSA Com o sistema no estado X n no tempo t n:

Gere dois números aleatórios r1 e r2 a partir da distribuição uniforme de intervalo de unidade U (0, 1).

Encontre o tempo até a próxima reação τ = 1 a 0 ln 1 r 1, onde a 0 (X n) = ∑ j = 1 M a j (X n).

Encontre a próxima reação j de ∑ m = 1 j - 1 a m (X n) & lt a 0 r 2 ≤ ∑ m = j M a m (X n).

Atualizar t n + 1 = t n + τ e X n+ 1 = X n + ν j.

O método direto requer dois números aleatórios recém-gerados em cada etapa de tempo. Embora existam outras implementações de SSA, como o Next Reaction Method [25] e o Optimized Direct Method [26], que podem ser mais econômicos, em geral o SSA é computacionalmente caro.

Τmétodo de salto

o τO algoritmo -leap salta ao longo do eixo histórico do SSA avaliando grupos de reações de uma só vez [10]. Isso significa significativamente menos cálculos, ou seja, menor tempo computacional, por simulação, mas a precisão da simulação está comprometida: não sabemos exatamente quantas reações ocorreram durante cada etapa de tempo, nem podemos dizer com mais precisão quando cada reação ocorre do que em qual intervalo de tempo. o condição de salto define um limite superior para o tamanho de cada passo de tempo τ: deve ser tão pequeno que as propensões não mudem significativamente durante sua duração, ou seja, a mudança de estado com o tempo t para t + τ é muito pequeno [10]. Desde a τ é pequena, a probabilidade uma(x)τ que uma reação ocorre durante t t + τ) também é pequeno, então o número de vezes K jcada reação disparada em um intervalo de tempo pode ser aproximada por P (a j (x) τ), uma variável aleatória de Poisson com média e variância uma j(x)τ[10]. O Euler τ-leap o algoritmo é o básico τ-leap método, e corresponde ao método de Euler para resolver EDOs ou o método de Euler-Maruyama para resolver SDEs.

Algoritmo 2. Euler τmétodo de salto Com o sistema no estado X nno tempo t n, e um intervalo de tempo τ:

Gerar M Números aleatórios de Poisson k j = P (a j (X n) τ).

Atualizar t n + 1 = t n + τ e X n + 1 = X n + ∑ j = 1 M ν j k j.

O Euler τ-leap tem ordem fraca um [13-15]. Embora um trabalho considerável tenha sido feito para melhorar o mecanismo de seleção dos passos de tempo τ[10-12] e eliminando etapas que resultariam em populações negativas [17, 27-29], isso não afeta a ordem do método, limitando sua precisão. Métodos com ordem superior são a única maneira de melhorar a precisão além de um certo ponto. Percebendo isso, Gillespie também propôs uma ordem superior τ-pelo método, o ponto médio τ-leap [10]. Isso é semelhante ao método do ponto médio para EDOs, onde a cada etapa uma estimativa é feita do gradiente de X no t n + τ/2. X n é então incrementado usando essas informações extras para fornecer uma aproximação mais precisa.

Algoritmo 3. Ponto médio τmétodo de salto Com o sistema no estado X nno tempo t n, e um intervalo de tempo τ:

Calcule X ′ = X n + 1 2 τ ∑ j = 1 M ν j a j (X n).

Gere M números aleatórios de Poisson k j = P (a j (X ′) τ).

Atualizar t n + 1 = t n + τ e X n + 1 = X n + ∑ j = 1 M ν j k j.

Embora sob a escala τ → 0 o ponto médio τ-leap tem a mesma ordem de precisão na média que o Euler τmétodo de salto, sob a escala de grande volume tem ordem dois fraca [15, 16]. Nossas simulações numéricas também sugerem que ele fornece aproximações de ordem superior para os dois primeiros momentos para sistemas lineares e não lineares (embora isso não esteja claro na literatura). No entanto, o erro de truncamento local de sua covariância é de primeira ordem [16].

Θ-trapezoidal τmétodo de salto

Com base no método SDE de Anderson e Mattingly [21], o θ-trapezoidal τ-leap [20] é um método de segunda ordem fraco. Consiste em duas etapas, uma etapa preditora com tamanho θτ e uma etapa corretora com tamanho (1 - θ)τ que visa cancelar quaisquer erros cometidos na primeira etapa.

Algoritmo 4. θ-trapezoidal τmétodo de salto Para um θ especificado ∈ (0,1), α 1 = 1 2 (1 - θ) θ, α 2 = (1 - θ) 2 + θ 2 2 (1 - θ) θ. Com o sistema no estado X nno tempo t n, e um intervalo de tempo τ:

Gerar M Números aleatórios de Poisson k j ′ = P (a j (X n) θτ), j = 1,…, M.

Calcule a etapa do preditor X ′ = X n + ∑ j = 1 M ν j k j ′.

Calcule l j = max α 1 a j (X ′) - α 2 a j (X n), 0.

Gerar M Números aleatórios de Poisson k j = P (l j (1 - θ) τ), j = 1,…, M.

Atualizar t n + 1 = t n + τ e X n + 1 = X ′ + ∑ j = 1 M ν j k j.

Especificamente, o θ-trapezoidal τO método -leap mostrou ter ordem de convergência fraca dois nos momentos, e um erro de truncamento local de O (τ 3 V - 1) para a covariância. τ = Vβ , 0 & lt β & lt 1 e na análise V → ∞, mas nas simulações o volume do sistema é mantido constante, portanto, parece que na prática isso também resulta em uma convergência de segunda ordem fraca na covariância [20].


Métodos

Extrapolação

A extrapolação de Richardson é uma técnica para aumentar a ordem de precisão de um método numérico, eliminando o (s) termo (s) de erro inicial em sua expansão de erro [23, 24]. Envolve resolver numericamente alguma função determinística Y(t) em um dado instante T=n τ usando o mesmo solucionador com tamanhos de etapas diferentes, onde definimos Y T τ como uma aproximação de Y(T) no tempo T usando o tamanho do passo τ. Y(T) pode ser escrito como

Onde ε g(τ) é o erro da solução aproximada em comparação com a verdadeira. Para um solucionador numérico geral, ε g(τ) pode ser escrito em termos de poderes do tamanho do passo τ:

onde o e k são vetores constantes e dependem apenas do tempo final e k 1& ltk 2& ltk 3,…. Eq. (2) nos diz que este método tem uma ordem de precisão k 1.

Essencialmente, a extrapolação de Richardson emprega a extrapolação polinomial das aproximações Y T τ q, q = 1, 2, ... e τ 1& gtτ 2& gt ..., para estimar Y T 0, ou seja, a solução numérica no limite do tamanho do passo zero, que corresponde a Y(T) (Figura 1). Cada extrapolação sucessiva remove o próximo termo de erro principal, que é a maior contribuição para o erro, aumentando assim a precisão da solução numérica e permitindo uma melhor estimativa Y(T).

Princípio de extrapolação de Richardson. Três soluções numéricas, com tamanhos de passo τ 1 = T, τ 2 = T 2, τ 3 = T 4 encontram estimativas cada vez mais perto da solução verdadeira Y (T) = Y T 0, ou seja, a solução numérica no limite de tamanho de passo zero. Eles podem ser extrapolados para encontrar uma estimativa muito próxima de Y T 0.

Para demonstrar isso, assuma um método numérico com tamanho de passo τ tem uma expansão de erro de

Por exemplo, o conhecido método de Euler para resolver EDOs tem essa expansão de erro. Agora em vez de τ, se usarmos um tamanho de passo τ/ 2, a expansão do erro é

Podemos tomar Y T τ, τ / 2 = 2 Y T τ / 2 - Y T τ, dando

O termo de erro inicial foi removido, resultando em uma aproximação de ordem superior. Isso pode ser repetido para obter uma ordem ainda maior de precisão usando mais aproximações iniciais Y T τ 1, ... Y T τ q, onde q pode ser qualquer número inteiro e τ 1& gtτ 2& gt… τ q. Definimos Y T τ 1, τ q como a solução extrapolada usando as aproximações iniciais Y T τ 1,… Y T τ q. A maneira mais fácil de visualizar isso é construir uma tabela de Neville (também chamada de tabela de Romberg) a partir das aproximações iniciais (Tabela 1).

A primeira coluna da tabela contém as aproximações numéricas iniciais. Em seguida, eles são extrapolados para encontrar a próxima coluna e assim por diante. Por exemplo, com três soluções iniciais YT τ, YT τ / 2, YT τ / 4, então YT τ, τ / 4 = 4 3 YT τ / 2, τ / 4 - 1 3 YT τ, τ / 2 (isto é facilmente calculado escrevendo primeiro uma fórmula semelhante à Eq. (3) para YT τ / 4, depois outra semelhante à Eq. (4) para YT τ / 2, τ / 4 e, mais uma vez, para YT τ, τ / 4 ) Em cada coluna subsequente, o próximo termo de erro inicial é cancelado, fornecendo uma solução de ordem ainda mais alta. Os coeficientes corretos para calcular cada novo termo da tabela de Neville podem ser encontrados em

Onde p=τ qr/τ qr+1 e k q é a ordem da solução na coluna q (c.f. Eq. (2)), e r=1,…,q-1. Isso pode ser generalizado para algum método de pedido com algum expansão de erro apropriada. A única condição para extrapolação é a existência de uma expansão de erro da forma na Eq. (2)

Método Bulirsch-Stoer

O método Bulirsch-Stoer é um solucionador de ODE preciso baseado na extrapolação de Richardson [25, 26]. Uma tabela de Neville é construída por extrapolação repetida de um conjunto de aproximações iniciais com tamanhos de etapas que são subintervalos diferentes de uma etapa geral maior τ, e é então usado para encontrar uma solução muito precisa. Isto acontece dentro de cada passo de tempo, permitindo τ para ser variado entre as etapas. Um método de ponto médio modificado (MMP, Algoritmo 3) é usado para gerar as aproximações iniciais na primeira coluna da tabela. Isso se presta bem a uma estrutura de extrapolação, já que o MMP subdivide cada etapa τ em n ̂ subetapas τ ̂ = τ / n ̂. Além disso, crucialmente, a expansão do erro do MMP contém apenas potências pares de τ ̂, resultando em convergência rápida [27].

Damos uma breve visão geral do método Bulirsch-Stoer determinístico aqui Ref. [28] tem uma excelente descrição do algoritmo, bem como um guia para sua implementação. Em cada etapa, uma coluna da mesa de Neville, k, no qual esperamos que as soluções aproximadas tenham convergido, bem como um tamanho de etapa τ são selecionados. A mesa de Neville é então construída executando k MMPs, com tamanhos de passo τ ̂ 1 = τ / 2,…, τ ̂ q = τ / n q, onde n q=2q,q=1,2,…,k e extrapolar sucessivamente as aproximações numéricas apropriadas. A convergência das soluções é avaliada com base na consistência interna da tabela de Neville, ou seja, a diferença entre a solução mais precisa em coluna k e isso na coluna k−1: da Tabela 1, isso é Δ Y (k, k - 1) = Y τ τ ̂ 1, τ ̂ k - Y τ τ ̂ 2, τ ̂ k. Conforme sucessivas aproximações iniciais Y τ τ ̂ q são adicionadas à primeira coluna, os resultados extrapolados em cada nova coluna convergem para a solução verdadeira e Δ Y(k,k-1) diminui. A aproximação final na coluna k é aceitável se e r r k≤1, onde e r r k é uma versão em escala de Δ Y(k,k-1) (consulte o Apêndice: algoritmo completo de Bulirsch-Stoer estocástico para obter mais detalhes). Se e r r k& gt1, a etapa é rejeitada e refeita com τ = τ 2.

Em uma implementação prática, a etapa inicial testa q=1,…,k m uma x, Onde k m uma x geralmente é definido como oito, a fim de estabelecer o k necessário para atingir a precisão exigida e garantir que o tamanho do passo seja razoável para os passos subsequentes, em seguida, teste a convergência apenas nas colunas k−1,k e k+1 [28]. Por causa de sua precisão, as etapas executadas pelo método Bulirsch-Stoer podem ser relativamente grandes em comparação com outros solucionadores numéricos. τ é alterado adaptativamente em cada etapa e é escolhido para minimizar a quantidade de trabalho realizado (ou seja, avaliações de função n ̂ + 1 do MMP) por tamanho da unidade. Desta forma, o método Bulirsch-Stoer adapta sua ordem e tamanho dos passos para maximizar a precisão e a eficiência computacional.

Método estocástico de Bulirsch-Stoer

O método SBS é baseado em sua contraparte determinística descrita na seção anterior. Existem algumas questões-chave que devem ser abordadas a fim de adaptá-lo com sucesso em um método estocástico. Os dois mais importantes estão interligados: primeiro, que quantidade deve ser calculada em cada etapa e, segundo, como a estocasticidade pode ser introduzida na imagem? O método determinístico de Bulirsch-Stoer calcula X n+1 a partir de X n usando o MMP para encontrar os estágios intermediários ao longo de [t,t+τ) No entanto, a estocasticidade não pode ser simplesmente adicionada a este esquema dentro ou fora do MMP, pois isso interferiria na extrapolação necessária para a tabela de Neville. Para atualizar o vetor de estado como na Eq. (1), devemos encontrar o número de reações por etapa.

Olhando para a fórmula de atualização para a trajetória de um processo de Markov de salto [29],

é claro que a quantidade que devemos calcular é ∫ t n t n + τ a (X (t)) dt, a fim de, então, tomar uma amostra de Poisson para a atualização (o τO método -leap se aproxima disso como uma(X(t n))τ) Assim, ao invés de calcular X n+1 usando diretamente o MMP, precisamos de uma maneira precisa de encontrar a integral das funções de propensão em cada etapa. Procedendo de maneira um tanto semelhante ao Algoritmo 3, chegamos ao Algoritmo 4: os estágios intermediários são encontrados usando o MMP, e as propensões são calculadas em cada estágio. Essas propensões intermediárias são então alimentadas em um método trapezoidal composto para fornecer uma estimativa precisa da integral.

Um ponto importante é que os estágios intermediários são resolvidos usando as equações de taxa de reação (Etapas 1, 3, 5 do Algoritmo 4), que fornecem a expectativa da trajetória estocástica ao longo de cada etapa, desde que ambos sejam iniciados no estado X n no tempo t n. Assim, encontramos o esperado ∫ t n t n + τ a (X (t)) dt usando integração de Romberg, e use isso para amostrar uma distribuição de Poisson a fim de incrementar X n. Este método é extremamente preciso para encontrar a média e totalmente estocástico, ou seja, cada simulação fornece uma realização estocástica diferente e a densidade de probabilidade total pode ser encontrada em um histograma de muitas simulações.

Chegamos agora à implementação do SBS. Primeiro, calculamos Δ a τ ̂ q (t n, t n + τ), a integral esperada das propensões sobre [t n,t n+τ), usando o Algoritmo 4 com vários tamanhos de etapas τ ̂ 1, τ ̂ 2,…. Em seguida, extrapolamos esses usando a tabela de Neville (Romberg) para chegar às soluções extrapoladas Δ uma extr (t n,t n+τ) isso é conhecido como integração Romberg. Uma vez que estes são suficientemente precisos, nós amostramos o número de reações como

Esta é a nossa aproximação da função de densidade de probabilidade subjacente em cada etapa. Combinado com o mecanismo de extrapolação descrito anteriormente e uma forma de adaptar o tamanho do passo, temos o método SBS completo.

O tamanho do passo é escolhido calculando a quantidade

Onde τ é o passo de tempo atual, τ k é o próximo passo de tempo hipotético para a coluna da tabela Romberg k e S 1 e S 2 são parâmetros de segurança, introduzidos no próximo parágrafo. Aqui e r r k é o erro local relativo a uma tolerância mista, com ordem O (τ 2 (k - 1) + 1) (ver Apêndice: algoritmo completo de Bulirsch-Stoer estocástico). Seu valor ideal é exatamente um: se for menor do que este, a etapa pode ter sido aumentada, e se for maior, significa que nosso limite de erro foi excedido e a etapa deve ser refeita usando um menor τ. Em cada etapa, um próximo passo de tempo candidato é selecionado para cada coluna da tabela de Neville k. Como sabemos alguma medida do trabalho realizado para cada k (o número de avaliações de função do MMP), podemos calcular a eficiência de cada um dos k tempos do candidato como o trabalho por unidade τ. Em seguida, selecionamos o candidato τ k que dá a maior eficiência. Para resumir, deixamos a descrição completa e a implementação passo a passo do SBS (Algoritmo 5) até o Apêndice.

O SBS usa vários parâmetros diferentes (consulte o Algoritmo 5), todos os quais têm algum efeito nos resultados. S 1 e S 2 são ambos fatores de segurança que redimensionam o próximo passo de tempo em alguma quantidade: quanto menores eles são, menor o passo de tempo e mais precisa a solução. Como sempre, no entanto, há um compromisso entre o tamanho da etapa e a velocidade, portanto, deve-se ter o cuidado de otimizar os parâmetros para obter a eficiência máxima. O mesmo também é verdadeiro para os vetores uma t o eu, a tolerância de erro absoluta e r t o eu, a tolerância de erro relativa. Eles são usados ​​para dimensionar o erro que é calculado a partir da consistência interna da tabela de Romberg. Eles geralmente são definidos bem baixos: cerca de 10-6 é comum. Há uma consideração adicional com o SBS, a saber, a da coluna de convergência, k. Mesmo quando os fatores de segurança são definidos como altos (significando passos de tempo maiores), o SBS pode atingir uma precisão muito alta simplesmente fazendo outra extrapolação e indo para uma coluna mais alta. Por este motivo, a relação entre os fatores de segurança e a precisão não é direta, sendo aconselhável verificar os passos de tempo e a coluna de convergência para cada novo conjunto de parâmetros.

Extensão: SBS-DA

Existe um esquema alternativo para a Eq. (6) para encontrar a atualização estocástica para o vetor de estado: este é o "grau de avanço", ou abordagem DA, e chamamos o método resultante de SBS-DA. Seu foco é o M× 1 processo aleatório Z j(t),j=1,…,M, o número de vezes que cada reação ocorre durante [0,t) [30, 31]. Z j(t) está relacionado ao vetor de estado X(t) por

Na verdade, X(t) é determinado exclusivamente por Z(t), Onde Z(t) é um M× 1 vetor [31]. Isso nos permite usar a abordagem DA para calcular o número de reações por etapa e, em seguida, retornar à abordagem da população para atualizar o vetor de estado, usando

onde nós definimos Z j([t,t+τ)) como o número de reações que ocorrem durante [t,t+τ) Observe que a Eq. (8) tem a mesma forma que as Eqs. (1) e (6). Na verdade, K j=Z j([t,t+τ)): no caso do SSA, o intervalo de tempo tende a ser muito pequeno e ocorre apenas uma reação, mas para o τ-leap e SBS é muito maior, então mais reações podem ocorrer. Da mesma forma que a Eq. (5), sabemos que [32]

A fim de encontrar a atualização do vetor de estado, devemos resolver para a média (e variância, veja abaixo) de K je amostra de acordo com a Eq. (9). As equações para a evolução da média e variância de K j,µ j(s) e V j(s), respectivamente (onde s executa apenas uma etapa [t,t+τ)), pode ser derivada de sua equação mestre e assumir a forma [19] (ver também [31])

Onde s ∈ [t,t+τ), X n é o valor do vetor de estado no início da etapa e fjj ′ (X n) = ∑ i = 1 N ∂ aj (X n) ∂ xi ν ij ′, j, j ′ = 1, ..., M são os elementos de um M×M matriz (note que nós apenas lidamos com seus elementos diagonais no caso da variância). Eqs. (10) e (11) devem ser resolvidos simultaneamente com as condições iniciais µ j(t)=V j(t) = 0 para encontrar µ j(t+τ) e V j(t+τ) Deve-se notar que eles são exatos apenas para sistemas com propensões lineares. No caso de propensões não lineares, as equações de momento contêm momentos mais altos e obtemos as Eqs. (10) e (11) por um argumento de fechamento padrão: nós Taylor expandimos as propensões e truncamos na primeira ordem [19]. Na verdade, a Eq. (11) só é necessário porque a Eq. (10) não é exato no caso geral. Para passos de tempo maiores, isso pode levar a um erro considerável, então devemos aproximar a verdadeira distribuição de Poisson de K j com um gaussiano cuja variância foi corrigida. Isso leva ao esquema de atualização

que agora substitui a Eq. (6). Aqui, ⌊ ⌉ denota o arredondamento para o inteiro mais próximo e o valor dez foi escolhido heuristicamente, pois acima desse valor uma amostra de Poisson pode ser bem representada por uma amostra gaussiana com a média e variância apropriadas.

Esta abordagem é um pouco semelhante ao imparcial τ-leap [19]. The key difference is that Ref.[19] uses this scheme in the context of a fixed-stepsize τ-leap method, basing the entire method around this scheme. In contrast, the SBS-DA is grounded in the Bulirsch-Stoer method, which it uses for its stepsize selection and combination of MMP and Richardson extrapolation to find the Poisson increment. The DA approach is only one part of the whole SBS-DA, and is only used as an alternative to Eq. (6), in order to also find the variance of the number of reactions occurring over each step.

The SBS and SBS-DA methods both calculate the parameter of the Poisson sample but they take different approaches to this. The two key differences are that (1) the SBS-DA attempts to correct using the variance of the sampled distribution in order to better approximate the true Poisson parameter when the stepsize is large, but (2) it sacrifices some of its performance because of the inherent inaccuracies of Eqs. (10) and (11) (see Results, Higher order of accuracy and robustness section).


EXACT STOCHASTIC METHODS

The so-called exact stochastic methods correctly account for inherent stochastic fluctuations and correlations. In addition, the discrete nature of the system is considered. Hence, they remain valid for very small particle numbers.

However, since they explicitly simulate each reaction event in the system, they have time complexity approximately proportional to the overall number of particles present in the system. Therefore, they are prohibitively slow on large systems.

This function determines the probability P(τ,μ|x,t)dτ that starting in state x at time-point t the next reaction in the system will occur in the time interval [t + τ,t + τ + dτ] and will be of type Rµ. The reaction times are exponentially distributed (homogeneous Poisson process).

By iteratively drawing random numbers according to that density function and updating the system state according to the chosen reaction's stoichiometry the system can be simulated over time, one reaction event after the other. The Direct Method and the First Reaction Method, as well as the other exact methods described below, are mathematically equivalent but differ in how they calculate samples of P(τ,μ|x,t)dτ.

Método direto

The Direct Method [ 5, 52] uses two random numbers per step to separately compute (i) the stochastic time step τ and (ii) the type of the next reaction Rµ. A similar method was already proposed in [ 53]. The algorithm proceeds as follows:

The corresponding reaction Rµ is realized, i.e. the number of the participating molecules is increased or decreased according to the stoichiometry, and the time is incremented by τ. The whole process (1–4) is repeated as many times as necessary to reach the desired simulation time.

First Reaction Method

The First Reaction Method [ 5, 52] uses M random numbers in each step to compute putative reaction times τµ para cada um dos M reactions Rµ in the system. The τµ are exponentially distributed with parameter umaµ. The reaction with the smallest reaction time τmin = min(τ1,…,τM), is executed and all putative reaction times are recalculated before the next step. Because of the wasteful use of random numbers and redundant recalculations the First Reaction Method is computationally inefficient. However, it can be numerically advantageous because, due to the machine-dependent resolution of floating-point numbers, the Direct Method cannot calculate rare reaction events amongst very frequent reactions in the system.

Next Reaction Method

Gibson and Bruck [ 54] improved the computational complexity of the First Reaction Method by the intelligent use of data structures:

A so-called dependency graph stores dependencies between the reactions → redundant recalculations of umaµ are avoided.

Absolute putative reaction times τµ (since the beginning of simulation) are used instead of relative ones (since the last reaction event). Random numbers are ‘recycled’ during a reaction time update according to ⁠ .

An indexed priority queue contains all reactions sorted according to their putative reaction times → the next reaction can always be found in the root of the tree in constant time.

In the Next Reaction Method recalculations of umaµ are minimized, and, asymptotically, only one random number per step is needed. Therefore it performs well, in particular, in the case of many reacting species and reactions.

Variants and extensions of the exact methods

Cao et al. (Optimized Direct Method [ 55]) and McCollum et al. (Sorting Direct Method [ 56]) reduce the computational cost of the Direct Method by sorting the reactions according to their propensities.

Li and Petzold (Logarithmic Direct Method [ 57]) use partial sums of propensities and a binary search for the next reaction in order to accelerate the computation.

Slepoy et al. [ 58] describe an algorithm that, in each iteration, chooses the next reaction event in constant-time independent of the number of different reactions.

The Direct Method has also been extended to accommodate delayed reactions often used to model, e.g. transcription or translocation processes [ 59–62].

Finally, there exist fast implementations that utilize special hardware, such as field-programmable gate arrays (FPGA) [ 63, 64] or computer clusters [ 65].


Discussão

The standard approach for estimating enzyme kinetic parameters even today continues to be based on the 100-year old MM equation (Eq. 1) 5,6 . However, when enzyme concentration is high, this approach can lead to biased estimation (Fig. 2). Even when enzyme concentration is relatively low, it may not be possible to identify kinetic parameters (Figs 3 and 4). To overcome the limitations of the canonical approach, we proposed an estimation method based on an alternative to the MM equation: the tQ model (Eq. 2), which is derived with the total QSSA 26,27,28,29 . Because the estimation procedure with the tQ model is not biased regardless of enzyme or substrate concentrations (Fig. 2), more accurate and precise estimations can be made when pooled data from different experimental conditions are used, unlike the canonical approach (Figs 5 and 6). It appears thus that the tQ model is especially appropriate for creating a consistent Bayesian inferential framework, which becomes more accurate as more data is used.

The canonical enzyme kinetic assay based on the MM equation generally requires a large excess of substrate over enzyme 42 . However, such conditions impose experimental limitations and cannot be always guaranteed and verified 15 . For instance, it is hard to generate a high concentration of barely soluble substrate 24 , and a low concentration of substrate is required for sensitive kinetic analysis, e.g., in the case of QD-FRET-based probes 39,40,41 . Importantly, to analyze na Vivo enzyme kinetics, where enzyme concentration is often high 16,17,18 , our approach, but not the canonical approach, can be used. For example, one needs to estimate the kinetic parameters underlying drug metabolism by CYP enzymes in the liver in order to predict the effects of drugs, as is essential for drug development 43 . Because of dosing requirements for potent drugs, the amount of CYP enzyme can greatly exceed the drug amount in the liver 44,45 . Another large area where our estimation method can be applied is in the development of nanobiosensors, which measure na Vivo activity of a specific enzyme for precise diagnostics, because such enzymes are often in large excess over biosensors 46,47 .

KinTek Explorer has been widely used to estimate enzyme kinetic parameters from the progress curves 48,49,50 . This software provides the confidence contours, which reveal the relationships between the estimated parameters. This approach recommends using multiple data sets to narrow down the confidence contours and thus improve precision of estimates and resolve the unidentifiability issue. Our finding (Fig. 7a,b) can provide the specific type of data sets required for the identifiability of k gato e K M, so that the KinTek Explore could perform parameter estimation more efficiently for the Michales-Menten type of enzyme reactions.

Since the initial velocity estimation with the MM equation is not accurate when enzyme concentration is high (Fig. 1), the standard initial velocity based on the MM equation would also be inaccurate 15 . On the other hand, the tQ model accurately captures the initial velocity for all conditions, and thus the modified initial velocity assay based on the tQ model is likely to be accurate over a wider range of conditions. To simplify such estimation procedures, an interesting future study could derive an analogous Lineweaver-Burk plot or the Hanes-Woolf plot 8,9,15,42 for the tQ model.

Even with relatively large noise in the data (Fig. S1), our proposed method leads to accurate estimation (Figs 5, 6 and 7), indicating its robustness against experimental noise and some minor inaccuracy of the tQ model in certain ranges of parameter observed in 14,30 . Furthermore, if there are departures from simple non-inhibitory enzyme kinetics (e.g. inhibition of enzyme by product) 51,52 , our method can be easily adjusted by modifying the tQ model (see 53,54 for the tQ model for other enzyme kinetics). Our work can also be used to improve the estimation of the kinetics underlying diverse biological functions, such as gene regulation 55,56 , cellular rhythms 57,58,59 , quorum sensing 60,61 , signal cascade 62,63 and membrane transport 64,65 , where the MM equation has been widely used.


1. Introdução

Biochemical systems have traditionally been modeled by a set of ordinary differential equations (ODEs). The general form of the Reaction Rate Equations (RREs) in that approach can be formulated as

para eu = 1, … n, where the state variables xeu represent the concentrations of involved species, and functions feu are inferred from the various chemical reactions of the system. This set of ODEs is usually solved with numerical methods packages such as DASSL and DASPK (Brenan et al. 1996), ODEPACK (Hindmarsh, 1983), or CVODE (Cohen and Hindamarsh, 1996) for example. An important feature of the equations in this approach is that the system is deterministic and continuous. For systems where all chemical species are present in large copy numbers, it is reasonable to model every species by its concentration and the traditional ODE approaches seem to work very well. However, if the system is small enough that the molecular populations of some of the reactant species are small, from one to thousands, discreteness and stochasticity may play important roles in the dynamics of the system. Such a case occurs often in cellular systems (McAdams and Arkin, 1997, Arkin et al., 1998, Fedoroff and Fontana, 2002), which typically involve copy numbers of one or two for the number of genes of a given protein, on the order of tens to hundreds for the corresponding RNAs and on the order of thousands for regulatory proteins and enzymes. In that case, Equation (1) cannot accurately describe the system's true dynamic behavior. Thus new modeling and simulation methods are needed to reflect the discrete and stochastic features of biochemical systems on a cellular scale.

To include discreteness, the most accurate way to simulate the time evolution of a system of chemically reacting molecules is to do a molecular dynamics simulation tracking the positions and velocities of all the molecules and the occurrence of all chemical reactions when molecules physically collide with each other. But molecular dynamics simulations are generally too expensive to be practical except in the case of a relatively small number of molecules and even then only for very short time scales. Instead we consider the case where the dynamics of biochemical systems can be approximated by assuming that the reactant molecules are “well-stirred” such that their positions become randomized and need not be tracked in detail. When that is true, the state of the system can be defined simply by the instantaneous molecular populations of the various chemical species. The chemical reactions can be defined as events that change the state of the system following biochemical rules, changing the molecular populations by integer numbers.

This chapter discusses numerical methods that can be applied to simulate such systems, taking into account discreteness and stochasticity at the molecular level. The focus will be on two major simulation methods: Gillespie's stochastic simulation algorithm (SSA) and the tau-leaping method. Finally, a merger of these two methods, the hybrid SSA/tau-leaping method will be described.


McAdams H, Arkin A. Stochastic mechanisms in gene expression. Proc Natl Acad Sci. 1997 94(3):814–9. http://www.pnas.org/content/94/3/814.

Fedoroff N, Fontana W. Small numbers of big molecules. Ciência. 2002 297(5584):1129–31. http://www.sciencemag.org/content/297/5584/1129.

Samoilov M, Plyasunov S, Arkin AP. Stochastic amplification and signaling in enzymatic futile cycles through noise-induced bistability with oscillations. Proc Natl Acad Sci U S A. 2005 102(7):2310–5. http://www.pnas.org/content/102/7/2310.

Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976 22(4):403–34.

Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977 81(25):2340–61.

Michaelis L, Menten ML. Die Kinetik der Invertinwirkung. Biochem Z. 1913 49:333–69.

Hill AV. The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J Physiol. 1910 40(Suppl):iv–vii.

Hattne J, Fange D, Elf J. Stochastic reaction-diffusion simulation with MesoRD. Bioinformática. 2005 21(12):2923–4.

Andrews SS, Bray D. Stochastic simulation of chemical reactions with spatial resolution and single molecule detail. Phys Biol. 2004 1:137–51.

van Zon JS, ten Wolde PR. Green’s-function reaction dynamics: A particle-based approach for simulating biochemical networks in time and space. J Chem Phys. 2005 123(23):234910. http://scitation.aip.org/content/aip/journal/jcp/123/23/10.1063/1.2137716.

Novère NL, Shimizu TS. StochSim: modelling of stochastic biomolecular processes. Bioinformática. 2001 17:575–6.

von Smoluchowski M. Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen. Annalen der Physik. 1906 326(14):756–80. http://dx.doi.org/10.1002/andp.19063261405.

Gardiner CW, McNeil KJ, Walls DF, Matheson IS. Correlations in stochastic theories of chemical reactions. J Stat Phys. 1976 14:307–31.

Nicolis G, Prigogine I. Self-organization in nonequilibrium systems : from dissipative structures to order through fluctuations. New York: A Wiley-Interscience Publication 1977. http://opac.inria.fr/record=b1078628.

doi M. Stochastic theory of diffusion-controlled reaction. J Phys A Math General. 1976 9(9):1479. http://stacks.iop.org/0305-4470/9/i=9/a=009.

Keizer J. Nonequilibrium statistical thermodynamics and the effect of diffusion on chemical reaction rates. J Phys Chem. 1982 86(26):5052–67. http://dx.doi.org/10.1021/j100223a004.

Fange D, Berg OG, Sjöberg P, Elf J. Stochastic reaction-diffusion kinetics in the microscopic limit. Proc Natl Acad Sci. 2010 107(46):19820–5. http://www.pnas.org/content/107/46/19820.

McQuarrie DA. Stochastic approach to chemical kinetics. J Appl Probab. 1967 4(3):413–78. http://www.jstor.org/stable/3212214.

Gillespie DT. A rigorous derivation of the chemical master equation. Physica A: Stat Mech Appl. 1992 188(1–3):404–25.

Baras F, Mansour MM. Reaction-diffusion master equation: a comparison with microscopic simulations. Phys Rev E. 1996 54:6139–48. http://link.aps.org/doi/10.1103/PhysRevE.54.6139.

Isaacson SA. The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target. SIAM J Appl Math. 2009 70(1):77–111.

Erban R, Chapman SJ. Stochastic modelling of reaction-diffusion processes: algorithms for bimolecular reactions. Phys Biol. 2009 6(4):046001. http://stacks.iop.org/1478-3975/6/i=4/a=046001.

Hellander S, Hellander A, Petzold L. Reaction-diffusion master equation in the microscopic limit. Phys Rev E. 2012 85:042901. http://link.aps.org/doi/10.1103/PhysRevE.85.042901.

Li S, Brazhnik P, Sobral B, Tyson JJ. A Quantitative Study of the Division Cycle of Caulobacter crescentus, Stalked Cells. PLoS Comput Biol. 2008 4(1):e9. http://dx.plos.org/10.1371%252Fjournal.pcbi.0040009.

Li S, Brazhnik P, Sobral B, Tyson JJ. Temporal controls of the asymmetric cell division cycle in Caulobacter crescentus. PLoS Comput Biol. 2009 5(8):000463. http://dx.doi.org/10.1371%252Fjournal.pcbi.1000463.

Subramanian K, Paul MR, Tyson JJ. Potential role of a bistable histidine kinase switch in the asymmetric division cycle of Caulobacter crescentus. PLoS Comput Biol. 2013 9(9):003221. http://dx.doi.org/10.1371%2Fjournal.pcbi.1003221.

Subramanian K, Paul MR, Tyson JJ. Dynamical localization of DivL and PleC in the asymmetric Division cycle of Caulobacter crescentus: a theoretical investigation of alternative models. PLoS Comput Biol. 2015 11(7):004348.

Collier J, Murray SR, Shapiro L. DnaA couples DNA replication and the expression of two cell cycle master regulators. EMBO J. 2006 25(2):346–56.

Collier J, Shapiro L. Spatial complexity and control of a bacterial cell cycle. Curr Opin Biotechnol. 2007 18(4):333–40. http://www.sciencedirect.com/science/article/pii/S0958166907000894.

Holtzendorff J, Hung D, Brende P, Reisenauer A, Viollier PH, McAdams HH, et al. Oscillating global regulators control the genetic circuit driving a bacterial cell cycle. Ciência. 2004 304(5673):983–7. http://www.sciencemag.org/content/304/5673/983.

Quon KC, Yang B, Domian IJ, Shapiro L, Marczynski GT. Negative control of bacterial DNA replication by a cell cycle regulatory protein that binds at the chromosome origin. Proc Natl Acad Sci. 1998 95(1):120–5. http://www.pnas.org/content/95/1/120.

McGrath PT, Iniesta AA, Ryan KR, Shapiro L, McAdams HH. A dynamically localized protease complex and a polar specificity factor control a cell cycle master regulator. Célula. 2006 124(3):535–47. http://www.sciencedirect.com/science/article/pii/S0092867406000663.

Jenal U, Fuchs T. An essential protease involved in bacterial cell-cycle control. EMBO J. 1998 17(19):5658–69.

Iniesta AA, McGrath PT, Reisenauer A, McAdams HH, Shapiro L. A phospho-signaling pathway controls the localization and activity of a protease complex critical for bacterial cell cycle progression. Proc Natl Acad Sci. 2006 103(29):10935–40. http://www.pnas.org/content/103/29/10935.

Tsokos C, Perchuk B, Laub M. A dynamic complex of signaling proteins uses polar localization to regulate cell-fate asymmetry in Caulobacter crescentus. Developmental Cell. 2011 20(3):329–41.

Isaacson SA. A convergent reaction-diffusion master equation. J Chem Phys. 2013 139(5):054101. http://scitation.aip.org/content/aip/journal/jcp/139/5/10.1063/1.4816377.

Cao Y, Gillespie DT, Petzold LR. The slow-scale stochastic simulation algorithm. J Chem Phys. 2005 122(1):014116.

Haseltine EL, Rawlings JB. Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. J Chem Phys. 2002 117(15):6959–69.

Hindmarsh AC. ODEPACK, a systematized collection of ODE solvers. IMACS Trans Sci Comput Amsterdam. 1983 1:55–64.

Petzold LR. A Description of DASSL: a differential/algebraic system solver, proceeding of the 1st IMACS World Congress. Montreal. 19821:65-68.

Financiamento

This work was partially supported by the National Science Foundation award DMS-1225160, CCF-0953590, CCF-1526666, and MCB-1613741. In particular, the publication of this paper is directly funded from CCF-1526666, and MCB-1613741.

Disponibilidade de dados e materiais

Authors’ contributions

FL initially realized the simulation error described in this paper. FL, MC and YC then designed the toy model and analyzed the numerical error caused by Hill function simulation. Later SW joined to help the implementation of hybrid simulation. FL and MC together drafted the manuscript, and YC gave critical revisions on the writing. All authors have read and approved the final manuscript.

Interesses competitivos

Os autores declaram não ter interesses conflitantes.

Consentimento para publicação

All authors consent to publish this work through BMC’17.

Ethics approval and consent to participate

Sobre este suplemento

Este artigo foi publicado como parte de BMC Systems Biology Volume 11 Supplement 3, 2017: Selected original research articles from the Third International Workshop on Computational Network Biology: Modeling, Analysis, and Control (CNB-MAC 2016): systems biology. The full contents of the supplement are available online at http://bmcsystbiol.biomedcentral.com/articles/supplements/volume-11-supplement-3.


Assista o vídeo: Cómo hacer una gráfica de Michaelis-Menten. Enzimas parte XVI (Outubro 2022).