ESTIMATIVA DE PARÂMETROS E VARIÁVEIS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/25/pemufrj2012msc... · Figura 2.2 - Abertura e fechamento dos canais iônicos (GUYTON

ESTIMATIVA DE PARÂMETROS E VARIÁVEIS DE ESTADO DO MODELO DE

HODGKIN-HUXLEY

Diego Cardoso Estumano

Dissertação apresentada ao Programa de Pós-

graduação em Engenharia Mecânica, COPPE,

da Universidade Federal do Rio de Janeiro,

como parte dos requisitos necessários à

obtenção do título de Mestre em Engenharia

Mecânica.

Orientadores: Helcio Rangel Barreto Orlande

Marcelo José Colaço

Rio de Janeiro

Outubro de 2012


HODGKIN-HUXLEY


DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DO INSTITUTO ALBERTO

LUIZ COIMBRA DE PÓS-GRADUAÇÃO E PESQUISA DE ENGENHARIA (COPPE)

DA UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOS

REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE EM

CIÊNCIAS EM ENGENHARIA MECÂNICA

Examinada por:

________________________________________________

Prof. Helcio Rangel Barreto Orlande, Ph.D.

________________________________________________

Prof. Marcelo José Colaço, D.Sc.

________________________________________________

Prof. Emanuel Negrão Macêdo, D.Sc.

________________________________________________

Haroldo Fraga de Campos Velho, D.Sc..

RIO DE JANEIRO, RJ - BRASIL

OUTUBRO DE 2012

iii

Estumano, Diego Cardoso

Estimativa de Parâmetros e Variáveis de Estado do

Modelo de Hodgkin - Huxley/ Diego Cardoso Estumano. –

Rio de Janeiro: UFRJ/COPPE, 2012.

XV, 143 p.: il.; 29,7 cm.



Dissertação (mestrado) – UFRJ/ COPPE/ Programa

de Engenharia Mecânica, 2012.

Referencias Bibliográficas: p. 136-143.

1.Problemas Inversos. 2. Filtros de Partículas 3.

Monte Carlo via Cadeia de Markov. 4. Hodgkin-Huxley. I.

Orlande, Helcio Rangel Barreto et al. II. Universidade

Federal do Rio de Janeiro, COPPE, Programa de

Engenharia Mecânica. III. Titulo.

iv

“A mente que se abre a uma nova ideia

jamais voltará ao seu tamanho original”.

Albert Einsten

A minha família.

v

AGRADECIMENTOS

Aos meus pais, Maria de Fátima Estumano e Miguel Estumano, por terem se

dedicado tanto para sustentar a família e por não medirem esforços quando o assunto

era educação, agradeço por todas as vezes que foram trabalhar na mercearia mesmo

estando cansados para que não faltasse nada aos meus irmãos Renata Estumano e

Rodrigo Estumano e a mim.

Ao meu irmão por ser especial em todos os aspectos e por isso foi a motivação

deste trabalho, pois o início deste estudo foi a tentativa de entender como ocorrem as

crises convulsivas e epilepsia.

Aos meus orientadores Helcio Rangel Barreto Orlande e Marcelo José Colaço pela

paciência e orientação neste trabalho.

Ao Prof. Emanuel Negão Macêdo por ter me ajudado e incentivado a fazer o

mestrado na COPPE.

A minha namorada Anna Tsukui por ter sido tão companheira e por ter me ajudado

em todas as vezes que tive dificuldades com qualquer coisa.

Ao Léo Tsukui e Luciana Sampaio por terem promovido vários finais de semana

divertidos durante este período.

Aos amigos do LTTC Bernard Lamien, Antônio do Nascimento, Kleber Lisbôa,

Leonardo Bermeo, Wellington Betencurte, José Mir, Martim Costa, Gino Andrade, Jean

Monteiro, Rodrigo Basto, Marcus Curi, Thiago Corrêa, Gustavo Nóbrega, Kelvin

Chen,Masoud Ghanbari, Ali Allahyarzadeh, Ivana Fernandes e Ivana Cerqueira pela

amizade durante todo esse tempo de mestrado.

Aos funcionários Vera Noronha, Evanise Silva, Júlio Ramos, Paulo Silva e Paulo

Veiga por terem dado suporte para o desenvolvimento do trabalho.

vi

Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos

necessários para a obtenção do grau de Mestre em Ciências (M.Sc.).


HODGKIN-HUXLEY


Outubro/2012



Programa: Engenharia Mecânica

Este trabalho tem por objetivo a solução de problemas inversos relacionados à

formulação de Hodgkin-Huxley, para problemas de potencial de ação em células

excitáveis. As aplicações estudadas aqui incluem os axônios e as fibras de Purkinje.

Os problemas de estimativa de parâmetros e de estimativa de variáveis de estado

foram resolvidos neste trabalho. Para estimar os parâmetros do modelo, foram feitas

análises dos coeficientes de sensibilidade, revelando alta dependência linear entre

muitos parâmetros e baixas magnitudes, para axônios e fibras de Purkinje. Assim,

utilizou-se o método de Monte Carlo com Cadeia de Markov (MCMC), através do

algoritmo de Metropolis-Hastings, para estimar a capacitância da célula no problema

do axônio, usando medidas simuladas do potencial elétrico. Com a finalidade de

estimar as variáveis de estado, aplicou-se os seguintes algoritmos do filtro de

partículas: filtro de Amostragem por Importância e Reamostragem Sequencial

(Sampling Importance Resampling - SIR) e filtro Auxiliar de Amostragem por

Importância e Reamostragem Sequencial (Auxiliary Sampling Importance Resampling

- ASIR). Os problemas de estimativa de estado foram resolvidos com incertezas

Gaussianas e aditivas para o modelo de evolução de estado e para o modelo de

observação. Medidas simuladas do potencial elétrico foram usadas na solução do

problema de estimativa de estado.

vii

Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the

requirements for the degree of Master of Science (M.Sc.)

ESTIMATION OF PARAMETERS AND STATE VARIABLES OF THE HODGKIN-

HUXLEY MODEL


October/2012

Advisors: Helcio Rangel Barreto Orlande.


Department: Mechanical Engineering

This work aims at the solution of inverse problems related to the Hodgkin-Huxley

model for the action potential in excitable cells. The applications examined in this work

included axons and Purkinje fibers. The parameter estimation problem and the state

estimation problem were solved in this work. For the parameter estimation problem, an

analysis of the sensitivity coefficients reveals strong linear dependence among many

parameters and small magnitudes, for axons and Purkinje fibers. Hence, the Markov

Chain Monte Carlo (MCMC) method, through the Metropolis-Hastings algorithm, was

used to estimate the cell capacitance for an axon, by using simulated measurements of

the electric potential. For the estimation of the state variables, the following algorithms

of the particle filters were applied: the Sampling Importance Resampling (SIR) filter and

the Auxiliary Sampling Importance Resampling (ASIR) filter. The state estimation

problems were solved with Gaussian and additive uncertainties, for the state evolution

model and for the observation model. Simulated measurements of the electric potential

were used for the solution of the state estimation problem.

viii

SUMÁRIO

CAPÍTULO 1 - INTRODUÇÃO ................................................................................. 1

1.1 MOTIVAÇÃO E OBJETIVOS .......................................................................... 1

1.2 ORGANIZAÇÃO DO TRABALHO ................................................................... 3

CAPÍTULO 2 - REVISÃO BIBLIOGRÁFICA ............................................................ 4

2.1 CÉLULAS ....................................................................................................... 4

2.2 NEURÔNIOS .................................................................................................. 7

2.3 FIBRAS DE PURKINJE ................................................................................ 12

2.4 PROBLEMAS INVERSOS ............................................................................ 17

2.5 INFERÊNCIA BAYESIANA ........................................................................... 19

CAPÍTULO 3 - PROBLEMA DIRETO ..................................................................... 24

3.1 MODELO DE HODGKIN-HUXLEY ............................................................... 24

3.2 MODELO DE HODGKIN-HUXLEY APLICADO para AXÔNIO ...................... 29

3.3 MODELO DE HODGKIN-HUXLEY APLICADO para FIBRAS DE PURKINJE

31

CAPÍTULO 4 - PROBLEMAS INVERSOS ............................................................. 35

4.1 ESTIMATIVA DE PARÂMETROS ................................................................. 35

4.2 ESTIMATIVA DE VARIÁVEIS DE ESTADO ................................................. 43

CAPÍTULO 5 - RESULTADOS E DISCUSSÕES ................................................... 52

5.1 VERIFICAÇÃO DA SOLUÇÃO DO PROBLEMA DIRETO ............................ 52

5.2 ANÁLISE DOS COEFICIENTES DE SENSIBILIDADE ................................. 58

5.3 ESTIMATIVA DE PARÂMETROS ................................................................. 74

5.4 ESTIMATIVA DAS VARIÁVEIS DE ESTADO ............................................... 81

CAPÍTULO 6 - CONCLUSÕES E SUGESTÕES DE TRABALHOS FUTUROS ... 134

ix

LISTA DE FIGURAS

Figura 2.1 - Estrutura de uma célula (GUYTON e HALL, 2006). ................................... 5

Figura 2.2 - Abertura e fechamento dos canais iônicos (GUYTON e HALL, 2006). ....... 7

Figura 2.3 - Desenho esquemático de um neurônio (BIANCHI, 2003). ......................... 8

Figura 2.4 - Fases de um potencial de ação (Adaptado JARBAS et al., 2008). ............. 9

Figura 2.5 - Sistema de Purkinje ( GUYTON E HALL, 2006). ...................................... 13

Figura 2.6 - Fases de um potencial de ação para fibras do coração (Adaptado JARBAS

et al., 2008). ................................................................................................................ 14

Figura 3.1 - Circuito elétrico representando uma membrana celular (Adaptado,

SORIANO et al., 2006). .............................................................................................. 25

Figura 3.2 - Circuito elétrico representando uma membrana celular de fibras de

Purkinje (Adaptado, NOBLE, 1962). ........................................................................... 32

Figura 4.1 -. Esquema de previsão e atualização de dados (KAIPIO e SOMERSALO,

2004). ......................................................................................................................... 44

Figura 4.2 - Algoritmo para o filtro de partículas SIR (Adaptado, Ristic (2004)). .......... 49

Figura 4.3 - Algoritmo para o filtro de partículas ASIR (Adaptado, Colaço(2011)). ...... 50

Figura 5.1 - Verificação dos resultados obtidos para mV para o axônio. ..................... 53

Figura 5.2 - Verificação dos resultados obtidos para m para o axônio. ...................... 54

Figura 5.3 - Verificação dos resultados obtidos para n para o axônio. ....................... 54

Figura 5.4 - Verificação dos resultados obtidos para h para o axônio. ....................... 55

Figura 5.5 - Verificação dos resultados obtidos para mV para fibras de Purkinje. ........ 56

Figura 5.6 - Verificação dos resultados obtidos para m para fibras de Purkinje. ........ 56

Figura 5.7 - Verificação dos resultados obtidos para n para fibras de Purkinje. ......... 57

Figura 5.8 - Verificação dos resultados obtidos para h para fibras de Purkinje. ......... 57

Figura 5.9 - Análise da perturbação para o axônio, (a) domínio completo e (b)

domínio específico. ..................................................................................................... 59

Figura 5.10 - Análise da perturbação para fibras de Purkinje , (a) domínio completo

e (b) domínio específico. ............................................................................................. 59

Figura 5.11 - Análise de mCX , máx

NaGX , máx

kGX ,

LGX , NaVX ,

kVX e LVX para o axônio (a)

domínio completo e (b) região de interesse. ............................................................... 61

Figura 5.12 - Análise dos coeficientes de sensibilidade dos parâmetros de m e m

para o axônio (a) domínio completo e (b) região de interesse. .................................... 62

x

Figura 5.13 - Análise dos coeficientes de sensibilidade dos parâmetros de h e

h


Figura 5.14 - Análise dos coeficientes de sensibilidade dos parâmetros de n e

n


Figura 5.15 - Análise da dependência linear entre mCX ,

mdX , hcX e

ndX para o axônio

(a) domínio completo e (b) região de interesse. .......................................................... 66

Figura 5.16 - Análise de ,m PCX ,

,máxNaG P

X , NalGX ,

,Na PVX e ,k PVX para fibras de Purkinje (a)

domínio completo e (b) região de interesse. ............................................................... 68

Figura 5.17 - Análise dos coeficientes de sensibilidade dos parâmetros de 1kG e

2kG

para fibras de Purkinje (a) domínio completo e (b) região de interesse. ...................... 69

Figura 5.18 - Análise dos coeficientes de sensibilidade dos parâmetros de ,m P e ,m P


Figura 5.19 - Análise dos coeficientes de sensibilidade dos parâmetros de ,h P e ,h P


Figura 5.20 - Análise dos coeficientes de sensibilidade dos parâmetros de ,n P e ,n P

para fibras de Purkinje. ............................................................................................... 72

Figura 5.21 - Análise da dependência linear entre ,k PVX ,

kcX , ,m PjX ,

,h PbX e ,n PbX em

uma região específica.para fibras de Purkinje. ............................................................ 73

Figura 5.22 - Análise de w , estimando mC F para axônio. .................................. 76

Figura 5.23 - Análise de w , estimando kV mV para fibras de Purkinje. ................... 76

Figura 5.24 - Estimativa de mC F para o axônio. .................................................. 79

Figura 5.25 - Estimativa de kV mV para fibras de Purkinje. ...................................... 81

Figura 5.26 - Estimativa de h para axônio usando o filtro de partículas SIR com 500

partículas. ................................................................................................................... 83


partículas. ................................................................................................................... 83


partículas. ................................................................................................................... 84


partículas. ................................................................................................................... 84

Figura 5.30 - Estimativa de h para axônio usando o filtro de partículas ASIR com 50

partículas. ................................................................................................................... 85

xi


partículas. ................................................................................................................... 86


partículas. ................................................................................................................... 87


partículas. ................................................................................................................... 87

Figura 5.34 - Estimativa de Pn para fibras de Purkinje usando o filtro de partículas SIR

com 500 partículas...................................................................................................... 89


com 750 partículas...................................................................................................... 90


com 1000 partículas. ................................................................................................... 90


com 5000 partículas. ................................................................................................... 91

Figura 5.38 - Estimativa de Pn para fibras de Purkinje usando o filtro de partículas

ASIR com 50 partículas. ............................................................................................. 92


ASIR com 75 partículas. ............................................................................................. 92


ASIR com 100 partículas. ........................................................................................... 93


ASIR com 500 partículas. ........................................................................................... 93

Figura 5.42 - Estimativa de mV de axônio para o caso 9. ............................................ 96

Figura 5.43 - Estimativa de m de axônio para o caso 9. ............................................ 97

Figura 5.44 - Estimativa de h de axônio para o caso 9............................................... 97

Figura 5.45 - Estimativa de n de axônio usando para o caso 9. ................................. 98

Figura 5.46 - Estimativa de NaG de axônio para o caso 9. .......................................... 98

Figura 5.47 - Estimativa de kG de axônio usando para o caso 9. ............................... 99

Figura 5.48 - Estimativa de NaI de axônio para o caso 9. .......................................... 99

Figura 5.49 - Estimativa de kI de axônio para o caso 9. .......................................... 100

Figura 5.50 - Estimativa de mV de axônio para o caso 10. ........................................ 101

Figura 5.51 - Estimativa de m de axônio para o caso 10. ........................................ 101

xii

Figura 5.52 - Estimativa de h de axônio para o caso 10........................................... 102

Figura 5.53 - Estimativa de n de axônio para o caso 10........................................... 102

Figura 5.54 - Estimativa de NaG de axônio para o caso 10. ...................................... 103

Figura 5.55 - Estimativa de kG de axônio para o caso 10. ....................................... 103

Figura 5.56 - Estimativa de NaI de axônio para o caso 10........................................ 104

Figura 5.57 - Estimativa de kI de axônio para o caso 10. ........................................ 104

















Figura 5.74 - Estimativa de mV de fibras de Purkinje para o caso 13. ....................... 115

Figura 5.75 - Estimativa de m de fibras de Purkinje para o caso 13. ........................ 116

Figura 5.76 - Estimativa de n de fibras de Purkinje para o caso 13. ......................... 116

Figura 5.77 - Estimativa de h de fibras de Purkinje para o caso 13. ......................... 117

Figura 5.78 - Estimativa de NaG de fibras de Purkinje para o caso 13. ..................... 117

Figura 5.79 - Estimativa de kG de fibras de Purkinje para o caso 13........................ 118

Figura 5.80 - Estimativa de NaI de fibras de Purkinje para o caso 13. ...................... 118

Figura 5.81 - Estimativa de kI de fibras de Purkinje para o caso 13. ........................ 119

xiii



















Figura 5.100 - Estimativa de n de fibras de Purkinje para o caso 16. ....................... 130

Figura 5.101 - Estimativa de h de fibras de Purkinje para o caso 16. ....................... 130

Figura 5.102 - Estimativa de NaG de fibras de Purkinje para o caso 16. ................... 131

Figura 5.103 - Estimativa de kG de fibras de Purkinje para o caso 16. ..................... 131

Figura 5.104 - Estimativa de NaI de fibras de Purkinje para o caso 16. .................... 132

Figura 5.105 - Estimativa de kI de fibras de Purkinje para o caso 16. ...................... 132

xiv

LISTA DE TABELAS

Tabela 2.1 - Algumas áreas que pode-se aplicar problemas de inversos. ................... 17

Tabela 3.1 - Parâmetros do modelo de HODGKIN-HUXLEY aplicados no axônio de

uma lula. ..................................................................................................................... 30

Tabela 3.2 - Condições iniciais para o axônio. ............................................................ 30

Tabela 3.3 - Parâmetros do modelo de HODGKIN-HUXLEY aplicados para fibras de

Purkinje. ...................................................................................................................... 33

Tabela 3.4 - Condições iniciais para fibras de Purkinje. .............................................. 34

Tabela 5.1 - Dados para analisar w para o axônio..................................................... 75

Tabela 5.2 - Casos analisados para a estimativa de mC para o axônio. ...................... 78

Tabela 5.3 - Estimativa de mC F para o axônio...................................................... 80

Tabela 5.4 - Estimativa de kV mV para fibras de Purkinje. ....................................... 80

Tabela 5.5 - Comparação entre os filtros SIR e ASIR para o caso do axônio. ............. 88

Tabela 5.6- Comparação entre os filtros SIR e ASIR para o caso de fibras de Purkinje.

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

Tabela 5.7- Casos analisados para estimar variáveis de estado. ................................ 95

Tabela 5.8 - Valores de RMS para as variáveis de estado para os casos 9, 10, 11 e

12. ............................................................................................................................ 114

Tabela 5.9 - Valores de erro RMS para as variáveis de estado para os casos 13, 14,

15 e 16...................................................................................................................... 133

xv

Lista de Símbolos

P vetor com parâmetros

Y vetor de medidas

I mA corrente elétrica

mC F capacitância da célula

mV mV diferença entre o potencial elétrico interno e externo

m , n , h probabilidade de abertura dos canais iônicos

G mS condutância

, 1s taxas de transferência de íons

W matriz de covariância das medidas

V matriz de covariância para os parâmetros

J coeficiente de sensibilidade

X coeficiente de sensibilidade reduzido mV

desvio-padrão

erro com distribuição normal e média zero

N número de estados na Cadeia de Markov

Sobrescrito e Subscrito

p Purkinje

máx máximo

med medidas

esti estimados

aN sódio

K potássio

1

CAPÍTULO 1 - INTRODUÇÃO

Neste capítulo inicial, serão apresentados a motivação, os objetivos e a forma de

organização deste trabalho.

1.1 MOTIVAÇÃO E OBJETIVOS

Para estudar determinadas doenças que estão relacionadas com os sistemas

nervoso e cardíaco é necessário descrever como ocorrem os fluxos de cargas pelos

canais iônicos e como pulsos elétricos são transmitidos através das células. Diversos

fenômenos biológicos são estimulados através de trocas iônicas que ocorrem em

determinados tipos de células, as quais possuem a capacidade de gerar impulsos

elétricos, por exemplo, os neurônios, células endócrinas e as que constituem fibras

musculares. Esses impulsos são conduzidos em alta velocidade, fazendo com que o

potencial elétrico das células seja alterado e retorne ao seu estado de repouso em um

curto intervalo de tempo. Este fenômeno é denominado potencial de ação e os íons de

sódio, Na+, e potássio, K+, são os que mais contribuem para o mesmo (JARBAS et al.

,2008) .

A permeabilidade dos íons de sódio e potássio alteram a seletividade da

membrana celular para os remédios e agentes químicos. Portanto, descrever como

ocorre o potencial de ação nas células ajudará a entender como acontece a

propagação de um impulso nervoso no sistema neurológico. Tal entendimento pode

ajudar no tratamento de algumas síndromes no sistema nervoso como, o mal de

Alzheimer, a Epilepsia e o mal de Parkinson, assim como no sistema cardíaco, como

certas arritmias. Descrevendo-se o potencial de ação das células que constituem as

fibras de Purkinje, localizadas no coração, pode-se melhorar o funcionamento de

2

marca-passos e fazer prevenções de algumas doenças. Enfim, descrever como

ocorrem os impulsos elétricos e o potencial de ação ajudará a entender melhor

síndromes e doenças que estejam relacionadas com a deficiência na propagação de

impulsos elétricos e, além disso, pode ser útil também na criação de equipamentos

para corrigirem algum distúrbio no organismo.

Para descrever o potencial de ação, Hodgkin e Huxley, em 1952, elaboraram uma

formulação matemática que é constituída por um sistema de equações diferenciais

ordinárias (E.D.O) acopladas (HODGKIN e HUXLEY, 1952). Tal formulação serviu de

base para diversas pesquisas que visaram descrever melhor como ocorre o potencial

de ação nas células.

Este trabalho tem por objetivo a solução de problemas inversos relacionados à

formulação de HODGKIN-HUXLEY. Tanto o problema de estimativa de parâmetros,

quanto o problema de estimativa de variáveis de estado, são resolvidos aqui. Com a

finalidade de estimar as variáveis de estado, aplicou-se as seguintes metodologias

Bayesianas do filtro de partículas: filtro de Amostragem por Importância e

Reamostragem Sequencial (Sampling Importance Resampling - SIR) e filtro Auxiliar de

Amostragem por Importância e Reamostragem Sequencial (Auxiliary Sampling

Importance Resampling - ASIR). Para estimar os parâmetros do modelo, utilizou-se o

método de Monte Carlo via Cadeia de Markov (MCMC), com o algoritmo de

Metropolis-Hastings.

3

1.2 ORGANIZAÇÃO DO TRABALHO

O capítulo 2 apresenta conceitos biológicos como células, neurônios, axônios,

fibras de Purkinje e potencial de ação. Além destes conceitos, é apresentada uma

revisão bibliográfica de estudos realizados em problemas que utilizaram o modelo de

HODGKIN-HUXLEY em diversas finalidades, bem como de trabalhos que utilizaram

filtros de partículas para estimar variáveis de estado e o método MCMC para estimar

parâmetros de modelos matemáticos.

O capítulo 3 apresenta o modelo de Hodgkin-Huxley. Tal modelo é então

particularizado para aplicações em axônios e fibras de Purkinje. Os dados de entrada

e a metodologia utilizada para resolver os problemas propostos também são

apresentados.

As técnicas usadas para as soluções dos problemas inversos de interesse neste

trabalho são descritas no capítulo 4. Neste capítulo são apresentados o método de

Monte Carlo via cadeia de Markov para estimativa dos coeficientes do modelo de

HODGKIN-HUXLEY aplicados em axônios e fibras de Purkinje e os filtros de

Amostragem por Importância e Reamostragem Sequencial (Sampling Importance

Resampling - SIR) e o filtro Auxiliar de Amostragem por Importância e Reamostragem

Sequencial (Axiliary Sampling Importance Resampling - ASIR) para estimar as

variáveis de estado do modelo de Hodgkin-Huxley.

O capítulo 5 apresentada os resultados obtidos para os problemas direto e

inverso, com as formulações e métodos descritos nos capítulos 3 e 4.

O capítulo 6 apresentada as conclusões e propostas para trabalhos futuros, tanto

do ponto de vista da solução do problema direto, quanto da solução do problema

inverso.

4

CAPÍTULO 2 - REVISÃO BIBLIOGRÁFICA

Diversas áreas de estudo tentam descrever com detalhe como ocorre o potencial

de ação em células excitáveis, as quais possuem a capacidade de transmitir impulsos

elétricos. Nos ramos da engenharia, física, matemática e biofísica, pesquisadores vêm

desenvolvendo trabalhos para descrever matematicamente este fenômeno. Nos ramos

da medicina, biologia, eletrofisiologia, bioquímica e farmácia, pesquisadores estudam

quais unidades biológicas, realizam a transmissão de impulsos elétricos. E como

entender este fenômeno é muito importante para compreender como algumas

doenças e síndromes ocorrem e, desta forma, estudar como preveni-las, ou fazer um

tratamento adequado. Apresenta-se neste capítulo uma revisão bibliográfica dos

trabalhos ligados ao tema em questão, que serviram de base para este estudo.

2.1 CÉLULAS

A descrição de como acontecem diversos fenômenos biológicos, desde um ato

relativamente simples como mover uma parte do corpo, até a cura ou a prevenção de

algumas doenças, deve-se entender como acontecem os fenômenos na unidade

celular. A célula é a unidade viva básica do organismo e cada órgão é um agregado de

muitas células diferentes, mantidas juntas por estruturas de suporte intracelular

(GUYTON e HALL, 2006). As células são constituídas basicamente por três partes:

núcleo, citoplasma e membrana celular. A Figura 2.1 ilustra uma célula com suas

estruturas básicas. Pode-se observar que o núcleo está limitado pela membrana

nuclear, já o citoplasma encontra-se entre o núcleo e a membrana celular e por fim a

membrana celular limita a célula do meio exterior.

5

Figura 2.1 - Estrutura de uma célula (GUYTON e HALL, 2006).

CRUZ (2006) descreveu em detalhe algumas funções básicas das principais

estruturas celulares. O núcleo celular é responsável pelo controle da célula, visto que

este contém as informações sobre seu funcionamento. O citoplasma é um meio

aquoso onde se encontram pequenas unidades funcionais denominadas de organelas,

como, por exemplo, o retículo endoplasmático granular e agranular, o complexo de

Golgi, mitocôndrias e lisossomos. A membrana celular tem como principais funções:

formar uma barreira permeável, criar e explorar gradientes químicos e elétricos

transmembrana e promover a transdução de sinal.

O fato de a membrana celular limitar o espaço da célula em relação ao meio

externo, a mesma é responsável pela diferença de cargas e de concentrações de

espécies químicas. No caso da diferença de cargas, isso gera um gradiente elétrico e

a movimentação das cargas é no sentido de igualar as cargas entre os meios extra e

intracelular. Já o gradiente de concentração é responsável pelo movimento de

espécies no sentido de igualar as concentrações entre os meios extra e intracelular.

(CAMPOS, 2008). Quando não houver mais diferença entre a concentração de

espécies químicas e cargas entre os meios extra e intracelular, é dito que a célula está

em equilíbrio. Nesta condição, o potencial elétrico entre os meios extra e intracelular é

denominado de potencial de repouso ou potencial de membrana.

6

O potencial elétrico de alguns tipos de células podem sofrer variações e

posteriormente voltar ao seu potencial de equilíbrio. Diversos tipos de células como,

por exemplo, neurônios, fibras musculares e células endócrinas, possuem a

capacidade de gerar impulsos elétricos. Estes impulsos são denominados de

potenciais de ação, e desempenham um papel importante na regulação da função

celular, bem como constituem um mecanismo biológico para a transmitir informações

entre as células (JACKSON, 2006).

O transporte iônico que ocorre através da membrana celular é realizado por meio

de poros denominados de canais iônicos. A membrana possui uma camada

fosfolipídica que é permeável a determinados tipos de íons. Para que as cargas

iônicas atravessem pelos canais iônicos, os canais devem estar abertos. A

condutância dos canais iônicos é na ordem de grandeza de picosiemens (pS)( RAUCH

e GAMBLE, 1989).

A abertura e o fechamento dos canais iônicos podem ser controlados de dois

modos: por variação de voltagem e por controle químico. No controle por variação de

voltagem, a conformação do canal reage ao potencial elétrico através da membrana

celular, e a abertura ou fechamento é realizada de acordo com a diferença de cargas

entre os meios interno e externo (GUYTON e HALL, 2006). Esse controle é realizado

com a finalidade de manter um equilíbrio de cargas, como pode ser observado na

Figura 2.2. No controle químico, algumas comportas dos canais dependem da ligação

de substâncias químicas com proteínas, ou de suas ligações químicas nos canais, que

faz com que as comportas sejam abertas ou fechadas.

7

Figura 2.2 - Abertura e fechamento dos canais iônicos (GUYTON e HALL, 2006).

Será dada ênfase neste trabalho ao estudo do potencial de ação em neurônios e

fibras musculares, as quais são descritos a seguir.

2.2 NEURÔNIOS

O sistema nervoso é capaz de processar informações, através de unidades que

possuem a capacidade de captar os estímulos do ambiente e interpretá-los, podendo

armazenar a informação resultante ou simplesmente gerar uma resposta, na forma de

impulso nervoso. Tais unidades são chamadas de neurônios (SOUZA, 2005). Os

neurônios podem apresentar diversas formas e tamanhos, mas possuem

determinados aspectos em comum (MALTA et. al, 2004). Os neurônios são formados

por basicamente três partes: corpo celular ou soma, dendritos e axônio. O corpo

celular contém o núcleo e um grande número de fibras finas que se prolongam a partir

desse corpo. Cada neurônio possui apenas uma fibra longa, que é chamada de

axônio. Estas fibras são linhas de comunicação entre neurônios, ou entre neurônios e

tecidos musculares. A Figura 2.3 ilustra um neurônio com seus componentes básicos

(BIANCHI, 2003).

8

Figura 2.3 - Desenho esquemático de um neurônio (BIANCHI, 2003).

Os dendritos funcionam como uma "antena" para os neurônios, sendo cobertos

por milhares de sinapses (região onde o axônio entra em contato com outros

neurônios e transmite informações). A membrana dendrítica possui moléculas de

proteínas chamadas de receptores, que são especializadas na detecção dos

neurotransmissores. Os axônios funcionam como um "fio de telégrafo", o qual envia

informações à longa distância. Para conseguir desempenhar suas funções, o axônio

possui algumas estruturas especiais, como as células de Schawann que são

responsáveis por formarem camadas de membrana que fazem o isolamento elétrico.

Essa camada é denominada de mielina, cuja função é acelerar a propagação de

impulsos nervosos ao longo do axônio (BEAR et al. 2006).

Um potencial de ação ocorre quando há alterações rápidas do potencial de

membrana, seguida por uma restauração do potencial de repouso (JARBAS et al.

2008). A importância de se estudar o potencial de ação está em se analisar como é

realizada a troca de carga nas células. As fases que descrevem este fenômeno para

um neurônio são apresentadas na Figura 2.4 e podem ser descritas como:

9

Figura 2.4 - Fases de um potencial de ação (Adaptado JARBAS et al., 2008).

Fase 1 - Repouso: nesta fase o potencial de membrana encontra-se praticamente

estável, sendo seu valor correspondente à ausência de estímulo externo.

Fase 2 - Despolarização: o estímulo faz com que a membrana se torne permeável

às cargas positivas que estão no meio extracelular. O íon de sódio tem maior

facilidade de penetrar na membrana devido a apresentar um dos maiores gradientes

de concentração. Com isso, o potencial no interior da membrana se torna bem mais

elevado. Os canais de sódio são constituídos por 2 tipos de comportas. As primeiras

são abertas quando ocorre o estímulo que provoca um aumento muito rápido do

potencial, liberando a entrada de íons de sódio. Devido à velocidade alta de abertura

das comportas, estes canais são chamados de canais rápidos de sódio. Quanto aos

canais de potássio, mais lentos, a variação da membrana para um potencial maior

provoca uma lenta abertura dos mesmos, que nesta fase não permite o fluxo do íon.

Nesta fase ocorre um grande aumento do potencial elétrico.

Fase 3 - Repolarização: Atingindo o valor de potencial interno máximo, os canais que

permitiam a entrada de sódio começam a se fechar e canais de potássio começam a

se abrir. Desta forma, se inicia um fluxo de potássio para o meio extracelular, por

haver gradientes favoráveis para esse fenômeno, e o potencial elétrico é reduzido.

10

Nesta etapa, as comportas de inativação de íons de sódio são fechadas, impedindo a

difusão de íons de sódio para o meio extracelular, ocorrendo a repolarização da célula.

Dependendo de como ocorra essa repolarização, pode-se ou não passar pela fase

seguinte, de pós-potencial hiperpolarizante. No final do processo, o potencial elétrico

se torna igual ao potencial de repouso.

Fase 4 - Pós-potencial hiperpolarizante: essa fase ocorre quando a repolarização é

rápida e o potencial elétrico acaba se tornando inferior ao potencial de repouso. A

razão para o aparecimento desta fase é o não fechamento de todos os canais lentos

de potássio, o que implica em um decréscimo do potencial elétrico da célula.

Fase 5 - Ação da bomba de sódio e potássio: Nesta fase final, as bombas de sódio

e potássio atuam com a função de restaurar as concentrações inicial de sódio e

potássio nos meios intra e extracelular.

Diversas pesquisas já foram realizadas para descrever como ocorre o potencial de

ação em células excitáveis. Um dos primeiros modelos matemáticos para analisar o

potencial de ação foi apresentado por HODGKIN e HUXLEY (1952). Este trabalho deu

origem a uma série de estudos relacionados a modelos matemáticos que possam

influenciar na forma como ocorre o potencial de ação. HODGKIN e HUXLEY (1952),

analisando nervos de lulas, descreveram os fluxos iônicos como sendo um circuito

elétrico. Neste circuito estão modeladas a permeabilidade e a condutância de cada

íon, bem como a fração dos canais iônicos por onde os íons atravessam a membrana

da célula. Esta formulação matemática é descrita por um sistema de quatro equações

diferenciais ordinárias acopladas. Por meio de experimentos, HODGKIN e HUXLEY

(1952) observaram que os íons de sódio e potássio são os que mais influenciam no

potencial de ação. Além disso, foi observado também que as condutâncias iônicas

variam com o tempo, enquanto que os potenciais de repouso dos íons, a resistência

da membrana e a condutância de alguns íons, são constantes.

PICKARD (1966) modificou o sistema descrito por HODGKIN e HUXLEY (1952)

para tentar minimizar duas dificuldades existentes na época: custo computacional

11

elevado, mesmo para casos simples, e a própria solução das equações diferenciais do

modelo de Hodgkin e Huxley. Em seu modelo matemático PICKARD (1966),

considerou que o axônio tinha a forma de um cilindro infinito. Em seu método de

solução foi utilizada a transformada de Fourier e a diferença de potencial elétrico foi

separada em duas partes, uma que depende do tempo e outra independente do

tempo. Assim, foi obtida a solução com um custo computacional bem inferior em

relação à solução obtida por HODGKIN e HUXLEY (1952).

STEIN e PEARSON (1971) analisaram a influência de alguns parâmetros sobre o

potencial de ação no axônio de lulas. Verificou-se que a amplitude do potencial de

ação é fortemente influenciada pela temperatura, tendo sido realizadas simulações a

6.3 ºC, 20 ºC e 32 ºC. Verificou-se que com o aumento da temperatura, o pico e o

tempo do potencial de ação reduzem consideravelmente. Outra análise realizada foi a

a descrição do nervo como sendo um cilindro, formulado através de um modelo

bidimensional, nas direções longitudinal e radial. Outra análise interessante foi verificar

que, quando se considera um axônio muito grande, a amplitude do potencial de ação

varia com a raiz quadrada do diâmetro do axônio.

LINDLEY et al. (1987) utilizaram um modelo semelhante ao desenvolvido por

HODGKIN e HUXLEY, onde foram adicionadas informações sobre o fluxo dos íons de

cálcio, e analisaram também a influência da temperatura na atividade dos neurônios.

Foi verificado que, para analisar a influência da temperatura, o modelo proposto por

HODGKIN e HUXLEY não se adequava bem aos dados experimentais obtidos para

neurônios do hipocampo. O modelo proposto por estes autores contém informações

sobre a temperatura na faixa de 35ºC a 40 ºC, tendo representado bem os dados

experimentais obtidos.

FRÖHLIC e JEZERNIK (2005), analisando o modelo de Hodgkin e Huxley,

formularam uma espécie de controladores para regular a abertura e fechamento dos

canais iônicos dos neurônios, mais precisamente os canais de sódio e, dessa forma,

controlar o potencial de ação dessas células. Os canais foram controlados na forma de

12

um circuito fechado e sem variar o estímulo externo. O objetivo deste trabalho foi

estimular engenheiros na área de controle a desenvolver aparelhos de estimulação

para reabilitação de pacientes que sofrem com disfunções do sistema nervoso.

SANTOS (2009) estudou modelos de atividades não sinápticas para descrever o

potencial trans-membrânico dos neurônios, utilizando o modelo de HODGKIN-

HUXLEY. Uma grande contribuição deste autor foi considerar a contribuição das

correntes iônicas provenientes dos cotransportadores, que são proteínas responsáveis

pelo transporte de íons e moléculas através da membrana celular.

BUHRY et al. (2012) utilizaram o método de evolução diferenciada para a

minimização da norma de mínimos-quadrados, a fim de obter vários parâmetros do

modelo de HODGIKIN-HUXLEY simultaneamente. Segundo os autores, métodos

utilizados anteriormente não eram capazes de estimar os parâmetros

simultaneamente, pelo fato de usarem métodos de otimização que conduziam a

função objetivo para um mínimo local e não para um mínimo global. Foram obtidos

resultados bastante satisfatórios, pois os parâmetros de valor absoluto pequenos

apresentaram erros bem pequenos em relação aos parâmetros teóricos. Por outro

lado, parâmetros como as condutâncias máximas de sódio e potássio apresentaram

erros de estimativa maiores que 10 %.

2.3 FIBRAS DE PURKINJE

O mau funcionamento do coração pode levar uma pessoa à morte. Isso se deve

ao fato de o coração ser responsável pelo bombeamento do sangue para todo

organismo e o bombeamento depende de um padrão sequencial de excitação e

contração (FONSECA, 2009).

O coração possui um sistema condutor e rítmico especial denominado de sistema

de Purkinje, que serve para gerar e conduzir rapidamente para todo o coração

13

impulsos elétricos rítmicos que provocam a contração do mesmo. Este sistema é

ilustrado na Figura 2.5 (GUYTON e HALL, 2006).

Figura 2.5 - Sistema de Purkinje ( GUYTON E HALL, 2006).

Segundo AIRES (1999), os mecanismos que descrevem o potencial de ação em

fibras do coração são semelhantes aos que ocorrem no axônio, mas há algumas

diferenças, por exemplo, os potenciais de ação cardíacos são mais longos e a fase de

repolarização é bem mais lenta. O comportamento do potencial de ação está

diretamente relacionado com a função que a célula ou a fibra desempenham. No caso

do axônio, a função básica é a transmissão de impulsos elétricos, e o potencial de

ação ocorre em períodos curtos. Já para o miocárdio, a função básica do potencial de

ação é garantir a propagação rápida e coordenada, pois este fenômeno está

relacionado com os ciclos de contração. Frequências ventriculares muito altas são

indesejáveis, pois reduzem o tempo do enchimento ventricular durante a diástole,

reduzindo assim a eficiência da bomba cardíaca.

14

Figura 2.6 - Fases de um potencial de ação para fibras do coração (Adaptado JARBAS

et al., 2008).

As fases para o potencial de ação para as células cardíacas, são apresentadas na

Figura 2.6. Uma diferença fundamental na fase de despolarização entre as células

nervosas e cardíacas é a maneira de como ocorre o fluxo de íons nessa etapa. Os

íons de potássio apresentam dois canais, que permitem sua passagem através da

membrana das células. Um canal é função do potencial elétrico da membrana e reduz

seu valor à medida que a célula é despolarizada, enquanto o outro canal começa a

aumentar este valor à medida que a célula cardíaca está sendo despolarizada

(GUYTON e HALL, 2006). O potencial de ação em células cardíacas forma um platô

na fase de despolarização. Este fato pode ser observado na Figura 2.6, e é devido à

forma como ocorre o fluxo de íons de potássio nesta fase. O platô que aparece no

potencial de ação no músculo cardíaco (Fibras de Purkinje), acontece pelo fato de que

a permeabilidade da membrana celular diminuir para os íons de potássio,

15

aproximadamente de cinco vezes. A seguir são descritas as fase do potencial de ação

em fibras de Purkinje (GUYTON e HALL, 2006).

Fase 1 - Repouso: ocorre da mesma maneira que em células nervosas, descritas na

seção 2.2.

Fase 2 - Despolarização: ocorre da mesma forma que em células nervosas, descritas

na seção 2.2.

Fase 3 - Repolarização: nas células cardíacas, não ocorre a repolarização imediata

como nas células nervosas, sendo esta etapa dividida em uma etapa rápida e outra

lenta. A etapa rápida é devido à inativação dos íons de sódio, sendo o início da

repolarização. Na etapa lenta, o potencial elétrico permanece em um platô por algum

tempo sendo este fenômeno causado pela velocidade baixa dos canais de potássio.

Fase 4 - Bomba de sódio e potássio: ocorre da mesma forma que em células

nervosas, descritas na seção 2.2.

Diversas pesquisas foram realizadas para descrever como ocorre o potencial de

ação nas fibras de Purkinje. NOBLE (1962) apresentou um modelo matemático

baseado no modelo HODGKIN-HUXLEY para descrever o potencial de ação nas fibras

de Purkinje. O modelo descreve a fluxo de íons de sódio da mesma forma

apresentada por HODGKIN e HUXLEY (1952). Porém, a grande diferença está na

forma como foi apresentado o modelo para descrever o fluxo do íon de potássio. Para

este íon, foi apresentada a condutância em dois termos, para representar o platô que

aparece no potencial de ação das fibras de Purkinje (ver Figura 2.5).

REUTER (1967) analisou a influência dos íons de cálcio em algumas propriedades

das fibras de Purkinje. Para esta análise foram testadas três concentrações de cálcio:

0, 1.8 e 7.2 mM. Foi verificado que a influência dos íons de cálcio sobre a resistência

da membrana é praticamente insignificante. Por outro, lado verificou-se

experimentalmente que os íons de cálcio tem uma influência significativa no potencial

de ação das fibras de Purkinje.

16

BASSANI (1979) estudou o potencial de ação das fibras de Purkinje analisando o

modelo de NOBLE (1952). Além de descrever como ocorre o potencial de ação nas

fibras de Purkinje, neste trabalho foram analisados pontos de equilíbrio na solução do

modelo de HODGKIN-HUXLEY. Para fazer esta análise, utilizou-se um método que

consiste na separação das variáveis em grupos, de acordo com suas constantes de

tempo. Através deste estudo, concluiu-se que há três pontos de equilíbrio no modelo

de HODGKIN-HUXLEY: o primeiro, próximo ao potencial de equilíbrio dos íons de

sódio;o segundo, próximo ao potencial de equilíbrio dos íons de potássio, sendo estes

dois primeiros pontos nós estáveis; o terceiro ponto de equilíbrio é um ponto de sela e

está situado em uma região entre os dois pontos de equilíbrio estáveis e separa os

tipos de trajetória diferentes.

VICTORRI et al. (1985) apresentaram um método híbrido para a solução do

modelo de HODGKIN-HUXLEY, aplicado para descrever o potencial de ação nas

fibras de Purkinje. Os autores notaram que métodos como Runge-Kutta apresentaram

erros elevados, sendo necessários incrementos bem pequenos na discretização do

tempo, resultando em um alto custo computacional.

WILLMS et al. (1999) usaram métodos de otimização para determinar alguns

parâmetros do modelo de HODGKIN-HUXLEY, usando um método que trata a função

como sendo a soma de várias funções da mesma forma, apenas com coeficientes

diferentes. Em seu modelo de inversão, foram utilizadas informações sobre a cinética

de abertura e fechamento dos canais iônicos por meio de medidas do potencial,

através da técnica de fixação de voltagem. Uma função de Boltzmann foi usada para

representar algumas funções no estado estacionário. Para fazer as estimativas, os

parâmetros foram definidos como tendo uma distribuição Gaussiana e os erros das

medidas foram considerados como sendo independentes. Foi verificado que, para

estimar vários parâmetros simultaneamente, os resultados obtidos não foram

satisfatórios. Uma metodologia criada para situações semelhantes a esta é a

formulação baseada em redes neurais. Um exemplo de aplicação desta técnica pode

17

ser observado no trabalho de HÄRTER et al. (2008), onde foram utilizadas as técnicas

de redes neurais e filtro de Kalman em um problema caótico, para fazer previsões em

meteorologia.

2.4 PROBLEMAS INVERSOS

Em vários campos da física e engenharia, observações diretas de algumas

variáveis são difíceis. Nestes casos, pode-se medir outras variáveis que sejam mais

fáceis e utilizar modelagens matemáticas adequadas para extrair informações da

variável que é difícil de ser medida. Estes tipos de problema de medição indireta são

denominados de problemas inversos (BECK,1999; KAIPIO e SOMERSALO, 2002).

A aplicação de problemas inversos pode ser feita em diversas áreas como, por

exemplo, nas engenharias: aeroespacial, biomédica, química, civil, computação,

elétrica e mecânica. A Tabela 2.1 a seguir mostra algumas aplicações ( McCORMICK,

2001).

Tabela 2.1 - Algumas áreas que pode-se aplicar problemas de inversos.

Engenharia Exemplos de aplicação

Aeroespacial Otimizar formas aerodinâmicas

Biomédica Imagem de tomografias

Química Concentração de traços químicos

Civil Localização de vazamentos em tubulações

Computação Reconhecimento de voz

Mecânica Estimar propriedades termofísicas

18

O estudo de problemas inversos faz parte de um novo paradigma de pesquisa,

onde as simulações computacional e experimental não são realizadas isoladamente,

mas sim de forma iterativa, a fim de que o máximo de informações sobre o problema

físico em questão seja obtido com as duas análises (ORLANDE et al., 2011;

BECK,1999).

Métodos sistemáticos para resolver problemas inversos foram bastante

pesquisados nos últimos vinte anos e tornaram esta técnica uma boa ferramenta para

análise e projetos de engenharia (ORLANDE, 2010).

A técnica de problemas inversos tem como uma de suas características a iteração

entre a simulação computacional e os dados obtidos de forma experimental, onde as

simulações computacional e experimental não são realizadas isoladamente, a fim de

que o máximo de informações sobre o problema físico analisado seja obtido com as

duas análises (BECK e ARNOLD, 1977).

Problemas inversos são matematicamente classificados como mal-postos, pois

sua solução não satisfaz as condições de existência, unicidade e estabilidade. Esta

característica faz com que, em alguns casos, exista a necessidade de usar técnicas de

regularização como as de Alifanov e Tikhonov (ORLANDE, 2010). A utilização de

técnicas de problemas inversos é capaz de estimar parâmetros e funções. A aplicação

de problemas inversos no estudo do modelo de HODGKIN-HUXLEY precisa da

aquisição de medidas de alguma variável. No caso deste estudo, esta variável é o

potencial elétrico das células. De fato, há equipamentos que são capazes de medir

esta variável, como por exemplo, o eletrocardiograma e o eletroencefalograma.

Diversos métodos para estimar parâmetros ou funções já foram desenvolvidos,

como os métodos de Levenberg-Marquardt, Gradiente Conjugado e de Monte Carlo

com Cadeia de Markov. Além destes métodos, outra metologia que está sendo

bastante estudada são os filtros de partículas, como os filtros de Kalman, SIR e ASIR.

Neste trabalho foram utilizados apenas o método de Monte Carlo via Cadeia de

Markov e os filtros de partículas SIR e ASIR.

19

2.5 INFERÊNCIA BAYESIANA

A informação que se tem sobre uma quantidade de interesse é fundamental na

estatística. O verdadeiro valor de é desconhecido e a ideia é tentar reduzir este

desconhecimento. Além disso, a intensidade da incerteza a respeito de pode

assumir diferentes graus. Do ponto de vista Bayesiano, estes diferentes graus de

incerteza são representados através de modelos probabilísticos para (EHLERS,

2007).

Em diversas situações a análise de dados envolve a estimativa de parâmetros ou

funções desconhecidas. Na maioria das aplicações, um conhecimento a priori sobre o

fenômeno é modelado. Este conhecimento permite formular a inferência Bayesiana

através do teorema de Bayes onde, através de distribuições a priori para parâmetros

desconhecidos e uma função de verossimilhança relacionada com os dados

observados, é possível determinar uma distribuição a posteriori para o evento

estudado (DOUCET et al., 2001).

A inferência Bayesiana é baseada em técnicas de simulação que utiliza amostras

da probabilidade a posteriori ( ) P | Y para extrair informações a seu respeito de P .

Como uma amostra é sempre um substituto parcial da informação contida em uma

densidade, métodos baseados em simulação são inerentemente aproximados e

devem apenas ser utilizados quando for constada a impossibilidade de extração

analítica de informação da posteriori, como é o caso do presente estudo (NAVEIRA-

COTTA, 2009). Para o cálculo da posteriori foram utilizados os métodos de Monte

Carlo via Cadeia de Markov e os filtros de partículas SIR e ASIR. A ideia básica destes

métodos é fazer um passeio aleatório no espaço de ( ) P | Y que converge para uma

distribuição estacionária, que é a distribuição de interesse do problema (NAVEIRA-

COTTA, 2009).

20

Esta seção apresenta uma breve revisão da literatura sobre Inferência Bayesiana

para análise de problemas inversos. Esta é uma técnica que pode ser aplicada em

diversas áreas como pode ser visto a seguir.

DOUCET et al.(2001) apresentaram em seu livro texto análises detalhadas do

método de Monte Carlo e de alguns filtros de partículas, além de apresentarem

diversas aplicações para os métodos explorados.

ARULAMPALAM et al. (2002) apresentaramm os algorítimos dos filtros de

partículas SIS (Sequential Importance Sampling), SIR, ASIR e filtros de Kalman

extendidos e compararam os erros obtidos de cada filtro. Os erros obtidos foram de,

5.35 (ASIR), 5.54 (SIR),6.09 (SIS), 23.8 (filtro de Kalman Extendido). Portanto para o

modelo analisado o método que se obteve menor erro foi o ASIR.

RISTIC et al . (2004) descreveram em seu livro-texto as metodologias e aplicações

de alguns tipos de filtros de partículas em diversos casos.

KAIPIO e SOMERSALO (2004) mostraram em seu livro, análises em metodologias

de problemas inversos em que o foco foi análises estatísticas dos erros gerados pela

modelagem.

ORLANDE et al. (2007) usaram uma aproximação para a função Likelihood

através de funções de bases radiais. Em seu trabalho mostrou-se que ao utilizar esta

técinica o tempo computacional para estimar parâmetros de um modelo que

representa a difusão e advecção de um sal é reduzido para um pouco menos da

metade.

PARTHASARATHY E BALAJI (2008), utilizando o método de Monte Carlo via

Cadeia de Markov, estimaram parâmetros como condutividade térmica e coeficiente

de convecção em uma placa bidimensional. O método empregado para resolver o

problema direto foi o de diferenças finitas. Neste estudo foi usado três tipos de

distribuição a priori: normal, log normal e Uniforme.

ORLANDE et al. (2008) estimaram o campo de temperatura para dois problemas

de condução de calor em regime transiente: um linear, onde foi aplicado o filtro de

21

Kalman; e outro com condutividade térmica e capacidade volumétrica variando com a

temperatura, onde foi utilizado o filtro de partículas SIR. Neste trabalho mostrou-se

que apesar de o desvio padrão para o modelo e para as medidas serem considerados

altos, obtiveram-se excelentes resultados para estimar o campo de temperatura com a

utilização dos filtros de Kalman e SIR nos casos em que o problema era linear e não

linear respectivamente. Além de mostrar que o custo computacional pa usar o filtro de

Kalman ser bem inferior em relação ao filtro de partículas SIR.

ORLANDE et al. (2011a) mostraram em seu livro diversos métodos de problemas

inversos, como inferência Bayesiana com o método de Monte Carlo via mostradores

de Gibbs ou com o algoritmo de Metrololis-Hastings, métodos de regularização de

Tikhonov entre outras. Apresentam também várias áreas de aplicação das técnicas

apresentadas.

ORLANDE et al. (2011b) apresentam as seguintes técnicas de problemas inversos

Levenberg-Marquardt, Gradiente Conjugado com problema Adjunto e Monte Carlo via

Cadeia de Markov. Além da teoria são apresentadas aplicações em problemas de

transferência de calor.

NAVEIRA-COTTA et al. (2011) utilizaram a Técnica da Transformada Integral

Generalizada, para resolver o problema direto de transferência de calor unidimensional

aplicado em material heterogêneo, e o método de Monte Carlo com Cadeia de Markov

(MCMC), para estimar simultaneamente a condutividade térmica, capacidade térmica

volumétrica, coeficiente de transferência de calor por convecção e fluxo de calor. Os

resultados estimados encontrados foram bem próximos dos valores teóricos.

CLEMENTI et al. (2011) utilizaram inferência Bayesiana, com o método MCMC,

para estimar uma distribuição de diâmetros de partículas de polímeros. Foram

estimadas distribuições multimodais, que se aproximaram consideravelmente das

distribuições verificadas experimentalmente.

BATISTA (2011) aplicou o filtro de partículas SIR para estimar função densidade

de tamanho das partículas com um modelo matemático de balanço populacional. O

22

problema direto foi resolvido através de transformada de Laplace com inversão

numérica.

COLAÇO et al. (2011) utilizaram filtros de partículas ASIR para estimar fluxo de

calor aplicado em problemas de convecção natural. Foram analisadas as influências

da frequência da aquisição de medidas e da quantidade de partículas utilizadas.

Verificaram que com mais partículas a estimativa do fluxo de calor foi melhor e com o

aumento da frequência de aquisição de dados o intervalo de confiança ficou maior.

SILVA et al. (2011) utilizaram os filtros de partículas SIR e ASIR em um problema

de solidificação, considerando um problema unidimensional, com as seguintes

finalidades: comparar os resultados obtidos com os filtros utilizados e estimar a

localização da frente de solidificação e fluxo de calor. Neste estudo verificou-se que é

melhor usar o filtro ASIR, pois este precisou de menos partículas para reproduzir

resultados semelhantes ao filtro SIR.

BEI et al.(2012), usando um modelo não linear, compararam a estimativa de

variáveis de estado usando o filtro de partículas SIR e um outro algorítimo que usa

uma combinação de partículas, de forma a diminuir o empobrecimento das partículas.

Para todos os modelos testados, este algoritmo proposto apresentou resultados

melhores que os obtidos com o filtro SIR.

LI et al.(2012), analisando um modelo de equações não lineares, apresentaram

uma metodologia para diversificar mais as partículas de modo a melhorar as

estimativas de variáveis de estado. Foi aplicado um método determinístico para

diversificar as partículas, apresentando resultados melhores que o filtro SIR.

A contribuição do presente estudo é a estimativa das variáveis de estado ou

parâmetros do modelo de Hodgkin-Huxley, com os algoritmos SIR e ASIR de filtros de

partículas e MCMC, respectivamente. Espera-se assim, desenvolver uma ferramenta

computacional que auxilie tratamentos de doenças estimuladas por descontrole de

fluxos iônicos em células excitáveis. Por outro lado, o presente trabalho se limita a

23

casos envolvendo apenas uma unidade biológica, por exemplo, um axônio ou uma

fibra de Purkinje.

24

CAPÍTULO 3 - PROBLEMA DIRETO

Neste capítulo, será apresentado o modelo proposto por Hodgkin e Huxley em

1952, com aplicações em axônio e fibras de Purkinje.

3.1 MODELO DE HODGKIN-HUXLEY

HODGKIN e HUXLEY (1952), ao analisarem o funcionamento de nervos de lulas

gigantes utilizaram uma técnica denominada de voltage clamp ou fixação de voltagem.

Nesta técnica dois eletrodos finos são inseridos ao longo do axônio. Como a

resistência dos eletrodos é baixa, pode-se considerar que todos os pontos do axônio

possuem o mesmo potencial elétrico e desta forma o potencial dependeria apenas do

tempo e não do espaço. Isto leva a uma modelagem que considera o axônio como

uma célula pontual. Um dos eletrodos é inserido na célula e outro fica no lado externo

da mesma. Dessa forma, é medida a diferença entre os potenciais elétricos no interior

e exterior da célula. Um terceiro eletrodo é inserido para controlar a diferença de

potencial elétrico entre o interior e exterior da célula (CRUZ, 2006).

HODGKIN E HUXLEY (1952), utilizando a técnica de fixação de voltagem

perceberam que algumas propriedades, como as condutâncias de sódio e potássio

variavam com a mudança de potencial no axônio. Por outro lado, há propriedades que

não variam com o tempo, como os potenciais de equilíbrio, a capacitância do axônio e

a condutância de outros íons.

No modelo de HODGKIN-HUXLEY, o axônio é modelado como sendo um circuito

elétrico, onde há uma corrente denominada de capacitiva, a qual é decorrente da

capacitância do axônio, e outra que é correspondente a soma das correntes de todos

25

os íons. Os íons de sódio e potássio são os que tem maior influência no potencial de

ação e estes, segundo o modelo de HODGKIN-HUXLEY, fluem através de

condutâncias que estão em paralelo, como pode ser observado na Figura 3.1. Além

dos canais de sódio e potássio, foi considerado um outro canal fictício para os íons

que influenciam menos o potencial de ação, como por exemplo, os íons de cloro.

A representação do sinal das correntes depende de como ocorre o fluxo dos íons.

No modelo é considerado como sendo positivos os íons que entram e negativos para

os íons que saem da célula.

Figura 3.1 - Circuito elétrico representando uma membrana celular (Adaptado,

SORIANO et al., 2006).

A corrente total I é calculada como sendo a soma das correntes capacitivas e as

correntes iônicas, isto é,

míons m

dVI I C

dt (3.1)

onde mC é a capacitância da célula no modelo de Hodgkin-Huxley.

26

A corrente total dos íons é calculada como sendo a soma das correntes dos íons

de sódio, potássio e íons fictícios, ou seja.

íons Na k LI I I I (3.2)

As correntes iônicas e capacitivas que atravessam a célula precisam ser

modeladas através das condutâncias de cada canal. Determinar seus valores é

necessário modelar como as condutâncias de cada canal podem ser representadas.

Para isso foram realizados experimentos utilizando a técnica de fixação de voltagem,

com as condutâncias dos íons de sódio e potássio, dadas pelas equações.(3.3)-(3.4).

Observou-se que a condutância do canal de potássio é diretamente proporcional à

(4n ), onde n representa a fração que está aberta deste canal e a condutância do

sódio é proporcional ao produto (3m ) com ( h ), onde m e h representam

respectivamente as probabilidades de o canal de sódio estar aberto e fechado

(HODGKIN E HUXLEY, 1952).

Segundo CRUZ (2006), as variáveis de ativação m e n presentes nas eqs.(3.3)-

(3.4) podem ser interpretadas como sendo a probabilidade de o canal estar aberto

para os canais de sódio e potássio respectivamente e a variável de inativação h seria

a probabilidade de o canal estar fechado para o canal de sódio. As condutâncias

referentes aos íons de sódio e potássio são mostradas nas equações 3.3 e 3.4.

3máx

Na NaG G m h (3.3)

4máx

k kG G n (3.4)

Após determinar as condutâncias, para encontrar os valores das correntes,

multiplica-se a condutância pela diferença entre o potencial do axônio e o de

referência para cada íon, na seguinte forma:

27

Na Na m NaI G V V (3.5)

k k m kI G V V (3.6)

L L m LI G V V (3.7)

Nas equações (3.3) e (3.4), máx

NaG , máx

kG e LG representam as condutâncias

máximas para os íons de sódio, potássio e íon fictício respectivamente onde mV é a

diferença de potencial entre os meios interno e externo da célula e NaV , kV e LV

representam respectivamente os potenciais de equilíbrio para os íons de sódio,

potássio e íon fictício.

Substituindo a equação (3.3) na equação (3.5) e a equação (3.4) na equação

(3.6), tem-se que:

3máx

Na Na m NaI G m h V V (3.8)

4máx

k k m kI G n V V (3.9)

Substituindo as equações (3.7)-(3.9) na equação (3.2), obtem-se:

3 4máx máx

íons Na m Na k m k L m LI G m h V V G n V V G V V (3.10)

Finalmente, substituindo a eq.(3.10) na eq.(3.1), obtém-se a primeira equação

diferencial ordinária do modelo de HODGKIN-HUXLEY.

3 4máx máxmm Na m Na k m k L m L

dVI C G m h V V G n V V G V V

dt (3.11)

28

A corrente I dada pela equação (3.11) representa o estímulo externo. Outras três

equações diferenciais ordinárias são necessárias para descrever a cinética de

abertura dos canais. Tais equações foram propostas por HODGKIN e HUXLEY (1952)

na seguinte forma:

1m m

dmm m

dt (3.12)

1h h

dhh h

dt (3.13)

1n n

dnn n

dt (3.14)

onde os coeficientes α's e β's dependem de mV (HODGKIN e HUXLEY, 1952).

Os coeficientes α's representam as taxas dos íons no sentido de fora para dentro

da célula, enquanto os coeficientes β's representam o contrário. As equações (3.11-

3.14) representam o modelo de HODGKIN-HUXLEY e os valores dos coeficientes α's

