1. Introduction

Investigation of the relations between changes in retinal structural and biomechanical properties and various diseases has been of increasing interest over the last decades. A large number of studies reported on the prospects of identifying retinal alterations in patients with degenerative or other cerebrovascular and cardiovascular diseases to define new biomarkers in order to improve clinical screening and monitoring [1–5]. Compared to other organs of the body, the eye has indeed the substantial advantage of being easily and non-invasively accessible in-vivo using optical coherence tomography (OCT). This imaging technique is already a gold standard in ophthalmology for the diagnostics of ocular diseases. While the different retinal layers are visualized in OCT intensity images, the phase data provide additional valuable information, especially for the analysis of the retinal vasculature. In several clinical applications, OCT angiography (OCTA), an extension of OCT, is increasingly used to examine the retinal vascular network [6]. In this label free method, the vascular information can be accessed either from the intensity and/or phase data by differentiating between static tissues and moving red blood cells. Doppler OCT, another extension of OCT, additionally allows to quantify blood flow from the phase difference between successive A-scans [7].

Phase-sensitive OCT (PhS-OCT) is not only a valuable technique to study fast motion in the vasculature itself but may also be applied for an analysis of slower rate length changes at the nanoscale in the eye [8,9] or for an investigation of tissue deformation in other tissues and organs [10–12]. PhS methods are commonly used in optical coherence elastography (OCE) as well, another extension of OCT that additionally provides access to tissue biomechanical properties, especially the tissue elasticity [13,14]. In OCE, tissue motion is usually accessed by determining differences in the amplitude and/or the phase between consecutive A-scans or B-scans [15–17] and inter-frame deformation can directly be obtained using the finite difference method (FDM) [18]. Phase gradient methods, that are based either on the least square method (LSM) [18] or on the vector method (VM) [19], are however more commonly applied for the estimation of the inter-frame strain in OCE, as the sensitivity of such methods is higher compared to the FDM [18]. Depending on the motion rate, unwrapping, which may be subject to error, may be required for most phase sensitive methods prior to strain estimation. This step is avoided in the VM as complex-valued signals are considered. The cumulative strain relative to the initial frame can finally be obtained by summation of the individual inter-frame data. The estimation of cumulative strain by such incremental method has been further optimized for slow rate tissue deformation imaging, as shown in recent studies [20,21]. In these works, rather than performing strain calculation for each pair of consecutive scans, the inter-frame time interval is adjusted according to the speed of the deformation. Another method used for the evaluation of tissue motion is based on the calculation of relative phase difference over time between two retinal layers and to compare this difference to the value of a fixed initial volume, as demonstrated for functional imaging in the retina using a full field (FF-) swept-source (SS-) OCT system [9,22]. The total issue displacement is then retrieved after temporal integration of the phase difference data. In non-moving (e.g. excised) tissue, a fixed reference is sometimes chosen for the phase difference calculation, as is done for example between A-scans in photo-thermal OCT [23].

Various studies investigated the use of PhS-OCT and OCE methods for ophthalmic applications. In OCE, investigations on corneal imaging, either in-vivo in animal models or ex-vivo [17,24–27] have mostly been reported. In recent years, in-vivo human corneal OCE imaging was also demonstrated [28–30], enabling to calculate the elastic modulus in human subjects in-vivo which was not possible with other existing techniques probing the mechanical properties of the cornea [24]. As non-invasive OCE measurements are challenging, reliable quantitative evaluation for clinical translation is up to now still limited [31]. OCE measurements in retinal tissue pose even greater challenges, and only ex-vivo samples [32,33] as well as in-vivo animal models [34] have so far been imaged. While full elastography requires information on the force exerted on the tissue, which is difficult to access in-vivo, passive elastography is a promising alternative technique used to access biomechanical properties of tissues without the need of an active external source. Passive excitation can result from the cardiac beatings or muscle contractions, among other examples. In-vivo OCT investigations in animal models have already shown the potential of the method for the evaluation of both corneal stiffness [35] and corneal tissue response [36] under natural pulsatile motion. While the first study focuses on the analysis of mechanical shear wave propagation, the second study is based on the calculation of tissue deformation and strain using the LSM over time. In-vivo heartbeat induced corneal strain measurements were also recently demonstrated in human eye using ultrasound elastography [37]. Similarly, several human in-vivo studies using PhS-OCT also focused on the investigation of deformation in corneal or retinal tissue under pulsation. Spahr et al. have shown retinal nanoscale tissue motion in the surrounding of retinal vessels over several cardiac cycles imaged with a FF-SS-OCT [38] (using the phase difference method [9,22]). Measurements of the trabecular meshwork motion [39,40] and the lamina cribrosa motion [41] under pulsation, particularly of interest for glaucoma research, were also conducted using PhS- spectral domain (SD-) OCT systems. In these two studies, tissue motion was retrieved from the phase difference calculated between adjacent B-scans only, providing therefore information on tissue deformation that is large and fast enough to induce measurable phase changes on a short time scale in the range of milliseconds.

