Transcript Modeling and Analysis of Atrial Fibrillation (ppt)

Modeling and Analysis
of Atrial Fibrillation
Radu Grosu
SUNY at Stony Brook
Joint work with
Ezio Bartocci, Flavio Fenton,
Robert Gilmour, James Glimm and Scott A. Smolka
Emergent Behavior in Heart Cells
EKG
Surface
Arrhythmia afflicts more than 3 million Americans alone
Modeling
CellExcite and Simulation
Tissue Modeling: Triangular Lattice
Communication by diffusion
CellExcite and Simulation
Tissue Modeling: Square Lattice
Communication by diffusion
Single Cell Reaction: Action Potential
Schematic Action Potential
Membrane’s AP depends on:
• Stimulus (voltage or current):
• Cell’s state
AP has nonlinear behavior!
• Reaction diffusion system:
u
t
 R (u )   ( D  u )
voltage
– External / Neighboring cells
Threshold
failed initiation
Resting potential
time
Behavior
In time
Reaction
Diffusion
Frequency Response
APD90: AP > 10% APm
DI90: AP < 10% APm
BCL: APD + DI
Frequency Response
APD90: AP > 10% APm
DI90: AP < 10% APm
BCL: APD + DI
S1-S2 Protocol: (i) obtain stable S1; (ii) deliver S2 with shorter DI
Frequency Response
APD90: AP > 10% APm
DI90: AP < 10% APm
BCL: APD + DI
S1S2 Protocol: (i) obtain stable S1; (ii) deliver S2 with shorter DI
Restitution curve: plot APD90/DI90 relation for different BCLs
Existing Models
• Detailed ionic models:
–
–
–
–
Luo and Rudi: 14 variables
Tusher, Noble2 and Panfilov: 17 variables
Priebe and Beuckelman: 22 variables
Iyer, Mazhari and Winslow: 67 variables
• Approximate models:
– Cornell: 3 or 4 variables
– SUNYSB: 2 or 3 variable
Stony Brook’s Cycle-Linear Model
Objectives
• Learn a minimal mode-linear HA model:
– This should facilitate analysis
• Learn the model directly from data:
– Empirical rather than rational approach
• Use a well established model as the “myocyte”:
– Luo-Rudi II dynamic cardiac model
HA Identification for the Luo-Rudi Model
(with P. Ye, E. Entcheva and S. Mitra)
• Training set: for simplicity 25 APs generated from the LRd
– BCL1 + DI2: from 160ms to 400 ms in 10ms intervals
• Stimulus: step with amplitude -80μA/cm2, duration 0.6ms
• Error margin: within ±2mV of the Luo-Rudi model
• Test set: 25 APs from 165ms to 405ms in 10ms intervals
Action Potential (AP) Phases
Stimulated
Identifying a Mode-Linear HA for One AP
u  E
u  P
u  F
Stimulated
u  U
son
u  R
soff  u  U
Identifying the Switching for one AP
Null Pts: discrete 1st Order deriv.
Infl. Pts: discrete 2nd Order deriv.
Thresholds: Null Pts and Infl. Pts
Segments: Between Seg. Pts
Problem: too many Infl. Pts
Problem: too many segments?
Identifying the Switching for one AP
Null Pts: discrete 1st Order deriv.
Infl. Pts: discrete 2nd Order deriv.
Thresholds: Null Pts and Infl. Pts
Segments: Between Seg. Pts
Problem: too many Infl. Pts
Problem: too many segments?
Solution: use a low-pass filter
-Moving average and spline LPF: not satisfactory
-Designed our own: remove pts within trains of inflection points
Identifying the Switching for all AP
Problem: somewhat different inflection points
Identifying the Switching for all AP
Solution: align, move up/down and remove inflection points
- Confirmed by higher resolution samples
Identifying the HA Dynamics for One AP
Modified Prony Method
u  E
u  P /
x i  ai
 ao