e β's são determinados de acordo com o tipo de célula que é analisada.

Segundo MALTA et al. (2004) uma grande limitação do modelo de HODGKIN-

HUXLEY está relacionada com o fato do mesmo ser baseado na técnica de fixação de

voltagem, onde os experimentos não levam em consideração a posição em que há o

registro do potencial elétrico. Outra grande desvantagem é o fato de o modelo de

HODGKIN-HUXLEY não levar em consideração as diferentes densidades de canais

iônicos ao longo do axônio.

29

3.2 MODELO DE HODGKIN-HUXLEY APLICADO PARA AXÔNIO

HODGKIN E HUXLEY (1952), determinaram as equações para as taxas α's e β's

para o axônio de uma lula, analisando dados experimentais. As expressões sugeridas

encontradas são apresentadas nas equações (3.15-3.20).

0.1 25

exp 0.1 25 1

m

m

m

V

V

(3.15)

4exp18

mm

V

(3.16)

0.07exp20

mh

V

(3.17)

1

exp 0.1 30 1h

mV

(3.18)

100.01

exp 0.1 10 1

m

n

m

V

V

(3.19)

0.125exp80

mn

V

(3.20)

onde mV é dado em milivolt.

As equações para representar as taxas α's e β's foram determinadas a partir de

ajustes de curvas aos dados experimentais. Os valores dos outros parâmetros

