Quantification of three-dimensional dynamics of intercellular geometry under mechanical loading using a weighted directional adaptive-threshold method

Capturing and quantifying dynamic changes in three-dimensional cellular geometries on fast time scales is a challenge because of mechanical limitations of imaging systems as well as of the inherent tradeoffs between temporal resolution and image quality. We have combined a custom high-speed two-photon microscopy approach with a novel image segmentation method, the weighted directional adaptive-threshold (WDAT), to quantify the dimensions of intercellular spaces of cells under compressive stress on timescales previously inaccessible. The adaptation of a high-speed two-photon microscope addressed the need to capture events occurring on short timescales, while the WDAT method was developed to address artifacts of standard intensity-based analysis methods when applied to this system. Our novel approach is demonstrated by the enhanced temporal analysis of the three-dimensional cellular and extracellular deformations that accompany compressive loading of airway epithelial cells. of the dynamics of this process. Thus the validity of the proposed mechanotransduction mechanism remains uncertain. Here we developed a means to capture and quantify the dynamic spatial changes in extracellular geometry of cultured cells that accompany compressive loading.


Introduction
Cells and tissues frequently experience changes in the mechanical state of their environment. Elucidating the mechanotransduction processes that convert forces into downstream signaling responses requires detailed knowledge of the temporal and spatial deformations that accompany mechanical loading. For example, airway constriction in situ and compressive stresses applied in vitro both activate the EGF receptor in airway epithelial cells [1,2]. In vitro imaging studies indicate that compressive stress predominantly affects the extracellular space separating neighboring airway epithelial cells, focusing attention on mechanotransduction occurring via autocrine EGFR signaling in a shrinking interstitial space [1]. However, imaging limitations have, until now, prevented a detailed examination of the dynamics of this process. Thus the validity of the proposed mechanotransduction mechanism remains uncertain. Here we developed a means to capture and quantify the dynamic spatial changes in extracellular geometry of cultured cells that accompany compressive loading.
Specifically, we modified a custom-built high-speed two-photon microscope [3,4] and developed a novel image analysis procedure that allows robust image segmentation, overcoming issues that occur when applying standard intensity-based segmentation algorithms to images with inhomogeneous intensity distributions, an inherent problem for many microscopic imaging techniques. These methods can be applied to a wide range of biological scenarios in which cells, their constituents, and their surrounding extracellular spaces undergo dynamic changes due to intrinsic or extrinsic factors.

Cell Culture
Normal human bronchial epithelial (NHBE) cells were obtained from Lonza Biosciences and cultured at an air-liquid interface as previously described [5]. In brief, cells were expanded on tissue culture-treated plastic (5% CO 2 , 37°C) in bronchial epithelial growth medium (BEGM, Clonetics) supplemented with bovine serum albumin (1.5μg/ml) and retinoic acid (50nM). Passage 3 cells were plated onto uncoated nucleopore membranes (0.4μm pore size, Transwell Clear, Costar) and grown to confluence with a 1:1 mixture of BEGM and Dulbecco's modified Eagle medium (DMEM, Life Technologies) applied both apically and basally. Once cells reached confluence they were grown at an air-liquid interface (basal feeding only) for 14 days.

