High-Throughput Analysis of Optical Mapping Data Using ElectroMap

1. Optical mapping data collection

1. Perform cardiac optical mapping using one of a wide range of experimental models including intact and isolated whole hearts^6,18, isolated atria^14,19, ventricular wedges²⁰, cardiac slices^21,22, and cellular monolayers²³. See associated references for experimental designs to collect raw optical mapping data from these preparations. Provided that data obtained can be converted to a tiff stack or saved in a .MAT file, it should be analyzable using ElectroMap. This includes data of varying dimensions (square/rectangular) and resolutions (maximum tested currently 2048 pixels x 2048 pixels).

2. Software installation and start-up

NOTE: Below are detailed the two methods for installing and running ElectroMap – either within MATLAB run from the source (.m) code or as a standalone executable file (.exe for windows). The final software and its functionality are invariant between the two setup options (other than a few differences in directory navigation). Therefore, the main considerations for choosing version to install are access to MATLAB and required toolboxes and whether access to source code is desired. Where possible, it is recommended to use the MATLAB version for faster start up times, shorter processing times, and easier error reporting.

1. Setup 1: Running electromap within MATLAB
   1. Install MATLAB. ElectroMap was designed in MATLAB 2017a, however, software has been tested for use in all subsequent releases of MATLAB (up to 2018b at time of writing). The following toolboxes are required: Image Processing, Signal Processing, Statistics and Machine Learning, and Curve Fitting.
   2. Download/clone all files from the latest ‘source code’ release of ElectroMap from the GitHub repository (https://github.com/CXO531/ElectroMap). Unzip the downloaded contents to a desired location.
   3. Open MATLAB and navigate to the folder location hosting the ElectroMap source code. Then, open the file ElectroMap.m and press run in the editor, or alternatively type ElectroMap in the command window and press RETURN. This will start the ElectroMap user interface, Figure 1A.
2. Setup 2: Standalone .exe file
   1. Download the installer file: https://drive.google.com/open?id=1nJyI07w9WIt5zWcit0aEyIbtg31tANxI.
   2. Follow the instructions in the installer, which will download required the MATLAB runtime from the web alongside ElectroMap software.
   3. Run ElectroMap.exe.
      NOTE: Start up time for the standalone version can be several minutes.

3. Image loading and pre-processing

1. Press Select Folder and navigate to the location of the data file(s) to be analyzed. This will populate the left-hand listbox with all files within that directory that are of the correct file type (.tif or .MAT). .MAT files must only contain the image stack variable.
   NOTE: Only folders and not individual files will appear as you navigate through the directory selector.
2. Select file to be loaded from within the interface and press Load Images.
   1. Once loaded, the first frame will appear, and the red outline will indicate automatic thresholding of the image. If needed, reload previously used ROIs by selecting Save/Load ROI. In this case, skip step 3.3.
   2. As default, thresholding is based on the pixel intensities in the first frame. If desired, modify this to a threshold based to the signal time course amplitude by changing the option in the Image for threshold dropdown menu. Please note that once the thresholding is selected, it is then applied for the whole image stack.
3. If desired, change the threshold option to manual, which will activate the slider to manually adjust the image threshold. Additionally, crop images (Crop Image) and/or draw a custom region of interest (Custom ROI) for analysis by selecting the appropriate tick box(es) below threshold options. Note that advanced options for region of interest selection such as number of areas are available from ROI Selection from the top menu.
4. Once an appropriate threshold has been applied, press Process Images to apply processing. Settings for processing are detailed below (step 3.4.1-3.4.5). At this point, ensure that the correct camera settings have been entered. These are Pixel Size in µm (IMPORTANT: this is the image pixel size, and not the size of the pixels that make up the chip or equivalent hardware in the imaging device) and Framerate in kHz.
   1. For signal inversion, tick the Invert Data checkbox to enable. If reported fluorescent signal is inversely proportional to parameter of interest (as with commonly used potentiometric dyes) the signal can be inverted.
   2. For spatial filtering, select Gaussian or Average from the kernel menu. The size of the spatially averaged area is controlled by the Size input adjacent to the Kernel dropdown menu (i.e. 3 results in 3 pixel x 3 pixel filter kernel). When applying a Gaussian filter, the standard deviation can also be set from the Sigma input.
   3. For baseline correction, select Top-hat²⁴ or polynomial (4^th or 11^th degree) correction²⁵ from the Baseline menu. Correction can be applied to each pixel individually (long processing time) or as an average of the entire image (quicker but assumes homogenous baseline alterations). Top-hat correction can also be modified by setting the Top-Hat Length in milliseconds, adjacent to the baseline selection dropdown menu. The length of the Top-Hat kernel should be greater than the timescale of the individual action potentials/calcium transients.
   4. For temporal filtering, select Savitzky-Goaly or infinite impulse (IIR) filtering from the Filtering menu.
      NOTE: Other than for the tissue averaged signal that appears in the bottom left, temporal filtering is applied to each pixel individually at time of parameter quantification from ensemble averaged image ranges. This has been implemented to reduce processing time by filtering small sections of data when required rather than entire files.
   5. For frame removal, note that if the Remove Frames option is selected, large peaks with amplitude greater than the signal of interest can be removed from the image set. This may be useful in optically paced datasets such as optogenetic pacing where depolarization is initiated by optical activation of opsins such as channelrhodopsin 2^11,12.
      NOTE: As frame removal will potentially introduce unphysiological step changes into the image signals, temporal filtering may introduce artifacts to the data and so is not recommended here.
5. Note that signal will be segmented once Process Images has been selected according to the options under Segmentation options, however this can quickly be changed without reprocessing the entire dataset (see section 4).

4. Data segmentation and ensemble averaging

NOTE: Once the file has been processed, peaks in the tissue averaged signal (bottom right trace, Figure 1A) will have been detected and labelled by red circles. Only peaks above a set threshold (blue line on trace that is set by Peak Threshold) are counted. Additionally, peaks are only counted if they are sufficiently delayed compared to the previous peaks, set by the Min Peak Distance input. Signal is then segmented based on the detected peaks. First, the effective cycle length (CL) of each peak is calculated by measuring the time between it and the next peak. If a number of peaks (set by Min Number of Peaks input) have similar CLs (threshold for which is set by Minimum Boundary input) then they are grouped and the average CL for those peaks calculated.

1. For further segmentation of the data, press Segment Signal. Sub-segmentation options are: None – all peaks with same CL grouped together; Tümü – Segments of n_peaks within the constant CL times (n_peaks is set by the Segment size input) are identified; Last – Final n_peaks before a CL change are identified and grouped, and all others are not analyzed; and Single Beat – This is equivalent to applying the Tümü segmentation with n_peaks = 1, and so no grouping or ensemble averaging (see 4.5) are applied. This can be applied by selecting the Single Beat button.
   1. Apply custom segmentation of the signal by zooming in on a time of interest and selecting Segment Signal. This will add an additional option entitled Zoomed Section to the section list box, corresponding to the time points selected.
2. The results of the segmentation will appear in the list-box adjacent to the tissue averaged signal, and will show section number and the estimated CL. All segmented time sections are denoted by different colors. Select a segment from the list-box to highlight that section in red. This will also automatically trigger analyses of this section, as if the Produce Maps button was selected (see section 5).
3. Analyses of grouped peaks will be performed on the ‘ensemble averaged’ data. This involves averaging the peaks in a segment together, with the reference times being the peaks identified in step 4.2. Update the time window to average by modifying the before and after inputs and pressing Segment Signal.

5. Action potential/calcium transient duration and conduction velocity analysis

1. Once images have been processed, the Produce Maps button will become active. Press Produce Maps to apply action potential duration (APD), activation time, conduction velocity and SNR analysis. By default, the analysis will be applied to the first signal segment. Select other segments from the list-box will apply analysis to chosen segment.
   NOTE: Results of analysis are displayed in results table, including mean, standard deviation, standard error, variance and 5th to 95th percentile analysis. Duration maps are termed ‘APD’ maps however, calcium signals processed using the same settings will measure calcium transient duration.
2. Select Get Pixel Info to see a detailed display of the signal from any pixel within the image, and Compare Pixels to simultaneously plot signals from up to 6 locations.
   1. Use the Signal Processing panel to adjust settings for duration analysis. These are: Duration – Time of percentage repolarisation/decay to measure from peak; ‘APD’ Baseline – Time period of signal that is defined as reference baseline for amplitude measurements; and ‘APD’ Start time – Start time for duration measurements. These are the same options for deciding the activation time for isochronal maps (discussed below) and are termed: Start ( d²F/dt²­­_max^­), Upstroke (dF/dt­­_max^­), Depolarisation midpoint (time of 50% amplitude), Peak (time of maximum amplitude). These definitions applied to mouse and guinea pig action potentials are shown in Figure 2A.
      NOTE: Changing any of these options will automatically update the duration map and the results table. Map scale and outlier removal options are also available.
3. Conduction velocity is also measured automatically within the main software interface. This is achieved using the multi-vector method of Bayly et al²⁶ from the isochronal map defined by the chosen activation measure (discussed in step 5.4). Press Activation Points to render a 3D representation of the activation map.
4. The multi-vector conduction velocity measurement method spatially segments the isochronal map into regions of n x n pixels. Set the value of n using the Local Window Size input, and set the range of activation times to apply analysis to using the Fitting activation times inputs.
   NOTE: For each local region, a polynomial surface, f, is fitted that best describes the relationship between activation time and spatial position, (x,y). The gradient vector, CV_local, of this surface is then calculated as:
   (1)
   where denotes the two-dimensional cartesian spatial differential operator²⁶.
5. For each pixel in the isochronal map, a local vector representing speed and direction of conduction is calculated. Select Isochronal Map with vectors from the display dropdown menu to view this analysis.
6. SNR is calculated as the ratio of the maximum amplitude compared to the standard deviation of the signal at baseline. This analysis is performed post all processing steps. Press SNR calculation in the top menu to edit settings for the period of the signal defined as baseline.

6. Conduction analysis module

1. Press Conduction to access more detailed analysis of conduction velocity. This opens a separate module where conduction can be quantified using the Bayly multi-vector method as in the main interface, single vector methods, and as an activation curve.
2. Press Single Vector to analyze conduction using the single vector method, where CV is calculated from the delay in activation time between two points. This can be done using Automatic or Manual methods, selectable below the Single Vector button.
   1. For automatic single vector method, select a distance and start point from which to measure conduction. The software will then perform a 360-degree sweep from the selected point, measuring the time delay and calculating the associated conduction velocity along all directions in 1-degree increments. The results of this analysis are displayed in the graph adjacent to the map, and the direction of slowest conduction is shown in red.
   2. For manual single vector method, choose both a start and end point from the isochronal map to calculate conduction velocity. To select a new start point, press Clear Start Point.
3. Press Local Vector to apply the multi vector method, with the settings matching those from the main interface. Within the conduction module, the distribution of conduction speeds, as well as the angular distribution of calculated vectors and angular dependence of conduction speed can be displayed.
4. Press activation curve to plot the percentage of tissue activated as a function of time. Time to 100% activation is automatically displayed, while custom values for minimum (blue) and maximum (red) activation percentages between which to measure can also be selected.

7. Additional analyses and modules

1. Aside from automatically performed duration and conduction velocity analyses, several other parameters can be quantified using ElectroMap. These analyses are selectable from the dropdown menu above the display map. Select one of these options to perform the analysis, and the results will appear in the 4^th row of the results table: 1) Diastolic Interval – Time from 90% repolarization to activation time of the next action potential; 2) Dominant Frequency – Frequency spectrum of each pixel is calculated using the fast Fourier Transform, and the frequency with the most power is defined as the dominant frequency. Advanced range and window settings for dominant frequency analysis are available by selecting Frequency Mapping; 3) Time to peak – The rise time between two user selected percentages (default 10 to 90%) of the depolarization phase of the action potential or the release of calcium. Percentage values can be changed by selecting TTP Settings; and 4) Relaxation constant (τ) – Relaxation constant is calculated by fitting a mono-exponential decay of the form of the form:
   (2)
   where the fluorescence level at time t depends on the peak fluorescence, F₀, and the subsequent decay (C is a constant)²⁷. The value between which to fit equation 2 are selectable within the main ElectroMap user interfaces, as well as a goodness of fit exclusion criteria based on the r² value.
