Newer
Older
\chapter{Cardiac Electrophysiology Solver}
The CardiacEPSolver is used to model the electrophysiology of cardiac
tissue, specifically using the monodomain or bidomain model. These models are
continuum models and represent an average of the electrical activity over many
cells. The system is a reaction-diffusion system, with the reaction term
modeling the flow of current in and out of the cells using a separate set of
ODEs.
\subsection{Bidomain Model}
The Bidomain model is given by the following PDEs,
\begin{align*}
g_{ix}\frac{\partial^2 V_i}{\partial x^2} + g_{iy}\frac{\partial^2 V_i}{\partial y^2} &= \chi \left[ C_m \frac{\partial(V_i-V_e)}{\partial t} + G_m(V_i-V_e) \right] \\
g_{ex}\frac{\partial^2 V_e}{\partial x^2} + g_{ey}\frac{\partial^2 V_e}{\partial y^2} &= -\chi\left[ C_m \frac{\partial(V_i-V_e)}{\partial t} + G_m(V_i-V_e) \right].
\end{align*}
However, when solving numerically, one often rewrites these equations in terms
of the transmembrane potential and extracellular potential,
\begin{align*}
\chi \left[ C_m \frac{\partial V_m}{\partial t} + J_{ion} \right] &= g_{ex}\frac{\partial^2 V_e}{\partial x^2} + g_{ey}\frac{\partial^2 V_e}{\partial y^2}\\
(g_{ix} + g_{ex})\frac{\partial^2 V_e}{\partial x^2} + (g_{iy} + g_{ey})\frac{\partial^2 V_e}{\partial y^2}
&= -g_{ix} \frac{\partial^2 V_m}{\partial x^2} - g_{iy} \frac{\partial^2 V_m}{\partial y^2}
\end{align*}
\subsection{Monodomain Model}
In the case where the intracellular and extracellular conductivities are
proportional, that is $g_{ix} = kg_{ex}$ for some $k$,
then the above two PDEs can be reduced to a single PDE:
\begin{align*}
\chi\left[ C_m \frac{\partial V_m}{\partial t} + J_{ion} \right] &= \nabla \cdot (\sigma \nabla V_m)
\end{align*}
\subsection{Cell Models}
The action potential of a cardiac cell can be modelled at either a biophysical
level of detail, including a number of transmembrane currents, or as a
phenomenological model, to reproduce the features of the action potential, with
fewer variables. Each cell model will include a unique system of ODEs to
represent the gating variables of that model.
A number of ionic cell models are currently supported by the solver including:
\begin{itemize}
\item Courtemanche, Ramirez, Nattel, 1998
\item Luo, Rudy, 1991
\item ten Tusscher, Panfilov, 2006 (epicardial, endocardial and
mid-myocardial variants)
\end{itemize}
Phenomological cell models are also supported:
\begin{itemize}
\item Aliev-Panfilov
\item Fitzhugh-Nagumo
\end{itemize}
It is important to ensure that the units of the voltage and currents from the
cell model are consistent with the units expected by the tissue level solver
(monodomain/bidomain). We will show as an example the Courtemanche, Ramirez,
Nattel, 1998 human atrial model.
The monodomain equation:
\begin{align*}
\chi \left[ C_m \frac{\partial V_m}{\partial t} + J_{ion} \right] &= \nabla \cdot (\sigma \nabla V_m)
\end{align*}
\begin{lstlisting}[style=BashInputStyle]
CardiacEPSolver session.xml
\end{lstlisting}
\section{Session file configuration}
\subsection{Solver Info}
\item \inltt{Eqtype} Specifies the PDE system to solve. The following values
are supported:
\begin{itemize}
\item \inltt{Monodomain}: solve the monodomain equation.
\item \inltt{BidomainRoth}: solve the bidomain equations using the Roth
formulation.
\end{itemize}
\item \inltt{CellModel} Specifies the cell model to use. Available cell
models are
\begin{center}
\begin{tabular}{l|l|l|l}
Value & Description & No. of Var. & Ref. \\
\inltt{AlievPanfilov} & Phenomological & 1 & \cite{AlPa96} \\
\inltt{CourtemancheRamirezNattel98} & Human atrial & 20 & \cite{CoRaNa98} \\
\inltt{FitzHughNagumo} & & & \\
\inltt{Fox02} & & & \\
\inltt{LuoRudy91} & Mammalian ventricular & 7 & \cite{LuRu91} \\
\inltt{PanditGilesDemir03} & & & \\
\inltt{TenTusscher06} & Human ventricular & 18 & \cite{TuPa06} \\
\inltt{Winslow99} & & & \\
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
\end{tabular}
\end{center}
\item \inltt{Projection} Specifies the Galerkin projection type to use. Only
\inltt{Continuous} has been extensively tested.
\item \inltt{TimeIntegrationMethod} Specifies the time integration
scheme to use for advancing the PDE system. This must be an IMEX scheme.
Suitable choices are: \inltt{IMEXOrder1}, \inltt{IMEXOrder2},
\inltt{IMEXOrder3}, \inltt{IMEXdirk\_3\_4\_3}. The cell model state
variables are time advanced using Forward Euler for the ion concentrations, and
Rush-Larsen for the cell model gating variables.
\item \inltt{DiffusionAdvancement} Specifies whether the diffusion is
handled implicitly or explicitly in the time integration scheme. The current
code only supports \inltt{Implicit} integration of the diffusion term. The
cell model is always integrated explicitly.
\end{itemize}
\subsection{Parameters}
The following parameters can be specified in the \inltt{PARAMETERS} section of
the session file. Example values are taken from \cite{NiKeBeBeBe11}.
\begin{itemize}
\item \inltt{Chi} sets the surface-to-volume ratio (Units:
$\mathrm{mm}^{-1}$).\\ Example: $\chi= 140 \mathrm{mm}^{-1}$
\item \inltt{Cm} sets the specific membrane capacitance (Units:
$\mathrm{\mu F\,mm}^{-2}$).\\ Example: $C_m= 0.01 \mathrm{\mu F\,mm}^{-2}$
\item \inltt{Substeps} sets the number of substeps taken in time
integrating the cell model for each PDE timestep.\\ Example: 4
\item \inltt{d\_min}, \inltt{d\_max}, \inltt{o\_min}, \inltt{o\_max}
specifies a bijective map to assign conductivity values $\sigma$ to
intensity values $\mu$ when using the \inltt{IsotropicConductivity}
function. The intensity map is first thresholded to the range
$[d_{\min},d_{\max}]$ and then the conductivity is calculated as
\begin{align*}
\sigma = \frac{o_{\max} - o_{\min}}{d_{\max}-d_{\min}} (1 - \mu) + o_{\min}
\end{align*}
\end{itemize}
\subsection{Functions}
The following functions can be specified inside the \inltt{CONDITIONS} section
of the session file. If both are specified, the effect is multiplicative.
Example values are taken from \cite{NiKeBeBeBe11}.
\begin{itemize}
\item \inltt{IsotropicConductivity} specifies the conductivity
$\sigma$ of the tissue. \\
Example: $\sigma =0.13341 \ \mathrm{mS\,mm}^{-1}$, based on $\sigma
= \frac{\sigma_i \sigma_e} {\sigma_i + \sigma_e}, \ \ \sigma_i=0.17, \sigma_e=0.62 \mathrm{mS\,mm}^{-1} $
The variable name to use is \inltt{intensity} since the conductivity may be
derived from late-Gadolinium enhanced MRA imaging. Example specifications are
\begin{lstlisting}[style=XmlStyle]
<E VAR="intensity" VALUE="0.13341" />
<F VAR="intensity" FILE="scarmap.con" />
\end{lstlisting}
where \inltt{scarmap.con} is a Nektar++ field file containing a variable
\inltt{intensity} describing the conductivity across the domain.
\item \inltt{AnisotropicConductivity} specifies the conductivity
$\mathbf{\sigma}$ of the tissue.
\end{itemize}
\subsection{Filters}
\begin{itemize}
\item \inltt{CellHistoryPoints} writes all cell model states over time at
fixed points. Can be used along with the \inltt{HistoryPoints} filter to
record all variables at specific points during a simulation.
\begin{lstlisting}[style=XmlStyle]
<FILTER TYPE="CellHistoryPoints">
<PARAM NAME="OutputFile">crn.his</PARAM>
<PARAM NAME="OutputFrequency">1</PARAM>
<PARAM NAME="Points">
0.00 0.0 0.0
</PARAM>
</FILTER>
\end{lstlisting}
\begin{itemize}
\item \inltt{OutputFile} specifies the filename to write history data to.
\item \inltt{OutputFrequency} specifies the number of steps between successive outputs.
\item \inltt{Points} lists coordinates at which history data is to be recorded.
\end{itemize}
\item \inltt{CheckpointCellModel} checkpoints the cell model. Can be
used along with the \inltt{Checkpoint} filter to record complete simulation
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
state and regular intervals.
\begin{lstlisting}[style=XmlStyle]
<FILTER TYPE="CheckpointCellModel">
<PARAM NAME="OutputFile"> session </PARAM>
<PARAM NAME="OutputFrequency"> 1 </PARAM>
</FILTER>
\end{lstlisting}
\begin{itemize}
\item \inltt{OutputFile} (optional) specifies the base filename to use.
If not specified, the session name is used. Checkpoint files are suffixed with the process ID and the extension `.chk`.
\item \inltt{OutputFrequency} specifies the number of timesteps between
checkpoints.
\end{itemize}
\item \inltt{Electrogram} Computes virtual unipolar electrograms at a
prescribed set of points.
\begin{lstlisting}[style=XmlStyle]
<FILTER TYPE="Electrogram">
<PARAM NAME="OutputFile"> session </PARAM>
<PARAM NAME="OutputFrequency"> 1 </PARAM>
<PARAM NAME="Points">
0.0 0.5 0.7
1.0 0.5 0.7
2.0 0.5 0.7
</PARAM>
</FILTER>
\end{lstlisting}
\begin{itemize}
\item \inltt{OutputFile} (optional) specifies the base filename to use. If
not specified, the session name is used. The extension `.ecg` is appended if not already specified.
\item \inltt{OutputFrequency} specifies the number of resolution of the
electrogram data.
\item \inltt{Points} specifies a list of coordinates at which electrograms
are desired. \emph{They must not lie within the domain.}
\end{itemize}
\item \inltt{Benchmark} Records spatially distributed event times for
activation and repolarisation (recovert) during a simulation, for
undertaking benchmark test problems.
\begin{lstlisting}[style=XmlStyle]
<FILTER TYPE="Benchmark">
<PARAM NAME="ThresholdValue"> -40.0 </PARAM>
<PARAM NAME="InitialValue"> 0.0 </PARAM>
<PARAM NAME="OutputFile"> benchmark </PARAM>
<PARAM NAME="StartTime"> 0.0 </PARAM>
</FILTER>
\end{lstlisting}
\begin{itemize}
\item \inltt{ThresholdValue} specifies the value above which tissue is
considered to be depolarised and below which is considered
repolarised.
\item \inltt{InitialValue} specifies the initial value of the
activation or repolarisation time map.
\item \inltt{OutputFile} specifies the base filename of activation and
repolarisation maps output from the filter. This name is appended
with the index of the event and the suffix `.fld`.
\item \inltt{StartTime} (optional) specifies the simulation time at
which to start detecting events.
\end{itemize}
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
\end{itemize}
\subsection{Stimuli}
Electrophysiological propagaion is initiated through the stimulus current
$I_{\mathrm{ion}}$. The \inltt{STIMULI} section describes one or more regions of
stimulus and the time-dependent protocol with which they are applied.
\begin{lstlisting}[style=XmlStyle]
<STIMULI>
...
</STIMULI>
\end{lstlisting}
A number of stimulus types are available
\subsubsection{Stimulus types}
\begin{itemize}
\item \inltt{StimulusRect} stimulates a cuboid-shaped region of the domain,
specified by two coordinates $(x_1,y_1,z_1)$ and $(x_2,y_2,z_2)$.
An additional parameter specifies the "smoothness" of the boundaries of the
region; higher values produce a sharper boundary. Finally, the maximum
strength of the stimulus current is specified in $\mu \mathrm{A} / \mathrm{mm}^3$
\begin{lstlisting}[style=XmlStyle]
<STIMULUS TYPE="StimulusRect" ID="0">
<p_x1> -15.24 </p_x1>
<p_y1> 14.02 </p_y1>
<p_z1> 6.87 </p_z1>
<p_x2> 12.23 </p_x2>
<p_y2> 16.56 </p_y2>
<p_z2> 8.88 </p_z2>
<p_is> 100.00 </p_is>
<p_strength> 50.0 </p_strength>
</STIMULUS>
\end{lstlisting}
\item \inltt{StimulusCirc} stimulates a spherical region of the domain, as
specified by a centre and radius. The smoothness and strength parameters are also specified as for `StimulusRect`.
\begin{lstlisting}[style=XmlStyle]
<STIMULUS TYPE="StimulusCirc" ID="0">
<p_x1> -15.24 </p_x1>
<p_y1> 14.02 </p_y1>
<p_z1> 6.87 </p_z1>
<p_r1> 12.23 </p_r1>
<p_is> 100.00 </p_is>
<p_strength> 50.0 </p_strength>
</STIMULUS>
\end{lstlisting}
\end{itemize}
\subsubsection{Protocols}
A protocol specifies the time-dependent function indicating the strength of the
stimulus and one such \inltt{PROTOCOL} section should be included within each
\inltt{STIMULUS}. This can be expressed as one of:
\begin{itemize}
\item \inltt{ProtocolSingle} a single stimulus is applied at a given start
time and for a given duration
\begin{lstlisting}[style=XmlStyle]
<PROTOCOL TYPE="ProtocolSingle">
<START> 0.0 </START>
<DURATION> 2.0 </DURATION>
</PROTOCOL>
\end{lstlisting}
\item \inltt{ProtocolS1} a train of pulses of fixed duration applied at a
given start time and with a given cycle length.
\begin{lstlisting}[style=XmlStyle]
<PROTOCOL TYPE="ProtocolS1">
<START> 0.0 </START>
<DURATION> 2.0 </DURATION>
<S1CYCLELENGTH> 300.0 </S1CYCLELENGTH>
<NUM_S1> 5 </NUM_S1>
</PROTOCOL>
\end{lstlisting}
\item \inltt{ProtocolS1S2} same as `ProtocolS1` except with an additional
single pulse applied at a different cycle length at the end of the train of S1 pulses.
\begin{lstlisting}[style=XmlStyle]
<PROTOCOL TYPE="ProtocolS1S2">
<START> 0.0 </START>
<DURATION> 2.0 </DURATION>
<S1CYCLELENGTH> 300.0 </S1CYCLELENGTH>
<NUM_S1> 5 </NUM_S1>
<S2CYCLELENGTH> 100.0 </S2CYCLELENGTH>
</PROTOCOL>
\end{lstlisting}
\end{itemize}