4
equazioni iconali. Una classificazione più utile si basa invece sulla scala
del modello [WEI97], distinguendo in:
• Modelli per cellule singole (modelli di membrana). Si occupano della
descrizione della dinamica elettrochimica cellulare e si fondano sulla
disponibilità di dati sperimentali, perciò la loro evoluzione procede con
l’aggiunta di nuove osservazioni. Questi modelli si basano sugli studi di
Hodgkin e Huxley [HOD52] sulle correnti ioniche di membrana nelle
cellule nervose. In tali modelli la membrana cellulare viene
rappresentata da una capacità lineare attraversata da una corrente pari
alla somma di una corrente dovuta a ioni e di un’eventuale corrente
imposta dall’esterno.
Oggi esiste un gran numero di tali modelli. A puro titolo di esempio
ricordo il modello di Earm-Noble [EAR90] di cellule atriali e quelli di
Beeler-Reuter [BEE77] e Luo-Rudy [LUO] per cellule ventricolari. Un
discorso a parte merita il modello di FitzHugh-Nagumo con le sue
varianti, che benché descriva l’attività di una singola cellula viene usato
solo nei modelli su scale superiori (i seguenti) alle quali la sua
semplicità è vantaggiosa in termini di costo computazionale [FIT60]
[FIT61], mentre è troppo riduttiva a questo livello.
• Modelli di tessuti. Lo scopo di questi modelli è caratterizzare la
propagazione dell’eccitazione nel miocardio in condizioni normali o in
presenza di sorgenti esterne, come nel caso della defibrillazione (la cui
applicazione - benché spesso efficace - è tuttora empirica, non essendo
noto il meccanismo con il quale agisce sull’attività cardiaca).
Tipicamente i modelli di questa classe sono formati dall’associazione di
un modello di propagazione ed un modello di membrana [COU96]
[KEE98]. Le ultime ricerche sono conformi alla concezione bidominio
5
del tessuto cardiaco, la quale afferma l’esistenza di due regioni che si
compenetrano: quella intracellulare e quella interstiziale. La differenza
tra i potenziali che si stabiliscono al loro interno ci dà il potenziale di
membrana. Un’altra caratteristica che ultimamente si cerca di introdurre
è l’anisotropia. L’anisotropia si manifesta in due aspetti diversi, ma
legati: la geometria delle fibre e la conducibilità variabile in funzione
della direzione di propagazione.
• Modelli integrali del cuore. Come dice il nome stesso questi modelli si
riferiscono all’intero volume cardiaco o a sue porzioni rilevanti, come i
ventricoli. Questi modelli tengono accuratamente conto della geometria
e dell’anatomia del cuore, incorporano un meccanismo per la
propagazione dell’attivazione nell’intero volume e mirano a riprodurre i
segnali di un elettrocardiogramma di superficie, previa inclusione di un
modello del torace umano.
Dato che in questo caso l’interesse è focalizzato su fenomeni a scala
globale, molti dei modelli di questa classe trascurano la dinamica
cellulare e per descrivere la propagazione dell’eccitazione si affidano ad
approssimazioni quali il principio di Huygens. Almeno questo è ciò che
è tradizionalmente accaduto. Ultimamente, infatti, si cerca di includere
in questi modelli i meccanismi che legano il potenziale di membrana
alle correnti ioniche [ALI96] [BER96] [PAN98]. Ovviamente un
modello così arricchito si comporta in maniera realistica in presenza di
una serie di situazioni importanti. L’inconveniente è un aumento della
complessità. Nel 1984 Barr e Plonsey [BAR84] espressero un giudizio
piuttosto scoraggiante sull’uso di un modello così complesso. Essi
stimarono che la simulazione di un ciclo d’attivazione per una schiera
bidimensionale di 106 nodi avrebbe richiesto 3000 anni. Ovviamente
negli anni il tempo necessario si è ridotto di parecchi ordini di
6
grandezza, ma ancora oggi è necessario accettare delle semplificazioni
per riportare ad un livello clinicamente pratico la quantità di risorse di
calcolo necessarie alla simulazione di modelli sofisticati del cuore.
Come ho già detto, data l’importanza pratica della questione, molte
ricerche mirano allo sviluppo di modelli integrali che mantengano le
proprietà importanti dei modelli a scala minore. I modelli ottenuti
consistono in gran parte dei casi in sistemi di PDE accoppiate, la cui
soluzione viene ottenuta per via numerica. Data la semplicità
d’implementazione, nella maggioranza dei casi viene adottato un metodo
alle differenze finite, spesso esplicito (Eulero). Questo allo scopo di evitare
di dover risolvere grossi sistemi di equazioni lineari, come nel caso dei
metodi semiimpliciti, o non lineari, per quelli impliciti. Purtroppo, dati i
vincoli imposti per ottenere la stabilità, l’utilizzo di un metodo esplicito
impone dei compromessi tra la durata della simulazione (dipendente dalla
discretizzazione temporale) e la risoluzione spaziale. Comunque proprio la
loro semplicità e la rapidità di realizzazione hanno fatto preferire questi
metodi rispetto ad alternative quali i metodi agli elementi finiti o quelli a
decomposizione di dominio (domain decomposition).
Ovviamente la ricerca continua per migliorare, soprattutto da un punto
di vista computazionale, le prestazioni delle implementazioni (si veda
[KEE98] [HAR97] [SAL98] solo per citarne alcuni).
L’approccio numerico soffre comunque della natura distribuita della
formulazione tramite PDE dei modelli. Infatti, il nostro è un problema di
evoluzione temporale su uno o più domini estesi nello spazio e quindi si
intuisce come mai gli studi sui modelli integrali vengano normalmente
condotti con l’ausilio di elaboratori massicciamente paralleli.
7
Negli ultimi anni si è assistito alla nascita e al diffondersi della
tecnologia delle reti neurali cellulari (CNN), circuiti analogici che al pari di
altri tipi di reti neurali realizzano il paradigma dell’elaborazione parallela e
distribuita, ma che a loro differenza presentano connessioni solo a livello
locale, il che permette una più semplice realizzazione hardware.
Sin da subito tra CNN e PDE si è ravvisata una similitudine che è stata
approfondita e sfruttata per usare questi dispositivi come alternativa alla
soluzione numerica di svariati problemi di PDE. Tipicamente si è affermato
che questi circuiti hanno rispetto agli approcci tradizionali il vantaggio di
operare in maniera intrinsecamente parallela e senza necessità di
discretizzazione temporale. Oltre a questo, trattandosi di circuiti analogici,
elaborano segnali i cui valori variano in intervalli del campo dei numeri
reali e quindi è assente l’errore di arrotondamento, inevitabile in ogni
elaborazione su calcolatore.
Alla luce di quanto detto, questo lavoro si propone di presentare le
CNN come alternativa agli strumenti solitamente usati per la simulazione
dei modelli cardiaci. Nonostante l’interesse per le CNN sia crescente, al
pari del loro campo d’applicazione, l’approccio qui esposto non è stato
finora tentato. E’ sembrato interessante quindi esaminare questo paradigma
di modellazione pensando che le sue potenzialità lo potrebbero rendere
un’importante alternativa all’uso di elaborazioni numeriche su calcolatori
paralleli.
Questo lavoro intende porsi al livello dei modelli di tessuti (si veda la
classificazione precedente), ma presentando alcune considerazioni utili al
fine dell’ottenimento di un modello integrale. Nel seguito si vedrà come è
possibile ottenere tramite CNN un modello dell’attività elettrica
ventricolare in un bidominio anisotropo e non omogeneo. Inoltre sarà presa
in esame l’influenza della geometria fibrosa sul tipo di maschere di
8
connessione usate per approssimare gli operatori differenziali, che saranno
sintetizzate tenendo conto di un nuovo metodo di approssimazione alle
differenze finite [SAL98].
Per quanto riguarda la modellazione della dinamica cellulare, saranno
descritti i casi dell’implementazione dei modelli di FitzHugh-Nagumo e di
Beeler-Reuter. Questi sono stati scelti non perché siano gli unici validi (la
validità dei modelli cellulari richiede considerazioni di fisiologia al di fuori
del campo di questo lavoro), ma semplicemente per illustrare la possibilità
di includere i modelli di membrana presentando un caso semplice ed uno
più complesso, ma direttamente legato alla fisiologia cellulare. Tra l’altro,
fin’ora non esistono implementazioni circuitali del modello di Beeler-
Reuter ( a differenza del modello di FitzHugh-Nagumo).
Infine i modelli di CNN proposti saranno utilizzati per delle
simulazioni al calcolatore di alcuni casi fisiologicamente importanti.
9
1. I MODELLI CARDIACI
Nel miocardio possiamo individuare due regioni coesistenti nello stesso
spazio: il dominio intracellulare e il dominio interstiziale, rispettivamente
all’interno della membrana cellulare e negli spazi tra le cellule. I due
domini sono ovunque separati dalla membrana cellulare e le loro proprietà
sono ricavate facendo delle medie sul volume da essi occupato.
Nonostante si tratti di domini non connessi, l’approssimazione bidominio ci
consente di trattare la propagazione dell’attivazione elettrica come se
questa avesse luogo in un mezzo continuo o, meglio, un cosiddetto sincizio.
Le cellule del miocardio, come tutte le cellule muscolari, hanno una
struttura fibrosa e quindi si estendono in una direzione molto più che nelle
altre; in tal modo è possibile associare ad ogni cellula una direzione
longitudinale. La conseguenza di questa particolare morfologia è
l’anisotropia della conducibilità elettrica. Per descrivere in termini
matematici questa proprietà associamo ad ogni punto r del miocardio un
tensore di conducibilità locale D(L)(r) ed un versore al (r) parallelo alla fibra
in r. Inoltre sia T(r) = [t1 (r), t2 (r), t3 (r)] una matrice le cui colonne sono i
versori di un riferimento cartesiano ortonormale e t3 (r) sia parallelo ad al
(r). Date le conducibilità g(L)l(r) e g(L)t(r) rispettivamente nella direzione
longitudinale e trasversale alla fibra, il tensore di conducibilità locale (in tre
dimensioni) è definito come una matrice diagonale tale che :
(1) D(L)(r) = diag [g(L)t(r), g(L)t(r), g(L)l(r)]
Il tensore di conducibilità globale è invece espresso da:
10
(2) D = TD(L)TT
Data la precedente espressione, D(r) è una matrice simmetrica e definita
positiva. Inoltre nel caso di modello bidominio abbiamo un tensore di
conducibilità per ciascun dominio, ognuno calcolato come sopra.
A causa delle reazioni chimiche che si svolgono all’interno delle cellule
e degli scambi ionici tra di esse, nel miocardio si stabiliscono due
potenziali elettrici, che indicheremo con Vi nel dominio intracellulare e con
Vo [mV] nel dominio interstiziale. Se indichiamo con ji e jo [mA⋅cm-2]
rispettivamente la densità di corrente intracellulare ed interstiziale,
otteniamo:
(3) oooiii VV ∇−=∇−= DjDj ,
da questa, per la legge di conservazione della carica elettrica espressa in
forma locale, consegue che:
(4) ⇒+−=+⋅∇ )],(),([)( tItI
sosioi rrjj
)],(),([ tItIVV
sosiooii rr +−=∇⋅∇+∇⋅∇⇒ DD
avendo indicato con Isi ed Iso [mA⋅cm-3] rispettivamente la stimolazione
dalla regione intracellulare e dalla regione interstiziale. Normalmente
si assume Isi = -Iso.
La membrana cellulare presenta delle strutture, dette canali ionici, che
regolano il passaggio in un senso o nell’altro di parecchie specie ioniche.
La somma delle correnti ioniche e della corrente dovuta alla presenza di
una capacità sulla membrana fornisce una corrente che è possibile mettere
11
in relazione con il potenziale di membrana, dato dalla differenza tra il
potenziale intracellulare e quello interstiziale:
(5) ),,(),( tItuj
t
uCI soionmm r−
+
∂
∂
= c
(5.a) u = Vi - Vo, v ≡ Vo
Nell’espressione precedente Cm [mF⋅cm-2] è la capacità di membrana, c
[cm-1] il rapporto superficie/volume, jion(u,t) [mA⋅cm-2]è la somma delle
correnti ioniche, mentre Im [mA⋅cm-3] è ovviamente la corrente
transmembrana. Per motivi di continuità tale corrente deve essere uguale
alla divergenza della densità di corrente interstiziale e quindi all’opposto
della divergenza della corrente intracellulare. Considerando la (5.a) e
introducendo il tensore Doi = Di + Do, otteniamo:
,
,
vVu
VVVV
oiooii
oioooiii
∇⋅−∇=∇⋅−∇=∇⋅∇
∇⋅−∇=∇⋅∇+∇⋅∇−∇⋅∇
DDD
DDDD
ed infine:
(6) ,0),(),( =∇⋅∇+−
+
∂
∂
vtItuj
t
uC osoionm Drc
(7) 0=∇⋅∇+∇⋅∇ vu oii DD
Il termine jion è stato precedentemente espresso come funzione del
potenziale di membrana e del tempo. In realtà la sua dipendenza dal tempo
non è diretta, ma sussiste a causa dell’esistenza delle cosiddette variabili di
12
attivazione/inattivazione, che rappresentano la maggiore o minore apertura
dei canali ionici. Tali variabili possono essere raggruppate nel vettore di
stato w ∈ RQ, Q ≥ 1, detto vettore di attivazione/inattivazione, la cui
evoluzione sarà a sua volta determinata da un’equazione differenziale
ordinaria. La forma di questa equazione e della funzione jion(u,w)
dipenderanno dal particolare modello cellulare scelto. In definitiva, il
modello bidominio nel caso generale consiste del seguente sistema di
equazioni differenziali:
(8)
[ ]
=
=∇⋅∇+∇⋅∇
+
∇⋅∇+
−=
∂
∂
),(
,0
,),(
),(1
wGw
DD
wD
u
dt
d
vu
uj
vtI
Ct
u
oii
ion
oso
m c
r
in cui G(u,w) ∈ RQ.
Il modello di Beeler-Reuter ed i suoi derivati
Come abbiamo già detto, i modelli di tessuti sono composti
dall’associazione di un modello di propagazione (che possiamo identificare
con il modello bidominio) e di un modello di membrana. I modelli di
membrana non fanno altro che mettere in relazione la tensione tra le sue
facce e la corrente che la attraversa. In condizioni fisiologiche tale corrente
è portata da ioni (Na+, K+, ioni organici negativi, ecc.), che attraversano la
membrana nelle due direzioni sotto la spinta di fattori quali la diversa
concentrazione all’interno ed all’esterno della cellula, la presenza di campi
13
elettrici e i meccanismi di pompaggio all’interno delle cellule.
L’andamento dei flussi ionici è regolato da strutture dette canali ionici che
modulano selettivamente la permeabilità della membrana di ogni tipo di
cellula rispetto ai vari elettroliti. La descrizione di tale attività nelle cellule
nervose è stata introdotta con il modello di Hodgkin e Huxley [HOD52] da
cui derivano i modelli ionici dell’attività elettrica ventricolare.
La maggiore o minore apertura dei canali ionici (e quindi la facilità con
cui avviene il passaggio degli ioni) viene quantificata per mezzo delle
variabili di attivazione-disattivazione. A causa della selettività della
membrana su di essa si stabilisce una differenza di potenziale poiché le
varie specie ioniche assumono concentrazioni diverse all’interno ed
all’esterno della cellula. Nello stato di riposo la permeabilità rispetto agli
ioni K+ e Cl- è massima e questi tendono a fluire rispettivamente verso
l’esterno e verso l’interno. Il processo di diffusione genera però un campo
elettrico che si oppone alla tendenza a diffondere nelle due direzioni.
All’equilibrio queste due forze si compensano a vicenda.
Durante il ciclo che dallo stato di riposo conduce alla generazione di un
potenziale d’azione, alla fase di recupero e quindi nuovamente
all’equilibrio, i canali ionici si aprono e si chiudono in funzione della
tensione e delle concentrazioni, producendo in diversi momenti del ciclo
impulsi di correnti ioniche di varia durata ed ampiezza. In ogni momento di
questo ciclo l’attività dei canali ionici è caratterizzata da una costante di
tempo ts con la quale ogni variabile di attivazione-inattivazione s tende ad
un valore asintotico s�. Ne consegue che possiamo descrivere molto
semplicemente la dinamica della variabile s con la seguente equazione
differenziale ordinaria:
14
(9)
s
ss
dt
ds
t
−
= ∞
Nella (9), apparentemente molto semplice, le grandezze s� e ts sono
funzioni del potenziale di membrana u. Infatti introducendo le cosiddette
costanti cinetiche (rate constants) as(u) e bs(u) è possibile esprimere tale
dipendenza nel seguente modo1:
(10) ts(u) =
)()(
1
uu ss ba +
,
s�(u) =
)()(
)(
uu
u
ss
s
ba
a
+
Le funzioni as(u) e bs(u) variano secondo il modello ionico a cui
appartengono e secondo la variabile di attivazione a cui si riferiscono.
Nell’ambito dei modelli ionici delle cellule ventricolari, quello di
Beeler-Reuter [BEE77] (successivamente indicato come BR) è
probabilmente il più citato e, comunque, quello da cui sono derivati gli altri
modelli, quali quello di Drouhard-Roberge [DRO86], quello di Ebihara-
Johnson [EBI80] e quello di Luo-Rudy [LUO91].
La densita di corrente ionica che attraversa una porzione di membrana
cellulare di 1 cm2 è, secondo BR:
(11)
sNaxKion iiiij +++= 11
1
Date queste definizioni l’equazione di stato delle generica variabile s viene spesso scritta anche così:
ds/dt = a – (a + b)s
15
Le correnti a secondo membro della (11) sono funzioni del potenziale di
membrana u, delle variabili di attivazione d,f,j,h,m,x e della concentrazione
intracellulare degli ioni calcio [Ca]i. Le espressioni di tali correnti sono le
seguenti:
iK1(u) =
−
+
+
+
−
+−++
+
)23(04.0)53(04.0)53(08.0
)85(04.0
1
23
2.0
1
435.0
uuu
u
e
u
ee
e
ix(u,x) =
[ ]
x
35)0.04(u
77)0.04(u
e
1e
8.0
+
+ −
iNa(u,m,h,j) = ( )( )NaNaCNa Eughjmg −+⋅ 3
is(u,d,f,[Ca]i) = )]]([[ iss CaEudfg −⋅
in cui le costanti hanno i seguenti valori:
09.050003.0,4 ====
sNaNaCNa gEgg
Le conduttanze sono in mS/cm2, le tensioni in mV, le correnti in mA/cm2.
La tensione Es ha la seguente dipendenza dalla concentrazione [Ca]i
(misurata in mole/l):
Es = -82.3 – 13.0287 ln[Ca]i
A sua volta la concentrazione [Ca]i varia secondo la seguente legge:
16
)][10(07.010
][ 77
is
i Cai
dt
Cad
−+⋅−= −−
Le restanti variabili sono governate da equazioni di stato della stessa
forma della (9), tenuto conto che le costanti cinetiche hanno la seguente
espressione generica (valida per a o b):
7
)(
54
)(
1
36
32 )(
Pe
PuPeP
PuP
PuP
+
++
=
+
+
a
Ogni costante cinetica ha i suoi rispettivi parametri P1, P2, … P7.
Per quanto riguarda gli altri modelli ionici citati all’inizio del paragrafo,
il modello di Luo-Rudy tiene in conto l’azione di diverse correnti e di
ulteriori concentrazioni ioniche per un totale di 14 variabili di stato
(potenziale interstiziale compreso) ed è molto complesso, per cui sembra
opportuno ometterlo; nei modelli di Ebihara-Johnson e di Drouhard-
Roberge, invece, si considera la sola corrente dovuta agli ioni Na, che è
quella principalmente responsabile della generazione di potenziali
d’azione. Di seguito sono riportate le equazioni relative al modello di
Drouhard-Roberge:
iNa(u,m,h) = ( )NaNa Euhmg −⋅ 3
)75.39(085.0
)65.42(22.
437.1
1
)65.42(09.0 +−
+−
=
−
+
= umum ee
u
ba
)5.20(095.
)65.79(193.0
1
7.1
1.0
+−
+−
+
==
uh
u
h
e
e ba
17
gNa = 15 mS/cm2, ENa = 40 mV.
Ricevi informazioni sui nostri servizi, sulle offerte e non perdere news e consigli su università e lavoro.