Successful investigation of tissue deformation in the eye, either with OCT or OCE, is highly dependent on phase stability. Different optimization approaches to achieve higher phase stability have been investigated by, for example, studying the instabilities of swept sources [42], increasing the imaging speed [43] or testing common path configurations [17,44]. In ophthalmic applications, the use of real time eye trackers, usually employing a scanning laser ophthalmoscope (SLO) integrated to an OCT system, has also greatly improved the stability of in-vivo measurements. Different tracking techniques have been implemented both in commercial and research OCT instruments, as seen in the literature [45–47]. Such technical improvements have fostered the use of OCT imaging for new applications.

In this work, we present the results of imaging and measuring deformation induced in the retinal nerve fiber layer (RNFL) by retinal vessels pulsating with the heartbeat in healthy human subjects. Measurements were done using a custom-built SD-OCT system with an integrated retinal tracker that stabilizes imaging with sub-speckle precision. To detect the subtle deformation induced changes, the phase information provided by OCT is used. Two methods were used for the evaluation of phase differences over time, both providing access to the temporal phase evolution (TPE) of data over several seconds. The phase data are converted into displacement data and displayed as maps that show tissue displacement vs. time at a given depth and lateral position within tissue. Thereby, tissue deformation on the scale of a few nm and over extended time periods of several seconds is demonstrated. We also analyze differences in tissue deformation induced by arteries and veins and report on the repeatability and precision of the method.

2. Methods

2.1 Temporal phase evolution OCT

In this study, we aim to measure sub-micron tissue motion over an extended timeframe, i.e. over several seconds. OCT based motion measurements typically compare phase signals recorded at the same position with a time lag Δτ between them. Depending on the motion speed to be measured, and depending on whether absolute or relative motions are to be analyzed, different evaluation strategies are used. The phase calculation can either be performed between A-scans (typically used for Doppler OCT [7], Fig. 1(a)), between B-scans (typically used for OCTA and also for OCE [24,48], Fig. 1(b)), or between volumes, and may consider a fixed reference, a moving reference, or referencing to all other frames simultaneously (as proposed by the extended Knox Thompson algorithm [49,50]). The standard Doppler method is not sensitive to bulk motion because it uses time lags in the order of ∼ 1–100 µs and therefore is only sensitive to rather fast motions that occurs in large retinal vessels. However, the range of measureable velocities can be increased using different scanning protocols as demonstrated in [51]. In OCTA, Δτ ∼ 1–100 ms is typically used for contrasting slow motion within capillaries This method is however more sensitive to both axial and lateral bulk motions, and quantitative information on slowly moving tissue components is difficult to achieve.

 

Fig. 1. Comparison between the TPE OCT methods (M1 and M2) used to access motion in the range of seconds and the usual PhS methods used for example in standard Doppler OCT and OCTA to evaluate motion in the range of microseconds and milliseconds, respectively.

Download Full Size | PDF

In order to analyze slower motion happening over several seconds, such as tissue deformation, we study the evolution of the optical phase signal over time with respect to a fixed reference frame (for which we introduce the term TPE-OCT). We used two different methods that also provide information on absolute tissue positions relative to a starting point. In the first method (M1), the phase differences between a reference B-scan and each of the subsequent B-scans of the acquired sequence, providing a Δτ of ∼ 0.1–10 s (Fig. 1 (c)), are calculated. As the SD-OCT system used for imaging has an integrated retinal tracker providing sub-speckle precision [46,52] (c.f. the preserved speckle structure in the averaged B-scans of Figs. 2(b), (e), (h)), the impact of lateral motions is greatly reduced in our measurements. In the second method (M2), the phase difference between consecutive B-scans is calculated followed by temporal integration relative to an initial reference frame. This second method has been used for example, in OCE and FF-OCT [36,38]. If phase unwrapping is performed for the data provided by M1, both methods are mathematically equivalent [53], however, each of them has slight advantages for some specific data presentation schemes, so we will demonstrate results obtained with both of them (in the majority of the presented results, we used M1).

 

Fig. 2. Overview of B-scans and their locations for the three study subjects. (b), (e) and (h) are the averaged intensity images of the traces marked by red dotted lines on the LSLO images (a), (d) and (g) for subjects S1, S2 and S3, respectively. The diameter of the vessels is calculated between the vessel walls as marked by the arrows in the zoomed images (c), (f) and (i). Some brighter walls, marked in red, can be more easily identified than others, marked in blue. The trace marked in orange in (d) is shown in Fig. 6.

Download Full Size | PDF

2.2 Measurement protocol

2.2.1 Experimental set-up

