SOLUÇÃO SÉRIE PARA O CAMPO ELÉTRICO INDUZIDO POR ESTIMULADORES
MAGNÉTICOS NA GEOMETRIA SLINKY GENERALIZADA
MARCÍLIO FEITOSA, EDUARDO FONTANA
Grupo de Fotônica, Dep. de Eletrônica e Sistemas, Centro de Tecnologia e Geociências, Universidade
Federal de Pernambuco. Av. Acadêmico Hélio Ramos, Cidade Universitária, Recife, PE, Brasil.
E-mails: [email protected],[email protected]
Abstract Magnetic stimulation is a technique to excite biological tissues by means of a time-varying magnetic field. This induced electric field can depolarize the cell membrane so as to evoke an action potential that propagates along neurons, eventually
being transmitted to other neurons or to a muscular cell. Design of a magnetic stimulator requires modeling of the impulse propagation along the nerve cell, as well as numerical simulations for coil design optimization to determine adequate excitation levels
as well as the degree of focalization on a given target cell. In this paper we report on a new methodology to calculate the stimulation field for the case of the traditional slinky coil geometry, that greatly reduces computation time, thus facilitating simulation
studies of the dynamics of electric impulse propagation along a nerve cell.
Keywords Biomagnetics, Biomedical applications of electromagnetic radiation, Biomedical engineering, Neuromuscular
stimulation.
Resumo A estimulação magnética consiste na excitação de tecidos orgânicos por um campo magnético variável no tempo. No
desenvolvimento dessa técnica, interesse crescente tem sido devotado na literatura ao emprego de estimuladores magnéticos na
configuração slinky (mola flexível de brinquedo), devido à possibilidade de melhor focalização magnética, relativamente a outras
configurações. Neste trabalho é desenvolvida uma solução série para o cálculo do campo elétrico induzido em uma configuração
Slinky generalizada válida para um número arbitrário de espiras do estimulador, computacionalmente eficiente para o projeto e
estudo da dinâmica do sistema.
Palavras-chave Efeitos biológicos da radiação eletromagnética, Bioeletromagnetismo, Cálculo de campos eletromagnéticos.
1
Introdução
A estimulação magnética de tecidos nervosos é uma
técnica que vem ganhando importância elevada para
médicos e demais profissionais de saúde, principalmente para aqueles que atuam na área fisioterápica
ou de reabilitação neuromuscular. É uma ferramenta
que, da mesma forma que a estimulação elétrica,
auxilia no diagnóstico e no tratamento de algumas
doenças que afetam o sistema nervoso [1] mas, devido a características particulares, tem se sobressaído
sobre esta técnica tradicional. Seu princípio consiste
em estimular tecidos através de um campo elétrico
induzido por um campo magnético externo, variável
no tempo. O campo elétrico induzido pode despolarizar a membrana celular e iniciar a propagação de
um potencial de ação, que pode atingir outra célula
nervosa ou uma célula muscular, causando uma contração da mesma. A característica que torna a estimulação magnética mais atrativa é o fato de o campo
magnético penetrar regiões eletricamente isoladas,
como as camadas de gordura e ossos, praticamente
sem sofrer atenuação [1]-[2]. Devido à alta impedância elétrica dessas regiões, a densidade de corrente
exigida para se conseguir uma excitação efetiva, através da estimulação com contato elétrico direto, é
elevada e capaz de produzir irritações na pele e
queimaduras.
O projeto de um estimulador magnético envolve
a modelagem da dinâmica de propagação do impulso
elétrico ao longo da fibra nervosa, necessária para se
estabelecer o nível mínimo de estímulo, mas também
envolve o estudo de como a geometria das bobinas
afeta a focalização do campo. É nesse contexto que o
presente trabalho se torna relevante, pois, através de
uma nova formulação matemática, diminui o esforço
computacional necessário à simulação do campo
gerado para cada geometria de bobina proposta.
2
Modelo Dinâmico da Membrana de uma Célula Nervosa
A dinâmica da propagação de um distúrbio elétrico
ao longo de um neurônio ou de uma fibra nervosa
pode ser descrita com o auxílio do modelo desenvolvido por Hodgkin-Huxley [3], indicado na Fig.1, no
qual o axônio, isto é, a fibra nervosa ao longo da
x
EL
EK ENa
gL
gK gNa
cm
ri
I(x+∆x)
∆x
I(x)
V(x)
V(x+∆x)
Fig1. Modelo de parâmetros distribuídos para a membrana de um
axônio.
qual o estímulo se propaga, é modelado como uma
linha de transmissão de parâmetros distribuídos. Nesse modelo os canais iônicos para o sódio, o potássio
e os demais íons envolvidos no processo são representados por condutâncias (gNa, gK e gL respectivamente) que são funções não-lineares da tensão entre
as faces da membrana [3]. As paredes da membrana
são modeladas eletricamente como as placas de um
capacitor cuja capacitância por unidade de comprimento é dada por cm e a solução iônica presente no
interior do axônio é representada por um resistor cuja
resistência por unidade de comprimento é ri. As diferentes concentrações iônicas nas soluções interna e
externa da célula produzem tensões de equilíbrio
para cada íon. Estas são representadas por fontes de
tensão contínuas ENa, EK, e EL, associadas aos gradientes iônicos produzidos pelos íons Na, K e demais
íons envolvidos.
Para o caso de uma fibra nervosa localizada ao
longo do eixo x, sob a ação de um campo elétrico Ex
induzido magneticamente, o potencial elétrico através da membrana V satisfaz a equação diferencial
Λ2
∂ 2V
∂x
2
−V = τ
∂E x
∂V
+ Λ2
,
∂t
∂x
(1)
onde as constantes de espaço e tempo são dadas respectivamente por
r
Λ =  m
 ri
