M Moodelização Matemática do Tremor Humano na Doença de … · As oscilações associadas a estes mecanismos de feedback, juntamente com as oscilações resultantes da própria

Modelização Matemática do Tremor Humano na Doença de Parkinson

Modelos dos Processos Fisiológicos no Homem 1

MMooddeelliizzaaççããoo MMaatteemmááttiiccaa ddoo TTrreemmoorr HHuummaannoo nnaa

DDooeennççaa ddee PPaarrkkiinnssoonn

ffeerrrreeiirraaeeddggaarr22000011@@hhoottmmaaiill..ccoomm EEddggaarrFFeerrrreeiirraa

JJooããoo DDuuaarrtteeRRiiccaarrddoo MMaarrttiinnss

TTiiaaggoo MMaarrççaall

EEnnggeennhhaarriiaa BBiioommééddiiccaa

MMooddeellooss ddooss PPrroocceessssooss FFiissiioollóóggiiccooss nnoo HHoommeemm



Índice Introdução …………………...……………………………….................................. 2 Caracterização da Doença de Parkinson ……..……………………...………….. 4

Perspectivas terapêuticas ……………………………………………….. 4 Modelos matemáticos do tremor Parkinsoniano ………………………………... 5

Análise do comportamento das unidades da rede neuronal quando estas têm uma resposta de acordo com o degrau de Heaviside ……………….

8

Simulação da rede neuronal para 1α = …………………………. 8 Simulação da rede neuronal para 0,8α = ………………………. 10 Simulação da rede neuronal para 0,5α = ………………………. 12

Análise do comportamento das unidades da rede neuronal quando estas têm uma resposta de acordo com uma função sigmóide………………...

14

Simulação da rede neuronal para 1α = …………………………. 14 Simulação da rede neuronal para 0,8α = ………………………. 16 Simulação da rede neuronal para 0,5α = ………………………. 18

Análise da Dinâmica da Rede Neuronal ………………………………………… 20 Conclusão …………………………………………………………………………. 28 Anexo……………………………………………………………………….............. 29 Bibliografia ………………………………………………………………………... 31



Introdução

Nos últimos anos tem-se vindo a verificar um interesse cada vez maior na actividade neuronal oscilatória do sistema nervoso central (SNC) e na importância que esta actividade pode ter no controlo motor.

O tremor é definido como sendo um movimento rápido, aproximadamente rítmico, de determinada parte do corpo, verificando-se normalmente nas extremidades do corpo.

Foi apenas por volta do século XIX que Horeley e Schofer demonstraram que o tremor teria uma origem neurológica.

Inicialmente pensava-se que os mecanismos associados à génese do tremor eram

resultado da contracção independente de cada músculo. Por outro lado, também era sugerido que o tremor se devesse à tendência natural de certas partes do corpo para oscilar a certas frequências. No entanto, estudos posteriores baseados na análise de electroencefalogramas realizados em animais anestesiados, sugeriram que o tremor se deve à ritmicidade da actividade neuronal do sistema nervoso central.

Devido à aplicação de novas técnicas como a magnetoencefalografia, apenas recentemente se conseguiu mostrar que a actividade oscilatória rítmica verificada no SNC se manifesta na periferia. Tem sido proposto que esta actividade tem um contributo importante no processamento de informação sensorial (Eckhorn et al) e de comandos motores (Farmer, 1998).

Pode classificar-se o tremor tendo em conta o estado de actividade da parte do corpo na qual se observa o tremor:

• Tremor de repouso: ocorre quando os músculos principais não estão activos, e a parte do corpo em estudo está completamente apoiada de forma a contrariar a gravidade.

• Tremor de acção: ocorre quando os músculos principais estão em actividade. No que diz respeito à causa do tremor, este pode ser classificado como tremor

normal ou patológico. O tremor normal (fisiológico), que resulta da actividade neuronal rítmica do SNC, é irregular, apresentando um espectro de frequências difuso e extenso.

As oscilações devidas aos mecanismos de feedback resultam de uma instabilidade inerente a estes ciclos (de feedback negativo). Um mecanismo de feedback negativo provoca, por si só, um processo oscilatório, dado que, com uma variação do sistema, este mecanismo provoca uma resposta de forma a contrariar essa mesma variação.

Existem três mecanismos de feedback responsáveis pelo controlo preciso do movimento, e que poderão induzir o tremor. Dois destes circuitos envolvem o transporte de informação pelas vias aferentes a partir de receptores sensíveis ao estiramento. Estas vias provêm dos músculos e vão fornecer informação acerca do estado dos mesmos, numa primeira fase, aos neurónios motores da medula (a resposta automática, ocorre sem pensamento consciente). Numa segunda fase essa informação é transmitida pelas vias aferentes às áreas do córtex responsáveis pelo controlo do movimento.



O terceiro ciclo envolve um mecanismo de feedback que transmite informação directamente dos neurónios para as áreas do córtex responsáveis pelo controlo do movimento. Este processo envolve apenas estruturas do SNC.

As oscilações associadas a estes mecanismos de feedback, juntamente com as oscilações resultantes da própria ressonância das estruturas mecânicas envolvidas no movimento (a ressonância resulta de quando essas estruturas mecânicas são perturbadas por forças externas ou internas), produzem um tremor denominado de tremor fisiológico.

Existem, igualmente, diversas circunstâncias extraordinárias (pontuais) que

podem induzir o tremor, tais como: estar sujeito a determinado tipo de medicação, situações de descontrolo emocional, ansiedade, medo, stress, raiva, fadiga física, hipotermia, dor, entre outros. Contudo, não é claro que estes medicamentos induzam alterações idênticas, quer no domínio do tempo, quer no domínio da frequência do tremor.