A custom built SD-OCT with an integrated retinal tracker is used for imaging [52]. It operates with a superluminescent diode (SLD) with a center wavelength of 860 nm and a bandwidth of 60 nm at an A-scan rate of 70 kHz and has a sensitivity of 98 dB. The resolution achieved is 4.2 µm axially (in tissue) and 20 µm laterally. The retinal tracker, essential feature for the purpose of this study, is composed of a line SLO (LSLO) with a laser light source at a center wavelength of 786 nm. A GPU-based algorithm is used for real time tracking as described in previous literature [46]. The power to the eye is 0.7 mW from the LSLO and 0.5 mW from the OCT beam, which satisfies the laser safety regulations.

2.2.2 Acquisition procedure

We are interested in measuring RNF tissue movement induced by retinal vessels that pulsate with the heartbeat, and to analyze the differences between deformation caused by the veins and the arteries. For each subject, we thus chose to acquire cross sections of RNF tissue in the vicinity of retinal vessels, including at least one artery and one vein. The B-scan traces, which can be of arbitrary orientation, were manually selected before imaging on the template LSLO image as shown in Figs. 2(a), (d), and (g). The coordinates of the selected points were then converted to input amplitudes for the x and y OCT scanners after linear interpolation at evenly spaced positions along the trace, as described in [54]. In order to obtain images with sufficient RNFL thickness and rather large vessels, ocular fundus regions nasally inferior or superior to the macula were preferably chosen. The lengths of the selected traces were about 1.6 to 2.2 mm and were adjusted to include sufficient RNFL tissue around the vessels while keeping a rather high sampling density and minimizing the effect of head rotation that will be more pronounced for longer traces. For each measurement, 1000 B-scans of 512 A-scans were recorded at the same position corresponding to an acquisition sequence of about 10 s.

2.2.3 Subjects’ selection

Three eyes of three healthy volunteers (aged 30, 34 and 32, respectively S1, S2 and S3) were included in this analysis. The study was approved by the university’s ethics committee and are in agreement with the tenets of the Declaration of Helsinki. The averaged intensity images of the traces as well as their location on the LSLO are shown in Fig. 2.

2.3 Data processing

2.3.1 Standard processing and registration

Standard Fourier Domain (FD-) OCT processing of the raw data was used to obtain the individual 1000 intensity and phase B-scans for each recorded data set. In order to correct for axial and lateral motions, registration was performed on the intensity scans, taking the first intensity B-scan as reference and using a similar method as described in [55]. The phase B-scans were registered by using the same shifts as for the intensity scans. Figure 3(a) shows an averaged intensity B-scan image after registration and flattening (with reference to the anterior RNFL surface), covering three cardiac cycles. The positions of two arteries and one vein are marked by arrows. The preserved speckle pattern provides a first impression of the high quality of tracking and registration. Figure 3(b) shows an M-scan (distance along B-scan trace vs. time) of the top boundary of the RNFL over time. Straight horizontal lines indicate individual bright speckles over time. Apart from a saccadic eye movement, the maximum lateral shift after registration was ∼ 2 pixels. Figure 3(c) shows the cross-correlation coefficient of each intensity B-scan with respect to the reference frame after registration. Values of ∼ 0.7 or larger can be considered as good.

 

Fig. 3. Illustration of image registration results. (a) Averaged intensity B-scan, covering three cardiac cycles, after registration and flattening. Two arteries (A1, A2) and one vein (V) are marked by arrows. (b) M-scan (vertical: distance along B-scan trace; horizontal: time) of the top surface of the RNFL. (c) Cross-correlation coefficient of intensity B-scans after registration (with respect to reference frame). The red circle marks the systolic peak at which images of Fig. 4 (a)-(c) are taken, the black circle marks the position where Fig. 4 (d) is taken.

Download Full Size | PDF

The diameter of each vessel was extracted from the averaged intensity scan, considering the number of pixels between the upper wall and lower wall in depth (see the arrows in Figs. 2(c), (f) and (i)). As the vessels may not be imaged perfectly perpendicularly and their walls cannot always be clearly distinguished, a margin of error is to be considered on the estimated diameter values. Standard Doppler images were calculated to identify the different heartbeat cycles from the phase variations in arteries over the full acquisition. For each cardiac cycle, the subset of intensity and phase B-scans were extracted from the entire sequence for the subsequent analysis. In the case of a Doppler angle close to 90°, the velocity could not be retrieved from the standard Doppler (see Section 4).

2.3.2 Processing of the phase data over time

In order to analyze RNF tissue changes over individual or several consecutive cardiac cycles, i.e. motion in the range of seconds, we used the two TPE methods, M1 and M2, introduced in section 2.1. For the evaluation of a single heartbeat sequence or a sequence of consecutive heartbeat cycles, several post-processing steps were followed as described below for M1:

Step 1. The first frame (phase and intensity data) of the heartbeat cycle was selected as reference for all subsequent steps. This corresponds to the beginning of the systolic phase. In the case of consecutive cardiac cycles, the reference frame is the first frame of the first heartbeat cycle for the processing of the entire sequence. The influence of the selection of the reference frame on the evaluated data is outlined in Section 4.