Instrumentation: High-Speed Two-Photon Microscope
To study deformation dynamics of cells and tissues under loading, an imaging modality with three-dimensional resolution and high frame rate was required. For this purpose we modified a previously described [3], custom-built high-speed two-photon microscope (HSTPM). The imaging speed of the microscope was enhanced by incorporating a polygonal mirror scanner composed of 50 aluminum-coated facets rotating at 3750 rpm to provide a line scan speed of 320μs/line in the x-axis (Lincoln Laser, Phoenix, AZ) ( Fig. 1). A slower galvanometer-driven scanner (Cambridge Technology, Cambridge, MA) deflected the line-scanning beam along the y-axis. To utilize the dead-time required for the flyback of the y-scan mirror, the image was also acquired during the reverse sweep. The specimen was raster scanned with a single excitation focus and photons were collected using high-sensitivity photomultiplier tubes (PMT, R3896, Hamamatsu, Bridgewater, NJ). The excitation light source was a femtosecond Ti:Sapphire laser (Mira 900, Coherent, Palo Alto, CA). To excite the fluorescein isothiocyanate (FITC) dextran conjugate (Sigma, Saint Louis, MO) an excitation wavelength of 780 nm was used. The optical setup was designed to interface with an upright microscope (Axioscope, Zeiss, Thornwood, NY). To perform a 3D volume scan, the objective was mounted on a computer-controlled piezoelectric objective translator (P-722.00, Physik Instrumente, Waldbronn, Germany). The axial translation of the objective along the z-axis yielded z-stacks of xy-plane images.
To image the cells and extracellular spaces in real-time an image transfer rate of 10 frames per second was used. For a given time-point, a typical z-stack involved about 30-40 individual images, meaning that the total acquisition time for one time-point was about 3-4 seconds. The voltage signal was measured using a 16-bit AD converter (NI 6251, National Instruments Corp., Austin, TX) in a common mode arrangement eliminating electric perturbations present in both the signal and ground lines. Images were displayed and saved in real-time by integrating the signal synchronized with the raster-scanning pattern using a custom-written LabView program. The field of view was approximately 80 μm × 95μm using a 63X long working distance, water immersion objective (Achroplan, 63X/0.9W, Zeiss). To avoid excitation saturation of chromophores and photodamage of the sample the average laser power incident upon the specimen was kept below 30mW.

Imaging Procedure
A custom stage was incorporated with the HSTPM to enable viewing of live NHBE cells during compressive loading. Individual transwell inserts were inverted (basal cell surface up, apical surface down) on the stage, and the objective of the upright microscope immersed in a drop of media on the basal surface of the transwell to maintain hydration (Fig 2). A 5% CO 2 controlled pressure source was connected via tubing to the otherwise sealed apical chamber of the transwell to allow specification of the transcellular pressure gradient (compressive stress). All experiments were performed at room temperature. To visualize the extracellular space surrounding cells, the basal media added when mounting the transwell insert on the microscope stage was loaded with a high molecular weight fluorescent FITC-dextran (either 70 or 150 kDa). These large dextran conjugates efficiently labeled the lateral intercellular space (LIS, Fig. 3), and exhibited minimal cellular uptake during the experiments relative to lower molecular weights dextrans (e.g. 10 or 20kDa, not shown). Prior to mechanical loading experiments, several minutes were allowed for diffusion and equilibration of the dextrans from the basal transwell surface into the LIS. In some experiments cells and the extracellular space were simultaneously labeled with CellTracker Green (Molecular Probes) and TexRed dextran respectively and visualized with standard two-channel epifluorescence microscopy.

