Em formação

Parâmetros para o modelo de potencial de ação de Fitzhugh-Nagumo


Sou um estudante do ensino médio que atualmente está cursando Biologia AP, então minha formação em biologia é bastante limitada. Estou interessado em modelar matematicamente os potenciais de ação.

Fazendo uma rápida pesquisa no Google, descobri o modelo Fitzhugh-Nagumo, um caso simplificado que modela a voltagem da célula durante um potencial de ação como uma equação diferencial de segunda ordem $ C frac {d ^ 2V} {dt ^ 2} + epsilon (V ^ 2-1) + frac {V} {L} = 0 $. Quero resolver essa equação numericamente, mas não consigo encontrar os valores numéricos de $ epsilon $, $ L $ ou $ C $ online. Quais são alguns parâmetros razoáveis ​​para este modelo para um neurônio humano?


Na verdade, este não é o modelo de Fitzhugh-Nagumo do neurônio, no entanto, está altamente relacionado a ele. Acredito que a equação que você tem aí é capaz de oscilações para quaisquer valores positivos de $ epsilon $, $ C $, $ L $. Eu gerei este gráfico com os valores dos parâmetros configurados em 1. Observe como ele é muito semelhante a uma onda senoidal.

Eu acredito que você estava tentando escrever a equação de Van der Pol aqui: $$ C frac {d ^ 2V} {dt ^ 2} + epsilon (V ^ 2-1) frac {dV} {dt} + frac { V} {L} $$ Observe a diferença sutil de multiplicar pela primeira derivada. Aqui está novamente uma imagem com todos os parâmetros definidos para 1

A equação de Fitzhugh-Nagumo pode ser encontrada na wikipedia aqui e na scholarpedia revisada por pares aqui. Para referência, é um modelo de dois estados descrito pela equação:

$$ frac {dV} {dt} = V-V ^ 3/3-W + I $$ $$ tau frac {dW} {dt} = V + a -bW $$

Um valor comum para os parâmetros são $ a = .7 $, $ b = .8 $, $ tau = 12,5 $ $ I $ pode ser qualquer coisa realmente, usei $ I = .5 $ para gerar o gráfico abaixo. Observe também que o modelo não aumentará se eu estiver abaixo de um determinado valor. A equação de Fitzhugh é um pouco mais geral, mas simplifica para o Van der pol sob a suposição de $ a = b = 0 $. Aqui está o gráfico com as constantes que descrevi.


Análise de parâmetros do modelo FitzHugh-Nagumo para uma simulação confiável

Derivado do modelo iônico pioneiro de Hodgkin-Huxley e devido à sua simplicidade e riqueza do ponto de vista da dinâmica não linear, o modelo FitzHugh-Nagumo (FHN) é um dos modelos de neurônios / células cardíacas simplificados de maior sucesso. Existem muitas variações do modelo FHN original. Embora esses modelos de tipo FHN ajudem a enriquecer a dinâmica do modelo FHN, os parâmetros usados ​​nesses modelos costumam estar em condições tendenciosas. Os resultados relacionados seriam questionáveis. Portanto, neste estudo, o objetivo é encontrar os limites dos parâmetros para um dos modelos FHN comumente usados, a fim de fornecer um melhor ambiente de simulação. Os resultados mostraram inicialmente que intervalo de tempo e tolerância de integração inadequados na solução numérica do modelo FHN podem fornecer alguns resultados tendenciosos que tornariam algumas publicações questionáveis. Em seguida, são apresentados os limiares dos parâmetros α, γ e ε. α controla a dinâmica global do FHN. α & gt 0, a célula está no modo refratário α & lt 0, a célula é excitável. ε controla a morfologia principal do potencial de ação gerado e tem relação com o período (P = 3,065 × αα, γ (-0,8275) + 4,397). Para mostrar oscilações de relaxamento com FHN, ε deve ser menor que 0,0085. 7 mal influencia o potencial de ação, apresentou relação linear com o período e duração do potencial de ação. Embora α & lt 0,1, ε & lt 0,0085, não haja um limite definido para γ, valores menores são recomendados.


Akcelik, V. (2002). Métodos multiescala de Newton-Krylov para propagação de onda acústica inversa. Ph.D. tese, Carnegie Mellon University, Pittsburgh, Pennsylvania.

Akcelik, V., Biros, G., Ghattas, O., Hill, J., Keyes, D., & amp van Bloemen Waanders, B. (2006). Algoritmos paralelos para otimização restrita por PDE. In: M. Heroux, P. Raghaven e H. Simon (Eds.), Fronteiras da computação paralela. SIAM.

Argentina, M., Coullet, P., & amp Krinsky, V. (2000). Colisões frontais de ondas em um sistema FitzHugh – Nagumo excitável: uma transição da aniquilação de ondas para o comportamento de ondas clássico. Journal of Theoretical Biology, 205, 47–52.

Balay, S., Buschelman, K., Gropp, W. D., Kaushik, D., Knepley, M. G., McInnes, L. C. & amp et al. (2007). PETSc Homepage. http://www.mcs.anl.gov/petsc.

Banks, H. T., & amp Kunisch, K. (1989). Técnicas de estimativa para sistemas de parâmetros distribuídos. Boston, MA: Birkhäuser.

Benzi, M., Haber, E., & amp Hanson, L. (2006). Algoritmos multinível para métodos de pontos interiores em grande escala em otimização limitada limitada. Relatório Técnico TR-2006-002-A, Departamento de Matemática e Ciência da Computação, Emory University. 16p.

Bernus, O., Verschelde, H., & amp Panfilov, A. V. (2002). Modelos iônicos modificados de tecido cardíaco para cálculos eficientes em grande escala. Física em Medicina e Biologia, 47, 1947–1959.

Biegler, L., Ghattas, O., Heinkenschloss, M., & amp Bloemen-Waanders, B. v. (Eds.) (2003). Otimização com restrição de PDE em grande escala. No Notas de aula em Ciência e Engenharia da Computação. Berlim: Springer-Verlag.

Biros, G., & amp Ghattas, O. (2005a). Métodos paralelos de Lagrange – Newton – Krylov – Schur para otimização restrita por PDE. Parte I: O solucionador Krylov-Schur. SIAM Journal on Scientific Computing, 27, 687–713.

Biros, G., & amp Ghattas, O. (2005b). Métodos paralelos de Lagrange – Newton – Krylov – Schur para otimização restrita por PDE. Parte II: O solucionador Lagrange – Newton e sua aplicação para o controle ideal de fluxos viscosos estáveis. SIAM Journal on Scientific Computing, 27, 714–739.

Brooks, D. H., Ahmad, G. F., MacLeod, R. S., & amp Maratos, G. M. (1999). Eletrocardiografia inversa por imposição simultânea de múltiplas restrições. IEEE Transactions on Biomedical Engineering, 46, 3–18.

Bub, G., Shrier, A., & amp Glass, L. (2002). Geração de onda espiral em meios excitáveis ​​heterogêneos. Cartas de revisão física, 88, 058101.

Chen, X. & amp Oshita, Y. (2006). Periodicidade e exclusividade de minimizadores globais de um funcional de energia contendo uma interação de longo alcance. SIAM Journal on Mathematical Analysis, 37, 1299–1332.

Cheng, L. K., Bodley, J. M., & amp Pullan, A. J. (2003). Comparação de formulações baseadas em potencial e ativação para o problema inverso da eletrocardiologia. IEEE Transactions on Biomedical Engineering, 50, 11–22.

Colli-Franzone P., & amp Pavarino, L. F. (2004). Um solucionador paralelo para sistemas de reação-difusão em eletrocardiologia computacional. Modelos e métodos matemáticos em ciências aplicadas, 14, 883–911.

Courtemanche, M., Skaggs, W., & amp Winfree, A. T. (1990). Circulação de potencial de ação tridimensional estável no modelo FitzHugh – Nagumo. Physica D, 41, 173–182.

Cox, S. J. (2006). Um método adjunto para localização de canal. Medicina Matemática e Biologia, 23, 139–152.

Cox, S. J., & amp Griffith, B. E. (2001). Recuperando propriedades quase ativas de neurônios dendríticos de gravações de potencial duplo. Journal of Computational Neuroscience, 11, 95–110.

Cox S. J., & amp Ji, L. (2003) Discerning iônicas correntes e sua cinética de dados de impedância de entrada. Boletim de Biologia Matemática, 63, 909–932.

Cox, S. J., & amp Wagner, A. (2004). Sobredeterminação lateral do sistema FitzHugh – Nagumo. Problemas inversos, 20, 1639–1647.

Dauby, P. C., Desaive, T., & amp Croisier, H. (2006). Ondas estacionárias no modelo FitzHugh – Nagumo de atividade elétrica cardíaca. Revisão Física E, 73, 021908.

Davidenko, J. M., Pertsov, A. V., Salomonsz, R., Baxter, W., & amp Jalefe, J. (1992). Ondas espirais estacionárias e flutuantes de excitação em músculo cardíaco isolado. Nature, 355, 349–351.

Dennis, Jr J. E., & amp Schnabel, R. B. (1996). Métodos numéricos para otimização irrestrita e equações não lineares. No Clássicos em Matemática Aplicada. Filadélfia: SIAM.

Eisenstat S. C., & amp Walker, H. F. (1994). Métodos de Newton inexatos globalmente convergentes, 4, 393–422.

Elmer, C. E., & amp Van Vleck, E. S. (2005). Equações de FitzHugh – Nagumo espacialmente discretas. SIAM Journal on Applied Mathematics, 65, 1153–1174.

Engl, H. W., Hanke, M., & amp Neubauer, A. (1996). Regularização de problemas inversos. Dordrecht: Kluwer.

FitzHugh, R. (1961). Impulsos e estados fisiológicos em modelos teóricos de membrana nervosa. Biophysical Journal, 1, 445–466.

Franzone, P. C., Pavarino, L. F., & amp Taccardi, B. (2005). Simulação de padrões de excitação, repolarização e duração do potencial de ação com modelos de bidomain e monodomínio cardíacos. Biociências matemáticas, 197, 35–66.

Gao, W., & amp Wang, J. (2004). Existência de frentes de onda e impulsos para as equações de FitzHugh – Nagumo. Análise não linear, 57, 667–676.

Haber, E., Ascher, U., & amp Oldenburg, D. (2000). Sobre técnicas de otimização para resolver problemas inversos não lineares. Problemas inversos, 16, 1263–1280.

Hansen, P. C., & amp O’Leary, D. P. (1993). O uso da curva L na regularização de problemas discretos mal colocados. SIAM Journal on Scientific Computing, 14, 1487–1503.

Hoffman, D. A., Magee, J. C., Colbert, C. M., & amp Johnston, D. (1997). K + Regulação do canal de propagação do sinal em dendritos de neurônios piramidais do hipocampo. Nature, 387, 869–875.

Isakov, V. (1998). Problemas inversos para equações diferenciais parciais. Nova York: Springer-Verlag.

Knoll, D. A., & amp Keyes, D. E. (2004). Métodos de Newton-Krylov livres de Jacobianos: Uma pesquisa de abordagens e aplicações. Journal of Comparative Physiology, 193, 357–397.

Krupa, M., Sandstede, B., & amp Szmolyan, P. (1997). Ondas rápidas e lentas na equação de FitzHugh – Nagumo. Journal of Differential Equations, 133, 49–97.

MacLeod, R. S., & amp Brooks, D. H. (1998). Progresso recente em problemas inversos em eletrocardiologia. Revista IEEE Engenharia em Medicina e Biologia, 17, 73–83.