Step 2. Registration was then performed: the shifts in x and z direction for each scan compared to the reference scan were obtained using the intensity data and applied to the phase data.

Step 3. Preliminary phase difference tomograms between the reference phase B-scan and each subsequent phase B-scan (at each recorded instant of the heartbeat cycle) were calculated.

Step 4. Bulk motion correction at each A-scan position was calculated on each preliminary phase difference scan to remove the remaining small amplitude motion. Two methods were implemented: either taking the mean phase difference value over the entire A-scan depth or taking the mean value over the pixels at the IS/OS junction only (as marked between the red dotted lines in Fig. 2(b)).

Step 5. The final bulk motion corrected phase difference tomograms between the reference phase B-scan and each subsequent phase B-scan were calculated.

For M2, the same processing steps are followed using a moving reference, i.e. taking the phase difference tomograms between adjacent B-scans, instead of a fixed reference (here the first frame of a heartbeat cycle). The sequence of deformation B-scans over time is obtained by temporal integration of the individual phase difference scans starting with the reference scan followed by conversion to tissue deformation in µm using Eq. (1).

Examples of the phase difference tomograms (obtained using M1) at the systolic peak of a cardiac cycle using each of the two bulk correction methods are given in Fig. 4(a) and (b). These phase difference tomograms correspond to the difference between the reference frame and the 15^th B-scan after the reference frame, corresponding to a time lag of Δτ = 0.14 s. For comparison, a phase difference tomogram between two adjacent B-scans (Δτ ⋍ 9 ms, corresponding to the procedure of OCTA imaging) is shown in Fig. 4(c). While Figs. 4(a), and (b) indicate absolute movements (with respect to the reference frame), the color coding of Fig. 4(c) is proportional to the speed of tissue motion. Compared to the latter case, large changes with several wrappings are seen in the RNF tissue surrounding the vessels over the entire scan taken with the long Δτ values. The phase difference data corresponding to the A-scans containing the large vessels are unreliable because of the low intensity in these parts of the images that prevent a reliable bulk phase correction. Therefore, these areas are shown in gray. By using only the data from the IS/OS junction for bulk motion correction, higher contrast is maintained in the RNFL phase difference pattern while the phase difference values at the IS/OS are about zero (Fig. 4(b)). This is expected as taking the entire depth range for bulk correction (Fig. 4(a)) compensates some of the motion of the RNFL itself. The IS/OS based bulk motion correction was used for the remainder of this study, assuming that the IS/OS phase information is stable over such period of time.

 

Fig. 4. Example of a phase difference scan at the systolic peak for S1 obtained with the M1 with the bulk correction using the entire A-scan depth (a) and only the pixels at the IS/OS junction (b). The scans correspond to the systolic peak marked in red in Fig. 3 (c). Large phase changes in the entire scan can be seen, with greater contrast in (b). For comparison, a phase difference scan between adjacent B-scans at the systolic peak is given in (c) (corresponds to one frame of M2), where only a slight increase in the values is seen around the vessels. Eye saccades result in de-correlated scans as shown in (d), corresponding to the saccade marked in black in Fig. 3 (b) and (c).

Download Full Size | PDF

Despite using the retinal tracker during acquisition and registration in the post-processing steps, eye saccades may be seen in some sequences. Such frames can easily be recognized and isolated as the phase correlation is lost and no pattern can be extracted from the phase difference scans (only the IS/OS green band is recognizable) as shown in Fig. 4(d).

2.3.3 Evaluation of the cardiac cycles sequences

After demonstration of the tissue motion in the individual consecutive phase difference scans, we aim to analyze and quantify the difference in RNF tissue deformation near veins and arteries over a full cardiac cycle. Such evaluation was done using RNFL deformation maps calculated in the following way for M1:

Step 1. The individual phase difference scans were flattened after segmentation of the upper boundary of the RNFL on the averaged intensity image.

Step 2. For each flattened phase difference tomogram, a band of five pixels along an A-scan (pixels i-2, i-1, i, i+1, i+2) at a given RNFL depth i was extracted and the data of the five pixels were averaged in complex space, similarly as done using the VM vector in OCE [56]. The data were additionally smoothened by averaging over the neighboring pixels (j-1, j+1) at a given A-scan j and over two consecutive bands in time.

Step 3. A wrapped temporal phase evolution OCT map (wrapped TPE map), i.e. an M-scan of phase evolution as a function of position along B-scan trace vs. time, of the smoothened phase difference corresponding to depth position i in the RNFL tissue, was calculated and plotted, covering one or several cardiac cycles.

Step 4. The corresponding deformation map (unwrapped TPE map) was obtained after unwrapping of the data, using the MATLAB unwrap function. The unwrapped map was smoothened by a Gaussian filter and converted to tissue deformation in µm using the relation below:

Step 5. Representative plots of the RNF tissue deformation were finally extracted close to each artery or vein by selecting a reference profile (tissue deformation over time for a given A-scan position) at the largest deformation and taking the mean of the 5 most correlated adjacent profiles to that reference profile.

For M2, similar steps were performed for the phase difference data between adjacent B-scans. The final deformation map is retrieved by temporal integration instead of unwrapping.

3. Results

3.1 Visualization of RNFL deformation

An example of a RNFL TPE map over three cardiac cycles, i.e. about 3.5 s, for S1 is given in Fig. 5 (corresponding to the intensity sequence of Fig. 3(b)). The maps were obtained at the RNFL band marked in red in (a) using M1. Large phase changes, including several phase wraps, are observed in the wrapped TPE map (b), indicating axial tissue motions beyond half a wavelength. In this map, the first frame of the first cardiac cycle is taken as reference for the calculation of the phase difference images to each of the subsequent scans. The wrapping pattern repeats itself over the three cardiac cycles. The TPE map given in (c) represents the tissue deformation after unwrapping and smoothing of the data displayed in (b). In the image displaying the phase difference between consecutive B-scans (d), changes in the RNFL are only seen close to the vessels during the systolic phase for each of the three heartbeat cycles as this data is related to the velocity of tissue. An integration of the map (d) along the temporal axis provides a comparable deformation map as given in (c). A video with the sequences of the individual intensity, phase difference scans (using M1) and deformation scans (using M2) over one cardiac cycle is provided in the supplementary material (Visualization 1). From a visual point of view, the changes in tissue expansion from frame to frame are better emphasized in the phase difference scans (using M1) while the depth extension of the deformation is better seen on the integrated deformation scans (using M2).

 

Fig. 5. Tissue motion analysis of RNFL over time. Wrapped (b) and unwrapped (c) RNFL TPE map (M-scan of phase/deformation evolution as a function of position along B-scan trace vs. time) at the location marked by red dotted line in (a) over three cardiac cycles for S1 are obtained using M1. The wrapping pattern repeats itself over the cardiac cycles in (b) corresponding to the deformations seen in (c). A similar TPE map as given in (c) is obtained after integration of the map (d) as done in M2 (not shown here). (e) Standard Doppler analysis extracted from artery A2 (A-scan 360) over the same temporal range. Colorbar: -π to π for (a), (b) and (d), <-0.25 µm to > 0.75 µm for (c).

Download Full Size | PDF

The acquired trace (B-scan) of Fig. 5 is crossing two arteries at A-scan locations 70 and 350 and a vein at location 200. While expansion starts in the RNF tissue around the arteries a few frames after the reference frame, a slight shrinkage is first seen on the vein pattern before expansion (marked by a white arrow in Fig. 5(b)). This may be explained by the fact that the tissue did not reach its fully contracted state at the time the reference frame was taken, i.e., there is a lag between arterial and venous pulsations. The behavior of the tissue deformation is further detailed in the next section 3.2, comparing the results for several subjects. An interesting pattern can be observed between the vein and the larger artery, i.e. between the A-scan positions 250 and 340. This pattern could be potentially caused by interference of mechanical tissue waves. The wrapped TPE and intensity maps, as a function of depth beneath the RNFL surface, recorded over one cardiac cycle, are presented in supplementary material (Visualization 2).

Figure 6 presents the RNFL deformation of another healthy subject over nine heartbeat cycles, i.e. almost 8 s, with the reference frame taken at the first scan of the first cycle. This example demonstrates the good phase stability over the acquisition sequence, which is ensured by the retinal tracker during acquisition. Despite some de-correlated scans, as seen for example at the beginning of the first cycle and probably due to an eye saccade, the phase stability is clearly retrieved in the next scans. This trace crosses both an artery and a vein, as observed in the center part of the averaged intensity B-scan. An additional small vessel can also be distinguished at the edge of the trace (at ∼ A-scan no. 6), causing additional RNF tissue deformation. A similar repetitive wrapping pattern under pulsation as in Fig. 5(b) is obtained (Fig. 6(b)). Less wrapping can however be observed, probably due to the smaller size of these vessels. The center artery and vein of Fig. 6 have diameters of about 67 µm and 58 µm, respectively, while the artery and the vein of Fig. 5 have diameters of about 104 µm and 145 µm, respectively. The amplitude of the deformation also seems to vary depending on the heartbeat cycle. In Fig. 6, the tissue returns to its initial state after each of the two first cycles. However, larger variations in the deformation pattern can be observed between the third and the seventh cycles, probably due to higher perfusion. The tissue seems to expand with the beginning of the next systolic phase before having time to first reach its initial state again. This phenomenon is clearly seen in the zoomed TPE maps on the center artery given in (c) and (d). The deformation map (d) was obtained using M2.

 

