1. Introduction

Every heartbeat is triggered by a rapidly propagating wave of electrical action potentials (conduction velocities up to 1 m/s) that synchronizes the contractions of the various chambers of the heart. Cardiac tissue consists of interconnected cardiac cells which are organized in muscle fibres. The electrical conductivity is much higher along than across fibres, making the myocardium highly anisotropic for electrical wave propagation [1]. The complex pattern of fibre organization throughout the myocardial wall allows for efficient contractions. Abnormal propagation of electrical activity severely compromises the mechanical function of the heart and can result in life-threatening cardiac arrhythmias [2]. Several computational studies have highlighted the complex three-dimensional (3D) activation patterns during normal and abnormal cardiac electrical activity [3–5]. Yet, so far only limited experimental information has been gained on the 3D spatiotemporal dynamics of wave propagation in cardiac tissue.

Optical imaging of cardiac electrical activity using voltage-sensitive dyes (VSDs) has become the method of choice to investigate arrhythmias experimentally at the tissue or whole heart level [6,7]. VSDs can be introduced through coronary flow without significant tissue damage and bind to the cardiac cell membranes. They respond to changes in transmembrane potential by changes in excitation and fluorescence spectra, which allow monitoring the cell’s electrical activity. Although recent advances have been made towards 3D optical imaging of cardiac electrical activity [8–10], surface epi-fluorescence imaging remains the most widely used technique in cardiac optical imaging [6, 7].

 

Fig. 1. A schematic representation of a typical cardiac epi-fluorescence experiment in isolated guinea pig hearts. In this setup the epicardial surface of the heart is uniformly illuminated by a light source (e.g. a tungsten-halogen lamp), whose light has been bandpass filtered at the appropriate wavelengths for excitation of the VSD. Fluorescence optical signals of cardiac electrical activity are then recorded using a CCD camera from the same epicardial area at the appropriate wavelength. A typical optical action potential as recorded from a single image pixel is shown on the right.

Download Full Size | PDF

A typical cardiac epi-fluorescence imaging system is schematically depicted in Fig. 1. Several laboratories have reported recently that fluorescent images recorded from the epicardium using such technique, contain contributions from deeper myocardial layers and that the magnitude of these contributions decreases with depth [11–15]. This effect is due to the photon scattering and absorption in tissue and can yield contributions to the optical signals from layers up to 1 mm below the epicardium. These studies indicated that the rising phase of the optical action potential, commonly termed the optical upstroke, is affected most by photon scattering: it can be up to 15 ms in duration, whereas the underlying electrical upstroke in a single cell is typically about 1 to 2 ms long [11].

 

Fig. 2. Subsurface wave front orientation and V_F*. The green area caricaturizes the subsurface volume from which contributions are made to the optical signal in the pixel of interest (dashed arrows). Panel A shows isochrones (white lines).of a wave front propagating towards the epicardium, whereas Panel B shows the opposite situation. The optical signal recorded from the pixel of interest on the epicardial surface is shown on top, as well as V_F* (black circle).

Download Full Size | PDF

Understanding this phenomenon helped in the interpretation of experimental data and led to the development of a novel technique to infer the subsurface wave front orientation [16, 17]. This technique relies on the observation that the morphology of the optical upstroke correlates to the subsurface propagation direction of electrical waves. Specifically, the fractional level at which the fluorescent signal has maximal time derivative V_F* (see Fig. 1), is low for waves propagating away from the epicardial surface, but high for waves moving towards the imaged surface.

Figure 2 illustrates schematically how V_F* correlates to subsurface wave front orientation. Panel A shows isochrones of a wave front propagating towards the epicardium. As the wave approaches the epicardium, it will gradually invade the subsurface volume from which contributions are made to the optical signal in the pixel of interest (green area). This results in an initially slow increase in V_F, as these regions only contribute modestly to V_F. As the wave approaches the epicardium, more substantial contributions are made to V_F, resulting in a faster increase in V_F. Maximal fluorescence is reached when the wave breaks through on the epicardium. Hence, the relative level of maximal slope during upstroke, i.e. V_F*, is located in the upper half of the upstroke, close to its maximum. The opposite is true for a wave propagating away from the epicardium as illustrated in Panel B. Intermediate values of V_F* are found for waves propagating at different angles with respect to the imaged surface [16,17].