requeridos pelo modelo (ver equações (3.11)-(3.14) são apresentados na Tabela (3.1)

(HODGKIN e HUXLEY, 1952)).

30

Tabela 3.1 - Parâmetros do modelo de HODGKIN-HUXLEY aplicados no axônio de uma lula.

Parâmetros Valores Parâmetros Valores

mC F 1 kV mV -12

máx

NaG S 120 máx

LG S 0.3

NaV mV 115 LV mV 10.6

máx

kG S 36

Para o caso do axônio foram utilizadas as condições iniciais utilizadas por

SORIANO et al. (2006). Estas podem ser verificadas na Tabela 3.2.

Tabela 3.2 - Condições iniciais para o axônio.

Parâmetros Axônio

mV mV -5

m 0

n 0.33

h 0.5

Após serem definidos os valores das taxas α's, β's, parâmetros e as condições

iniciais para resolver o modelo de HODGKIN-HUXLEY é necessário definir a corrente

de estímulo I . Neste trabalho este valor é de 6 mA para tempo entre 25 e 27 ms e 0

mA para os demais intervalos de tempo (SORIANO et al., 2006).

31

3.3 MODELO DE HODGKIN-HUXLEY APLICADO PARA FIBRAS

DE PURKINJE

NOBLE (1962) adaptou o modelo de HODGKIN-HUXLEY para fazer uma

