1.

Introduction

Drosophila melanogaster has been a useful myocardial model for investigating human heart disease.^1–4 Recently, many approaches have been proposed for assessing cardiac functions in Drosophila, such as using electropacing stress,⁵ optical coherence tomography (OCT),⁶ multielectrode array systems,⁷ semi-automatic optical heartbeat analysis (SOHA),⁸ atomic force microscopy,⁹ and multiple sensor electrocardiography.¹⁰

Among these, OCT is a powerful tool for producing noncontact and noninvasive tomographic images of biological tissues. By measuring the echo time delay and the intensity of backscattered light from a specimen, OCT uses inherent differences in the index of refraction in tissue, rather than enhancement with dyes, to differentiate various tissue types. The extended imaging depth in the scattering medium enables OCT to successfully obtain in vivo images of the adult heart in Drosophila.^6,11–15

In previous studies using time-domain (TD) OCT, given the limitation of imaging speed, a number of heartbeat parameters have been obtained by manual counting; this involves tracking heart activity through a limited numbers of beats from M-mode OCT records.^6,11,12 Since the recently developed frequency-swept lasers have improved sensitivity, making real-time B-mode and Doppler OCT imaging possible,^13–15 high-throughput studies in Drosophila are now feasible. This requires an automated analyzing method for obtaining reliable measurements of dynamic heart function parameters in large samples or over a long period of time.

Therefore, in this study we proposed a rapid OCT-based method for obtaining cardiac parameters from adult flies. Our method is very useful, particularly in experiments with a large sample size. We adopted the random walker algorithm for semiautomatic cardiac chamber segmentation in B-mode OCT images. Random walker is a semiautomatic segmentation method based on graph theory by Leo Grady;¹⁶ it requires that the user give some labels as input seed points for the cardiac chamber area and background. Greater detail of random walker segmentation will be described in the next section. First, we created an initial set of seed points for the random walker algorithm on the first frame. In the next frame, the cross-sections of the heart tube are automatically segmented in a total of 2000 frames, and the size of the inner margin is represented by area for each Drosophila. Then, depending on the histogram distribution of the changing area during each heartbeat cycle, an optimum threshold can be chosen automatically by an iterative process.¹⁷ Thus, each heart period (HP) is subdivided into two meaningful sections, which are used to represent diastolic and systolic intervals (DI and SI), respectively. These two parameters are important for detailed heartbeat analysis.

The proposed method combines high-speed swept-source OCT (SS-OCT), for optical recording of beating hearts, with robust, semiautomatic analysis to rapidly detect and quantify, on a beat-to-beat basis, not only heart rate (HR) and HP but also DI and SI, and end-diastolic (EDA) and end-systolic area (ESA). These parameters can be used to determine dynamic parameters of heart function, including the arrhythmia index (AI) and percent fractional shortening (%FS). We also demonstrate age-related alterations observed in fly heart-function using this method. The convenience of the algorithm, and the relevance of our results to the human body, is also discussed.

2.

Methods

2.1.

SS-OCT Imaging Technique

Drosophila ($w^{1118}$), used as wild-type strain, were maintained on regular yeast glucose media at 25°C. For image assessment, flies were first anesthetized by exposure to Fly Nap® (Carolina Biological Supply Company, Burlington, North Carolina). Flies were then immobilized on a plastic petri dish using jelly glue, with the dorsal side facing the OCT probe. For each time point, a group of 10 flies of different ages, including one-, three-, five-, and nine-week-old flies, were examined using the SS-OCT system (OCM1300SS; ThorlabsInc., Newton, New Jersey).

Two-dimensional OCT images were obtained in the longitudinal direction to identify the conical chamber (CC) of the heart [Fig. 1(a)]. Then, the OCT image orientation was rotated by 90 deg, such that the CC was acquired in the transverse plane of the Drosophila. Figure 1(b) and 1(c) shows representative transverse two-dimensional OCT images of the CC during diastole and systole, respectively, in a three-week-old wild-type Drosophila.

Fig. 1

Representative B-mode OCT images. The arrow shows the conical heart chamber (CC) (a) in longitudinal orientation and (b) and (c) in transverse orientation. The dorsal side is at the top of the images. Squares show the heart wall in the (b) diastolic and (c) systolic stage. (d) and (e) The process of implementation of the autosegmentation algorithm. (d) Original image, (e) Enhanced image, (f) Related probability distribution, (g) Red circle indicates the heart segmentation result.