V_F* maps of surface activity provided the first tool to gain information on the three-dimensional nature of cardiac electrical activity utilizing VSDs. Simulated V_F* maps showed good agreement with experimental V_F* maps in various species [17,18]. In all these studies the photon transport model is essential in determining the actual subsurface wave front orientation. Utilizing the diffusion equation in thick slabs of cardiac tissue (up to 1 cm thick), a strong linear correlation between V_F* and the optically averaged subsurface wave orientation was shown. The photon diffusion approximation is valid in the high scattering regime for tissue samples at least one order of magnitude thicker than the photon attenuation length, which for typical VSDs lies in the range of 1–2 mm. It is therefore unclear how previous computational results can be applied to thinner cardiac preparations, such as the widely used guinea pig hearts, where the diffusion approximation is expected to break down. Specifically, the linear relationship between V_F* and the subsurface wave front orientation in such thin hearts remains to be shown.

The major goal of this study is to perform a comprehensive analysis of the sensitivity of the simulated epi-fluorescence cardiac action potential to the photon transport model in slabs of various thicknesses. We compare optical action potentials obtained using the photon diffusion equation to a Monte Carlo based model. Monte Carlo transport models are well-established in other applications, but their use in combination with computed cardiac activity is fairly limited and has not yet been validated against experimental data [10,19]. We focus on the optical upstroke duration and morphology, as these are the features of the optical action potential most affected by optical blurring. All computational results are compared to experimental data obtained from intact guinea pig hearts as a mean of validation. We also determine, for the various models and tissue thicknesses, the relationship between V_F* and the optically averaged subsurface wave front orientation. Finally, subsurface electrical activity is reconstructed form experimental epi-fluorescence recordings using the various photon transport models and the differences and similarities between the various models are discussed.

2. Methods

2.1 Hybrid electro-optical models

Optical signals of cardiac electrical activity can be simulated using hybrid electro-optical models, which can be applied to different types of acquisition modes, tissue geometries, and boundary conditions [8–10, 19–21]. These hybrid models consist of: (1) an electrophysiological model to simulate the three-dimensional propagation of electrical waves in cardiac tissue and (2) a photon transport model to solve the forward optical problem and generate surface optical signals produced by voltage-sensitive dyes.

Briefly, let V_m(r⃗, t) be the cellular transmembrane potential produced by an electrical wave in cardiac tissue at location r⃗=(x, y, z) and time t; let β be the quantum yield of the voltage-sensitive dye; let Φ_e(r⃗) be the photon density (or fluence) of the excitation light in the tissue, and let Γ(z′,ρ) be the point-spread function describing the flux Γ of photons emitted by the fluorescent dye at point (x′,y′,z′) and exiting the tissue at point (x,y), with ρ the radial coordinate on the imaged xy surface, centered around (x’,y′), i.e. $\rho =\sqrt{{(x-x^{\prime })}^2+{\left(y-y^{\prime }\right)}^2}$. The fluorescence optical signal VF recorded from the surface point (x,y) can then be obtained by convolving V_m (r⃗, t) with Φ_e and Γ:

where V is the volume of the cardiac tissue.

We will investigate the sensitivity of Eq. (1) to the photon transport model by calculating Φ_e and Γ using both a Monte Carlo approach and the photon diffusion equation in slabs of ventricular tissue. We also investigate, for each transport model, the sensitivity to the slab thickness L.

2.2 Electrophysiological model

Electrical wave propagation in cardiac tissue was simulated using a well-established and validated guinea pig electrophysiological model [22]. This was achieved by using reaction-diffusion equations in order to obtain the spatial and temporal distribution of the transmembrane potential V_m(r⃗, t) (see for example [23]):

where C_m is the membrane capacitance, and D_E is the electrical diffusivity tensor. I_ion represents the total transmembrane ionic current and is obtained through the solution of several additional differential equations describing ion channel kinetics and intracellular ion concentrations. Here, we used the dynamic Luo-Rudy II model specifically created for guinea pig ventricular myocytes [22]. Propagation of electrical waves is much faster along than across muscle fibers. Hence, the diffusivity tensor D_E was scaled to produce steady-state conduction velocities of 60 cm/s in longitudinal and 20 cm/s in transverse direction, consistent with our experimental measurements. The fibers were assumed to rotate at a linear rate with depth. The rate of transmural fiber rotation was set to 12°/mm in counterclockwise direction when going from epi- to endocardium yielding a total fiber rotation of 120° in a 10 mm thick slab [1]. Note that we showed recently that the rate of transmural fiber rotation only mildly affects the epi-fluorescence activation patterns following epicardial stimulation [18]. Electrical activity following a twice threshold point stimulus on the epicardium was simulated.

Propagation of electrical waves was simulated in rectangular slabs of 20×20 mm and varying thicknesses L=2.5, 5.0, 7.5, and 10.0 mm. To integrate Eq. (2) we have used an explicit finite difference scheme as described elsewhere using a time step of 0.01 ms and a space step of 0.1 mm. Simulations were carried out on a parallel cluster consisting of 16 dual AMD Athlon MP2200+ processors running at 1.8 GHz. We used the MPI library and a “domain slicing” algorithm to parallelize the code [24]. A typical simulation took approximately 12 hours central processing unit time for the 10mm thick slab.