2. Press Single File Analysis to open a dedicated module for high-throughput duration and conduction analysis of each identified segment in a file. Analysis can be performed on either the whole image (duration, conduction and activation time) or on selected regions or points of interest (currently duration only). Results are outputted to a .csv file.
   NOTE: For APD values from the whole image, the first column in the .csv file is the mean, while the second column is the standard deviation.
3. Press Alternans to initiate a standalone module for dedicated analysis and mapping of beat-to-beat variability. See O’Shea et al. 2019¹³ for details on alternans processing and analysis options. Specifically, this module is designed to identify two period oscillations, known as alternans. Both duration and amplitude alternans are calculated and outputted.
   NOTE: Duration alternans are measured by comparing the duration measurement from one peak to the next; i.e. if peak one and two and APD₁ and APD₂ respectively, then the duration alternan (ΔAPD) is calculated as
   (3)
   The duration measurement is performed using the settings in the main interface. Meanwhile, amplitude alternans can be quantified and mapped across multi-cellular preparations as absolute change (defined as a percentage where 0% = same amplitude between one beat and the next). Furthermore, the effects of phenomena such as calcium load can be further investigated by measuring and comparing Load and Release alternans, as has been previously reported²⁸. If L is defined as the peak amplitude of the large beats (i.e. where the amplitude is greater than the previous beat), S the amplitude of the small beats, andD the diastolic load of the small beats, the release alternans () are defined as:
   (4)
   Conversely, load alternans () are defined as:
   (5)
   Alternans measurements can be made across the whole tissue, and the results of the analysis are displayed in the bottom right of the module. When first using the module, the analysis is performed across the entire experimental file, and the results displayed are an average beat-beat difference across the whole file. However, analysis can be restricted to specific times in the file by de-selecting Hold Zoom, zooming in on a specific time period, and selecting Analyse Zoomed Section. This will update the results panel to show analysis from the selected time period.
   1. Select play to show a beat-to-beat video of the alternans analysis. Additionally, select Create Mean Map to export a map of the alternans behavior averaged from the selecting time points, which are set in the pop-up menu when using this feature.
4. Press Phase Map to initiate the phase mapping module. A Hilbert transform is performed to calculate the instantaneous phase (between -π and +π) of the signals at each timepoint. Press play or drag slider to visualize phase behavior over time and click on a pixel to render a phase diagram.

8. Exporting data

1. Data is exported from ElectroMap in a variety of forms. Press Export Values to save the values of the currently displayed map in the main used interface. Measured values can be saved as either a map (preserving pixel locations) or condensed into a single list, and can be saved as .csv,.txt or .MAT files.
2. Press Export Map to bring up a pop-up containing the currently displayed map, which can then be saved in a variety of image formats. Display options for the map are controlled by selecting Map Settings but can also be edited once Export Map has been selected. For example, a color bar can be added by selecting this icon from the top menu, and the scale can be set by selecting Edit > Colormap.
3. Press Activation Video to render an animation of the activation sequence, which can be saved as an animated .gif file.
4. Press Segment Video to save a .avi video file of the currently displayed parameter of each identified segment.