modelagem matemática para ser utilizada nas fibras de Purkinje. Uma diferença entre

os potenciais de ação de axônios e fibras de Purkinje deve-se ao fato de que, nas

fibras de Purkinje há a presença de canais que possuem a cinética de abertura lenta;

estes são responsáveis pela presença do platô que se verifica no potencial de ação

que pode ser observado na Figura 2.6. Devido à presença do platô no potencial de

ação das fibras de Purkinje, a condutância de potássio foi dividida em duas partes,

uma para representar os canais rápidos e outro para representar os canais lentos. De

fato, segundo GUYTON e HALL (2006), a permeabilidade aos íons de potássio na

membrana celular miocárdica diminui aproximadamente cinco vezes, efeito que não

ocorre para os neurônios.

NOBLE (1962) representou a condutância dos canais lentos e rápidos de potássio

como sendo uma associação em paralelo com as condutâncias dos demais íons,

conforme pode ser observado na Figura (3.2).

Outra diferença em relação ao axônio está relacionada com a condutância dos

íons de sódio. Para as fibras de Purkinje, NOBLE (1962) observou que a condutância

de sódio tem um comportamento análogo ao descrito por HODGKIN e HUXLEY (1952)

mais um valor fixo. Em relação aos outros íons, além de sódio e potássio, foi

observado que não influenciavam muito no potencial de ação. Por isso, não foram

incluídos no modelo descrito. Por fim, a corrente de estímulo, I , é nula para este

caso, visto que estes tipos de células são auto excitáveis, isto é, não precisam de

estímulos para ocorrer a propagação do potencial de ação.

32

Figura 3.2 - Circuito elétrico representando uma membrana celular de fibras de

Purkinje (Adaptado, NOBLE, 1962).

As equações apresentadas por NOBLE (1962) correspondentes as taxas α's e β's,

para as correntes e condutâncias dos íons de sódio e potássio em fibras de Purkinje

são apresentadas nas equações (3.21)-(3.32).

3

, , 0.14máx

Na P Na P P pG G m h (3.21)

, ,

1

90 901.2exp 0.015exp

50 60

m P m P

k

V VG

(3.22)

4

2 1.2k PG n (3.23)

, 1 2k P k kG G G (3.24)

, , , ,Na P Na P m P Na PI G V V (3.25)

, , , ,k P k P m P k PI G V V (3.26)

33

,

,

,

480.1

48exp 1

15

m P

m P

m P

V

V

(3.27)

,

,

,

0.12 8

8exp 1

5

m P

m P

m P

V

V

(3.28)

,

,

900.17exp

20

m P

h P

V

(3.29)

1

,

,

42exp 1

10

m P

h P

V

(3.30)

4

,

,

,

10 50

50exp 1

10

m P

n P

m P

V

V

(3.31)

,

,

900.002exp

80

m P

n P

V

(3.32)

Além de definir as equações para as taxas α's e β's, assim como no caso do

axônio, é necessário definir os parâmetros para usar o modelo de HODGKIN-

HUXLEY. Estes parâmetros são apresentados na Tabela (3.3) enquanto a Tabela (3.4)

mostra os valores usados para condição inicial das variáveis de estado.

Tabela 3.3 - Parâmetros do modelo de HODGKIN-HUXLEY aplicados para fibras de

Purkinje.

Parâmetros Valores

,m PC F 12

,

máx

Na PG mS 400

,Na PV mV 40

,k PV mV -100

I mA 0

34

Para o caso das fibras de Purkinje foram utilizadas as condições iniciais utilizadas

por SORIANO et al. (2006), estas podem ser verificadas na Tabela 3.4.

Tabela 3.4 - Condições iniciais para fibras de Purkinje.

Parâmetros Fibras de Purkinje

,m PV mV 70

Pm 0.079

Pn 0.323

Ph 0.602

35

CAPÍTULO 4 - PROBLEMAS INVERSOS

Neste capítulo serão apresentadas as técnicas de problemas inversos que foram

utilizadas neste trabalho. Para estimar parâmetros, foi utilizado o método de Monte

Carlo com Cadeia de Markov (MCMC) e para estimar as variáveis de estado, foi

utilizado o filtro de partículas, nas formas de seus algoritmos SIR e ASIR.

4.1 ESTIMATIVA DE PARÂMETROS

Considerando o modelo de HODGKIN-HUXLEY para a formulação do problema

direto, será representado como P o vetor que contém os parâmetros da formulação e

Y o vetor que contém as medidas na seguinte forma:

1 2 3, , ..., NP P P PTP (4.1)

1 2 3, , ..., NY Y Y YTY (4.2)

onde N representa o número de parâmetros que aparecem na formulação do problema

Supondo que os erros de medidas sejam Gaussianos, com média zero, conhecida

a matriz de covariância W , assim o fato dos erros serem aditivos e independentes do

vetor de parâmetros P , a função de verossimilhança pode ser expressa como (BECK

e ARNOLD, 1977):

1

1221

| 2 | | exp2

IT

Y P

m mW Y V P W Y V P (4.3)

36

onde ,1 ,2 ,, ,....,m m m m I V P V P V P V P é obtido da solução do problema direto

dado pelo modelo de HODGKIN-HUXLEY nos tempos de medida it , 1,....,i I , com

um conjunto de parâmetros P .

A função de verossimilhança representa a densidade de probabilidade condicional

das medidas. Uma solução bastante comum de problemas inversos de estimativa de

parâmetros é maximizar a função de verossimilhança. Esta maximização pode ser

obtida por meio da minimização do expoente de função verossimilhança, equação

(4.3) (BECK e ARNOLD, 1977; ORLANDE, 2010), ou seja:

1T

S m mP Y V P W Y V P (4.4)

onde W é a matriz de covariância das medidas. A minimização da função objetivo,

equação (4.4), resulta em estimativas de máxima verossimilhança. Supondo que os

erros são, não-correlacionados, com distribuição normal, media zero e desvio-padrão

constante, a matriz 1

W torna-se uma matriz diagonal com elementos dados pelo

inverso das covariâncias dos erros (OZISIK e ORLANDE, 2000).

Supondo ainda válidas as hipóteses que os erros das variáveis medidas sejam

não correlacionados, com distribuição normal, média zero e desvio-padrão constante,

a norma de mínimos-quadrados torna-se uma função objetivo que fornece variância

mínima. Deste modo a norma dos mínimos quadrados pode ser escrita como (BECK e

ARNOLD, 1977; OZISIK e ORLANDE, 2000):

T

S m mP Y V P Y V P (4.5)

Na abordagem Bayesiana é possível levar-se em conta as informações

qualitativas a respeito dos parâmetros que aparecem na formulação do problema,

disponíveis antes das medidas e codificadas em termos de distribuições de

37

probabilidade, juntamente com as informações das medidas experimentais, para

aprimorar o processo de inferência estatística ou de tomada de decisão.

Assim, a solução de problemas inversos dentro de uma abordagem Bayesiana se

dá na análise sobre a densidade de probabilidade a posteriori, que é a distribuição de

probabilidade condicionais dos parâmetros dadas as medidas experimentais (KAIPIO

e SOMERSALO, 2004; BECK e ARNOLD, 1977).

Para a solução de um problema inverso dentro de uma abordagem Bayesiana,

são usadas as seguintes hipóteses básicas (KAIPIO e SOMERSALO, 2004):

1 - Todas as variáveis incluídas no modelo são modeladas como randômicas.

2 - A aleatoriedade descreve o grau de informação ou incerteza a respeito da

realização dos experimentos.

3 - O grau de informação a respeito das variáveis é codificado em distribuição de

probabilidade.

4 - O mecanismo para combinar as informações apontadas pelas medidas com

aquelas supostas para os parâmetros antes da realização das medidas é representado

pelo teorema de Bayes:

posteriori

||

P Y PP P Y

Y (4.6)

onde posteriori P é a densidade de probabilidade a posteriori, P é a densidade a

priori, | Y P é a função de verossimilhança e Y é a densidade de

probabilidade marginal das medidas, que representa, uma constante de normalização.

Supondo que os parâmetros P tenham uma distribuição a priori Gaussiana com

média μe matriz de covariância V a distribuição a priori para os parâmetros pode ser

escrita como (BECK e ARNOLD, 1977; KAIPIO e SOMERSALO, 2004):

38

1

1221

2 | | exp2

NT

P V P μ V P μ (4.7)

onde μ e V são respectivamente a média e a matriz de covariância para os

parâmetros.

Neste caso, podem ser obtidas estimativas pontuais para os parâmetros através

da maximização da distribuição a posteriori, que é obtida através da minimização da

seguinte função objetivo:

1 1T T

MAPS m mP Y V P W Y V P μ P V μ P (4.8)

O modelo de HODGKIN-HUXLEY, quando aplicado em axônio apresenta 26

parâmetros e quando aplicado em fibras de Purkinje apresenta 37 parâmetros.

Portanto para estimar parâmetros deste modelo deve-se usar um método que seja

capaz de estimar vários parâmetros e que seja capaz de calcular a distribuição a

posteriori quando esta apresentar integrais complexas para serem resolvidas.

A teoria dos processos estocásticos, onde a cadeia de Markov esta inserida é

geralmente definida como a parte dinâmica da teoria das probabilidades, onde se

estuda um conjunto de variáveis aleatórias. O método de Monte Carlo via Cadeia de

Markov é uma versão iterativa dos métodos de Monte Carlo tradicionais. A ideia é

obter uma amostra da distribuição a posteriori e calcular estimativas amostrais das

características desta distribuição usando técnicas de simulação iterativa, baseadas em

Cadeias de Markov (EHLERS, 2007, NAVEIRA-COTTA, 2009). Devido ao seu caráter

geral, que permite inclusive o tratamento de distribuições a posteriori não- analíticas,

utilizou-se neste trabalho o Método de Monte Carlo com Cadeia de Markov (MCMC).

39

Uma cadeia de Markov é um processo estocástico 0 1, ,..., nP P P , tal que a

distribuição de tP , dados todos os valores anteriores 0 1 1, ,..., tP P P depende apenas

de 1tP . Matematicamente temos (KAIPIO e SOMERSALO, 2004, EHLERS, 2007):

0 1 1 1| , ,..., |t t t tp p P A P P P P A P (4.9)

Uma cadeia de Markov é mais precisamente definida pela sua probabilidade de

transição ,p i j p i j , a qual define a probabilidade de que o processo,

estando no estado is mova-se para o estado js em um único passo, conforme segue:

1, |t j t ip i j p i j p s s P P (4.10)

O método de Monte Carlo via Cadeia de Markov requer, para que se obtenha uma

única distribuição de equilíbrio, que a cadeia de Markov seja (EHLERS, 2007):

homogênea, isto é, as probabilidades de transição de um estado pera outro

sejam invariantes;

irredutível, isto é, cada estado pode ser atingido a partir de qualquer outro em

um número finito de iterações;

aperiódica, isto é, não haja estados absorventes;

Assim uma condição suficiente para se obter uma única distribuição estacionária é

que o processo atenda à seguinte equação de balanço:

| |i i j jp i j p p j i p P Y P Y (4.11)

40

onde |i ip P Y e |j jp P Y são as probabilidades dos estados distintos da

distribuição de interesse.

Uma questão importante, de ordem prática, é como os valores iniciais influenciam

o comportamento da cadeia. A ideia é que conforme o número de iterações aumente,

a cadeia gradualmente converge para uma distribuição de equilíbrio. Assim, em

aplicações práticas é comum que os estados iniciais sejam descartados, como se

formassem uma amostra de aquecimento. O problema então consiste em construir

algoritmos que gerem cadeias de Markov cuja distribuição converge para a distribuição

de interesse (EHLERS, 2007).

Neste trabalho foi utilizado o método de Monte Carlo via Cadeia de Markov com o

algoritmo de Metropolis-Hastings para estimar os parâmetros do modelo de

HODGKIN-HUXLEY.

Este algoritmo é baseado no método de aceitação-rejeição, onde são gerados

valores candidatos *

P que pertencem a uma distribuição de proposta 1* |t

p

P P .

Neste trabalho foram utilizadas distribuições Gaussianas. O algoritmo de Metropolis-

Hastings é descrito a seguir (KAIPIO e SOMERSALO, 2004):

1. Inicializa-se o contador de iterações da cadeia 0i e arbitra-se um valor inicial

0P ;

2. Gera-se um candidato P* da distribuição 1* |t

p

P P ;

3. Calcula-se a probabilidade de aceitação 1 *|t

P P do valor candidato na forma:

1* *

1 *

1 1*

| || min 1,

| |

t

t

t t

p

p

P Y P PP P

P Y P P (4.12)

41

4. Gera-se um numero randômico u com distribuição uniforme, isto é ~ 0,1u U .

5. Se 1 *|t

u

P P , aceita-se o novo valor e faz-se 1 *i

P P . Caso contrário

faz-se 1i i

P P .

6. Incrementa-se o contador de i para 1i e volta-se ao passo 2.

Um parâmetro importante de ser observado neste método é a taxa de aceitação.

Uma consequência de ter uma baixa taxa de aceitação é deixar a cadeia de Markov

por várias iterações no mesmo estado. Por outro lado, para altas taxas de aceitação,

apesar de a cadeia de Markov precisar de menos iterações para evoluir, a

convergência ocorre de forma muito lenta. A taxa de aceitação pode ser controlada

pela escolha da distribuição da proposta 1* |t

p

P P . Para isso, pode-se fazer uma

análise da evolução da cadeia de Markov.

O coeficiente de sensibilidade ijJ , definido na equação (4.13), é uma medida da

sensibilidade do potencial estimado ,m iV em relação às variações do parâmetro jP .

Pequenos valores de magnitude de ijJ indicam que grandes variações em jP causam

pequenas mudanças em ,m iV . Nestes casos as estimativas dos parâmetros jP podem

ser difíceis, porque um mesmo valor de mV pode ser obtido para uma grande faixa de

valores de jP . Além da magnitude grande para o coeficiente de sensibilidade, é

necessário analisar quais parâmetros são linearmente dependentes, pois não é

possível estimar simultaneamente parâmetros que tenham os coeficientes de

sensibilidade linearmente dependentes (BECK et al., 1985).

,m i

ij

j

VJ

P

(4.13)

42

Em problemas que envolvem parâmetros com diferentes ordens de magnitude, os

coeficientes de sensibilidade com respeito aos vários parâmetros podem ter diferentes

ordens de grandeza, criando assim dificuldades na comparação e identificação da

dependência linear. Esta dificuldade pode ser aliviada através de uma análise dos

coeficientes de sensibilidade reduzidos, definidos como (BECK et al., 1985).:

,i j

iP j

j

VmX P

P

(4.14)

A aproximação por diferenças finitas pode ser empregada na determinação dos

coeficientes de sensibilidade, aproximando as derivadas de primeira ordem que

aparecem na própria definição dos coeficientes de sensibilidade, equação (4.13). Se

uma aproximação por diferença avançada for usada, tem-se os coeficientes de

sensibilidade aproximados segundo a equação (4.15). Se a aproximação de primeira

ordem não for suficientemente acurada, o coeficiente de sensibilidade pode ser

aproximado por diferença centrada na forma dada pela equação (4.16), esta foi a

formulação utilizada para o cálculo dos coeficientes de sensibilidade neste trabalho

(OZISIK e ORLANDE, 2000):