v xo V
P
u  x&i  x&o  I s
x&v
bV
xF
i
i i
x&o  b o x o
Stimulated
u  U
son
u  R
soff  u  U
Summarizing all HA
u   E (d i )
u   P (di ) /
x i  ai (d i )
x o  a o (d
( d i ))
u
P
i
u  x&i  x&o  I s
x&iub i( dF i(d
) xi i)
x&o  b o ( d i ) x o
Stimulated
u  U (d i )
u   R (d i )
/t  0
son / d i  t
soff  u  U (di )
Finding Parameter Dependence on DI
Solution: apply mProny once again on each of the 25 points
Cycle Linear
Summarizing all HA
u   E (d i )
u   P (di ) /
x i  ai (d i )
x o  a o (d
( d i ))
u
P
i
b i (d i )   i 1 e
 i1 di
b o (d i )   o 1 e
 o1 di
Stimulated
u  U (d i )
  i2e
 i 2 di
  o2e
 o 2 di
u  x&i  x&o  I s
x&iub i( dF i(d
) xi i)
x&o  b o ( d i ) x o
u   R (d i )
/t  0
son / d i  t
soff  u  U (di )
Frequency Response on Test Set
AP on test set: still within the accepted error margin
Restitution on test set: follows very well the nonlinear trend
Cornell’s Nonlinear Minimal Model
Objectives
• Learn a minimal nonlinear model:
– This should facilitate analysis
• Approximate the detailed ionic models:
– Rational rather than empirical approach
• Identify the parameters based on:
– Data generated by a detailed ionic model
– Experimental, in-vivo data
Switching Control



R (u , u s1 , u s 2 )  



u  u s1
0
u  u s1
u s 2  u s1
S ( k s (u  u s )) 
1 e
u  us2
1
 0

H (u  u s )  
 1
else
1
u  us
u s  0 .5
u  us
ks  16
 ks (u  u s )
Cornell’s Minimal Model
u
  ( D  u )  ( J fi  J si  J so )
J fi   H (u   v )(u   v )(u u  u )v / 
voltage
fi
Diffusion
input input
Slow output
Laplacia Fast Slow
current
currentcurrent
n
Cornell’s Minimal Model
u

JJ fifi 
J si 
Activation
Heaviside
Fast input
Threshol
 ( D  Slow
u )  ( Input
J fi  (step)
J si Slow
J so ) Output
GatePiecewise
Resistance
d
Gate
Gate
Nonlinear
Time Cst
 H (u


)(u


)(u

u
)v
/

