1. Introduction

The field of mechanobiology studies the role of material properties and cell forces in physiological and disease processes such as cancer, wound healing, and development. It has been shown that mechanical cues, such as substrate stiffness or local stresses, have bearing on cellular behavior in many contexts [1–6]. For example, it has been shown that extracellular matrix (ECM) stiffness can promote the onset and progression of cancer [7–10], and that metastatic cancer cells exert greater traction forces than their healthy counterparts [11]. As the field advances, mechanobiology will offer a more complete understanding of cell-ECM interactions, and may support the development of novel cell physics-based models and diagnostics or therapeutics [12, 13].

Traction force microscopy (TFM) is a method used to measure the forces exerted by cells on their environment [14–17]. This is achieved by measuring cell-induced deformations in a substrate, such as a micropillar array, hydrogel, or ECM, followed by reconstruction of the exerted forces by inverting a stress-strain constitutive relation using the substrate material properties in an analytical model or finite element method (FEM) solver. TFM has traditionally been performed on samples consisting of cells grown on 2D substrates. One area of interest is that of collective cellular behavior, where 2D TFM studies have revealed several collective cellular behaviors including wave patterns in cell migration fronts [18], space-filling behaviors [19], tendencies to minimize intercellular shear stresses [20], and jamming transitions [5, 21]. Recent work has led to the development of TFM in 3D cell culture, with confocal microscopy taking the lead as the modality of choice to enable imaging of 3D samples and deformations [15, 17, 22–26]. It remains an open area of research how the collective cell behaviors observed in 2D settings will translate within 3D environments.

Mechanical signals, and the responses they invoke, span a wide range of spatiotemporal scales—from individual molecules to entire tissues and organs, and from collective cell migration over hours, days, or weeks, to cellular mechanotransduction events on the order of seconds to minutes [3, 6, 27–29]. To understand how all of these systems interact, it is crucial that TFM and other methods in mechanobiology capture information over a large range of spatiotemporal scales. Moreover, as cellular behavior can strongly depend on the composition and dimensionality of surrounding substrates, the ability to observe cellular behavior within 3D models of the ECM is of high priority [3, 6, 27, 28]. High spatiotemporal coverage—the ability to capture large volumes at high resolution, with rapid acquisition that may be maintained for extended times—is therefore a key trait required for the advancement of 3D TFM research. To date, 3D TFM has addressed these needs primarily through the use of confocal microscopy. However, this imaging modality is limited to a depth range of only a few hundred micrometers, placing a limit on the volumetric field-of-view that may be observed. Furthermore, total imaging time and maximum temporal sampling rate are limited by photobleaching and phototoxicity [15, 17]. To circumvent the limitations imposed by confocal microscopy, while addressing the needs of TFM studies, we propose the use of optical coherence tomography as a novel platform for 3D TFM.

Optical coherence tomography (OCT) is an interferometric imaging modality that enables non-invasive, label-free imaging in scattering biological samples. OCT can achieve rapid, volumetric imaging with micrometer level resolutions and imaging depths on the order of hundreds of micrometers to a few millimeters [30, 31]. Optical coherence microscopy (OCM) is a variation of OCT methods, characterized by a high transverse resolution [32]. At resolutions approaching 1 µm, OCM can enable the capture of cellular-level features. Finally, as a volumetric extension of holographic imaging methods (which are otherwise typically considered as confined to optically thin samples [33]), OCT/OCM provides access to phase information of the backscattered optical field, which can be used to enhance tomographic image formation [34, 35] and displacement sensitivity [36]. These attributes make OCM a promising candidate to offer new imaging capabilities in mechanobiology research.

Prior work has established OCT as a viable imaging modality for the study of both cell migration dynamics and cell-material interactions. OCT has been used to image and track the 3D migration dynamics of cells embedded in scattering substrates, with temporal resolutions on minute timescales maintained over multi-hour experiments [37–39], and to analyze cell-material interactions of cells migrating on substrates with different micro-topographies [37]. Among the fastest variations of OCT/OCM imaging are full-field OCM, which has been used to capture millisecond-scale cellular dynamics [40, 41], and swept-source OCM, which can boast line scan rates in the multi-MHz range [42], enabling acquisition of volumes at speeds and sizes more often seen in light sheet microscopy [43, 44]. Recent work has shown that diffusion-sensitive OCT can be used to quantify pore size (fiber spacing) in the ECM, a physical feature which plays a crucial role in the behavior of invasive cancer [45]. Another class of methods, optical coherence elastography (OCE), enables the measurement of mechanical properties of tissue by measuring small displacements from artificially induced (as opposed to cell-induced) forces [36, 46, 47].

We propose traction force optical coherence microscopy (TF-OCM) as a new technique to image cell force – an important biophysical parameter studied in the field of cell mechanics. Previous cellular-resolution OCM research focused on directly imaging cell behavior, whereas TF-OCM seeks to infer cell forces based on precise imaging of (cell-induced) substrate deformations. We demonstrate that volumetric data acquired with OCM can be used to measure the cell-induced ECM deformations vital to TFM applications in fully 3D environments.

2. Methods

2.1 Cell culture and sample preparation