2.3 Monte Carlo simulations

An established Monte Carlo model of light transport as developed by Wang and co-workers [25, 26] was used to calculate the excitation fluence and emission point-spread functions. The input parameters required by the main Monte Carlo simulation program, MCML, are the photon absorption coefficient, µa, the photon scattering coefficient, µ_s, anisotropy factor, g, the refractive index of both the tissue, n_tissue, and surrounding medium n_surround (saline solution, n_tissue/n_surround=1.05), and the thickness of the tissue, L. The values used for µ _a, µ _s and g, were those provided by Ding and co-workers [13] for the 488 nm excitation and 669 nm emission of di-4-ANEPPS in heart tissue.

Table 1. Optical parameters for excitation (488 nm) and emission (669 nm).

View Table | View all tables in this article

To obtain the excitation fluence Φ_e as a result of uniform excitation (illumination) of the heart surface, we first obtained the tissue response to an infinitely narrow photon beam normally incident on the heart surface (MCML), assuming tissue optical properties at the di-4-ANEPPS excitation wavelength of 488 nm (see Table 1). The resulting tissue photon scattering and absorption response was then subject to a convolution procedure (using a second program, CONV) to produce the desired tissue response to a uniform photon beam with finite radius and normal incidence. The radius of the uniform, circularly flat, finite beam was set to 1.5 cm, large enough to obtain a fluence (in joules/cm²) versus tissue depth (in cm) at the origin (at r=0) similar to that of a uniform photon beam with infinite radius.

Next, we simulated the voltage-sensitive fluorescent photons that arise from inside the heart tissue and are eventually emitted from the upper heart surface, defined in MCML and CONV as the diffuse reflectance, R _d. To accomplish this, point photon sources were placed at regularly spaced depths, z, inside the tissue. Photon packets were launched isotropically from each of these point source locations. Isotropy was achieved by setting the anisotropy factor, g, to zero for the first propagation step only. For all subsequent photon packet propagation steps in the tissue, we used the value of the anisotropy factor, g, characteristic of the heart tissue at the 669 nm emission wavelength (Table 1). The other optical characteristics, µ _a, and µ _s, were as described in Table 1 for the di-4-ANEPPS 669 nm emission wavelength [13]. We then recorded the resulting emission, i.e., R _d (in photons/cm²) from the tissue surface. The line passing through a point source that is normal to the heart surface is hereafter referred to as the ‘optical axis.’ The value R _d at a distance ρ on the surface from the optical axis and for a point source at depth Z is nothing else than the value of the point spread function Γ(Z,ρ) in Eq. (1). The simulations were performed at spatial resolutions of 0.01 mm in the radial direction and 0.1 mm in the Z direction. The reason for the higher radial resolution is that the point-spread functions for sources near the surface are very narrow and require finer resolution to fully characterize their radial distribution.

All computations were performed on a desktop with an Athlon XP 3200+ processor. It took 63 minutes CPU to compute the point-spread function of a source at 0.7 mm depth in a 5 mm thick slab and using 20,000,000 photon packets.

2.4 Diffusion approximation

Optical diffusion theory has been used in several recent studies to calculate cardiac optical signals [16–18, 20–21, 23, 25]. Using the diffusion approximation one can derive analytical expressions for Φ_e and Γ in the slab geometry and for given boundary conditions. Under uniform illumination of the z=0 surface (further referred to as the epicardial surface) and using Robin or partial-flux boundary conditions, the excitation fluence Φ_e only depends on depth z [8]:

where Φ₀ is the density of incident light, δ_e is the attenuation length of excitation light, and d_e is the extrapolation distance for excitation.

The point-spread function for emission from a point source located at (x’,y’,z’) is given by

where δ is the attenuation length and d the extrapolation distance for emission [8]. The attenuation lengths are determined as $\sqrt{\frac{D}{{\mu }_{\alpha }}}$ with the photon diffusion coefficient D=(3µ_a+3µ_s’)^-1 (see Table 1 with µ_s’=(1 - g) µ_s). The tissue-saline interface was modeled following an approach developed by Aronson to calculate extrapolation distances (d_e=0.65 mm and d=1.27 mm). Note that we have also investigated the case of Dirichlet or zero boundary conditions that have been used previously in various computational studies of cardiac optical imaging. In this case the extrapolation distances were all set to 0.

2.5 Experiments