v  )
v  
u )(u  u fi
Piecewise
H (u 
(u
)v
/

v
v
u
fi
Piecewise
Bilinear
Nonlinear
 H (u   w ) w s /  si
J so  (1  H (u   w )) (u  u o ) /  o  H (u   w ) /  so
Piecewise


 (1  H (u


))
(v

v
)
/


H
(u


)v
/

Sigmoid
v

v
v
v Linear


w&  (1  H (u(s-step)
  w ))(w   w ) /  w Nonlinear
H (u   w )w /  w
v
s&
 (S (2 k s (u  u s ))  s ) /  s
J fi   H (u   v )(u   v )(u u  u )v / 
fi
Time Constants and Infinity Values
v

 (1  H (u   v ))  v1  H (u   v )  v 2



s
 (1  H (u   w ))  s1  H (u   w )  s 2
Piecewise
Constant

 o  (1  H (u   o ))  o1  H (u   o )  o 2
w






 w   w 1  ( w 2   w 1 ) S (2 k w (u  u w ))
Sigmoidal
 so   so1  ( so 2   so1 ) S (2 k so (u  u so ))
w

v   (1  H (u   v ))
Piecewise
Linear
w   (1  H (u   o )) (1  u /  w  )  H (u   o ) w 
*
 so 
(1  H (u   o ))  o1
 H (u   o )  o 2
Single Cell Action Potential
Cornell’s Minimal Model
v  u
u&   ( D  u )  (u   v )(u u  u )v / 
fi
 ws / 
fi
 1 /  so

v&   v /  v

w&   w /  w
w  u  v
s&  (S (2 k s (u  u s ))  s ) /  s 2
u&   ( D  u )  w s /  si  1 /  so

o  u  w
v&   v /  v 2

w&   w /  w
u&   ( D  u )  u /  o 2
s&  (S (2 k s (u  u s ))  s ) /  s 2
v&   v /  v 2
u  o

w&  (w   w ) /  w 1
*
u  v
u  w

s&  (S (2 k s (u  u s ))  s ) /  s
u  o
u   v  0 .3
u   w  0. 13
u&   ( D  u )  u /  o1

u   o   v  0 .0 0 6

v&  (1  v ) /  v1

w&  (1  u /  w   w ) /  w
s&  (S (2 k s (u  u s ))  s ) /  s
Partition with Respect to v
v  vc
u  v
u   v  0 .3
u  w
u   w  0. 13
u  o
u   o   v  0 .0 0 6

Partition with Respect to v
v  vc
v  u
u&   ( D  u )  (u   v )(u u  u )v / 
fi
 ws / 
fi
 1 /  so

v&   v /  v

w&   w /  w
s&  (S (2 k s (u  u s ))  s ) /  s 2
(  v  u )  (v  v c )
u&   ( D  u )  w s / 
fi
 1 /  so

v&   v /  v

w&   w /  w
s&  (S (2 k s (u  u s ))  s ) /  s 2
u  v
u   v  0 .3
u  w
u   w  0. 13
u  o
u   o   v  0 .0 0 6

Superposed Action Potentials
HA for the Model
(  v  u )  (v  v c )
u&   ( D  u )  (u   v )(u u  u )v / 
u  v
fi
 ws / 
fi
 1 /  so
u  v


v&   v /  v
v  vc
w&   w /  w

v  vc

w  u  v
s&  (S (2 k s (u  u s ))  s ) /  s 2
(  v  u )  (v  v c )
u  v
u&   ( D  u )  w s /  si  1 /  so
u&   ( D  u )  w s / 

v&   v /  v 2
 1 /  so

v&   v /  v

w&   w /  w

w&   w /  w
u  w
s&  (S (2 k s (u  u s ))  s ) /  s 2
o  u  w
u&   ( D  u )  u /  o 2
s&  (S (2 k s (u  u s ))  s ) /  s 2
u  o

u  w
fi
v&   v /  v 2

w&  (w   w ) /  w 1
*
s&  (S (2 k s (u  u s ))  s ) /  s
u  o
u&   ( D  u )  u /  o1

v&  (1  v ) /  v1

w&  (1  u /  w   w ) /  w
u  o
s&  (S (2 k s (u  u s ))  s ) /  s
Analysis of Sigmoidal Switching






 w   w 1  ( w 2   w 1 ) S (2 k w (u  u w ))
 so   so1  ( so 2   so1 )S (2 k so (u  u so ))
s& 

(S (2 k s (u  u s ))  s ) /  s




 w  (1  H (u  u w )) w 1  H (u  u w ) w 2
s  (rs R (u ,  v )  s ) /  s
Superposed Action Potentials
Current HA of Cornell’s Model
(  v  u )  (v  v c )
u  v
u&   ( D  u )  (u   v )(u u  u )v / 

fi
 ws / 
fi
 1 /  so
u  v


v&   v /  v
v  vc
v  vc

w&   w /  w
w  u  v
s&  ((u   v ) / (2 rs u s )  s ) /  s 2
u&   ( D  u )  w s /  si  1 /  so
v&   v / 
u  v

v2
u&   ( D  u )  w s / 

 1 /  so



w&   w /  w

uw  u  w
u  uw
s&  ((u   v ) / (2 rs u s )  s ) /  s 2
u&   ( D  u )  u /  o 2
u  w
fi
v&   v /  v
u  w
w&   w /  w
s&   s /  s 2
(  v  u )  (v  v c )

 o  u  uw

v&   v /  v 2
w&  (w   w ) /  w 2
u&   ( D  u )  u /  o1
s&   s /  s1
v&  (1  v ) /  v1
*

u  o

w&  (w   w ) /  w 1
u&   ( D  u )  u /  o1
s&   s /  s1
v&  (1  v ) /  v1
*
u  u

w
u  o



w&  (1  u /  w   w ) /  w 1
u  o
s&   s /  s1
Analysis of 1/τso ?
 so   so1  ( so 2   so1 )S (2 k so (u  u so ))
J so  (1  H (u   w ))(u  u o )  H (u   w ) /  so
Cubic Approximation of 1/τso ?
 so   so1  ( so 2   so1 )S (2 k so (u  u so ))
J so  (1  H (u   w ))(u  u o )  H (u   w ) /  so
Superposed Action Potentials
Very sensitive!
Summary of Models
• Both models are nonlinear
– Stony Brook’s: Linear in each cycle
– Cornell’s: Nonlinear in specific modes
• Both models are deterministic
• Both models require identification
– Stony Brook’s: On a mode-linear basis
– Cornell’s: On an adiabatically approximated model
Modeling Challenges
• Identification of atrial models
– Preliminary work: Already started at Cornell
• Dealing with nonlinearity
– Analysis: New nonlinear techniques? Linear approx?
• Parameter mapping to physiological entities
– Diagnosis and therapy: To be done later on
Analysis
Atrial Fibrillation (Afib)
• A spatial-temporal property
– Has duration: it has to last for at least 8s
– Has space: it is chaotic spiral breakup
• Formally capturing Afib
– Multidisciplinary: CAV, Computer Vision, Fluid Dynamics
– Techniques: Scale space, curvature, curl, entropy, logic
Spatial Superposition
• Detection problem:
– Does a simulated tissue
contain a spiral ?
• Specification problem:
– Encode above property as a
logical formula?
– Can we learn the formula?
How? Use Spatial Abstraction
Superposition Quadtrees (SQTs)
 !m  {s,u,p,r}. p l (m ) = 1
p i (m ) =
1
4
p

4
ij
(m j )
j= 1
Abstract position and compute PMF p(m) ≡ P[D=m]
Linear Spatial-Superposition Logic
Syntax
Semantics
The Path to the Core of a Spiral
Root
1
1
1
Click the core to
determine the quadtree
1
2
2
3
3
4
1
2
3
2
3
4
2
3
4
4
4
Overview of Our Approach
Emerald: Learning LSSL Formula
Emerald: Bounded Model Checking
Curvature Analysis
N -
N orm al
T -
T angent
 T - C urvature
N
N
T
T
T
• Some properties of the curvature:
– The curvature of a straight line is identical to 0
– The curvature of a circle of radius R is constant
– Where the curve undergoes a tight turn, the curvature is large
• Measuring the curvature:
– Adapting Frontier Tool [Glimm et.al]: MPI code on Blue Gene
– Also corrects topological errors
Edge Detection
Scalar field
Front wave
Canny algorithm
Normal Vectors Computation
Compute the Gradient
Tangent Vectors Computation
Based on the Gradient
The Curl of the Tangent Field
Curl = infinitesimal rotation of a vector field (circulation density of a fluid)
Verification Setup
• Models are deterministic with one initial state:
– A spiral: induced with a specific protocol
• Verification becomes parameter estimation/synthesis:
– In normal tissue: no fibrillation possible
– Diseased tissue: brute force gives parameter bounds
– Parameter space search: increases accuracy
• Parameters are mapped to the ionic entities:
– Obtained mapping: used for diagnosis and therapy
Possible Collaborations
• Pancreatic cancer group:
– Spatial properties: also a reaction diffusion system
– Nonlinear models: approximation, diff. invariants, statistical MC
– Parameter estimation: information theory, statistical MC
• Aerospace / Automotive groups:
– Monitoring & Control: low energy defibrillation, stochastic HA
– Machine learning: of spatial temporal patterns