NIH-3T3 fibroblasts (ATCC, CRL-1658) were maintained in media consisting of a 100:10:1 solution of Dulbecco’s Modified Eagle Medium (Life Technologies), bovine calf serum (Life Technologies), and penicillin-streptomycin (Life Technologies). 3D cell cultures were prepared for imaging by first suspending fibroblasts in chilled (4°C) media at a density of 200,000 cells/mL. To add scattering contrast for use in ECM deformation tracking, 0.5 µm diameter polystyrene beads in aqueous solution (Sigma-Aldrich) were added to the cell suspension to achieve a concentration of approximately 1 × 10⁹ beads/mL. This suspension was then combined in a ratio of 1:1 mL/mL with Matrigel (Corning). The resulting mixture was deposited in 100 µL aliquots onto the surface of glass-bottomed petri dishes with a 10 mm diameter circular microwell (MatTek) and left to polymerize for 30 minutes in an incubator. Samples were then covered with culture media, left overnight in an incubator (~12 hours), and imaged the following day. In addition to these 3D cell cultures, imaging phantoms for displacement tracking validation were prepared using an identical protocol, excluding the addition of cells.

2.2 OCM system and data acquisition

Samples were imaged using a custom-built spectral-domain (SD)-OCM system. The illumination source was a Ti:Sapphire laser (Femtolasers, INTEGRAL Element) with a central wavelength of 800 nm and a full-width-at-half-maximum bandwidth of 160 nm. A-scans were acquired using a spectrometer (Wasatch Photonics) with a 300 nm full bandwidth, coupled to a 4096 pixel 12-bit CMOS line scan camera (Teledyne Dalsa, Piranha4). This yielded an axial sampling period of 0.75 µm per pixel (in aqueous media) after image reconstruction. The system exhibited a sensitivity of 85 dB at 350 μm optical path delay from the zero optical path delay position, and a sensitivity fall-off of −10 dB/mm of optical path length. We believe that the sensitivity was limited by the alignment of the line scan camera, which had square 10.56 μm pixels, rather than ‘tall’ rectangular pixels which provide superior tolerance to misalignment. The objective lens used in this study was an Olympus XLUMPlanFl 20 × /0.95 W ∞/0. Samples were imaged in an inverted configuration, with illumination and collection occurring through the bottom surface of the glass-bottomed petri dishes. In this configuration, the axial and transverse resolutions of the system were approximately 2.4 and 1.5 µm, respectively.

Volumes acquired under the time-lapse imaging protocol consisted of 1024 × 1024 A-scans, acquired at a line scan rate of 30 kHz with an exposure time of 12 µs. The dimensions (z × x × y) of each full reconstructed volume were 1536 × 420 × 420 µm in aqueous media (2048 × 1024 × 1024 voxels). However, due to constraints imposed by the limited depth-of-field of the objective lens, only a 225 × 420 × 420 µm (300 × 1024 × 1024 voxels) volume surrounding the imaged cell/s was reconstructed and saved for further analysis using a high-speed SD-OCM processing method for depth-selective reconstruction [48]. To achieve this, the imaged cells and the focal plane of the system were colocalized prior to image acquisition at a position of 1 mm below the zero optical path delay position, enabled by custom-built acquisition software with a real-time live display and computational adaptive optics enabled [49]. Using these volume dimensions enabled the use of a high transverse spatial sampling frequency (at 410 nm/pixel, to accommodate the high transverse resolution) while maintaining a field-of-view which captured a large region surrounding each imaged cell, all while maintaining a manageable data set size for analysis in MATLAB with limited computer memory (128 GB RAM on our system). Approximately 150 GB of raw SD-OCM data was acquired for each time-lapse imaging experiment.

2.3 Experimental design and time-lapse imaging protocol

Physiological temperature, humidity, and pH were maintained throughout imaging using an incubation stage (Okolab, UNO-PLUS). To demonstrate the viability of using OCM imaging to quantify the dynamics of 3D cell-induced ECM deformations, cells in the previously described 3D culture were exposed to either cytochalasin D (a well-characterized actin cytoskeleton destabilizing agent that inhibits cell contractility [50]) dissolved in dimethyl sulfoxide (DMSO), or pure DMSO as a control. The imaging and exposure protocol was as follows.

To identify and select a fibroblast cell for imaging, the sample was translated until a cell (which appears as a ~30-100 μm scattering structure embedded in the optically clear Matrigel) became visible in the live display. As we sought to image cells that were actively interacting with the ECM, the candidate cell was inspected for the presence of filopodia protruding into the ECM, and for a generally elongated morphology, as is typical for healthy NIH-3T3 cells in 3D culture [51]. If these features were not present (i.e., the cell lacked protrusions or had a roughly spherical morphology), the cell was rejected and the searching process continued. After identifying what appeared to be a suitable cell or cell pair, the sample and imaging field-of-view were fixed in place for the remainder of the experiment.

Each sample was imaged every 5 minutes for a total time of 90 minutes. Baseline cell activity (without reagent exposure) was recorded during the first 30 minutes. After 30 minutes, samples were exposed to the appropriate reagent (cytochalasin D solution or pure DMSO). Cytochalasin D was dissolved into DMSO at a concentration of 0.5 mM. The solution was then added to the culture media in which the sample was immersed, resulting in a final cytochalasin D concentration of 1 µM, and a final DMSO concentration of 0.2 vol%. In the control experiment, pure DMSO was added to the culture media to achieve a final concentration of 0.2 vol%. Cell activity was then monitored for an additional 60 minutes following the introduction of these compounds, yielding a final total of 19 time-points per sample.

2.4 Reconstruction of volumetric OCM data sets and image sequences for displacement tracking