All experimental protocols conformed to the Guide for the Care and Use of Laboratory Animals. Guinea pigs (n=3) were injected with heparin (550 U/100g), and euthanized by sodium pentobarbital (7.5 ml/kg), after which the heart was immediately isolated and placed in ice-cold cardioplegic solution (CPS), composed of 13.44 KCl, 12.6 NaHCO3, 280 Glucose, and 34 Mannitol (in mmol/l). The aorta was cannulated and perfusion of oxygenated Tyrode’s solution at 37°C (composed of (in mmol/l) 130 NaCl, 24 NaHCO3, 1.2 NaH2PO4, 1.0 MgCl2, 5.6 Glucose, 4.0 KCl, 1.8 CaCl2; buffered to a pH of 7.4) was achieved using a constant pressure (80 mm Hg) Langendorff apparatus. The heart was put into a special transparent chamber and superfused with the same Tyrode solution at a rate of 30–40 ml/min. Perfusion and superfusion temperatures were continuously monitored and kept at 37±0.5°C by using two sets of glass heating coils and heated-refrigerated circulators. Electrodes were sutured to whole hearts to monitor the electric activity. All preparations were continuously paced in the mid-portion of the left ventricular free wall through a tungsten electrode, at a basic cycle length of 300 ms and at twice the excitation threshold. After a 30 min stabilization period, the excitation-contraction uncoupler diacetylmonoxime (DAM) was added to the Tyrode’s solution (15 mmol/l) to stop contractions of cardiac tissue. Finally the preparation was stained through a bolus injection into the perfusion flow of the voltage-sensitive dye DI-4-ANEPPS (from Molecular Probes, concentration 5 µg/ml).

The optical setup consisted of a 14-bit CCD video camera (“Little Joe”, SciMeasure, Decatur, GA) with a Computar H1212FI lens (focal length 12 mm, 1:1.2 aperture ratio; CBC Corp., Commack, NY). We used a collimated beam from a 250 W tungsten halogen lamp to uniformly illuminate the left ventricular free wall. The light was heat filtered and then passed through a 520±40 nm bandpass excitation filter. Fluorescence was recorded at 640±50 nm. The optical images (80x80 pixels) were acquired from a 20 mm×20 mm area of the preparation at 2000 frames per second. The background fluorescence was subtracted from each frame to obtain the voltage-dependent signal.

In all experiments, 50 optical action potentials were ensemble averaged to reduce noise and avoid the excessive use of spatial and temporal filters that might affect optical upstroke morphology. The alignment error of the optical action potentials was as little as one-half frame (0.25 ms).

2.6 Data analysis

To characterize the optical upstroke morphology, we used the relative fractional level V_F*, at which the time derivative of the voltage-sensitive fluorescence is maximal (see Fig. 1). For example, V_F*=0 means that the slope of the upstroke is maximal right at its onset, while V_F*=1 means that the slope is maximal at the end of the upstroke. To obtain accurate estimates of V_F*, we used cubic spline interpolation (5 additional points per interval) between actual data points in our fluorescence signals V_F from each pixel [16–18].

Previous studies have indicated that V_F* correlates to the optically-weighted subsurface wave propagation direction. The optically-weighted mean angle of propagation with respect to the epicardial surface, ϕ_F, was determined as in previous studies [16,17]: first, we computed the spatial gradient in transmembrane voltage, ∇V_m, over the entire tissue volume at each time step in the computer simulation. The location of the wave front in each frame was determined by identifying the voxels for which ∇V_m had reached its maximum value, ∇V_m (r⃗,t)_max, in time. Next, the local angle of the wave front with respect to the surface, ϕ_E(r⃗), was determined by computing the angle of the vector ∇V_m (r⃗,t)_max with respect to the surface. For a point (x,y) on the epicardial surface, we then computed the optically-weighted mean angle ϕF(x,y) of the propagating wave front beneath that point by convolution of the local wave front angle ϕ_E(r⃗) with the optical weighting function, Φ_e ·Γ, over the entire tissue volume; which was achieved by replacing V_m with ϕ_E in Eq. (1):

3. Results

3.1 Excitation fluence and intramural optical weight

 

Fig. 3. Excitation fluence and intramural optical weight. Panels A and B show the intramural excitation fluence Φ_e for two different tissue thicknesses (2.5 and 5mm respectively). The data was normalized to the photon density on the epicardial surface (z=0). Panel C and D show the relative contribution of an intramural layer at depth z to the total fluorescence for a uniform source distribution for two different tissue thicknesses (2.5 and 5mm respectively).

Download Full Size | PDF