1 2 1 2, ,...., ,...., , ,...., ,....,m j j Npar m j Npar

ij

j

V P P P P P V P P P PJ

P

(4.15)

1 2 1 2, ,...., ,...., , ,...., ,....,

2

m j j Npar m j j Npar

ij

j

V P P P P P V P P P P PJ

P

(4.16)

43

4.2 ESTIMATIVA DE VARIÁVEIS DE ESTADO

Diversos problemas na ciência requerem a estimativa de variáveis de estado de

sistemas que são alteradas ao longo do tempo (ARULAMPALAM, 2002). Problemas

de estimativa de estado, também são designados problemas inversos não

estacionários. Em tais tipos de problemas, os dados disponíveis são usados junto com

o conhecimento prévio sobre os fenômenos físicos e os dispositivos de medição, a fim

de produzir estimativas sequenciais das variáveis dinâmicas desejadas (KAIPIO e

SOMERSALO, 2004). Isto é realizado de tal maneira que o erro seja minimizado

estatisticamente (MAYBECK, 1979).

Para fazer a estimativa das variáveis de estado, utilizou-se neste trabalho o

método do filtro de partículas através dos algoritmos Sampling Importance

Reasampling (SIR) e Auxiliary Sampling Importance Reasampling (ASIR) (DOUCET et

al., 2001). Estes métodos utilizam informações de alguma variável que possa ser

medida, junto com um modelo estocástico de evolução das variáveis de estado, para

fazer a estimativa das mesmas. A estimativa sequencial é feita de modo a diminuir

estatisticamente os erros, conforme ilustrado na Figura 4.1, onde . representa uma

distribuição estatística. A estimativa usando filtros de partículas segue a metodologia

apresentada por DOUCET et. al (2001).

44

Figura 4.1 -. Esquema de previsão e atualização de dados (KAIPIO e SOMERSALO,

2004).

Considerando um modelo de evolução para variáveis de estado x apresentado na

equação 4.17:

,kk k-1 k

x f x v (4.17)

onde f é em geral uma função não linear que depende de x e apresenta incertezas

representadas por v , o vetor x possui as variáveis de estado que devem ser

estimadas e por fim k representa o tempo.

As observações descritas pelo modelo dependem das variáveis de estado

contidas em x e dos valores das medidas z que estão relacionados por meio de uma

função h , conforme apresentado na equação 4.18 a seguir.

,k k k k

z h x n (4.18)

onde n representa a incerteza das medidas.

45

As equações (4.17)-(4.18) representam, respectivamente, os modelos de evolução

das variáveis de estado e das observações. Estes modelos são baseados nas

hipóteses a seguir:

(a) A sequência kx para .....,3,2,1k é um processo Markoviano, isto é:

| | k 0 1 k-1 k k-1

x x ,x , .......,x x x (4.19)

(b) A sequência kz para .....,3,2,1k é um processo Markoviano, isto é:

| | k 0 1 k k k

z x ,x , .......,x z x (4.20)

(c) A sequência para kx

depende apenas das medidas através das variáveis de

estado do tempo imediatamente anterior, isto é,

| | k k-1 1 k-1 k k-1

x x ,z , ......., z x x (4.21)

onde | a b denota a probabilidade condicional de ocorrer o evento a dado que b

ocorreu.

Para os ruídos das variáveis de estado e de observação, as seguintes suposições

são feitas (KAIPIO e SOMERSALO, 2004):

(a) Para i j , os ruídos de iv e jv , assim como in e j

n são mutuamente

independente e são também mutuamente independente do estado inicial 0x ;

(b) O ruído dos vetores iv e jn são mutuamente independente para todos os valores

de i e j .

46

Diferentes problemas podem ser considerados como modelos de evolução e

observação (RISTIC, 2004), conforme listado abaixo:

i) O problema de previsão, onde o objetivo é estimar |k 1:k-1

x z (4.22-a)

ii) O problema de filtragem, onde o objetivo é estimar |k 1:k

x z (4.22-b)

iii) O problema de suavização com retardo fixo, onde o objetivo é estimar

|k 1:k+p

x z , onde 1p é o retardo fixo (4.22-c)

iv) O problema de suavização com domínio inteiro, onde o objetivo é estimar

| Kk 1:

x z , onde , 1,...,K i K 1: i

z z é toda a sequência de medidas (4.22-d)

O problema estudado neste trabalho é um problema de filtragem. Assumindo que

| 0 0 0

x z x é uma informação conhecida (a priori) no instante t0, a distribuição

a posteriori |k 1:k

x z , é o objetivo do problema de filtragem. Esta distribuição pode

ser obtida por meio dos filtros Bayesianos através de duas etapas: previsão e

atualização como é mostrado na Figura 4.1, que ilustra o processo (KAIPIO e

SOMERSALO, 2004).

O filtro de Kalman é o mais conhecido e utilizado filtro Bayesiano. Porém, sua

aplicação é limitada a modelos lineares com ruídos Gaussianos e aditivos (DOUCET

et al., 2001). Extensões do filtro de Kalman foram desenvolvidas para casos menos

restritivos, usando técnicas de linearização. Do mesmo modo, os métodos de Monte

Carlo foram desenvolvidos de forma a representar a densidade a posteriori, em termos

de amostras aleatórias e pesos associados. Tais métodos de Monte Carlo,

normalmente designados por filtros de partículas, não requerem as hipóteses

restritivas do filtro de Kalman. Assim os filtros de partículas podem ser aplicados a

modelos não lineares e com erros não Gaussianos. Como o modelo de HODGKIN-

47

HUXLEY é não linear, foram usados os filtros de partículas para estimar as variáveis

de estado neste trabalho.

A ideia principal do filtro de partículas é a de representar a função de densidade a

posteriori requerida, através de um conjunto de amostras aleatórias com pesos

associados de modo a calcular as estimativas baseadas nestas amostras e pesos.

Tomando , 0,...,i Ii

0:kx com os respectivos pesos , 0,...,i Ii

kw e considerando

, 0,...,j j I 0:k

x x como o conjunto de todos estados até kt , onde I é o número

de partículas e os pesos são normalizados, isto é, 1

1I

i

k

i

w

. Então, a densidade a

posteriori em kt pode ser aproximada da seguinte forma (DOUCET, 2001; ORLANDE,

2008):

1

|I

i

k

i

w

i

0:k 1:k-1 0:k 0:kx z x x (4.23)

onde . é a função delta Dirac. Assumindo as hipóteses em (4.22-a,d) a equação

(4.23) pode escrita como:

1

|I

i

k

i

w

i

k 1:k-1 k kx z x x (4.24)

Um problema comum do método de filtro de partículas é o fênomeno da

degeneração, onde depois de alguns estados as partículas podem ter pesos

insignificantes. A degeneração implica em grande esforço computacional, porque a

densidade a posteriori torna-se praticamente nula. Este problema pode ser resolvido

aumentando o número de partículas ou, de forma mais eficiente, selecionando-se as

melhores partículas através de uma técnica de reamostragem.

48

A reamostragem geralmente envolve um mapeamento de medidas aleatórias

, wi i

k kx com medidas aleatórias 1, N i

kx com pesos uniformes. Isto pode ser

realizado se o número de partículas com pesos grandes ficar inferior a um

determinado número limitante. Uma alternativa é efetuar reamostragens

indiscriminadamente em todos os tempos kt , como é realizado nos filtros de SIR e

ASIR (DOUCET, 2001; ORLANDE, 2008).

Embora a etapa de reamostragem reduza os efeitos de degeneração, isso pode

levar a uma perda de diversidade e a amostra resultante pode conter uma grande

quantidade de partículas repetidas. Este problema é conhecido como empobrecimento

da amostra. Isso pode ser grave nos casos em que se tenha pequenos ruídos no

modelo de evolução de estado. Neste caso ocorre o colapso de todas as partículas em

uma única partícula em intervalos de tempo pequenos (DOUCET, 2001; ORLANDE,

2008).

Segundo RISTIC (2004), o método ASIR foi desenvolvido para melhorar o método

SIR. A ideia básica deste método é a reamostragem no tempo 1k (usando medidas

avaliadas no tempo k ), antes das partículas avançarem para o tempo k .

As variáveis de estado a serem estimadas neste trabalho são as apresentadas no

modelo de HODGKIN-HUXLEY: mV , m , n e h . Após a estimativa dessas variáveis é

possível estimar as propriedades NaG , kG , NaI e kI através das equações 3.3, 3.4,

3.8, 3.9, 3.21, 3.22, 3.24, 3.25 e 3.26, bem como suas incertezas utilizando-se:

2 2

2 2 2

Na

Na NaG m h

G G

m h

(4.25-a)

2

2 2

k

kG n

G

n

(4.25-b)

2 2

2 2 2

Na

Na NaI m h

I I

m h

(4.25-c)

49

2

2 2

K

KI n

I

n

(4.25-d)

onde i é a incerteza associada a variável i .

Os algoritmos para os filtros de partículas SIR e ASIR são apresentados nas

Figuras (4.3) e (4.4), respectivamente.