Image Segmentation
In order to quantify morphological changes during loading, a segmentation algorithm is needed to process microscopy images and generate a simplified image representation [6]. The objective for our algorithm was to produce results comparable to those that would be obtained via manual segmentation, but in a systematic and automated manner. For our system, the simplified representation designates pixel as either LIS or intracellular space. Since the cells form a confluent monolayer, every pixel in the image corresponds to one of these categories. The number of pixels corresponding to the LIS region can be used to estimate the total LIS area. By repeating the segmentation procedure for all Z slices over multiple time points, the LIS volume dynamics can be quantified.
Initially, we applied two standard intensity-based thresholding methods to our systemintensity thresholding and adaptive thresholding. While alterations to the two-photon microscope enabled observations on shorter timescales, they also resulted in image properties which are challenging to analyze by standard intensity-based segmentation methods. The use of a fast-scanning polygonal mirror resulted in a non-uniform intensity distribution in the horizontal direction of the image due to partial clipping of the excitation beam by the edges of the polygonal mirror facets. This non-uniform illumination characteristic was a challenge for simple, standard intensity-based segmentation approaches.
Prior to the segmentation, all images were preprocessed using a 3×3 median filter -a standard de-noising method that removes shot-noise without blurring shape boundaries -a known problem for linear filters [7]. A previously applied method for identifying the LIS boundary is the application of a flat intensity threshold [1,8,9]. However, after alterations to enable fastsampling on the two-photon microscope, we found that thresholding was complicated by the non-uniform illumination characteristic of the images (Fig 4A-B). No choice of thresholding parameters could adequately capture the LIS geometry across images with heterogeneous intensity distributions (Fig. 5). Adaptive thresholding [8][9][10] algorithms are known to be more robust to variations in illumination intensity than flat threshold approaches. In the standard adaptive thresholding procedure, the local average of each pixel is subtracted prior to thresholding. The implicit assumption of the standard adaptive thresholding procedure is that the local average intensity is approximately proportional to the local background intensity. A relative threshold is applied to the intensity difference between the pixel of interest and its local average. While the use of adaptive thresholding significantly improved the quality of the segmentation, (Fig. 4C, Fig. 6), segmentation artifacts were present, particularly in LIS-dense regions of the image, which exhibit artificially high local averages. In such regions, pixels can be incorrectly classified as intracellular space because the pixel intensity is not sufficiently greater than the high local average intensity (Fig. 6A,B). Attempting to compensate for this problem by lowering the relative threshold results in "islands" of intracellular space incorrectly labeled as LIS (Fig. 6C). Thus, we sought a simple segmentation method that would be more robust than these standard intensity-based methods for our system.
To address the limitations of standard intensity-based segmentation methods, we developed the Weighted Directional Adaptive-Threshold (WDAT) method. Like adaptive thresholding, WDAT utilizes local averages to estimate local changes in background intensity. WDAT differs from standard adaptive thresholding in that it utilizes a set of directional averaging filters to calculate multiple weighted, local averages (Fig. 7). The use of multiple directional averages is particularly important near the LIS-dense regions. Combining multiple segmentation estimates based on directional averaging filters utilizes the line-like characteristic of the LIS to ensure proportional sampling of non-LIS background pixels in at least one of the directions.
The implementation details of the WDAT method are as follows. Raw images are preprocessed using a 3×3 median filter to remove shot noise (Fig. 7, step (1)). Local averages are weighted using a normalized Gaussian function centered on the pixel of interest. Specifically, for each pixel, estimates of the local average are taken by averaging along one of four directionsvertical, horizontal, and the two diagonals (Fig. 7, step (2)): (1) Where , correspond to the local intensity averages at (x,y) in the vertical, horizontal, and two diagonal directions, respectively. I(x, y) is the image intensity at (x,y) and Q is a fixed normalizing factor. N and σ are parameters that correspond to the size of the local averaging window, and the standard deviation of the Gaussian used for the weighting function, respectively. By applying a smooth weighting function rather than a flat average, the local intensity average estimates vary smoothly over the image. A pixel is classified as LIS if its intensity is sufficiently greater (as determined by a relative threshold) than the local average calculated by any of the four directional averaging filters (Fig. 7, step (4)): (2) Where LIS(x,y) is 1 if the pixel at (x,y) that corresponds to LIS and 0 if (x,y) corresponds to intracellular space. T is the relative threshold intensity parameter.
Finally, segmentations were processed with heuristics designed to remove regions uncharacteristic of a LIS space (Fig 7. step (5)). Spurious, single-pixel regions in the segmentation were eliminated with a 3×3 median filter. Since LIS spaces are contiguous, regions in the segmentation with a small number of pixels are removed from the segmentation. Small, disconnected region in the segmentation most often correspond to non-LIS intracellular regions where the local average intensity is low. The parameters of the WDAT method are the local averaging filter size (N), the width of Gaussian weighting function (σ), the relative threshold intensity (T), and the minimum LIS region size (A). WDAT parameters were calibrated using a preliminary set of images and then applied to all subsequent data analysis. To improve the robustness of the relative threshold intensity (T) across different data sets, it is specified in terms of the intensity standard deviation of the data, rather than an absolute intensity. A local filter size of 15μm, Gaussian width of 4μm, minimum LIS region size of 20μm 2 , and relative threshold intensity of 0.3σ I (where σ I is the intensity standard deviation of the image stack corresponding to the first time point in a dataset), were found to yield a good segmentation for a calibration image. These parameters values were fixed and applied to all subsequent data analysis.

