Reconstructing subsurface electrical wave orientation from cardiac epi-fluorescence recordings : Monte Carlo versus diffusion approximation

The development of voltage-sensitive dyes has revolutionized cardiac electrophysiology and made optical imaging of cardiac electrical activity possible. Photon diffusion models coupled to electrical excitation models have been successful in qualitatively predicting the shape of the optical action potential and its dependence on subsurface electrical wave orientation. However, the accuracy of the diffusion equation in the visible range, especially for thin tissue preparations, remains unclear. Here, we compare diffusion and Monte Carlo (MC) based models and we investigate the role of tissue thickness. All computational results are compared to experimental data obtained from intact guinea pig hearts. We show that the subsurface volume contributing to the epi-fluorescence signal extends deeper in the tissue when using MC models, resulting in longer optical upstroke durations which are in better agreement with experiments. The optical upstroke morphology, however, strongly correlates to the subsurface propagation direction independent of the model and is consistent with our experimental observations. ©2008 Optical Society of America OCIS codes: (170.3880) Medical and Biological Imaging; (170.3660) Light propagation in tissues References and links 1. D. Streeter, Handbook of Physiology, (Bethesda, MD, American Physiological Society, 1979). 2. R. A. Gray, A. M. Pertsov and J. Jalife, "Spatial and temporal organization during cardiac fibrillation," Nature 392, 75-78 (1998). 3. O. Bernus, M. Wellner and A. M. Pertsov, "Intramural wave propagation in cardiac tissue: asymptotic solutions and cusp waves," Phys. Rev. E 70, 061913 (2004). 4. F. Fenton and A. Karma, “Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation,” Chaos 8, 20-47 (1998). 5. K.H.J.W. Ten Tusscher, R. Hren and A.V. Panfilov, “Organization of ventricular fibrillation in the human heart,” Circ. Res. 100, e87-e101 (2007). 6. D. S. Rosenbaum and J. Jalife, Optical mapping of Cardiac excitation and arrhythmias, (Armonk, N Y, Futura Publishing Company, Inc. 2001). 7. I. R. Efimov, V. P. Nikolski and G. Salama, "Optical imaging of the heart," Circ. Res. 95, 21-33 (2004). 8. V. D. Khait, O. Bernus, S. Mironov and A. M. Pertsov, “Method for the three-dimensional localization of intramyocardial excitation centers using optical imaging,” J. Biomed. Opt. 11, 34007 (2006). 9. M. Wellner, O. Bernus, S. F. Mironov and A. M. Pertsov, “Multiplicative optical tomography of cardiac electrical activity,” Phys. Med. Biol. 51, 4429-46 (2006). 10. E. M. C. Hillman, O. Bernus, E. Pease, M. B. Bouchard and A. M. Pertsov, “Depth-resolved optical imaging of transmural electrical propagation in perfused heart,” Opt. Express 15, 17827-17841 (2007). 11. S. D. Girouard, K. R. Laurita, and D. S. Rosenbaum, "Unique properties of cardiac action potentials recorded with voltage-sensitive dyes," J. Cardiovasc. Electrophysiol. 7, 1024–1038 (1996). 12. W. Baxter, S. F. Mironov, A. V. Zaitsev, A. M. Pertsov and J. Jalife, "Visualizing excitation waves in cardiac muscle using transillumination," Biophys. J. 80, 516–530 (2001). #97989 $15.00 USD Received 26 Jun 2008; revised 11 Aug 2008; accepted 18 Aug 2008; published 21 Aug 2008 (C) 2008 OSA 1 September 2008 / Vol. 16, No. 18 / OPTICS EXPRESS 13758 13. L. Ding, R. Splinter and S. B. Knisley, "Quantifying spatial localization of optical mapping using Monte Carlo simulations," IEEE Trans. Biomed. Eng. 48, 1098-1107 (2001). 14. D. L. Janks and B. J. Roth, "Averaging over depth during optical mapping of unipolar stimulation," IEEE Trans. Biomed. Eng. 49, 1051-1054 (2002). 15. M. A. Bray and J. P. Wikswo, “Examination of optical depth effects on fluorescence imaging of cardiac propagation,” Biophys. J. 85, 4134-4145 (2003). 16. C. J. Hyatt, S. F. Mironov, M. Wellner, O. Berenfeld, A. K. Popp, D. A. Weitz, J. Jalife and A. M. Pertsov, "Synthesis of voltage-sensitive fluorescence signals from three-dimensional myocardial activation patterns," Biophys. J. 85, 2673-2683 (2003). 17. C. J. Hyatt, S. F. Mironov, F. J. Vetter, C. W. Zemlin and A. M. Pertsov, "Optical action potential upstroke morphology reveals near-surface transmural propagation direction," Circ. Res. 97, 277-284 (2005). 18. C.W. Zemlin, O. Bernus, A. Matiukas, C.J. Hyatt and A.M. Pertsov, “Extracting Intramural Wavefront Orientation From Optical Upstroke Shapes in Whole Hearts,” Biophys. J. 95, 942-950 (2008). 19. M.J. Bishop, G. Bub, A. Garny, D.J. Gavaghan and B. Rodriguez , “An investigation into the role of the optical detection set-up in the recording of cardiac optical mapping signals: A Monte Carlo simulation study,” Physica D (to be published). 20. O. Bernus, M. Wellner, S. F. Mironov and A. M. Pertsov, "Simulation of voltage-sensitive optical signals in three-dimensional slabs of cardiac tissue: application to transillumination and coaxial imaging methods," Phys. Med. Biol. 50, 215-229 (2005). 21. M. J. Bishop, B. Rodriguez, J. Eason, J. P. Whiteley, N. Trayanova, and D. J. Gavaghan, “Synthesis of voltage-sensitive optical signals: application to panoramic optical mapping,” Biophys. J. 90, 2938–2945 (2006). 22. G. M. Faber and Y. Rudy, “Action potential and contractility changes in Nai overloaded cardiac myocytes: a simulation study,” Biophys. J. 78, 2392–2404 (2000). 23. O. Bernus, K. S. Mukund and A. M. Pertsov, "Detection of intramyocardial scroll waves using absorptive transillumination imaging," J. Biomed. Opt. 12, 14035 (2007). 24. R. Zaritsky and A. M. Pertsov, “Simulation of 2-D spiral wave interactions on a Pentium-based cluster,” in Proc. of Neural, Parallel, and Scientific Computations, M. P. Bekakos, G. S. Ladde, N. G. Medhin, and M. Sambandham, eds., (Dynamic Publisher, Atlanta, 2002). 25. L.-H. Wang, S. L. Jacques, and L.-Q. Zheng, "MCML – Monte Carlo modeling of photon transport in multilayered tissues," Comput. Methods Programs Biomed. 47, 131-146 (1995). 26. L.-H. Wang, S. L. Jacques, and L.-Q. Zheng, "CONV – Convolution for responses to a finite diameter photon beam incident on multilayered tissues," Comput. Methods Programs Biomed. 54, 141-150 (1997). 27. S. T. Flock, M. S. Patterson, B. C. Wilson and D. R. Wyman, “Monte Carlo modeling of light propagation in highly scattering tissues – I: model predictions and comparison with diffusion theory,” IEEE Trans. Biomed. Eng. 36, 1162-1168 (1989). 28. B. J. Roth, “Photon density measured over a cut surface: implications for optical mapping of the heart,” IEEE Trans. Biomed. Eng. 55, 2102-2104 (2008).


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][4][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][9][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.
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][12][13][14][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).ofa 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).
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 threedimensional 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 wellestablished 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.

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][9][10][19][20][21].These hybrid models consist of: (1)  with Φ e and Γ: (1) 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.

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 reactiondiffusion equations in order to obtain the spatial and temporal distribution of the transmembrane potential ) , ( t r V m (see for example [23]): (2) 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 x 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.

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.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 2 ) 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 2 ) 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.

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 Φ 0 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 ( 4) where δ is the attenuation length and d the extrapolation distance for emission [8].The attenuation lengths are determined as a D μ / 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.

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 x 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).

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][17][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] 1): (5) 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.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).

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.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.66ms across the whole volume of the slab.Consistent with previous studies [11,[16][17][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.  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.

V F * 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.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).).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.  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).

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.

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][12][13][14][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][9][10][16][17][18][19][20][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 finitedifference 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 V F * (>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.

,
: first, we computed the spatial gradient in transmembrane voltage, m V ∇ , 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 in time.Next, the local angle of the wave front with respect to the surface, the surface.For a point (x,y) on the epicardial surface, we then computed the opticallyweighted mean angle φ F (x,y) of the propagating wave front beneath that point by convolution of the local wave front angle ) r ( E φ with the optical weighting function, Φ e •Γ, over the entire tissue volume; which was achieved by replacing V m with φ E in Eq. (

Fig. 3 .
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).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

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

Fig. 5 .
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.

Fig. 6 .
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).

Fig. 7 .
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.

Figure 7 (
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

Fig. 8 .
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.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 Diff Zero Diff Robin MC A 1.59 1.53 1.09 B -3.72 -3.57-2.22 R -0.94 -0.93 -0.98By 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 Table4.The absolute value of the correlation coefficient R is, for all models, close to 1 and indicates the strong linear relationship between

Fig. 9 .
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.

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

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