Volumetric data sets were reconstructed in MATLAB R2014b using the standard operations of background spectrum subtraction, computational compensation of spectrometer nonlinearity and dispersion, and transformation from frequency domain to space domain [52]. The high resolution objective resulted in rapidly increasing defocus artifacts with distance from the focal plane. Therefore, defocus was compensated using computational adaptive optics (CAO) [35]. The technique used consisted of a plane-by-plane 2D deconvolution, which was performed in the spatial frequency domain using multiplication with a phase-only correction kernel [35, 49]. The kernel parameters were determined with the use of a peak intensity image metric. This metric was used to obtain an initial estimate of the proper defocus correction at each depth in a given volumetric data set. The final applied kernel parameters were obtained from a linear fit of the initial estimates, as the severity of the defocus phase aberration increases approximately linearly with distance from the focal plane. This fitting operation helped to reduce noise in the metric-driven correction kernel calibration. The underlying assumption is that phase aberrations vary slowly with depth, such that adjacent or nearby depths must undergo similar corrections. The linear relationship of the defocus phase aberration in depth is an assumption used in prior work, and is based upon axial beam propagation [49].

To study local average ECM deformations in the transverse and vertical dimensions of the sample, maximum intensity projections were obtained from the reconstructed volumes. This operation collapses the 3D displacement tracking problem into a 2D displacement tracking problem, which may be performed independently within cross-sectional slices of the imaged volume along both the transverse and vertical orientations. Projections in the transverse orientation were obtained from the full 225 × 420 × 420 µm volumes previously described (forming a 420 × 420 µm projection in the xy-plane). Vertically-oriented projections were obtained by projecting along a 40 µm thick vertical slice intersecting the imaged cell body/bodies (forming a 225 × 420 µm projection in the zx-plane).

Over the 90 minute imaging span, some bulk sample drift was present between acquisitions. This bulk drift between time-points was removed using 2D cross-correlation performed over the full field-of-view of the maximum intensity projection image sequences. This operation acts as a high-pass spatial filter on the computed cell-induced ECM displacement fields, and assumes that the imaged cell does not cause a bulk displacement of the substrate it is embedded in. This operation has practical benefits for computation time which will be discussed in Section 2.6. As the vertical slices exhibited intensities which rapidly decreased with distance from the focal plane, the mean intensity at each depth was normalized prior to any displacement tracking.

2.5 Manual tracking of individual particles

To obtain a basic picture of cell-induced ECM deformation over time, polystyrene beads embedded in the Matrigel substrate at selected regions of interest were identified and tracked manually in the transverse image sequences. For this proof-of-concept work, this manual tracking was a simple alternative to automating the identification of cells, their orientation, and their behavior. It also served as a simple means to assess the reliability of the automated displacement tracking discussed in Section 2.6. In principle, automating the identification of cell features in combination with our automatic tracking could eliminate this step, but was not explored in this work.

The manual tracking consisted of two steps. First, the approximate location of a chosen bead was manually identified at each time-point. Second, the bead centroid was computed using the image intensity. The centroid location was assumed to correspond to the location of the bead. For both the control and experimental cases, four beads were selected from four different regions of interest in the vicinity of the imaged cell body/bodies. The selection of each specific bead was arbitrary, as the manual tracking data are only meant to represent the general behavior of the neighborhood surrounding each bead. The regions of interest for bead selection included the neighborhood near the leading edge of the cell (in the general direction of migration, determined from time-lapse animations, such as Visualization 1 and Visualization 2), the neighborhood of the trailing edge of the cell (opposite the leading edge), the neighborhood immediately adjacent to the cell body (neither leading nor trailing, and less than 100 μm away), and a neighborhood far from the cell body (at least 100 μm) that does not lie in the leading or trailing paths. The first three neighborhoods are intended to demonstrate cell-induced ECM deformations, while the fourth is meant to be a region of little to no deformation activity. In the event that a suitable particle could not be found in a given region, a close alternative was selected instead. For example, the control sample shown in Fig. 3(a) contained two adjacent cells with opposite migration paths, so that no definite trailing edge could be identified. As a result, the trailing edge measurement was replaced with another adjacent edge measurement.

2.6 Automated 3D displacement tracking

ECM displacement maps were obtained in a fully automated fashion using a digital image correlation (DIC) based technique to process the registered maximum intensity projection time sequences described in Section 2.4 [17, 53–55]. For a given time series of 2D maximum intensity projections, a 2D grid of interrogation points was established to specify where displacement data would be computed. In the xy-planes, the interrogation point grid spacing was 4.1 µm (10 pixels). In the zx-planes, interrogation point grid spacings were 7.5 µm (10 pixels) and 4.1 µm (10 pixels) along the z and x directions, respectively. For each interrogation point, a reference image was windowed from the current time-point image, centered at the interrogation point. This reference image was 41.4 × 41.4 µm (101 × 101 pixels) in the xy-planes, and 75.8 × 41.4 µm (101 × 101 pixels) in the zx-planes. This reference image was cross-correlated with a deformed-state image, obtained from the next time-point image and centered at the same interrogation point. This deformed-state image was 33.2 × 33.2 µm (81 × 81 pixels) in the xy-planes, and 72.8 × 33.2 µm (97 × 97 pixels) in the zx-planes. It should be noted that the difference in size between the reference and deformed-state images only allows for the measurement of displacements that do not exceed ± 10 pixels in the xy-planes or ± 2 pixels in the zx-planes. While this limits the range of local displacements that may be detected between time points, it helps to speed up computation over the many cross-correlation operations required across a time-lapse data set. Any bulk drift present between time points would further reduce the maximum detectable local displacement, if the reference and deformed state image windows were not increased to compensate. This makes the image registration step described in Section 2.4 vital for efficient performance of displacement tracking in this work.