The SS-OCT system had a median wavelength of 1310 nm; an axial and transverse resolution of around 9 and 15 μm, respectively, in tissue; total power of 10 mW; and an A-scan rate of 16 kHz. For each measurement, around 2000 frames (each covering an area of $0.75\times 3{\mathrm{mm}}^2$, corresponding to $128\times 512\text{pixels}$ in the $Y-Z$ direction) was obtained at 92 frames per second.

2.2.

Algorithm

We used a combination of two algorithms, written in Matlab, to automatically measure several heartbeat parameters in Drosophila. First, after contrast-enhancement of the original image, as shown in Fig. 1(d) and 1(e), an initial set of seed points are created by the user to indicate the inner margin (green point) of the heart chamber for the random walker algorithm. The red point is automatically provided by the algorithm, showing seed points at fixed positions to indicate the outer margin of the heart chamber and background. When the random walker moves from any pixel, it calculates a vector of probabilities, based on the correlation of both grey-level and distance for each pixel, to reach the specified seeds. As demonstrated in Fig. 1(f) and 1(g), the probability of every unlabeled point reaching the labeled points will be computed by this algorithm. Then, segmentation results can be obtained by performing thresholding on the probability map. The algorithm is summarized below, and more details can be found in Ref. 16.

Based on graph theory, a graph consists of a pair of points, $G=(V,E)$, with vertices set $v\in V$ and edge set $e\in E$. The weight of an edge $e_{ij}$ has a weight $w_{ij}$ between neighboring vertices $v_i$ and $v_j$, with intensity $g_i$ and $g_j$. The weighting function is given by

Eq. (1)

$$w_{ij}=\exp [-\beta {(g_i-g_j)}^2],$$

where $\beta$ is the constant for adjustment of $w_{ij}$. The discrete Dirichlet integral in a region $x$ is defined as

Eq. (2)

$$D[x]=\frac{1}{2}{x}^{T}Lx,$$

where the Laplacian matrix $L$ is defined as

Eq. (3)