Results and Discussion
To visualize the extracellular space (LIS) of cultured normal human bronchial epithelial cells bounded apically by tight junctions, laterally by neighboring cells, and basally by a porous substrate, we introduced fluorescent FITC-dextran into the underlying basal media (Fig. 3). We acquired stacks of x-y images separated by 1 μm along the apico-basal axis using the highspeed two-photon microscope described above, spanning the entire depth of the LIS (about 10-20 μm). While conventional laser scanning microscope can achieve 1 frame/second resolution, it is insufficient to measure the change in LIS dimension that requires 3D imaging with about 30 planes axially that corresponds to an imaging rate of about 2 stack/min. Clearly, dynamics changes of LIS that has time scale on the order of a minute cannot be accurately resolved. The current system has temporal resolution of about 10 Hz and can image an image stack within 5 seconds allowing much more accurate measurement of the dynamics of LIS dimensional changes. The use of the fast-scanning polygonal mirror-based two-photon microscope allowed us to observe and quantify the three-dimensional dynamic collapse of the LIS under mechanical loading.
Initially, we applied standard intensity-based approaches (thresholding and adaptive thresholding) to systematically quantify LIS volume dynamics. However, alterations to the imaging system implemented to optimize image acquisition time also results in image properties that are problematic for these segmentation approaches (see "methods" for details). We therefore developed the weighted directional adaptive-threshold method (WDAT) to identify LIS regions in images with inhomogeneous intensity distributions. WDAT properly segmented the entire field and captured the LIS architecture, allowing for a more detailed and accurate systematic characterization of LIS dynamics than standard intensity-based segmentation algorithms (Fig. 4-6).
To examine the sensitivity of our image analysis methods to the type of labeling used, we applied WDAT to a two-channel standard fluorescence image (Fig. 8) of NHBE cells labeled intracellularly with CellTracker Green (green pseudocolor) and extracellularly with TexRed dextran (red pseudocolor). Segmentation of each channel was performed independently and for the intracellular label we applied the inverse algorithm to provide a filtered image comparable to that resulting from the extracellular label. Based on the agreement between the LIS architectures obtained using either intra-or extracellular labeling (Fig. 8), we performed subsequent experiments utilizing only extracellular fluorescent dextran as the LIS marker.
The enhanced high-speed two photon microscopy approach allowed us to image dynamic loading-induced changes in extracellular geometry on the timescale of seconds, a timescale that was previously inaccessible. After introducing 150kDa FITC-dextran into the LIS until an equilibrium concentration was reached, a constant physiological pressure gradient of 30cmH 2 O was applied continuously for 600 seconds [1]. A stack of 2D sections was obtained for the initial, pre-pressure state (time point 0) followed by acquisition of stacks for time points 20, 120, 300, and 600 seconds after the application of pressure (Fig. 9). More frequent sampling was possible, but was unnecessary given the kinetics of the deformations we observed. For clarity, in Fig. 9 we show only three matched raw consecutive images from stacks obtained at time points 0, 20, and 600s seconds and the corresponding WDAT analyzed images. Complete LIS image stacks for 0 and 600s are shown in Fig. 10. Counting the pixels segmented as LIS (i.e. all of the white pixels in the analyzed images) in a given stack and multiplying by the slice thickness (1 □m) yielded an estimate of the LIS volume for a given time point. Normalizing to the initial, pre-collapse (t = 0) volume allowed us to quantify the dynamic LIS collapse curve (Fig. 9). Time-matched control experiments showed little drift in LIS volume in the absence of external loading.
Previous measurements of signaling induced by compressive stress focused on the time frame of 5 to 30 minutes after the onset of compression, with peak signaling activity shown to occur at ∼30 minutes [1,2]. Whereas previous imaging of deformations under compressive loading were limited to several minute after application of pressure [1], here we show that the decrease in LIS volume occurs on the order of seconds to minutes (Fig. 9). These finding suggests that the deformations and key mechanotransduction events in this system happen on a time scale shorter than previously appreciated. Our results thus serve to focus future efforts on understanding the dynamics of mechanotransduction events that occur in synchrony with the measured deformations. These quantitative imaging results also provide an essential input for mathematical models [11] that relate changes in local geometry to autocrine signaling, helping to refine our understanding of the dynamics of extracellular mechanotransduction. Most importantly, the improved imaging and image analysis methods detailed here offer an enhanced capacity to capture and quantify dynamic, three dimensional cellular deformations under loading.

