COPPE/UFRJ
DINÂMICA NÃO-LINEAR E CAOS EM RITMOS CARDÍACOS
Sandra Regina Freitas da Silva Morgado de Gois
Tese de Doutorado 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 Doutor em Engenharia
Mecânica.
Orientador: Marcelo Amorim Savi
Rio de Janeiro
Outubro de 2010
DINÂMICA NÃO-LINEAR E CAOS EM RITMOS CARDÍACOS
Sandra Regina Freitas da Silva Morgado de Gois
TESE 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 DOUTOR EM CIÊNCIAS EM
ENGENHARIA MECÂNICA.
Examinada por:
________________________________________________
Prof. Marcelo Amorim Savi, D.Sc.
________________________________________________
Prof. Luciano Luporini Menegaldo, D.Sc.
________________________________________________
Profa. Djenane Cordeiro Pamplona D.Sc.
________________________________________________
Prof. Alberto Paiva, D.Sc.
________________________________________________
Prof. Daniel Alves Castello, D.Sc.
RIO DE JANEIRO, RJ - BRASIL
OUTUBRO DE 2010
Gois, Sandra Regina Freitas da Silva Morgado de
Dinâmica Não- linear e Caos em Ritmos Cardíacos/
Sandra Regina Freitas da Silva Morgado de Gois. – Rio de
Janeiro: UFRJ/COPPE, 2010.
XXII, 88 p.: il.; 29,7 cm.
Orientador: Marcelo Amorim Savi
Tese (doutorado) – UFRJ/ COPPE/ Programa de
Engenharia Mecânica, 2010.
Referencias Bibliográficas: p.75-81.
1. Dinâmica não-linear. 2. Ritmos Cardíacos. 3. EDD.
4. Biomêcanica I. Savi, Marcelo Amorim. II. Universidade
Federal do Rio de Janeiro, COPPE, Programa de
Engenharia Mecânica. III. Titulo.
i
Ao meu marido Jorge Audrin, sem sua ajuda nada seria possível.
À minha filha Maria Júlia, você que chegou no meio deste trabalho e deixou minha
vida mais iluminada.
ii
O meu agradecimento vai para aqueles que acreditaram que eu era capaz quando
nem eu mais acreditava...
O meu agradecimento vai para minha princesa, para ela que com poucos dias de vida
soube me entender, que foi minha amiguinha, que dormia quando eu precisava
estudar, e ainda me presenteava todos os dias com olhinhos de “te amo mamâe”.
Filha depois que você nasceu não me lembro de ter passado um dia sequer sem ter
sorrido...Obrigada, mamãe te ama!
O meu agradecimento ao meu anjo, protetor, amigo e incentivador, à meu grande
amor, ao meu marido Audrin. Sem palavras para expressar meu agradecimento e meu
amor por você...conseguimos!
O meu agradecimento vai para aqueles que mesmo longe se fizeram presentes
sempre me incentivando, me motivando e ainda compreenderam minhas poucas
visitas no decorrer deste tempo: aos meus pais, irmãos, cunhadas e sobrinhas.
O meu agradecimento especial vai para meu orientador, Marcelo Savi, responsável
direto pela missão que agora se cumpre. As indicações, as dicas, as correções, as
discussões, os conselhos e até mesmo alguns entreveros, enfim, tudo isto compôs
uma somatória fundamental não só para a construção do pensamento que se traduz
nas páginas deste texto, mas como para a maturidade de toda uma vida a seguir:
antes de tudo, este momento se dedica a este grande ser humano com carinho.
Aos meus amigos que sempre estiveram comigo mesmo quando eu não estava com
eles, que compreenderam a minha ausência em casamentos, nascimentos,
aniversários e no corriqueiro dia-a-dia...
Aos parentes que sentiram meu "sumiço" e ouviram muitos "não". "Não posso ir". "Não
vai dar"." Vou faltar". "Dessa vez eu não vou" "Avisa aí que não deu". Bom sentir
vocês todos sempre na torcida e acreditando.
Por fim, mas não por último, Àquele que ainda não descobri ao certo, mas sei que
inquieta meu espírito, instiga o coração e me leva à sempre buscar: agradeço a Deus
por todos os momentos maravilhosos que tenho tido em minha vida.
iii
Resumo da Tese apresentada à COPPE/UFRJ como parte dos requisitos necessários
para a obtenção do grau de Doutor em Ciências (D.Sc.)
DINÂMICA NÃO- LINEAR E CAOS EM RITMOS CARDÍACOS
Sandra Regina Freitas da Silva Morgado de Gois
Outubro/2010
Orientador: Marcelo Amorim Savi
Programa: Engenharia Mecânica
Fenômenos rítmicos representam uma das mais fascinantes manifestações
do comportamento dinâmico de sistemas biológicos. A compreensão dos mecanismos
reponsáveis pelos ritmos biológicos é essencial para entender a dinâmica da vida.
Ritmos naturais podem ser periódicos ou irregulares no tempo e no espaço. Este
trabalho avalia os ritmos cardíacos a partir de duas perspectivas: modelagem
matemática a partir de modelos de análogos físicos; e análise de séries temporais. A
modelagem matemática considera o acoplamento de três osciladores de Van der Pol
modificados com defasagem no tempo. Assim, a dinâmica cardíaca é representada
por um sistema de equações diferencias com defasagem. Simulações numéricas são
conduzidas apresentando uma boa concordância com o comportamento geral dos
ritmos cardíacos. De uma maneira geral, repreduzem-se ritmos normais e patológicos
representados pelo eletrocardiograma (ECG) a partir da variação dos parâmetros do
modelo. Depois da modelagem e simulação numérica, apresenta-se uma análise a
partir de séries temporais. A idéia é efetuar a reconstrução do espaço de estado e
verificar a caoticidade das séries a partir dos expoentes de Lyapunov. Esta análise é
realizada para sinais obtidos a partir do modelo matemático e depois, para sinais
experimentais. A contaminação por ruído torna a análise difícil impondo o uso de
filtros. O trabalho faz uma breve discussão sobre os filtros e a sua utilização na análise
de ECGs.
iv
Abstract of Thesis presented to COPPE/UFRJ as a partial fulfillment of the
requirements for the degree of Doctor of Science (D.Sc.)
NONLINEAR DYNAMICS AND CHAOS IN CARDIAC RHYTHMS
Sandra Regina Freitas da Silva Morgado de Gois
Octubre/2010
Advisor: Marcelo Amorin Savi
Department: Mechanical Engineering
Rhythmic phenomena represent one of the most striking manifestations of
dynamic behavior in biological systems. Understanding the mechanisms responsible
for biological rhythms is essential for the comprehension of the dynamics of life. Natural
rhythms could be either periodic or irregular over time and space. This work analyzes
cardiac rhythms from two different perspectives: mathematical modeling based on
physically similar systems; and time series analysis. In the mathematical modeling the
system is composed of three coupled Van der Pol oscillators with time delay between
them. Thus, the heart dynamics is represented by a system of differential-delay
equations. Numerical simulations are carried out and good agreement of the results to
the general heart rhythm is achieved. By variation of parameter values in the system,
normal and pathological rhythms represented by the electrocardiogram (ECG) are
reproduced. Later, a time series analysis is presented, where a state-space
reconstruction is carried out in order to verify the caoticity in the system, by means of
Lyapunov exponents. First, this analysis is conducted on signals obtained from the
mathematical model and, later, on experimental signals. Degradation, caused by noise,
makes the analysis much more difficult, what makes necessary the use of filters. A
brief discussion on filters and their application in ECG analysis is done.
v
SUMÁRIO
SUMÁRIO ...................................................................................................................................................... vi
1.
INTRODUÇÃO ....................................................................................................................................... 1
1.1.
O RGANIZAÇÃO DO T RABALHO ............................................................................................... 6
2. O CORAÇÃO .............................................................................................................................................. 8
2.1. ATIVIDADE E LÉTRICA C ARDÍACA ................................................................................................ 9
2.2. POTENCIAL DE ATIVAÇÃO .........................................................................................................11
2.3. ELETROCARDIOGRAMA (ECG) .................................................................................................14
3. MODELAGEM MATEMÁTICA ................................................................................................................20
3.1. EQUAÇÕES DIFERENCIAIS COM DEFASAGEM .............................................................................24
3.2. MÉTODO DE SOLUÇÃO : MÚLTIPLOS PASSOS .............................................................................26
3.3. SIMULAÇÃO NUMÉRICA ............................................................................................................27
4. ANÁLISE DE SÉRIES TEMPORAIS .........................................................................................................40
4.1. MÉTODO DA I NFORMAÇÃO MÚTUA MÉDIA .................................................................................42
4.2. MÉTODO DOS F ALSOS VIZINHOS ..............................................................................................44
4.3. EXPOENTES DE L YAPUNOV ......................................................................................................46
4.4. ANÁLISE DE ECG S G ERADOS A PARTIR DO MODELO MATEMÁTICO ............................................50
5. ANÁLISE DE ECGS EXPERIMENTAIS CONSIDERANDO PROCEDIMENTOS DE FILTRAGEM........61
5.1. REDUÇÃO L INEAR DE RUÍDO ....................................................................................................61
5.2. REDUÇÃO N ÃO -LINEAR DE RUÍDO ............................................................................................64
5.3. ANÁLISE DE ECG S UTILIZANDO REDUÇÃO NÃO -LINEAR DE RUÍDO ............................................66
6. CONCLUSÃO ............................................................................................................................................72
6.1. T RABALHOS F UTUROS .............................................................................................................73
7. REFERÊNCIAS ..........................................................................................................................................75
vi
1.
INTRODUÇÃO
Os fenômenos naturais apresentam ritmos que podem ser tanto periódicos quanto
irregulares em relação ao tempo e ao espaço. A dinâmica regular ou irregular de
sistemas biomédicos pode estar associada a funcionamentos fisiológicos normais ou
patológicos. Uma dinâmica extremamente regular pode representar doenças, como
por exemplo, certas anomalias respiratórias, doenças sanguíneas cíclicas, epilepsia,
espasmos neurológicos e tremores. Por outro lado, há fenômenos onde uma dinâmica
regular reflete comportamento saudável, como os ciclos do sono e menstrual. Ritmos
irregulares também podem refletir doenças tais como arritmias cardíacas e diferentes
desordens neurológicas (Savi, 2006, 2005; Glass, 2001; Ferriere & Fox, 1995).
Mudanças rítmicas da pressão sangüínea, da freqüência cardíaca e de outras
medidas cardiovasculares indicam a importância dos aspectos dinâmicos na
compreensão desses ritmos. Assim, para se obter uma compreensão dos sistemas
cardíacos e seus ritmos, faz-se necessário um estudo da dinâmica desses sistemas
biológicos.
Por desempenhar um papel fundamental na fisiologia dos seres vivos e pelo fato
de apresentar um comportamento rítmico complexo, o coração é um órgão
amplamente estudado por meio de modelos que abordam sua dinâmica sob os pontos
de vista quantitativos ou qualitativos.
Muitos estudos apontam para o fato de que certas arritmias cardíacas estão
associadas a uma resposta caótica (Witkowski et al., 1995; Radhakrishna et al., 2000).
Isso pode ser útil para sugerir diferentes estratégias terapêuticas, mudando as
abordagens clássicas. Dentre as arritmias cardíacas, destacam-se as fibrilações atrial
e ventricular, a bradicardia, e a taquicardia.
As arritmias clínicas que têm maior potencial para aplicações terapêuticas da teoria
do caos são as taquiarritmias aperiódicas, das quais fazem parte as fibrilações atrial e
1
ventricular. Garfinkel et al. (1992, 1995) discutem a aplicação de técnicas de controle
de caos a fim de evitar respostas cardíacas arrítmicas, abordagem que pode, por
exemplo, ser incorporada a marca-passos para evitar a fibrilação ventricular. Como
respostas caóticas podem ser controladas de modo eficiente através do método OGY
(método proposto por Ott, Gebrogi e Yorke, 1990) ou de suas variantes, isto pode
inspirar algumas abordagens interessantes a fim de estabilizar órbitas instáveis
associadas ao ritmo cardíaco normal (de Paula & Savi, 2007; Savi et al., 2006;
Pereira-Pinto et al., 2004).
As primeiras observações da dinâmica não-linear e caótica em sistemas
fisiológicos se deram no final da década de 1970 e começo da de 1980 por Mackey &
Glass (1977) enquanto Winfree (1980) aplicava métodos geométricos da dinâmica
não-linear às oscilações biológicas, especialmente aos ritmos circadianos e aos ritmos
cardíacos.
Diversos estudos realizados nos anos 1980 evidenciaram a natureza não-linear
dos processos cardíacos (Ritzenberg et al., 1984). Babloyantz & Destexhe (1988)
foram os primeiros a estudar o ECG em repouso, com ritmo sinusal – aquele ditado
pelo nódulo sino-atrial – normal, utilizando técnicas de dinâmica caótica. Eles
descobriram que a dimensão de correlação de um ECG varia entre 3,6 e 5,2 e,
portanto, sugere que a dinâmica ocorre em um atrator estranho determinístico de
baixa dimensão. Ravelli & Andolini (1990) também calcularam a dimensão de
correlação da atividade cardíaca humana, variando o ritmo sinusal normal para
fibrilação ventricular. Eles constataram que, conforme se evolui do ritmo sinusal para
uma completa fibrilação ventricular ocorre um aumento na complexidade da dinâmica
do sinal de ECG correspondente.
Goldsberger et al. (1990) propuseram que a variação da freqüência cardíaca
(VFC) era um processo dinâmico não-linear, de natureza fractal, cujo comportamento
complexo poderia ser explicado a partir de modelos caóticos. O controle do ritmo
cardíaco é perturbado por alterações na função do sistema nervoso autônomo em um
2
número importante de síndromes clínicas incluindo a morte súbita cardíaca, a falha
congestiva, a intoxicação com drogas, o sofrimento fetal, o envelhecimento fisiológico
entre outras. Essas condições estão associadas com a perda da complexidade da
dinâmica do batimento cardíaco. Tais alterações, que não são detectáveis usando
técnicas convencionais, podem ser quantificadas usando métodos derivados da
análise não-linear. Portanto, uma variedade de enfermidades que alteram a função
autônoma leva a uma perda da complexidade fisiológica e, portanto, a uma maior
regularidade. A dinâmica é tipicamente menos complexa em pacientes doentes do que
em sãos.
Deve-se enfatizar que a perda desse tipo de complexidade não-linear não
necessariamente vem acompanhada de uma diminuição da variação calculada por
métodos tradicionais. A medida da complexidade da VFC pode ser um novo indicador
de diagnóstico e prognóstico de maior exatidão e utilidade do que aqueles utilizados
atualmente (Glass, 1993).
Durante os últimos anos tem aumentado o número de estudos da VFC com o
enfoque não-linear e se tem comparado os resultados obtidos com aqueles obtidos
através de métodos lineares (Voss et al.,1994), observando que ambas as análises
proporcionam informações distintas. Também se tem proposto diversos métodos nãolineares para estimar a complexidade de séries temporais (Voss et al.,1994, 1995;
Lombardi et al.,1996; Abarbanel, 1996). Muitos trabalhos com enfoque em
metodologias de técnicas não-lineares tais como dimensão do atrator, a entropia, os
índices de informação mútua, os expoentes de Lyapunov, os mapas de Poincaré,
dinâmica simbólica, entre outras, têm sido particularmente úteis para caracterizar
diversas patologias cardíacas, observando-se que esses índices têm valores inferiores
em enfermos cardíacos, e que em uma mesma patologia, valores de dimensão
menores caracterizam maior gravidade da cardiopatia (Carvajal et al., 1996,1997;
Zebrowski et al., 1997, 1997; Signorini et al., 1994).
3
Há diversas formas de se avaliar indiretamente o funcionamento do coração.
Muitas delas a partir da medição não invasiva de alguma variável ligada ao processo
cardíaco. O eletrocardiograma (ECG) é uma dessas formas, onde são gravados os
impulsos elétricos relacionados ao funcionamento cardíaco sob a forma de ondas que
representam a corrente elétrica em diferentes áreas do coração. O sinal captado pode
ser utilizado para medir a freqüência e a regularidade dos batimentos cardíacos, assim
como o tamanho e posição das câmaras cardíacas.
A idéia de tratar o coração como um sistema de osciladores não-lineares
acoplados data de 1928 com o artigo clássico de Van der Pol & Van der Mark onde
descrevem o batimento cardíaco por meio de sistemas elétricos acoplados exibindo
osciladores de relaxação. Mais tarde, as equações de Van der Pol, originalmente
concebidas para descrever osciladores de relaxação em sistemas eletrônicos, foi
freqüentemente usada em modelos teóricos das funções cardíacas. Em 1952, Hodkin
& Huxley descrevem o comportamento do batimento cardíaco considerando diferentes
modelos. Recentemente, Grudzinski & Zebrowski (2004) propuseram uma variação do
oscilador de Van der Pol clássico a fim de capturar alguns aspectos importantes,
considerando o potencial de ação do coração. Santos et al. (2004) utilizam osciladores
acoplados a fim de descrever os aspectos gerais dos ritmos cardíacos, discutindo
diferentes termos de acoplamento. Campbell & Wang (1998) apresentam uma
discussão lidando com algumas características desses osciladores, destacando a
importância de se considerar os acoplamentos com defasagem.
Outra abordagem para a análise dos ritmos cardíacos é a partir de séries
temporais (Abarbanel, 1996). Takens (1981) afirmou que é possível recriar um
retrato equivalente da topologia do sistema multidimensional usando observações
de uma única variável do sistema. A idéia básica é que o efeito de todas as outras
variáveis não observáveis se reflete na variável observada. As regras que
governam o comportamento do sistema original podem ser reconstruídas a partir
da saída, isto é, vetores em um novo espaço, o espaço de reconstrução ou espaço
4
de imersão que pode ser construído a partir de valores em atraso da série temporal
observada.
Este trabalho propõe a análise dos ritmos cardíacos a partir de duas abordagens
distintas: modelagem matemática a partir de modelos de análogos físicos; e análise de
séries temporais.
A modelagem matemática usa o fato de que as características qualitativas do
potencial de atuação cardíaco apresentam grande semelhança com a resposta
dinâmica do oscilador de Van der Pol clássico (VdP). Desta forma, esse oscilador é
considerado como ponto inicial da modelagem. Propõe-se então, um modelo inédito
com três osciladores de VdP modificados conectados considerando acoplamentos
defasados a fim de reproduzir sinais de ECG. Portanto, a dinâmica cardíaca é
representada por um sistema de equações diferenciais com defasagem. Simulações
numéricas são realizadas reproduzindo ritmos normais e patológicos ao se tomar
variações apropriadas tanto nos termos de acoplamento quanto no forçamento. Com
isso, pretende-se chegar a um modelo que possa ser utilizado como ferramenta de
análise ou mesmo de predição do comportamento cardíaco, podendo assim verificar a
influência de diversos fatores sobre o ritmo cardíaco.
A análise de séries temporais, por sua vez, considera sinais de ECG simulados e
experimentais. Em linhas gerais, inicia-se a análise a partir de sinais gerados através
do modelo matemático proposto. Depois disso, passa-se a tratar sinais experimentais
onde, por meio de uma comparação desses sinais com os ECGs gerados pelo modelo
proposto, pretende-se fazer uma validação do modelo e mesmo obter uma melhor
compreensão do significado físico-biológico de cada um de seus termos.
O método de coordenadas defasadas é empregado para estabelecer o espaço de
estado reconstruído. Parâmetros de defasagem são calculados a partir dos métodos
da informação mútua média e do método dos falsos vizinhos. Os expoentes de
Lyapunov são utilizados para avaliar a caoticidade do sinal, sendo estimados através
5
do algoritmo de Kantz. Os resultados desses métodos são baseados no pacote
TISEAN 3.0 (Hegger et al., 1999).
Desta forma, através do uso de ferramentas para a análise de sistemas nãolineares, são estabelecidos parâmetros que permitam uma comparação qualitativa das
dinâmicas dos sistemas modelado e real. Uma dificuldade encontrada nesta
comparação é a forte presença de ruído em eletrocardiogramas reais, pois tal
perturbação causa variações nos parâmetros calculados para a caracterização da
dinâmica do sistema elétrico cardíaco. Portanto, antes de ser feita a caracterização do
ECG real é necessário submetê-lo a um processo de filtragem. Este trabalho
apresenta uma discussão sobre processos de filtragem visando a sua aplicação em
ECGs.
1.1.
ORGANIZAÇÃO DO TRABALHO
O Capítulo 1 apresenta uma exposição da motivação do trabalho bem como de
seus objetivos, sendo também apresentada uma revisão bibliográfica sobre o tema.
No Capítulo 2 são discutidos os aspectos gerais do comportamento do elemento
central deste tema – o coração – assim como do eletrocardiograma (ECG).
O Capítulo 3 propõe um modelo matemático da dinâmica do coração assim como
diversas simulações numéricas do modelo proposto. Utilizam-se ECGs experimentais
como referência, indicando o comportamento normal e diferentes patologias
identificadas a partir dos ECGs.
A análise por meio das séries temporais obtidas a partir do comportamento elétrico
do coração é apresentada nos Capítulos 4 e 5. No Capítulo 4 apresentam-se algumas
das ferramentas básicas para a análise de séries temporais. Além disso, desenvolvese a análise de séries obtidas a partir do modelo matemático. A análise de sinais
experimentais reais também é realizada, mostrando as suas principais dificuldades. A
6
partir dos resultados obtidos, apresenta-se o Capítulo 5 que trata da aplicação de
filtros para realizar a análise de séries temporais.
Finalmente o Capítulo 6 apresenta as conclusões do trabalho, sendo também
discutidas possibilidades de trabalhos futuros.
7
2. O CORAÇÃO
O coração é um órgão oco e musculoso, cuja função é bombear sangue para os
pulmões, órgãos e tecidos. Localiza-se anatomicamente um pouco à esquerda do
centro do tórax, no sentido ântero-posterior, sendo que o ápice do coração é
acentuadamente deslocado para a esquerda.
O coração consiste em duas partes distintas: direita e o esquerda. Uma bomba
propulsiona o sangue através dos pulmões para as trocas de oxigênio e dióxido de
carbono (a circulação pulmonar), e a outra bomba propulsiona o sangue para outros
tecidos do corpo (a circulação sistêmica). Cada uma dessas partes impele sangue
para os pulmões e para os tecidos. Cada parte do coração é dividida em duas, o que
faz com que o coração tenha quatro compartimentos: os átrios direito e esquerdo
(parte superior) e os ventrículos direito e esquerdo (parte inferior), mostrados na
Figura 1. O sangue retorna do corpo pelo sistema circulatório para o átrio direito e dali
vai para o ventrículo direito. Ele é então ejetado do ventrículo direito para os pulmões,
onde é oxigenado, seguindo dos pulmões para o átrio esquerdo, e dali para o
ventrículo esquerdo. Finalmente o sangue é bombeado através da válvula aórtica para
a artéria aorta e o sistema circulatório. O átrio funciona como uma bomba para o
ventrículo, enchendo o ventrículo de sangue para que seja bombeado para adiante
pelo sistema circulatório. Os ventrículos são os grandes responsáveis por impulsionar
o sangue para a circulação pulmonar e sistêmica.
O fluxo de sangue através do coração é unidirecional sendo conseguido por
distribuição adequada de válvulas. Embora o débito cardíaco seja intermitente, o fluxo
contínuo para os tecidos corporais (periferias) ocorre pela distensão da aorta e suas
ramificações durante a contração ventricular (sístole) e pela retração elástica das
paredes das grandes artérias, com a propulsão do sangue, para frente, durante o
relaxamento ventricular, a diástole (Malmivuo & Plonsey, 1995).
8
Figura 1 – Figura esquemática do coração.
As paredes do coração são compostas por músculo cardíaco, chamado miocárdio.
Este tipo de tecido muscular é estriado, assim como os músculos esqueléticos, porém
sua contração é involuntária.
A célula muscular do coração (miócito) tem uma ativação elétrica que acontece por
meio de um fluxo de entrada de íons de sódio através da membrana celular. Uma fase
estacionária segue a despolarização cardíaca e, depois disso, ocorre a repolarização
que é conseqüência do fluxo de saída de íons de potássio. A contração mecânica do
coração é associada à ativação elétrica das células musculares cardíacas.
2.1. ATIVIDADE ELÉTRICA CARDÍACA
Há dois séculos, Galvani e Volta demonstraram que fenômenos elétricos estavam
envolvidos nas contrações espontâneas do coração. Kölliker e Müller (1855)
descobriram que quando eles colocavam o nervo de uma preparação de músculo
9
esquelético inervado em contato com a superfície de um coração de rã batendo, o
músculo esquelético contraía-se a cada contração cardíaca. Assim os pesquisadores
concluíram que a excitação espontânea do coração havia gerado atividade elétrica
suficiente para excitar as fibras nervosas motoras e estimular o músculo esquelético.
Uma distinção importante entre tecido muscular cardíaco e tecido muscular
esquelético é que a ativação do músculo cardíaco pode se propagar de uma célula
para a outra em qualquer direção. Portanto, as frentes de onda de ativação são
bastante complexas. A única exceção é o limite entre átrios e ventrículos, por onde a
onda de ativação normalmente não pode passar a não ser através de um sistema
condutor especial, já que há uma barreira de tecido fibroso isolante.
Os eventos elétricos que normalmente ocorrem no coração iniciam a contração
cardíaca. Distúrbios na atividade elétrica podem levar a distúrbios sérios e, às vezes,
letais no ritmo cardíaco. Todos os ritmos cardíacos são governados por potenciais de
atuação que são controlados por células especializadas. O nódulo sinusal (sinoatrial
ou nódulo SA) é localizado no átrio direito, junto à veia cava superior (Figura 1). As
células do nódulo SA são auto-excitatórias, células marca-passo, que geram o
potencial de ação numa freqüência de 70 pulsos/minuto. A ativação se propaga a
partir do nódulo SA através do átrio, mas não pode se propagar diretamente através
da interface entre átrios e ventrículos. O nódulo atrioventricular (nódulo AV) é
localizado na fronteira entre os átrios e os ventrículos (Figura 1). Ele tem uma
freqüência intrínseca de cerca de 50 pulsos/minuto, no entanto, este nódulo segue
freqüências mais altas quando excitado. Em um coração normal, o nódulo AV é o
único caminho condutor dos átrios para os ventrículos. Portanto, sob condições
normais, ventrículos podem ser excitados apenas por impulsos que passem através
dele. A propagação do nódulo AV para os ventrículos é feita através de um sistema
especializado de condução. Junto ao nódulo, este sistema é composto por um feixe
comum, chamado de feixe de His. Mais adiante, ele se separa em dois ramos,
continuando ao longo de cada lado do septo, formando os ramos direito e esquerdo do
10
feixe. Ainda mais adiante, os feixes se ramificam, tornando-se as fibras de Purkinje,
que divergem para o lado interno das paredes ventriculares.
Como a freqüência intrínseca do nódulo SA é a mais alta, ela dita a freqüência de
ativação de todo o coração. Se a conexão do átrio ao nódulo AV falhar, este último
adota sua freqüência intrínseca. Se o sistema de condução do feixe de His falhar, os
ventrículos irão se contrair numa freqüência determinada por eles próprios.
2.2. POTENCIAL DE ATIVAÇÃO
A atividade elétrica do coração está intimamente relacionada aos ritmos celulares
regulados pela bomba de Na-K. A concentração dos íons de sódio (Na+) é cerca de 10
vezes mais alta fora da membrana celular do que dentro, enquanto a concentração de
íon de potássio (K+) é cerca de 30 vezes mais alta do lado de dentro se comparada ao
lado de fora. Estimulando a membrana, o potencial transmembrana sobe cerca de 20
mV e atinge o limite, isto é, quando a tensão elétrica muda de −70 mV para cerca de
−50 mV, a permeabilidade da membrana ao sódio e ao potássio iônicos muda.
A Figura 2 apresenta um esboço do fluxo iônico que gera o potencial de ação.
Primeiramente, a permeabilidade ao íon de sódio aumenta muito rapidamente,
permitindo aos íons de sódio fluir do lado de fora para dentro da célula, aumentando o
potencial positivo dentro, atingindo cerca de +20 mV. Em seguida, há um lento
crescimento da permeabilidade aos íons de potássio, que fluem de dentro para fora,
causando o retorno do potencial intracelular ao seu valor de repouso.
O pico máximo do valor de tensão elétrica na membrana durante a ativação é de
cerca de 100 mV; a duração do impulso é aproximadamente 300ms. Durante o
repouso, depois da ativação, a bomba de Na-K restaura as concentrações iônicas
dentro e fora da membrana aos seus valores originais.
11
Figura 2 – A bomba de Na-K e o potencial de ativação.
Uma vez que a ativação foi iniciada, a membrana não é mais sensível a novos
estímulos, não importa sua magnitude. Esta fase é chamada de período refratário
absoluto. O processo de ativação envolve fatores como correntes, potenciais,
condutividades, concentrações e fluxo de íons. O termo impulso de ação descreve
todo o processo. Medições bioelétricas focalizam a diferença de potencial através da
membrana. Portanto, a medida elétrica do impulso de ação é chamada de potencial de
ação que descreve o comportamento do potencial de membrana durante a ativação.
Diferentes grupos de células apresentam diferentes configurações de canais
iônicos, o que leva a diferentes velocidades do impulso de ação. Portanto, as células
são classificadas como as de potencial de ação de resposta lenta e de resposta rápida
que apresentam formas diferentes para o sinal do potencial de ação, representado na
Figura 3. O primeiro grupo – resposta lenta – inclui os músculos atriais e ventriculares,
12
o feixe de His e as fibras de Purkinje, enquanto o segundo - resposta rápida –
corresponde aos nódulos SA e AV.
Figura 3 – Forma da curva do Potencial de Ação: resposta lenta (lado esquerdo) e
resposta rápida (lado direito).
Os ritmos celulares produzem atividade elétrica que é relacionada às ondas de
despolarização e repolarização. O fluxo de sódio para dentro da célula é acionado pelo
crescimento grande e rápido da permeabilidade ao sódio. Esta corrente é confinada à
zona chamada de zona de despolarização. A corrente devido ao fluxo de saída é a
corrente do circuito local, que inicialmente despolariza o tecido em repouso e avança
na direção de propagação.
Quando uma célula é despolarizada, ela produz um campo elétrico que inicia o
fenômeno de despolarização em outra célula próxima a ela e então se despolariza.
Portanto, a despolarização pode ser compreendida como uma onda que se propaga
pelo tecido cardíaco. A natureza da onda de repolarização é diferente daquela da onda
de despolarização.
Diferentemente da despolarização, a repolarização não é um fenômeno que se
propaga, porém, ao se examinar a localização das células onde está ocorrendo
repolarização em instantes de tempo consecutivos, é possível aproximar a
repolarização por uma onda que se propaga pelo tecido.
13
2.3. ELETROCARDIOGRAMA (ECG)
Os fenômenos elétricos que se originam durante a atividade cardíaca podem ser
registrados por meio de microeletrodos denominados eletrocardiógrafos. Desse modo,
é possível registrar na superfície do tórax o potencial elétrico gerado pela atividade
elétrica do coração. A medida do comportamento elétrico extra-celular do tecido
muscular cardíaco é chamada eletrocardiograma.
Em 1888, Waller, realizando estudos experimentais em cães, variou a posição de
dois eletrodos na superfície cardíaca obtendo uma diferença de potencial que era
máxima quando eletrodos situavam-se na base e na ponta do coração. Esse órgão
pôde então ser comparado a uma pilha cujos pólos geradores se encontravam
naqueles pontos. Já em 1902, Einthoven, considerado “pai” do eletrocardiograma,
idealizou um método para captar energia elétrica na superfície corpórea empregando
um galvanômetro de corda, criando assim o sistema de derivações clássicas, que
trouxe grandes contribuições a uma área, até então, pouco desenvolvida.
A propagação da onda do sinal elétrico cardíaco através do corpo apresenta uma
forma muito complicada devido aos valores diferentes da condutibilidade dos tecidos
que compõem o corpo humano e sua distribuição espacial complexa.
Figura 4 - Sistema de ECG com 10 e 3 eletrodos respectivamente.
14
Adicionalmente, como os fenômenos da despolarização e da repolarização se
propagam através do sistema de condução do coração, o campo elétrico tem uma
resposta periódica. Muitos sistemas de medida do sinal cardíacos já foram propostos
para medir a diferença potencial entre pontos na superfície do corpo (Figura 4).
Conseqüentemente, dependendo da configuração das ligações, diferentes
medidas são obtidas, mas, geralmente, o sinal contém as seguintes ondas:
• Onda P: é a primeira onda registrada no ECG, representa a ativação do
átrio logo após a estimulação do nódulo sinus. Dura normalmente entre 60 e
90ms em adultos, e tem uma forma arredondada com amplitude máxima entre
0,25 e 0,30mV.
•
Intervalo PQ: É medido normalmente a partir da onda P até o começo
do complexo QRS e dura 90ms.
•
Complexo QRS: Corresponde à ativação ventricular e é medido do
começo da primeira onda (não importa que a onda seja Q ou R), à última onda
(seja ela R ou S). Em adultos normais, o complexo dura aproximadamente
80ms e apresenta uma forma pontiaguda, por causa das altas freqüências do
sinal. Sua forma varia muito, dependendo do sistema de leitura usado.
•
Intervalo ST: Vai do fim do complexo QRS ao começo da onda T e
corresponde à parte do processo de repolarização ventricular
•
Onda
T: Representa
a ativação
ventricular,
tem um formato
arredondado e uma amplitude em torno de 0.60mV. Nem sempre aparece no
ECG normal.
15
A diferença de potencial no tecido do coração pode ser interpretada como um
dipolo elétrico assim como a propagação do pulso pode ser interpretado como um
dipolo em movimento. O sinal medido é a projeção do dipolo ao longo da linha definida
pelas ligações.
Na Figura 5 pode-se acompanhar o deslocamento do dipolo móvel (seta amarela)
bem como sua trajetória (linha verde tracejada) e os correspondentes sinais captados
por um sistema de três eletrodos (linha preta). Assim seguindo a legenda:
(a) o coração sai do seu estágio original de polarização;
(b) surge o primeiro pulso no nódulo sinoatrial, que se propaga pelo átrio, gerando
assim o dipolo;
(c) em seguida, o átrio esta completamente estimulado;
(d) através do nódulo atrioventricular o pulso é transmitido para o ventrícul,
gerando novamente o dipolo;
(e) o pulso é transmitido pelos feixes de His, passando então por ambos os lados
dos feixes;
(f) a partir daí é transmitido através das fibras de Purkinje, que se espalham pelo
ventrículo;
(g) o pulso atravessa as paredes ventriculares;
(h) neste ponto, o ventrículo se encontra completamente estimulado e contraído,
sendo que o dipolo é novamente nulo;
(i) a repolarização ocorre do lado externo para o interno do ventrículo
(j) por fim, o coração retorna ao seu estado original de polarização, estando
pronto para recomeçar o ciclo.
16
Figura 5 – Formação do ECG e o dipolo.
17
Tabela 1 – Eventos Elétricos no Coração
Localização no coração
Evento
Emissão do
Nódulo AS
ECG
Velocidade de
(ms)
Terminologia
condução [m/s]
0
Pulso
Átrio Direito
Tempo
Despolarização
0,05
5
0,8-1,0
P
Átrio Esquerdo
Despolarização
Chegada do
Nódulo AV
85
50
Pulso
0,8-1,0
Intervalo
P-Q
Partida do
Nódulo AV
0,02-0,05
125
Pulso
Feixe de His
Ativado
130
1,0-1,5
Ramificação do feixe
Ativado
145
1,0-1,5
Fibras de Purkinje
Ativado
150
3,0-3,5
Despolarização
175
Despolarização
190
Despolarização
225
Septo
Endocárdio
Ventrículo
Esquerdo
Ventrículo
Esquerdo
0,3
QRS
0,8
Ventrículo
Direito
Despolarização
250
Epicárdio
Ventrículo
Esquerdo
Repolarização
400
Ventrículo
Direito
Endocárdio
Ventrículo
Esquerdo
Repolarização
Repolarização
18
T
600
0,5
Algumas características dos eventos elétricos cardíacos e de suas ondas estão
relacionadas na Tabela 1. A partir de tais características é possível indicar a forma de
cada onda e também o resultado de sua combinação em termos de amplitude do
dipolo, mostrados na Figura 6. O ECG é na verdade a projeção de tal sinal de acordo
com um ângulo variável (a orientação do dipolo), sendo chamada de derivação, cada
uma das projeções em cada diferente sistema de medição.
Figura 6 – Diferentes formas de ondas relacionadas a cada tipo de célula
especializada do coração
O sinal do eletrocardiograma é medido entre dois eletrodos, sendo os sistemas
mais comuns o de dez e o de três eletrodos, mostrados na Fig. 4. Os eletrodos são
posicionados sobre pontos pré-definidos do corpo, formando pares que definem linhas
sobre este mesmo corpo. As derivações do eletrocardiograma são então a projeções
do dipolo elétrico cardíaco sobre cada uma dessas linhas.
19
3. MODELAGEM MATEMÁTICA
A equação de Van der Pol (VdP) foi originalmente empregada na descrição de
osciladores de relaxação em circuitos eletrônicos, tendo sido freqüentemente usada
em modelos teóricos do ritmo cardíaco. Esse tipo de equação é muito útil na
modelagem fenomenológica de sistemas naturais, especialmente dos batimentos
cardíacos, já que ele apresenta características presentes em sistemas biológicos,
como ciclo limite, sincronização e caos. A forma geral desta equação é apresentada
abaixo:
(
)
&x& + a 1 − bx 2 x& + cx = Γ(t )
(1)
onde a, b e c são parâmetros do sistema e Γ(t) é um forçamento externo.
As características de um oscilador de VdP isolado apresentam grande semelhança
com as características do potencial de atuação cardíaco. Os dois tipos de resposta do
potencial de ação (rápido e lento) podem ser facilmente reproduzidos pelo oscilador de
VdP. Grudzinski & Zebrowski (2004) propuseram um oscilador modificado a fim de
simular importantes características fisiológicas do potencial de ação. Em geral, a nova
equação tem dois pontos fixos e um termo de dissipação assimétrico relativo à tensão
elétrica:
&x& + a ( x − w1 )( x − w2 ) x& −
x( x + d )( x + e)
= Γ(t )
ed
(2)
aqui a, d, e, w1 e w2 são parâmetros do sistema e Γ(t) é um forçamento externo.
Como já mencionado, o ritmo cardíaco normal é gerado primariamente pelo nódulo
SA, que é considerado o marca-passo normal, sendo o nódulo AV um segundo marca-
20
passo. Cada uma dessas instâncias apresenta um potencial de atuação que é
fundamental para a dinâmica do coração, mas não necessariamente a mais
expressiva ao compor o sinal do ECG. Cada ativação (despolarização seguida de
repolarização) corresponde a uma região diferente do coração e, por conseguinte, a
uma massa de tecido diferente, gerando correntes de diferentes magnitudes. Portanto,
a combinação de ondas de ativação oriundas de cada região do coração é
responsável pela forma do ECG, não tendo de corresponder exatamente ao esboço da
Figura 6. Alguns desses sinais podem ser preponderantes nessa composição, como
as ondas originadas nos átrios e ventrículos. Por outro lado, como essas regiões
seguem a ativação dos nódulos SA e AV, a assinatura no sinal de ECG é
representativa de sinais de marca-passo e é possível associar estes sinais ao átrio e
ao ventrículo, respectivamente.
A partir dessas hipóteses espera-se que osciladores acoplados possam
representar uma dinâmica geral do batimento cardíaco, onde cada um dos osciladores
representa o sinal de uma região diferente do coração. Normalmente, a dinâmica do
coração é modelada a partir de dois osciladores que representam os nódulos SA e AV.
No entanto, observa-se que esses dois osciladores não são suficientes para reproduzir
o sinal de ECG, pois o sinal do primeiro oscilador corresponde à ativação do nódulo
SA e do átrio, e o sinal do segundo oscilador corresponde apenas à despolarização
ventricular. Assim, é possível reproduzir a curva P, mas não o complexo QRS, pois
este intervalo corresponde principalmente à repolarização ventricular.
Esta observação motiva a inclusão de um terceiro oscilador que represente a
propagação do pulso através dos ventrículos, o que fisiologicamente representa o
complexo de His-Purkinje, composto pelo feixe de His e as fibras de Purkinje. A Figura
7 representa o modelo conceitual mostrando tanto os osciladores quanto os
acoplamentos entre eles. A fim de construir um modelo geral, considera-se que todos
os acoplamentos são assimétricos e bidirecionais. Além disso, excitações externas
são incorporadas ao sistema, considerando-se um termo periódico em cada oscilador,
21
Γi(t). Este modelo conceitual pode ser representado por um conjunto de equações
diferenciais como mostrado a seguir.
Γ1
Nódulo
SA
Γ2
Complexo
Nódulo
HP
AV
Γ3
Figura 7 – Modelo conceitual com três osciladores acoplados.
x&1 = x2
x& 2 = −aSA x2 ( x1 − wSA1 )( x1 − wSA2 ) − x1 ( x1 + d SA )( x1 + eSA ) + ρ SA sin(ω SAt ) +
+ k SA− AV ( x1 − x3 ) + k SA−HP ( x1 − x5 )
x&3 = x4
x& 4 = −a AV x4 ( x3 − wAV1 )( x3 − wAV2 ) − x3 ( x3 + d AV )( x3 + e AV ) + ρ AV sin(ω AV t ) +
+ k AV −SA ( x3 − x1 ) + k AV −HP ( x3 − x5 )
(3)
x&5 = x6
x&6 = −aHP x6 ( x5 − wHP1 )( x5 − wHP2 ) − x5 ( x5 + d HP )( x5 + eHP ) + ρ HP sin(ω HPt ) +
+ k HP−SA ( x5 − x1 ) + k HP− AV ( x5 − x3 )
Este conjunto de equações pode ser reescrito em uma forma compacta, como:
(4)
&y& = f(y, y& ) + Γ(t ) + Ky
22
y = [x1
x3
x5 ]
T
y& = [x 2
x4
x6 ]
T
Γ = [ρ SA sin(ω SA t ) ρ AV sin( ω AV t ) ρ HP sin(ω HP t )]
T
 k SA− AV + k SA− HP