We investigated the differences in excitation fluence and emission point-spread functions for the various models and different tissue thicknesses L. The results are presented in Fig. 3. Panels A and B show the intramural excitation fluence Φ_e for two different tissue thicknesses (2.5 and 5 mm respectively). We see that changing the type of boundary condition in diffusion theory (Zero vs. Robin) does not substantially affect the excitation fluence. An exponential fit yields a penetration depth (i.e. the exponential decay constant) of 0.58 mm for both models. Using the MC model, we found a two-fold increase in depth penetration of the excitation light and a typical subsurface maximum at 0.15 mm below the epicardial surface. Comparison between the panels A and B indicates minor dependence of the excitation fluence on slab thickness. Similar results were obtained in the 7.5 and 10.0 mm thick slabs (not shown).

Panel C and D show the relative contribution of an intramural layer at depth z to the total fluorescence for a uniform source distribution. It is obtained by integrating the point-spread function in the radial direction for various depths and calculating the relative contribution made to the total fluorescence (i.e. integration over the whole slab volume). Interestingly, we observe little differences between the MC model and the diffusion theory using Robin boundary conditions. However, diffusion theory with zero BC shows substantial deviations as its total emission has a faster drop as function of depth. This is a result of enforcing the zero BC at z=2.5 mm. This boundary effect is most pronounced in thin slabs, where the tissue thickness is comparable to the attenuation length of fluorescence. The differences between the models are greatly attenuated in the thicker preparations as illustrated in Panel D for the 5 mm thick slab, as all curves approach zero for depths larger than 3 mm.

 

Fig. 4. Subsurface optical integration volume for the different models (see text for details).

Download Full Size | PDF

The optical properties of tissue for both excitation and emission determine the subsurface volume that contributes to the optical signal in a single pixel, and thus the degree of optical “blurring”. We investigated the size and shape of this optical integration volume (OIV) for the various photon transport models by calculating the relative contribution of each voxel to the total signal recorded from the pixel (0,0) on the epicardium. To do so, we considered a uniform transmembrane potential distribution (V_m=1) in Eq. (1) and calculated Φ_e·Γ/V_F(0,0) for each voxel. Because of the cylindrical symmetry with respect to the optical axis, and for ease of representation, we further integrated all contributions within annular regions of radius r_k and at constant depths z, as shown on the left panel of Fig. 4. The results are shown on the right for the 2.5 mm thick slab: dark red indicates a small contribution to the total optical signal (<0.1%), whereas yellow stands for a large contribution (~2%). The white lines show isosurfaces obtained at 0.1% intervals. The green isosurface encloses a volume that contributes 75% of the total optical signal in the pixel of interest, and will be our definition of the OIV. The diffusion model with zero BC yields a substantially smaller OIV than both the diffusion model with Robin BC and the MC model. The difference between the OIVs of the two diffusion models is more pronounced in the lateral direction (r_k) than in depth (z): both models yield contributions up to 1mm depth, but the zero BC OIV is only 2 mm at its widest vs. 5.2 mm for Robin BC. When comparing the OIVs of the diffusion model with Robin BC and the MC model, we find that they have similar lateral extents (5.2 vs. 4.8 mm respectively), but that the MC model yields contributions up to 1.5 mm depth, consistent with the deeper penetration of excitation fluence as shown in Fig. 3. Similar results were found for the thicker slabs, with slight increases (<5%) in the size of OIV for the diffusion model with zero BC and the MC model and no noticeable effect in the case of diffusion with Robin BC (not shown).

3.2 Optical action potential upstroke

Having calculated the excitation fluence and point-spread function, and characterized the OIV for the various models, we proceed with calculating voltage-sensitive optical signals for the various photon transport models utilizing the hybrid electro-optical equations as described in the methods section. We focus on the optical action potential upstroke, as previous studies have indicated that this is the part of the optical action potential that is most sensitive to optical blurring effects.

Activation maps are a widely used tool to investigate spatiotemporal dynamics of cardiac activity [7]. They are directly related to the optical upstrokes, as they are obtained by measuring in each pixel the time at which the normalized V_F reaches a preset level during its upstroke. Figure 5 shows optical activation maps obtained for epicardial point stimulation (white cross) in the 5 mm thick slab. For comparison we show an activation map obtained experimentally for mid free wall pacing of the left ventricle of an isolated guinea pig heart. The three left panels show qualitative agreement between the different models, with fast propagation (short activation times indicated by purple-blue colors) along the epicardial fiber orientation (vertical direction), and slow propagation in a perpendicular direction. A similar activation pattern is observed experimentally. Careful analysis reveals, however, small differences between the different models: the earliest activation occurs at 4.73 ms, 5.22 ms, and 9.9 ms after stimulation, for diffusion with zero BC, diffusion with Robin BC, and MC respectively. The latest activation occurs at 50.32 ms, 50.19 ms, and 47.75 ms, respectively. For comparison, the simulated epicardial electrical activity V_m yields earliest activation at 0.8 ms and latest activation at 51.49 ms after applying the stimulus. The differences between optical and electrical activation times reflect the increased duration of the optical upstroke due to optical smearing. This is further illustrated in the right panel, where optical upstrokes recorded from the black asterisk are compared for the different models. Comparison with an experimental optical upstroke recorded from a similar location shows that the MC model yields the most realistic optical upstrokes in this particular pixel.

 