Moreau-Villéger, V., Delingette, H., Sermesant, M., Ashikaga, H., Faris, O., & amp McVeigh, E., et al. (2006). Criação de mapas de condutividade aparente local do epicárdio com um modelo eletrofisiológico 2D do coração. Transações IEEE em Engenharia Biomédica 53(8), 1457–1466.

Murillo, M., & amp Cai, X.-C. (2004). Um algoritmo paralelo totalmente implícito para simular a atividade elétrica não linear do coração. Álgebra Linear Numérica com Aplicações, 11, 261–277.

Murray, J. D. (1993). Biologia matemática. Berlim: Springer

Nagumo, J., Arimoto, S., & amp Yoshizawa, S. (1962). Uma linha de transmissão de pulso ativa que simula o axônio do nervo. Anais do Instituto de Engenheiros de Rádio, 50, 2061–2070.

Nii, S. (1997). Estabilidade de viajar soluções de ondas frontais (múltiplas traseiras) das equações de FitzHugh – Nagumo. SIAM Journal on Mathematical Analysis, 28, 1094–1112.

Nocedal, J. & amp Wright, S. J. (1999). Otimização numérica. Nova York: Springer-Verlag.

Patel, S. G., & amp Roth, B. J. (2005). Solução aproximada para as equações de bidomain para problemas de eletrocardiograma. Revisão Física E, 72, 051931.

Pennacchio, M., Savare, G., & amp Franzone, P. C. (2006). Modelagem multiescala para a atividade bioelétrica do coração. SIAM Journal on Mathematical Analysis, 37, 1333–1370.

Pernarowski, M. (2001). Controlabilidade de sistemas excitáveis. Boletim de Biologia Matemática, 63, 167–184.

Petersson, J. H. (2005). Sobre a existência global de sistemas parabólicos semilineares. Análise não linear, 60, 337–347.

Riccio, M. L., Koller, M. L., & amp Gilmour, P. F. (1999). Restituição elétrica e organização espaço-temporal durante a fibrilação ventricular. Circulation Research, 84, 955– 963.

Roth, B. J. (2004). Art Winfree e o modelo de bidomain de tecido cardíaco. Journal of Theoretical Biology, 230, 445–449.

Saad, Y. (2003). Métodos iterativos para sistemas lineares esparsos, 2ª ed. Filadélfia: SIAM.

Scott, A. C. (1975). A eletrofísica de uma fibra nervosa. Resenhas de Física Moderna, 47, 487–533.

Shahidi, V., Savard, P., & amp Nadeau, R. (1994). Problema direto e inverso da eletrocardiografia: Modelagem e recuperação de potenciais epicárdicos em humanos. IEEE Transactions on Biomedical Engineering, 41, 249–256.

Sneyd, J., Dale, P. D., & amp Duffy, A. (1998). Ondas viajantes em sistemas protegidos: aplicações às ondas de cálcio. SIAM Journal on Applied Mathematics, 58, 1178–1192.

Suckley, R., & amp Biktashev, V. N. (2003). Comparação de assintóticos de excitabilidade do coração e nervo. Revisão Física E, 68, 011902.

Tsai, J.-C., & amp Sneyd, J. (2005). Existência e estabilidade de ondas viajantes em sistemas protegidos. SIAM Journal on Applied Mathematics, 66, 237–265.

Vanier, M. C., & amp Bower, J.-M. (1999). Uma pesquisa comparativa de métodos automatizados de busca de parâmetros para modelos neurais compartimentais. Journal of Computational Neuroscience, 7, 149–171.

Vogel, C. R. (1996). Não convergência do método de seleção do parâmetro de regularização da curva L. Problemas inversos, 12, 535–547.

Vogel, C. R. (2002). Métodos computacionais para problemas inversos. In Frontiers in Matemática Aplicada. Filadélfia: SIAM.

Weiss, J. N., Garfinkel, A., Karagueuzian, H. S., Qu, Z., & amp Chen, P. S. (1999). Caos e a transição para a fibrilação ventricular - uma nova abordagem para avaliação de drogas antiarrítmicas. Circulação, 99, 2819–2826.

Willms, A. R., Baro, D. J., Harris-Warrick, R. M., & amp Guckenheimer, J. (1999). Um método de estimativa de parâmetro aprimorado para modelos de Hodgkin-Huxley. Journal of Computational Neuroscience, 6, 145–168.

Winfree, A. T. (1990). Soluções semelhantes a partículas estáveis ​​para as equações de ondas não lineares de meios excitáveis ​​tridimensionais. Revisão SIAM, 32, 1–53.

Yamada, H. & amp Nozaki, K. (1990). Interação de pulsos em sistemas dissipativos. Equações de FitzHugh – Nagumo. Progresso da Física Teórica, 84, 801–809.


Resultados

Modelo de monodomínio de FitzHugh-Nagumo

Nesta seção, alguns exemplos numéricos são dados para verificar a eficiência do método GS proposto. Nos neurônios, vários tipos de canais iônicos podem influenciar o potencial de membrana. Os canais iônicos dependentes de voltagem são controlados pelo potencial de membrana, enquanto o potencial de membrana é influenciado por esses mesmos canais iônicos, o que causa loops de feedback que permitem dinâmicas temporais complexas, incluindo oscilações e eventos regenerativos, como AP. Primeiramente, o método GS é usado para resolver um modelo bidimensional de FitzHugh-Nagumo em regiões retangulares e circulares 40,53. O modelo FitzHugh-Nagumo é


Compreenda a dinâmica do modelo FitzHugh-Nagumo com um aplicativo

Em 1961, R. Fitzhugh (Ref. 1) e J. Nagumo propuseram um modelo para emular o sinal atual observado em células excitáveis ​​de um organismo vivo & # 8217s. Isso ficou conhecido como o modelo FitzHugh-Nagumo (FN) de neurociência matemática e é uma versão mais simples do modelo Hodgkin-Huxley (HH) (Ref. 2), que demonstra as correntes de pico nos neurônios. Na postagem do blog de hoje & # 8217s, examinaremos a dinâmica do modelo FN criando um aplicativo interativo no software COMSOL Multiphysics®.

Potencial de ação de uma célula nervosa e # 8217s

As células nervosas são separadas da região extracelular por uma membrana de bicamada lipídica. Quando as células não estão conduzindo um sinal, há uma diferença potencial de cerca de -70 mV através da membrana. Esta diferença é conhecida como potencial de repouso da célula & # 8217s. Íons minerais, como sódio e potássio, e íons de proteína carregados negativamente, contidos na célula, mantêm o potencial de repouso. Quando a célula recebe um estímulo externo, seu potencial aumenta em direção a um valor positivo, um processo conhecido como despolarização, antes de cair novamente para o potencial de repouso, chamado repolarização.


Gráfico do potencial de ação de uma célula.

Em um exemplo, a concentração dos íons de sódio em repouso é muito maior na região extracelular do que dentro da célula. A membrana contém canais fechados que permitem seletivamente a passagem de íons através deles. Quando a célula é estimulada, os canais de sódio se abrem e há um fluxo de íons sódio para dentro da célula. Esta corrente de sódio & # 8220 & # 8221 aumenta o potencial da célula, resultando na despolarização. No entanto, como os portões do canal são acionados por voltagem, os portões de sódio fecham depois de um tempo. Os canais de potássio então se abrem e uma corrente de potássio de saída flui, levando à repolarização da célula.

Hodgkin e Huxley explicaram este mecanismo de geração de potencial de ação por meio de equações matemáticas (Ref. 2). Embora tenha sido um grande sucesso na modelagem matemática de fenômenos biológicos, o modelo completo de Hodgkin-Huxley é bastante complicado. Por outro lado, o modelo de FitzHugh-Nagumo é relativamente simples, consistindo em menos parâmetros e apenas duas equações. Um é para a quantidade V, que imita o potencial de ação, e o outro é para a variável C, que modula V.

Hoje, vamos nos concentrar no modelo FN, enquanto o modelo HH será um tópico de discussão em um momento posterior.

Definindo o modelo FitzHugh-Nagumo

As duas equações do modelo FN são

O parâmetro I corresponde a uma excitação enquanto uma e b são os parâmetros de controle do modelo. A evolução de C é mais lento do que a evolução de V devido ao parâmetro ε multiplicar tudo no lado direito da segunda equação. Os pontos fixos das equações do modelo FN são as soluções do seguinte sistema de equações,

o V-nullcline e C-nullcline são as curvas V-V ^ 3/3-W + I e V + a-bW, respectivamente, no plano VW. Observe que o V-nullcline é uma curva cúbica no plano VW e o C-nullcline é uma linha reta. A inclinação da reta V + a-bW é controlada de tal forma que as nulas se cruzam em um único ponto, tornando o sistema & # 8217s único ponto fixo.

O parâmetro I simplesmente muda V-nullcline para cima ou para baixo. Assim, a mudança de I modula a posição do ponto fixo de modo que diferentes valores de I façam com que o ponto fixo fique à esquerda, no meio ou à direita da curva V-V ^ 3/3-W + I.

Compreendendo a dinâmica do modelo FN

Para simular o que acontece quando o ponto fixo está em cada região, podemos usar o DAE global interface incluída no pacote básico do COMSOL Multiphysics.

o V-nullcline é mostrado na figura abaixo em uma cor verde. Na região acima deste nullcline frac

& lt0, enquanto na região abaixo é positivo. o C-nullcline é mostrado em vermelho na região à direita desta linha reta, frac
& gt0, e à esquerda, frac
& lt0.

Vamos primeiro examinar o que acontece se o ponto fixo estiver no lado direito, região 3, do V-nullcline. Diremos isso quando t, representando o tempo, é igual a zero, ambos V e C também estão em zero. Neste caso, ambos frac

e frac
são positivos no ponto de partida e próximo a ele e, portanto, ambos mudam com o passar do tempo. Mas desde V evolui mais rápido do que C, V aumenta rapidamente enquanto C permanece praticamente inalterado. Na figura, podemos ver que isso resulta em uma parte quase horizontal da curva V-W.

Conforme a curva se aproxima do V-nullcline, a taxa de mudança de V desacelera e C torna-se mais proeminente. Desde frac

ainda é positivo, C deve aumentar, e a curva se move para cima. O ponto fixo então atrai esta curva e a evolução termina no ponto fixo.


Gráfico do plano VW quando o ponto fixo está no lado direito do V-nullcline.

Se o ponto fixo está no meio, região 2, então o que discutimos até agora é verdadeiro. A diferença é que uma vez que a curva vai além do joelho direito do V-nullcline, frac

torna-se negativo e V decai rapidamente. Enquanto se move para a esquerda, a curva cruza a nuliclina vermelha da direita para a esquerda. Deste ponto em diante, enquanto ambos V e C diminuir, a evolução de V domina e a curva torna-se horizontal novamente.

Isso continua até que a curva atinja a parte esquerda do V-nullcline. A curva começa a abraçar o V-nullcline e inicia uma jornada lenta para baixo. Quando toca o joelho esquerdo do V-nullcline, ele se move rapidamente em direção à parte direita do V-nullcline. Observe que esse movimento nunca atinge o ponto fixo e, portanto, continua se repetindo, o que podemos ver no gráfico abaixo.


Gráfico do plano VW quando o ponto fixo está na região central do V-nullcline.

Isso nos deixa com um último caso para discutir & mdash quando o ponto fixo está na parte esquerda, Região 1, do V-nullcline. Os resultados devem ser semelhantes ao gráfico a seguir. Observe que as análises que realizamos anteriormente são transferidas.


Gráfico do plano VW quando o ponto fixo está no lado esquerdo do V-nullcline.