1
2
 ,

τ = rm c m ,
(2)
(3)
com rm representando a resistência transmembrana equivalente para o modelo. Essa resistência é
uma função não linear do potencial de membrana.
Podemos observar que (1) tem, a princípio, a estrutura de uma equação de difusão devido a sua segunda
derivada no espaço e primeira derivada no tempo. De
fato, se a intensidade do estímulo (representado pelo
gradiente longitudinal da componente na direção x
do campo elétrico induzido) for inferior a determinado valor, o distúrbio de tensão produzido em um
ponto da membrana vai se difundir longitudinalmente, mas sem se propagar. Esse distúrbio causará apenas um pequeno desvio na tensão de equilíbrio da
membrana a qual se re-estabelece logo em seguida.
Existe, no entanto, um valor limiar para a tensão de
membrana que, se for superado devido à intensidade
do estímulo, um potencial de ação irá se propagar
como uma onda não amortecida cuja velocidade depende de parâmetros intrínsecos ao nervo [1].
3
Estudo da Geometria Slinky Generalizada
Como pode ser observado em (1), o termo responsável pela estimulação é a derivada em x da componente x do campo elétrico externo. Esse termo, que representa a função de ativação, depende da configuração geométrica da bobina de estimulação, além da
distribuição temporal do pulso de corrente produzido
pela descarga de um banco de capacitores e controlado por um dispositivo semicondutor de alta potência. Faz-se necessário então um estudo detalhado
sobre qual configuração possibilita alcançar estímulos mais efetivos e mais focalizados, utilizando-se
um mínimo de energia. Em algumas pesquisas o uso
da bobina circular clássica [4]-[5] é bastante difundido (Fig.2), mas outros resultados da literatura a respeito do tema indicam que as bobinas do tipo slinky
apresentam melhor rendimento [2]. Nessa geometria,
as espiras são todas de mesmo raio e têm um ponto
em comum. A geometria estudada no presente trabalho, mostrada na Fig.3 é uma generalização da configuração slinky, em que os raios, correntes e inclinações relativas das espiras são arbitrários. O nome
slinky foi adotado por Richard James em 1943 para
designar um conhecido brinquedo por ele criado, que
tem a forma de uma mola flexível. A configuração
Slinky traz como benefício a diminuição do espalhamento lateral do campo usado para estimulação,
além de minimizar a possibilidade de estímulo de
outros nervos que se encontrem próximo ao nervo
alvo [2].
Em estudos anteriores [6] mostramos que, para a
geometria slinky de duas espiras, se mantivermos o
ápice da bobina a uma distância fixa do nervo, obtemos uma distribuição do gradiente longitudinal do
campo elétrico com uma largura rms mínima para um
ângulo de aproximadamente 47º entre o plano das
espiras e a superfície de trabalho (Fig.4). Nesses
mesmos estudos, observamos que essa bobina, para
esse mesmo ângulo, concentra o gradiente do campo
em uma região aproximadamente 35% menor que a
região excitada por uma bobina circular convencional [6]-[7].
Fig.2. Aplicação do pulso magnético no nervo mediano, utilizandose uma bobina circular clássica.
z
i3
(a)
zk
X X’
X
z
r3
θk
X’
ik
rk
ik
rk
yk
Rk
Xk
i2
r2
r1
αk
y
i1
y
x
x
x
Fig.3. Geometria da configuração Slinky generalizada.
4
âzk
Modelo Teórico
Em busca de uma configuração de bobinas que permita a geração de um campo de estimulação mais
eficaz, foram realizados estudos sobre como a geometria das bobinas afeta esse campo. As variáveis
envolvidas nesses estudos foram o número de espiras
na bobina, o ângulo que cada espira faz com a superfície de trabalho (plano paralelo àquele onde está
posicionado o nervo alvo do estímulo), o raio de cada espira e a intensidade e a forma de onda da corrente que percorre cada espira. Esses parâmetros
estão detalhados na Tabela I. Cada nova espira adicionada à bobina é definida de forma uniformizada
de acordo com os parâmetros citados com relação a
um sistema de coordenadas rotacionado em torno do
eixo x como indicado na Fig.5a. Essa uniformização
permitiu uma redução no trabalho de se descrever
cada bobina individualmente em função do plano
cartesiano original.
Foi considerada a determinação do campo elétrico induzido devido a uma configuração de N espiras
circulares, onde a k-ésima espira possui um raio rk e
é percorrida por uma corrente ik. Todas as espiras
Relative rms width
1.80
1.70
Variable apex
1.60
1.50
1.40
(b)
âz
âyk
αk
âφk
αk
ârk
ây
φk
âx
Fig.5. (a) Geometria e parâmetros de definição de uma dada espira
na configuração Slinky generalizada. (b) Vetores de definição de
uma espira na configuração Slinky generalizada.
têm um ponto em comum na origem do sistema do
laboratório xyz. O sistema xykzk representa o sistema
de coordenadas da k-ésima espira, desviada de um
ângulo αk do plano xy.
A formulação convencional para o cálculo do
campo elétrico em um determinado ponto do espaço
devido a uma dada distribuição de corrente, que emprega integrais elípticas, exige um grande esforço
computacional o que introduz um longo atraso nas
simulações numéricas da dinâmica do potencial de
ação. Devido a isso, foi desenvolvida uma expressão
analítica para o campo elétrico induzido bem como
para sua derivada espacial usada como fonte de estímulo em (1).
Em uma aproximação quase estática, o campo
elétrico induzido por uma espira percorrida por uma
corrente é dado por [8]
G
G
∂A
(4)
E=−
∂t
G
onde A representa o vetor potencial magnético, que
pode ser escrito da forma
Stationary apex
1.30
0
10
20
30
40
50
60
70
80
G µ
A= 0
4
90
Angle(degrees)
Fig.4. Largura rms relativa do gradiente do campo elétrico induzido ao longo do nervo.
com
N
∑ i F (R , u
k
k =1
k
k
k
)âφk ,
(5)
TABELA I. PARÂMETROS PARA A CONFIGURAÇÃO SLINKY
GENERALIZADA.
Parâmetro
G
Xk
G
X
αk
rk
ik
xyz
xykzk
Rk, θk, φk
â x , â y , âz
Definição
Localização do centro da k-ésima espira.
Ex =
Vetor posição do ponto de observação.
Ângulo entre o plano da espira e o plano xy.
Raio da k-ésima espira.
Corrente que percorre a k-ésima espira.
Sistema de coordenadas do laboratório.
Sistema de coordenadas da k-ésima espira
Coordenadas esféricas no sistema xykzk
Vetores unitários do sistema xyz.
onde
Vetores unitários do sistema xykzk
â x , â yk , â zk
Vetores unitários (radial e azimutal), no
plano xyk.
ârk , âφk
µ0
4π
di k rk ( y cos α k + zsinα k − rk )
Fk
1/ 2
dt
R 1− u 2
N
∑
k =1
(
k
k
(17)
)
− ysinα k + z cos α k
.
(18)
Rk
A derivada de (17) pode ser posta na forma
uk =
di k rk ( y cos α k + zsinα k − rk )H k
N
µ
∂E x
= −x 0
4
∂x
∑ dt
(
R k2 1 − u k2
k =1
)
, (19)
com
2π
Fk =
∫
0
em que
rk cos φdφ
G G ,
X − X'
(6)
G G
X − X ' = R k â Rk − rk â r 'k .
Assim,
G
µ
E=− 0
4π
N
∑
k =1
di k
Fk (R k , u k )âφk .
dt
(8)
â rk = cosφ k â x + sin φ k â yk ,
(10)
â yk = cos α k â y + sinα k â z .
(11)
O vetor posição do ponto de observação no sistema xykzk pode ser obtido em função do sistema de
coordenadas xyz por
G
G G
Rk = X − X k .
Utilizando a relação
(12)
(XG − XG )• (XG − XG ) ,
k
k
F2k =
(
G
X = xâ x + yâ y + zâ z ,
)
(14)
(15)
obtém-se
2
2
Rk = x + y + z
2
∫ R
0