O tremor patológico resulta de uma distorção ou amplificação da actividade oscilatória do SNC e caracteriza-se por ser mais regular, apresentando por isso um espectro de frequências com picos bem definidos (mais concentrados).

Assim, de entre as várias doenças às quais está associado o tremor patológico (tais como a coreia de Sydenham, paralisia cerebral, doença de Parkinson), visaremos neste trabalho a modelização matemática de uma rede neuronal (e o consequente tremor por esta provocado), envolvida no controlo da actividade motora, e das alterações na sua actividade resultantes do evoluir da doença de Parkinson num indivíduo.

CÓRTEX CEREBRAL

(Estriado, substância nigra,

núcleos sub-talâmicos)

FEEDBACK CENTRAL

CONJUNTO DE

NEURÓNIOS MOTORES

RECEPTORES DE ESTIRAMENTO

MÚSCULO

MOVIMENTO (Tremor fisiológico)

Propriedades da unidade motora (ritmos não sincronizados e sincronização não oscilatória)

Factores mecânicos

Oscilações de grande latência devido aos mecanismos de feedback

Oscilações de baixa latência devido aos mecanismos de feedback

Oscilações resultantes da actividade do SNC

Oscilações devidas ao mecanismo de feedback central

Figura 1 Diagrama ilustrativo dos vários pontos ao longo da via motora onde o tremor fisiológico, podeter origem.



Figura 2 Representação dos núcleos da base do hemisfério esquerdo.

Caracterização da Doença de Parkinson

A doença de Parkinson é uma doença neurodegenerativa que aparece, normalmente, depois dos 55 anos de idade (não é contagiosa nem hereditária), causada pela perda de neurónios dopaminérgicos, existentes na substância nigra. Consequentemente, existe uma diminuição da produção de dopamina (neurotransmissor) nesta região, afectando de forma directa a actividade do estriado.

O estriado e os restantes núcleos da base (núcleos subtalâmicos, substância nigra – ver Figura 2) são importantes na organização e coordenação dos movimentos e da postura, bem como na diminuição do tónus muscular e na inibição da actividade muscular indesejável. Assim sendo, as perturbações ocorridas na substância nigra, causadoras da doença de Parkinson, provocam uma diminuição da capacidade de controlar os movimentos e de se movimentar (bradicinésia), os movimentos tornam-se mais rígidos (bruscos), ocorre uma diminuição da capacidade para manter a postura e há perda de expressão facial. Surge nestes doentes um tremor característico, que consiste em movimentos circulares e transversais das pontas do polegar e do indicador. A maioria dos pacientes com doença de Parkinson apresenta um tremor de repouso com uma frequência entre os 4 e 6 Hz, o qual é inibido pela activação voluntária dos músculos e aumentado quando o indivíduo se movimenta ou está distraído. Estes indivíduos (cerca de 60%) apresentam também um tremor postural (tremor de acção associado à manutenção da postura) que tem uma frequência entre os 5 e os 12 Hz e que pode ser observado mesmo quando não se verifica tremor de repouso. Perspectivas terapêuticas No sentido de minorar os efeitos da doença de Parkinson diversas estratégias têm sido desenvolvidas em variadas áreas científicas.

Numa vertente farmacológica tem tentado diminuir-se os sintomas da doença de Parkinson através da utilização da Levodopa (L-dopa), precursora da dopamina, ou, ainda mais eficazmente, pela combinação de L-dopa e carbidopa (a qual impede que a L-dopa seja absorvida por tecidos que não sejam os encefálicos).

Foi descoberta uma proteína, o factor neurotrópico derivado da linha de células gliais (GDNF), que promove a sobrevivência selectiva dos neurónios secretores de dopamina.

Numa perspectiva electrofisiológica, a estimulação eléctrica contínua do núcleo intermédio ventral a uma frequência superior a 100 Hz pode diminuir o tremor. Esta mesma estimulação, quando aplicada ao globus pallidus interno e aos núcleos subtalâmicos, tem como efeito o controlo do tremor e a redução dos sintomas relacionados com perturbações do movimento.

Numa perspectiva médico-cirúrgica, os transplantes de células da medula adrenal ou de células da substância nigra de fetos têm vindo a ser testados em doentes

com Parkinson como potenciais fontes de dopamina, de forma a colmatar o decréscimo da produção de dopamina nestes indivíduos.



Modelos matemáticos do tremor Parkinsoniano Ao longo do tempo diversos investigadores propuseram vários modelos que

pretendiam ser representativos do tremor Parkinsoniano. Contudo, embora esses modelos reproduzissem de forma bastante aproximada determinadas características do tremor, o mesmo não acontecia com outras características, sendo por isso modelos incompletos. Desta forma, as tentativas realizadas para criar modelos representativos do tremor na doença de Parkinson não lhes permitiu determinar as características dinâmicas do funcionamento das redes neuronais envolvidas na doença.

Actualmente chegou-se a um consenso generalizado entre investigadores de que o tremor Parkinsoniano está associado a grupos de neurónios, localizados no córtex cerebral, que geram potenciais de acção sincronizados a frequências semelhantes à frequência do tremor. Investigações desenvolvidas por Edwards, Beuter e Glass levaram ao desenvolvimento de um modelo matemático representativo de uma rede neuronal inibitória (pois a dopamina é um neurotransmissor na maioria dos casos inibitórios), rica em conexões entre os diferentes grupos de neurónios.

Neste modelo é possível, como será explicado mais adiante, variar determinados parâmetros. Esta variação representa o decréscimo da eficácia sináptica dos neurónios dopaminérgicos da via nigro-estriada (diminuição da eficácia de transmissão de informação da substância nigra para o estriado – característica da doença de Parkinson). Neste modelo é possível simular a transição de um tremor fisiológico normal (dinâmica irregular) para um tremor Parkinsoniano (dinâmica regular) à medida que a eficácia sináptica diminui. Isto implica que neste modelo se considere que o tremor Parkinsoniano e o tremor fisiológico têm origem no mesmo circuito motor, embora resultem de diferentes comportamentos desse circuito.