Fig. 5. Optical activation maps and optical upstrokes. The three left panels show simulated activation maps obtained using the different photon transport models. An experimental activation map is shown on the right, as well as a comparison of optical upstrokes obtained from the pixel indicated by an asterisk * in the activation maps.

Download Full Size | PDF

We further quantified the differences in optical upstroke durations by measuring them in each image pixel for the different photon transport models and different thicknesses of the slabs of cardiac tissue. The results are presented in Table 2 as mean ± standard deviation (SD) over the whole field of view. As expected we find that diffusion with zero BC yields the shortest optical upstrokes and MC the longest, independent of tissue thickness. The diffusion model with zero BC and the MC model are the most sensitive to the tissue thickness, as the upstroke duration shows a 10% increase when increasing the slab thickness from 2.5 to 5 mm thick in both models. The diffusion model with Robin BC shows in this case only a 5% increase in upstroke duration. Further increasing the slab thickness (up to 10 mm) has in all models little or no effect on upstroke duration. For comparison, the electrical upstroke duration was 2.88±0.66 ms across the whole volume of the slab. Consistent with previous studies [11,16–18], the optical upstroke duration is longer than electrical upstroke duration due to optical blurring. The MC model yields the longest upstrokes due to its larger OIV. In our experiments on guinea pig left ventricles we obtained upstroke durations of 12.86±0.65, 11.13±1.75, and 11.36±1.96 ms in the three different hearts.

Table 2. Optical upstroke durations (mean ± SD ms) for various tissue thicknesses L.

View Table | View all tables in this article

Table 2 presents a qualitative comparison between the different models as upstroke durations have been averaged over the whole field of view. However, the subsurface wave front orientation and propagation velocity can change dramatically as the wave moves away from the pacing site [3]. We should therefore expect the differences between the three models to vary with location on the imaged surface due to their different OIVs. Figure 6 shows color maps of the differences in upstroke duration between the diffusion and MC models. Black indicates maximal difference, whereas yellow-white stands for minimal difference. Both maps are symmetrical with respect to the stimulus location. They indicate that the largest differences in upstroke duration are found at about 5 mm from the stimulus site, essentially in the direction of slow propagation (compare with activation maps in Fig. 5). The smallest differences are found at about the same distance from the stimulus location, but in the direction of fast propagation. The differences in upstroke duration can further be appreciated from representative optical upstrokes in both areas (a and b, respectively) as shown on the right.

 

Fig. 6. Differences in optical upstroke durations in the various models. The left panels show color-coded maps of the differences in optical upstroke durations in each pixel between the diffusion models and the MC model. The tracings on the right compare the optical upstrokes of the three models recorded from two different areas (a and b).

Download Full Size | PDF

3.3 VF* vs. subsurface wave front orientation

Our simulation results have indicated that Monte Carlo models yield more realistic optical upstroke durations when compared to the diffusion based model. In this subsection we investigate the correlation between V_F* and the optically averaged subsurface wave front orientation for the different models and tissue thicknesses, and show how it can be applied to reconstruct subsurface wave front orientation from experimental recordings.

Table 3. Maximal and minimal values of VF* for the different models and tissue thicknesses.

View Table | View all tables in this article

Figure 7(A) shows V_F* maps obtained after epicardial point stimulation in the 5mm thick slab. For all models we observe the typical pattern that was recently reported in experimental studies in both porcine slabs of ventricular tissue and guinea pig hearts: a vertical strip centered around the stimulation site of low V_F* (<0.3) is observed, corresponding to regions of fast epicardial propagation (compare with Fig. 5). In the direction of slow propagation (horizontal direction), V_F* rapidly increases to values above 0.7. The maps obtained with the diffusion models are almost undistinguishable from each other. The MC model results in a slightly wider strip of low V_F*, and a broader range of values over the same field of view. This is further quantified in Table 3 where the maximal and minimal values for V_F* are shown for the different models and different tissue thicknesses. For comparison, we show a V_m* map, which represents the relative level of maximal electrical upstroke velocity obtained from the simulated electrical activity on the epicardial surface. The V_m* values all lie close to 0.45, clearly indicating that the typical V_F* patterns are a result of subsurface contribution. An experimental V_F* map obtained for mid free wall pacing in the isolated guinea pig heart is shown in Fig. 7(C).

 