The cross-correlation of the reference and deformed-state images was upsampled by a factor of 11 via zero-padding in the Fourier domain to allow measurement of subpixel displacements. The coordinates of the peak of the upsampled cross-correlation were taken to provide the 2D translation between the reference and deformed-state images. This process was repeated for all interrogation points between all pairs of images adjacent in time. The output of this process is a local ECM displacement field at all interrogated locations in the image sequence. A median filter was applied to reduce noise artifacts in the computed displacement fields. The 2D filter used spanned 10 × 10 interrogation points. Using a method based on simulated sinusoidal deformations of an imaged phantom [22], the spatial resolution of the displacement tracking algorithm (that is, the full-width-at-half-maximum of the algorithm response to an idealized point displacement) was found to be approximately 35 µm and 51 µm in the transverse and axial directions, respectively.

To determine the algorithm’s sensitivity to small displacements, displacement tracking was performed on sequential images of a stationary Matrigel phantom (described in Section 2.1). This measurement yielded a displacement noise floor of approximately 110 nm and 60 nm between time points in the xy- and zx-planes, respectively. Results and discussion of more generalized sensitivity measurements using different correlation window and median filter sizes may be found in Fig. 1 and Section 3.1.

 

Fig. 1 Measurement noise floor of the DIC-based ECM displacement tracking algorithm (described in Section 2.6) in the (a) xy-plane and (b) zx-plane, respectively. The variable ‘correlation window side length’ refers to the side length of the windowed deformed-state image. The different curves represent the use of median filters of different sizes following the cross-correlation operation. Displacement noise floors corresponding to the final parameters used in the tracking of cell-induced displacements are denoted by black arrows.

Download Full Size | PDF

To obtain the cumulative displacements of the ECM over time, the computed displacements between time points must undergo a cumulative integration in time. However, the process described above computes displacements at interrogation points which are defined at fixed locations with respect to the lab coordinate frame. A cumulative summation in time will not accurately represent the displacement experienced by individual particles unless performed in the material coordinate frame. Let the vector $r$ denote a fixed position, defined with respect to the origin of the lab coordinate frame, and let $n$ denote the integer time step ranging from 0 to $N-1$, where $N$ is the total number of time-points acquired. In this work, a total of $N=19$ time steps were acquired for each sample, so $n$ is taken to hold integer values from 0 to 18. Now, we define two displacement fields. The first is $u_n\left(r\right)$, defined as the computed displacement undergone by a particle between time $n-1$ and time $n$, given that its position in the lab coordinate frame is $r$ at time $n-1$. The second is $U_{0,n}\left(r\right)$, defined as the cumulative displacement undergone by a particle between time 0 and time $n$, given that its position in the lab coordinate frame is $r$ at time 0. In this work, the cumulative particle displacement was computed using the relation

3. Results

3.1 Sensitivity of automated displacement tracking

Images of stationary Matrigel phantoms (discussed in Section 2.1) were used to test the displacement measurement noise floor of our fully automated displacement tracking algorithm under varying conditions, including both changes in cross-correlation window size and median filter size, in order to study how these parameters affect our algorithm’s performance. The results of this testing are summarized in Fig. 1. The correlation window side length values correspond to the size of the windowed deformed state image. (The final window side lengths used in this work were 81 pixels and 97 pixels in the transverse and vertical planes, respectively.) In all cases, the side length of the windowed reference-state was maintained at a constant 20 pixels greater than the deformed-state window side length. In this way, each case was able to detect displacements over a fixed range of ± 10 pixels. After the initial correlation-based displacement tracking was performed, median filters of different sizes were applied across the data obtained at each interrogation point. (The final median filtering used in this work took place over 10 × 10 interrogation points.) The noise floor of the automated displacement tracking algorithm was taken to be the standard deviation of the resulting computed displacement field. As stated previously, the displacement sensitivity achieved under the processing conditions used in the remainder of this work was 110 nm (~0.25 transverse pixels) and 60 nm (~0.1 axial pixels) in the xy- and zx-planes, respectively, comparable to that reported in prior DIC work in OCE literature [53].

There are two trends evident in these results. First, as the correlation window size increases, displacement tracking noise declines. This is likely due to the fact that larger window sizes capture larger numbers of scattering beads embedded in the sample substrate. The added structure and signal makes the cross-correlation more robust to noise in reconstructed OCM images. However, this reduction in noise comes at the cost of a larger window, which acts as spatial low-pass filter when measuring the displacement field. Second, as median filter size increases, displacement tracking noise declines. The median filter acts to suppress rapid changes across space in the measured displacement field, and so also exhibits low-pass filter-like behavior. Though both of these trends are related to the principles of low-pass filtering, it should be noted that both the cross-correlation-based tracking and the median filtering are nonlinear operations on the space domain displacement field signal. As a result, though the general trends discussed here may be useful when designing experiments and selecting processing parameters, the methods remain vulnerable to potentially unpredictable and inconsistent behaviors that may depend on the images being used.

Another feature that affects displacement measurement sensitivity is CAO (performed here as a digital refocusing across depth, as described in Section 2.4). To demonstrate the impact of CAO, out-of-focus en face planes (located 30 μm below the focus) were extracted from the volumetric images of the stationary phantoms described previously. Displacement tracking was performed on both raw and digitally refocused versions of these en face planes in a manner almost identical to that described in the Methods sections, with the only difference being that these en face planes were extracted directly from the image volume and not through a maximum intensity projection along depth. The resulting displacement noise floor in the raw and digitally refocused cases were 1270 nm and 340 nm, respectively. The use of CAO therefore improved sensitivity of the displacement tracking algorithm by nearly a factor of 4. We attribute this to the signal-to-noise ratio and resolution improvement that can accompany CAO [35]. The decrease in performance of these cases relative to the 110 nm transverse displacement noise floor previously reported is likely due to a reduction in both signal content and strength, as the en face planes were not obtained through maximum intensity projection, and therefore lack additional signal contributions (information) from nearby depths.