Neste modelo considera-se uma unidade como sendo um grande conjunto de neurónios a actuar simultaneamente, e no qual a actividade de cada unidade afecta a actividade das restantes.

Assim, cada unidade desta rede neuronal constituída por N unidades pode ser representada pela seguinte equação:

1( )

Ni

i ij j ij

dy y w g ydt

τ=

= − + −∑

Fisiologicamente wij representa a eficácia sináptica da unidade j que actua em i,

enquanto iτ representa o limiar de potencial de acção da unidade i. A função g(yj) é a função resposta de determinada unidade j aos seus estímulos de entrada.

A função resposta g pode tomar as seguintes formas: • Função sigmóide:

1 tanh( )

( )2

jj

yg y

β+= , em que β é a inclinação (ganho) da função sigmóide.

• Função degrau de Heaviside:

1 , 0

( )0 , 0

jj

j

yg y

y≥= <



A equação apresentada é capaz de representar um número infinito de sistemas neuronais diferentes, pois podemos variar o número de unidades, bem como o tipo e quantidade de ligações entre estas.

Para efectuar o nosso estudo considera-se a seguinte rede neuronal constituída por seis unidades e em que as diferentes unidades se encontram relacionadas da seguinte forma:

Nesta representação da rede neuronal considerada, as ligações que não apresentam o parâmetro α representam ligações cuja eficácia sináptica não é alterada durante as simulações (sendo permanentemente igual a 1). Por sua vez, as ligações que possuem o parâmetro α têm a eficácia sináptica alterada durante a simulação. Dado que o parâmetro α é menor que 1, a eficácia sináptica diminui e este modelo representa a evolução da doença de Parkinson. Para as situações em que não se estabelece uma ligação directa entre duas unidades, considera-se nula a eficácia sináptica da ligação entre essas duas unidades. Deste modo obtém-se uma matriz W que contém as eficácias sinápticas, wij, das diversas ligações da rede neuronal apresentada:

0 1 0 0 0 10 0 0 1 0 10 0 0 1 01 0 0 0 0 11 1 0 0 0 00 0 0 0

Wα

α α

= −

As eficácias sinápticas consideram-se negativas, pois nas regiões do SNC que

este modelo pretende representar, a conexão entre a região da substância nigra e do estriado, a dopamina (principal neurotransmissor envolvido na transmissão de informação nesta região), tem um efeito predominantemente inibitório.

Pelo facto de as unidades se localizarem em regiões próximas entre si, considera-se que estas são constituídas por neurónios do mesmo tipo. Assim, toma-se o limiar de potencial de acção das diversas unidades como sendo igual, tendo o valor de

1,5− (como sugerido no livro Nonlinear Dynamics in Physiology and Medicine, Anne Beuter).

1

4

6

2

35α

αα

Figura 3 Diagrama ilustrativo da rede neuronal constituída por seis unidades.



As condições iniciais consideram-se aleatórias, tendo sido escolhidos os seguintes valores:

1

2

3

4

5

6

(0) 0,125(0) 0, 200(0) 0,019(0) 0, 234(0) 0,195(0) 0, 257

yyyyyy

= = − = − = − =

= −

Deste modo é possível escrever as seis equações diferenciais que descrevem o comportamento em cada uma das seis unidades ao longo do tempo:

11 12 2 16 6 1 1 2 6

22 24 4 26 6 2 2 4 6

33 34 4 35 5 3 3 4 5

44 41 1 46 6 4

( ) ( ) 1 ( ) 1 ( ) 1,5

( ) ( ) 1 ( ) 1 ( ) 1,5

( ) ( ) 1 ( ) ( ) 1,5

( ) ( )

dy y w g y w g y y g y g ydtdy y w g y w g y y g y g ydt

dy y w g y w g y y g y g ydtdy y w g y w g ydt

τ

τ

τ α

τ

= − − ⋅ − ⋅ − = − − ⋅ − ⋅ +

= − − ⋅ − ⋅ − = − − ⋅ − ⋅ +

= − − ⋅ − ⋅ − = − − ⋅ − ⋅ +

= − − ⋅ − ⋅ − 4 1 6

55 51 1 52 2 5 5 1 2

66 63 3 65 5 5 6 3 5

1 ( ) 1 ( ) 1,5

( ) ( ) 1 ( ) 1 ( ) 1,5

( ) ( ) ( ) ( ) 1,5

y g y g y

dy y w g y w g y y g y g ydtdy y w g y w g y y g y g ydt

τ

τ α α

= − − ⋅ − ⋅ +

= − − ⋅ − ⋅ − = − − ⋅ − ⋅ +

= − − ⋅ − ⋅ − = − − ⋅ − ⋅ +

Análise do comportamento das unidades da rede neuronal quando estas têm uma resposta de acordo com o degrau de Heaviside Como referido anteriormente, o degrau de Heaviside é definido como:

1 , 0( )

0 , 0j

jj

yg y

y≥= <

De forma a implementar o modelo da rede neuronal constituída por seis

unidades, para este caso, criou-se um subsistema representativo da função g(yj):

1Out1

>=

RelationalOperator Product

1Constant10

Constant

1In1

Figura 4 Diagrama do subsistema que simula a função g(yj).