Algoritmo SIR do Filtro de Partículas 1. Início 1.1. Tome um conjunto de partículas de distribuição inicial )( 0x e obtenha

Niwx ii ,......1;, 00

1.2. Faça t=1; 2. Avaliação dos Pesos: Avance os estados no tempo t-1 para o tempo t usando o modelo de evolução de estado

2.1. Calcule os novos pesos: 1|i i

t t tw z x

2.2. Normalize os pesos:

N

i

i

t

i

ti

t

w

ww

1

ˆ

3. Reamostragem 3.1. Construa a soma dos pesos acumulativos, sendo computada por

i

tii wcc 1 para i = 2,...N ,com 01 c

3.1.2 Tome i = 1 e gere 1u de uma distribuição uniforme 1,0 NU

3.1.3 Para j = 1,,....,N faça

11

1 jNuu j .

Enquanto ij cu faça i = i + 1

Designar as partículas i

t

j

t xx

Designar os pesos para 1 Nw j

k

4. Cálculo da estimativa estado atual

1

ˆ|N

i i

t t t t

i

x z x w

Faça 1t t , se 1 finaltt , então pare.

Com as novas partículas, retorne ao segundo passo.

Figura 4.2 - Algoritmo para o filtro de partículas SIR (Adaptado, Ristic (2004)).

50

Algoritmo ASIR do Filtro de Partículas

Passo 1

Para 1,...,i N fazer partículas i

kx a partir da densidade a priori

1| i

k kx x e calcular alguma característica de kx , dado 1kx ,

como por exemplo a média 1|i i

k k kE x x . Usando a função

de verossimilhança calcule 1|i i i

k k k kw z w

Passo 2

Calcular o peso total i

kit w , normalizar o peso das

partículas e então calcular o peso para cada partícula, 1i i

k kiw t w

Passo 3

Reamostragem das partículas

Construir a soma dos pesos acumulativos pela relação,

1

i

i i kc c w , para 1,...,i N , com 0 0c

Para 1i começar com 1u com uma distribuição uniforme

10,U N

Para j = 1,...,N

Fazer 1

1 1ju u N j

Enquanto j iu c , fazer 1i i

Registrar amostras j i

k kx x

1j

kw N

.i j i

Passo 4

Para j = 1,...,N

Gerar partículas j

kx a partir da densidade de probabilidade a

priori 1| ij

k kx x e, usando a função de verossimilhança,

calcular os pesos correspondentes.

|

|

j

k kj

k ij

k k

z xw

z

Passo 5

Calcular o peso total i

kit w e normalizar o peso das

partículas e então calcular o peso para cada partícula, 1ij j

k kw t w

Figura 4.3 - Algoritmo para o filtro de partículas ASIR (Adaptado, Colaço(2011)).

Neste trabalho não foram realizados experimentos. Por isso, para as medidas de

,m medV , foi utilizado a seguinte metodologia:

51

,m med m medV V (4.26)

onde mV é obtido no modelo direto com parâmetros conhecidos,

med é o desvio-

padrão das medidas e é uma variável randômica com distribuição normal, média

zero e desvio-padrão unitário.

As estimativas realizadas serão analisadas com nível de 99% de confiança e, para

isso será utilizada a equação (4.27)

2.576esti esti

esti x x (4.27)

onde esti

x pode ser uma variável de estado ou um parâmetro e esti é o desvio-padrão

do que foi estimado.

.

52

CAPÍTULO 5 - RESULTADOS E DISCUSSÕES

Neste capítulo serão apresentados os resultados e discussões em relação a cada

caso analisado. Serão apresentados a verificação do problema direto, a estimativa dos

parâmetros com o método MCMC e a estimativa das variáveis de estado com os filtros

SIR e ASIR.

5.1 VERIFICAÇÃO DA SOLUÇÃO DO PROBLEMA DIRETO

A solução do problema foi obtida utilizando um computador com processador I5 de

2.3 GHz e 6 Gb de memória RAM. Os códigos foram implementados no software

Matlab® e para a solução do sistema de equações diferenciais ordinárias acopladas,

foi utilizada a subrotina ode15s.

Nos problemas relacionados ao potencial de ação nos axônios foi utilizada a

subrotina ode15s, que é baseada no método de Gear. Os valores das tolerâncias

absoluta e relativas para este caso foram de 610. Para o caso em que se analisou as

fibras de Purkinje foi utilizada também a subroutina ode15s, mas para este caso foi

necessário aumentar as tolerâncias absoluta e relativa para 510, pois com valores

menores apareceram mensagens de erro na simulação.

A verificação da solução do problema direto é uma etapa muito importante, já que

o mesmo será usado repetidas vezes nas soluções dos problemas inversos de

interesse neste estudo. Este trabalho teve como referência para a verificação do

problema direto os resultados obtidos por SORIANO et al.(2006), onde para fazer a

verificação da solução utilizou-se o software pega ponto para obter resultados obtidos

por SORIANO et al.(2006).

53

SORIANO et al.(2006) resolveram o modelo de HODGKIN-HUXLEY para os casos

do axônio e das fibras de Purkinje, com os dados de entrada apresentados nas

Tabelas (3.2) e (3.4). SORIANO et al.(2006) usaram a subrotina ode45 a qual utiliza o

método de Runge-Kutta na forma explícita para a solução do modelo de HODGKIN-

HUXLEY.

Uma comparação entre os resultados obtidos neste estudo e os de SORIANO et

al. (2006) são apresentados nas Figuras (5.1-5.4) para o caso do axônio.

0 10 20 30 40 50 60

0

20

40

60

80

100

Tempo (ms)

Vm

(m

V)

Presente trabalho

Soriano et al. (2006)

Figura 5.1 - Verificação dos resultados obtidos para mV para o axônio.

54

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tempo (ms)

Ativação m

Presente trabalho


Figura 5.2 - Verificação dos resultados obtidos para m para o axônio.

0 10 20 30 40 50 60

0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.7

0.75

Tempo (ms)

Ativação n

Presente trabalho


Figura 5.3 - Verificação dos resultados obtidos para n para o axônio.

55

0 10 20 30 40 50 60

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

Tempo (ms)

Inativação h

Presente trabalho


Figura 5.4 - Verificação dos resultados obtidos para h para o axônio.

Pode-se verificar nas Figuras (5.1-5.4) que as soluções obtidas para o axônio

estão em excelente concordância em relação aos resultados obtidos por SORIANO et

al. (2006).

Após ter verificado as soluções obtidas para o axônio, verifica-se o modelo de

HODGKIN-HUXLEY em fibras de Purkinje. Uma comparação dos resultados obtidos

por SORIANO et al.(2006) e pelo presente trabalho são apresentados nas Figuras

(5.5-5.8).

56

0 200 400 600 800 1000 1200 1400 1600 1800 2000

-80

-70

-60

-50

-40

-30

-20

-10

0

10

20

Tempo (ms)

Vm

(m

V)

Presente trabalho


Figura 5.5 - Verificação dos resultados obtidos para mV para fibras de Purkinje.

0 200 400 600 800 1000 1200 1400 1600 1800 2000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tempo (ms)

Ativação m

Presente trabalho


Figura 5.6 - Verificação dos resultados obtidos para m para fibras de Purkinje.

57

0 200 400 600 800 1000 1200 1400 1600 1800 2000

0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.7

Tempo (ms)

Ativação n

Presente trabalho


Figura 5.7 - Verificação dos resultados obtidos para n para fibras de Purkinje.

0 200 400 600 800 1000 1200 1400 1600 1800 20000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tempo (ms)

Inativação h

Presente trabalho


Figura 5.8 - Verificação dos resultados obtidos para h para fibras de Purkinje.

58

Pode-se verificar nas Figuras (5.5-5.8) que as soluções obtidas aqui para as fibras

de Purkinje também estão em excelente concordância com os resultados obtidos por

SORIANO et al.(2006).

5.2 ANÁLISE DOS COEFICIENTES DE SENSIBILIDADE

Antes de dar sequência à estimativa dos parâmetros com o método MCMC na

seção 5.3, procede-se nesta seção à análise dos coeficientes de sensibilidades com

relação a cada um dos parâmetros do modelo de HODGKIN-HUXLEY, para os casos

do axônio e das fibras de Purkinje.

Os coeficientes de sensibilidade foram calculados aqui por diferenças finitas,

usando a equação (4.16). Assim deve ser avaliado o valor da perturbação que deve

ser utilizado neste cálculo. De fato, ao escolher um valor muito pequeno para , os

valores do numerador e do denominador da equação (4.16) se tornam muito

pequenos. Desta forma, quando é efetuado o cálculo para o coeficiente de

sensibilidade há o aparecimento de erros numéricos. Por outro lado, o uso de grandes

valores resulta em erros de aproximação da fórmula de diferenças-finitas para as

derivadas contínuas dos coeficientes de sensibilidade.

As Figuras (5.9-5.10) mostram os resultados obtidos para os coeficientes de

sensibilidade reduzidos para os parâmetros mC no caso do axônio e kV no caso das

fibras de Purkinje, respectivamente. Nos dois casos foram analisados os seguintes

valores de : 310,

410e

510.

59

0 10 20 30 40 50 60-700

-600

-500

-400

-300

-200

-100

0

100

200

XC

m

(m

V)

Tempo (ms)

Coeficiente de sensibilidade reduzido para Cm

Vm

=10-3

=10-4

=10-5

(a)

26 27 28 29 30 31 32

-20

0

20

40

60

80

100

120

140

160

XC

m

(m

V)

Tempo (ms)

Coeficiente de sensibilidade reduzido para Cm

Vm

=10-3

=10-4

=10-5

(b)

Figura 5.9 - Análise da perturbação para o axônio, (a) domínio completo e (b)

domínio específico.

0 50 100 150 200 250 300 350 400 450 500-9000

-8000

-7000

-6000

-5000

-4000

-3000

-2000

-1000

0

1000

X V

k (m

V)

Tempo (ms)

Coeficiente de sensibilidade reduzido para Vk

Vm

=10-3

=10-4

=10-5

(a)

300 310 320 330 340 350 360 370 380 390 400

-800

-700

-600

-500

-400

-300

-200

-100

0

100

X V

k (m

V)

Tempo (ms)

Coeficiente de sensibilidade reduzido para Vk

Vm

=10-3

=10-4

=10-5

(b)

Figura 5.10 - Análise da perturbação para fibras de Purkinje , (a) domínio completo

e (b) domínio específico.

Ao analisar as Figura 5.9 e 5.10 pode-se observar que para valores de iguais a

410, os coeficientes de sensibilidade reduzidos das variáveis analisadas convergem,

em termos do incremento usado no cálculo de diferenças finitas. Embora não sejam

apresentados nestas figuras, valores de incremento menores que 10 -5 resultaram em

erros numéricos. Portanto, o valor para a perturbação usada para calcular as

derivadas através de diferenças finitas para o axônio será de 410.

Sendo definido o valor da perturbação a ser usado no cálculo dos coeficientes de

sensibilidade por diferenças finitas, procede-se agora às análises de dependência

linear e magnitude dos mesmos.

60

Os parâmetros do modelo de HODGKIN-HUXLEY aplicados ao axônio foram

apresentados nas equações 3.11 e 3.15-3.20. As equações a seguir são escritas com

os parâmetros onde serão aplicadas incertezas.

0.1 25

exp 0.1 25 1 exp

m m m m

m m

m m m m m

V a b V

V c d V e

(5.1)

4exp exp18

m mm m m

m

V Vf

g

(5.2)

0.07exp exp20

m mh h h

h

V Va

b

(5.3)

1 1

exp 0.1 30 1 exph h

m h h m hV c d V e

(5.4)

100.01

exp 0.1 10 1 exp

m n m

n n n

m n n m n

V b Va

V c d V e

(5.5)

As Figuras 5.11.a - 5.14.a apresentam os coeficientes de sensibilidade reduzidos

para o axônio. Verificou-se que os coeficientes apresentam magnitude elevada,

apenas entre 27 e 31 ms. Por isso, os coeficientes de sensibilidade também foram

apresentados apenas neste intervalo de tempo nas figuras 5.11.b - 5.14.b.

61

0 10 20 30 40 50 60

-600

-400

-200

0

200

400

600m

V

Tempo (ms)

Coeficiente de Sensibilidade Reduzido

Vm

XC

m

XG

máx

Na

XG

máx

k

XG

L

XV

Na

XV

k

XV

L

(a)

27.5 28 28.5 29 29.5 30 30.5 31 31.5 32

-100

-50

0

50

100

150

mV

Tempo (ms)


Vm

XC

m

XG

máx

Na

XG

máx

k

XG

L

XV

Na

XV

k

XV

L

(b)

Figura 5.11 - Análise de mCX , máx

NaGX , máx

kGX ,

LGX , NaVX ,

kVX e LVX para o axônio (a)

domínio completo e (b) região de interesse.

62

0 10 20 30 40 50 60

-5000

-4000

-3000

-2000

-1000

0

1000

2000

mV

Tempo (ms)


Vm

Xa

m

Xb

m

Xc

m

Xd

m

Xe

m

Xfm

Xg

m

(a)

27.5 28 28.5 29 29.5 30 30.5 31

-400

-200

0

200

400

600

800

1000

mV

Tempo (ms)


Vm

Xa

m

Xb

m

Xc

m

Xd

m

Xe

m

Xfm

Xg

m

(b)

Figura 5.12 - Análise dos coeficientes de sensibilidade dos parâmetros de m e m

para o axônio (a) domínio completo e (b) região de interesse.

63

0 10 20 30 40 50 60-100

0

100

200

300

400

500

600

700

mV

Tempo (ms)


Vm

Xa

h

Xb

h

Xc

h

Xd

h

Xe

h

(a)

27 27.5 28 28.5 29 29.5 30 30.5 31

-100

-80

-60

-40

-20

0

20

40

60

80

100

mV

Tempo (ms)


Vm

Xa

h

Xb

h

Xc

h

Xd

h

Xe

h

(b)

Figura 5.13 - Análise dos coeficientes de sensibilidade dos parâmetros de h e h


64

0 10 20 30 40 50 60

-1500

-1000

-500

0

500

1000

1500

2000

2500

mV

Tempo (ms)


Vm

Xa

n

Xb

n

Xc

n

Xd

n

Xe

n

Xfn

Xg

n

(a)

27.5 28 28.5 29 29.5 30 30.5-600

-500

-400

-300

-200

-100

0

100

200

300

400

mV

Tempo (ms)


Vm

Xa

n

Xb

n

Xc

n

Xd

n

Xe

n

Xfn

Xg

n

(b)

Figura 5.14 - Análise dos coeficientes de sensibilidade dos parâmetros de n e n


65

Nas Figuras 5.11.b-5.14.b pode-se perceber que todos os coeficientes de

sensibilidade de cada grupo analisado são linearmente dependentes. Foram então

selecionados os parâmetros que apresentaram maiores magnitudes de cada grupo,

para fazer uma análise entre os mesmos a fim de escolher quais têm maiores

magnitudes e se há independência linear entre eles. Esta comparação é apresentada

na Figura 5.15.

0 10 20 30 40 50 60

-5000

-4000

-3000

-2000

-1000

0

1000

2000

mV

Tempo (ms)


Vm

XC

m

Xd

m

Xc

h

Xd

n

(a)

66

27 27.5 28 28.5 29 29.5 30 30.5

-500

0

500

1000

mV

Tempo (ms)


Vm

XC

m

Xd

m

Xc

h

Xd

n

(b)

Figura 5.15 - Análise da dependência linear entre mCX ,

mdX , hcX e

ndX para o axônio

(a) domínio completo e (b) região de interesse.

Analisando a Figura 5.15 pode-se notar a dependência linear entre todos os

parâmetros. Pode-se também verificar que o parâmetro que apresenta maior

magnitude para o coeficiente de sensibilidade reduzido é mdX . Por outro lado, optou-

se neste trabalho pela estimativa da capacitância, mC , por apresentar uma magnitude

considerável e principalmente por dar uma informação de como está o acúmulo de

cargas nas superfícies da membrana celular.

O mesmo procedimento realizado para a análise dos coeficientes de sensibilidade

para o axônio será feito para aplicação em fibras de Purkinje. Os parâmetros do

modelo de HODGKIN-HUXLEY aplicados em fibras de Purkinje foram apresentados

nas equações 3.21-3.32 e as equações a seguir são escritas com os parâmetros que

serão aplicados incertezas.

67

3 3

, , ,0.14l

máx máx

Na P Na P P p Na P P p NaG G m h G m h G (5.6)

, , , ,

1

90 901.2exp 0.015exp exp exp

50 60

m P m P m P k m P k

k k k

k k

V V V c V fG a e

d g

(5.7)

4 4

2 1.2k P k PG n b n (5.8)

, , ,

, ,

, , ,

,

,

480.1

48exp 1 exp

15

m P m P m P

m P m P

m P m P m P

m P

m P

V V ba

V V ce

d

(5.9)

, , , ,

,

, , ,

,

,

0.12 8

8exp 1 exp

5

m P m P m P m P

m P

m P m P m P

m P

m P

V f V g

V V hj

i

(5.10)

, , ,

, ,

,

900.17exp exp

20

m P m P h P

h P h P

h P

V V ba

c

(5.11)

11

, , ,

, ,

,

42exp 1 exp

10

m P m P h P

h P h P

h P

V V eg

f

(5.12)

4

, , , ,

,

, , ,

,

,

10 50

50exp 1 exp

10

m P n P m P n P

n P

m P m P n P

n P

n P

V a V b

V V ce

d

(5.13)

, , ,

, ,

,

900.002exp exp

80

m P m P n P

n P n P

n P

V V gf

h

(5.14)

As Figuras 5.16 - 5.21 apresentam em todo o tempo analisado os coeficientes de

sensibilidade reduzidos para fibras de Purkinje.

68

0 50 100 150 200 250 300 350 400 450 500

-5000

-4000

-3000

-2000

-1000

0

1000

mV

Tempo (ms)

Vm,P

X C

m,P

X G

máx

Na,P

X G

Nal

X V

Na,P

X V

k,P

(a)

310 320 330 340 350 360 370 380 390 400-800

-600

-400

-200

0

200

400

mV

Tempo (ms)

Vm,P

X C

m,P

X G

máx

Na,P

X G

Nal

X V

Na,P

X V

k,P

(b)

Figura 5.16 - Análise de ,m PCX ,

,máxNaG P

X , NalGX ,

,Na PVX e ,k PVX para fibras de Purkinje (a)

domínio completo e (b) região de interesse.

69

0 50 100 150 200 250 300 350 400 450 500

-1500

-1000

-500

0

500

1000

1500

2000

2500

3000

mV

Tempo (ms)

Vm,P

X a

k

X b

k

X c

k

X d

k

X e

k

X fk

X g

k

(a)

300 310 320 330 340 350 360 370 380 390

-300

-200

-100

0

100

200

300

400

mV

Tempo (ms)

Vm,P

X a

k

X b

k

X c

k

X d

k

X e

k

X fk

X g

k

(b)

Figura 5.17 - Análise dos coeficientes de sensibilidade dos parâmetros de 1kG e 2kG

para fibras de Purkinje (a) domínio completo e (b) região de interesse.

70

0 50 100 150 200 250 300 350 400 450 500

-1

-0.5

0

0.5

1

1.5

2

x 104

mV

Tempo (ms)

Vm,P

X a

m,P

X b

m,P

X c

m,P

X d

m,P

X e

m,P

(a)

310 320 330 340 350 360 370 380 390

-4000

-3000

-2000

-1000

0

1000

2000

3000

4000

mV

Tempo (ms)

Vm,P

X a

m,P

X b

m,P

X c

m,P

X d

m,P

X e

m,P

(b)

Figura 5.18 - Análise dos coeficientes de sensibilidade dos parâmetros de ,m P e ,m P


71

0 50 100 150 200 250 300 350 400 450 500

-2500

-2000

-1500

-1000

-500

0

500

1000

mV

Tempo (ms)

Vm,P

X a

h,P

X b

h,P

X c

h,P

X d

h,P

X e

h,P

X fh,P

X g

h,P

(a)

280 300 320 340 360 380 400

-1500

-1000

-500

0

500

1000

mV

Tempo (ms)

Vm,P

X a

h,P

X b

h,P

X c

h,P

X d

h,P

X e

h,P

X fh,P

X g

h,P

(b)

Figura 5.19 - Análise dos coeficientes de sensibilidade dos parâmetros de ,h P e ,h P


72

0 50 100 150 200 250 300 350 400 450 500

-400

-300

-200

-100

0

100

200

300

mV

Tempo (ms)

Vm,P

X a

n,P

X b

n,P

Xc

n,P

X d

n,P

X e

n,P

X fn,P

X g

n,P

X h

n,P

Figura 5.20 - Análise dos coeficientes de sensibilidade dos parâmetros de ,n P e ,n P

para fibras de Purkinje.

Nas Figuras 5.16.b-5.20.b pode-se perceber que todos os coeficientes de

sensibilidade de cada grupo analisado são linearmente dependentes. Foram, então,

selecionados os parâmetros que apresentaram maiores magnitudes de cada grupo,

para fazer uma análise entre os mesmos a fim de escolher quais têm maiores

magnitudes e se há independência linear entre eles. Esta comparação é apresentada

na Figura 5.21.

73

290 300 310 320 330 340 350 360 370 380 390

-1500

-1000

-500

0

500

1000

1500

mV

Tempo (ms)

Vm,P

X c

k,P

X jm,P

Xb

h,P

X b

n,P

X V

k,P

Figura 5.21 - Análise da dependência linear entre ,k PVX ,

kcX , ,m PjX ,

,h PbX e ,n PbX em

uma região específica.para fibras de Purkinje.

Analisando as Figuras 5.16-5.21 pode-se verificar que o parâmetro que apresenta

maior magnitude para o coeficiente de sensibilidade reduzido é ,h PbX . Por outro lado,

optou-se por estimar o parâmetro relacionado à capacitância, ,k PV , por apresentar uma

magnitude considerável e principalmente por dar uma informação de como está o

potencial de repouso do íon de potássio. Por este motivo para o modelo de HODGKIN-

HUXLEY aplicado em fibras de Purkinje o parâmetro estimado seria o potencial de

repouso do íon de potássio, ,k PV .

Além da análise de magnitude e dependência linear dos coeficientes de

sensibilidade reduzidos, é interessante observar o intervalo de tempo em que há

grandes magnitudes para os coeficientes de sensibilidade. Caso este intervalo de

tempo seja muito curto, pode ser que poucas medidas tragam informação para a

74

estimativa dos parâmetros de interesse. Isto, de fato, se passou com o parâmetro ,k PV

em fibras de Purkinje, como será apresentado na seção seguinte.

5.3 ESTIMATIVA DE PARÂMETROS

Nesta seção, serão apresentados os resultados obtidos para a estimativa dos

parâmetros do modelo de HODGKIN-HUXLEY, com o método MCMC descrito no

capítulo 4.

Para a proposta de candidatos na cadeia de Markov foi usado um modelo de

passeio aleatório com uma distribuição Gaussiana, na seguinte forma:

1* .t

w

P P ε (5.15)

onde ~ε tem média zero, é não-correlacionado e os elementos da matriz de covariância

dada por:

2

,teórico,

se cov

0 se

i ii j

P i j

i j

ε (5.16)

Ao estimar os parâmetros com o método MCMC, uma etapa importante é

determinar qual valor do passo de procura, w , será usado. Escolhendo-se um valor

muito baixo para w , a convergência para estimativa será bastante lenta. Caso se

escolha um valor alto para w , serão necessárias várias iterações para o método

aceitar novos parâmetros. Estes fatos podem ser verificados analisando-se a taxa de

aceitação.

A análise de w foi realizada tomando distribuições Gaussianas com desvio-

padrão constante para as distribuições a priori e com os demais dados apresentados

75

na Tabela 5.1. Foram analisados os seguintes valores para o passo de procura:

0.0005, 0.001 e 0.0015, para o caso de aplicação em axônios, e 0.0001, 0.003 e

0.0005 para o caso de aplicação em fibras de Purkinje. Os resultados obtidos para as

cadeias de Markov de mC nestes casos são mostrados nas Figuras 5.22-5.23.

Tabela 5.1 - Dados para analisar w .

AXÔNIO PURKINJE

Parâmetros Valores Parâmetros Valores

esti

priori

P 0.2 mC esti

priori

P 0.2 kV

prioriP(i)

0.01 iP prioriP(i)

0.01 iP

med ,0.05 m máxV med ,0.05 m máxV

N (Axônio) 2. 105 N (Purkinje) 104

Onde esti

priori

P representa o desvio-padrão para a distribuição a priori do parâmetro a

ser estiado, priori

P(i) é o desvio-padrão para a distribuição a priori dos demais

parâmetros, med é o desvio-padrão para as medidas simuladas e N é a quantidade

de estágios na cadeia de Markov.

76

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 105

0.9

0.95

1

1.05

1.1

1.15

1.2

1.25

N (estado)

Cm

(

F)

Teórico

w = 0.0015

w = 0.001

w = 0.0005

Figura 5.22 - Análise de w , estimando mC F para axônio.

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000-120

-115

-110

-105

-100

-95

-90

Vk (

mV

)

N (estado)

Exato

w = 0.0001

w = 0.0003

w = 0.0005

Figura 5.23 - Análise de w , estimando kV mV para fibras de Purkinje.

77

Para análise realizada em axônio com 0.0005w , 0.001w e 0.0015w foram

encontradas as taxas de aceitação 12.55%, 4.27% e 4.17%, respectivamente. Para

análise realizada em fibras de Purkinje com 0.0001w , 0.0003w e 0.0005w

foram encontradas as taxas de aceitação 7.15%, 6.20% e 3.07%, respectivamente.

Além da taxa de aceitação reduzir com o aumento do passo de procura, nota-se a

redução do número de estados de cadeias de Markov necessários para a

convergência. Pode-se verificar nas Figuras 5.22 e 5.23 que, para passos de procura

maiores, são necessários menos estados para a convergência da cadeia de Markov.

Pode-se perceber que os estados de aquecimento para o caso em que se analisou

axônio, para passos de procura 0.0005w , 0.001w e 0.0015w , são de 1.6x105,

6x104 e 2x104 estados, respectivamente. Pode-se perceber que os estados de

aquecimento para o caso em que se analisou fibras de Purkinje, para passos de

procura 0.0001w , 0.003w e 0.0005w , são de 4500, 6600 e 700 estados,

respectivamente

Os casos analisados para estimar o parâmetro mC para o axônio e kV para fibras

de Purkinje são apresentados na Tabela 5.2, os quais envolvem diferentes desvios

padrão para as medidas e para a distribuição a priori dos parâmetros a serem

estimados.

78

Tabela 5.2 - Casos analisados para a estimativa de parâmetros.

Aplicação Caso priori med

Axônio 1 0.1 mC ,0.1 m máxV




Fibras de Purkinje 5 0.1 kV ,0.1 m máxV




Para os demais parâmetros foram utilizados desvios padrão de 1 % dos valores

teóricos, com intuito apenas de levar em consideração pequenas incertezas

associadas a cada um. Para o primeiro estado da cadeia de Markov, foram utilizados

os respectivos valores teóricos, para cada parâmetro, exceto para mC , que para o

axônio foi inicializado com um valor 20% acima do valor teórico e para kV , que para

fibras de Purkinje foi inicializado com um valor 20% inferior do valor teórico. Os

resultados obtidos para os 8 casos analisados são apresentados nas Figura 5.24 e

5.25.

79

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 105

0.8

0.9

1

1.1

1.2

1.3

N (estados)

Cm

Teórico

Caso 1

Caso 2

Caso 3

Caso 4

Figura 5.24 - Estimativa de mC F para o axônio.

Analisando as Figura 5.24 pode-se verificar que o caso 3 não resulta em uma boa

estimativa, pois os desvios-padrão para as medidas e para o parâmetro são 20% do

potencial máximo e de mC , respectivamente. A convergência é atingida com 40000,

60000 e 100000 para os casos 1, 4 e 2 respectivamente. Os valores medidos

estimados para o parâmetro mC são apresentados na Tabela 5.3, juntamente com

seus intervalos de 99% de confiança e com erro relativo referente à medida estimada.

Analisando a Tabela 5.3, pode-se verificar que o caso 3 resulta na pior estimativa.

Por outro lado, a melhor estimativa foi encontrada para o caso 2, pois para este caso

os desvios para as medidas e para o parâmetro são os menores. Os erros relativos

encontrados foram 2.00%, 0.00%, 20.00% e 1% para os casos 1, 2, 3 e 4

respectivamente. Além de verificar os erros relativos é importante observar o intervalo

de confiança. Analisando este, pode-se notar que apesar de o caso 3 não apresenta

uma boa estimativa, por outro lado os casos 2 e 4 o valor teórico está no intervalo de

80

confiança, analisando o erro relativo e o intervalo de confiança de 99%, pode-se

concluir que para estes casos obteve-se uma boa estimativa para o parâmetro mC .

Analisando a Tabela 5.4, pode-se verificar que os casos em que foram analisadas

as fibras de Purkinje o valor teórico do parâmetro estimado não está no intervalo de

confiança e os erros relativos encontrados foram de 6.00, 6.25, 8.20 e 9.20 %. Tal fato

também fica evidenciado pela análise da Figura 5.25. Pode-se concluir que as

estimatvas de ,k PV não foram boas, por isso não foi possível fazer estimativas de

parâmetros para o modelo de HODGKIN-HUXLEY aplicado em fibras de Purkinje.

Tabela 5.3 - Estimativa de mC F para o axônio.

Caso priori med ( )mC teórico ( )mC estimado Erro relativo (%)

1 0.1 mC ,m x0.1 m áV 1 0.98±0.03 2.00

2 0.1 mC ,m x0.05 m áV 1 1.00±0.02 0.00

3 0.2 mC ,m x0.1 m áV 1 1.20±0.01 20.00

4 0.2 mC ,m x0.05 m áV 1 1.01±0.03 1.00

Tabela 5.4 - Estimativa de kV mV para fibras de Purkinje.

Caso priori med ( )kV teórico ( )kV estimado Erro relativo (%)

5 0.1 kV ,0.1 m máxV -100 -106.48±0.05 6.00

6 0.1 kV ,0.05 m máxV -100 -106.25±3.86 6.25

7 0.2 kV ,0.1 m máxV -100 -108.20±0.02 8.20

8 0.2 kV ,0.05 m máxV -100 -109.20±0.03 9.20

81

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000-120

-115

-110

-105

-100

-95

-90

N (estado)

Estimativa para Vk

Exato

Caso 5

Caso 6

Caso 7

Caso 8

Figura 5.25 - Estimativa de kV mV para fibras de Purkinje.

5.4 ESTIMATIVA DAS VARIÁVEIS DE ESTADO

As estimativas para as variáveis de estado foram realizadas com as seguintes

metodologias Bayesianas do filtro de partículas: filtro de Amostragem por Importância

e Reamostragem Sequencial (Sampling Importance Resampling - SIR) e filtro Auxiliar

de Amostragem por Importância e Reamostragem Sequencial (Axiliary Sampling

Importance Resampling - ASIR). Uma análise a ser feita aqui será a comparação entre

estes filtros, aplicados para estimar as variáveis de estado no caso de aplicação do

modelo de HODGKIN-HUXLEY em axônios e fibras de Purkinje.

Neste trabalho, as incertezas dos modelos de evolução de estado e de

observação serão consideradas aditivas e Gaussianas.

82

As variáveis de estado serão estimadas a cada 1 s , nos instantes de tempo que

coincidem com os de medidas, para o axônio. Para as fibras de Purkinje se trabalhará

com intervalos de 10 s .

A comparação entre os filtros será realizada utilizando desvios-padrão de 5% do

potencial para as medidas e também de 5% da condição inicial para o modelo.

Para escolher qual tipo de filtro de partículas é melhor para ser usado, um

parâmetro utilizado para a comparação foi o erro RMS , cuja expressão é apresentada

a seguir:

2

1

1 M

estimado teóricoiRMS X i X i

M (5.17)

onde M é a quantidade de medidas.

Os filtros de partículas são métodos que heurísticos e por isso os resultados

podem não ser o mesmo quando forem realizadas repetições, pelo fato de os filtros

apresentarem essa característica foram realizadas 100 repetições em todos os casos

estudados e os resultados foram considerados como sendo a média entre os valores

obtidos nas 100 repetições.

Ao usar os filtros de partículas é importante avaliar o número mínimo de partículas

necessário para representar adequadamente as distribuições a posteriori das variáveis

de estado do problema proposto. As Figuras (5.26-5.29) apresentam as estimativas do

coeficiente de inativação h para o axônio usando 500, 750, 1000 e 5000 partículas,

respectivamente, utilizando o filtro de partículas SIR. Nestas Figuras são apresentadas

as médias e os intervalos de confiança de 99% da distribuição a posteriori das

variáveis de estado em cada tempo.

83

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


partículas.

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


partículas.

84

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


partículas.

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


partículas.

85

Analisando as Figuras (5.26-5.29), verifica-se que quanto maior for a quantidade

de partículas utilizadas, as estimativas das variáveis de estado se tornam melhores.

Além dos intervalos de confiança ficarem melhor definidos. Uma análise análoga é

agora realizada utilizando o filtro de partículas ASIR. As Figuras (5.30-5.33) mostram

os resultados obtidos para a estimativa de h usando 50, 75, 100 e 500 partículas.

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


partículas.

86

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


partículas.

87

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


partículas.

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


partículas.

88

Analisando as Figuras (5.30-5.33), de forma semelhante ao que aconteceu com o

filtro de partículas SIR, verifica-se que quanto maior for o mínimo de partículas

utilizadas, as estimativas das variáveis de estado se tornam melhores, além dos

intervalos de confiança ficarem melhor definidos.

A comparação das Figuras (5.26-5.29) com as Figuras (5.30-5.33) revela

claramente que o filtro ASIR requer um número bem menor de partículas, para que

sejam obtidos resultados semelhantes ao do filtro SIR com um grande número de

partículas. Embora esta análise tenha sido realizada para a variável de estado, h ,

resultados semelhantes são observados para as outras variáveis de estado. Além dos

gráficos apresentados, são apresentados na Tabela 5.5 os erros RMS e entre

parênteses o desvio padrão da estimativa de cada uma das variáveis de estado, com

os filtros SIR e ASIR, e os diferentes números de partículas analisados com os

mesmos desvios-padrão de 5% do potencial para as medidas e também de 5% da

condição inicial para o modelo. Nesta Tabela também são apresentados os tempos

computacionais de cada caso.

Tabela 5.5 - Comparação entre os filtros SIR e ASIR para o caso do axônio.

Filtro N Tempo

Cpu (s) mRMS V mV RMS h RMS m RMS n

SIR 500 137.237 1.624 (0.049) 0.019 (0.002) 0.017 (0.001) 0.010 (0.001)

SIR 750 198.436 2.779 (0.090) 0.029 (0.003) 0.020 (0.001) 0.016 (0.001)

SIR 1000 267.217 1.557 (0.061) 0.015 (0.001) 0.010 (0.001) 0.010 (0.001)

SIR 5000 1730.591 1.361 (0.008) 0.010 (0.001) 0.010 (0.001) 0.010 (0.001)

ASIR 50 29.339 0.818 (0.079) 0.014 (0.003) 0.008 (0.002) 0.007 (0.001)

ASIR 75 44.258 1.300 (0.109) 0.010 (0.002) 0.012 (0.002) 0.006 (0.002)

ASIR 100 58.819 1.456 (0.130) 0.018 (0.004) 0.012 (0.001) 0.007 (0.001)

ASIR 500 313.231 1.145 (0.040) 0.012 (0.001) 0.012 (0.001) 0.007 (0.001)

89

Analisando a Tabela 5.5, pode-se verificar que, de fato o filtro de partículas ASIR,

com um número de partículas dez vezes menor e tempo computacional bastante

inferior é capaz de apresentar erros RMS para as variáveis de estado semelhantes ou

até mesmo menores que os do filtro de partículas SIR, em relação aos desvios-padrão

observou-se que os dois filtros apresentaram baixos desvios-padrão para os casos

analisados.

A comparação entre os filtros SIR e ASIR é agora realizada para as fibras de

Purkinje. Nesta comparação são usadas incertezas com desvios-padrão de 5% das

variáveis de estado na condição inicial e 5% do potencial máximo obtido na solução

direta para as medidas, e aquisição de dados a cada 10 s .

As Figuras (5.34-5.37) apresentam as estimativas de n para fibras de Purkinje

usando 500, 750, 1000 e 5000 partículas utilizando o filtro de partículas SIR.

100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança


com 500 partículas.

90

100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança



100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança



91

100 200 300 400 500 600 700 800 900 10000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança



Analisando as Figuras (5.35-5.37), verifica-se que quanto maior for a quantidade

de partículas utilizadas, as estimativas das variáveis de estado se tornam melhores,

além dos intervalos de confiança ficarem melhor definidos. A mesma análise é agora

realizada utilizando-se o filtro de partículas ASIR. As Figuras (5.38-5.41) mostram os

resultados obtidos para a estimativa de n usando 50, 75, 100 e 500 partículas.

92

100 200 300 400 500 600 700 800 900 10000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n


ASIR com 50 partículas.

100 200 300 400 500 600 700 800 900 10000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n



93

100 200 300 400 500 600 700 800 900 10000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n



100 200 300 400 500 600 700 800 900 10000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n



94

Analisando as Figuras (5.38-5.41), de forma semelhante ao filtro de partículas SIR,

verifica-se que quanto maior for o número de partículas utilizadas, as estimativas das

variáveis de estado se tornam melhores, além dos intervalos de confiança ficarem

melhor definidos. Pode-se verificar que, para os casos em que foram usadas 50 e 75

partículas, os valores teóricos não estão no intervalo de confiança.

A comparação das Figuras (5.34-5.37) com as Figuras (5.38-5.41) revela

claramente que o filtro ASIR requer um número bem menor de partículas, para que

sejam obtidos resultados semelhantes ao do filtro SIR com um grande número de

partículas. Embora esta análise tenha sido realizada para a variável de estado, n ,

resultados semelhantes são observados para as outras variáveis de estado.

Apresenta-se na Tabela 5.5 os erros RMS e desvios-padrão para as estimativa de

cada uma das variáveis de estado, com os filtros SIR e ASIR, e os diferentes números

de partículas analisados. Nesta Tabela também são apresentados os tempos

computacionais de cada caso.

Tabela 5.6- Comparação entre os filtros SIR e ASIR para o caso de fibras de Purkinje.

Filtro Np Tempo Cpu (s)

mRMS V mV RMS h RMS m RMS n

SIR 500 662.49 1.265 (0.076) 0.016 (0.002) 0.011 (0.002) 0.024 (0.005)

SIR 750 1629.4 1.280 (0.095) 0.017 (0.003) 0.010 (0.001) 0.021 (0.004)

SIR 1000 1337.4 1.272 (0.094) 0.017 (0.002) 0.011 (0.001) 0.020 (0.002)

SIR 5000 6985.2 1.288 (0.122) 0.017 (0.001) 0.011 (0.002) 0.019 (0.002)

ASIR 50 163.67 1.127 (0.042) 0.016 (0.002) 0.011 (0.001) 0.041 (0.015)

ASIR 75 246.16 1.228 (0.033) 0.018 (0.002) 0.010 (0.001) 0.036 (0.012)

ASIR 100 321.79 1.138 (0.032) 0.017 (0.002) 0.010 (0.001) 0.034 (0.013)

ASIR 500 1656.1 1.370 (0.024) 0.018 (0.001) 0.010 (0.001) 0.029 (0.010)

95

Analisando a Tabela 5.5, pode-se verificar que o filtro de partículas ASIR, com um

número de partículas dez vezes menor e tempo computacional bastante inferior, é

capaz de apresentar RMS para as variáveis de estado semelhantes ao encontrado

com o filtro de partículas SIR.

Após realizar comparações entre os filtros de partículas SIR e ASIR para estimar

as variáveis de estado do modelo de HODGKIN-HUXLEY com aplicação em axônios e

fibras de Purkinje, conclui-se que é melhor usar o filtro de partículas ASIR para os dois

casos. Assim, serão analisados agora casos com o uso do filtro de partículas ASIR e

500 partículas, com desvios-padrão, 1% e 5% do valor máximo de mV para as

medidas e 1% e 5% da condição inicial para o modelo. Os casos analisados são

apresentados na Tabela 5.7.

Tabela 5.7- Casos analisados para estimar variáveis de estado.

Casos Aplicação med modelo Np t ms

9 axônio ,0.01 m máxV 00.01x 500 1

10 axônio ,0.01 m máxV 00.05x 500 1

11 axônio ,0.05 m máxV 00.01x 500 1

12 axônio ,0.05 m máxV 00.05x 500 1

13 fibras de Purkinje ,0.01 m máxV 00.01x 500 10




As Figuras (5.42-5.49) apresentam os resultados para a estimativa das variáveis

de estado obtidos para o caso 9. Como esperado, a concordância entre as variáveis

96

de estado estimadas e exata é excelente para este caso, já que as incertezas

associadas aos modelos de evolução e de observação são pequenas.

0 10 20 30 40 50 60-20

0

20

40

60

80

100

120

Tempo (ms)

Vm

(m

V)

Exata

Medidas

Estimada

99% Confiança

Figura 5.42 - Estimativa de mV de axônio para o caso 9.

97

0 10 20 30 40 50 600

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação m

Exata

Estimada

99% Confiança

Figura 5.43 - Estimativa de m de axônio para o caso 9.

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança

Figura 5.44 - Estimativa de h de axônio para o caso 9.

98

0 10 20 30 40 50 600.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança

Figura 5.45 - Estimativa de n de axônio usando para o caso 9.

0 10 20 30 40 50 60-10

0

10

20

30

40

50

60

Tempo (ms)

GN

a (m

S)

Exata

Estimada

99% Confiança

Figura 5.46 - Estimativa de NaG de axônio para o caso 9.

99

0 10 20 30 40 50 600

2

4

6

8

10

12

Tempo (ms)

Gk

(mS

)

Exata

Estimada

99% Confiança

Figura 5.47 - Estimativa de kG de axônio usando para o caso 9.

0 10 20 30 40 50 60-1200

-1000

-800

-600

-400

-200

0

200

Tempo (ms)

I Na

(mA

)

Exata

Estimada

99% Confiança

Figura 5.48 - Estimativa de NaI de axônio para o caso 9.

100

0 10 20 30 40 50 60-100

0

100

200

300

400

500

600

700

800

900

Tempo (ms)

I k (m

A)

Exata

Estimada

99% Confiança

Figura 5.49 - Estimativa de kI de axônio para o caso 9.


de estado obtidos para o caso 10.

101

0 10 20 30 40 50 60-20

0

20

40

60

80

100

120

Tempo (ms)

Vm

(m

V)

Exata

Medidas

Estimada

99% Confiança


0 10 20 30 40 50 600

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação m

Exata

Estimada

99% Confiança


102

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


0 10 20 30 40 50 600.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança

Figura 5.53 - Estimativa de n de axônio para o caso 10.

103

0 10 20 30 40 50 60-10

0

10

20

30

40

50

60

Tempo (ms)

GN

a (m

S)

Exata

Estimada

99% Confiança


0 10 20 30 40 50 600

2

4

6

8

10

12

Tempo (ms)

Gk

(mS

)

Exata

Estimada

99% Confiança

Figura 5.55 - Estimativa de kG de axônio para o caso 10.

104

0 10 20 30 40 50 60-1200

-1000

-800

-600

-400

-200

0

200

Tempo (ms)

I Na

(mA

)

Exata

Estimada

99% Confiança


0 10 20 30 40 50 60-100

0

100

200

300

400

500

600

700

800

900

Tempo (ms)

I k (m

A)

Exata

Estimada

99% Confiança


105

A diferença entre os casos 9 e 10 está no desvio padrão do modelo de evolução

de estado, onde, para o caso 9 é de 1% da condição inicial e para o caso 10 é de 5%

da condição inicial. Pode-se observar que para o caso 10 Figuras (5.50-5.57) que o

intervalo de confiança é maior e não é tão regular como no caso 9 Figuras (5.42-5.49)

para a estimativa das variáveis de estado. Isto se dá particularmente para as variáveis

de estado para as quais não há medições, como, m , n , h , NaG ,

kG , NaI e kI .


de estado obtidos para o caso 11. A comparação da Figura 5.42 com as Figuras 5.50

e 5.58 mostram as maiores incertezas das medidas de mV para o caso 11.

0 10 20 30 40 50 60-20

0

20

40

60

80

100

120

Tempo (ms)

Vm

(m

V)

Exata

Medidas

Estimada

99% Confiança


106

0 10 20 30 40 50 600

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação m

Exata

Estimada

99% Confiança


0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


107

0 10 20 30 40 50 600.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança


0 10 20 30 40 50 60-10

0

10

20

30

40

50

60

Tempo (ms)

GN

a (m

S)

Exata

Estimada

99% Confiança


108

0 10 20 30 40 50 600

2

4

6

8

10

12

Tempo (ms)

Gk

(mS

)

Exata

Estimada

99% Confiança


0 10 20 30 40 50 60-1200

-1000

-800

-600

-400

-200

0

200

Tempo (ms)

I Na

(mA

)

Exata

Estimada

99% Confiança


109

0 10 20 30 40 50 60-100

0

100

200

300

400

500

600

700

800

900

Tempo (ms)

I k (m

A)

Exata

Estimada

99% Confiança


Pode-se observar na Figura (5.42) que as estimativas de mV se aproximam mais

que os valores medidos dos valores teóricos. Analisando as Figuras (5.58-5.65) pode-

se verificar que os valores teóricos estão dentro do intervalo de confiança de 99% para

a estimativa para as variáveis de estado. Apesar das maiores incertezas das medições

neste caso, quando comparadas ao caso 10 Figuras (5.50-5.57) nota-se nas Figuras

(5.58-5.65) uma excelente concordância das variáveis de estado estimadas com seus

valores exatos.



110

0 10 20 30 40 50 60-20

0

20

40

60

80

100

120

Tempo (ms)

Vm

(m

V)

Exata

Medidas

Estimada

99% Confiança


0 10 20 30 40 50 600

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação m

Exata

Estimada

99% Confiança


111

0 10 20 30 40 50 600

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


0 10 20 30 40 50 600.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança


112

0 10 20 30 40 50 60-10

0

10

20

30

40

50

60

Tempo (ms)

GN

a (m

S)

Exata

Estimada

99% Confiança


0 10 20 30 40 50 600

2

4

6

8

10

12

Tempo (ms)

Gk

(mS

)

Exata

Estimada

99% Confiança


113

0 10 20 30 40 50 60-1200

-1000

-800

-600

-400

-200

0

200

Tempo (ms)

I Na

(mA

)

Exata

Estimada

99% Confiança


0 10 20 30 40 50 60-100

0

100

200

300

400

500

600

700

800

900

Tempo (ms)

I k (m

A)

Exata

Estimada

99% Confiança


114

Analisando as Figuras (5.66-5.73) pode-se verificar que o intervalo de confiança

de 99 % de confiança aumenta significativamente no caso 12 em relação aos casos

anteriores. Por outro lado, nota-se nas Figuras (5.66-5.73) uma excelente

concordância entre as medidas estimadas para as variáveis de estado e seus valores

exatos, apesar do caso 12 envolver grandes incertezas nos modelos de evolução de

estado e de observação.

Os erros RMS e entre parênteses os desvios-padrão para a estimativa de cada

uma das variáveis de estado, para os casos 9 a 12, são apresentados na Tabela 5.8.

A análise desta Tabela revela que as estimativas são mais afetadas pelas incertezas

nas medidas experimentais do que nos modelos de evolução de estado.

Tabela 5.8 - Valores de RMS para as variáveis de estado para os casos 9, 10, 11 e

12.

Variável Caso 9 Caso 10 Caso 11 Caso 12

NaG S 0.047 (0.015) 0.238 (0.135) 0.786 (0.071) 0.249 (0.057)

kG S 0.021 (0.006) 0.130 (0.049) 0.095 (0.013) 0.149 (0.020)

NaI mA 3.311 (0.394) 5.579 (4.137) 10.589 (0.618) 10.892 (1.038)

kI mA 1.995 (0.544) 4.947 (3.229) 10.210 (1.626) 10.356 (1.633)

mV mV 0.210 (0.008) 0.600 (0.011) 1.002 (0.070) 1.145 (0.040)

h 0.002 (0.001) 0.009 (0.001) 0.003 (0.001) 0.012 (0.001)

m 0.002 (0.001) 0.004 (0.001) 0.018 (0.001) 0.023 (0.001)

n 0.001 (0.001) 0.004 (0.001) 0.003 (0.001) 0.007 (0.001)

Pode-se observar na Tabela 5.8 que os valores de erro RMS são os menores

para o caso 9. Este fato já era esperado já que este caso tem as menores incertezas

115

para as medidas e o modelo de evolução de estado. Por outro lado, o caso 12

apresentou os maiores valores de erro RMS , já que envolve os maiores valores para

as incertezas de medidas e do modelo de evolução de estado.

A mesma análise realizada para o modelo de HODGKIN-HUXLEY aplicado em

axônios é agora feita para fibras de Purkinje. As Figuras (5.74-5.81) apresentam os

resultados para as estimativas das variáveis de estado obtidas para o caso 13.

0 100 200 300 400 500 600 700 800 900 1000-100

-80

-60

-40

-20

0

20

Tempo (ms)

Vm

(m

V)

Exata

Medidas

Estimada

99% Confiança

Figura 5.74 - Estimativa de mV de fibras de Purkinje para o caso 13.

116

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação m

Exata

Estimada

99% Confiança

Figura 5.75 - Estimativa de m de fibras de Purkinje para o caso 13.

0 100 200 300 400 500 600 700 800 900 10000.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança

Figura 5.76 - Estimativa de n de fibras de Purkinje para o caso 13.

117

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança

Figura 5.77 - Estimativa de h de fibras de Purkinje para o caso 13.

0 100 200 300 400 500 600 700 800 900 1000-4

-2

0

2

4

6

8

10

12

14

16

Tempo (ms)

GN

a (m

S)

Exata

Estimada

99% Confiança

Figura 5.78 - Estimativa de NaG de fibras de Purkinje para o caso 13.

118

0 100 200 300 400 500 600 700 800 900 10000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Tempo (ms)

Gk

(mS

)

Exata

Estimada

99% Confiança

Figura 5.79 - Estimativa de kG de fibras de Purkinje para o caso 13.

0 100 200 300 400 500 600 700 800 900 1000-200

-150

-100

-50

0

50

Tempo (ms)

I Na

(mA

)

Exata

Estimada

99% Confiança

Figura 5.80 - Estimativa de NaI de fibras de Purkinje para o caso 13.

119

0 100 200 300 400 500 600 700 800 900 100010

20

30

40

50

60

70

Tempo (ms)

I k (m

A)

Exata

Estimada

99% Confiança

Figura 5.81 - Estimativa de kI de fibras de Purkinje para o caso 13.

Pode-se observar nas Figuras (5.74-5.81) que para o caso 13, as estimativas das

variáveis de estado estão em excelente concordância com os valores exatos. Pode-se

perceber também que os intervalos de confiança estimados são pequenos. Assim,

sendo as estimativas tem pequenas incertezas. Assim como aconteceu com os

resultados obtidos para o caso 9, onde o modelo de HODGKIN-HUXLEY é aplicado

em axônios, isto se deve às pequenas incertezas associadas às medidas e ao modelo

de evolução neste caso.



120

0 100 200 300 400 500 600 700 800 900 1000-100

-80

-60

-40

-20

0

20

Tempo (ms)

Vm

(m

V)

Exata

Medidas

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação m

Exata

Estimada

99% Confiança


121

0 100 200 300 400 500 600 700 800 900 10000.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


122

0 100 200 300 400 500 600 700 800 900 1000-4

-2

0

2

4

6

8

10

12

14

16

Tempo (ms)

GN

a (m

S)

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 10000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Tempo (ms)

Gk

(mS

)

Exata

Estimada

99% Confiança


123

0 100 200 300 400 500 600 700 800 900 1000-200

-150

-100

-50

0

50

Tempo (ms)

I Na

(mA

)

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 100010

20

30

40

50

60

70

Tempo (ms)

I k (m

A)

Exata

Estimada

99% Confiança


124

Pode-se observar nas Figuras (5.82-5.89) que para o caso 14, as estimativas de

variáveis de estado ainda se encontram em boa concordância com os valores exatos.

Pode-se perceber também que o intervalo de confiança das estimativas é maior do

que aquele para o caso 13.



0 100 200 300 400 500 600 700 800 900 1000-100

-80

-60

-40

-20

0

20

Tempo (ms)

Vm

(m

V)

Exata

Medidas

Estimada

99% Confiança


125

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação m

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 10000.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança


126

0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 1000-4

-2

0

2

4

6

8

10

12

14

16

Tempo (ms)

GN

a (m

S)

Exata

Estimada

99% Confiança


127

0 100 200 300 400 500 600 700 800 900 10000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Tempo (ms)

Gk

(mS

)

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 1000-200

-150

-100

-50

0

50

Tempo (ms)

I Na

(mA

)

Exata

Estimada

99% Confiança


128

0 100 200 300 400 500 600 700 800 900 100010

20

30

40

50

60

70

Tempo (ms)

I k (m

A)

Exata

Estimada

99% Confiança



variáveis de estado estão em excelente concordância com os valores exatos, apesar

das maiores incertezas das medições.

A diferença entre os casos 13 e 15 é apenas o desvio padrão para as medidas,

lembrando que ,0.01med m máxV e ,0.05med m máxV , respectivamente, para os casos

13 e 15.



129

0 100 200 300 400 500 600 700 800 900 1000-100

-80

-60

-40

-20

0

20

Tempo (ms)

Vm

(m

V)

Exata

Medidas

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Ativação m

Exata

Estimada

99% Confiança


130

0 100 200 300 400 500 600 700 800 900 10000.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Tempo (ms)

Ativação n

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 10000

0.2

0.4

0.6

0.8

1

Tempo (ms)

Inativação h

Exata

Estimada

99% Confiança


131

0 100 200 300 400 500 600 700 800 900 1000-4

-2

0

2

4

6

8

10

12

14

16

Tempo (ms)

GN

a (m

S)

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 10000.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Tempo (ms)

Gk

(mS

)

Exata

Estimada

99% Confiança


132

0 100 200 300 400 500 600 700 800 900 1000-200

-150

-100

-50

0

50

Tempo (ms)

I Na

(mA

)

Exata

Estimada

99% Confiança


0 100 200 300 400 500 600 700 800 900 100010

20

30

40

50

60

70

Tempo (ms)

I k (m

A)

Exata

Estimada

99% Confiança


133


variáveis de estado estão em boa concordância com os valores exatos. Porém, para

este caso os intervalos de confiança são os maiores dentre os analisados, e devido as

grandes incertezas nas medições e nos modelos de evolução de estado.

Apresentam-se na Tabela 5.9 os erros RMS para as variáveis de estado

estimadas, nos casos 13 a 16.

Tabela 5.9 - Valores de erro RMS para as variáveis de estado para os casos 13, 14,

15 e 16.

Variável Caso 13 Caso 14 Caso 15 Caso 16

NaG S 0.035 (0.020) 0.041 (0.021) 0.101 (0.011) 0.403 (0.132)

kG S 0.006 (0.001) 0.015 (0.002) 0.013(0.001) 0.033 (0.003)

NaI mA 0.875 (0.250) 0.929 (0.264) 1.327 (0.118) 5.082 (1.596)

kI mA 0.228 (0.009) 0.639 (0.143) 0.460 (0.012) 1.208 (0.055)

mV mV 0.254 (0.003) 0.299 (0.008) 0.740 (0.010) 1.370 (0.024)

h 0.004 (0.001) 0.009 (0.002) 0.008 (0.001) 0.018 (0.001)

m 0.002 (0.001) 0.005 (0.001) 0.009 (0.001) 0.010 (0.001)

n 0.005 (0.0020) 0.024 (0.007) 0.006 (0.001) 0.007 (0.010)

Tal como para os casos envolvendo o axônio (ver Tabela 5.8), de um modo geral

as estimativas são mais afetadas pelas incertezas nas medidas do que pelas

incertezas no modelo de evolução de estado. De qualquer forma, tanto para o axônio

quanto para as fibras de Purkinje, os erros RMS aumentam com maiores incertezas

nas medidas ou maiores incertezas no modelo de evolução de estado. Por outro lado,

a concordância entre valores estimados e valores exatos é muito boa para todos os

casos analisados.

134

CAPÍTULO 6 - CONCLUSÕES E SUGESTÕES DE

TRABALHOS FUTUROS

Neste trabalho foi explorada a técnica de Monte Carlo via Cadeia de Markov para

estimar parâmetros, assim como os filtros de partículas SIR e ASIR para estimar as

variáveis de estado do modelo de HODGKIN-HUXLEY.

A análise dos coeficientes de sensibilidade revelou não ser possível estimar

simultaneamente todos os parâmetros do modelo de HODGKIN-HUXLEY, por terem

os coeficientes de sensibilidade linearmente dependentes. Usando o método de Monte

Calo via Cadeia de Markov, foi possível estimar o parâmetro que representa a

capacitância elétrica das células no modelo de HODGKIN-HUXLEY quando o modelo

foi aplicado em axônios.

A estimativa de variáveis de estado do modelo de HODGKIN-HUXLEY foi

realizada utilizando os filtros de partículas SIR e ASIR. Comparativamente ao filtro

SIR, obteve-se com o filtro ASIR melhores resultados para a estimativa das variáveis

de estado do modelo de HODGKIN-HUXLEY, aplicado tanto em axônios quanto em

fibras de Purkinje, com um número reduzido de partículas e com menores tempos

computacionais.

Espera-se que a estimativa de variáveis de estado do modelo de HODGKIN-

HUXLEY realizada neste trabalho possa auxiliar na identificação de fluxos iônicos

descontrolados e, assim, ajudar em prevenções e tratamentos de doenças como

Epilepsias, mal de Alzheimer, mal de Parkinson e arritmias cardíacas.

Como sugestões para trabalhos futuros temos o uso de modelos mais refinados

do que o de HODGKIN-HUXLEY, que representem a transmissão do potencial de

ação através das sinapses, a influência de neurotransmissores e a dependência

135

espacial do potencial de ação, além da estimativa do erro do modelo utilizado. A

sugestão em relação às estimativas de parâmetros e variáveis de estado é a utilização

de um algoritmo capaz de estimar simultaneamente parâmetros e variáveis de estado

e não de forma separada como foi realizado neste trabalho.

136

REFERÊNCIAS BIBLIOGRÁFICAS

AIRES, M.M., Fisiologia. 2 ed. Rio de Janeiro, Guanabara Koogan S.A, 1999.

ARULAMPALAM, M.S., MASKELL, S., GORDON, N., CLAPP, T., "A Tutorial on

Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking", IEE

Transactions on Signal Processing, vol. 50, pp.174-188, 2002

BASSANI, J.W.M, 1979, A Excitabilidade Celular: uma Formalização Matemática

Baseada no Modelo de Hodgkin-Huxley e sua Extensão às Células Cardíacas , Tese

de Mestrado, UNICAMP, CAMPINAS, SP, Brasil.

BATISTA,C.S, Solução de Equações de Balanço de Massa Populacional Usando a

Técnica da Transformada de Laplace e Filtro de Partículas, Tese de Doutorado,

PRODERNA, Belém, Pará, Brasil, 2011.

BEAR, M.F., CONNORS, B.W., PARADISO, M.A., Neurociências Desvendando o

Sistema Nervoso, 2 ed. São Paulo, Artmed S.A, 2006.

BECK, J.V., "Sequential Methods in Parameter Estimation", 3rd International

Conference on Inverse Problems in Enginnering, Tutorial Session, Port Ludlow, WA,

1999.

BECK, J.V., BLACKWELL, B., CLAIR, C.R., Inverse Heat Conduction: ILL-Posed

Problems, Wiley Interscience, New York, 1985.

137

BECK, J.V., ARNOLD, K.J., "Parameter Estimation in Engineering and Science",

Wiley-Interscience, New York, 1977.

BEI, C., CAI-WEN, M., ZHEN-TAO, L., " Improved Particle Filter Based on Fine

Resampling Algorithm", The Journal of China Universities of Posts and

Telecommunications, vol. 19, pp. 100-106, 2012.

BIANCHI, A.G.C., 2003, Caracterização, Modelagem e Simulação Matemático-

Computacional da Dinâmica do Crescimento e Conexões de Células Neurais, Tese,

USP, SÃO CARLOS, SP, Brasil.

BUHRY, L., PACE, M., SAÏGHI, S. " Global Parameter Estimation of an Hodgkin-

Huxley Formalism Using Membrane Voltage Recording: Application to Neuro-Mimetic

Analog Integrates Circuits", Neurocomputing, vol. 81, pp. 75-85, 2012.

CAMPOS, F.O., 2008, Modelagem Computacional da Eletrodifusão Cardíaca: O

Desenvolvimento de um Novo Modelo para Células de Camundongos e a Avaliação de

Novos Esquemas Numéricos, Dissertação, UFJF, JUIZ DE FORA, MG, Brasil.

CLEMENTI, L.A., VEGA, J.R., GUGLIOTTA, L.M., ORLANDE, H.R.B., "A Bayesian

Inversion Method For Estimating The Particle Size Distribution of Latexes from

Multiangle Dynamic Light Scattering Measurement", Chemometrics and Intelligent

Laboratory Systems, vol.107, pp. 165-173, 2011.

COLAÇO, M.J., ORLANDE, H.R.B, SILVA, W.B., DULIKRAVICH,G.S, “Aplication of

Bayesian Filter to Estimation Unknow Heat Fluxes in a Natural Convection Problem”,

International Design Engineering Technical Conferences & Computers and

Information in Engineering Conference, 2011.

138

CRUZ, E.F.,2006, Modelos Markovianos para Canais Iônicos em Membranas

Celulares. Dissertação, IFUSP/USP, SÃO PAULO, SP, Brasil.

DOUCET, A., FREITAS, N., GORDON, N., Sequential Monte Carlo Methods in

Practice, Springer, New York, 2001.

EHLERS, R.S., Introdução à Inferência Bayesiana, Departamento de Estatística,

Universidade Federal do Paraná, 2007.

FONSECA, A.V.S.,2009, Estimulação Multidirecional de Células Cardíacas:

Instrumentação e Experimentação . Dissertação, UNICAMP, SÃO PAULO, SP, Brasil.

FRÖHLIC, F., JEZERNIK, S., " Feedback Control of Hodgkin-Huxley Nerve Cell

Dynamics", Control enginnering pratice, vol.13, pp. 1195-1206, 2005.

GUYTON, A.C., HALL, J.E., Tratado de Fisiologia Médica. 11 ed. Rio de Janeiro,

Elsevier Editora Ltda, 2006.

HÄRTER, F.P., CAMPOS VELHO, H.F., REMPEL, E.L., CHIAN, A.C.L., "Neural

networks in Auroral Data Assimilation" Journal of Atmospheric And Solar-Terrestrial

Physics vol 70, pp.1243-1250, 2008.

HASTINGS, W.K., "Monte Carlo Sampling Methods Using Markov Chains and Their

Applications" Biometrika, vol 57, pp.97-109, 1970.

HODGKIN, A.L., HUXLEY, A.F. "A Quantitative Description of Membrane Current and

It's Application to Conduction and Excitation in Nerve", Journal of Phisiology, vol.

117, pp. 500-544, Mar. 1952.

139

JACKSON, M. B., Molecular and Cellular Biophysics. 1 ed. New York, Cambridge

University Press, 2006.

JARBAS, O., WACHTER, P.H., AZAMBUJA, A.A., BIOFÍSICA para ciências

biomédicas. 3 ed. Porto Alegre, edipurcs, 2008.

KAIPIO, J.P., SOMERSALO, E., " Estimating Anomalies From Indirect Observations",

Journal of Computational Physics, vol.181, pp-398-406, 2002.

KAIPIO, J.P., SOMERSALO, E., " Computational and Statistical Methods for Inverse

Problems", Springer, Berlin, 2004.

LI, T., SATTAR, T.P., SUN, S., " Deterministic Resampling: Unbiased Sampling to

Avoid Sample Impoverishment in Particle Filters", Signal Processing, vol. 92,

pp.1637-1645, 2012.

LIMA, M.F., 2011, "Processos Estocásticos não-Markovianos em Difusão Anômala"

Tese de doutorado, IF/UFAL, Alagoas, AL, Brasil.

LINDLEY, M.A., THOMAS, M.P., HOROWITZ, J.M., " A Simulation of the Effects of

Temperature on Hippocampal Neurons", computer methods and programs in

biomedicine, v. 25, pp. 3-12, 1987 .

MALTA, C., PINTO, R.D., LENZI, H.L., ROMANHA, W.S., MACHADO, M.P."

Modelagem Matemática do Comportamento Elétrico de Neurônios e Algumas

Aplicações", Sociedade Brasileira de Matemática Aplicada e Computacional, São

Carlos, vol. 13, Sp, Brasil. 2011.

140

MAYBECK, P., " Stochastic Models, Estimation and Control", Academic Press, New

York, 1979.

McCORMICK, N.J., "Inverse Problems: Methods and Applications", Tendências em

Matemática e Computacional, vol. 2, pp. 1-12. 2001.

MOTA, C.A.A., "Estimativa Simultânea do Fluxo de Calor e de Propriedades

Termofísicas de Materiais em Altas Temperaturas", Tese de Doutorado, UFRJ, Rio de

Janeiro, Rj, Brasil, 2007.

NAVEIRA-COTTA, C.P., "Problemas Inversos de Condução de Calor em Meios

Heterogêneos: Análise Teórico-Experimental Via Transformação Integral, Inferência

Bayesiana e Termografia por Infra Vermelho", Tese de Doutorado, UFRJ, Rio de

Janeiro, Rj, Brasil, 2009.

NAVEIRA-COTTA, C.P., COTTA.,R.M., ORLANDE, H.R.B.,"Inverse Analysis With

Integral Transformed Temperature Fields: Identification of Thermophysical Properties

in Heterogeneous Media", International Journal of Heat and Mass Transfer, vol.54,

pp.1506-1519, 2011.

NOBLE, D. A. “Modification of Hodgkin- Huxley Equations Applicable to Purkinje Fibre

Action and Pace – Maker Potentials”, Journal of Physiology, vol. 160, pp. 317 – 352.

(1962).

ORLANDE, H.R.B., "Inverse Problems in Heat Transfer: New Trends on Solution

Methdologies and Applications",14 th International Heat Transfer Conference.

141

ORLANDE, H.R.B., COLAÇO, M.J., DULIKRAVICH, G.S., "Approximation of the

Likelihood Function in the Bayesian Technique for the Solution of Inverse Problem",

Inverse Problems, 2008.

ORLANDE, H.R.B., FUDYM, O., MAILLET, D., COTTA, R.M., "Thermal Measurements

and Inverse Techniques", Taylor & Francis, 2011a.

ORLANDE, H.R.B., DULIKRAVICH, G.S., COLAÇO, M.J., "Application of Bayesian

Filters to Heat Conduction Problems", International Conference on Engineering

Optimization, 2008.

ORLANDE, H.R.B., COLAÇO, M.J., COTTA, C.P.N., GUIMARÃES, G., BORGES,

V.L.B., " Problemas Inversos em Transferência de Calor", sociedade brasileira de

matemática aplicada e computacional, São Carlos, vol. 51, Sp, Brasil. 2011b.

ÖZISIK, M.N., ORLANDE, H.R.B., " Inverse Heat Transfer", New York : Taylor &

Francis, 2000.

PARTHASARATHY, S., BALAJI, C., " Estimation of Parameters in Multi-Mode Heat

Transfer Problems Using Inference - Effect of Noise and a Priori", International

Journal of Heat and Mass Transfer, vol. 51, pp.2313-2334, 2008.

PICKARD, W.F., " On the Propagation of the Nervous Impulse Down Medullated and

Unmedullated Fibres", Journal of Theoretical Biology, vol.11, pp. 30-45, 1966.

RAUCH, G., GAMBLE, F. "Computerized Analysis of Ionic Channel Multi-State

Condutance in Planar Lipid Bilayers", Computer Methods And Program In

Biomedicine, vol.30, pp.297-303. 1989.

142

REUTER, H. “The Dependence of Slow in Purkinje on the Extracellular Calcium-

Concentration”, Journal of Physiology, vol. 192, pp. 479 – 492., Abr. 1967.

RISTIC, B., ARULAMPALAM, S., GORDON, N., Beyond the Kalman Filter, Artech

House, Boston, 2004.

SANTOS, D.O.C., 2009, Padrões de Disparo e Sincronização em um Modelo para

Atividade Epileptiformes Não-Sinápticas, Tese, UFMG, MINAS GERAIS, MG, Brasil.

SILVA, W. B. ORLANDE, H. R. B., COLAÇO, M. J., FUDYM, O.,"Application of

Bayesian Filters to a One-Dimensional Solidification Problem", 21nd Brazilian

Congress Mechanical Engineering , Brazil. 2011.

SORIANO, D.C., RICARDO, R.A., BASSANI., J.W. "Resolvendo o Modelo de Hodgkin-

Huxley: uma proposta didática", XX Congresso brasileiro de engenharia biomédica,

pp. 29-32, 2006.

SOUZA, G., 2005, SPINS - Um Simulador Neural Para Visualização De Aspectos De

Aprendizado Utilizando Neurônios Spiking. Dissertação, UFSC, SANTA CATARINA,

SC, Brasil.

STEIN, R.B., PEARSON, K.G.," Predict Amplitude and Form of Action Potentials

Recorded From Unmyelinated Nerve Fibres", Journal of Theoretical Biology, vol.32,

pp. 539-558, 1971.

143

VICTORRI, B., VINET, A., ROBERG, F.A., DROUHARD, J.P., "Numerical Integration in

the Reconstruction of Cardiac Action Potentials Using Hodgkin-Huxley Type Models",

Computers And Biomedical Research, vol. 18, pp. 10-23, 1985.

WILLMS, A.R., BARO, D.J., WARRICK, R.M.H., GUCKENHEIMER, J., " An Improved

Parameter Estimation Method For Hodgkin-Huxley Models ", Journal of

Computational Neuroscience, vol. 6, pp. 145-168, 1999.

Documents

ESTIMATIVA DE PARÂMETROS E VARIÁVEIS DE …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/25/pemufrj2012msc... · Figura 2.2 - Abertura e fechamento dos canais iônicos (GUYTON