Conclusion
The quantification of the dynamic decrease of the extracellular geometry of mechanically loaded epithelial cells presented a challenge both in terms of instrumentation and subsequent image analysis. To overcome these challenges we combined a custom high-speed two-photon microscopy approach with a newly developed automated image segmentation method. The custom high speed two-photon microscope allowed us, for the first time, to acquire real-time 3-D images of the collapsing lateral intercellular space (LIS) separating neighboring epithelial cells. Our novel automated image processing method, termed the weighted directional adaptive-threshold (WDAT) method, enabled us to systematically analyze and quantify the dynamics of the LIS volume. This experimental setup represents a significant advance over previous imaging experiments with mechanically loaded airway epithelial cells [1], permitting measurements on shorter timescales, and is applicable to other systems for which quantification of cellular geometries on short time scales is of interest. Applying our imaging-analysis system to mechanically loaded human bronchial epithelial cells revealed that key mechanotransduction events occur on timescales shorter than previously accessible. The imaging-analysis system allows for additional studies aimed at further elucidating mechanotransduction events in this short timescale regime. Schematic of the high-speed two-photon microscope imaging system: Infrared excitation pulsed beam (optical path depicted in red) is scanned through a polygonal mirror (fast, horizontal) and a galvanometric mirror (slow, vertical) to generate raster scanning on a sample. The emitted fluorescence (optical path depicted in green) is collected by the objective lens and detected by PMT after spectral filtration by a dichroic mirror and a barrier filter. Schematic showing the interface of the two-photon imaging system with transwell cell culture insert. The HSTPM in interfaced to an upright microscope, hence the transwell is inverted on a custom stage. A controlled pressure source provides pneumatic pressure via a tube going through the plug on the apical side of the cells. The bottom of the transwell (area ∼0.33cm 2 ) is covered with 100μL of media containing fluorescent FITC-dextran, which diffuses into the LIS of the cells (bounded by tight junctions apically). For convenience, we define the z-axis as the baso-apical depth axis, and the x and y axes as the corresponding planar axes.    Detailed comparison of adaptive thresholding and WDAT. The raw intensity image of a ∼30×20μm region is shown on the left of each segmentation image for comparison. (a) Adaptive thresholding with a high relative threshold parameter results in patchy, overly-sparse LIS structures. (b) Adaptive thresholding with a moderate relative threshold is superior to standard thresholding, but the artificially high local average intensities at LIS intersections results in misclassified pixels at these regions. (c) Adaptive thresholding with a low threshold results in small patches intracellular regions being incorrectly classified as LIS space (d) The WDAT method addresses both of these LIS segmentation issues. Outline of the WDAT image analysis procedure. Images are preprocessed with a median filter to remove shot noise (step 1), local averages are calculated using weighted directional averaging filters (step 2), a threshold is applied relative to each local average (step 3), segmentation results are combined (step 4), and removal of regions uncharacteristic of LIS regions (small, disconnected regions) is performed on the combined segmentation (step 5).  Timecourse of changes in LIS dimensions. Analyzed and raw images for 3 consecutive 1micron sections for t=0 (pre-collapse) and 20, 600 seconds after application of a constant 30cmH 2 O transcellular pressure gradient; the white region is LIS, determined using WDAT algorithm. The corresponding LIS volume change plot is shown, normalized to the initial LIS volume, obtained by summing all of the white pixels (LIS) for a stack at a given time point and normalizing to the t=0 value. Full LIS image stacks at 0 and 600 seconds. Top row: WDAT analyzed images for consecutive 1micron sections (white is LIS). Middle row: original raw images. Bottom row: combined image (original and inverted analyzed image, black is LIS) to demonstrate the effectiveness of WDAT algorithm in identifying the LIS.