Neste subsistema a entrada yj é comparada de forma contínua com a constante 0. Sempre que yj for maior ou igual que 0, isto é, sempre que a condição apresentada for verdadeira, o Relational Operator devolve o boleano “1”. Caso contrário, a condição apresentada é falsa, e o Relational Operator devolve o boleano “0”. De seguida efectua-se o produto do boleano obtido por 1, de forma a convertê-lo num número que possa ser manipulado algebricamente pelas operações posteriores no modelo. Fisiologicamente, o degrau de Heaviside indica que apenas quando a saída de determinada unidade da rede neuronal é superior ou igual a 0, esta actua nas unidades com as quais tem ligações estabelecidas. Tal indica que, nesta situação, nem sempre as unidades estão a exercer influência entre si. As unidades exercem influência apenas em determinados períodos de tempo.

• Simulação da rede neuronal para 1α = :

Ao tomar-se 1α = considera-se, neste modelo, que as eficácias sinápticas das ligações entre as unidades não se encontram afectadas. Deste modo, a actividade de determinada unidade é transmitida de forma inalterada às unidades com as quais tem ligações estabelecidas. Assim, considerando 1α = obtém-se o comportamento da rede neuronal que origina o tremor fisiológico normal.

De seguida, é apresentada a implementação desta rede neuronal, tendo em conta os parâmetros acima referidos, bem como a saída de uma das unidades da rede, a unidade 1, pois as restantes 5 unidades têm um comportamento que embora não sendo exactamente igual é equivalente:

In1Out1

g(y6)

In1Out1

g(y5)

In1Out1

g(y4)

In1Out1

g(y3)

In1Out1

g(y2)

In1Out1

g(y1)

y1

To Workspace1

In1

In2

Out1

Subsistema 6

In1

In2Out1

Subsistema 5

In1

In2

Out1

Subsistema 4

In1

In2

Out1

Subsistema 3

In1

In2Out1

Subsistema 2

In1

In2

Out1

Subsistema 1

Scope

1

Alfa

1 Alfa

1Out1

1s

Integrator

-1.5

Constant

2In2

1In1

Figura 5 Diagrama de blocos representativo da simulação da rede neuronal para α=1. O bloco Workspacecria um vector y1, no Matlab, contendo as coordenadas da simulação. A caixa ao lado representa os subsistemas 1 a 6. Nestes subsistemas, Constant é o bloco que representa o limiar do potencial de acção.



Analisando a resposta da unidade 1, obtida no Scope, verifica-se que esta

apresenta um comportamento bastante irregular. Desta forma, a irregularidade da resposta desta unidade, que se supõe estar envolvida no controlo e organização de determinadas funções motoras, reflecte-se no aparecimento de um tremor com um comportamento igualmente irregular. Este tremor irregular é o denominado tremor fisiológico. De forma a fazer uma análise mais completa da resposta da saída da unidade 1, determinou-se o seu espectro de frequências. Para tal recorreu-se a um algoritmo descrito no anexo 1, tendo-se obtido o espectro de frequência a seguir apresentado:

0 0.5 1 1.5 2 2.5 30

5

10

15

20

25

30

35

40Espectro de Frequências de y1

Frequência (Hz)

Am

plitu

de

Figura 7 Espectro de frequências da saída da unidade 1 para α=1.

0 5 10 15 20 25 30 35 40 45 50-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

Tempo (s)

Am

plitu

de

Saída da unidade 1

Figura 6 Gráfico da saída da unidade 1 da rede neuronal para α=1.



Ao analisar o espectro de frequências obtido verifica-se que a saída da unidade 1 é constituída por uma gama considerável de frequências, sendo 0,24Hz a frequência de maior amplitude no espectro. A gama de frequências obtidas reflecte o comportamento irregular da saída da unidade 1 visualizada no Scope.

• Simulação da rede neuronal para 0,8α = :

Ao considerar-se 0,8α = , as eficácias sinápticas de determinadas ligações entre unidades ficam afectadas. Assim, a actividade de determinada unidade é transmitida de forma alterada, menos eficaz, às unidades com as quais tem ligações estabelecidas.

De seguida, apresenta-se a implementação desta rede neuronal, e pelo motivo referido anteriormente, apresenta-se apenas a saída de uma das unidades da rede, a unidade 1:

In1Out1

g(y6)

In1Out1

g(y5)

In1Out1

g(y4)

In1Out1

g(y3)

In1Out1

g(y2)

In1Out1

g(y1)

y1

To Workspace1

In1

In2

Out1

Subsistema 6

In1

In2

Out1

Subsistema 5

In1

In2

Out1

Subsistema 4

In1

In2

Out1

Subsistema 3

In1

In2

Out1

Subsistema 2

In1

In2

Out1

Subsistema 1

Scope

0.8

Alfa

0.8 Alfa

Figura 8 Diagrama de blocos representativo da simulação da rede neuronal para α=0,8. O bloco Workspacecria um vector y1, no Matlab, contendo as coordenadas da simulação.



Analisando a saída da unidade 1, obtida através do Scope, verifica-se que o

gráfico ainda apresenta um comportamento bastante irregular, sendo no entanto, mais regular do que o obtido para 1α = .

Recorrendo ao algoritmo descrito no anexo 1 obteve-se para esta saída o seguinte espectro de frequências:

0 5 10 15 20 25 30 35 40 45 50-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

Tempo (s)

Am

plitu

de

Saída da unidade 1

Figura 9 Gráfico da saída da unidade 1 da rede neuronal para α=0,8.

0 0.5 1 1.5 2 2.5 30

2

4

6

8

10

12

14


Frequência (Hz)

Am

plitu

de

Figura 10 Espectro de frequências da saída da unidade 1 para α=0,8.



Neste gráfico, a gama de frequências envolvidas na saída da unidade 1 é ligeiramente menor do que em 1α = e apresenta um pico mais definido quando a frequência é aproximadamente 0,3Hz. Os restantes picos presentes no espectro têm, relativamente aos obtidos quando 1α = , menor amplitude.