Fig. 6. RNFL TPE map recorded over an extended time of ∼ 8 s (b) at the top of the RNFL of S2 at position marked by red-dashed line in (a). Good phase stability is observed over most of the recording period, with exception of few de-correlated scans such as the ones marked by arrows. The reference scan was selected at the beginning of the first cardiac cycle for the entire sequence of nine cycles. A zoomed TPE map centered on the artery A is given in (c) with its corresponding deformation map in (d) obtained using M2. The variations in the extent of the tissue deformation over the different cycles is especially visible on the deformation map. Colorbar: -π to π for (b) and (c), <-0.25 µm to > 0.75 µm for (d).

Download Full Size | PDF

3.2 Investigation of deformation properties

The RNF tissue deformation near arteries and veins under pulsation are compared for the three healthy subjects of Fig. 2 (traces marked in red dotted lines on the LSLO images). For each subject, the three best-quality heartbeat cycles obtained in the 10 s recording time (higher coefficient of correlation over the sequence duration and lower amount of eye saccades) were considered in this analysis. The reference frame was taken at the beginning of each cardiac cycle, considering the three repetitions separately.

Figure 7 gives examples of the wrapped temporal phase evolution maps and unwrapped deformation maps of one cycle for each volunteer using M1. Despite the difference in maximum deformation observed, similar patterns are seen for all subjects. The expansion in the tissue surrounding the arteries is starting earlier during the systolic phase as compared to the expansion surrounding the veins. The plots given in Figs. 7(i), (j), and (k) confirm this behavior. A double hump pattern can be recognized for the tissue expansion near the arteries, probably corresponding to the peak of the systolic phase and the dicrotic notch. Only a single, broad expansion peak is observed around the veins, reaching significantly higher deformation as compared to the arteries for S1 and S3.

 

Fig. 7. Comparison of the RNF tissue deformation over one cardiac cycle for the three study subjects evaluated with the M1. For S1, S2 and S3 respectively: (a), (b) and (c) are the RNFL wrapped TPE maps at the top boundary of the RNFL, (d), (e) and (f) are the corresponding RNFL TPE deformation maps after unwrapping, (i), (j) and (k) are the deformation profiles extracted from (d), (e), and (f), respectively, close to each artery and vein. Similar tissue deformation patterns are observed in the three subjects. The expansion of the veins is delayed compared to the one of the arteries.

Download Full Size | PDF

The results obtained in the three repetitive cycles are listed in Table 1 providing vessel diameters and RNFL tissue deformation parameters relative to the reference frame at beginning of the heartbeat cycle. A maximum expansion of 1.58 ± 0.10 µm was obtained at the vein of S1. The average precision for a single pixel (standard deviation) of tissue deformation measurements is ∼ 30 nm. It is obtained from the analysis of the noise values in an area without pulsation on the wrapped TPE scan of S1. Also temporal measurements (such as the time at which this maximum is reached within a cardiac cycle) shows a good repeatability. From Table 1, we calculated the temporal delays between tissue responses around arteries and veins. Delays of 0.35 s, 0.15 s and 0.27 s (corresponding to about 29%, 19% and 30% of the total time of the cycle) were obtained between the maximum deformation reached at the arteries and the one reached at the veins for S1, S2 and S3, respectively. Delays between the beginning of the tissue expansion at the arteries and at the veins were found to be 0.05 s, 0.04 s and 0.05 s for S1, S2 and S3, respectively. Similarly, delays of 0.08 s, 0.04 s and 0.08 s are calculated for the maximum of the gradient during the tissue expansion (maximum tissue expansion speed (du/dt)_max).

Table 1. Vessel diameters and RNFL tissue deformation parameters (Δu and times relative to the reference frame) in the vicinity of vessels observed in three healthy volunteers. Mean values ± standard deviation for three cardiac cycles. Right before the arriving of the pulse wave, the tissue reaches its maximum contraction Δu_min. While the pulse wave propagates through the tissue, the maximum tissue expansion speed (du/dt)_max is then reached. The tissue expands up to a maximum value Δu_max before it starts retracting. NC = not calculable.

View Table

4. Discussion

We demonstrated the use of TPE OCT to analyze RNF tissue deformation induced by retinal vessels pulsating with the heartbeat. Using a custom-built SD-OCT system with an integrated retinal tracker, we acquired phase stable measurements over several heartbeat cycles in different healthy subjects. After processing of the TPE data, a reproducible phase pattern was observed in the vicinity of the retinal vessels over the course of the different cardiac cycles. From the unwrapped deformation maps, we could quantify and compare the tissue deformation near arteries and veins with an average precision for a single pixel of about 30 nm. The theoretical value, calculated as introduced in [57], is about 15 nm (based on an average SNR of 10 dB computed from the ratio of the average signal in the RNFL and the average noise).