3.2 Three-dimensional ECM displacements of control vs. force-inhibited cells

The results of fully automated 3D displacement tracking are summarized in Fig. 2, Visualization 1, and Visualization 2. Figure 2 depicts ECM displacements at a single time-point, whereas the visualizations depict displacements across all time-points as animations. Though every sample underwent displacement tracking, only one control sample and one contractility-inhibited sample are fully illustrated in Fig. 2 and the visualizations. These samples will be discussed in detail in the following sections as representative results of our algorithm. The remaining imaged samples and their behaviors will be discussed at the end of Section 3.3.

 

Fig. 2 Automated tracking of 3D deformations induced by NIH-3T3 fibroblasts cultured in a Matrigel-derived ECM. These images represent the deformations accumulated over a 90 minute imaging time, with reagents introduced to the sample after the first 30 minutes of imaging. Cells were exposed to pure DMSO (a-c), or cytochalasin D solution (d-f). Each subfigure depicts displacements in the en face (upper panels) and vertical (lower panels) orientations. (a,d) Superposition of the initial (t = 0 minutes, red channel) and final (t = 90 minutes, green channel) states of the sample, obtained from the initial and final registered maximum intensity projection images described in Section 2.4. (b,e) Cumulative displacement magnitude of the extracellular matrix (in µm) from a given initial location. (c,f) Cumulative displacement field depicting the direction and relative magnitude of ECM displacement (with arrow lengths exaggerated for visibility), superimposed on the initial maximum intensity projection images. Scale bars = 50 µm. Refer to the text for a discussion of the arrows in (a), (d), and (e).

Download Full Size | PDF

The control sample (visible in Figs. 2(a)-2(c) and Visualization 1) contains two adjacent fibroblasts (indicated by white arrows in Fig. 2(a)). Though the cells are distinct during the final 10 minutes of the experiment, it was unclear whether the cells were fully independent, or were a pair of daughter cells in an ongoing cell division. These cells appear to pull on the ECM in opposing directions toward the top and bottom of the en face image frame, respectively. In the path of the upper cell lies a thin scattering spindle (indicated by the yellow arrow in Fig. 2(a)), which may be a protein fiber, dead cell, or other structure. As the upward-migrating cell reaches this structure, ECM deformation in the vicinity of this object increases and extends over a large area (see Visualization 1). In fact, Figs. 2(a)-2(c) show that both cells influence ECM deformation as far as 100 µm or more away from their bodies, consistent with prior literature [56]. The areas of greatest displacement appear along the leading edges of cells, where cellular protrusions can be seen dynamically changing structure and interacting with the ECM in Visualization 1. It should be noted that the calculated displacement magnitudes and fields (Figs. 2(b), 2(c), 2(e), 2(f)) may have greater error near the edge of the field-of-view or in the vicinity of objects with variable scattering structure. This may be due to boundary effects of CAO, cross-correlation, and median filtering, or due to speckle decorrelation resulting from fast sub-resolution changes in object structure. For example, see the region indicated by the white arrow in Fig. 2(e). This region of apparent high displacement overlaps with a round cell intersecting the edge of the field-of-view. As can be seen in Visualization 2, the large local displacements reported by the tracking algorithm do not appear to be truly occurring in the neighborhood of the cell. This may be due to the constantly changing speckle pattern formed by the internal structure of the round cell. This speckle pattern, which changes quickly relative to the time-lapse sampling rate, could cause the assumptions underlying DIC (in this case, assuming small changes in scattering structure between time points) to break down, resulting in incorrect reported displacements.

As shown by the images and computed displacement fields, the deformations are not confined to a single plane, but extend in all three dimensions radiating away from the cellular protrusions into the ECM, forming a pattern reminiscent of a dipole field [13]. The physiological significance of such long-range, 3D deformations is an open area of investigation. As cells have been shown to respond to mechanical stimuli, the study of these deformation fields may shed further light on mechanical communication and collective behaviors among neighborhoods of cells [57].

The sample exposed to cytochalasin D (visible in Figs. 2(d)-2(f) and Visualization 2) exhibits differing temporal behavior compared to the control case. Early in the experiment, the cell (indicated by the white arrow in Fig. 2(d)) can be seen deforming the ECM in a similar fashion to the control sample discussed previously. It should be noted that the baseline rate of cell-induced deformation (apparent in Visualization 1 and Visualization 2) differs between the cytochalasin D and DMSO experiments. We attribute this to natural variability in cellular activity. Such variability is apparent both between and within the control and contractility-inhibited experimental groups, as is discussed in Section 3.3 and can be seen in Figs. 4 and 5. Future studies that seek to identify novel biological responses under uncharacterized conditions will require large numbers of imaged cells and large quantities of raw imaging data which are far in excess of that presented in this initial proof-of-concept work.

Following the baseline activity of the first 30 minutes, Visualization 2 shows that shortly after exposure to the contractility inhibitor at the 30 minute time-point, the cell and embedded scattering beads in both the en face and cross-sectional planes can be seen to relax back toward their initial locations. This relaxation is not complete, suggesting plastic deformation of the substrate, ECM remodeling, and/or some lingering contractility of the cell. Shortly before the end of the sequence, the cell can be seen to begin contraction again, though further observation of this event was cut short by the end of the experiment, as the behavior was not noticed until the data was studied in post-processing. Overall, the observed particle displacements show that cell contractility inhibition by cytochalasin D is measurable through the analysis of OCM data.