Fig. 7. V_F* and optically weighted subsurface angle ϕ_F. Panel A shows V_F* maps obtained after epicardial point stimulation in the 5mm thick slab. Panel B shows the optically averaged subsurface wave front orientation ϕ_F with respect to the epicardial surface. The electrical V_m* and ϕ_E maps are shown for comparison. Panel C shows an experimental V_F* map and Panel D shows simulated 3D isochronal surfaces following an epicardial point stimulus.

Download Full Size | PDF

Figure 7(B) shows the optically averaged subsurface wave front orientation ϕ_F with respect to the epicardial surface (see methods section). An angle of 90 denotes a plane wave propagating away from the epicardium, whereas an angle of -90 denotes a plane wave propagating towards the epicardial surface. The three models yield similar results, and agree qualitatively with the electrical wave front orientation ϕ_E at the epicardial surface. The optical maps are however blurred due to subsurface contributions. In all case, we observe positive angles in the direction of fast propagation (vertical direction), and a gradual change towards negative angles in the other directions. This phenomenon can be explained by looking at the actual 3D wave propagation pattern in the slab as shown in Panel D. This figure is a 3D activation pattern showing isochronal surfaces at 10 ms intervals. The wave originates in the middle of the epicardial (top) surface after applying the stimulus. The wave expands rapidly along the epicardial fibers (x-axis), but more slowly in the perpendicular direction (y-axis), resulting in an ellipsoidal shape. The transmural fiber rotation, however, affects wave front orientation such that, in the direction of slow propagation, the wave front gradually rotates and is eventually oriented towards the epicardial surface. In the direction of fast propagation the wave front remains oriented away from the epicardial surface. This readily explains the typical ϕ_F patterns presented in Panel B.

 

Fig. 8. Linear regression analysis of V_F* vs ϕ_F for the different photon transport models. Each dot in the scatter plot indicates the results obtained for a single pixel and the red line shows the linear fit to the data.

Download Full Size | PDF

Table 4. Linear regression analysis of subsurface angle ϕ_F vs V_F* in a slab of 5 mm thickness. ϕ_F=A+B·V_F* and R is the correlation coefficient

View Table | View all tables in this article

By comparing Figs. 7(A) and 7(B) a correlation between V_F* and ϕ_F can already be observed independent of the model: high values of V_F* indicate waves approaching the epicardial surface, while the opposite is true for low values of V_F*. Figure 8 shows scatter plots of V_F* vs. ϕ_F for the data presented in Fig. 7. The red line is the result of a linear regression analysis which results are shown in Table 4. The absolute value of the correlation coefficient R is, for all models, close to 1 and indicates the strong linear relationship between V_F* and ϕ_F, especially for the MC model. The main difference between models lies in the slope B, with larger values obtained for the diffusion based models. Interestingly, the linear relationship between V_F* and ϕ_F did not show significant dependence on tissue thickness for any model (<2% variations in parameters A, B and R).

3.4 Reconstructing subsurface wave front orientation from experimental epi-fluorescence recordings

The linear relationship between V_F* and ϕ_F allows for a straightforward reconstruction of the average subsurface wave front orientation from experimental data. Figure 9 shows the results for such reconstruction in the case of two different activation patterns: mid free wall pacing (top) and sinus rhythm (bottom). The left panels show the experimentally obtained V_F* maps, whereas the right panels show the reconstructed ϕ_F for the different models. For each type of activity, the reconstructed patterns look qualitatively similar, because they are all a linear transformation of the same V_F* map. The reconstructed maps using the MC model are less sharp than those obtained with the diffusion based models. This is primarily due to the smaller slope B for this model, again reflecting the larger OIV in this case and more pronounced optical blurring. For paced activity (top) the major differences are found in the direction of slow propagation with maximal differences in ϕ_F of about 20°. The MC model yields in general somewhat larger values of ϕ_F, especially in the basal region of the ventricles.

 

Fig. 9. Reconstructing subsurface wave front orientation from V_F* maps for epicardial pacing (top) and sinus rhythm (bottom). For each row, the V_F* map is shown on the left and the reconstructions using the different models are shown on the right.

Download Full Size | PDF

4. Discussion