Probing the biomechanical properties of ocular tissues has been of increasing interest, especially with the development of OCE and its capability to access tissue elastic modulus. Non-invasive OCE measurements are however still challenging in the eye and have essentially been limited to the cornea in human subjects in-vivo [24]. Compared to standard OCE which uses an external pulse source, our method relies on the natural cardiac pulsation to investigate tissue deformation, similarly to recent work reported on passive elastography in the cornea in humans [37] and in animal models [36]. While such method does not provide information on the elastic modulus (the force exerted on the tissue is unknown), it gives interesting insights on retinal tissue motion over the range of several seconds. In a further step, strain estimation could be performed based either on the LSM [18], as in heartbeat OCT [36], or the VM [19] and optimized VM [20,21].

Another PhS-OCT method that provides access to motion on a larger time scale was presented by Spahr et al., using a FF-SS-OCT system [38]. This method uses a very fast area camera to record retinal OCT volumes at a volume rate of 2000 volumes/s. From these volumes, 2D en-face images corresponding to two retinal layers (RPE and RNFL) are segmented, and the phase difference between these layers is obtained. From the temporal evolution of these phase differences, tissue deformation over time is derived by temporal integration of phase difference changes. An advantage of this method is the access to tissue deformation over a 2D en face area. The technique is, however, more sensitive to multiple scattered photons and artifacts caused by fast motion within vessels may deteriorate the image quality. Using FF-OCT, the choice of the imaging parameters is also less flexible compared to standard scanning OCT systems, which can be a limiting factor for various measurements. Additionally, scanning systems are already commonly used in clinical applications, while clinical translation of FF systems is still pending.

In order to retrieve motion over the range of seconds, high phase stability is essential during acquisition. Despite the moderate imaging speed of our SD-OCT system (A-scan rate of 70 kHz) compared to recent ultra-fast OCT systems [43], phase stability is ensured by the sub-speckle precision of the integrated retinal tracker [46,52]. In this study, we demonstrated phase stable measurements over a sequence of nine cardiac cycles, corresponding to about 8 s. During eye saccades, the phase correlation may be lost over a few scans but is retrieved in the subsequent scans. While we only acquired sequences of about 10 s (1000 B-scans in total), the phase stability could potentially be preserved over longer time, after adaptation of the acquisition software. In addition to the phase stability of the system, compliance of the imaged subject is essential for such long measurements. Over longer timeframes, rotation and other misalignments during acquisition may deteriorate the phase stability. The method was so far only applied to a few healthy volunteers, with good fixation skills and selecting cross-section traces of about 2 mm length. Imaging of a larger set of healthy subjects and/or patients with reduced fixation skills may require further adaptation of the system. For example, the robustness of the imaging method could be improved by implementing rotation correction on the retinal tracker (the current version only correcting translation).

After the recording of phase stable measurements, an essential post-processing step for the purpose of our analysis is the bulk motion correction. In standard Doppler OCT and OCTA, various methods have been used, as shown for example in [58,59]. In our analysis we tested two different methods: correcting the phase difference scans by the mean phase difference value over the entire depth at each A-scan position or by only considering the mean of the pixels corresponding to the IS/OS junction. We observed that the second method provides higher contrast in the TPE scans, as the first method is incorrectly compensating some of the RNF tissue motion. We chose to use the bulk motion correction taking only the IS/OS data, under the assumption that no or low motion will happen in this layer in the timeframe of our analysis (a few seconds) as compared to the larger RNFL motion induced by vessel pulsation. All quantitative measurements were performed on single and separately analyzed heartbeat cycles taking the reference frame at the beginning of each cycle. For a quantitative evaluation over longer time, for example taking the reference frame at the beginning of a first cycle and analyzing several subsequent cycles as in Fig. 6, the impact of slower phase changes at the photoreceptor layer will need to be further investigated.

Another important step of the data processing is the choice of the reference frame for the TPE analysis. As our study focuses on motion induced by vessel pulsation, we chose to take the reference scan at the beginning of a cardiac cycle. The choice of the reference frame will influence both, the computation of each phase difference tomogram using M1 and the temporal integration using M2. In this work, we determined the reference frame from the standard Doppler OCT data derived at the artery location. However, in the case of a Doppler angle close to 90° or when imaging smaller arteries, the Doppler signal may be noisier or not accessible. In these cases (such as S3 in this study), the beginning of the heartbeat cycle was determined by manually shifting the reference frame until obtaining a phase pattern similar to those usually observed near the retinal vessels (cf., e.g., Fig. 5(b)). Four other healthy volunteers were imaged additionally to the three study subjects, as we tested different acquisition parameters, and similar phase patterns at the arteries and veins have been observed in each case.