3.3 Time-lapse particle displacements at selected positions around cell bodies

Single particle tracking was used to assess the temporal dynamics of both samples. The results of both manual tracking of individual particles and automated DIC-based displacement tracking are depicted in Figs. 3, 4, and 5. The manual tracking can be considered the ‘true’ displacement. Discrepancies between the manual and automated tracking may arise from several factors affecting the automatic tracking, including the low-pass nature of the cross-correlation operation, measurement noise, interpolation error, and median filtering. However, both cases exhibit comparable results, and may both be used to distinguish the behaviors of the control and contractility-inhibited samples.

 

Fig. 3 Automated and manual tracking of embedded polystyrene bead cumulative displacement magnitudes in time (top) at varying locations around the cells (bottom) exposed to (a) pure DMSO, or (b) cytochalasin D dissolved in DMSO. All displacement magnitudes are defined with respect to the initial location of a given bead. Solid curves depict results of manual single particle tracking; dashed curves depict results of automated DIC-based displacement tracking. The vertical dashed lines in the displacement plots mark the time at which the DMSO or cytochalasin D was added to the samples. Scale bars = 50 µm.

Download Full Size | PDF

 

Fig. 4 Automated and manual tracking of cumulative displacement magnitudes undergone by embedded polystyrene beads at various selected locations around fibroblasts exposed to the control conditions (DMSO). Each subplot (a-h) depicts the results obtained from independent trials of the experimental protocol. The first subplot (a) depicts the same data discussed in Fig. 3(a). All displacement magnitudes are defined with respect to the initial location of a given bead. Solid curves depict results of manual single particle tracking; dashed curves depict results of automated DIC-based displacement tracking. The vertical dotted lines mark the time at which the samples were exposed to DMSO.

Download Full Size | PDF

 

Fig. 5 Automated and manual tracking of cumulative displacement magnitudes undergone by embedded polystyrene beads at various selected locations around fibroblasts exposed to the contractility inhibiting conditions (cytochalasin D + DMSO). Each subplot (a-h) depicts the results obtained from independent trials of the experimental protocol. The first subplot (a) depicts the same data discussed in Fig. 3(b). All displacement magnitudes are defined with respect to the initial location of a given bead. Solid curves depict results of manual single particle tracking; dashed curves depict results of automated DIC-based displacement tracking. The vertical dotted lines mark the time at which the samples were exposed to cytochalasin D solution.

Download Full Size | PDF

As shown in Fig. 3, in both the control and contractility inhibition samples (the same samples discussed in Section 3.2), particles far from the cell bodies exhibit little displacement from their initial positions (less than 2 µm) over the 90 minute duration of both experiments. This helps confirm that the displacements measured in this study are true cell-induced forces and are not due to other factors. In contrast, particles near cells or placed in the leading or trailing paths of cell migration exhibited much greater total displacements (with magnitudes up to 10 µm). In the control experiment, the cells continually deform the ECM over time, causing particles to displace progressively further from their points of origin. The sample exposed to cytochalasin D exhibits distinctly different ECM deformation dynamics. Ten minutes after introduction of the contractility inhibitor, particle displacements attain a maximum, shortly followed by motion back toward their points of origin. After the 70-minute time-point, this ECM relaxation appears to level off at approximately half of the peak measured displacements. As stated previously, this partial relaxation may be due to plastic deformation/remodeling of the ECM or remaining cell forces not inhibited by the cytochalasin D during the experiment.

The trends discussed in these example cases are shared throughout the experimental groups, which are summarized in Figs. 4 and 5. For both the control and contractility-inhibited groups, a total of 8 separate samples were prepared, imaged, and analyzed for each case, using the same procedure as that used to obtain results presented in Fig. 3. As shown in Fig. 4, cells exposed to DMSO tend to continually deform the ECM in an approximately linear fashion throughout the duration of the imaging experiment. Though each sample exhibited this trend, the rate of ECM deformation is variable. Particles located near the leading edges of cells exhibited average displacements rates of 76 ± 29 nm/min (where the error bounds here, and in the remainder of this section, denote plus or minus one standard deviation). We attribute this variation in displacement rates to both biological variability in cell behavior, and to variability in the selected tracking location with respect the cell boundary and orientation.

Similar to the control case, samples exposed to cytochalasin D (as shown in Fig. 5) also exhibit the approximately linear ECM deformation behavior (with average leading edge deformation rates of 102 ± 45 nm/min.) until sometime after the introduction of cytochalasin D. As the compound takes effect, particles embedded in the ECM relax back toward their initial positions. The features of this behavior manifest cell-to-cell variability. These include the time that passes until relaxation begins (16 ± 8 min.), the average rate of relaxation (−56 ± 60 nm/min.), and the amount (percentage) of relaxation that occurs before reaching a steady state or the end of the experiment time (relaxation to 46 ± 30% of the peak displacement). Such variations may not only result from biological variability, but from other factors such as the time required for cytochalasin D to diffuse to the imaged cell.

4. Discussion

4.1 Sample preparation for effective cross-correlation-based displacement tracking