• Simulação da rede neuronal para 0,5α = :

Quando se considera 0,5α = as eficácias sinápticas de determinadas ligações

entre unidades encontram-se reduzidas a metade. Para estas condições, a actividade de algumas unidades é transmitida a outras unidades de forma muito pouco eficaz, reflectindo-se à posteriori na actividade das restantes unidades. Isto verifica-se pelo facto de a actividade de qualquer das unidades ser dependente, de forma mais ou menos directa, da actividade das restantes.

De seguida é apresentada a implementação em Simulink desta rede neuronal:

Pela análise da saída da unidade 1, obtida através do Scope, verifica-se que esta apresenta um comportamento regular, que se repete ao longo do tempo. A regularidade da saída desta unidade, bem como das restantes unidades, reflecte-se no surgimento de um tremor com um comportamento igualmente regular, característicos de indivíduos com a doença de Parkinson.

Recorrendo ao algoritmo descrito no anexo 1 obteve-se para esta saída o seguinte espectro de frequências:

0 5 10 15 20 25 30 35 40 45 50-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4Saída da unidade 1

Tempo (s)

Am

plitu

de




Analisando o espectro de frequências obtido para a resposta da unidade 1

verifica-se que este apresenta um único pico bem definido para uma frequência de 0,34Hz, devido ao comportamento periódico da resposta da unidade 1. Com a análise dos gráficos obtidos, verifica-se que a diminuição do valor do parâmetro α origina uma evolução da saída da unidade 1 (e, consequentemente, do tremor) de um comportamento irregular (aleatório, imprevisível) para um comportamento regular, o qual se apresenta periódico. O comportamento irregular inicialmente registado ( 1α = , situação fisiológica normal) deve-se em parte à própria natureza da função de saída das diversas unidades às suas entradas. A função de saída considerada neste caso é a função de Heaviside, a qual origina uma saída binária, do tipo “tudo ou nada”, ou seja, o seu valor apenas pode tomar os valores “0” ou “1”. Assim, em cada unidade neuronal tem-se uma variação brusca da resposta à sua entrada, pelo que todas as unidades do sistema serão afectadas desta variação, devido ao facto de esta variação ser transmitida a todas elas.

Fisiologicamente, este tipo de comportamento das unidades neuronais corresponde a uma tentativa permanente de ajustar a posição real de determinada parte do corpo à posição desejada. No entanto, neste caso, devido à natureza da função de saída de cada unidade, este ajuste é feito de modo grosseiro (tipo “tudo ou nada”), pelo que nunca se atinge uma posição de equilíbrio definitiva.

Para 0,5α = , e como referido anteriormente, verifica-se que as unidades apresentam um comportamento regular ao longo do tempo. Fisiologicamente, este comportamento implica que, embora possa haver uma predisposição por parte de um indivíduo para manter determinada parte do corpo numa posição fixa, tal não ocorre, uma vez que as ligações entre as diversas unidades têm uma eficácia sináptica bastante reduzida. Isto implica que as ligações entre as unidades envolvidas no controlo motor se encontrem enfraquecidas. Verifica-se permanentemente este tipo de comportamento por parte das unidades da rede neuronal, independentemente da vontade do indivíduo, provocando o tremor Parkinsoniano de forma contínua.

0 0.5 1 1.5 2 2.5 30

10

20

30

40

50


Frequência (Hz)

Am

plitu

de




(1+tanh(4*u))/2

Fcn

Análise do comportamento das unidades da rede neuronal quando estas têm uma resposta de acordo com uma função sigmóide

Como foi referido anteriormente, a função sigmóide é definida por:

1 tanh( )

( )2

jj

yg y

β+=

De forma a implementar em Simulimk a função sigmóide para 4β = (valor

referido na bibliografia), recorre-se a um bloco Fcn no qual se insere 1 tanh(4 )2

u+ e no

qual u é a entrada do bloco. Para melhor compreender os resultados obtidos com este modelo, apresenta-se a

representação gráfica da função sigmóide:

• Análise da rede neuronal para 1α = :

Como referido anteriormente, ao considerar-se 1α = , pretende-se simular a situação em que a eficácia sináptica das ligações entre as unidades não se encontra alterada (situação fisiológica normal).

A implementação deste modelo no Simulink recorre a uma estrutura semelhante à utilizada anteriormente, com excepção das funções g(yj). Através do Scope obtém-se a representação da actividade da unidade 1:

Figura 13 Representação gráfica da função sigmóide para 4β = .

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0

0.2

0.4

0.6

0.8

1

Função sigmóide para Beta = 4

Entrada

Saí

da



f(u)

g(y6)

f(u)

g(y5)

f(u)

g(y4)

f(u)

g(y3)

f(u)

g(y2)

f(u)

g(y1)

y1

To Workspace

In1

In2

Out1

Subsistema 6

In1

In2

Out1

Subsistema 5

In1

In2Out1

Subsistema 4

In1

In2

Out1

Subsistema 3

In1

In2

Out1

Subsistema 2

In1

In2

Out1

Subsistema 1

Scope

1

Alfa

1Alfa

Figura 14 Diagrama de blocos representativo da simulação da rede neuronal para α=1. O bloco Workspacecria um vector y1, no Matlab, contendo as coordenadas da simulação.

Figura 15 Gráfico da saída da unidade 1 da rede neuronal para α=1.

0 5 10 15 20 25 30 35 40 45 50-0.05

0

0.05

0.1

0.15

0.2

0.25

0.3


Tempo (s)

Am

plitu

de