$$L_{ij}=\{\begin{array}{ll}{d}_i,& \text{if}i=j\\ -w_{ij},& \text{if}{v}_i\text{and}{v}_j\text{are adjacent nodes}\\ 0,& \text{otherwise}\end{array},$$

in which $L_{ij}$ is indexed by vertices $v_i$ and $v_j$ and $d_i$ is the degree of vertex $v_i$. All pixels of the same region have a similar intensity, so that differentiating $D[x]$ approximates zero. Given the labeled points $x_M$, the unlabeled points $x_U$ can be solved by the discrete Dirichlet integral function

Eq. (4)

$$D[x_U]=\frac{1}{2}[x_M^T{x}_U^T]\left[\begin{array}{cc}{L}_M^T& B\\ B^T& L_U^T\end{array}\right]\left[\begin{array}{c}{x}_M\\ x_U\end{array}\right].$$

The unlabeled area, $x_U$, can be obtained by solving a sparse linear system in the differentiating $D[x_U]$ as follows:

Eq. (5)

$$L_u{x}_U=-B^T{X}_M.$$

It is only necessary to position seeds on one frame; thereafter, the edges of the CC are automatically detected on every frame. Figure 2 demonstrates the area inside the boundary of the heart tube in 100 frames, where peak and trough positions represent EDA and ESA of the Drosophila. Here, we define $\mathrm{FS}=(\mathrm{EDA}-\mathrm{ESA})/\mathrm{EDA}$, which represents the extent of change in the cardiac chamber area during systole; it can provide an estimate of the contractility of the heart tube. HP is quantified as the time elapsed between the ends of two consecutive systoles (i.e., at the troughs of the filtered data, as indicated by the red line in Fig. 2). AI is defined as the standard deviation of HP, normalized to the median HP. This parameter provides quantitative representation of the arrhythmicity of the hearts of individual flies, which is useful for comparing groups of individuals.

Fig. 2

Change in area inside the segmented boundary of the heart tube during heart contraction. Blue: original data, red: filtered data for determining positions of troughs only. HP is the time between two troughs.

Furthermore, in order to automatically estimate the length of diastole and systole, we adopted an automatic thresholding-algorithm.¹⁷ First, we used the average value of the maximum and minimum areas on the histogram distribution during each heartbeat cycle as an initial estimate of the threshold $T$. Then, we partitioned the area in each HP into lower and higher groups, R1 and R2, using the threshold $T$. A new threshold is then created by calculating the mean values ${\mu }_1$ and ${\mu }_2$ of the two partitions, using the equation

Eq. (6)

$$T=\frac{1}{2}({\mu }_1+{\mu }_2).$$

These three steps are repeated until the mean values, ${\mu }_1$ and ${\mu }_2$, do not change in successive iterations; the optimum threshold can then be chosen automatically as a result of this iterative process. In the process, each HP is subdivided into two meaningful regions; Fig. 3 demonstrates two heartbeats, where SI is quantified as the interval during contraction (systole), which is, in turn, defined as the time occupied when a cardiac area smaller than the mean ${\mu }_1$ values of the lower group R1 are encountered. DI is then quantified as HP minus SI.

Fig. 3

Change in area inside the segmented boundary of the heart tube in two heart periods. Horizontal lines represent the automatic thresholding result. T, optimum threshold; ${\mu }_2$, the average of the upper group; ${\mu }_1$, the average of the lower group.

3.

Results

We first tested how inputting different numbers as seed points and in different regions could affect the probability distribution and the final segmentation. The representative result is shown in Fig. 4. For comparison, Fig. 4(a) shows the segmentation results manually measured by our expert. Figure 4(b) to 4(d) shows the segmentation results from inputting one, two, and three seed points, respectively, to indicate the inner margin. The right column of Fig. 4 shows the enlarged images. For each image, we input seed points in 20 different randomly selected regions. Random walker segmentation resulted in the same mean area size ($5419{\mu \text{m}}^2$), irrespective of whether one, two, or three seed points were provided. The coefficients of variation (CV), which are defined as the standard deviation of area in percentage, normalized to the mean area, for seed points provided in different regions were 0.13%, 0%, and 0% for inputting one, two, or three seed points, respectively. This indicates that the algorithm has very high precision in segmenting the Drosophila heart tube in OCT images.

Fig. 4

Segmentation results from (a) manual segmentation by our expert, (b) to (d) semi-automated method, by inputting one, two, and three seed points, respectively. Right column shows the enlarged images.

Thereafter, we randomly selected 28 images (each with two seed points) and compared the size of the segmentation area with that obtained by manual segmentation. As shown in Fig. 5(a), the correlation coefficient between automated method and manual segmentation was 0.993, which demonstrated that the random walker algorithm is as accurate as the manual method in segmenting the Drosophila heart tube in OCT images.

Fig. 5

(a) Correlation of the segmentation results of area ${\mu \text{m}}^2$ between the automated and manual methods. (b) Canonical correlation results of DI and SI using the automated thresholding and manual methods. R, correlation coefficient.

In order to assess whether SI and DI detection, based on the automated thresholding algorithm, differs from that defined by manual interpretation by an expert user, manual measurements for DI and SI were obtained by the user adjusting the thresholds on the histogram distribution of a segmented area (e.g., in Fig. 2, where five representative DI and SI were measured for each fly). Using the contraction pattern provided by M-mode OCT as an objective assessment of heart wall motion, the user could then finely adjust the thresholds. We randomly selected four movies taken of one-, three-, five-, and nine-week-old flies and compared the manual measurements obtained in this manner with the output of our algorithm. Both DI and SI periods in the heart of Drosophila were measured manually and using the automated thresholding algorithm. We evaluated the correlation of the manual and automated methods using canonical correlation, a statistical analysis method used to assess correlation of two sets of variables with more than one dimension.¹⁸ Figure 5(b) shows that the maximal canonical correlation found was 0.9989, which indicates that the algorithm defined SI and DI as accurately as did manual measurements.

Figure 6 shows the M-mode contraction pattern, in which we can observe the phenomenon of prolonged relaxation time and irregular heart rhythms that occur with increasing age. In particular, flies in their ninth week of age displayed various forms of arrhythmia. Figures 7 to 10 summarize eight heartbeat parameters of male Drosophila $w^{1118}$ derived and quantified by our automated algorithm, in their first, third, fifth, and ninth week, respectively. Data points represent the mean [±standard error of the mean (SEM)] for 10 files per datum point. Figure 7 shows a significant difference both in HP and HR between one- and three-week-old [analysis of variance (ANOVA), $P=0.003$)], one- and five-week-old (ANOVA, $P=0.001$), one- and nine-week-old (ANOVA, $P=0.001$), three- and five-week-old (ANOVA, $P=0.001$), and three- and nine-week-old flies (ANOVA, $P=0.001$). It was also quite clear from our data (Fig. 8) that DI in particular increased with age, as compared to SI, with a significant difference in DI between three- and five-week-old flies (ANOVA, $P=0.001$). Moreover, the variation in DI was more obvious when comparing flies in their fifth and ninth weeks; this was consistent with Fig. 9, which shows the average AI, reflecting the age-related increase in arrhythmicity in flies (ANOVA, $P=0.007$). However, the average size for ESA remains similar across age, while a small but statistically insignificant increase in EDA is shown in Fig. 10(a). For this reason, no significant decline in FS was observed, as shown in Fig. 10(b).