The 3D displacement tracking shown in this work is derived from performing 2D cross-correlation-based displacement tracking in transverse and vertical planes. Full 3D displacement tracking throughout an imaged volume is possible, in principle, using digital volume correlation and/or single particle tracking [15, 17, 22, 23]. A key parameter for future work will be the density of scattering particles present in a sample. Each particle acts as a sampling point of the local deformation field. Successful tracking using correlation-based methods benefits from the use of windows that capture several particles at once, so that the arrangement can form a distinct structure that can be uniquely tracked over time. As such, higher particle density allows for the capture of the deformation field with a higher spatial resolution. In this work, the bead concentration resulted in a mean particle spacing of approximately 12.5 µm. By comparison, recent work in TFM makes use of mean particle spacings in the range of 1.5 µm to 25 µm, though high sampling frequency methods tend to keep within the 1.5 – 3 µm regime [25].

A key assumption of cross-correlation based methods is that local deformations are the result of purely translational motion. Other forms of deformation can degrade or corrupt results (if not accounted for with additional processing algorithms) [22, 23]. However, the high temporal sampling frequency offered by OCM could help reduce the severity of this problem, by capturing more intermediate steps between deformation states. Future work will have to investigate the limits of OCM for detecting small and variable deformations in artificial samples. More advanced algorithms may be adopted from the existing TFM literature to build ‘smarter’ and more efficient algorithms than those discussed so far [23]. Looking toward the future, tracking native scattering features or speckle patterns may offer an avenue toward new processing techniques and investigations of cell traction forces in vivo. Such algorithms are already in development for application within the optical coherence elastography community [53, 54].

4.2 Computing cellular traction forces from measured displacement fields

Cellular tractions are related to ECM deformations through the constitutive relationship between stress and deformation in the ECM. Reconstructing cell traction forces from raw displacement data is a formidable challenge for 3D TFM. Reconstruction is an ill-posed problem, and requires extensive knowledge of local material properties, boundary conditions, and reliable regularization techniques [22, 58]. The most sophisticated method of solving this inverse problem is to use finite element analysis, but even the use of this approach to date has not accounted for the complexities of biological systems. Many TFM studies have used simple linear elastic, isotropic hydrogels as a substrate for easy experimentation and modeling. However, both natural and artificial ECM can violate many of these conditions. Biological materials are often viscoelastic, nonlinear, anisotropic, and have varying properties at different spatial scales [56, 59, 60]. Furthermore, as cells are capable of modifying their environment, even successful measurement of these properties in the neighborhood of a cell will only be valid over a limited timeframe. For these reasons, TFM research must include a thorough accounting of assumptions and limitations of the approach that is used to compute cellular forces. If cell forces cannot be calculated, other metrics of cell behavior, such as mean deformation mechanics [61] or ECM strain energy [62], may still extract useful information from an experiment. In this work, the raw displacement data yielded results that could be used to reveal differences in cell behavior between the control and experimental groups. If cell traction forces or a quantification of ECM remodeling are desired, the field of TFM and mechanobiology at large will require new methods to capture the complex properties of biological materials.

4.3 The ‘big data’ problem, TF-OCM system design, and data acquisition/reconstruction strategies

Volumetric TF-OCM requires cellular resolution over a cuboidal field-of-view with side dimensions of ~10² μm for single/isolated cell studies, or side dimensions of ~10³ μm for studies on cell populations/networks. When coupled with the requirement of high temporal resolution, in order to provide sensitivity to fast cellular mechano-transduction events on the minute time-scale or shorter, this will require the acquisition, and (ideally real-time) processing, of large (‘big data’) volumetric data sets. A major challenge to the practical implementation of TF-OCM will be to obtain as much relevant biological information as possible per volumetric acquisition, while minimizing the amount of raw data that must be stored and processed to reconstruct images.

In the present study, we employed CAO to increase the usable OCM depth range from ~75 μm (limited by degraded transverse resolution due to defocus) to ~200 μm (limited by photon collection of the confocal response). Furthermore, we pre-aligned the location of isolated cells to be near the focal plane at the beginning of each experiment, so as to improve signal strength and reduce the distance from the focal plane that needed to be reconstructed with CAO. Finally, analysis of the resulting 150 GB of raw spectral data per experiment (~2.5 TB in total) was streamlined to ‘ignore’ depths (see Section 2.2 for details) that were outside the usable CAO-OCM depth range. As a result, only 1.4 TB of useful space domain data needed to be reconstructed, as opposed to the 19.2 TB that would have been required for a full volumetric OCT reconstruction using standard methods (assuming double-precision complex numbers in MATLAB). Using these combined methods, a greater amount of useful information could be obtained from each data set acquired, and computational resources were not wasted reconstructing data that contained no useful information to contribute toward the results of this work.

In the more general case, such as the 3D tracking of ECM displacements due to multiple isolated cells or cellular networks, or when imaging multiple cell cultures to obtain statistically significant biological results, employing additional strategies to address the ‘big data’ problem will become essential. One approach, in the case of sparse or low density samples, is the use of compressed sensing [63]. Another approach is to maximize the information content of each volumetric CAO-OCM data set, for example through alternative hardware designs. In particular, the rapid degradation of signal with distance from focus that currently limits the usable depth range could be improved through the use of an astigmatic optical system to equalize photon collection vs. depth, while maintaining diffraction-limited focal-plane resolution throughout the volume with CAO [35]. This could be coupled with real-time volumetric visualization to optimize sample pre-alignment at the beginning of each TF-OCM experiment.