+ rk2
(
(
(cosφ )2
2
k

 , (20)
1/ 2 

)
+ rk2
)
cos φ 

(
1/ 2
− 2rk Rk 1 − u k2
3/ 2
3/ 2
u k = cos θ k .
dφ , (21)
dφ (22)
(23)
As equações (19)-(23) representam uma generalização da formulação integral para o cálculo do
campo de estimulação gerado por uma bobina do
tipo slinky. Essa formulação demanda um elevado
esforço computacional no cálculo da distribuição do
gradiente pois, para cada ponto no espaço, é necessário o cálculo de três integrais numéricas.
A formulação integral pode ser evitada pelo uso
de uma expansão em harmônicos esféricos [8] para o
inverso do denominador de (7).
De acordo com o teorema da adição para harmônicos esféricos [8],
∞
l
∑∑ 2l1+ 1 rr
l
<
l +1
>
l =0 m = −l
* π
Υlm
( , φ k ' )Υlm (θ k , φ k ) , (24)
2
com Ylm representando o harmônico esférico de índices (l,m). Inserindo (24) em (6) e aplicando-se as
propriedades de ortogonalidade dos harmônicos esféricos [9], obtém-se
Fk (Rk , u k ) = rk
− 2rk ( y cosα k + zsinα k ) . (16)
∞
∑
l =0
Substituindo-se (9) em (8), com o emprego de
(11) e considerando-se apenas a componente de interesse no campo elétrico induzido, teremos
)
cosφ 

1/ 2
− 2rk Rk 1 − u k2
+ rk2
e
G G = 4π
X − X'
G
X k = rk cos α k â y + sinα k â z
1
F1k − rk F2k  + Fk

Rk 1 − u k2
2π
1
e com

2
k
(13)
e com base na Fig.5,
1/ 2
cosφ
∫ R
0
(9)
)
2π
F1k =
âφk = −sinφ k â x + cos φ k â yk ,
Rk =
(
(7)
Os vetores unitários mostrados na Fig.5b obedecem às relações
com


H k =  Rk 1 − u k2


onde
al ≡
al
R<2l +1
R>2l + 2
P21l +1 (u k )
(−1)l (2l − 1)!!
,
2l (l + 1)!
(25)
(26)
com P12l+1 representando o polinômio associado de
Legendre de índices (2l+1,1).
Substituindo-se (25) em (5), obtemos a expressão para o campo elétrico induzido de (4). O gradiente longitudinal desse campo é obtido após algumas
manipulações algébricas, podendo ser expresso no
sistema xyz na forma
dE x
µ
=− 0
dx
4
di k x( y cos α k + zsinα k − rk )
Gk ,
R k4
N
∑ dt
k =1
(27)
∞
∑

∑ a  Rr
l =1
∞
r
al  k
 Rk
l
k
k
l =1




2l + 2




2l +1
H 1k , rk < Rk
,
(28)
H 2 k , rk > Rk
H 1k = − (2l + 3)R k P12 l +1 + ( ysin α k − z cos α k )P 2 2 l +1 , (29)
e
H 2 k = 2lR k P12l +1 + ( ysin α k − z cos α k )P 2 2l +1 . (30)
Para se evitar problemas computacionais no cálculo das derivadas dos polinômios de Legendre, particularmente para valores específicos de x e elevados
valores do índice l, nas equações (29) e (30) as derivadas primeira e segunda dos polinômios foram expressas como combinações dos polinômios de Legendre [9] da seguinte forma:
P1l (u k ) =
l
u k2
−1
e
P2l (u k ) =
(
u k2
l
)
−1
2
[u k Pl (u k ) − Pl −1 (u k )]
{[(l − 1)u
5
2
k
]
uma matriz, como indicado a se" I max k
" Sinal k
t0 k
"
σk
"
" Raio k
" Ang k
" I max N 

" Sinal N 
t0 N 
"

σN 
"

" Raio N 
" Ang N 
(33)
Onde, para cada uma das N espiras, temos:
onde temos



Gk = 



como colunas em
guir.
 I max1

 Sinal1
 t
0
LOOP =  1
 σ1

 Raio1
 Ang
1

(31)
}
− (l + 1) Pl (u k ) + 2 xPl −1 (u k ) . (32)
Discussão
Foram realizadas simulações computacionais no ambiente Mathcad 13 para calcular a distribuição espacial do campo elétrico induzido bem como do seu
gradiente, tanto para a formulação integral como para
a formulação série baseada nas expansões de Legendre. Nesses cálculos foram utilizadas as funções disponíveis no ambiente para o cálculo das integrais e
das funções de Legendre.
Apesar de o algoritmo desenvolvido aceitar
qualquer distribuição temporal para os pulsos de
corrente aplicado às espiras, optamos por pulsos
gaussianos para as simulações. Os parâmetros de
definição de cada espira são informados ao algoritmo
Imax
- Amplitude máxima do pulso gaussiano de
corrente;
Sinal
- Sinal dependente do sentido da corrente;
- Centróide do pulso gaussiano de corrente;
t0
σ
- Desvio padrão do pulso de corrente;
Raio
- Raio da espira;
Ang
- Ângulo de inclinação da espira com relação ao plano xy.
O algoritmo inicialmente define as funções que
serão utilizadas: os polinômios de Legendre e suas
derivadas primeira e segunda, os polinômios associados de Legendre e a derivada temporal do pulso
gaussiano de corrente. Em seguida são definidas as
coordenadas esféricas do sistema de coordenadas
rotacionado em função do sistema cartesiano fundamental. Só então o campo elétrico é calculado bem
como a sua derivada longitudinal. Esta função é utilizada como entrada em (1) e então o potencial de
membrana é calculado.
Para o caso da expansão polinomial, as expressões em série foram truncadas quando a precisão do
resultado chegou a 0.1%. Essa foi a mesma tolerância adotada no cálculo das integrais numéricas realizado pelo Mathcad. Observou-se que a solução série
permitiu uma redução no tempo de computação de
um fator de pelo menos 10 relativamente à formulação integral. Esta vantagem se torna ainda maior
quando o número de espiras na bobina aumenta. De
fato, para um elevado número de espiras de corrente
(>10), a formulação integral torna-se praticamente
inviável.
Foi observada uma desvantagem no uso da formulação por expansão de Legendre, na qual cada
termo da expansão do campo elétrico, apesar de ser
uma função contínua, possui derivada descontínua
para a k-ésima espira quando a distância do centro da
espira até um ponto no espaço é igual ao raio da espira. Obviamente, se todos os termos da série fossem
incluídos, este problema não ocorreria. No entanto,
quando a série é truncada, a derivada do campo se
torna descontínua nessa condição. Para o caso onde
todas as espiras se encontram no hemisfério superior
e o nervo alvo se encontra no hemisfério inferior,
paralelo ao eixo x, no plano y = 0, esse problema não
ocorre. O problema de descontinuidade pode ser resolvido expandindo-se a derivada de (6) com respeito à variável x. Esse procedimento, no entanto, re-
quer o emprego dos polinômios ultra-esféricos de
Gegenbauer [9], e ainda está em desenvolvimento.
6
Conclusões
Foi desenvolvido um método computacional eficiente para o cálculo do gradiente espacial do campo
elétrico induzido por estimuladores na configuração
Slinky generalizada. O método consiste no cálculo
do campo, e do seu gradiente, através de uma expansão em série baseada nos polinômios de Legendre, ao
invés do cálculo por integrais elípticas. Apesar de
alguns problemas terem sido encontrados, como o
aparecimento de algumas descontinuidades no cálculo do gradiente em situações específicas, o procedimento se mostrou bem mais eficiente, reduzindo o
tempo de simulação em até 10 vezes. Este método
está sendo utilizado em estudos de otimização e projeto de estimuladores magnéticos com boa focalização de campo para aplicações em nervos periféricos.
O comportamento em tempo real da dinâmica das
fibras nervosas sob influência do campo induzido
está sendo estudado e novos resultados serão publicados em breve.
Referências Bibliográficas
[1] J. Malmivuo, R. Plonsey, “Bioelectromagnetism:
principles and applications of bioelectric and
biomagnetic fields”, Oxford University, Oxford,
1995.
[2] C. Ren, P. P. Tarjan, and D. B. Popovic, “ A
Novel Electric Design for Electromagnetic
Stimulation – The Slinky Coil,” IEEE
Transactions on Biomedical Engineering, vol.
42, no. 9, pp. 918-925, Sep. 1995.
[3] A. L. Hodgkin and A. F. Huxley, "A quantitative
description of membrane current and its
application to conduction and excitation in
nerve", Journal of Physiology, no. 117, pp.500544, 1952.
[4] M. George, E. M. Wassermann and R. M. Prost,
"Transcranial
magnetic
stimulation:
a
neuropsychiatric tool for the 21st century",
Journal of Neuropsychiatric and Clinical
Neurosciences, vol. 8, pp.373-382, 1996.
[5] J.
Ruohonen,
“Transcranial
magnetic
stimulation: modeling and new techniques”,
Helsinki University of Technology, 1998.
[6] M. A. F. Feitosa and Eduardo Fontana,
“Prospects for the Development of a Magnetic
Stimulation Device for Human Tissue,”
Proceedings of the 2005 International
Microwave and Optoelectronics Conference.
Brasília 2005. v. 1. p. 1-4.
[7] Kai-Hsiung Hsu, Srikantan S. Nagarajan, and
Dominique M. Durand, “Analysis of Efficiency
of Magnetic Stimulation,” IEEE Transactions on
Biomedical Engineering, vol.50, no. 11, pp.
1276-1285, Nov. 2003.
[8] J. D. Jackson, “Classical Electrodynamics”, 3rd.
edition. New York, John Wiley & Sons Inc.,
1999, pp. 174-242.
[9] Arfken and Weber, “Mathematical Methods for
Physicists”, 5th. edition. United States,
Academic Press, 2001, pp. 745-820.