Fig. 6

A representative M-mode OCT images in male Drosophila $w^{1118}$ at one, three, five, and nine weeks of age.

Fig. 7

Cardiac parameters in male Drosophila $w^{1118}$ at one, three, five, and nine weeks of age showing (a) HP and (b) HR.

Fig. 8

Cardiac parameters in male Drosophila $w^{1118}$ at one, three, five, and nine weeks of age showing DI and SI.

Fig. 9

Cardiac parameters in male Drosophila $w^{1118}$ at one, three, five, and nine weeks of age showing AI.

Fig. 10

Cardiac parameters in male Drosophila $w^{1118}$ at one, three, five, and nine weeks of age showing (a) EDA and ESA and (b) FS.

4.

Discussion

Drosophila is an important model for investigating human heart diseases. Automated quantification of potentially relevant parameters from the model can replace labor-intensive manual work and generate more objective, accurate, and comprehensive metrics. A widely used approach for cardiac image analysis is to reduce complexity by using M-mode, which is the monodimensional projection of the image along a line. M-mode images are often used in echocardiography for the analysis of heart motion dynamics.¹⁹ Previous methods^6,11,12 using M-mode OCT images were based on continuous in-depth scans in the midline of a fly’s cardiac chamber over time; cardiac parameters, including HR, end-diastole diameters, and end-systole diameters, are manually measured and derived by averaging a few cardiac cycles during a normal, regular rhythm.

Automatic segmentation of cardiac M-mode images for identifying end-systole and end-diastole can be a challenge, although many methods have been proposed, such as gray-level histograms,²⁰ fuzzy set,²¹ Markov random field models, neural networks, binarization, morphology,²² and many others. In these methods, borders are not well defined and the gray level is not uniform. Thus, Bertelli et al.²³ proposed a semi-automated method, which segments cardiac M-mode images by using a multiclassifier for computing cardiac parameters, such as the ratio between the diameters in end-diastolic and end-systole, or the ratio between diastolic and systolic durations. The required computational time for the process is about 25 s. Fink et al.²⁴ used a combination of two movement detection algorithms to track movement of heart edges from a semi-intact Drosophila heart. The user needs to adjust the thresholds and filters to finely tune the algorithm output, using M-mode as an objective assessment of heart wall movements.

In this study, instead of following the segmentation of cardiac M-mode images, cardiac parameters are measured on the basis of changes in the whole cross-sectional area of the cardiac chamber, rather than by measuring the brightness (i.e., gray level change) between the edge of the superior and inferior walls during mid-diastole and mid-systole in M-mode OCT images. Our proposed semi-automated segmentation method, based on the random walker algorithm, can provide a unique, quality solution that is robust to weak/noisy object boundaries.¹⁶ This is because, according to the theory of random walker segmentation, the performance of the algorithm depends on the weighting function in Eq. (1). Noise is usually considered to be a disturbance in a local region. This means that a region with noise has a greater intensity difference between adjacent pixels than a region without noise. The weighting function determines the influence of intensity difference between adjacent pixels. It decreases the effect of noise and increases the boundary effect. Further details and proofs can be found in Ref. 25.