Cardiac optical imaging has become the method of choice to visualize propagating electrical waves in the heart during both normal and abnormal activity [6, 7]. The effects of photon scattering and absorption on the optical signals have been acknowledged for over a decade [11–15], but it is only recently that hybrid electro-optical models have been utilized to quantify these effects and carefully reproduce cardiac optical imaging experiments in silico [8–10,16–21]. Even today, most of these models are based on the photon diffusion equation, either with zero BC or Robin BC [8–9, 17–18, 20–21]. One of the main advantages of the diffusion equation is its versatility as it can be solved on any type of geometry using finite-difference or finite-element schemes [18, 21, 23], or even yield analytical solutions when considering, for example, the slab geometry [8, 20]. Diffusion based models have successfully been applied to study cardiac optical signals in thick slabs of cardiac tissue, and recently led to the development of a novel technique to reconstruct subsurface wave front orientation utilizing V_F* maps [16, 17]. We recently showed that experimental V_F* maps obtained in thick porcine ventricles [17] are similar to those observed in small guinea pig hearts [15]. However, the photon diffusion equation is believed to fail in such thin hearts. It therefore remained unclear how previous computational tools developed to calculate the subsurface wave angle and obtained utilizing the diffusion equation can be applied to smaller tissue samples.

The goal of this paper was to address this question by using a hybrid electro-optical model based on a MC photon transport model and compare the simulated cardiac optical signals with those obtained using diffusion models. Moreover, we investigated the sensitivity of the optical signal to tissue thickness, to cover the broad range of cardiac optical mapping experiments being routinely performed [6,7].

We have found that MC based models yield the longest optical upstroke durations, and are in good agreement with experimentally recorded values. Previous studies using the diffusion equation indicated moderately prolonged optical upstrokes using different models of electrical activity, consistent with our findings [16, 17, 20]. The ionic model used here has been validated against experimental data obtained specifically from guinea pig ventricular cells [22]. Our results therefore indicate that MC models are the most realistic in terms of optical upstroke duration.

The differences in upstroke durations between the different models depend however on the location where these are recorded with respect to the stimulus and thus on subsurface wave orientation (Fig. 6). In the region of maximal difference (point a in Fig. 6), the wave front propagates at relatively low velocity towards the epicardial surface (Fig. 7(D)). Hence, the optical upstrokes will be relatively long and differences in depth extent of the OIV between different models will be emphasized. This is the reason why MC model yields substantially longer upstroke in that region, while both diffusion based models yield similar durations. In the region of minimal difference (point b in Fig. 6) however, wave propagation is much faster and away from the epicardial surface (Fig. 7(D)), yielding smaller differences between the different models. Our simulation study indicates also that the effects of tissue thickness on the optical signal are relatively minor. Its most distinct effect was a modest prolongation of the upstroke (~10%) when increasing the thickness from 2.5mm to 5mm.

Our most important finding is that V_F* shows a strong linear correlation with the optically averaged wave front orientation, independent of the model and tissue thickness. As such V_F* is a good indicator of sub-surface wave propagation, as is demonstrated for two distinct types of activity, i.e. epicardial pacing and sinus rhythm. During sinus rhythm, the ventricles are activated through a specialized conduction network of Purkinje fibers. These fibers emerge in the ventricular muscle primarily on its endocardial surface. Hence, we expect the wave front to propagate primarily towards the epicardial surface (negative angle ϕF) resulting in high VF* (>0.4) [15]. Our experimental V_F* map and the reconstructed ϕ_F maps shown in Fig. 9 are in good agreement with this general activation picture. The parameters obtained from the linear correlation between V_F* and ϕ_F, and necessary for such reconstruction, do however depend on the photon transport model (Fig. 8 and Table 4). This is to be expected since the models yield quite different optical integration volumes. Previous computational studies utilizing the diffusion equation should therefore be revisited accordingly [13, 15, 18], but without affecting their main conclusions. Our study highlights that utilizing V_F* maps to infer subsurface wave front orientation from cardiac epi-fluorescence signals is a robust method, which can be applied in a wide range of cardiac optical mapping experiments across species.

It is worth noting that the formulation of the photon diffusion equation utilized in this study and previous work [17–18,20–21], assumes that all incoming photons immediately diffuse into the tissue. This leads to monotonically decreasing excitation fluences as a function of depth as shown by Eq. (2). A different approach consists in treating the first few scattering events of the excitation photons separately through the diffusion equation’s source function [27]. In such case, a subsurface maximum is observed in the excitation fluence, similar to what is obtained with the Monte Carlo model. This approach was recently used by Roth [28] to investigate excitation fluences in a cardiac wedge model, but has not yet been applied in combination with cardiac electrophysiological models. Such hybrid models could potentially resolve some of the discrepancies observed in this paper between the diffusion and Monte Carlo based models of optical upstroke durations, and would be a worthwhile addition to the hybrid electro-optical models.