Simulando o modelo FN usando o construtor de aplicativos

Para explorar a rica dinâmica do modelo FN descrito acima, precisamos alterar repetidamente várias entradas sem alterar o modelo subjacente. Como tal, é desejável uma interface de usuário que nos permite alterar facilmente os parâmetros do modelo, realizar a simulação e analisar os novos resultados sem ter que navegar na estrutura de árvore do Model Builder para realizar essas várias ações.

Para fazer isso, podemos recorrer ao Application Builder. Esta plataforma nos permite criar um aplicativo de simulação fácil de usar que expõe todos os aspectos essenciais do modelo, enquanto mantém o resto nos bastidores. Com este aplicativo, podemos alterar rapidamente os parâmetros por meio de uma interface amigável e estudar os resultados usando figuras estáticas e animações. O aplicativo também torna mais fácil para os alunos compreenderem a dinâmica do modelo FN & # 8217s sem ter que se preocupar em criar um modelo.

Os parâmetros importantes do modelo FN, ou seja, a, b, ε, e eu, são exibidos no aplicativo & # 8217s Parâmetros do modelo seção. Os painéis gráficos exibem várias quantidades de interesse, como a forma de onda para V e C. Exibimos o diagrama do plano de fase no painel superior direito junto com o V& # 8211 e C-nullclines. A posição do ponto fixo é facilmente identificável nesse gráfico. Assim que a simulação for concluída, podemos animar as trajetórias de tempo escolhendo a opção de animação na barra de ferramentas da faixa de opções. Para obter um resumo dos parâmetros e resultados da simulação, podemos selecionar o Relatório de Simulação botão.


App que mostra a dinâmica do modelo FitzHugh-Nagumo quando o ponto fixo está na Região 2.

Podemos reproduzir facilmente os casos descritos na seção anterior com nosso aplicativo. A imagem acima, por exemplo, mostra o que acontece quando o ponto fixo está na região 2. Podemos facilmente mover o ponto fixo para a região 1 ou 3, tornando o atual 0,1 ou 2,5, respectivamente. Observe que quaisquer outros parâmetros no aplicativo também podem ser alterados para ver se outras tendências interessantes surgem.

O aplicativo que apresentamos aqui é apenas um exemplo do que você pode criar com o Application Builder. O design do seu aplicativo, desde o layout até os parâmetros incluídos, depende de você. A flexibilidade do Application Builder permite adicionar tanta complexidade quanto necessário, em parte graças aos métodos do Method Editor for Java®. Em uma postagem do blog de acompanhamento, criaremos um aplicativo para ilustrar a dinâmica do modelo mais complicado de HH. Fique ligado!

Referências

Mais recursos sobre como usar o Application Builder

  • Veja como combinar modelagem baseada em equações e aplicativos nesta postagem do blog: Implementando a forma fraca com um aplicativo COMSOL®
  • Para descobrir como construir seu próprio aplicativo, confira nossa série de vídeos de introdução ao construtor de aplicativos

Oracle e Java são marcas registradas da Oracle e / ou de suas afiliadas.


Discussão

Neste artigo, propusemos uma nova técnica de modelagem denominada modelo ressonante (RM) para descrever o comportamento elétrico de células biológicas observadas in vivo ou in vitro por meio da representação simplificada in silico de sua morfologia de potencial de ação. O modelo prevê propriedades experimentais e características de células biológicas, incluindo (1) morfologia SAN AP humana e de coelho, (2) efeitos de bloqueio de corrente, (3) arrastamento de frequência em células autorítmicas, (4) morfologia AP de miócito cardíaco humano e (5) ) Morfologia AP do neurônio gigante de lula. Validamos a morfologia AP produzida pelo RM da célula autorrítmica humana e do miócito cardíaco humano. Além disso, a frequência de arrastamento exibida pelo RM de SAN humano foi validada em relação aos resultados experimentais. Os dois componentes do RM são o gerador de forma de onda e o controlador de estado. O gerador Waveshape produz a forma de onda AP básica e consiste em osciladores. Cada oscilador representa um harmônico da série de Fourier. RM é capaz de produzir formas de onda AP muito diferentes e complexas com a mesma estrutura de modelagem que vimos no caso de células e neurônios SAN. Ambas as células possuem características de AP marcadamente diferentes, com um neurônio disparando muito rápido com pico acentuado em comparação com a célula SAN de ritmo lento. No entanto, RM capturou essas morfologias AP apenas variando o número de osciladores. Isso adiciona certeza computacional, pois o comportamento numérico de cada oscilador será o mesmo. Para incorporar várias propriedades dinâmicas de uma célula, o controlador de estado foi projetado com a ajuda de um stateflow ® que inclui várias funções. Essas funções são avaliadas apenas uma vez na mudança de alguma condição de trabalho e isso resulta em um aumento mínimo no custo computacional. Os resultados (Fig. 9) mostram que com o aumento do número de harmônicos, a precisão do RM também aumenta. No entanto, mais osciladores resultarão em um aumento no tempo de execução e, portanto, há uma compensação entre fidelidade e tempo de execução. Além disso, a variação paramétrica da morfologia AP pode ser usada para modelar o comportamento de novos medicamentos. A administração de um novo medicamento pode levar ao bloqueio de uma ou mais correntes iônicas. Em relação a γ na Eq (4). γ poderia ser proporcional à concentração do novo medicamento. Embora a droga possa influenciar vários canais iônicos, o impacto na morfologia da forma de onda pode ser parametrizado por uma única variável (γ).

Para fins de comparação, analisamos quantitativamente a complexidade de modelos eletrofisiológicos detalhados e genéricos com relação ao MR. Os resultados são apresentados na Tabela 3. Os resultados mostram o número médio de pontos de dados computados por segundo pela Unidade Central de Processamento (CPU) para diferentes modelos. Todos os modelos foram simulados em tempo discreto no MATLAB em várias etapas de tempo oitenta vezes por 5000 ms. Os modelos detalhados de eletrofisiologia [21, 22] com mais de 27 variáveis ​​de estado são vinte e cinco vezes complexos do que o RM com doze termos de Fourier (FT). Além disso, os modelos de eletrofisiologia detalhados requerem o uso de tamanhos de etapa menores nos cálculos numéricos em comparação com o RM. Reconhecemos o fato de que os modelos detalhados de eletrofisiologia são capazes de capturar comportamentos complexos. No entanto, para aplicações específicas, como simulações de tecido em grande escala envolvendo milhões de pontos de grade, uma técnica alternativa que resulte em acelerações de alta velocidade é necessária. Os resultados indicam que o RM é duas vezes mais rápido que o modelo mais simples de três variáveis ​​do potencial de ação cardíaco [38]. Esta comparação mostra que o RM melhora a tratabilidade computacional enquanto reproduz com precisão a maioria das características do potencial de ação.


3. Resultados