Inicialmente a actividade da unidade 1 apresenta grandes oscilações, contudo, a sua amplitude tende a diminuir rapidamente. A partir de determinado instante, a unidade 1, bem como as restantes, atinge uma actividade praticamente constante (0,1046 para um tempo suficientemente grande). Este comportamento da unidade 1 deve-se ao facto de a função resposta considerada (função sigmóide) ser uma função contínua, o que faz com que a saída das diversas unidades não varie tão bruscamente como na situação em que se utilizou a função de Heaviside (para 1α = ). O facto de a variação não ser tão brusca permite que seja atingido um ponto a partir do qual a actividade (saída) das unidades é constante.

Ao aplicar o algoritmo da transformada de Fourier (Anexo 1) à saída da unidade 1, obtém-se o seguinte espectro de frequências:

Analisando o espectro de frequências obtido verifica-se que este apresenta

apenas dois picos de frequência significativos. Um destes picos, correspondente a 0Hz, tem grande amplitude pois, durante a maioria do tempo, a saída da unidade 1 é constante. O outro pico, correspondente a 0,24Hz, tem uma amplitude reduzida e resulta das oscilações iniciais da saída da unidade 1.

Do ponto de vista fisiológico, pode-se entender este tipo de comportamento das unidades como resultante de um controlo mais fino e mais preciso, por parte das componentes envolvidas, do movimento e da manutenção da posição de algumas partes do corpo (que sofreram uma variação da posição). A actividade das unidades varia cada vez menos, à medida que a posição de certas partes do corpo se ajusta à pretendida.

• Análise da rede neuronal para 0,8α = :

Para 0,8α = pretende simular-se a situação em que a eficácia sináptica de

algumas das ligações entre as diferentes unidades da rede neuronal se encontra algo enfraquecida.

De seguida apresenta-se a implementação desta rede neuronal no Simulink e a saída da unidade 1:

Figura 16 Espectro de frequências da saída da unidade 1 para α=1.

0 0.5 1 1.5 2 2.5 30

10

20

30

40

50

60


Frequência (Hz)

Am

plitu

de



Ao analisar-se a representação da actividade (saída) da unidade 1, verifica-se

que esta apresenta inicialmente oscilações de grande amplitude que tende a diminuir rapidamente com o tempo. No entanto, relativamente ao caso anterior ( 1α = ), esta diminuição é menos acentuada. Também para esta situação, a saída da unidade 1 atinge, em determinado momento e devido à natureza da função sigmóide, um valor


f(u)

g(y6)

f(u)

g(y5)

f(u)

g(y4)

f(u)

g(y3)

f(u)

g(y2)

f(u)

g(y1)

y1

To Workspace

In1

In2Out1

Subsistema 6

In1

In2Out1

Subsistema 5

In1

In2

Out1

Subsistema 4

In1

In2Out1

Subsistema 3

In1

In2

Out1

Subsistema 2

In1

In2

Out1

Subsistema 1

Scope

0.8

Alfa

0.8Alfa


0 5 10 15 20 25 30 35 40 45 50-0.1

-0.05

0

0.05

0.1

0.15

0.2

0.25

0.3


Tempo (s)

Am

plitu

de



praticamente constante (0,0607). Este valor é, portanto, superior ao valor obtido na situação anterior ( 1α = ).

Procedendo-se à análise das frequências da saída da unidade 1, obtém-se a seguinte representação:

O espectro obtido apresenta dois picos para as frequências de 0Hz e 0,2599Hz.

O pico de 0Hz corresponde ao período de tempo em que a saída da unidade 1 é constante. O outro pico (de pequena amplitude mas, ainda assim, superior à registada no caso de 1α = , para o pico equivalente) é provocado pelo curto período inicial de oscilações.

Do ponto de vista fisiológico, esta situação corresponde a um estado intermédio entre a situação fisiológica normal e a situação resultante da doença de Parkinson. Verifica-se, assim, a necessidade de um maior período de tempo para que a actividade da unidade 1 atinja um valor constante, ao qual corresponde o atingir de uma posição de equilíbrio da parte do corpo considerada.

• Análise da rede neuronal para 0,5α = :

Para 0,5α = pretende simular-se a situação em que a eficácia sináptica de algumas ligações está bastante diminuída.

Após a implementação em Simulink obtém-se, no Scope, a representação da saída da unidade 1:


0 0.5 1 1.5 2 2.5 30

5

10

15

20

25


Frequência (Hz)

Am

plitu

de



Através da observação da saída da unidade 1, verifica-se que esta apresenta

