Automatic Detection of Larval Zebrafish ECG: Computational Tool for High-throughput Cardiac Activity Analysis

Automatic analysis of larval zebrafish electrocardiographs (ECG) is essential for high-throughput measurements in environmental toxicity assays, cardiotoxicity measurements and drug screening. We have developed a MATLAB based software is built on methods that have previously been used to analyze human ECG, such as the Pan-Tompkins algorithm and Wavelet. For the first time these sophisticated tools have been applied to larval zebrafish ECG to automatically characterize the heart-beat waveforms. The ability of the automated algorithm to detect the QT interval for both normal and pharmacologically altered larval ECG is found and compared to previously used software that is built into LabChart® (AD Instruments). Using wavelet transforms it is shown that the typical larval ECG features are within the frequency range of 1 to 31 Hz. It is also shown that the automated software is capable of detecting QTc (heartrate corrected heartbeat interval) changes within pharmacologically altered zebrafish larval ECG. The automated process is a significant improvement on the approaches that were previously applied to the zebrafish ECG. The automated process described here that is based on established techniques of analyzing ECG can sensitively measure pharmacologically induced changes in the ECG. The novel, automated software is faster, more sensitive at detecting ECG changes and less dependent on user involvement, thus minimizing user error and human bias. The automated process can also be applied to human ECG.

Introduction 6 118 From this average view, the user defines the position of the Q-start point, Q-end point, T-peak and T-end.
119 From this input the software automatically calculates the QTc (based on the Bazett formula) and other 120 heartbeat characteristic. This information was exported as a table into Microsoft Excel for further analysis.
121 Performing the analysis on a section of data takes at least 5 minutes.
122 This process was performed on each section of larval ECG to determine the heartbeat characteristics before 123 and after the introduction of the pharmacological agent.

124
125 Applying the Automated process 126 The automated process was applied to each data section and the outputted QTc was then exported to Excel.
127 The program measured the QTc for each 40 consecutive heartbeats and the output was tabulated as shown 128 in Table I 177 The original architecture of the Pan-Tompkins algorithm is divided into three processes which can be 178 thought of as the Initial Learning Phase (ILP), Secondary Learning Phase (SLP), and Detection Phase (DP).
179 The ILP initializes detection thresholds based upon the size of the "signal" and "noise" peaks detected. The 180 SLP uses two full heart cycles to determine the average RR interval and then set the limit of the possible 9 181 RR-interval for the rest of the ECG trace. To allow for any adaptation of the recorded ECG signal, for 182 example due to pharmacologically induced change in heart rate, these thresholds are adjusted periodically.
183 The DP processes the ECG and generates a pulse for each QRS. It then uses thresholds based on both the 184 filtered ECG and the processed signal to detect the QRS waves. Their detection thresholds were set to just 185 above the "noise" that is sensed by the algorithm. This approach reduced the number of false positives 186 caused by noise that mimics the QRS characteristics, which is always a problem in human ECG due to 187 electromyography artefacts. In zebrafish these artefacts could also be a problem due to sporadic motion 188 artefacts caused by twitching or other motion. The automated process conserves and builds on this 189 architecture. 2) The second stage of the algorithm is to differentiate the signal to accentuate the turning points. If 210 all artefacts have been removed then the turning points can only be associated with biological events, which 211 for a normal ECG signal occur at the P-wave, R-wave and T-wave. By differentiating, these features 212 become amplified. The amount that they are amplified is proportional to their frequency as higher 213 frequency activity is changing at a faster rate. 11 214 215 As R-waves have the highest frequency of the ECG features it is accentuated the most in the differentiated 216 output (Fig 2A).

217
218 The differentiation is applied via the following iteration, Where is the differentiated output and is the filtered ECG signal. In MATLAB this can again be 221 applied using an inbuilt function for which the code is given below, 252 253 Where is the integrated ECG output. If Z is roughly equal to the width of the QRS peak then it has 254 the effect of creating large peaks in the region of the QRS, whilst ensuring that everything else is roughly 255 zero. The output of this algorithm is shown in Fig. 2A. 268 The algorithm first searches through the integrated waveform and an R wave is 'detected' every time a 269 peak is found above the established threshold level. If the detected peak is below threshold, it is treated as a 270 noise peak. If an R wave is not detected within 166% of the previously measured RR interval the program 271 performs a search back using the second, lower threshold. Every time a peak is detected the algorithm 272 determines if it is an R-wave using the previously established thresholds and updates the thresholds by the 273 following algorithm: 274 If the peak that has been found is an R-wave then the new running estimate of the R peak height (R PK ) is 275 updated as, 276 R PK = 0.125*PEAK + 0.875*R PK 277 Where PEAK is the height of the detected peak.
278 If the peak that has been found is below threshold and therefore a noise peak, then the new running average. This enables the mean voltage at every position in the heart beat, relative to the R wave to be 309 determined using the formula below, For the automated process the average waveform is calculated for a period that is equal to the measured RR 312 interval that runs from 0.15 s before the R peak. This is shown in Fig. 2D. 317 To find Q the software searches for the local minima in the region that is 0.15 s before the R peak. This 318 window is guaranteed to contain a Q-wave if the R-wave frequency is greater than 6 Hz.
319 To find the T Peak the software searches for the local maxima that has the largest amplitude after the R 320 peak. Once the T-peak has been found the program searches for the next point at which the ECG changes

334
A wavelet transform decomposes a signal into its fundamental parts that have a well-defined time 335 and frequency localization. Through a convolution, the continuous wavelet transform (CWT) finds the 336 inner product between the inputted signal and the analyzing wavelet that has a well-defined time-duration 337 and frequency band. In a CWT the signal is compared to the analyzing wavelet that is time-shifted and 338 scaled to yield coefficients that correspond to a measurement of the ECG constituents within the section 339 and frequency band. In essence, the wavelet transform provides information about the ECG frequency at