A fim de obter um melhor entendimento de como o algoritmo de gradiente atua, suas vantagens e limitações, nós o aplicamos a três cenários distintos: o disparo de um único potencial de ação em um sistema monoestável, o início de disparos repetitivos em um sistema biestável e a supressão de disparos repetitivos em um sistema biestável. Prosseguimos com dois dos sistemas de modelos neuronais mais clássicos: o modelo Hodgkin-Huxley como nosso sistema monoestável e o sistema FitzHugh-Nagumo como o sistema biestável. Embora o sistema FitzHugh-Nagumo tenha sido usado principalmente em sistemas neuronais, ele também tem aplicações mais amplas em outros sistemas biológicos (Aliev & # x00026 Panfilov 1996 Kawato & # x00026 Suzuki 1980).

3.1 Modelo Hodgkin-Huxley de Excitação Neuronal

Como podemos ver na Figura 2, o algoritmo de gradiente começa com um estímulo gerado aleatoriamente e, em algumas gerações, a forma aproximada do estímulo ideal é vista. Em 30 iterações, o estímulo ideal é revelado, com poucas melhorias adicionais na norma L 2 após as iterações subsequentes. Uma observação interessante que notamos é que a primeira iteração geralmente produz um estímulo com um índice de desempenho inferior. Isso se deve ao fato de que, como estamos gerando o estímulo aleatoriamente, ele muito provavelmente falha em atender à condição terminal. Assim, o algoritmo primeiro altera o estímulo para estar mais perto de atender a condição terminal, antes de começar a melhorar o índice de desempenho.

O algoritmo de gradiente molda os estímulos aleatórios em direção a uma forma de onda ideal. Aqui, mostramos a progressão das soluções à medida que o algoritmo de gradiente começa com um estímulo de ruído branco e encontra a solução mais eficiente energeticamente em 100 iterações. O painel superior mostra a trajetória da norma L 2 em 100 iterações. Os seis painéis mostram a forma de onda em evolução na 1ª, 5ª, 10ª, 20ª e 30ª iteração do algoritmo. Há muito pouca melhora entre a 30ª e a 100ª iteração, conforme visto na trajetória da norma L 2 (painel superior). Por convenção, a corrente positiva despolariza

A Figura 3 mostra a forma de onda do estímulo do algoritmo de gradiente em 100 iterações e o potencial de ação que causou. Observe que o potencial de membrana V na Figura 3 é deslocado do modelo HH por & # x0221260 mV, que também é consistente com o uso moderno do modelo Hodgkin-Huxley. O modelo original definiu arbitrariamente o potencial de repouso em 0 mV, enquanto os neurônios, na verdade, descansam em um estado hiperpolarizado. O resultado do algoritmo de gradiente tinha uma norma L 2 de 15,5 & # x000b5J / cm 2. Como um ponto de comparação, usamos um estímulo de amplitude constante de 25 ms e reduzimos a amplitude para baixo até que apenas criasse um potencial de ação. Descobrimos que, usando uma amplitude de 2,255 & # x000b5A / cm 2, o neurônio dispararia um potencial de ação a 11 ms. Em seguida, reduzimos a duração do estímulo até que ele não disparasse mais um potencial de ação e descobrimos que quando a amplitude era 2,255 & # x000b5A / cm 2 e a duração era 9,7 ms, o neurônio mal dispararia um potencial de ação. A norma L2 deste pulso retangular mal supra-limiar foi de 49 & # x000b5J / cm 2. Assim, a forma de onda gerada com esta técnica mostrou uma grande redução da energia necessária para fazer com que o neurônio disparasse um único potencial de ação.

O estímulo ideal derivado do algoritmo de gradiente (parte superior) dispara um único potencial de ação no modelo de Hodgkin-Huxley (parte inferior)

À medida que buscamos obter uma melhor compreensão da aplicação do algoritmo & # x02019s, examinamos o efeito de que a duração do estímulo afetaria a norma de estímulo & # x02019 L 2 ideal. Nós mudamos tf, mas manteve o resto dos parâmetros do sistema constantes. Isso significava que não apenas o estímulo diminuía em duração, mas também o potencial de ação ocorria mais cedo. A Figura 4 mostra o gráfico que ilustra como a duração do estímulo afetou o estímulo ideal & # x02019 L 2 -norm. Como podemos ver no gráfico, quando a duração do estímulo foi aumentada, menores L 2 -norms puderam ser alcançados, até certo ponto. Além de 25 milissegundos, as melhorias na norma L 2 são mínimas.

A duração do estímulo mais longa fornece estímulos mais eficientes energeticamente. Uma vez que a duração do estímulo se estende além de um certo ponto, não há nenhuma melhoria adicional na eficiência energética

A Figura 5 mostra a forma das formas de onda ideais sob as diferentes condições de quando o sistema cruza o limite especificado. É interessante notar que quando um potencial de ação é desejado mais cedo (menos de 7 ms), a forma de onda ideal é monofásica, ao passo que tempos de potencial de ação posteriores (mais de 7 ms) são atingidos de maneira ideal com formas de onda bifásicas. Em um estudo recente, Clay, Forger e Paydarfar (2012) explicaram que a fase de hiperpolarização dos estímulos ótimos é útil para remover uma pequena quantidade de inativação de sódio, permitindo assim que uma fase de despolarização menos energética ainda elicie um potencial de ação. Esta figura mostra que isso é realmente verdade, mas apenas quando o tempo desejado para o potencial de ação é longo o suficiente. Se o tempo para o potencial de ação for mais curto, a fase de hiperpolarização é reduzida e ela desaparece completamente se o tempo para o potencial de ação for muito curto. De uma perspectiva de projeto, isso pode sugerir que em sistemas neuronais que requerem elicitação rápida de um potencial de ação, as correntes pós-sinápticas excitatórias seriam muito mais prevalentes do que em sistemas neuronais que são mais suscetíveis à elicitação retardada de potenciais de ação.

As formas de onda ideais mudam de um estímulo monofásico para um estímulo bifásico conforme a quantidade de tempo para o potencial de ação aumenta

Examinamos a robustez do algoritmo de gradiente para estimativas iniciais geradas aleatoriamente do estímulo ideal. Como tal, criamos um conjunto de sementes geradas aleatoriamente com amplitudes variadas (em resolução de 0,1 ms) e rastreamos para cada iteração a norma L 2 e a distância em que o ponto final calculado estava do estado terminal desejado. A Figura 6 mostra um gráfico da trajetória da norma L 2 ao longo das primeiras 30 iterações. Embora cada trajetória tenha começado em locais diferentes, todas terminaram em soluções com aproximadamente a mesma norma L 2. Existe alguma variação no erro devido às condições terminais, mas o sistema é extremamente sensível e, portanto, quaisquer pequenas alterações no estímulo resultarão em pequenas variações na distância do ponto final às condições terminais.

O algoritmo de gradiente é robusto ao estímulo inicial. As trajetórias de energia de estímulo de 100 sementes diferentes com amplitudes diferentes são mostradas

3.2 O Modelo FitzHugh-Nagumo: Quiescência para Disparos Repetitivos

Os valores da norma L 2 para todas as 10 execuções de cada um dos 68 estudos computacionais são mostrados na Figura 7. Para cada uma dessas execuções, o sistema começa no ponto fixo quiescente e termina em uma fase especificada, & # x003d5, conforme definido anteriormente. O algoritmo de gradiente falhou em convergir dentro do limite definido de 1000 iterações para algumas das execuções e nós as marcamos de acordo. Este gráfico revela um agrupamento de valores baixos da norma L 2 em torno de & # x003d5 = 0,58, uma falta de convergência no intervalo 0,05 & # x0003c & # x003d5 & # x0003c 0,5 e uma multiplicidade de soluções onde & # x003d5 & # x0003e 0,9 ou & # x003d5 & # x0003c 0,05.

Apenas certas regiões de fase convergiram durante a transição de quiescência para disparos repetitivos. Os pontos que convergiram são marcados em preto, enquanto aqueles que falharam são marcados em vermelho. Para uma definição da fase, consulte os Métodos 2.4 e 2.5

Para começar, notamos que a transição ideal do estado quiescente para o estado de disparo repetitivo é quando a condição terminal está em torno de & # x003d5 = 0,58. Na verdade, há um pequeno intervalo, 0,55 & # x0003c & # x003d5 & # x0003c 0,65, onde o algoritmo de gradiente produz valores da norma L 2 consistentemente baixos. Isso mostra uma trajetória particular com a & # x0201centrance & # x0201d mais ideal no estado de disparo repetitivo a partir do estado quiescente. Nessas fases, parece haver apenas um extremo no espaço de solução.

Em segundo lugar, existe uma grande variedade de fases, nas quais o algoritmo gradiente é incapaz de convergir para uma solução com o número predeterminado de iterações. Aumentamos o número de iterações para 5000, mas o algoritmo de gradiente ainda não conseguiu convergir para uma solução no intervalo 0,05 & # x0003c & # x003d5 & # x0003c 0,4, o que corresponde a mudanças rápidas nas variáveis ​​de estado durante o potencial de ação. Como estamos usando um algoritmo de gradiente de primeira ordem, conforme o algoritmo se aproxima dessas soluções, pequenas mudanças no estímulo podem fazer com que o algoritmo ultrapasse sua estimativa do ótimo, levando à perda de convergência. Acreditamos que existe uma solução porque a condição terminal está em um ciclo limite de estado estacionário. Se pegarmos a solução quando & # x003d5 = 0,58 e preenchermos o final com zeros, deveremos ser capazes de encontrar um estímulo que nos leve a um ângulo de fase entre 0,05 e 0,4. Isso sugere que talvez o conjunto de estimativas iniciais que permitem a convergência para uma solução ótima local, também conhecido como a região de convergência, seja pequeno em comparação com a região de convergência que vimos para o Hodgkin-Huxley monoestável ou o FitzHugh biestável -Nagumo quando a condição do terminal estiver na faixa de 0,55 & # x0003c & # x003d5 & # x0003c 0,65.

Finalmente, notamos vários ótimos locais encontrados quando a condição terminal no ciclo limite está perto do pico de x1 (& # x003d5 & # x0003c 0.05 e & # x003d5 & # x0003e 0.9). A Figura 8 mostra um exemplo de duas formas de onda de estímulo que resultaram em uma transição de quiescência para disparo repetitivo usando exatamente os mesmos pontos inicial e final. Conforme ilustrado aqui, essa multiplicidade ocorre porque a condição terminal está em um ciclo limite oscilatório. Neste exemplo, dois estímulos ótimos distintos resultam porque pode haver uma multiplicidade de oscilações subliminares antes que a trajetória salte para o ciclo de limite estável. Neste exemplo, o algoritmo convergiu para duas formas de onda de estímulo ideais que causaram uma ou três oscilações subliminares antes de induzir um salto para o ciclo limite.O estímulo que induz um salto após uma oscilação subliminar tem duração mais curta e amplitude maior, em comparação com o estímulo ideal que transita mais gradualmente. Essa descoberta correspondeu a um dos resultados que encontramos no modelo de Hodgkin-Huxley: se quisermos fazer a transição de estados mais rapidamente, é necessária mais energia.

O algoritmo de gradiente encontra várias soluções ótimas que induzem uma transição de quiescência para disparos repetitivos. A condição inicial (ponto fixo quiescente) e a condição terminal (& # x003d5 = 0) do algoritmo de gradiente são as mesmas para ambas as tentativas. A única diferença é o estímulo inicial gerado aleatoriamente que é dado ao algoritmo de gradiente. O painel superior mostra os dois estímulos ideais (verde e azul), enquanto o painel inferior mostra a resposta do x1 variável no modelo FitzHugh-Nagumo aos estímulos (verde e azul correspondentes). Para uma definição da fase, consulte os Métodos 2.4 e 2.5

Para fornecer uma comparação, calculamos novamente o pulso retangular ideal para mudar o modelo de FitzHugh-Nagumo de quiescência para disparo repetitivo e descobrimos que o melhor estímulo tinha uma amplitude de 0,11 e uma duração de 21,102, resultando em uma norma L 2 de 0,26. Os resultados do algoritmo de gradiente variaram de 0,0038 na melhor das hipóteses e 0,0097 na pior. Aqui podemos ver uma melhoria substancial em relação ao estímulo de pulso retangular padrão.

Como nossos resultados no modelo de Hodgkin-Huxley, a Figura 9 mostra os caminhos em direção à otimização resultante de diferentes sementes geradas aleatoriamente. Embora as diferentes sementes tenham uma ampla gama de normas L 2 e distâncias para a condição terminal, quase todas convergem para um dos dois clusters. Descobrimos que das soluções que convergiram, 78% das sementes geradas aleatoriamente convergiram para a solução maior da norma L 2 e 22% convergiram para a solução menor da norma L 2. Notamos na Figura 9 que há alguns pontos que derrubam a norma L 2 muito acentuadamente com iterações sucessivas, quase como se atraídos pelo ótimo local inferior, mas então voltam a subir, estabelecendo-se no ótimo local superior. Este fenômeno se deve ao fato de que uma única iteração pode fazer com que a norma L 2 caia muito, mas pode aumentar o erro nas condições terminais. Na próxima iteração, o algoritmo gradiente corrige o erro na condição terminal, levando ao retorno no valor da norma L 2.

Trajetórias de energia de estímulo por meio de 200 iterações de algoritmo de gradiente mostram convergência para duas formas de onda. Executamos 100 diferentes sementes geradas aleatoriamente com diferentes fatores de escala. As cores verde e azul correspondem às respectivas formas de onda verde e azul vistas na Figura 8

3.3 O Modelo FitzHugh-Nagumo: Disparo Repetitivo para Quiescência

A Figura 10 mostra as normas L 2 de 10 execuções de cada um dos 68 estudos computacionais de cada um dos pontos ao longo do ciclo de limite de disparo repetitivo até o ponto quiescente. Em todos os casos, a duração do estímulo foi fixada em 8 ms. A partir da figura, podemos ver que há novamente uma janela ótima onde & # x003d5 está entre 0,4 e 0,6 para começar a transição do ciclo de limite de disparo repetitivo para o ponto fixo quiescente. Isso corresponde ao local mais próximo do ponto fixo quiescente. Em torno desta janela, os melhores valores da norma L 2 ficam em torno de 0,012. Em comparação, a forma de onda de estímulo constante mais ótima tem uma norma L 2 de 0,064, por uma duração de 7,87 ms e uma amplitude de 0,09 dada em & # x003d5 = 0,72. É interessante notar que, ao usar uma forma de onda de estímulo constante, o estímulo discreto falha em suprimir disparos repetitivos em uma grande variedade de ângulos de fase. Em contraste, o método do gradiente permite a descoberta de novas formas de onda que induzem uma transição de uma ampla gama de fases em torno do ciclo limite para o estado quiescente.

O algoritmo de gradiente revela multiplicidade na transição de disparos repetitivos para quiescência em diferentes fases. O painel superior mostra a energia do estímulo dos diferentes estímulos otimizados (8 ms de duração) que induzem transições de diferentes fases de disparo repetitivo para o ponto fixo quiescente. Exemplos específicos de diferentes soluções são mostrados no painel inferior, rotulado uma Através dos d. Observe que b e c são soluções com as mesmas fases iniciais. Para uma definição da fase, consulte os Métodos 2.4 e 2.6

Havíamos discutido como a multiplicidade ocorria no exemplo anterior de estímulos que induzem o neurônio FitzHugh-Nagumo à transição do estado quiescente para o estado de disparo repetitivo, devido à natureza cíclica do sistema. Aqui, podemos ver que mesmo quando a duração do estímulo é menor que um ciclo, podemos encontrar multiplicidade com o sistema FitzHugh-Nagumo. Escolhemos o exemplo em que & # x003d5 = 0,678. A Figura 11 mostra a convergência de um conjunto de condições iniciais geradas aleatoriamente para os dois estímulos diferentes. A Figura 12 mostra dois estímulos que começam no mesmo ponto no ciclo de limite de disparo repetitivo e terminam no ponto quiescente.

Os estímulos aleatórios convergem para dois estímulos ótimos diferentes que induzem uma transição do disparo repetitivo para a quiescência. As trajetórias da energia do estímulo por meio de 100 iterações do algoritmo de gradiente são mostradas aqui. Todos esses estímulos começam na mesma fase inicial (& # x003d5 = 0,678). Para uma definição da fase, consulte os Métodos 2.4 e 2.6

O algoritmo de gradiente revela diferentes mecanismos de supressão de disparos repetitivos. O painel superior mostra os estímulos, o painel do meio o x1 resposta aos estímulos. O painel inferior mostra toda a resposta do espaço de estado a ambos os estímulos. Como podemos ver, um estímulo (verde) suprime o sistema rapidamente, enquanto o outro estímulo (tracejado em azul) faz com que o sistema funcione mais rapidamente em torno do estado oscilatório antes de entrar no estado quiescente. A linha preta na figura é marcada para mostrar o primeiro curso do sistema sem qualquer estimulação. Isso marca o ciclo limite do estado de disparo repetitivo. Para uma definição da fase, consulte os Métodos 2.4 e 2.6

Uma implicação interessante sobre esse resultado em particular é que há potencial para aplicar o algoritmo de gradiente à mudança de fase de um sistema oscilatório. O que podemos ver aqui na Figura 12 é que uma das soluções faz a transição quase imediatamente para o ponto quiescente, enquanto a outra faz um loop em torno do ciclo limite antes de entrar no ponto quiescente. Isso é interessante porque o período normal de um ciclo é de 12,84 ms. Capturamos uma instância em que o estímulo em 8 ms percorreu todo o ciclo limite. Poderíamos teoricamente configurar o algoritmo de gradiente para viajar de um ponto no ciclo limite para um ponto diferente do ciclo limite, exigindo que essa mudança de fase ocorra dentro de uma fração do período do ciclo limite. Desta forma, o algoritmo de gradiente pode ser usado para encontrar estímulos de deslocamento de fase ideais (Dean, Forger, & # x00026 Klerman 2009 Forger & # x00026 Paydarfar 2004 Serkh & # x00026 Forger 2014).


Bem vindo a CESE

Simulogic Inc. fornece produtos e serviços comerciais personalizados para a plataforma CESE, incluindo modelos prontos para uso para CESE, modelos desenvolvidos sob medida, extensões para o ambiente CESE e suporte e treinamento adicionais.

Suporte completo para MacOS X além de Windows e Linux, display dividido inovador para comparação de vários traços de dados simulados, importação / exportação de dados em MS Excel, arquivos de texto Axon e tabelas ASCII, cursores poderosos para selecionar regiões de interesse e realizar medições online e análise estatística, VirtuClamp para simular os protocolos de pinça de tensão e corrente e muito mais.

Analise os parâmetros do potencial de ação com o plugin de análise do potencial de ação para CESE Plus.

A pinça de potencial de ação está agora a um clique de distância com o plug-in de estímulo CESE Plus e HEKA.

O Cell Electrophysiology Simulation Environment (CESE) é uma estrutura abrangente projetada especificamente para realizar simulações eletrofisiológicas computacionais, por exemplo, simulações da atividade elétrica do miócito cardíaco. CESE é útil para simulações de potenciais de ação, correntes iônicas individuais e mudanças nas concentrações iônicas.

CESE é um programa de plataforma cruzada, ele roda em qualquer sistema que tenha o Java Runtime Environment (JRE) versão 1.5 (5.0) ou superior. Ele foi testado em Windows, Linux, Solaris, MacOS X e AIX.

Este documento descreve as principais características da CESE. Eles são apresentados para os dois domínios principais: Usuários e Desenvolvedores.

Usuários CESE

CESE é um ambiente integrado para a realização de simulações computacionais usando uma variedade de modelos eletrofisiológicos.

Nesta fase, a CESE permite a criação e execução dos modelos unicelulares (contendo as formulações da corrente Hodgkin-Huxley (HH) e da corrente Markoviana). Modelos de atividade elétrica de miócitos cardíacos com código-fonte estão incluídos na distribuição CESE. Esperamos estender o número de modelos disponíveis e adicionar alguns modelos neuronais no futuro.

A principal força da CESE está em sua uniformidade e uma interface de programa permanece a mesma para diferentes tipos de modelos. Você pode alternar facilmente entre os modelos e comparar os resultados da simulação. Os parâmetros do modelo podem ser modificados, selecionados para saída e / ou fixados da mesma forma padrão.

CESE estende o significado eletrofisiológico convencional do "grampo de tensão". Você pode fixar virtualmente qualquer variável do modelo, incluindo tensão (potencial de membrana), correntes iônicas totais ou individuais, concentrações iônicas, temperatura, variáveis ​​de bloqueio, etc. Os comandos de fixação podem ser funções complexas por peça, definidas individualmente para a variável do modelo de interesse . Isso abre possibilidades infinitas para a exploração do comportamento do modelo complexo.

CESE fornece visualizações de dados simples, mas eficientes. Os resultados da simulação podem ser apresentados na forma gráfica e tabulada. Os gráficos podem ser personalizados e as regiões de interesse ampliadas.

Mesmo que a CESE não tenha sido projetada para ser uma ferramenta de análise de dados, você pode gerar relações corrente-tensão (I-Vs) e calcular parâmetros estatísticos para um determinado sinal dentro do programa. Você pode exportar seus dados para os formatos ASCII, Axon Text File (ATF) e NetCDF para continuar a análise em seu pacote favorito.

Desenvolvedores CESE

A CESE foi criada desde o início para incorporar as melhores práticas de programação disponíveis para desenvolvedores Java, tanto em termos de consistência da interface do usuário quanto de clareza e reutilização de código. Sempre que possível, o CESE conta com APIs Java disponíveis (por exemplo, Java2D, JavaBeans, JAXP) para simplificar o código.

A criação do modelo requer uma série de funções de manutenção a serem codificadas & mdash, incluindo integradores ODE, rotinas para lidar com os parâmetros do modelo, salvar / restaurar o estado do modelo, visualizar os resultados da simulação, etc. A CESE fornece implementação para essas rotinas, portanto, você pode concentre-se em escrever o código para corrente (s) iônica de concreto, e a CESE cuidará do resto.

CESE não está tentando criar estruturas de programação complicadas por conta própria & mdash, em vez disso, utiliza APIs Java básicas. Por exemplo, os modelos são componentes Java em conformidade com a especificação JavaBeans. Usamos XML para especificar comandos de fixação e serialização de objetos Java para salvar / restaurar os parâmetros do modelo.

CESE é um projeto de código aberto criado com a reutilização de código em mente. Nós o encorajamos a estender a CESE e contribuir com o projeto. Envie um e-mail para [email protected] com comentários e perguntas. Visite nossa página SourceForge para obter acesso ao banco de dados de bugs, listas de discussão, repositório CVS, arquivos e fóruns.

Referências

  • Vascotto S, Missan S, La C
    CESE Plus: uma abordagem computacional para acelerar a pesquisa em eletrofisiologia.
    Laboratory Focus, 13 (4): 8-11, 2009.
  • Missan S, McDonald TF
    CESE: Ambiente de Simulação de Eletrofisiologia Celular.
    Appl Bioinformatics, 4 (2): 155-6, 2005.

Se você usou a CESE para sua pesquisa, por favor nos avise e nós adicionaremos seu artigo à lista.


Parâmetros para o modelo de potencial de ação de Fitzhugh-Nagumo - Biologia

Podemos usar o FPGA para fazer integração numérica rápida para resolver modelos de equações diferenciais de neurônios. Esta página descreve alguns modelos de neurônios e sua solução por meio de técnicas de DDA. O Digital Differential Analyzer (DDA) é um dispositivo para calcular diretamente a solução de equações diferenciais. Consulte a página DDA para obter mais detalhes sobre algoritmos e numéricos.

Modelo Neural de Izhikevich:
O neurônio
Eugene Izhikevich desenvolveu um modelo simples e semi-empírico de neurônios corticais. As propriedades de cada neurônio são controladas por 4 parâmetros, mais uma entrada de corrente constante. Seu site inclui programas matlab e descrições detalhadas do modelo de célula. Existem duas variáveis ​​de estado: v é o potencial de membrana e você é o recuperação de membrana variável. As equações diferenciais, configurações de parâmetros e exemplos de saída do modelo são mostrados abaixo em uma figura do site de Izhikevich. Uma versão eletrônica da figura e as permissões de reprodução estão disponíveis gratuitamente em www.izhikevich.com .

  • Modelo Soma (corpo celular e gerador de pico)
    • Constante de tempo de 16 etapas de tempo (intervalo de tempo nominal de 1/16 ms)
    • Dinâmica da lei quadrada (conforme explicado no História seção abaixo)
    • Parâmetros definíveis a, b, c, d e I (conforme explicado nos diagramas acima)
    • Atrasos de 1 a 64 etapas de tempo para simular o movimento de um pico ao longo de uma conexão de axônio.
    • Fonte de corrente constante controlada por um potencial de ação pré-sináptica.
    • Constante de tempo de queda exponencial configurável em 4, 8, 16, 32, 64, 128, 256 intervalos de tempo.
    • Peso sináptico ajustável entre -1 (inibitório) e +1 (excitatório).
    • Entrada para corrente em cadeia de outro módulo de sinapse ou de uma polarização de corrente constante.
    • Corrente determinada por uma condutância e a diferença de voltagem entre duas células.
    • Versões retificadoras (para frente ou para trás) e simétricas disponíveis.
    • Entrada para corrente em cadeia de outro módulo de sinapse ou de uma polarização de corrente constante.
    • O módulo altera o peso de uma conexão sináptica dependendo se o pico pós-sináptico segue (causal) ou conduz (não causal) o pico de entrada pré-sináptica. Conexões causais tornam-se mais fortes e as não causais mais fracas.
    • Constante de tempo de interação de pico exponencial configurável em 4, 8, 16, 32, 64, 128, 256 etapas de tempo
    • Os pesos delta para mudanças causais e não causais são configuráveis, assim como um peso inicial e um peso máximo.

    Dois pares recíprocos com sinapse elétrica entre eles.
    O módulo de nível superior (e arquivo do projeto) instancia 4 neurônios com as seguintes conexões:

    1. A saída do pico do neurônio 1 está conectada ao
      1. Neurônio 2 por meio de uma sinapse com peso -0,016
      2. Neurônio 3 através de uma sinapse elétrica com retificação e uma condutância de 1/64. O fluxo atual ocorre apenas se v1 & gtv2.
      3. LED 1
      1. Neurônio 1 por meio de uma sinapse com peso -0,016
      2. LED 2
      1. Neurônio 4 por meio de uma sinapse com peso -0,016
      2. Neurônio 1 através de uma sinapse elétrica com retificação e uma condutância de 1/64. O fluxo atual ocorre apenas se v1 & gtv2.
      3. LED 3
      1. Neurônio 3 por meio de uma sinapse com peso -0,016
      2. LED 4

      Par recíproco com sinapse modificada por STDP para um terceiro neurônio
      O módulo de nível superior (e arquivo do projeto) instancia 3 neurônios com as seguintes conexões:

      1. A saída do pico do neurônio 1 está conectada ao
        1. Neurônio 2 através de uma sinapse com peso -0,016
        2. Neurônio 3 por meio de uma sinapse com STDP e um peso inicial de zero. Os incrementos de peso STDP são ajustados para que o disparo do neurônio 3 eventualmente sincronize com o neurônio um.
        3. LED 1
        1. Neurônio 1 por meio de uma sinapse com peso -0,016
        2. LED 2

        As três imagens abaixo mostram as tensões iniciais não sincronizadas (neurônio 1 na parte inferior, neurônio 3 na parte superior), um estado intermediário e o estado final convegado gerado pelo módulo verilog acima. Inicialmente, ambos os neurônios estão espontaneamente ativos, mas com peso zero de conexão sináptica entre eles. A sinapse Hebbian é ajustada para que, após alguns milhares de picos, o acoplamento entre as células se desenvolva lentamente. No estado intermediário, você pode ver que todas as outras explosões do neurônio 1 estão disparando um pico no neurônio 3, mas também pode ver a pequena voltagem gerada por um acoplamento subliminar. O estado final de equilíbrio mostra um pico no neurônio 3 para cada explosão no neurônio 1.

        História (sequência de desenvolvimento de código)
        O módulo verilog de nível superior inclui o seguinte código para construir uma célula. Os parâmetros da célula e variáveis ​​de estado são redefinidos quando o botão KEY3 é pressionado. A variável de contagem é um prescaler de clock para desacelerar o cálculo por um fator de 4096 para que possa ser enviado através do codec de áudio. Observe que v e u são reduzidos por um fator de 100 para que a faixa dinâmica seja de -1 a +1, conforme exigido pelo sistema aritmético e pelo conversor D / A. Em comparação com as equações acima, as equações em escala têm todos os termos constantes divididos por 100. As constantes a e b não são escaladas (porque são multiplicadas pelas variáveis ​​em escala), mas c e d são divididos por 100. A constante (0,04 ) no termo v 2 acima é ampliado em 100 porque elevar ao quadrado (o escalado) v reduz o valor em 10.000. As equações finais em escala são:

        v1 (n + 1) = v1 (n) + dt * (4 * (v1 (n) ^ 2) + 5 * v1 (n) +1,40 - u1 (n) + I)
        u1 (n + 1) = u1 + dt * a * (b * v1 (n) - u1 (n))

        • dt = 1/16 (deslocar para a direita 4 vezes)
        • a no intervalo (0 a 0,1)
        • b no intervalo (0 a 0,3)
        • c no intervalo (-0,70 a -0,50)
        • d no intervalo (0,0005 a 0,08)

        A primeira passagem do código é mostrada abaixo.

        Duas imagens de osciloscópio mostram a tensão no traço superior e a variável de recuperação da membrana na parte inferior. Os traços mostrados usam parâmetros adequados para um neurônio em burst.

        Todo o projeto está compactado aqui.

        Dois neurônios com construção modularizada
        Todo o modelo neural, incluindo integradores, pode ser abstraído em um módulo de modo que o módulo de nível superior tenha apenas as instanciações de neurônios e as interconexões. As interconexões são feitas conectando a saída de pico de uma célula à entrada de corrente de outra. A corrente de polarização nesta versão pode ser definida a partir dos interruptores e o relógio e os sinais de reinicialização têm uma aparência um pouco mais limpa. Cada um dos neurônios pode ter seu trem de pico (binário) trazido para a porta paralela.

        O módulo neurônio encapsula todas as coisas numéricas. Observe que o Espigão a saída tem apenas uma amostra (etapa de integração) de largura. Cada neurônio usa 6 multipilers, então o número máximo de neurônios que podem rodar em paralelo é cerca de 11. Para mais isso, o sistema terá que ser serializado ou simplificado. O módulo produz a variável de tensão e uma saída de pulso. Os próximos 4 parâmetros controlam a dinâmica da célula.O próximo parâmetro é a entrada atual para a célula.

        A imagem do osciloscópio abaixo, que mostra as tensões dos dois neurônios, foi obtida com a corrente de entrada para cada um deles definida como 18'h0_2000.

        Otimizando o módulo neurônio e adicionando uma sinapse
        O modelo neural acima usa 6 de 70 multiplicadores de hardware para cada neurônio. Olhando para o intervalo de parâmetros fornecido na imagem do site de Izhikevich, sugere que os parâmetros aeb não precisam ser muito precisos para que o modelo funcione. Os dois multiplicadores por constantes foram substituídos por deslocamentos para aproximar as constantes em um fator de dois. Com essa mudança, o número de multiplicadores cai para dois / neurônio de forma que é possível colocar cerca de 35 neurônios modelo no FPGA. O módulo de nível superior foi otimizado ainda mais com a introdução de um modelo mais realista de interconexão de neurônios (veja a referência de Ruben Guerrero-Rivera abaixo) chamado de sinapse. O novo módulo de sinapse obtém uma saída de pulso de outro neurônio e usa o impulso (ou Espigão) como uma porta para carregar a condição inicial de uma unidade de decaimento exponencial. O rastreamento do osciloscópio abaixo mostra a variável de voltagem de um neurônio, junto com a saída da unidade de sinapse impulsionada pela mesma saída de pico do neurônio. Você pode ver que o peso negativo (18'sh3_8000) causa uma explosão de corrente negativa no momento de cada potencial de ação.

        No código do módulo de nível superior mostrado abaixo, a taxa de clock foi reduzida ao ponto em que picos individuais podem ser vistos em LEDs piscando. Alterar a corrente girando os interruptores causa uma variedade de comportamentos dinâmicos. Uma configuração de switch de cerca de 16'h0180 produz rajadas alternadas. Configurações mais altas causam disparos síncronos de alta taxa.

        O módulo sinapse é apenas outro integrador numérico que resolve a equação itertativa
        I (n + 1) = I (n) - I (n) / 16 + espiga1 * peso1 + espiga2 * peso2 + espiga3 * peso3
        A divisão por 16 resulta em uma constante de tempo razoável, mas pode ser transformada em variável. O módulo de sinapse produz uma variável de saída atual. Os próximos seis parâmetros são pares de de modo que cada unidade de sinapse tem um fan-in de três de outras células. O módulo de neurônio otimizado usa apenas dois multiplicadores de hardware, mas torna os parâmetros aeb configuráveis ​​apenas em uma potência de dois. Melhor resolução pode ser obtida usando outra entrada para definir um segundo valor de deslocamento. o Modelo FitzHugh-Nagumo com um oscilador semelhante a um neurônio
        O modelo FitzHugh-Naugumo é uma versão simplificada do modelo Hodgkin-Huxley (HH) de produção do potencial de ação do nervo. É também uma variação da equação de Van der Pol. Existem duas variáveis ​​de estado: v é o variável de ativação e w é o variável de inativação. Em termos de HH, v é alguma combinação de voltagem de membrana e ativação de sódio, enquanto w é uma combinação de inativação de sódio e ativação de potássio. Um código matlab para resolver o sistema é mostrado abaixo. Com a corrente constante mostrada, o sistema oscila. Do ponto de vista numérico, este sistema é interessante porque uma das variáveis ​​dinâmicas é ao cubo, o que requer alguns escalonamentos e multiplicadores. Converter isso para um ponto fixo de 18 bits exigiu um pouco de escala e mexeu com constantes. A versão implementada com dt = 2 -6 foi

        v (i + 1) = v (i) + dt * (v (i) - v (i) ^ 3 - w (i) + corrente)
        w (i + 1) = w (i) + dt * (1/16 * (v (i) - w (i)))

        Observe que as constantes foram arredondadas para quase uma potência de dois. O estouro durante a soma sugeriu dimensionar as variáveis ​​antes de adicioná-las. O código verilog no módulo de nível superior era

        sinal_multo v_2 (v2, v, (v & gt & gt & gt1)) // forma v-quadrado / 2
        sinal_mult v_3 (v3, v2, (v & gt & gt & gt1)) // forma v-cubed / 4
        atribuir vnew = v + (((v & gt & gt & gt2) -v3- (w & gt & gt & gt2) + (c & gt & gt & gt2)) & gt & gt & gt4) // dimensionar vars para corresponder a v3 (deslocamento total 6)
        atribuir wnew = w + (((v & gt & gt & gt1) - (w & gt & gt & gt1)) & gt & gt & gt9) // mult por um fator de 1/16 * dt (deslocamento total 10)

        A taxa de atualização foi dimensionada para colocar a forma de onda na região de frequência de áudio para que o codec de áudio pudesse ser usado para a saída. Um visor de osciloscópio é mostrado abaixo. O traço superior é v e o inferior é w.

        Todo o projeto está compactado aqui.

        FitzHugh_Naugumo com dois osciladores acoplados
        Conectar dois osciladores entre si pela conexão cruzada de suas variáveis ​​de ativação, mas com sinal negativo, faz com que os osciladores travem a fase 180 graus fora de fase. Há um estado metaestável de amplidude inferior se as condições iniciais para os dois osciladores estiverem exatamente equilibradas. Uma entrada de corrente para a primeira célula foi usada para quebrar a simetria. A parte principal do código do módulo de nível superior é mostrada abaixo. O acoplamento cruzado ocorre porque v1 / 4 é subtraído do feedback positivo para o oscilador 2 e v2 / 4 é subtraído do feedback positivo para o oscilador 1.
        O estado de bloqueio de fase de 180 graus é mostrado abaixo. Biologicamente, isso corresponde a sinapses inibitórias recíprocas entre os dois osciladores, forçando-os a alternar.

        O estado em fase metaestável é mostrado abaixo. Observe que a amplitude é menor porque cada oscilador inibe o outro. A próxima figura mostra a mudança da fase em fase para a fase 180 para uma diferença muito pequena nas condições iniciais. Observe que ambos os osciladores aumentam de fase no primeiro ciclo, mas no terceiro ciclo para dois estão bloqueados em 180 fases. O traço plano à esquerda do gatilho representa o estado de reinicialização.

        Todo o projeto está compactado aqui.

        • Ruben Guerrero-Rivera, et al, Kits de construção de lógica programável para modelagem neuronal em tempo hiper-real, Neural Computation, Volume 18, Issue 11 (novembro de 2006) Páginas: 2651 - 2679
        • Eugene M. Izhikevich, Simple Model of Spiking Neurons, IEEE Transactions on Neural Networks (2003) 14: 1569-1572


        Modelagem e simulação do potencial de ação em tecidos cardíacos humanos pelo método dos elementos finitos.

        As doenças cardíacas são uma das principais causas de morte no mundo e muito trabalho tem sido feito para elucidar as causas e os mecanismos dos problemas cardíacos. Os modelos de computador tornaram-se ferramentas valiosas para o estudo e compreensão dos fenômenos complexos da eletrofisiologia cardíaca. Os modelos têm desempenhado um papel importante neste campo [1] e apoiam os testes de novos fármacos, o desenvolvimento de novos dispositivos médicos e técnicas de diagnóstico não invasivas.

        A atividade elétrica é responsável pela contração periódica e ciclo de relaxamento do coração que impulsiona o sangue por todo o corpo. Portanto, a atividade elétrica é essencial para que o coração execute sua função. A maioria dos problemas cardíacos graves está, de fato, relacionada a distúrbios na atividade elétrica do coração [2].

        Modelos eletrofisiológicos do coração descrevem como a eletricidade flui pelo coração, controlando sua contração. Os modelos nos quais estamos interessados ​​consistem em sistemas de equações diferenciais. Os modelos da eletrofisiologia de uma célula são governados por sistemas de equações deferenciais ordinárias (ODEs), e os modelos da eletrofisiologia do tecido são governados por uma ou mais equações deferenciais parciais (PDEs). Normalmente, um modelo de PDE acoplado a um modelo de ODE para simular o tecido cardíaco que consiste em uma rede de células, os ODEs modelam a atividade elétrica nas células e os PDEs modelam a propagação da atividade elétrica pela rede como um todo.

        A atividade elétrica do coração como um todo é, portanto, caracterizada por uma estrutura multiescala complexa, que vai desde a atividade microscópica dos canais iônicos na membrana celular até as propriedades macroscópicas da propagação anisotrópica das frentes de excitação e recuperação em todo o coração. O modelo mais completo de tal configuração complexa é o modelo anisotrópico Bidomain [3], que consiste em um sistema de duas equações de difusão de reação parabólica degenerada que descrevem os potenciais intra e extracelulares no músculo cardíaco, juntamente com um sistema de equações deferenciais comuns que descrevem as correntes iônicas que fluem através da membrana celular, que estão associadas ao termo de reação não linear. Este modelo é computacionalmente muito caro devido ao envolvimento de diferentes escalas de espaço e tempo, e um modelo de tecido simplificado é o sistema de monodomínio anisotrópico, que consiste em uma equação de difusão de reação parabólica que descreve a propagação do potencial transmembrana acoplado a um modelo iônico, que tem tem sido amplamente utilizado na literatura [4, 5].

        O objetivo do presente trabalho é modelar e simular o potencial de ação no tecido cardíaco humano como ferramentas poderosas no estudo das arritmias cardíacas. E construir uma programação de computador para resolver os fenômenos de propagação de excitação usando o método dos elementos finitos. O resto do artigo está organizado da seguinte forma. Na seção 2, formulamos o modelo Bidomain. Na seção 3, formulamos o modelo de monodomínio. Na seção 4, descrevemos os tensores de condutividade. A seção 5 descreve as condições de contorno. A seção 6 descreve as correntes iônicas e os modelos de membrana FHN modificados. Finalmente, na seção 7 é apresentada a simulação numérica, variando o modelo de monodomínio do tecido cardíaco e o modelo iônico (FHN).

        O tecido cardíaco pode ser classificado em dois grupos: intracelular e extracelular. Para explicar os efeitos das diferenças de potencial através da membrana celular, o modelo Bidomain trata esses dois grupos como dois domínios separados. Cada ponto no coração é considerado em ambos os domínios, que podem ser considerados como sobrepostos um ao outro [6]. Cada ponto possui um potencial elétrico definido em cada domínio [2].

        O modelo Bidomain de tecido cardíaco é baseado no fluxo de corrente, distribuição de potencial elétrico e conservação de carga e corrente [3]. A descrição de cada domínio é baseada em uma versão generalizada da lei de Ohm que define a relação entre o campo elétrico E, derivado do potencial [phi] (V), da densidade de corrente J e do tensor de condutividade D.

        Considerando os espaços intracelulares e extracelulares especificamente, temos:

        [J.sub.i] = - [D.sub.i] [nabla] [[phi] .sub.i] [J.sub.e] = - [D.sub.e] [nabla] [[phi ] .sub.e] (2)

        Onde [J.sub.i] e [J.sub.e] são as densidades de corrente intra e extracelular, [D.sub.i] e [D.sub.e] são tensores de condutividade intra e extracelular, respectivamente, e [[phi] .sub.i] e [[phi] .sub.e] são o potencial elétrico nos espaços intracelulares e extracelulares.

        [nabla]. [J.sub.i] = [I.sub.m], [nabla]. [J.sub.i] = [I.sub.m] [nabla] ([J.sub.i] + [J.sub.e]) = 0 (3)

        Onde [I.sub.m] é a corrente transmembrana por unidade de volume [7], que é composta por um componente capacitivo e um componente iônico [I.sub.ion] resultante do fluxo de corrente através de canais de íons, bombas e trocadores no membrana celular.

        [I.sub.m] = [[beta] .sub.m] ([C.sub.m] [[derivada parcial] [V.sub.m] / [derivada parcial] t] + [I.sub. ion] + [I.sub.app]) (4)

        Onde [[beta] .sub.m] é a razão entre a área de superfície e o volume de uma célula cardíaca, [C.sub.m] a capacitância da membrana celular especificada, [I.sub.app] é a corrente de estímulo e [ V.sub.m] é a tensão transmembrana que é dada por:

        [V.sub.m] = [[phi] .sub.i] - [[phi] .sub.e] (5)

        Combinando as Equações (2) - (5), obtemos:

        [nabla]. [D.sub.i] ([nabla] [V.sub.m] + [nabla] [[phi] .sub.e]) = [[beta] .sub.m] ([C. sub.m] [[derivada parcial] [V.sub.m] / [derivação parcial] t] + [I.sub.ion] + [I.sub.app]) (6)

        [nabla]. (([D.sub.i] + [D.sub.e]) [nabla] [[phi] .sub.e]) = - [nabla]. ([D.sub.i] [ nabla] [V.sub.m]) (7)

        As equações parabólica (6) e elíptica (7) são as equações governantes do modelo de biodomínio do tecido cardíaco.

        O modelo monodomínio é uma simplificação do modelo Bidomain que é mais fácil de analisar e menos exigente computacionalmente [2]. Supondo que a anisotropia dos espaços intracelulares e extracelulares seja a mesma, isto é, que a condutividade no espaço extracelular é proporcional à condutividade intracelular.

        Onde, [lambda] é um escalar, representando a razão entre a condutividade dos espaços intracelulares e extracelulares. A escolha do valor de [lambda] pode determinar a precisão fisiológica, mas não é trivial selecionar um valor que forneça os melhores resultados [2]. O custo computacional de usar o modelo de monodomínio é cerca de metade a um décimo do custo de usar o modelo de Bidomain [9,5], dependendo da complexidade do modelo de célula usado.

        Substituir a equação (8) na equação (7) dá:

        [nabla]. (([D.sub.i] + [lambda] [D.sub.i]) [nabla] [[phi] .sub.e]) = - [nabla]. ([D.sub. i] [nabla] [V.sub.m]) (9)

        [nabla]. ([D.sub.i] [[phi] .sub.e]) = - [1/1 + [lambda]] [nabla]. ([D.sub.i] [nabla] [V .sub.m]) (10)

        Substituir a equação (10) na equação (6) dá:

        [nabla]. [1/1 + [lambda]] [D.sub.i] [nabla] [V.sub.m] = [[beta] .sub.m] ([C.sub.m] [[ derivada parcial] [V.sub.m] / [derivada parcial] t] + [I.sub.ion] + [I.sub.app]) (11)

        Se introduzirmos uma condutividade efetiva D = [[lambda] / (1 + [lambda]) [D.sub.i]], [1] obtemos o modelo de monodomínio do tecido cardíaco como:

        [nabla] .D [nabla] [V.sub.m] = [[beta] .sub.m] ([C.sub.m] [[derivada parcial] [V.sub.m] / [derivada parcial] t] + [I.sub.ion] + [I.sub.app]) (12)

        [[derivada parcial] [V.sub.m] / [derivada parcial] t] = -1 / [C.sub.m] ([I.sub.ion] + [I.sub.app]) + (1 /[C.sub.m][[beta].sub.m]) [nabla] .D [nabla] [V.sub.m] (13)

        O tensor de condutividade, d na equação (13) é definido em grande parte pela estrutura do coração. As células cardíacas são agrupadas em fibras musculares, e as fibras musculares são agrupadas em camadas de fibras [3]. Veja a Fig. 1. Essas estruturas influenciam o fluxo de condutividade elétrica é maior ao longo das fibras ao invés de através delas [2].

        Em cada ponto em cada domínio um tensor de condutividade local, denotado [D.sub.nfo], o tensor de condutividade para o caso em que não há orientação da fibra que pode ser representada por uma matriz simétrica 3 x 3 e definida na base formada como .

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (14)

        Onde [D.sub.x], [D.sub.y] e [D.sub.z] o valor escalar do tensor de condutividade na direção, x, y e z, respectivamente. O tensor de condutividade global D, que inclui a orientação da fibra, pode ser determinado usando uma transformação que envolve [D.sub.nfo] e uma matriz de rotação R

        Onde [R.sup.T] é a transposta da matriz de rotação R. Se cada folha é perpendicular a apenas uma das globais, consulte a Fig. 2 (por exemplo, o eixo x em um sistema de coordenadas (x, y, z)), então a base local para cada folha é fornecida por.

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (16)

        A fim de especificar as condições de contorno, deixe a superfície fechada [S.sub.H] ser um limite separando a região de bidomain H e o condutor de volume circundante B, e deixe a superfície fechada [S.sub.B] ligar a região B. n denota a unidade normal externa para [S.sub.H] e [S.sub.B]. [[phi] .sub.o] é o potencial cardíaco extra, [D.sub.i] é a condutividade escalar isotrópica fora de H en é a superfície externa normal do coração [9].

        [[phi] .sub.o] = [[phi] .sub.e], n. [D.sub.o] [nabla] [[phi] .sub.o] = n. ([D.sub. i] [nabla] [V.sub.m] + [D.sub.i] [nabla] [[phi] .sub.e]) em [S.sub.H] (17)

        n. [D.sub.o] [nabla] [[phi] .sub.o] = 0, em [S.sub.H] (18)

        Uma vez que as fontes em H estão relacionadas à presença de meio intracelular, que está ausente em B, podemos assumir que o vetor [D.sub.i] [nabla] [V.sub.m] na Eq. (17) é tangente à superfície SH, isso se torna as condições de contorno Numan (sem fluxo) [10] são:

        n. ([D.sub.i] [nabla] [V.sub.m]) = 0 em [S.sub.H] (19)

        n. [D.sub.e] [nabla] [[phi] .sub.e] = 0 em [S.sub.H] (20)

        6 Modelos de Membros e Correntes Iônicas

        O primeiro modelo de membrana para correntes iônicas que aparecem no monodomínio e no modelo Bidomain depende da escolha do modelo de membrana para a condutividade celular. O modelo mais antigo apareceu no trabalho sobre o potencial de ação do nervo por Hodgkin e Huxley [11]. Modelos deste tipo foram amplamente estudados e desenvolvidos por fisiologistas para o potencial de ação cardíaca sob a suposição de uma célula equipotencial, a variação no tempo do potencial de membrana [V.sub.m] para uma única célula é regida em tais modelos por uma equação diferencial ordinária [12].

        [[derivada parcial] [V.sub.m] / [derivada parcial] t] = - ([[I.sub.ion] + [I.sub.stim] / [C.sub.m]]) (21 )

        Onde [I.sub.ion] é a soma de todas as correntes iônicas transmembrana, [I.sub.stim] é a corrente de estímulo externa, e [C.sub.m] é a capacitância total da membrana. Tais modelos são naturalmente incluídos no quadro de trabalho acima, onde a corrente de estímulo é fornecida pelo termo de difusão nas configurações de Bidomain e monodomínio, modelando a propagação do estímulo através de células vizinhas.

        6.1 O modelo de célula Fitzhugh-Nagumo modificado

        FitzHugh e Nagumo reduziram a equação de Hodgkin-Huxley (quadridimensional) a um sistema bidimensional, extraindo a excitabilidade da dinâmica na equação original de Hodgkin-Huxley [13]. Doravante, o modelo iônico mais simples é o FitzHugh-Nagumo (FHN), que consiste em uma corrente iônica e uma variável de disparo. Supondo que o potencial V seja zero em repouso, a corrente iônica usa apenas uma variável de recuperação:

        [I.sub.ion] = cu (u - [alfa]) (u - 1) + w (22)

        [[derivada parcial] w / [derivada parcial] t] = [épsilon] (u - [gama] w) (23)

        Aqui w é uma variável de recuperação e u é o potencial transmembrana normalizado, definido por

        u = [[V.sub.m] - [V.sub.rest]] / [[V.sub.p] - [V.sub.rest]] (24)

        Onde [V.sub.m] é o potencial transmembrana, [V.sub.rest] é o potencial de repouso, e [V.sub.p] é o potencial de platô [14]. O potencial limite normalizado, a é definido de forma semelhante

        [alpha] = [[V.sub.th] - [V.sub.rest]] / [[V.sub.p] - [V.sub.rest]] (25)

        Onde, [V.sub.th] é o potencial limite.

        7 Método Numérico e Resultado da Simulação

        O método dos resíduos ponderados é uma técnica que pode ser usada para obter soluções aproximadas para equações diferenciais lineares e não lineares. Se usarmos este método, as equações de elementos finitos podem ser derivadas diretamente das equações diferenciais governantes do problema, sem qualquer necessidade de conhecer o funcional [15]. Derivamos as equações dos elementos finitos usando a abordagem de Galerkin.Um método de elemento finito de Galerkin foi desenvolvido para resolver o potencial de ação de um tecido cardíaco usando o modelo de monodomínio que se acoplou ao modelo modificado de FitzHugh-Nagumo (FHN) em um domínio geral com isotropia igual e sem orientação de fibra. Embora a implementação suporte problemas bidimensionais (2-D) e tridimensionais (3-D), para simplificar, apenas soluções 2-D são discutidas neste artigo.

        7.1 Consideração Numérica

        Ao resolver um sistema de EDOs ou PDEs, pode ser ineficiente usar um método numérico para cada parte do sistema [2]. Por exemplo, alguns componentes do sistema podem ser resolvidos de forma mais eficiente com um método numérico e outras partes do sistema resolvidas de forma mais eficiente com outro método numérico. Em vez de resolver tal sistema com um método numérico e aceitar as consequências da ineficiência, muitas vezes é melhor usar um método de divisão [6]. Um método de divisão usa uma estratégia de dividir e conquistar para resolver o sistema, dividindo o sistema em partes que podem ser resolvidas de forma eficiente com um método específico.

        7.2 Equação de Elementos Finitos

        Nesta seção, descrevemos um método de elemento finito (FEM) do modelo de monodomínio que se acoplou ao modelo de célula FitzHugh-Nagumo (FHN) modificado. A equação do modelo de monodomínio é

        [[derivada parcial] [V.sub.m] / [derivada parcial] t] + 1 / [C.sub.m] ([I.sub.ion] - [I.sub.app]) = ([1 /[C.sub.m][[beta].sub.m]]) [nabla] .D [nabla] [V.sub.m] (26)

        O RHS da equação (26) é chamado de parte de reação e o LHS da equação (26) é chamado de parte de difusão. No caso de um modelo de célula única, a parte de difusão é inexistente. Para resolver a equação (26), primeiro a expandimos.

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (27)

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (28)

        Estamos aplicando o método de Euler para derivada de tempo e a aplicação do método de Galerkin ao termo de difusão apenas sobre todo o domínio [OMEGA] para a equação (28).

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (29)

        Onde, W é a função de ponderação e a variável de campo de aproximação [[??]. Sub.m] é a seguinte.

        Onde, <[V.sub.m]>, o vetor do elemento nodal, i: índice do nó da grade e N a função de interpolação.

        Substituindo a equação (30) na equação (29), obtém-se:

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (31)

        Por integração por partes, o termo de difusão na equação (31).

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (32)

        Substituindo a equação (30) na equação (32), obtém-se:

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (33)

        Onde S são as condições de contorno de [OMEGA], Substituindo a equação (19) na equação (33) dá:

        [EXPRESSÃO MATEMÁTICA NÃO REPRODUÍVEL EM ASCII] (34)

        Substituindo as equações (34) e (30) na equação (32), aplicamos o método de Euler para a derivada de tempo. Podemos obter a equação algébrica resultante representada pelo problema de matriz usando uma notação de tensor.

        [([V.sub.m]). Sup.n + ​​1.sub.i] = [([V.sub.m]). Sup.n.sub.i] - ([[M] .sup. -1] [K] [([V.sub.m]). Sup.n.sub.i] + [S.sup.n.sub.i]) [DELTA] t (35)

        Onde sobrescrito n é o índice de tempo, M é o FEM, matriz de massa concentrada e K é a matriz de enrijecimento do FEM

        7.3 A Divisão do Operador

        Etapa 1 - Por uma etapa de tempo [DELTA] t / 2, resolvendo a equação de difusão

        [[derivada parcial] [V.sub.m] / [derivada parcial] t] = ([1 / [C.sub.m] [[beta] .sub.m]]) [[D.sub.x] ([[[derivada parcial] .sup.2] [V.sub.m]] / [[derivada parcial] [x.sup.2]]) + [D.sub.y] ([[[derivada parcial] .sup.2] [V.sub.m]] / [[derivada parcial] [y.sup.2]])] (36)

        [([V.sub.m]). Sup.n + ​​1.sub.i] = [([V.sub.m]). Sup.n + ​​1 / 2.sub.i] - ([[M ] .sup.-1] [K] [([V.sub.m]). sup.n + ​​1 / 2.sub.i]) [DELTA] t (37)

        Etapa 2- Para uma etapa de tempo [DELTA] para resolver a equação de reação

        [[derivada parcial] [V.sub.m] / [derivada parcial] t] = -1 / [C.sub.m] ([I.sub.ion] + [I.sub.app]) (38)

        [([V.sub.m]). Sup.n + ​​3 / 2.sub.i] = [([V.sub.m]). Sup.n + ​​1.sub.i] + [S.sup .n + 1.sub.i] [DELTA] t (39)

        Etapa 3- Para uma etapa de tempo [DELTA] t / 2 resolvendo a equação de difusão

        [[derivada parcial] [V.sub.m] / [derivada parcial] t] = ([1 / [C.sub.m] [[beta] .sub.m]]) [[D.sub.x] ([[[derivada parcial] .sup.2] [V.sub.m]] / [[derivada parcial] [x.sup.2]]) + [D.sub.y] ([[[derivada parcial] .sup.2] [V.sub.m]] / [[derivada parcial] [y.sup.2]])] (40)

        Aqui [V.sub.m] em RHS é o resultado da etapa 2.

        [([V.sub.m]). Sup.n + ​​2.sub.i] = [([V.sub.m]). Sup.n + ​​3 / 2.sub.i] - ([[M ] .sup.-1] [K] [([V.sub.m]). sup.n + ​​3 / 2.sub.i]) [DELTA] t. (41)

        Nesta seção, apresentamos os resultados de nossas simulações numéricas obtidas pelo código MATLAB. Consideramos domínios bidimensionais e rodamos modelos de células únicas e simulações de monodomínio de modelos de células acopladas, com o FHN modificado. nas simulações numéricas foi utilizado o parâmetro mostrado na tabela 1 para o modelo FitzHugh-Nagumo modificado.

        O tecido tem um potencial de repouso negativo, normalmente v = -90mV. A propagação do sinal no coração assume a forma de despolarização, onde o potencial aumenta rapidamente e atinge um potencial de pico positivo após alguns milissegundos. Segue-se uma fase de platô em que o potencial permanece positivo por pouco mais de 100 ms. Finalmente, a célula se repolariza, retornando ao seu potencial de repouso negativo. Todo o processo é denominado potencial de ação e normalmente dura cerca de 300ms. A Figura 1 mostra o potencial de ação de uma célula normal.

        Apresentamos aqui uma abordagem para simular a propagação da excitação nos tecidos cardíacos, baseada em modelos não lineares do tipo reação-difusão, considerando a abordagem de monodomínio. As correntes iônicas são expressas pelo modelo FHN modificado simples, especialmente projetado para tecidos humanos. Simulações numéricas em duas dimensões são fornecidas para mostrar o comportamento da propagação da excitação e da fase de repolarização.

        [1.] D. Noble. Modeling the heart insights, fracassos e progresso, Bioessays 24-1155-63, 2002.

        [2.] J. Sundnes, G.T. Lines, X.Cai, B. F. Nielsen, K.- A. Mardal e A. Tveito., Computing the Electrical Activity in the Heart, Springer-Verlag, Berlin, 2006.

        [3.] C. S. Henriquez, Simulando o comportamento elétrico do tecido cardíaco usando o modelo de Bidomain, Crit. Rev. Biomedical. Engineering, 21: 1-77, 1993.

        [4.] J. M. Rogers e A. D. McCulloch, A colocação-Galerkin modelo de elemento finito de propagação de potencial de ação cardíaca, IEEE Trans. Biomédico. Engenharia. 41: 743-57,1994.

        [5.] Yu Zhang, Ling Xia, Yinglan Gong, Ligang Chen, Guanghuan Hou e Min Tang, Solução Paralela na Simulação de Propagação Anisotrópica de Excitação Cardíaca, FIMH, LNCS 4466, 170-179, 2007.

        [6.] E. J. Vigmond, R. W. dosSantos, A. J. Prassl, M. Deo e G. Plank, Solvers for the Heart Bidomain equations, Progress in Biophysics Molecular Biology. 96 (1-3): 3-18, 2008.

        [7.] Barr, R. C. e Joseph D. Bronzino, Basic Electrophysiology, The Biomedical Engineering Handbook: Second Edition. Boca Raton: CRC Press LLC, 2000.

        [8.] J. Sundnes, BF Nielsen, K.- A. Mardal, X. Cai, GT Lines e A. Tveito, On the computational complex of the Bidomain and the monodomain models of electrophysiology, Annals Biomedical Engineering, 34 ( 7), 1088-1097, 2006.

        [9.] Yu Zhang, Ling Xia e Guanghuan Hou, Solução Eficiente de Equações de Bidomain na Simulação de Propagação Anisotrópica de Excitação Cardíaca, ICIC, LNBI 4115, 571-581, 2006.

        [10.] RH Clayton, O. Bernus, EM Cherry, H. Dierckx, FH Fenton, L. Mirabella, Modelos de eletrofisiologia do tecido cardíaco Progress, Challenges and open questions, Progress in Biophysics Molecular Biology, 104-22: 48, 2011 .

        [11.] Xile Wei, Jiang Wang e Bin Deng, Apresentando modelo interno para sincronização de saída robusta de neurônios FitzHugh-Nagumo em estimulação elétrica externa, Commun. Não linear. Sci. Numer. Simulat., 14: 3108-3119, 2009.

        [12.] K. H. W. J. tenTusscher, A.V.Panfilov, Alternans and spiral breakup in a human ventricular tissue model, American J. of Phys., 291: H1088-1100, 2006.

        [13.] Qingyun Wang, Qishao Lu, GuanRong Chen, Zhaosheng feng e LiXia Duan, Bifurcação e sincronização de modelos FHN acoplados sinapticamente com atraso de tempo, Caos, Solitons e Fractais, 39: 918-925, 2009.

        [14.] A. J. Pullan, M. L. Buist e L. K. Cheng, Mathematically Modeling the Electrical Activity of the Heart: From Cell to Body Surface and Back Again, World Scientific, New Jersey, 2005.

        [15.] Singiresu S. RAO, The Finite Element Method In Engineering, Fourth Edition, Elsevier Science & amp Technology Books, dezembro de 2004.

        [16.] Maria Murillo e Xiao-Chuan Cai, Um algoritmo paralelo totalmente implícito para simular a atividade elétrica não linear do coração, Numer. Lin. Algebra Appl., 11: 261-277, 2004.

        Recebido: 14 de agosto de 2011 / Aceito: 19 de agosto de 2011

        S. M. Shuaiby * M. A. Hassan * M. El-Melegy

        Departamento de Engenharia Mecânica, Universidade Assiut, Faculdade de Engenharia, Assiut, Egito.


        Assista o vídeo: The Action Potential (Dezembro 2021).