algumas oscilações iniciais de grande amplitude e que esta diminui com o decorrer do tempo. Após estas oscilações de maior amplitude verificam-se, durante um grande período de tempo, oscilações regulares que diminuem lentamente até se atingir um instante (após um enorme período de tempo) em que se tornam nulas (a saída da unidade 1 torna-se constante.


f(u)

g(y6)

f(u)

g(y5)

f(u)

g(y4)

f(u)

g(y3)

f(u)

g(y2)

f(u)

g(y1)

y1

To Workspace

In1

In2

Out1

Subsistema 6

In1

In2

Out1

Subsistema 5

In1

In2

Out1

Subsistema 4

In1

In2

Out1

Subsistema 3

In1

In2

Out1

Subsistema 2

In1

In2Out1

Subsistema 1

Scope

0.5

Alfa

0.5Alfa


0 5 10 15 20 25 30 35 40 45 50-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

0.25


Tempo (s)

Am

plitu

de



Procedendo-se à análise das frequências da saída da unidade 1, obtém-se a seguinte representação:

Neste caso obtém-se um pico bem definido para uma frequência de 0,2599Hz, o

que reflecte o comportamento regular da saída da unidade 1, como foi referido anteriormente. Verifica-se, também, um pico de pequena amplitude para uma frequência de 0Hz, indicativo do curto período de tempo em que a saída apresenta um valor constante.

Do ponto de vista fisiológico, o comportamento oscilatório da saída da unidade 1 (com tendência para atingir um valor constante), não se verifica, pois é necessário um grande período de tempo até que se atinja o valor constante. Como a posição corporal é alterada constantemente, não se dispõe do tempo necessário para ajustar de forma fina (precisa) a posição da parte do corpo considerada. Como consequência, regista-se uma actividade regular e oscilatória das diversas unidades, o que se reflecte no surgimento de um tremor permanente e também regular, o temor Parkinsoniano.

Análise da dinâmica da rede neuronal

De forma a descrever o comportamento da rede neuronal considerada, e consequentemente o tremor a ela associada, utilizam-se seis equações diferenciais.

Deste modo não é possível fazer a análise da dinâmica deste sistema recorrendo ao programa PPlane6, como vinha sendo feito nas aulas práticas.

Assim, de forma a analisar-se a evolução do comportamento da rede neuronal para diverso valores de α e diferentes funções resposta, escolheram-se aleatoriamente conjuntos de saídas de 3 unidades diferentes. A representação destas foi feita recorrendo ao Matlab através da função plot3(…), após ter sido executada a simulação em Simulink, o qual guardou, recorrendo a um bloco adequado, os resultados numa variável.


0 0.5 1 1.5 2 2.5 30

0.5

1

1.5

2

2.5

3

3.5

4

4.5Espectro de Frequências de y1

Frequência (Hz)

Am

plitu

de



Unidades da rede neuronal com uma resposta de acordo com o degrau de Heaviside

• α=1 Considerando y1, y2, y3 :

-0.3-0.2

-0.10

0.10.2

0.30.4

0.5

-0.3-0.2

-0.10

0.10.2

0.30.4

0.50.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

y1

Dinâmica da rede neuronal para alfa=1

y2

y3

Figura 23: Representação da trajectória resultante da evolução da rede neuronal, escolhendo as saídas y1, y2, y3.

-0.3-0.2

-0.10

0.10.2

0.30.4

0.50.6

-0.4-0.3

-0.2-0.1

00.1

0.20.3

0.40.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

y2


y5

y6

Figura 24 Representação da trajectória resultante da evolução da rede neuronal, escolhendo as saídas y2, y5, y6.

Para qualquer um dos conjuntos de saídas das unidades a que se recorreu, não é possível identificar qualquer tipo de singularidade. Este resultado reflecte o comportamento irregular e aleatório da saída das unidades para estes parâmetros, como referido atrás.



• α=0.8

Considerando y1, y2, y3

-0.4-0.3

-0.2-0.1

00.1

0.20.3

0.4

-0.4-0.3

-0.2-0.1

00.1

0.20.3

0.40.5

0.6-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

y1

Dinâmica da rede neuronal para alfa=0.8

y2

y3



-0.4-0.3

-0.2-0.1

00.1

0.20.3

0.40.5

0.6

-0.4-0.2

00.2

0.40.6

0.81

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

y2


y5

y6

Figura 26 Representação da trajectória resultante da evolução da rede neuronal, escolhendo as saídas y2, y5, y6. Também neste caso para qualquer um dos conjuntos de saídas escolhidas, não é possível identificar qualquer tipo de singularidade, resultado do comportamento irregular, aleatório anteriormente verificado para a saída das unidades, consequentemente do tremor, para estes parâmetros. Contudo parece evidenciar-se uma tendência para o sistema evoluir atingindo um ciclo limite em torno de um centro.



• α=0.5 Considerando y1, y2, y3

-0.4-0.3

-0.2-0.1

00.1

0.20.3

0.4

-0.4-0.3

-0.2-0.1

00.1

0.20.3

0.40.5

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

y1


y2

y3



-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5

0

0.2

0.4

0.6

0.8

1-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

y2


y3

y4


Para este caso é possível verificar que o sistema, partindo das suas condições iniciais, evolui de tal forma que descreve uma trajectória aproximadamente rectangular (ciclo limite) em torno de um centro. Tal deve-se ao comportamento periódico das saídas das diversas unidades para estes parâmetros.



Unidades da rede neuronal com uma resposta de acordo com uma função sigmóide ( β=4 )

• α=1


-0.050

0.050.1

0.150.2

0.250.3

0.35

-0.2-0.1

00.1

0.20.3

0.40.5

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

y1


y2

y3


Considerando y2, y4, y5:

-0.2-0.1

00.1

0.20.3

0.40.5

-0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

0.25

-0.5

0

0.5

y2


y4

y5

Figura 30 Representação da trajectória resultante da evolução da rede neuronal, escolhendo as saídas y2, y4, y5. Para estes casos é possível identificar uma singularidade. A trajectória resultante da evolução da rede neuronal tende para o foco oscilando, por isso designa-se este tipo de singularidade de foco estável. No caso em que se analisou y1, y2 e y3 o foco é



(y1=0.10456,y2=0.10456, y3=0.10456), quando se analisou y2, y4 e y5 o foco é (y2=0.10456, y4=0.10456,y5=0.10456).

• α=0.8

Considerando y1, y2 e y3 tem-se:

-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 -0.2

0

0.2

0.4

0.6

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

y2

y1

Dinâmica da rede neuronal para Alfa=0.8

y3


Considerando y2, y4 e y6, tem-se:

-0.2-0.15

-0.1-0.05

00.05

0.10.15

0.20.25

0.3

-0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

-0.4

-0.2

0

0.2

0.4

y2


y4

y6


Novamente para estes casos é possível identificar uma singularidade. A trajectória resultante da evolução do comportamento da rede neuronal, evolui em cada



um dos casos para um foco oscilando. Designa-se por isso este tipo de singularidade por foco estável. Para a análise de y1, y2, y3 o foco atingido é (y1=0.069653,y2=0.069653,y3=0.17506), para a análise de y2,y4,y6 o foco é (y2=0.069653,y4=0.069653,y6=0.16906).

• α=0.5

Considerando y1, y3, y5 tem-se:

-0.2-0.15

-0.1-0.05

00.05

0.10.15

0.20.25

0.3

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0

0.2

0.4

0.6

0.8

y3

Dinâmia da rede neuronal para Alfa=0.5

y1

y5


Considerando y2, y4 y6:

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

0.25 -0.25-0.2

-0.15-0.1

-0.050

0.050.1

0.150.2

-0.5

0

0.5

1

y4


y2

y6




À semelhança dos casos anteriores, também aqui é possível identificar uma singularidade. Analisando as representações gráficas obtidas é possível visualizar que a trajectória descrita por cada uma delas tende para um foco oscilando (foco estável). Quando se considera y1, y3, y5 o foco atingido é (y1=0.0051416,y3=0.50029,y5=0.47944). Para y2, y4, y6 o foco para o qual a trajectória tende é (y2=0.0051416,y4=0.0051416,y6=0.51954).

De notar que para chegar a estes resultados para α=0.8 e α=0.5, foi necessário efectuar a simulação para um maior período de tempo (10000) visto o sistema evolui devagar, e soa assim foi possível visualizar uma trajectória que nos permitisse chegar a estes resultados.



Conclusão

Após a realização deste trabalho, e depois da análise dos resultados, conclui-se que este modelo representa, de forma aceitável, o comportamento da rede neuronal envolvida no controlo do tremor humano, bem como as alterações do mesmo à medida que a doença de Parkinson evolui. No entanto os resultados obtidos por este modelo apresentam algumas imprecisões relativamente aos esperados.

Estas imprecisões dizem respeito sobretudo à periodicidade dos sinais de saída obtidos nas unidades (consequentemente no tremor), reflectindo-se nos espectros de frequência obtidos. A forma destes foi a esperada, contudo os picos verificados ocorrem para frequências bastante inferiores às esperadas (como referido no inicio do trabalho). Estas imprecisões poder-se-ão dever ao número reduzido de unidades considerado e/ou a limitações associadas ao Simulink.



Anexo

• Descrição do script usado na determinação dos espectros de frequência

Na determinação do espectro de frequência da saída da unidade 1, nas diversas situações ao longo do trabalho, recorreu-se ao seguinte script: Y=fft(y1); Pyy = Y.* conj(Y)/5001; f = 100*(0:150)/5001; plot(f,Pyy(1:151)) title('Espectro de Frequências de y1') xlabel('Frequência (Hz)') ylabel('Amplitude') Explicação do script: Y=fft(y1); Determina, recorrendo à função do Matlab fft (fast Fourier transform), a transformada discreta de Fourier dos pontos do vector y1, o qual contém o resultado da simulação efectuada no Simulink. Pyy = Y.* conj(Y)/5001;

Determina a magnitude de cada um dos elementos do vector Y, normalizando à dimensão do vector Y.

f = 100*(0:150)/5001; Cria um vector linha que irá servir de suporte para definir a escala de

frequências do espectro de frequência determinado posteriormente. Escolhe-se um vector linha com 151 elementos pois este intervalo de frequência é suficiente para a análise dos sinais em questão. Divide-se por 5001 para normalizar à dimensão do vector Y de forma a ficar coerente com Pyy. Multiplica-se por 100 (1/0.01), pois na simulação em Simulink fazem-se contagens em intervalo de tempo de 0.01 unidades, para ter a frequência em Hz. plot(f,Pyy(1:151))

Faz a representação do espectro de frequência de y1, ou seja, representa a magnitude da transformada de Fourier de y1 em função da frequência em Hz.

title('Espectro de Frequências de y1') xlabel('Frequência (Hz)') ylabel('Amplitude')

Introduzem formatações ao nível do título, eixo das abcissas, quer ao nível do eixo das ordenadas, no gráfico obtido no comando anterior.



Teste do script: De forma a testar o script utilizado, implementou-se o seguinte modelo em Simulink:

y1

To Workspace

Sine Wave1Sine Wave

frequência=2*pi*2 rad /s= 2Hz

Sine Wavefrequência=2*pi*1 rad /s

= 1Hz

Add

Figura 35 Diagrama de blocos representativo da simulação feita para testar o script que foi utilizado para determinar o espectro de frequências dos sinais ao longo do trabalho. O bloco Workspace cria um vector y1, no Matlab, contendo as coordenadas da simulação.

Somam-se duas funções seno, uma com frequência de 2π rad/s = 1Hz e outra com frequência 4π rad/s = 2Hz. Como esperado, obtêm-se no espectro de frequências 2 picos, um na frequência 1Hz e outro na frequência 2Hz.

0 0.5 1 1.5 2 2.5 30

200

400

600

800

1000

1200


Frequência (Hz)

Am

plitu

de

Figura 36 Espectro de frequências do sinal utilizado como teste.



Bibliografia

Berne, Robert M. et al; Physiology; 5ª edição; Mosby, 2004 Beuter, Anne et al; Nonlinear Dynamics in Physiology and Medicine; Springer, 2003

Haykin, S., Veen, B. Van; Signals and Systems; 2ª edição; John Wiley & Sons, Inc.; 2003

Mackay, William A.; Neurofisiologia Sem Lágrimas; Edições Gulbenkian; 1999 Seeley, Rod R. et al; Anatomia & Fisiologia; 1ª edição, Lusodidacta, 1997 McAuley, J. H.; Marsden, C. D.; Physiological and pathological tremors and rhythmic central motor control; Brain, Vol. 123, No. 8, 1545-1567, 2000

Documents

M Moodelização Matemática do Tremor Humano na Doença de … · As oscilações associadas a estes mecanismos de feedback, juntamente com as oscilações resultantes da própria