K =  − k AV − SA
 − k HP − SA
− k SA− AV
(5)


− k AV − HP

k HP − SA + k HP − AV 
− k SA− HP
k AV − SA + k AV − HP
− k HP − AV
Note que as equações que governam o problema têm a forma geral f (y, y& )
(relacionada ao oscilador de VdP modificado), Γ(t) é um termo de forçamento e K
representa a matriz de acoplamento.
A defasagem na transmissão de sinais é inevitável e, já que até mesmo pequenas
defasagens podem alterar a dinâmica do sistema, é necessário entender como a
defasagem na condução muda o comportamento de osciladores acoplados. A inclusão
da defasagem em equações diferenciais pode causar mudanças drásticas e até
mesmo o surgimento do caos em um sistema que em outra situação apresentaria um
comportamento regular (Campbell & Wang, 1998). Deste modo, o modelo matemático
proposto é modificado a fim de incluir aspectos de defasagem temporal nos termos de
acoplamento. Portanto, as equações de governo são modificadas para:
x&1 = x2
x&2 = −aSA x2 ( x1 − wSA1 )( x1 − wSA2 ) − x1 ( x1 + d SA )( x1 + eSA ) + ρ SA sin(ωSAt ) +
+ kSA− AV ( x1 − xτ3 SA− AV ) + kSA− HP ( x1 − xτ5 SA−HP )
x&3 = x4
x&4 = −aAV x4 ( x3 − wAV1 )( x3 − wAV2 ) − x3 ( x3 + d AV )( x3 + eAV ) + ρ AV sin(ωAV t ) +
+ k AV − SA( x3 − x1τ AV −SA ) + k AV − HP ( x3 − xτ5 AV −HP )
x&5 = x6
x&6 = −aHP x6 ( x5 − wHP1 )( x5 − wHP2 ) − x5 ( x5 + d HP )( x5 + eHP ) + ρ HP sin(ωHPt ) +
+ kHP− SA( x5 − x1τ HP−SA ) + kHP− AV ( x5 − xτ3 HP− AV )
23
(6)
onde xiτ = xi (t − τ ) e τ representa a defasagem. Note que, na verdade, há diferentes
defasagens dependendo do tipo de conexão.
A idéia geral destes osciladores acoplados é a de que o sinal do ECG é formado
pela composição destes sinais:
(7)
X = ECG = α 0 + α1 x1 + α 3 x3 + α 5 x5
Analogamente, é possível definir:
d ( ECG)
X& =
= α1 x 2 + α 3 x 4 + α 5 x 6
dt
(8)
A solução de sistemas de equações diferenciais com a forma mostrada em Eq. 6
não se dá pela simples integração como em um sistema de equações diferenciais
ordinárias (EDOs), portanto, alguns métodos de solução são discutidos a seguir.
3.1. EQUAÇÕES DIFERENCIAIS COM DEFASAGEM
Equações diferenciais que possuem dependência de instantes defasados são
chamadas equações diferenciais com defasagem (delay-differential equations), cuja
forma geral é mostrada em Eq. 9,
, , (9)
onde yt = { y( τ ): τ ≤ t }. Tais equações podem se apresentar em diferentes formas,
contendo uma defasagem contínua ou discreta (Eq. 10), sendo neste último caso
24
chamadas de equações diferenciais com defasagem, EDD (difference-differential
equations).
, , τ , τ , . . ., τ com τ1 >...> τn ≥ 0 (10)
O sistema de equações apresentado em Eq. 6 tem a forma apresentada em Eq.
10, com defasagens discretas, assim, o modelo proposto trata-se de um conjunto de
EDDs. Isso faz com que seja necessário o conhecimento não apenas das condições
iniciais, mas também de um histórico temporal das variáveis defasadas. Se
consideradas no domínio da freqüência ao invés do tempo, é como se os sinais
correspondentes às variáveis defasadas tivessem fase diferente do restante.
Sistemas de equações diferenciais com defasagem ocupam atualmente um lugar
de grande importância em diversas áreas, particularmente nas ciências biológicas (por
exemplo, na dinâmica de populações e em epidemologia). O interesse em tais
sistemas surge quando modelos pontuais tradicionais são substituídos por modelos
distribuídos mais realistas, como por exemplo, em um modelo presa-predador onde a
taxa de reprodução não depende apenas da população atual, mas também de um
período anterior, correspondente ao período de gestação.
Neste trabalho, o modelo proposto do sistema elétrico cardíaco, de natureza
distribuída, onde ocorrem acoplamentos entre os elementos do modelo na forma de
realimentação dos estados, levam à ocorrência de defasagens. Isto devido ao tempo
finito de transmissão do sinal elétrico pelos tecidos do coração, levando a defasagens
finitas e, portanto, a equações diferenciais mistas com equações com defasagem,
EDDs.
Como passam a ser envolvidos valores dos estados em instantes passados, é
necessário fornecer uma função que descreva o histórico temporal inicial do sistema, a
fim de definir o valor dos estados antes do instante inicial do intervalo de tempo de
25
interesse. Em muitos modelos, por não se conhecer os valores reais dos estados em
instantes anteriores, a função histórico é definida arbitrariamente como um vetor
constante ou como uma função variável. Por este motivo, por vezes, a EDD e a função
histórico são incompatíveis, isto é, para alguma de suas derivadas, normalmente a
primeira, ocorrem descontinuidades. Quaisquer descontinuidades na função histórico
ou em suas derivadas necessitam de tratamento próprio, pois elas podem se propagar
no tempo.
Visando evitar descontinuidades, a partir do sistema de equações diferenciais com
defasagem, busca-se aqui uma função histórico y0 equivalente à solução do sistema e
que, para instantes de tempo 0 < t < τ, forneça valores correspondentes a instantes de
tempo anteriores ao instante inicial de observação (Boucekkine et al., 1997):
y 0 (t ) = y(t − τ )
(11)
A determinação de y0(t) pode ser feita utilizando-se uma expansão por série de
Taylor, aproximando a solução exata da EDD por uma EDO, da seguinte forma
(Cunningham, 1954):
y (t − τ ) ≅ y (t ) − τ y& (t ) +
τ2
2
&y&(t )
(12)
3.2. MÉTODO DE SOLUÇÃO : MÚLTIPLOS PASSOS
Os métodos de solução para EDDs tendem reduzir as EDDs a EDOS, cuja
solução já é bem explorada. Deste modo, o integrador é baseado no método clássico
de Runge-Kutta de quarta ordem.
No método dos múltiplos passos o intervalo de integração é considerado como
sendo dividido em subintervalos de duração τ. É utilizada a aproximação por EDOs
26
proposta em Eq. 12 apenas para gerar o histórico do primeiro subintervalo, ou seja,
obter valores aproximados da função no intervalo de −τ a 0. A partir daí tais valores
são alimentados nas EDDs, que passam a ser tratadas como EDOs, empregando-se o
método de Runge-Kuta. Quando t torna-se maior que τ, o histórico temporal para o
próximo subintervalo já está disponível, não sendo necessário recorrer-se à
aproximação de Eq. 12 para gerar o hitórico temporal, funcionando de modo
seqüencial, como mostrado na Figura 8.
f ( t, y, y0)
-τ
f ( t, y, yτ)
0
y = yτ
f ( t, y, y2τ)
f ( t, y, y3τ)
t
3τ
τ
2τ
y = y2τ
y = y3τ
Figura 8 – Método dos múltiplos passos.
3.3. SIMULAÇÃO NUMÉRICA
Esta seção trata de simulações numéricas realizadas com o modelo de ritmo
cardíaco proposto. Valores de parâmetros sugeridos por Grudzinski & Zebrowski
(2004) e Santos & Viana (2004) são usados como referência do oscilador do nódulo
SA para o ECG normal. Os outros valores são ajustados a fim de atingir uma
concordância qualitativa com sinais reais de ECG obtidos a partir da segunda
derivação,
os
quais
estão
disponíveis
(http://library.med.utah.edu/kw/ecg/image_index),
no
no
“ECG
banco
learning
aberto
de
center”
dados
“Physionet” (www.physionet.org) e na “ECG library” (http://www.ecglibrary.com). Está
além do escopo desta contribuição uma determinação automática dos parâmetros do
sistema. Aqui, a escolha de parâmetros é feita de um modo artesanal, com o objetivo
de entender os ritmos cardíacos sob um ponto de vista dinâmico e, portanto, o
interesse recai essencialmente sobre a análise qualitativa do sistema.
27
Deste modo, adotam-se os parâmetros da Tabela 2 para representar o
funcionamento normal do coração, bem como algumas patologias cardíacas,
destacando-se o flutter ventricular, a bradicardia sinusal e a fibrilação ventricula.
Em relação a aspectos do acoplamento, considera-se que o coração normal tem
um acoplamento unidirecional do nódulo SA para o nódulo AV e também do nódulo AV
para o HP. Além disso, é necessário estabelecer valores apropriados para as
defasagens nos acoplamentos. Se não há acoplamento, a defasagem não importa já
que ela representa o tempo decorrido quando o sinal viaja de um marca-passo a outro.
Nódulo
SA
Nódulo
Complexo
AV
HP
(a)
(b)
3
mV
2
1
0
-1
4
5
6
7
Tempo (s)
(c)
8
9
10
Figura 9 – ECG normal: (a) modelo conceitual, (b) ECG experimental, (c) ECG
simulado.
28
Nessas condições, efetua-se a simulação do ECG e seus resultados são
apresentados na Figura 9, onde é mostrado o funcionamento esquemático do coração
e o histórico do sinal resultante X = α 0 + α1 x1 + α 3 x3 + α 5 x5 . Além disso, apresenta-se
um ECG experimental medido a partir da segunda derivação, referindo-se a um par de
eletrodos de um sistema de dez eletrodos, com o negativo fixado no membro superior
direito e o eletrodo positivo fixado no membro inferior esquerdo. Esta é uma das
derivações mais características no estabelecimento de diversos diagnósticos, portanto,
bastante representativa do comportamento elétrico cardíaco.
Comparando simulações numéricas com o sinal real é possível observar que eles
apresentam uma boa concordância. A Figura 10 apresenta uma comparação entre o
ECG real e o simulado, mostrando que o modelo captura o comportamento do ECG
real, e que as simulações numéricas apresentam boa concordância com as medidas
reais. Essas conclusões se tornam mais claras observando-se uma ampliação do ECG
normal simulado apresentado na Figura 11, mostrando que ele captura suas
características gerais, apresentando as ondas mais importantes: P, QRS, T.
ECG (mV)
2
1
0
simulado
real
-1
3
3.5
4
4.5
5
5.5
Tempo (s)
6
6.5
7
7.5
Figura 10 – Comparação entre sinais de ECG normais reais e simulados.
29
8
1.5
R
ECG (mV)
1
0.5
P
0
T
Q
S
-0.5
4
4.1
4.2
4.3
4.4
Tempo (s)
4.5
4.6
4.7
Figura 11 – Vista detalhada do ECG normal simulado.
A dinâmica do sistema pode ser melhor compreendida observando o espaço de
fase. Como X& = d ( ECG ) / dt = α1 x2 + α 3 x4 + α 5 x6 é possível observar na Figura 12
que o espaço de fase do ECG normal é essencialmente periódico. Na mesma figura
também são mostradas projeções bi-dimensionais do espaço de fase hexadimensional, representando a dinâmica de cada oscilador (SA, AV e HP). A
característica periódica do comportamento do sistema é novamente perceptível.
30
6
4
d(ECG)/dt
2
0
-2
-4
-6
-0.5
0
10
0.5
ECG
1
1.5
(a)
3
10
2
5
1
0
0
X6
0
X4
X2
5
-1
-10
-2
-15
-5
-10
-3
-2
-1
X1
0
1
2
-5
-3
-2.5
-2
-1.5
(b)
X3
-1
-0.5
0
-20
-3
(c)
-2
-1
X5
0
1
2
(d)
Figura 12 – Espaços de fase relacionados ao ECG normal: (a) ECG; (b)
1ooscilador; (c) 2º oscilador; (d) 3º oscilador.
A análise a seguir trata de diferentes termos de acoplamento a fim de simular
algumas patologias cardíacas identificadas a partir do ECG. Primeiramente considerase o flutter ventricular, caracterizado principalmente pela falta de funcionamento do
nódulo sinoatrial. Neste caso, não há comunicação entre os nódulos SA e AV, o que é
simulado eliminando-se o acoplamento entre o primeiro e o segundo osciladores
(kAV-SA = 0), enquanto todos os outros parâmetros permanecem os mesmos do ECG
normal. Desta forma, o único termo de acoplamento não nulo é kHP-AV = 20. A Figura
31
13a apresenta o modelo esquemático dos acoplamentos existentes no caso desta
patologia. Sob tais condições, o sistema é conduzido pelo nódulo AV, apresentando
uma freqüência mais alta em relação ao ECG normal. Essa patologia é chamada de
flutter ventricular, sendo caracterizada pelo fato da resposta dos ventrículos passar a
ser diferente, não apenas em freqüência – devido ao ritmo mais rápido do nódulo AV –
mas também na forma do sinal, cujo padrão característico é representado pelo ECG
mostrado na Figura 13b, que apresenta grande similaridade com a simulação
apresentada na Figura 13c. Observe que o ECG tem um comportamento regular, o
que é confirmado nos espaços de fase e em algumas de suas projeções bidimensionais apresentadas na Figura 14.
Nódulo
SA
Nódulo
Complexo
AV
HP
(a)
(b)
mV
2
1
0
-1
4
5
6
7
Tempo (s)
8
9
10
(c)
Figura 13 – ECG sem o acoplamento SA-AV (flutter ventricular): (a) Modelo
conceitual; (b) ECG experimental e (c) ECG simulado.
32
8
6
d(ECG)/dt
4
2
0
-2
-4
-6
-1
0
1
ECG
2
(a)
1
15
8
10
6
0.5
5
0
0
X6
2
X4
X2
4
-5
0
-2
-10
-0.5
-4
-15
-6
-3
-2
-1
X1
(b)
0
1
-1
-0.3
-0.2
-0.1
X3
0
0.1
0.2
-20
-3
(c)
-2
-1
X5
0
1
2
(d)
Figura 14 - Espaços de fase relacionadas ao flutter ventricular: (a) ECG; (b)
1ooscilador; (c) 2º oscilador; (d) 3º oscilador.
O acoplamento AV-SA é investigado a seguir, considerando que o único termo de
acoplamento não nulo é kAV-SA=5, o que significa que kHP-AV se anula. Sob tais
condições, a influência dos ventrículos é diferente e o ritmo cardíaco é representado
pelo ECG mostrado na Figura 15. Este tipo de patologia é chamado de bradicardia
sinusal, ocorrendo quando é perdida a comunicação do nódulo átrio-ventricular para o
33
feixe de His. Ela apresenta um comportamento regular, como se pode observar nos
espaços de fase e em algumas de suas projeções bi-dimensionais apresentadas na
Figura 16. É verificada uma boa concordância qualitativa entre o ECG experimental e
o simulado, visto que todas as características topológicas do primeiro são
reproduzidas no segundo.
Nódulo
SA
Nódulo
Complexo
AV
HP
(a)
(b)
mV
2
1
0
-1
4
5
6
7
Tempo (s)
8
9
10
(c)
Figura 15 – ECG sem o acoplamento HP-AV (bradicardia sinusal). (a) Modelo
conceitual; (b) ECG experimental e (c) ECG simulado.
34
2
d(ECG)/dt
1.5
1
0.5
-2.2 -2.1
-2 -1.9
ECG
-1.8
(a)
0
8
2
1
2
0
0
-0.01
X6
4
X4
X2
6
-0.02
-1
-2
-0.03
-2
-4
-6
-3
-2
-1
X1
(b)
0
1
-3
-2.5
-2
-1.5
X3
-1
-0.5
-0.04
-7
(c)
-6.9
-6.8
X5
-6.7
(d)
Figura 16 – Espaços de fase relacionadas a bradicardia sinusal: (a) ECG; (b)
1ooscilador; (c) 2º oscilador; (d) 3º oscilador.
Neste ponto, considera-se mais uma vez o ECG normal, no entanto, marca-passos
externos estão excitando o sistema. Portanto, são considerados os mesmos
parâmetros relacionados ao ECG normal (Figura 9-11) e diferentes parâmetros de
forçamento: ρSA = 1, ρAV = 1, ρHP = 20, ωSA = ωAV = ωHP = 2π/(60/70) (ver Tabela 2).
Nessas condições, induz-se uma fibrilação ventricular. Esta resposta patológica é
causada por diferentes estimulações ao ventrículo, sendo caracterizada por um ECG
irregular com uma resposta QRS rápida. A Figura 17 mostra a simulação numérica
relacionada a esta condição juntamente com o sinal de ECG correspondente.
35
Γ1
Nódulo
SA
Γ2
Γ3
Complexo
Nódulo
HP
AV
(a)
(b)
mV
2
1
0
-1
4
5
6
7
Tempo (s)
8
9
10
(c)
Figure 17 – ECG com fibrilação ventricular. (a) Modelo conceitual; (b) ECG
experimental e (c) ECG simulado.
Na Figura 18 são apresentados espaços de fase e algumas de suas projeções bidimensionais. A característica irregular deste tipo de comportamento é perceptível já
que a resposta não está relacionada a uma curva fechada. Analisando as seções de
Poincaré, traçadas segundo a freqüência de amostagem de 2π/(60/70) Hz,
correspondente
à
freqüência
cardíaca
normal,
observa-se
uma
estrutura
aparentemente caótica (Figura 19). Já que a resposta parece ser caótica, esta
conclusão deve ser confirmada pela análise de algum invariante do sistema como os
expoentes de Lyapunov.
36
6
d(ECG)/dt
4
2
0
-2
-4
-6
-8
-1
0
1
ECG
2
(a)
4
20
10
2
10
5
0
0
X2
X4
X6
15
0
-2
-10
-5
-4
-20
-10
-3
-2
-1
X1
(b)
0
1
2
-6
-3
-2
-1
X3
0
1
-30
-3
-2
(c)
-1
X5
0
1
(d)
Figura 18 – Espaços de fase relacionados à fibrilação ventricular: (a) ECG; (b)
1ooscilador; (c) 2º oscilador; (d) 3º oscilador.
37
2
8
6
d(ECG)/dt
4
2
0
-2
-4
-6
-8
-0.5
0
0.5
ECG
1
1.5
(a)
4
20
2
10
0
0
10
0
-5
-3
-2
-1
X1
(b)
0
1
2
X6
X2
X4
5
-2
-10
-4
-20
-6
-3
-2
-1
X3
0
1
-30
-3
(c)
-2
-1
X5
0
1
2
(d)
Figura 19 – Seções de Poincaré relacionadas à fibrilação ventricular: (a) ECG; (b) 1º
oscilador; (c) 2º oscilador; (d) 3º oscilador
38
Tabela 2 – Parâmetros de Simulação
Variáveis
ECG Normal
Fibrilado
Bradicardia
Flutter
Primeiro oscilador
aSA
1/(mV.s)
3
3
3
3
wSA1
(mV)
0.2
0.2
0.2
0.2
wSA2
(mV)
-1.9
-1.9
-1.9
-1.9
dSA
(mV)
3
3
3
3
eSA
(mV)
4.9
4.9
6
4.9
Segundo oscilador
aAV 1/(mV.s)
3
3
3
3
wAV1
(mV)
0.1
0.1
0.1
0.1
wAV2
(mV)
-0.1
-0.1
-0.1
-0.1
dAV
(mV)
3
3
3
3
e AV
(mV)
3
3
3
3
Terceiro oscilador
aHP 1/(mV.s)
5
5
5
5
wHP1
(mV)
1
1
1
1
wHP2
(mV)
-1
-1
-1
-1
dHP
(mV)
3
3
3
3
eHP
(mV)
7
7
7
7
Acoplamentos ( 1/s2 )
KSA-AV
0
0
0
0
KAV-HP
0
0
0
0
KHP-SA
0
0
0
0
KAV-SA
-5
-5
-5
0
KHP-AV
-20
-20
0
-20
KSA-HP
0
0
0
0
Defasagem (s)
0
0
0
0
τSA-AV
0
0
0
0
τSA-HP
0.8
0.8
0.8
0.8
τAV-AS
0
0
0
0
τAV-HP
0
0
0
0
τHP-AS
0.1
0.1
0.1
0.1
τHP-AV
2
Amplitudes (ρ – mV/s ) e frequencias (ω – Hz ) das entradas
0
1
0
0
ρAS
0
1
0
0
ρAV
0
20
0
0
ρHP
ωSA
0
0
0
2π/(60/70)
ωAV
0
0
0
2π/(60/70)
ωHP
0
0
0
2π/(60/70)
Ainda: α0 = 1 mV e α1 = 0,1; α3 = 0,05; α5 = 0,4 adimensionais
A 2ª parcela de cada oscilador é multiplicada por um termo unitário de unidades 1/(mV2.s2)
apenas para compatibilizar unidades.
39
4. ANÁLISE DE SÉRIES TEMPORAIS
Sistemas dinâmicos podem ser classificados em contínuos ou discretos. No
primeiro caso a sua dinâmica obedece à Eq. 13, onde a observação é feita de modo
contínuo, aplicando-se a muitos sistemas presentes na natureza.
(13)
Já no segundo caso, a dinâmica se desenvolve de forma discreta, o sinal de
interesse é observado por meio de amostras. Um sistema contínuo pode ser
observado por meio de amostras de suas variáveis o que o torna um sistema discreto,
conforme a seguir:
∆ ∆ 1 (14)
Na Eq. 14 as amostras são tomadas em intervalos de tempo iguais, isto é,
segundo um período de amostragem ∆t. Assim, esta dinâmica pode ser referenciada
por um índice
n,
relativo à ordem da amostra dentro do conjunto tomado. A um
conjunto de dados tomados ao longo do tempo, a partir da observação de algum
sistema, dá-se o nome de série temporal.
Deve-se observar, no entanto, que ao serem efetuadas medições sobre um
sistema, os valores medidos podem não corresponder diretamente aos estados do
sistema, mas sim a um valor dado por uma função s(.), que mapeia os pontos do
espaço de estados do sistema em números reais, os quais correspondem às variáveis
observadas pelo aparato de medição.
40
Assim, ao se efetuar um conjunto medidas a partir de um experimento, obtém-se
uma série temporal Sn , composta por elementos s(n), dados por Eq. 15, onde ξ(n)
representa uma perturbação na medição, ou seja, um ruído.
s (15)
No caso do eletrocardiograma, são medidos sinais de potencial elétrico a partir da
superfície do corpo humano, resultantes do dipolo elétrico formado por seus nódulos
de células marca-passo. Por este motivo, o valor medido corresponde ao valor do
campo resultante em pontos específicos do corpo humano e não diretamente aos
estados do sistema.
Devido ao uso de equipamentos digitais para o armazenamento dos dados, o ECG
é comumente tratado como um sinal discreto, formando uma série temporal. Por tratarse de um sistema de medição de sinais elétricos, sinais externos podem facilmente
influenciar a medição, inclusive sinais na forma de ondas eletromagnéticas, visto que o
corpo humano capta várias dessas ondas, como uma antena.
No contexto de análise dos ritmos cardíacos por meio do ECG, não
necessariamente tem-se que conhecer detalhes sobre a dinâmica do coração para
avaliar os seus ritmos. A análise de séries temporais é uma abordagem que pode ser
bastante útil neste contexto.
A idéia básica da reconstrução do espaço de estados é a de que, um sinal de
saída do sistema, contém informação sobre variáveis não-observáveis, que pode ser
usada para estimar o estado atual. Portanto, uma série temporal escalar Sn pode ser
utilizada para construir uma série temporal vetorial cuja topologia seja equivalente à da
dinâmica do sistema em análise.
Para se reconstruir o espaço de estados do sistema é necessário formar um
conjunto de pontos ordenado que capture a estrutura das órbitas deste espaço. Isto
41
pode ser feito por meio da defasagem da série temporal, formando séries defasadas
S n +δ i , onde δi é a defasagem da série que formará a componente i do espaço
reconstruído (Takens, 1981). Tomando-se uma defasagem adequada, é possível gerar
um conjunto de defasagens para criar um espaço de dimensão De como mostrado
abaixo:
u = {S n −( De−1)δ , S n −( De−2 )δ , . . . , S n−δ , S n }
(16)
A literatura relata muitos métodos empregados para determinar o tempo de
defasagem. O Método da Informação Mútua (Fraser & Swinney, 1986) apresenta bons
resultados, o que disseminou seu uso. Para a determinação da dimensão de imersão
De também há diversos métodos, tendo o Método dos Falsos Vizinhos (False Nearest
Neighbors - FNN) (Kennel et al., 1992) se apresentado como uma boa opção.
4.1. MÉTODO DA INFORMAÇÃO MÚTUA MÉDIA
Fraser & Swinney (1986) estabelecem que a defasagem δ está associada ao
primeiro mínimo local da função de informação mútua média I(δ), a qual é definida
como:
 P(s(n), s(n + δ )) 
N −δ
I (δ ) =
∑ P(s(n), s(n + δ ))log  P(s(n) ) P(s(n + δ ))
2
(17)
n =1
Para um melhor entendimento do significado da informação mútua, cabe aqui uma
explicação mais detalhada de cada um dos termos da função que a define: P(s(n)) é a
probabilidade de que a medida s(n) apresente determinado valor. Considerando-se
um domínio, possivelmente o intervalo sobre o qual se distribuem os valores de Sn ,
42
esta probabilidade pode ser calculada pela razão entre o número de ocorrências do
valor de s(n) na série e N, o número total de elementos da série. Do mesmo modo,
P(s(n+δ)) é a probabilidade da ocorrência de determindao valor na medição s(n+δ),
onde este é um termo qualquer da série, localizado δ amostras depois de s(n) dentro
da ordenação de Sn . O termo P(s(n),s(n+δ)) refere-se à probabilidade conjunta das
medidas s(n) e s(n+δ), ou seja, a probabilidade de um dado valor observado no
instante n seja observado novamente em um instante n+δ.
Desta forma, I(δ) expressa a probabilidade de ocorrência simultânea de um valor
em uma posição n e na correspondente posição n+δ, independente de qual seja o
valor, pois foi adimensionalizado pelas respectivas probabilidades absolutas.
Assim, a expressão da informação mútua fornece, para um dado valor de
defasagem, o somatório das probabilidades da repetição de valores após um intervalo
δ, assemelhando-se a uma forma de autocorrelação do sinal. Verifica-se que quando
as medidas s(n) e s(n+δ) são completamente independentes tem-se I(δ) = 0. Por outro
lado, quando s(n) e s(n+δ) são iguais, I(δ) é máximo.
A Figura 22 apresenta uma representação esquemática de modo a representar o
significado físico da informação mútua. Existem duas curvas para auxiliar na análise,
com defasagens δ1 e δ2, cada uma. No esboço toma-se um ponto s(j) da curva,
verificando-se que para uma defasagem δ1 o correspondente s(j+δ1) não apresenta
valor muito próximo ao seu. Quando é analisado o ponto s(j+δ2), verifica-se que seu
valor e o do ponto s(j) são coincidentes, indicando um máximo de similaridade,
indicando haver mais informação em comum entre as duas curvas, ou seja, uma maior
informação mútua entre Sn e Sn +δ 2 . O que se deseja é recompor o espaço de estados
do sistema, o que é feito pala geração de suas componentes a partir da aplicação de
diferentes defasagens à série temporal produzida pelo sistema em análise. O espaço
de estado reconstruído guarda certa quantidade de informação, sendo tanto mais
43
semelhante ao espaço de estado real do sistema, quanto maior for a quantidade de
informação que se consegue extrair da série temporal.
δ2
Sn
δ1
j
j+δ1
j+δ2
n
Figura 20 – Representação esquemática da informação mútua.
Vê-se então que a minimização da quantidade de informação mútua (repetida)
entre os sinais defasados leva a uma maior quantidade de informação nova no espaço
reconstruído, causando uma maior similaridade sua com espaço real. Por isso,
escolhe-se uma defasagem que leve a um mínimo de informação mútua entre a série
e sua cópia defasada. Toma-se o menor valor de defasagem que leve a um mínimo
local da informação mútua, pois o aumento da defasagem tende a diminuir a
quantidade de informação contida no espaço reconstruído, repetida ou não. Isto ocorre
porque, uma maior defasagem faz com que trechos menores da série sejam tomados
para formar as componentes do espaço reconstruído, carregando consigo menos
informação.
4.2. MÉTODO DOS FALSOS VIZINHOS
O algoritmo dos Falsos Vizinhos (False Nearest Neighbors – FNN) foi
originalmente desenvolvido para determinar o número de coordenadas defasadas
necessárias para se recriar uma dinâmica autônoma, mas ele é extendido para
examinar o problema da determinação da dimensão de imersão apropriada.
44
Em geral não se consegue medir diretamente os valores dos estados do
sistema, mas apenas uma série temporal resultante de um mapeamento deles em um
espaço de menor dimensão. Em uma dimensão de que seja muito baixa para
representar o atrator formado pela dinâmica real do sistema, nem todos os pontos que
se situam próximos uns aos outros são realmente vizinhos, por causa dos processos
dinâmicos do sistema, que não são ali corretamente retratados. Alguns pontos se
localizam na verdade longe uns dos outros e simplesmente aparecem como vizinhos
porque a estrutura geométrica do atrator foi projetada em um espaço de menor
dimensão (Kennel et al., 1992).
De = 1
De = 3
De = 2
Figure 21 – Falsos vizinhos.
Basicamente, o método dos falsos vizinhos consiste na análise do percentual de
falsos vizinhos que surgem ao se passar de uma dimensão para a imediatamente
superior. Após o cálculo da defasagem adequada para a reconstrução, são testadas
diferentes dimensões, buscando-se a mais baixa delas que forneça um pequeno
percentual de falsos vizinhos.
Para implementar o método, considera-se um espaço reconstruído com uma
dimensão D, adotando-se índices sobrescritos para indicar a dimensão utilizada.
Assim, o ponto uD(n) tem uD(r) como seu vizinho mais próximo na dimensão D, de
modo que a fração de falsos vizinhos passa a ser dada por:
45
 u D +1 (n ) − u D +1 (r )
 σ
 U  − u D (n ) − u D (r ) 
U
−
R
D
D
 u (n ) − u (r )
 R



N − D −1
σ

U  − u D (n ) − u D (r ) 
R

n =1
N − D −1
∑
F ( D, R ) =
n =1
∑
(18)
onde U(.) é a função degrau unitário e σ é o desvio padrão dos valores da série, aqui
utilizado para previamente eliminar pontos muito afastados, que não poderiam ser
vizinhos. A função degrau é utilizada para que os somatórios da equação acima se
tornem contadores.
Calcula-se a distância de cada ponto uD(n) ao seu vizinho mais próximo uD(r).
Passa-se então da dimensão D para a D+1, isto é, reconstrói-se o espaço por meio de
coordenadas defasadas para esta nova dimensão. São tomadas as projeções destes
mesmos pontos no novo espaço, respectivamente uD+1(n) e uD+1(r), calculando-se a
distância entre elas. Se a razão entre as distâncias nas dimensões D e D+1
ultrapassar um valor limite R, uD(r) é considerado um falso vizinho.
Apenas quando esta distância mudar radicalmente de uma dimensão para a outra
é que um ponto será considerado um falso vizinho, sendo o limite desta mudança
dado pelo fator R. No pacote Tisean é sugerido um valor para R igual a 2, ou seja,
para que um ponto seja considerado um falso vizinho, a distância dele ao ponto de
referência teria que dobrar após o aumento de dimensão.
4.3. EXPOENTES DE L YAPUNOV
Os expoentes de Lyapunov são a taxa exponencial média de divergência ou
convergência de órbitas próximas no espaço de estados de um sistema, de modo que
terminam por avaliar a sensibilidade do sistema a variações de condições iniciais. Os
sinais dos expoentes de Lyapunov fornecem uma visão qualitativa da dinâmica do
46
sistema, estando valores positivos associados à divergência das órbitas e os negativos
à sua convergência.
A existência de expoentes de Lyapunov positivos indica direções de instabilidade
local na dinâmica do sistema (Franca & Savi, 2003), levando a órbitas divergentes.
Como órbitas próximas representam estados similares, a divergência das órbitas
significa que, após certo tempo, comportamentos bastante diferentes surgirão a partir
de cada uma. Assim, perde-se a capacidade de predizer o comportamento do sistema
quando houver pelos menos um expoente positivo, sendo o sistema classificado como
caótico.
Para o cálculo do expoente de Lyapunov, Wolf et al. (1985) propõem a seguinte
expressão:
 ui (t ) − uiε (t )
1
λi = lim log 2 
t →∞ t
 u ( 0) − u ε ( 0)
i
 i




(19)
onde o índice i refere-se à i-ésima componente de um ponto do espaço de estados,
situando-se o estado u(t) sobre uma trajetória de referência e uε(t) pode ser
considerado este estado após uma pequena perturbação. Da equação anterior vê-se
que o expoente de Lyapunov quantifica a velocidade com que a informação é criada
ou destruída pela dinâmica do sistema, sendo expresso em informação por unidade de
tempo.
De modo a esclarecer este conceito, pode-se considerar um exemplo simples:
estabelecida uma condição inicial com uma incerteza da ordem de 10-6 vezes o seu
valor, em quanto tempo esta informação estaria completamente degradada pela
dinâmica do sistema? A incerteza inicial pode ser vista como uma perturbação
(denominador igual a 10-6 em Eq. 19) e para que se degrade completamente a
informação, a incerteza deve atingir a ordem de grandeza da mesma (numerador igual
47
a 1 na mesma equação). Assim, o intervalo de tempo necessário para ocorrer a
degradação seria igual a 19,93 / λi.
Deve-se notar que o valor do expoente depende da base logarítmica utilizada, o
que pode variar entre autores, bastando utilizar sempre a mesma base para manter
uma coerência de escalas. No entanto, os sinais dos expoentes são independentes da
base, pois dependem apenas da divergência ou convergência das órbitas, sendo
ditados pela topologia do espaço de estados.
Por outro lado, o algoritmo de Wolf pode apresentar resultados errôneos,
especialmente por considerar que há de fato divergência exponencial. O algoritmo
proposto por Kantz (1994) utiliza as mesmas idéias apresentadas, porém, toma o
atrator reconstruído e examina a divergência orbital em escalas de comprimento. O
método monitora a evolução em longo prazo de um único par de órbitas próximas e
está apto a estimar expoentes de Lyapunov não-negativos. Kantz (1994) considera
que a taxa de divergência das trajetórias flutua ao longo das mesmas, sendo esta
flutuação dada pelo espectro de expoentes de Lyapunov efetivos. A média dos
expoentes de Lyapunov efetivos ao longo da trajetória fornece o verdadeiro valor do
máximo expoente de Lyapunov. Para estimá-lo, toma-se como uma boa aproximação
o logarítimo L da média dos afastamentos, dado por:
∆ ∑!
# $%&'
! ( -
( %
∑'+&'( |* ∆ ∆|,
(20)
Nesta implementação toma-se como referência um ponto u(n0) do espaço de
imersão e em torno dele é estabelecida uma vizinhança U, de raio ε, equivalente a
uma perurbação na referência e tão pequeno quanto possível. Tendo a última
componente de cada ponto u(n) em U como base, toma-se um ponto ∆n amostras à
frente na série, calculando-se então a média das distâncias em relação à referência
u(n0), a partir da última componente.
48
Daí troca-se a referência u(n0) e o processo é repetido para N elementos da série,
possivelmente todos, e depois é calculanda a média dos logarítmos. Partindo de Eq.
19 verifica-se que, para realmente calcular os expoentes, é necessário dividir uma
razão de deformação pelo intervalo de tempo (∆n), no qual se permitiu que o sistema
evoluísse.
Paralelamente, para que a equação anterior exprima uma razão de deformação, é
necessário que o afastamento médio seja divido pelo afastamento inicial. No entanto,
se o raio da vizinhança for suficientemente pequeno, pode-se considerar que a
distância inicial para todos os pontos é igual a ε. Isto acarreta apenas uma translação
da curva, de valor –ln(ε), ao longo do eixo vertical.
Vê-se assim que, diferentemente do que ocorre no algoritmo de Wolf, aqui são
obtidos logarítmos de afastamentos médios. Se o crescimento de tais afastamentos
obedecer a uma lei exponecial, então a curva de seus logaritmos pelo intervalo ∆n se
aproxima de uma reta, onde o coeficiente angular é o máximo expoente de Lyapunov.
É importante notar que, desta forma, a translação não altera o valor dos expoentes.
Figura 22 – Expoentes de Lyapunov segundo o algoritmo de Kantz.
49
Também deve ser destacado que é neste ponto que os algoritmos de Kantz e Wolf
diferem, pois neste último, parte-se da premissa de que existe uma lei exponencial de
crescimento, sendo que ele já calcula diretamente o suposto expoente. Já o algoritmo
de Kantz utiliza o logaritmo da média dos afastamentos, tendo o usuário que analisar
as curvas no gráfico L x ∆n e verificar a presença de uma lei exponencial de
crescimento do afastamento.
Fisicamente é como se, a partir de um espaço de imersão, fosse tomado um ponto
de referência u(n0), onde um conjunto de pontos u(n) de órbitas próximas seriam seu
estado perturbado. A partir daí as órbitas evoluem por um intervalo de tempo
correspondente a ∆n amostras. O logaritmo da média das distâncias entre os novos
pontos e a nova referência aproxima um fator de deformação causada pelo sistema,
sendo feita uma média deste valor com diferentes pontos de referência para evitar a
influência de fatores locais.
4.4. ANÁLISE DE ECGS GERADOS A PARTIR DO MODELO MATEMÁTICO
Esta seção trata os sinais de ECG como séries temporais geradas através do
modelo matemático proposto. Basicamente, consideram-se dois sinais diferentes:
ECG normal e ECG associado à fibrilação ventricular. A ideia é estabelecer uma
comparação entre os dois comportamentos. Os expoentes de Lyapunov são usados
para caracterizar a dinâmica do sistema sendo úteis para identificar os tipos de
patologias dos ritmos cardíacos.
Inicialmente considera-se um sinal associado a um ECG normal com 20.000
amostras. Os parâmetros de reconstrução são analisados através dos métodos da
informação mútua e dos falsos vizinhos. A Figura 23 apresenta os resultados desta
análise, mostrando que a defasagem deve ser δ = 0.4s (400 amostras) e a dimensão
de imersão deve ser De = 3. Apesar do primeiro mínimo não estar claramente definido,
50
é utilizado um mínimo global da primeira região, definida pelo primeiro máximo local.
Assim, o método das coordenadas defasadas pode ser empregado para a
reconstrução do espaço de estados. A Figura 24 apresenta o espaço reconstruído
junto ao espaço obtido diretamente a partir de simulação numérica, mostrando uma
boa concordância, visto que ambos apresentam as mesmas características
topológicas. A Figura 25 apresenta o maior expoente de Lyapunov avaliado pelo
algoritmo de Kantz, indicando um valor nulo, o que é esperado para respostas
periódicas.
5
4.5
Informaçao Mútua
4
3.5
3
2.5
2
1.5
1
0.5
0
0
500
1000
1500
2000
2500
3000
Defasagem (amostras)
(a)
0.7
Fração de Falsos Vizinhos
0.6
0.5
0.4
0.3
0.2
0.1
0
1
2
3
4
5
6
Dimensão de Imersão (De)
7
8
9
(b)
Figura 23 – Análise da defasagem (a) e dimensão de imersão (b) para o ECG
normal simulado.
51
1.5
E C G (t-0 .0 8 )
1
0.5
0
-0.5
-0.5
0
0.5
ECG(t)
1
0.5
ECG
1
1.5
(a)
6
4
d(ECG)/dt
2
0
-2
-4
-6
-0.5
0
1.5
(b)
Figure 24 – Espaço de fase do ECG normal utilizando o modelo do coração:
reconstruído a partir da série gerada pelo modelo (a) e traçado a partir de seus
estados (b).
52
Logarítmo do Fator de Deformação
-6
-6.5
-7
De=3
De=4
De=5
De=6
-7.5
-8
0
10
20
30
40
50
Número de Iterações
Figura 25 – Expoentes de Lyapunov do ECG normal simulado.
A análise a seguir é dedicada à análise de um ECG associado a uma fibrilação
ventricular. Os parâmetros utilizados foram os mesmos da simulação representada na
Figura 17. Mais uma vez foram geradas por simulação com 20.000 amostras.
O sinal de ECG é utilizado para executar a reconstrução do espaço de estado. A
análise dos parâmetros de defasagem é apresentada na Figura 26, mostrando que a
defasagem de reconstrução deve ser δ = 0,100s (42 amostras) e que a dimensão de
imersão é De = 10. Embora o primeiro mínimo não esteja claramente definido, é
utilizado o mínimo global da primeira região, definido pelo primeiro máximo local. Esse
resultado indica o aumento da complexidade do sinal. Com o emprego desses
parâmetros, o método das coordenadas defasadas pode ser empregado para
reconstruir o espaço de estado. A Figura 27 apresenta o espaço reconstruído,
indicando uma boa concordância com aquele obtido a partir de simulação numérica. A
Figura 28 apresenta o maior expoente de Lyapunov calculado a partir do algoritmo de
53
Kantz, indicando um valor λ=0.2286, confirmando que a fibrilação ventricular está
relacionada ao caos.
2.5
Informação Mútua
2
1.5
1
0.5
0
0
50
100
150
200
250
300
350
400
450
500
Defasagem (amostras)
(a)
1
0.9
Fração de Falsos Vizinhos
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
1
2
3
4
5
6
7
8
Dimensão de Imersão (De)
9
10
11
12
(b)
Figura 26 – Análise dos parâmetros de defasagem para a fibrilação ventricular do
ECG (a) Informação Mútua e (b) Número de Falsos Vizinhos.
54
2
1.5
2
1.5
1
ECG(t-tao)
ECG(t-2tao)
1
0.5
0
0.5
0
-0.5
-0.5
-0.5
-2
ECG(t-tao)
0
2
2
0.5
1
ECG(t)
1.5
0
1
1.5
ECG(t)
-0.5
0
0.5
2
1.5
2
1
0
0.5
-0.5
ECG(t-9tao)
1.5
2
1
1.5
ECG(t-6tao)
ECG(t-4tao)
2
0.5
0
1
0.5
0
-0.5
-0.5
-0.5
-0.5
0
0.5
ECG(t-2tao)
1
1.5
2
1.5
2
0
0.5
1
-0.5
-1
0
0.5
0
1
ECG(t)
1
1.5
2
2
ECG(t-2tao)
ECG(t-4tao)
Figura 27 – Espaço de estado reconstruído a partir de um ECG simulado com
fibrilação ventricular.
Logarítmo do Fator de Deformação
-4
-4.4
-4.8
De=3
De=4
De=5
De=6
-5.2
-5.6
5
10
15
20
25
30
35
40
45
50
Número de Iterações
Figura 28 – Expoentes de Lyapunov p/ ECG simulado com fibrilação ventricular.
55
4.5. ANÁLISE DE ECGS EXPERIMENTAIS
Com base nos resultados da análise dos ECGs gerados a partir do modelo
matemático, pretende-se conduzir a análise de ECGs experimentais reais. A princípio
seria de se esperar que os resultados obtidos fossem similares, porém, a presença de
ruído nos sinais experimentais causa, em geral, alteração nas características dos
sinais.
As análises consideram séries temporais de eletrocardiogramas retiradas do banco
de dados aberto Physionet (endereço http://www.physionet.org/physiobank/database/).
Para o ECG normal utiliza-se uma série com um total de 7680 amostras. A Figura 29
ECG (mV)
apresenta o ECG normal utilizado na análise.
2
0
-2
0
1
2
3
4
5
6
7
8
9
10
Tempo (s)
Figura 29 – ECG normal experimental.
A Figura 30 apresenta os resultados da análise dos parâmetros de imersão:
informação mútua média (esquerda) e falsos vizinhos (direita), mostrando indicando os
seguintes resultados: δ = 0.18s (45 amostras) e De= 20. A partir daí, passa-se a avaliar
os expoentes de Lyapunov. O cálculo do máximo expoente de Lyapunov através do
algoritmo de Kantz está apresentado na Figura 31 para diferentes valores da
dimensão de imersão. Note que existe uma tendência para um expoente nulo, no
entanto, seria difícil afirmar isso sem nenhuma informação a respeito da série.
56
1.4
1.2
Informação Mútua
1
0.8
0.6
0.4
0.2
0
0
100
200
300
400
500
600
700
800
900
1000
Defasagem (amostras)
(a)
0.8
Fração de Falsos Vizinhos
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
5
10
15
20
25
Dimensão de Imersão (De)
(b)
Figura 30 – ECG normal experimental: informação mútua (a) e falsos vizinhos (b).
57
-1
Logaritmo do fator de Deformação
-1.5
-2
De=1
De=2
De=3
De=4
De=5
De=7
De=8
De=9
De=10
De=11
De=12
De=13
-2.5
-3
-3.5
-4
-4.5
-5
5
10
15
20
25
30
Número de Iterações
35
40
45
50
Figura 31 – Máximo expoente de Lyapunov para o ECG normal experimental.
Em seguida é analisado um sinal de ECG diagnosticado como uma fibrilação
ECG (mV)
ventricular, notando-se sua semelhança com a curva mostrada na Figura 17.
2
0
-2
50
52
54
56
58
60
62
64
66
68
70
Tempo (s)
Figura 32 – ECG fibrilado experimental.
Na Figura 35 apresenta-se a análise dos parâmetros de imersão, tempo de
defasagem e dimensão de imersão, respectivamente através da informação mútua
média e os falsos vizinhos, obtendo-se valores de δ = 0,84s (210 amostras) e De = 20
para este eletrocardiograma. O expoente de Lyapunov máximo está apresentado na
Figura 34. Mais uma vez é difícil identificar o valor do expoente, mas existe uma
58
tendência que sugere um expoente positivo, conforme pode ser observado na linha
inclinada.
2.5
Informação Mútua
2
1.5
1
0.5
0
0
500
1000
1500
2000
2500
3000
Defasagem (amostras)
(a)
1
0.9
Fração de Falsos Vizinhos
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
2
4
6
8
10
12
14
Dimensão de Imersão(De)
16
18
20
(b)
Figura 33 – ECG experimental com fibrilação ventricular: informação mútua (a) e
falsos vizinhos (b).
59
-1
Logaritmo do Fator de Deformação
-2
-3
De=1
De=2
De=6
De=7
De=8
De=9
De=10
De=11
De=12
De=13
De=14
De=15
De=16
-4
-5
-6
-7
0
5
10
15
20
25
30
Número de Iterações
35
40
45
50
Figura 34 – Máximo expoente de Lyapunov para o ECG experimental com
fibrilação ventricular.
Tendo em vista a dificuldade em se avaliar os expoentes de Lyapunov a partir
de dados experimentais, ferramentas auxiliares devem ser empregadas. Parte da
divergência entre o ECG real e o simulado deve-se à presença de perturbações
durante a medição representada pelo ruído.
Desta forma, torna-se necessário uma melhor caracterização do sitema a partir
dos valores medidos no ECG. Para tanto o próximo capítulo apresenta uma discussão
sobre processos de filtragem aos dados experimentais, o que permite diminuir a
influência do ruído na análise dinâmica.
60
5.
ANÁLISE
DE
ECGS
EXPERIMENTAIS
CONSIDERANDO
PROCEDIMENTOS DE FILTRAGEM
Todo dado experimental é contaminado por ruído. O nível de ruído presente em
um sinal pode prejudicar a sua análise, impossibilitando conclusões a respeito de sua
dinâmica. Neste contexto, torna-se necessário a separação do sinal e de seu ruído o
que é determinado a partir das naturezas do sistema observado e do ruído.
O ruído pode apresentar uma correlação com a dinâmica do sistema observado ou
ser independente. Enquadra-se neste último caso o chamado ruído de medição que se
refere à corrupção das medidas durante o processo de observação por fatores
independentes da dinâmica do sistema. Em geral, o ruído de medição ocorre pela
presença de fatores externos ao sistema que afetam o aparato de aquisição de dados,
estando fora do controle ou da percepção do usuário.
A filtragem é o processo de separação entre grupos de dados. Várias abordagens
podem ser utilizadas com este objetivo e um modo de tentar separar o ruído e o sinal
de interesse é estimar o valor real deste último, estabelecendo uma distinção do sinal
de ruido.
Este capítulo apresenta técnicas de filtragem, selecionando a mais adequada ao
sinal experimental do ECG. Em seguida aplica-se o filtro ao sinal experiemental,
procedendo-se a seguir à análise do sistema.
5.1. REDUÇÃO LINEAR DE RUÍDO
Em geral, não é possível aplicar um processo linear de filtragem diretamente a
uma série temporal no domínio do tempo com o intuito de obter uma estimativa
. do estado discretizado X(n). Isto pode ser visto ao tomar-se a Eq. 17 e multiplicála por um fator θa, o que é equivalente a aplicar um filtro linear de ganho θa sobre o
61
sinal observado. Caso o ruído apresente uma magnitude muito menor que a
observação, ela pode ser desprezada, como mostrado em Eq. 21.
/0 /0 1 /0 2 /0 1
(21)
Em uma implementação real isso ocorre porque o aparato de medição possui uma
dada resolução, que sendo maior que a amplitude do ruído após a filtragem, faz com
que ele não seja adquirido na leitura. Esse processo ocorre numa fase analógica da
medição. Em um segundo momento, o processamento digital torna possível multiplicar
o sinal por um segundo ganho θd, de modo que a série temporal recupere sua
amplitude original sem correr o risco de recuperar o ruído.
/ /0 (22)
Essa abordagem é válida apenas no caso do ruído ser de amplitude muito menor
que o sinal observado. Caso o sinal de interesse apresente alguma oscilação da
mesma magnitude do ruído, há perda de informação.
Outra abordagem a ser investigada é a filtragem linear no domínio da freqüência.
Aplicando-se a transformada de Fourier à série temporal s(n) obtém-se a descrição
̃ 4 desses elementos no domínio da freqüência, onde N é número de amostras
consideradas.
̃ 4 1
!
7 8 9:;⁄!
(23)
√6 -
Assim, aplicando-se a filtragem linear no domínio da freqüência, com o objetivo de
obter os valores estimados da série no domínio da freqüência =.> 4, a Eq. 24 é
obtida.
62
=>.4 /; ̃ 4
(24)
Calculando-se o erro de estimação e, entre a saída do sistema e a série temporal
medida, obtém-se:
@
8 ∑ A
@
;-B
A
%̃? 4 ̃ 4%
@
∑A
%/; =>4 4 =>4%
@
;-B
A
(25)
!/
∑;-B!//; 1 |=>4| /; |4|
Na busca do filtro que minimize o erro de estimação, deriva-se o erro em relação a
θk, igualando-o a zero. Assim, encontra-se o valor dos parâmetros do filtro linear ótimo
para estimar o estado real no domínio da freqüência.
O filtro ótimo dado pela Eq. 25 pode promover uma melhora na qualidade do sinal,
diminuindo a quantidade de ruído presente caso o estado X(n) e o ruído ξ(n) possuam
bandas de freqüência diferentes. No caso de ruídos de medição, é comum que as
faixas de freqüência sejam largas, mas uma vez que a dinâmica do sistema se
desenvolva sobre uma faixa de frequências relativamente estreita, é ainda possível
diminuir o ruído.
No entanto, quando o sistema é caótico, sua faixa de frequências abrange todo o
espectro, o que torna impossível tentar extrair o ruído. Portanto, fica claro que o uso
de filtragem linear é ineficiente na extração do ruído quando se trata de sistemas
caóticos. Portanto, a redução de ruído não-linear não pode depender da informação de
freqüência a fim fazer a separação entre o sinal e o ruído.
63
5.2. REDUÇÃO NÃO -LINEAR DE RUÍDO
A filtragem de sinais oriundos de sistemas não-lineares exige o uso de métodos
especiais já que as estruturas usuais de filtros espectrais ou lineares podem interagir
desfavoravelmente com a estrutura não-linear.
Dentre as várias abordagens para se fazer a filtragem de séries temporais
oriundas de sistemas não-lineares, uma delas é por meio da estimação do valor do
sinal real. Muitas vezes a estimação é feita por meio de algoritmos de predição, o que
nada mais é que obter-se uma estimativa do valor real do sinal em um instante futuro a
partir de seus valores passados. No intuito de fazer uma filtragem, o valor a ser
estimado é no instante atual.
Em um esquema de filtragem onde simplesmente substitui-se o sinal por uma
predição sua, ela já estaria contaminada, pois foi baseada em valores passados
contaminados. Se o sistema em questão for caótico, o erro, na média, tende a
aumentar, causando uma amplificação do ruído.
Portanto um segundo ingrediente é necessário para que a redução de ruído seja
bem sucedida. Uma alternativa é o ajuste à trajetória estimada da dinâmica do
sistema.
Diferentes métodos de redução de ruído foram propostos na literatura, diferindo no
modo em que a dinâmica é aproximada. Neste trabalho explora-se a estrutura no
espaço de fase reconstruído.
O algoritmo de redução não-linear de ruído mais simples substitui a coordenada
central de cada vetor de imersão pela média local desta coordenada. Isso atinge uma
aproximação localmente constante da dinâmica, sendo baseado na suposição que a
dinâmica é contínua.
Este esquema da redução de ruído é executado de modo simples, como mostra
Eq. 26. Primeiramente uma dimensão de imersão De e um passo de reconstrução têm
que ser escolhidos, para construir o espaço de imersão. Para cada ponto deste
64
espaço u(n) = { s(n-(De-1)), s(n-(De-2)),..., s(n) } são formadas vizinhanças Vn,
contendo todos os pontos u(p) tais que | u(p) - u(n) | < ζ. Para cada ponto u(p) da
vizinhança é tomada a coordenada de posição média s(p-De/2), sendo computada sua
média sobre a vizinhança Vn :
̂ EF /2 1
7 I EF /2
|H |
'J+KL
(26)
Deve-se ter o cuidado de fazer o raio ζ das vizinhanças, grande o bastante para
englobar as perturbações causadas pelo ruído, mas ao mesmo tempo deve ser menor
do que um típico raio de curvatura presente no sinal. Em um sinal qualquer, o raio de
curvatura quantifica as oscilações locais, descrevendo numericamente a magnitude
destas oscilações.
Se as oscilações induzidas localmente pelo ruído tiverem magnitude maior, ou
mesmo próxima, à magnitude das oscilações locais do sinal não ruidoso, a verdadeira
natureza da dinâmica do sinal pode ficar distorcida. Assim, para o filtro proposto em
Eq. 26, é necessário que sua área de atuação seja grande o bastante para abranger
algumas oscilações induzidas pelo ruído, mas muito menor que a área correspondente
a uma oscilação devida à dinâmica do sistema.
Estas restrições podem não ser sempre satisfeitas simultaneamente, neste caso
se tem que repetir o processo com diferentes escolhas de raio da vizinhança e
parâmetros de imersão, avaliando-se os resultados. Se o nível de ruído for
substancialmente menor do que os raios de curvatura presentes no sinal, isto é, o raio
das vizinhanças for cerca de 2 a 3 vezes maior que a amplitude do ruído, melhores
resultados são obtidos (Schreiber, 1993).
Depois de uma varredura completa sobre a série temporal todas as medidas são
substituídas pelos valores corrigidos. Naturalmente, para o primeiro e último pontos,
nenhuma correção está disponível. A correção média pode ser tomada como um raio
65
para a nova vizinhança na iteração seguinte. Deve-se estar atento para o detalhe de
que a vizinhança de cada ponto contém pelo menos o próprio ponto. Se ele é o único
membro da vizinhança, a filtragem fornece simplesmente a medida não corrigida e
nenhuma mudança é feita. Assim, podem-se executar iterações múltiplas, diminuindo
os valores de raio de vizinhança até que nenhuma mudança adicional seja feita.
5.3. ANÁLISE DE ECGS UTILIZANDO REDUÇÃO NÃO -LINEAR DE RUÍDO
Esta seção considera a análise de ECGs experimentais considerando uma
redução não-linear de ruído. Considera-se a aplicação do filtro definido na Eq. 26.
Inicialmente considera-se um ECG normal. A Figura 35 apresenta o sinal experimental
junto com o sinal filtrado que possui uma curva mais suave.
4
Ruidoso
Filtrado
3
ECG (mV)
2
1
0
-1
-2
10
10.5
11
11.5
12
12.5
Tempo (s)
13
13.5
14
14.5
15
Figura 35 – ECG normal experimental e sua filtragem.
A aplicação do procedimento de filtragem exige que se conheça os parâmetros de
imersão do sistema que são determinados através da informação mútua média, que
define um tempo de defasagem igual a 0,048s (12 amostras) e do método dos falsos
66
vizinhos, que define uma dimensão de imersão igual a 6, conforme mostrado na
Figura 36, onde é trabalhado um sinal ruidoso.
0.8
0.7
Informação Mútua
0.6
0.5
0.4
0.3
0.2
0.1
0
0
20
40
60
80
100
120
140
160
180
200
Defasagem (amostras)
(a)
0.8
Fracao de Falsos Vizinhos
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
2
4
6
8
10
12
14
16
18
Dimensao de Imersao (De)
20
(b)
Figura 36 – Análise dos parâmetros de imersão do ECG normal experimental:
Informação mútua (a) e falsos vizinhos (b)
Agora, considera-se o cálculo do maior expoente de Lyapunov correspondente à
série, através do algoritmo de Kantz. A dificuldade de analisar o expoente a partir do
sinal sem filtragem pode ser verificada nas Figuras 31 e 34.
67
0
Logarítmo do Fator de Deformação
-1
-2
-3
-4
e
e
e
e
e
-5
=
=
=
=
=
5.32x10E-04
6.87x10E-04
8.88x10E-04
1.14x10E-03
1.48x10E-03
De=
De=
De=
De=
De=
3
3
3
3
3
-6
-7
-8
5
10
15
20
25
30
35
40
45
50
Número de Iterações
Figura 37 – Máx. expoente de Lyapunov para o ECG normal experimental filtrado.
A Figura 37 mostra as curvas utilizadas para se obter o expoente de Lyapunov do
sinal filtrado. Consideram-se diferentes valores para o raio ε da vizinhança em torno da
trajetória de referência. Note que é possível obter curvas razoavelmente bem
comportantadas, onde rapidamente os valores passam a seguir uma reta horizontal,
mostrando a periodicidade do sinal de eletrocardiograma normal.
A seguir, considera-se um ECG diagnosticado como uma fibrilação ventricular. A
Figura 38 mostra uma comparação entre o sinal original e o obtido após a aplicação do
filtro.
68
4
ECG Fibrilado
ECG Fibrilado Filtrado
ECG (mV)
3
2
1
0
-1
-2
54
55
56
57
58
59
Tempo (s)
Figura 38 – ECG experimental associado à fibrilação ventricular e sua filtragem.
Mais uma vez, a utilização do processo de filtragem necessita do conhecimento
dos parâmetros de imersão cuja análise está apresentada na Figura 39. A informação
mútua média indica um tempo de defasagem δ = 0,84s (210 amostras) e os falsos
vizinhos indicam uma dimensão De = 20.
69
60
2.5
Informação Mútua
2
1.5
1
0.5
0
0
500
1000
1500
2000
2500
3000
Defasagem (amostras)
(a)
1
0.9
Fração de Falsos Vizinhos
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
2
4
6
8
10
12
14
Dimensão de Imersão(De)
16
18
20
(b)
Figura 39 – ECG experimental com fibrilação ventricular: informação mútua (a) e
falsos vizinhos (b).
70
-2
Logarítmo do Fator de Deformação
-2.2
-2.4
-2.6
-2.8
e=
e=
e=
e=
e=
e=
e=
e=
e=
e=
-3
-3.2
-3.4
-3.6
1.18 E-001
1.41 E-001
1.41 E-001
1.41 E-001
1.68 E-001
1.68 E-001
1.68 E-001
1.68 E-001
2.00 E-001
2.00 E-001
De=
De=
De=
De=
De=
De=
De=
De=
De=
De=
4
5
6
7
6
8
9
10
7
8
-3.8
-4
5
10
15
20
25
30
35
40
45
50
Número de Iterações
Figura 40 – Máximo expoente de Lyapunov para o ECG experimental filtrado com
fibrilação ventricular.
A análise do maior expoente de Lyapunov a partir do sinal filtrado está
apresentada na Figura 40. Utilizam-se diferentes valores dos parâmetros de filtragem
que não alteram significativamente o resultado. Deve-se obervar que existe uma
inclinação correspondente a um expoente λ entre 0,06 e 0,026, o que comparado ao
valor de 0,023 obtido para o máximo expoente, a partir de um sinal simulado de ECG
fibrilado, mostra que os valores são coerentes.
A interpretação do expoente de Lyapunov gerados a partir do algortimo de Kantz
continua sendo difícil. Deve-se notar que a partir de certo tamanho de janela de tempo
os afastamentos estabilizam em torno de um valor, levando a uma reta horizontal. Isto
ocorre devido à saturação das distâncias, pois a janela passou a ser tão grande, que o
afastamento médio passou a ser do tamanho médio do atrator. Após a filtragem, as
curvas de logaritmo do afastamento médio ficaram mais definidas, mostrando mais
claramente o trecho linear, o que facilita o cálculo do expoente máximo de Lyapunov.
No entanto, deve-se tomar cuidado para que a filtragem não descaraterize a dinâmica
do sistema.
71
6. CONCLUSÃO
Este trabalho apresenta a análise dos ritmos cardíacos a partir de duas
abordagens
diferentes.
Inicialmente,
apresenta-se
um
modelo
matemático
considerando-se três osciladores de Van der Pol modificados, conectados por
acoplamentos com defasagem. Cada oscilador representa um marca-passo natural do
coração: nódulo sinoatrial (SA), nódulo atrioventricular (AV) e complexo de HisPurkinje (HP). Desta forma, o modelo matemático se apresenta sob a forma de um
conjunto de Equações Diferenciais com Defasagem – EDD.
O método dos múltiplos passos é utilizado para solução das EDDs. É construida
uma função histórico para o sistema com defasagem por meio da expansão em série
de Taylor. Esta função irá fornecer o histórico temporal do sistema antes do intervalo
de simulação, podendo-se a partir daí empregar o método de Runge-Kutta de quarta
ordem para fazer a integração numérica das equações do modelo.
Simulações numéricas são realizadas mostrando que o modelo proposto é capaz
de capturar o comportamento geral da dinâmica dos batimentos cardíacos,
representando a forma do ECG normal com as curvas P, QRS e T. A seguir, alguns
ritmos patológicos
são
estudados estabelecendo-se diferentes
situações de
acoplamento. Basicamente, consideram-se interrupções de comunicação no sistema
elétrico cardíaco. Sob tais condições, é possível simular cardiopatologias como flutter
ventricular e bradicardia sinusal. Além disso, considera-se a excitação de marcapassos externos, representando outros comportamentos patológicos como fibrilação
ventricular. Esses resultados mostram que o modelo de três osciladores captura o
comportamento geral dos ritmos cardíacos representados pelo sinal de ECG e podem
levar a identificação de diferentes ritmos patológicos. Além disso, mostra-se que a
abordagem dinâmica pode ser útil para identificar características do sistema, o que
pode ser destacado através da análise do espaço de fase e de mapas de Poincaré.
72
A segunda abordagem utilizada neste trabalho é a análise de séries temporais de
ritmos cardíacos representados pelo sinal de ECG. Inicialmente, consideram-se sinais
gerados a partir do modelo matemático proposto. A análise utiliza a reconstrução do
espaço de estado e o cálculo dos expoentes de Lyapunov. A reconstrução do espaço
de estado emprega o método das coordenadas defasadas utilizando a informação
mútua média e o método dos falsos vizinhos para determinar os parâmetros de
imersão. Os expoentes de Lyapunov são avaliados a partir do algoritmo de Kantz.
Consideram-se dois tipos de ECGs: normais e com fibrilação ventricular. A análise
confirma que o ECG normal está relacionado a um comportamento periódico,
enquanto a fibrilação ventricular está relacionada ao caos. Depois disso, as técnicas
são aplicadas em sinais experimentais. A inevitável presença de ruído nos sinais
experimentais torna a análise mais difícil o que motiva a utilização de processos de
filtragem. Emprega-se então uma técnica não-linear de redução de ruído, onde é feita
uma suavização local do sinal, semelhante a uma média móvel. Esta filtragem causa
alguma melhora nas curvas empregadas no cálculo do máximo expoente de
Lyapunov, no entanto, ainda é difícil identificar um sinal onde não se conhece os
detalhes de suas características.
6.1. TRABALHOS FUTUROS
Vislumbra-se como continuação deste trabalho o desenvolvimento de um estudo
das órbitas periódicas instáveis (OPIs) presentes em sinais de ECG. Tais órbitas,
imersas em atratores caóticos, são fundamentais para o entendimento da dinâmica
caótica associada ao sistema. Uma aplicação particularmente importante destas
órbitas é o controle de sistemas caóticos em que o primeiro passo essencial
geralmente é a determinação de órbitas periódicas.
Sugere-se aplicar técnicas de detecção de órbitas instáveis aos sinais de ECG, de
modo a se atingir uma maior compreensão das características do sistema. O pior
cenário para a detecção de tais órbitas é a situação em que se deseja extrair OPIs a
73
partir de uma série temporal obtida experimentalmente. Isso advém do fato de que, em
geral, esse tipo de série apresenta ruído observacional (ruído branco).
O modelo matemático e os sinais obtidos a partir de sua integração numérica
serão utilizados para calibrar os procedimentos utilizados. Essencialmente, pretendese inserir ruído branco em séries geradas computacionalmente e, desta forma
estabelecer uma comparação com a série entre os sinais ideal e com ruído. A
existência de uma infinidade de OPIs em uma resposta caótica é a base para as
técnicas de controle de caos. Dessa forma, pode-se utilizar essas técnicas para
estabilizar a dinâmica do coração em uma órbita instável, não caótica, mais próxima
do comportamento normal. Isto corresponderia a retirá-lo de uma situação de caos,
correspondente possivelmente a uma fibrilação, trazendo para uma órbita periódica
estável, ou seja, batimentos cardíacos normais.
74
7. REFERÊNCIAS
Abarbanel, H. D. I., 1996, “Analysis of observed chaotic data”. Springer.
Babloyantz A., Destexhe, A., 1988, “Is the normal heart a periodic oscillator?”,
Biological Cybernetics, v. 58, pp. 203-211.
Boucekkine, R., Licandro, O. & Paul, C., 1997, “Differential-difference equations in
economics: On the numerical solution of vintage capital growth models”, Journal of
Economic Dynamics and Control, v. 21, pp. 347-362.
Campbell, S. R. & Wang, D., 1998, “Relaxation oscillators with time delay coupling”,
Physica D, v.111, pp.151-178.
Carvajal, R., Vallverdú, M., Baranowski, R., Chojnowska, L., Rydlewska Sadowska, W., Jané, R., Caminal, P., 1996, “Nonlinear analysis of heart rate variability
in patients with hiper-trophic cardiomyopathy”. Proceedings of the. 18th annual
conference IEEE-EMBS, CD-ROM paper 278.
Carvajal, R., Vallverdú, M., Zebrowski, J. J., Baranowski, R., Caminal, P., 1997, ”
Clasificación de pacientes con cardiomiopatía hipertrófica utilizando la dimensión de
correlación de la variabilidad del ritmo cardiaco”, Memorias del XV Congreso de la
SEIB, pp. 226-229.
Cunningham, W. J., 1954, “A nonlinear differential-difference equation of growth”
Proceedings of the National Academy of Sciences of the United States of America,
v.40, n.8, pp.708-713.
75
de Paula, A.S. & Savi, M.A., 2009, “A multiparameter chaos control method based
on
OGY approach”,
Chaos,
Solitons
&
Fractals,
v.40, n.3,
pp.1376-1390.
doi:10.1016/j.chaos.2007.09.056
Dubin, D., 1996, “Interpretação rápida do ECG”, Editora de Publições Biomédicas –
EPUB, Rio de Janeiro.
Einthoven, W. (1902), “Ein neus galvanometer”, Ann. Physik., v.12, pp.1059.
Ferriere, R. & Fox, G. A., 1995, “Chaos and evolution”, Trends in Ecology and
Evolution, v.10, n.12, pp.480-485.
Franca, L.F.P. & Savi, M.A., 2003, “Evaluating noise sensitivity on the time series
determination of Lyapunov exponents applied to nonlinear pendulum”, Shock and
Vibration, v.10, n.1, pp.37-50.
Fraser, A. M. & Swinney, H. L., 1986, “Independent coordinates for strange
attractors from mutual information”, Physical Review A, v.33, pp.1134–1140.
Garfinkel, A., Spano, M. L., Ditto, W. L. & Weiss, J. N., 1992, “Controlling cardiac
chaos”, Science, v.257, pp.1230-1235.
Garfinkel, A., Weiss, J. N., Ditto, W. L. & Spano, M. L., 1995, “Chaos control of
cardiac arrhythmias”, Trends in Cardiovascular Medicine, v.5, n.2, pp.76-80.
Glass L., Kaplan D., 1993, “Time series analysis of complex dynamics in physiology
and medicine”, Med. Progress through Technol., v. 19, pp. 115-128.
76
Glass, L., 2001, “Synchronization and rhythmic processes in physiology”, Nature,
v.410, pp.277-284.
Gois, S.R.F.S.M. & Savi, M.A., 2009, “An analysis of heart rhythm dynamics using
a three coupled oscillator model”, Chaos, Solitons & Fractals, v.41, n.5, pp. 2553-2565.
doi: 10.1016/j.chaos.2008.09.040.
Goldberger, A. L., 1992, “Fractal mechanisms in the electrophysiology of the heart”,
Transactions of the IEEE on Engineering in Medicine and Biology, pp. 47-52.
Goldberger, A. L., Rigney, D. R., West, B.J., 1990, “Chaos and fractals in human
physiology”, Scientific American, v. 262, pp. 42-49.
Grudzinski, K. & Zebrowski, J.J., 2004, “Modeling cardiac pacemakers with
relaxation oscillators”, Physica A, v.336, pp.153-162.
Guyton, A. C. & Hall, J.D., 1997, “Tratado de fisiologia médica”, 9 ed., Editora
Guanabara Koogan.
Hegger, R. Kantz H. & Schreiber T., 1999, “Practical implementation of nonlinear
time
series
methods:
The
TISEAN
package”,
Chaos,
v.
9,
pp.413.
DOI:10.1063/1.166424.
Hodkin, A. L. & Huxley, A. F., 1952, “A quantitative description of membrane
current and its application to conduction and excitation in nerve”, Journal of Physiology,
v.117, pp.500-544.
Jenkins, D. & Gerred S. “ECG library”, available in: <http://www.ecglibrary.com>.
77
Kantz, H., 1994, “A robust method to estimate the maximal Lyapunov exponent of a
time series”, Physics Letters A, v.185, pp. 77–87.
Koelliker, A & Mueller, H., 1855, “Nachweis der negativen schwankung des
muskelstromes am naturlich sich contrahieden muskel verhandel.” Journal Phys Med
Gesellsch., v. 6, pp.528.
Lombardi F., Sandrone G., Mortara A., Torzillo D., La Rovere M.T., Signorini M.G.,
Cerutti S., Malliani A., 1996, “Linear and nonlinear dynamics of heart rate variability
after acute myocardial infarction with normal and reduced left ventricular ejection
fraction”. J. Cardiol., v. 77, pp. 1283-1288.
Mackey, M. C., Glass, L., 1977, “Oscillation and chaos in physiological system”,
Science, v. 197, pp. 287-289.
Malmivuo, J. & Plonsey, R., 1995, “Bioelectromagnetism - principles and
applications of bioelectric and biomagnetic fields”, Oxford University Press, New York.
Moffa, P. J. & Sanches, P. C. R., 2001, “Eletrocardiograma normal e patológico”,
Editora Roca.
Pereira-Pinto, F.H.I., Ferreira, A.M & Savi, M.A., 2004, “Chaos control in a
nonlinear pendulum using a semi-continuous method”, Chaos, Solitons & Fractals,
v.22, n.3, pp.653-668.
Radhakrishna, R. K. A., Dutt, D. N., Yeragani, V. K., 2000, “Nonlinear measures of
heart rate time series: Influence of posture and controlled breathing”, Autonomic
Neuroscience: Basic and Clinical, v.83, pp.148–158.
78
Ritzenberg, A. L., Adam, D. R., Cohen, R. J., 1984, “Period multiplying evidence for
nonlinear behavior of the canine heart”, Nature, v.114, pp. 321-327.
Rosenstein, M. T., Collins, J. J., and De Luca, C. J., 1993, “A practical method for
calculating largest Lyapunov exponents from small data sets”, Physica D, v.65, pp.
117–134.
Sano, M. and Sawada, Y., 1985, “Measurement of the Lyapunov spectrum from a
chaotic time series”, Physical Review Letters, v.55, pp. 1082–1085.
Santos, A. M., Lopes, S. R. & Viana, R. L., 2004, “Rhythm synchronization and
chaotic modulation of coupled Van der Pol oscillators in a model for the heartbeat”,
Physica A, v.338, pp.335-355.
Savi, M. A., 2005, “Chaos and order in biomedical rhythms”, Journal of the Brazilian
Society of Mechanical Sciences and Engineering, v.27, n.2, pp.157-169.
Savi, M. A., 2006, “Dinâmica não linear e caos”, Editora E-papers.
Savi, M.A., Pereira-Pinto, F.H.I. & Ferreira, A.M., 2006, “Chaos control in
mechanical systems”, Shock and Vibration, v.13, n.4/5, pp.301-314.
Schreiber, T., 1993, “Extremely simple nonlinear noise reduction method”, Phys.
Rev. E, v.47, pp.2401.
Signorini M.G., Cerutti S., 1994, “Lyapunov exponents calculated from heart rate
variability time series”, IEEE-EMBS, pp. 119-120.
79
Takens, F., 1981, “Detecting strange attractors in turbulence”', Lecture Notes in
Math., v.898, Springer, New York.
Van Der Pol, B. & Van Der Mark, J., 1928, “The heartbeat considered as a
relaxation oscillator and an electrical model of the heart”, Philosophical Magazine, v.6 ,
pp. 763.
Van Der Pol, B., 1926, “On relaxation oscillations”, Philosophical Magazine, v.2,
pp.978.
Voss A., Kurths J., Kleiner H.J., Wessel N., 1995, “Improved analysis of heart rate
variability by methods of nonlinear dynamics”. J. Electrocardiol, vol. 28, pp. 81-88.
Voss A., Kurths J., Kleiner H.J., Witt A., Saparin P., Dietz R., Fiehring H., Wessel
N., 1994, “Nonlinear dynamics versus traditional methods of heart rate variability
analysis”, Computers in Cardiology, pp. 125-128.
Waller, A. D., 1888, “Introductory address on the electromotive properties of the
human heart “, Br. Med J. v.2, pp.751–754.
Winfree, A. T., 1980, “The geometry of biological time”, Springer.
Witkowski, F. X., Kavanagh, K. M., Penkoske, P. A., Plonsey, R., Spano, M. L.,
Ditto, W. L. & Kaplan, D. T., 1995, “Evidence for determinism in ventricular fibrillation”,
Physical Review Letters, v.75, n.6, pp.1230-1233.
80
Witkowski, F. X., Leon, L. J., Penkoske, P. A., Giles, W. R., Spano, M. L., Ditto,
W. L. & Winfree, A. T., 1998, “Spatiotemporal evolution of ventricular fibrillation”,
Nature, v.392, pp.78-82.
Wolf, A., Swift, J. B., Swinney, H. L. & Vastano, J. A., 1985, “Determining Lyapunov
exponents from a time series”, Physica D, v.16, pp.285-317.
Yanowitz,
F.G.
“ECG
learning
center
in
cyberspace”,
available
in:
http://library.med.utah.edu/kw/ecg/image_index.
Zebrowski J.J., Buchner T., Baranowski R., Poplawska W., 1996, “Nonlinear
analysis of heart dynamics”, Med. Biol. Eng. Comput., v. 34, pp. 379-380.
Zebrowski J.J., Poplawska W., Baranowski R., Buchner T., 1997, “Nonlinear
analysis, theory, methods and applications”, Proc. 2nd World Congress of Nonlinear
Analysts, v.30, pp. 1007-1017.
81