Other OCT systems designs can provide alternative approaches to address the ‘big data’ problem. The use of FF-OCT [40, 41, 54], where the scanning of the reference arm can be coupled with focus tracking, can provide highly efficient data acquisition spanning only the depth range of interest. This approach offers a simpler image reconstruction procedure (relative to CAO-OCM), together with enhanced focal-plane photon collection. Another approach offering enhanced depth-dependent photon collection for volumetric imaging is full-field swept-source (FF-SS)-OCM [64]. This approach can offer ultra-high volumetric acquisition rates, but the use of spatially-coherent full-field illumination makes FF-SS-OCM more susceptible to cross-talk/multiple scattering than beam-scanning OCM systems [65, 66]. An alternative fast imaging scheme with enhanced depth-dependent photon collection is parallel Fourier domain OCT, which uses line illumination to collect A-scans in parallel, with reduced cross-talk effects compared to FF-SS-OCM [67]. Photon collection in beam-scanned OCM can be enhanced via the use of Bessel beam illumination, that can provide spatially-invariant transverse resolution and uniform photon collection over an extended depth range [68, 69]. This approach requires use of separate illumination and collection beams (that can be difficult to align in practice) to enhance SNR [70, 71], and is yet to demonstrate cell imaging over a depth range greater than ~200-300 μm.

4.4 Potential synergy with optical coherence elastography

TF-OCM may benefit from ongoing research in the field of optical coherence elastography (OCE). OCE is a technique used to infer sample mechanical properties by measuring material deformations resulting from a known perturbation. The scientific problem addressed by OCE is very similar to that of TFM [36, 46, 47]. Both forms of microscopy make use of the fundamental relationship between forces, ECM properties, and deformations. They simply differ in which part of this relationship is treated as an unknown quantity, and from what source the deformations originate (externally applied vs. cell-induced). In fact, OCE could benefit TF-OCM by providing access to the ECM mechanical properties required for reconstructing cellular forces in mechanically heterogeneous media. We envision that, compared to existing approaches for TFM, an OCM platform could provide a more complete measurement of cell-ECM interactions, where OCE is used to quantify material properties of a substrate or tissue, followed by TF-OCM to monitor how cells behave and modify the environment. The ability to perform both OCE and TFM in parallel, to monitor dynamic changes in mechanical properties and cellular forces in 3D environments, would constitute a significant advance that could support studies into the dynamics of biological systems that are not possible with current methods used in mechanobiology research.

4.5 Advantages of an OCM-based platform for cell mechanics research

The use of OCM in TFM studies may offer several advantages over current confocal microscopy-based methods. Spectral-domain and swept-source OCT collect signal across all image depths simultaneously, allowing volumetric acquisition with only a 2D transverse raster scan [30]. While OCM images do suffer from aberrations, defocus, and a decrease in signal strength away from the focal plane, the development of hardware- and computation-based methods to mitigate these issues is an active and promising area of research in the field. For example, alternative signal collection schemes like Bessel beam illumination, dynamic focusing [71], and Gabor domain OCM [72] help maintain resolution and signal collection over large depth ranges. Computational methods like interferometric synthetic aperture microscopy [34] and CAO [35, 73–75] allow for the mitigation of defocus and aberrations, respectively. CAO may also enable the use of novel illumination schemes to enhance signal collection in depth [35]. Current techniques and continuing developments in these areas allow the imaging depth range of OCT and OCM to match or surpass that of confocal microscopy. In fact, the samples made for TFM may form an ideal setting for application of computational adaptive optics, as particles added to a sample for displacement tracking can also serve as ‘guide stars’, enabling optimization of resolution and image quality in volumetric data sets, post-acquisition [73, 74]. Alternatively, as done in this work, particle intensities can be used for image-metric-based CAO algorithms.

OCM is based on scattering contrast using near-infrared illumination, effectively eliminating the concern of photobleaching and phototoxicity common to fluorescence microscopy [15, 17]. This feature, in combination with rapid acquisition, enables an unprecedented temporal sampling frequency for volumetric imaging that may be maintained for extended experiment times in both low- and highly-scattering biological samples. As such, OCM may enable the study of large 3D cell populations and/or collectives for both long term processes like migration, division, and ECM remodeling, as well as short term dynamics such as environmental probing with filopodia and mechanotransduction events [27]. Further down the line, OCM may enable novel mechanobiology research in in vivo settings using endogenous scattering contrast.

The high temporal sampling enabled by (especially swept-source and full-field) OCT may also allow for reliable feature/particle tracking with less computation-intensive algorithms. As the time between acquisitions decreases, deformations between images will also decrease. This means that particle or feature tracking algorithms will have smaller search spaces to analyze and may return more robust results due to less structure decorrelation between time-points [23]. It is also possible that new dynamic processes on the minute-to-second timescales may be revealed by these high spatiotemporal resolution acquisitions [3, 27]. Finally, the high phase sensitivity of frequency domain OCT methods may allow the detection of subwavelength (down to subnanometer) displacements at time scales up to the order of the system line scan rate [36]. It is not currently known if any such displacements will be present, but the possibility of detection nevertheless exists.

5. Conclusion

We have demonstrated the first usage of OCM to quantify dynamic cell-induced 3D displacements for use in TFM research, marking a first step in the development of TF-OCM as a new imaging technique in the field of mechanobiology. By imaging and analyzing cells exposed to both control and contractility-inhibiting conditions, we have shown that our technique allows for comparisons of cell traction forces in an in vitro system of interest to cell mechanics research. The high-speed volumetric acquisition capabilities of TF-OCM offers several advantages over existing techniques, and may create opportunities for new discoveries in cell mechanics at previously inaccessible spatiotemporal scales, in both low- and highly-scattering biological media. With the potential to leverage the existing foundation and future advances in both traditional TFM and elastography, TF-OCM is primed for future development and application as a technique for advancing mechanobiology.