To further evaluate the influence of the reference scan on the results and to improve the interpretation of the data, we investigated the use of different reference frames for the B-scan trace of S1 as shown in Fig. 8. Additionally to the usual wrapped TPE map (c) with the reference frame at the beginning of the systolic phase, we reprocessed the data after choosing different reference frames: during the previous diastolic phase (0.09 s before the start of the systolic phase) (b), at half of the systolic peak (d) and at the systolic peak (e). When setting the reference frame at the end of the previous diastolic phase, we can see a larger tissue shrinkage with several phase wraps near the vein as well as a slight shrinkage near the arteries before expansion. When setting the reference frame at the beginning of the systolic phase and then at half of the systolic peak, this shrinkage is not seen any more close to the arteries. In these cases, the expansion of the arteries started after 0.1 s with the reference frame at the beginning of the systolic phase and right away with the reference frame at half of the systolic peak. Similarly, only a reduced tissue retraction is seen close to the vein before expansion on these images. Taking the reference frame at the systolic peak, the tissue expands directly near the vein while the usual double hump wrapping pattern is lost near the arteries. Only a slight tissue expansion can be distinguished before the RNF tissue shrinks again at the top and bottom arteries. Some more distorted phase wrapping patterns can still be seen near the middle artery, possibly due to the pattern of tissue wave interference between this artery and the vein. The quantitative analysis of tissue deformation is thus greatly dependent on the choice of the reference frame.

 

Fig. 8. Influence of shifting the reference frame on the measured phase difference. The TPE maps of S1 over the time frame of one cardiac cycle were calculated setting the reference frame at (b) the end of the previous diastolic phase, (c) the beginning of the systolic phase, (d) half of the systolic peak, (e) the systolic peak. The maps (b)-(d) are plotted until the end of the cardiac cycle and the map (e) until the systolic peak of the next cardiac cycle. Time points of the selected reference frames are marked in the standard Doppler OCT plot (derived at the position of the artery at A-scan 360) in (a).

Download Full Size | PDF

Data evaluation method M1 seems especially useful for analyzing subtle details of tissue deformation, visualized by phase difference contours that immediately show patterns of sub-wavelength deformation. To quantitatively analyze larger tissue deformation (larger than half the wavelength of the OCT sampling beam), phase unwrapping is used to convert the phase data into absolute distances. We used here the built-in MATLAB unwrapping function. The RNF tissue deformation measured in this study showed peak values up to about 1.6 µm, which is in agreement with literature data [38]. In the case of noisier data (as it was the case for some of the data of the additionally imaged subjects), the built-in algorithm failed to correctly unwrap the phase difference information of M1. The development of a custom algorithm might be beneficial that is less sensitive to the varying quality of datasets and that excludes unwrapping errors, for example caused by eye saccades. As only small displacements are induced between consecutive B-scans, unwrapping and therefore unwrapping errors are avoided using M2, which seems therefore especially suited for quantification of larger tissue deformation, as, e.g., for visualization of the depth extension of tissue deformation as shown in Visualization 1. In the case of higher rate motions, as for example seen in OCE [60], an additional unwrapping step may however be required. Eye saccades may nonetheless affect the performance of the temporal integration in M2. Especially in areas where only little motion occurs, inaccuracies in the phase difference values may have a stronger influence on the cumulated data. Optimization of the processing steps will permit a larger quantitative analysis and be particularly important for future work in order to evaluate patient data, which may be of lower quality. Implementation of the extended Knox-Thompson algorithm for the computation of the phase difference data, may additionally improve the overall quality of the processed data [50].

Our study showed that TPE OCT is able to provide measurement of tissue deformation with resolution in the 10 nm range over extended times. Possible applications comprise, e.g., the investigation of pulse wave propagation, the temporal analysis of tissue deformation, as well as time lags between motions induced by arterial and venous pulsations, under normal physiologic conditions, but also in various ocular disorders that are known to be associated with altered blood circulation, like glaucoma, diabetic retinopathy, or retinal vein occlusions. Moreover, our method might also be useful for analyzing cardiovascular diseases, e.g. by looking for changes in tissue deformation patterns near vessels that are affected by various types of plaques. On a longer term (and more speculative), by using the retina as a “window into the brain”, our method might be useful to identify patients at risk of brain damage from vascular complications. Other possible applications could comprise, e.g., the analysis of perfusion changes under visual stimulus, which could complement other recent optophysiologic and opto-retinographic studies [61–63] and thereby support various physiologic studies in the field of vision research.

5. Conclusion

In this work, we investigated the TPE method for measuring deformations in retinal tissue over the range of seconds. We evaluated this deformation in the RNFL in the vicinity of retinal vessels under cardiac pulsation in healthy human subjects. The integrated retinal tracker of the OCT system used for imaging, providing sub-speckle precision, ensured phase stability during acquisition. Reproducible results with a precision of about 30 nm were obtained and differences between deformation induced by veins and arteries of different caliber could be evaluated. In future work, TPE analysis of tissue deformation may provide valuable new information not only for the monitoring of ocular disorders but also other diseases affecting the blood circulation such as cardiovascular ones.