There are some other well-known segmentation methods, such as active contours or graph cuts. Segmentation based on active contouring²⁶ has large time complexity and poses a local energy minimization problem, so that the method does not always present a correct solution, especially when segmenting images with a tortuous boundary or high noise levels. Figure 1 shows images with a markedly tortuous boundary and high noise level. Another well-known algorithm, graph cuts,²⁷ has been developed for automated and semi-automated segmentation. For automated segmentation, the performance of graph cuts is noticeably better than other methods. However, if an insufficient number of seeds is provided in the semi-automated segmentation, the “small cut” problem will decrease the correctness of segmentation. Grady¹⁶ showed that random walker segmentation can cope with images with a tortuous boundary or high noise levels, has no local energy minimization problem, or “small cut” drawback, even if a small number of seeds are used. Our data confirmed that even in the presence of a degraded boundary, such as the cardiac tube of Drosophila, the algorithm could perform segmentation successfully (with a very high correlation to manual segmentation, $R=0.993$). We also proved that when inputting more than two seed points, different regions of input would not affect the final segmentation result ($\mathrm{CV}=0$).

The automated thresholding algorithm is a simple and efficient statistical method for computing a threshold value automatically between two bimodal distributions, simply according to their mean value. In this study, we observed that the area distribution in diastole and systole followed a bimodal distribution. Thus, DI and SI, based on our algorithm, could be automatically and independently detected in each HP. This method was validated against manual interpretation by our expert, and the algorithm generated the SI and DI as accurately as did manual measurements (the methods were highly correlated, $R=0.998$). It should be noticed that there are periods of time when the speckle pattern is fixed on M-mode OCT images (Fig. 6); the rate of change of the heart segmented area at these times is zero, so that our automated thresholding algorithm detects “long” diastoles or systoles.

The quantitative heartbeat parameters measured using our proposed method show that as flies age, HP lengthens [Fig. 7(a)], resulting in an age-related decline in HR [which is the inverse of HP; Fig. 7(b)]. We also found that the age-dependent decrease in HR may be due to a disproportionate increase in DI, compared to SI, with age (Fig. 8). Moreover, the current algorithm not only measures HR but also quantitatively expresses the age-dependent increase of heartbeat arrhythmicity (i.e., AI), which may be due to large variations in DI associated with age, particularly when comparing flies in their fifth and ninth weeks of age (Fig. 9). These findings are consistent with results obtained with a movement detection algorithm by Fink et al.²⁴ Ocorr et al.⁷ also suggested that age-related cardiac dysfunction in the fly heart is due to a decrease in the efficacy of cardiac relaxation. Since age-dependent decreases in intrinsic HR, and increases in the incidence of AI, have also been documented in humans,^28–31 our observation of the age-related changes in heart function of flies has correlates in humans.

Large variations existed in both EDA and ESA within the tested groups, which may be due to differences in size between flies; however, use of %FS can eliminate the effect of body size. Fink et al.²⁴ have reported age-dependent decreases in FS. However, in our study, no significant declines in %FS in fly hearts were observed with age, suggesting that muscle contractility remained similar across age. Previous studies using different cardiovascular stress methods to estimate end-systolic and end-diastolic dimensions also failed to detect significant alterations associated with aging.^32,33

In this study, because the eight heartbeat parameters were automatically detected by manually positioning two seeds only on one frame, it only requires around 0.5 s to process a single fly (including a minimum of 100 cardiac cycles) with the proposed algorithm; this can be an efficient and useful approach for high-throughput gene-screening studies in Drosophila. Moreover, solving the solution of the random walker segmentation according to Eq. (5) is a linear system. Douglas et al.³⁴ have demonstrated the design of analog circuit implementation from a linear system. This direct correspondence with analog electrical circuits opens the possibility for hardware implementation of the random walker segmentation algorithm. As the random walker algorithm can also operate in three dimensions, further research, combining faster SS-OCT scanning and more control of fly movement, may allow heart volume segmentation.

5.

Conclusion

Here, we presented an algorithm for cardiac image segmentation with high precision and high accuracy, based on a small set of pre-labeled pixels. Furthermore, the technique we described here allows an optimum threshold to be chosen automatically as a result of an iterative process. By applying the semi-automated segmentation and auto-threshold algorithm, we established an easier and more accurate quantification method for measuring heart parameters in Drosophila. Both algorithms involve the use of a few empirical parameters and are expected to perform robustly in the presence of noise. Using our approach, we were able to show age-dependent changes in fly heart-function, including HP increases; disproportionate increases in DI and SI, increases in AI, and a similar %FS across age. The combination of high-speed and successive B-mode imaging features of SS-OCT gives this technique potential for automatically assessing large numbers of individuals over significant periods of time. OCT-based automatic and rapid quantification of cardiac parameters can facilitate future high-throughput studies in Drosophila.