Microrheological quantification of viscoelastic properties with photonic force optical coherence elastography

: Photonic force optical coherence elastography (PF-OCE) is a new approach for volumetric characterization of microscopic mechanical properties of three-dimensional viscoelastic medium. It is based on measurements of the complex mechanical response of embedded micro-beads to harmonically modulated radiation-pressure force from a weakly-focused beam. Here, we utilize the Generalized Stokes-Einstein relation to reconstruct local complex shear modulus in polyacrylamide gels by combining PF-OCE measurements of bead mechanical responses and experimentally measured depth-resolved radiation-pressure force profile of our forcing beam. Data exclusion criteria for quantitative PF-OCE based on three noise-related parameters were identified from the analysis of measurement noise at key processing steps. Shear storage modulus measured by quantitative PF-OCE was found to be in good agreement with standard shear rheometry, whereas shear loss modulus was in agreement with previously published atomic force microscopy results. The analysis and results presented here may serve to inform practical, application-specific implementations of PF-OCE, and establish the technique as a viable tool for quantitative mechanical microscopy. a perfect lateral alignment between the of the PF and the of the beads, initial of the displacement of the beads (e.g., via during the measurements. Nevertheless, our results show that quantitative PF-OCE is in good agreement with standard shear rheometry for the characterization of the stiffness of the PAM gels.


Introduction
Optical micromanipulation was first demonstrated by Ashkin in 1970, when he showed that radiation pressure from a low numerical aperture (NA) beam was able to accelerate neutral microparticles in aqueous suspension [1]. This work paved the way for his invention (more than a decade later) of the now ubiquitous "optical tweezers" (OTs) [2]. OTs have since become the cornerstone for numerous discoveries in single-molecule biophysics and nanoscale sciences [3][4][5]. In addition, OTs are also widely adopted in both passive and active microrheology [6,7] for quantification of microscopic viscoelasticity of soft materials, including biological hydrogels and living cells [8][9][10][11]. Meanwhile, optical micromanipulation by low-NA radiation pressure has led to relatively fewer applications in the life sciences [12][13][14].
Using optical coherence tomography (OCT) to monitor the dynamics of accelerated beads in real-time, we recently performed the 'OCT-version' of Ashkin's seminal experiment and measured depth-resolved radiation-pressure force profile from a low-NA Gaussian beam on polystyrene micro-beads in various viscous fluids [15]. We found that radiation-pressure force on the order of 0.1 pN/mW could be achieved depending on the characteristics of the forcing beam as well as the properties of the beads and the medium. Notably, the motions of the beads induced by the radiation-pressure force extended over hundreds of micrometers about the focal plane of the forcing beam. The combination of radiation pressure from a low-NA beam and OCT detection of particle displacements has the potential to open up new modes of quantitative optical micromanipulation and sensing with large depth coverage.
techniques that utilizes spatially localized mechanical excitation provided by a focused optical beam, and makes use of embedded microparticles similar to prior active microrheological approaches [8][9][10]. Other OCE techniques that utilize focused optical excitation include nanobomb OCE [20] and laser-induced surface acoustic wave OCE [21], however, they are fundamentally based on different principles of mechanical excitation and subsequent reconstruction of mechanical properties. These techniques use absorption (as opposed to radiation pressure arising from the transfer of photon momentum) of a short-pulse focused laser to generate localized heating, which generates propagating mechanical waves resulting from rapid thermo-elastic expansion. Furthermore, absorption-mediated photothermal effects, including the thermo-optic effect and thermal expansion, are the physical basis for photothermal OCT (PT-OCT) [22][23][24], a functional variant of OCT imaging that detects optical path length changes associated with the presence of absorbers in the sample. Meanwhile, magnetomotive OCE also performs microrheological measurements on embedded particles, but is based on displacements of embedded magnetic nanoparticles in response to a modulated external magnetic field [25]. Beyond elastographic measurements, particle tracking by OCT has also been adopted for passive microrheological characterization of fluid viscosity [26] and pore sizes of polymer networks [27].
Our PF-OCE technique utilizes harmonically modulated radiation-pressure force to apply localized mechanical excitation on beads randomly distributed over an extended depth range in viscoelastic medium, and leverages phase-sensitive OCT to detect the resulting bead oscillations (on the order of picometers to a few nanometers) [16]. After compensating for the accompanying photothermal response (the same effects utilized in PT-OCT) of the medium [16], we can isolate the complex mechanical response of each bead from the measured changes in optical path length. With the potential to reconstruct viscoelastic properties from bead mechanical response, our preliminary experiments demonstrated PF-OCE as a new approach for performing high-throughput volumetric microrheological measurements of 3D engineered biological systems [16]. Nevertheless, the quantification of viscoelastic properties from bead mechanical response and the challenges associated with that process have not been explicitly addressed in our previous study.
In this paper, we combined our PF-OCE approach [16] with depth-resolved radiationpressure force measurement [15] to quantitatively reconstruct microscopic complex shear modulus around each probed bead in polyacrylamide (PAM) gels, and compare our quantitative PF-OCE results to shear rheometry. In order to obtain reliable measurements of bead mechanical response for quantitative PF-OCE, we analyzed measurement noise at key processing steps and identified data exclusion criteria based on three noise-related parameters. We also discuss modifications that have been made to the previously presented technique [16] in order to improve the performance of PF-OCE for quantitative reconstruction of viscoelastic properties. The analysis and results presented here can serve as a guide to inform future implementations of PF-OCE for specific applications, such as the study of cell-induced microscopic spatiotemporal variation in extracellular matrix (ECM) mechanical properties [9,28] in the rapidly growing field of mechanobiology [29][30][31][32].

Experimental setup
The optical setup for PF-OCE ( Fig. 1(a)) was based on a spectral domain (SD)-OCT system with an additional forcing beam (PF beam) for providing radiation-pressure force, F rad (z), excitation, combined in free-space with the OCT sample arm beam, as previously described [15,16]. The SD-OCT system used a broadband superluminescent diode source (Thorlabs, LS2000B) with a center wavelength of 1300 nm and a full-width-half-maximum (FWHM) bandwidth of 200 nm. Spectral data was acquired with a spectrometer (Cobra 1300, Wasatch Photonics) with a bandwidth of 245 nm and a line scan camera (GL2048, Sensors Unlimited) with 2048 pixels. The axial resolution of the OCT system was 3.7 µm in air. The PF beam was provided by a laser diode at the wavelength of 789 nm (Frankfurt Laser Company, FLU0786M250, HI780 fiber output). Both the OCT and the PF beams were focused by an air objective lens (Olympus, LCPLN20XIR). The OCT beam had an effective NA of 0.28 with a transverse resolution of 2.3 µm. The PF beam had an effective NA of 0.18 and a FWHM focal spot size of 2.2 µm. Fig. 1. Experimental setup and sample configuration. (a) The optical setup consisted of an SD-OCT system and a PF beam combined in free space with the OCT sample arm beam such that the foci of both beams were co-aligned in 3D space. Current-control input from a function generator modulated the power of the PF beam. SLD: superluminencent diode, LD: laser diode, PR: photoreceiver, LP: long-pass dichroic filter, BCM: beam control module, XY: twoaxis galvanometer, OBJ: objective lens. (b) Profile of F rad (z) as a function of depth, z, w.r.t. the focal plane of the PF beam at 78-mW PF beam power. (c) Peak F rad as a function of PF beam power, P, at the sample. Linear best-fit line is shown. In (b) and (c), F rad (z) was measured on 1.7-µm polystyrene beads in 10% glycerol-water mixture (refractive index 1.3469). The preparation and the naming convention of chemically-crosslinked PAM gels was based on Tse et al [33]. First, acrylamide monomer (40% Acrylamide Solution, Bio-Rad), bisacrylamide crosslinker (2% Bis Solution, Bio-Rad), and deionized water were mixed at desired concentration (Table 1). Then, aqueous suspensions of 1.7-µm diameter (Spherotech, PP-15) and 0.1-µm diameter (Sigma-Aldrich, LB1) polystyrene beads were added to the mixture at concentrations of 28 µL/mL (mean particle separation of 15 µm) and 1.2 µL/mL (mean particle separation of 2 µm) mixture volume, respectively. The 1.7-µm bead size was selected for optimal bead oscillation amplitude when forced by the harmonically-modulated 2.2-µm PF beam. The 0.1-µm beads provided background scattering signals in the medium and served as photothermal reporters [16]. After careful stirring, the monomer-crosslinker-bead solution was desiccated for 15 minutes. Then, a portion of the solution was set aside for bulk mechanical testing by shear rheometry.

Sample preparation
To the remaining solution, redox initiator, 10% ammonium persulfate (APS, Bio-Rad), and catalyst, tetramethylethylenediamine (TEMED, Bio-Rad), were added at the concentrations of 10 µL/mL and 1 µL/mL, respectively. After mixing, 100 µL of the activated solution was pipetted into the well in a glass-bottomed imaging petri dish. Another glass coverslip was placed on top of the well, concealing the PAM gel, and the edges were sealed with mineral oil to prevent evaporation. The gel was allowed to polymerize for 60 minutes prior to the PF-OCE measurements. One sample was prepared per PAM gel concentration.

System alignment procedure and data acquisition
Prior to each experiment, co-alignment between the foci of the OCT beam and the PF beam in 3D space was verified by comparing the detected OCT signal against the confocal response of the PF beam acquired with a photoreceiver (Newport, 2051-FS), as previously described [16]. Beam alignment was adjusted by the beam control module ( Fig. 1(a)) as needed. A total of 13 M-mode OCT data sets were acquired in each sample. The first 9 data sets were acquired with a 1.7-µm bead located at the focal planes of the PF and OCT beams, and aligned to the optical axis of these beams-these provided PF-OCE measurements on 9 different beads per PAM gel concentration. The remaining 4 data sets were acquired without a 1.7-µm bead directly in the beam path-these served as additional data for depth-dependent photothermal response measurements.
Each M-mode data set consisted of 10 frames and 32000 A-scans per frame, acquired at a line-scan rate of 10 kHz and an exposure time of 16 µs. The function generator (Tektronix, AFG3051C) harmonically modulated the output power of the PF beam at the modulation frequency of 20 Hz and peak power of 78 mW at the sample ( Fig. 1(a)), corresponding to a peak force of 12 pN at the focal plane ( Fig. 1(c)). The modulation was synchronized to the acquisition of A-scans such that the PF beam power remained at 0 mW until the modulation began at the 321st A-scan in each frame. The sample was first positioned axially such that the focal plane of the OCT beam (co-aligned with the PF beam focal plane) was located roughly at the center of the PAM gel (i.e., mid-way between the petri dish bottom and the coverslip placed on top of the well). Then, the sample was only translated laterally to change the imaging location between data sets. For the 9 data sets containing a 1.7-µm bead at the focal plane, the bead was positioned in the OCT beam path (co-aligned with the PF beam path) by laterally translating the sample to maximize the M-mode OCT intensity of the bead at the focal plane. This procedure was repeated for each of the 9 beads measured per sample.

Data processing
M-mode OCT images, S(z,t), were reconstructed according to a standard SD-OCT reconstruction procedure (background subtraction, spectrum resampling, dispersion correction, and inverse Fourier transformation). The complex M-mode OCT images are given by: where t denotes the time at each A-scan and z denotes the depth w.r.t. the focal plane of the PF and OCT beams at each pixel depth. The following processing steps were performed independently at each pixel depth to obtain the oscillatory response induced by the harmonically modulated PF beam. First, average OCT signal-to-noise ratio (SNR) was calculated for each pixel depth, as previously described [16]. Then, pixel depths with OCT SNR < 10 dB were removed and excluded from further processing. The rationale supporting this exclusion criteria will be discussed in Section 3.2. The raw optical path length changes, ΔOPL(z,t), induced by the harmonic modulation was extracted from the raw phase of the complex OCT signal. This phase must be registered to the reference phase obtained from the petri dish bottom surface in order to eliminate any systematic phase drifts in the system [16]. The registered phase angle was extracted from the product of the complex OCT signal at each pixel depth and the complex conjugate of the OCT signal at depth z ref , corresponding to the pixel depth with maximum average OCT intensity on the petri dish bottom surface: where λ 0 and n denote the center wavelength of the OCT beam and the refractive index of the medium, respectively. The resulting phase-registered ΔOPL response was averaged across all 10 frames to suppress random noise in the phase of the complex OCT signal. From the phaseregistered ΔOPL(z,t), the noise amplitude, δA(z), of the ΔOPL response was calculated from the root-mean-square of the frequency spectrum adjacent to the response peak at the 20-Hz modulation frequency [16]. Then, the ΔOPL response at the modulation frequency was isolated via multiplication with a brick-wall filter (20-Hz center frequency and ± 0.4-Hz passband) in the frequency domain. The complex filtered ΔOPL response can be expressed as the product of a depth-dependent phasor and the PF drive waveform: is the PF drive waveform expressed as a complex exponential of unit amplitude with a phase shift φ drive . A(z) and φ(z) denote the amplitude and phase shift w.r.t. h(t) of the ΔOPL response at each pixel depth, respectively. Note that with this notation, φ < 0 corresponds to a phase delay w.r.t. h(t). For brevity, any mention of phase shift from this point on shall imply a shift w.r.t. h(t), according to Eq. (3). Given that the sample is assumed to behave linearly, and the modulation was provided at a single frequency, the oscillation amplitude and phase shift are expected to be constant in time. Thus, A(z) and φ(z) were obtained by taking the median of the amplitude and the phase shift of ΔOPL(z,t) over the duration that the modulation was applied in each frame (approximately 3 s): where med{…} t denotes median operation along t. After obtaining ΔOPL(z,t), δA(z), A(z), and φ(z), the data was further excluded based on oscillation amplitude-to-noise ratio, A/δA, where pixel depths with A/δA ≤ 6 were removed and excluded from further processing. Then, the remaining data was segmented into a photothermal response data region (corresponding to the 0.1-µm beads and the coverslip surface) and a total response data region (corresponding to the 1.7-µm beads) based on OCT SNR thresholds, as previously described [16]. However, the specific threshold values were adjusted to account for differences in the noise performance of the system, acquisition schemes, and bead sizes; values used in the presented experiments were 10 dB ≤ OCT SNR < 20 dB for the 0.1-µm beads and OCT SNR ≥ 22 dB for the 1.7-µm beads.
After the segmentation, the photothermal response data region from all 13 data sets acquired per sample was combined to obtain a single depth-dependent photothermal response amplitude, A PT (z), and phase shift, φ PT (z), for each PAM gel concentration. A PT (z) and φ PT (z) were obtained via nonlinear least-square curve-fitting of the theoretically simulated depthdependent photothermal response amplitude and phase shift curves [16] to the measured A(z) and φ(z) data (in the photothermal response data region), respectively. The obtained best-fit curves represented A PT (z) and φ PT (z) for each PAM gel concentration.
Mechanical response (oscillation) of the 1.7-µm beads induced by the harmonically modulated F rad was isolated from the filtered ΔOPL(z,t) in the total response data region after compensating for the photothermal response by: where ΔOPL mech (z,t), ΔOPL tot (z,t), and ΔOPL PT (z,t) denote the complex mechanical, total, and photothermal responses at each pixel depth, respectively. ΔOPL mech (z,t) can also be expressed as the product of a depth-dependent phasor and the PF drive waveform: where the mechanical response amplitude, A mech (z), and phase shift, φ mech (z), were obtained from Eqs. (4) and (5), respectively. In order to assess the reliability of the measured A mech (z) and φ mech (z), the standard deviation of the mechanical response amplitude and phase shift, denoted by σ A (z) and σ φ (z), were computed over the duration that the modulation was applied in each frame. This matter will be discussed fully in Section 3.2.
After obtaining A mech (z), φ mech (z), σ A (z), and σ φ (z), the results were further excluded based on the stability in the time domain of the ΔOPL mech (z,t) response, where pixel depths with σ A /A mech ≥ 0.1 or σ φ ≥ π/16 were removed and excluded from further analysis. The remaining pixel-wise A mech (z) and φ mech (z) results from the pixel depths comprising each 1.7-µm bead were grouped together and the median values were taken to represent the mechanical response amplitude and phase of each bead [16]. These bead-wise A mech and φ mech results were used to reconstruct the local complex shear modulus, G G iG * ′ ′ ′ = + , of the PAM medium around each bead via the Generalized Stokes-Einstein Relation (GSER) [34,35] where z b denotes the depth at the center of the beads w.r.t. the focal plane of the PF beam and a denotes the radius of the beads. Refer to Section 2.5 for the experimental measurement of depth-resolved F rad (z). We reiterate that φ mech is the phase shift of ΔOPL mech w.r.t. h(t), where a negative value of φ mech implies that bead oscillation lags behind the mechanical excitation provided by the modulated F rad . In other words, the mechanical 'phase delay' (as referred to in rheometry) is given by -φ mech .

Radiation-pressure force measurement
Depth-resolved measurement of the axial profile of F rad (z) with OCT followed a previously reported method [15]. Briefly, 1.7-µm beads were dispersed in a solution of 10% glycerol in distilled water (approximating the refractive index of the PAM gels). Then, the beads were accelerated by the PF beam operating at a constant power, during which the axial trajectories of the beads were tracked in real-time via M-mode OCT imaging. Instantaneous velocity and acceleration of each bead were computed from its trajectory. Then, the axial equation of motion for a spherical bead in viscous fluid was solved to obtain the instantaneous magnitude of F rad at each depth. The axial profile of F rad (z) at the PF beam power of 78 mW at the sample is shown in Fig. 1(b). At the focal plane, the peak force per unit PF beam power was 0.15 pN/mW (Fig. 1(c)).

Shear rheometry
Bulk complex shear modulus of each PAM sample was measured with a parallel-plate shear rheometer (TA Instruments, DHR-3) using a 20-mm diameter plate. APS solution and TEMED were added to the monomer-crosslinker solution immediately before each measurement. After mixing well, 200 µL of the activated polymer solution was pipetted onto the bottom plate of the rheometer before the top plate was lowered down to achieve a gap of 500 µm. Excess polymer solution was wiped away, then, the sample was sealed around the gap with mineral oil to prevent evaporation. Polymerization progress was monitored in a time-sweep oscillatory test at oscillation frequency of 1 rad/s and shear strain of 0.5%, measuring the complex shear modulus as a function of polymerization time. All samples were left to polymerize for 60 minutes; stabilization of shear moduli within this time frame confirmed complete polymerization of the PAM gels. After polymerization, frequency-sweep oscillatory test was performed at the oscillation frequency ranging from 1 to 50 Hz and applied torque of 10 µN⋅m. Two different samples were tested per PAM concentration. The test was repeated 2 times on one sample and 3 times on the other, totaling 5 measurements per PAM gel concentration.

Statistical analysis
For the comparison of complex shear modulus measured by quantitative PF-OCE to shear rheometry measurements, Pearson's linear correlation test was performed on the median value of G ′ and G ″ measured on all beads in each PAM gel concentration by PF-OCE against the mean value of G ′ and G ″ from 5 measurements in each PAM gel concentration by shear rheometry.

Bead response to harmonically-modulated radiation pressure
An example of the reconstructed space-domain M-mode OCT image shows a 1.7-µm bead located at the focal plane of the OCT beam, which was co-aligned to the focal plane of the PF beam ( Fig. 2(a)). Elsewhere inside the PAM gel, lower intensity background signals from the 0.1-µm beads and the intrinsic scattering of the PAM gel can be observed. The petri dish bottom, which was used for phase registration, appears at the top of the image due to the inverted imaging configuration ( Fig. 1(a)). Similarly, the glass coverslip, which was placed on top of the petri dish well, appears at the bottom of the image. The signal from the surface of the coverslip in contact with the PAM gel encoded the total cumulative photothermal response.
The raw phase of the complex OCT signals, φ OCT , contained a systematic drift over time that was presented at all pixel depths. The oscillatory response of the 1.7-µm bead, ΔOPL, was obtained after registering its phase to that of the petri dish bottom surface (Fig. 2(b)). Note that this ΔOPL represents the total response, which still contains the contribution of the photothermal response of the PAM medium. The 20-Hz response, induced by the harmonically modulated PF beam, is evident in the frequency-domain representation of ΔOPL (Fig. 2(c)). In this data set, the phase noise amplitude, δA, adjacent to the 20-Hz response was 23 pm. After isolating the 20-Hz response via a brick-wall filter, the oscillation amplitude and phase shift were extracted according to Eqs. (4) and (5), respectively. Both the amplitude and the phase shift exhibited variations over time, with the median values of A and φ, and the standard deviations of σ A and σ φ , respectively (Fig. 2(d)). As shown in the inset of Fig. 2(d), a phase lag in the ΔOPL response w.r.t. to h(t) corresponds to a negative-valued phase shift.

Analysis of noise and criteria for quantitative PF-OCE
In order to obtain reliable bead mechanical responses for the quantitative reconstruction of complex shear modulus, PF-OCE data was excluded based on noise-related parameters at several processing steps, as described in Section 2.4. A flowchart summarizes key processing steps and associated exclusion criteria (Fig. 3(a)). These exclusion criteria were based on three factors that contribute to the errors in the PF-OCE measurements of A mech and φ mech : 1) SNR of the OCT signal, 2) oscillation amplitude-to-noise ratio of the ΔOPL response, and 3) instability of the ΔOPL response in the time domain. In the results that follow, measurements from all pixel depths (regardless of the OCT SNR thresholds) in a PAM gel sample are included for analysis purpose.
First, the SNR of the OCT signal has a direct effect on the level of noise that is presented in φ OCT [36,37]. Phase noise amplitude at 20 Hz, δA, in the raw ΔOPL response was plotted as a function of the average OCT SNR at each pixel depth ( Fig. 3(b)). Remarkably, for pixel depths with OCT SNR above 10 dB, δA followed the theoretical shot-noise limited prediction [36]. However, at pixel depths with OCT SNR below 10 dB, the measured ΔOPL contained over an order of magnitude larger δA than the theoretical shot-noise limits, suggesting that the measured signals at these SNR levels were dominated by other sources of phase noise. Based on the trends observed in Fig. 3(b), it was determined that the 20-Hz oscillatory response induced by the modulated F rad could only be reliably obtained at the pixel depths with OCT SNR ≥ 10 dB, where the phase noise obeyed the theoretical shot-noise limited noise floor [36].
Second, as a direct consequence of SNR-dependent δA, the oscillation amplitude-to-noise ratio, A/δA, of the measured ΔOPL response is also influenced by the OCT SNR. In order to accurately measure the amplitude and phase shift of the oscillatory response to the modulated F rad , the oscillation amplitude must be sufficiently larger than the noise amplitude at the modulation frequency. Two data regimes, with higher and lower A/δA, can be clearly observed in the plot of A/δA as a function of average OCT SNR (Fig. 3(c)), where the boundary between these two regimes also predominantly coincide with the OCT SNR cutoff at 10 dB. Based on the two regimes observed in Fig. 3(c), it was determined that only the ΔOPL response with A/δA > 6, corresponding to a relative power of approximately 15 dB, would facilitate reliable calculations of A and φ. Third, any variation in time of the mechanical response amplitude and phase shift, despite the constant modulation frequency, is an indication of uncertainties in the isolated mechanical response. Potential sources of this instability over time include the measurement phase noise, power fluctuation of the PF beam, or any drifting of the 1.7-µm beads relative to the PF beam while each frame is being acquired. As a result, an unstable ΔOPL mech response in time may imply that the corresponding A mech and φ mech results cannot be accurately extracted by Eqs. (4) and (5), respectively. To analyze the instability of the ΔOPL response over time, standard deviations of its amplitude and phase shift, σ A and σ φ , were calculated over the duration of the applied modulation (see Section 2.4) and plotted as a function of the corresponding median values, A and φ (Figs. 3(d) and 3(e)). Larger values of σ A /A and σ φ , indicating a higher degree of instability in the ΔOPL response over time, correspond to erroneous calculated values of A and φ that are spread out over a wide range. Notably, most of these data points also coincide with those that exhibit OCT SNR < 10 dB and A/δA ≤ 6 (light blue markers). This observation may suggest that the instability in the time domain is largely a result of the SNR-dependent phase noise in the measured ΔOPL response. In addition, this observation also supports the effectiveness of the two previous exclusion criteria in selecting the pixel depths with stable response that allowed for reliable calculation of A and φ. Based on the analysis in Figs. 3(d) and 3(e), only A mech and φ mech results from the isolated mechanical response with σ A /A mech < 0.1 or σ φ < π/16 were considered reliable for the quantitative PF-OCE reconstruction of complex shear modulus.

Depth-dependent photothermal responses
In order to isolate the mechanical response of the 1.7-µm beads from the measured total response, the depth-dependent photothermal response of the PAM medium was obtained for each PAM gel concentration. Pixel depths in the photothermal response data region were excluded based on the OCT SNR and the A/δA criteria prior to fitting to the theoretical curves to obtain A PT (z) and φ PT (z) (Fig. 4). Although all pixel depths in the photothermal data region are displayed, only data points from the pixel depths that have passed the first two exclusion criteria (dark blue markers) contributed to the curve fit in Figs. 4(b) and 4(c). For the 4T1C gel shown in Fig. 4, the photothermal response amplitude and phase shift at the focal plane (z = 0) were 1.1 ± 0.2 nm and 2.55 ± 0.12 rad, respectively. Note that the positive-valued phase shifts in Fig. 4(c) does not imply that the photothermal response is in advance of the PF drive waveform; they reflect the decrease in optical path length as the temperature increases [16,38] as a result of the dominating effect of the negative-valued thermo-optic coefficient of water [39], which produces an almost out-of-phase response w.r.t. h(t).

b) and (c)
A PT and φ PT as a function of z, respectively, where z is defined w.r.t. to the focal plane of the OCT beam. Blue markers are data points from the photothermal response data region. Black solid lines are the best-fit curves obtained from fitting the data to the theoretical curves. Red dotted lines indicate ± 1 median absolute difference between the best-fit curves and the data. In (a)-(c), light blue indicates pixel depths that were excluded from analysis based on the first two exclusion criteria. Dark blue indicates pixel depths that were used for the curve fit to produce A PT (z) and φ PT (z).

Reconstruction of complex shear modulus and comparison to shear rheometry
Bead mechanical response amplitude, A mech , and phase shift, φ mech , were obtained for each measured 1.7-µm bead in the PAM gels after compensating for the photothermal response according to Eq. (6), using the respective A PT (z) and φ PT (z) curves (Figs. 4(b) and 4(c)) for each PAM concentration, and enforcing the time-domain instability exclusion criteria described in Section 3.2 ( Fig. 5(a)). A decreasing trend from 3T1C-6T1C was apparent for A mech , consistent with the expected increasing trend in stiffness of the PAM gels as the polymer concentrations increased. For consistency with rheometry convention, the mechanical phase delay, -φ mech , is plotted instead of the calculated phase shift. Under this definition, a completely elastic solid is characterized by the phase delay of 0 rad while a viscous fluid is characterized by the phase delay of π/2 rad. The phase delay was close to 0 rad for all PAM gel concentrations, with the 3T1C gel exhibiting a relatively more viscous behavior than the rest at a median phase delay of π/8 rad.
Complex shear modulus was reconstructed from A mech and φ mech of each 1.7-µm bead via GSER [34,35] using to Eq. (8) (Fig. 5(b)). This represents the microscopic local viscoelastic properties of the medium around each bead. In this study, the PAM gels were assumed to be microscopically homogeneous at the length scale of the 1.7-µm bead size. This assumption is supported by [35,40]. An increasing trend from 3T1C-6T1C could be observed for both G ′ and G ″ . Quantitative PF-OCE results showed statistically significant correlation to bulk shear rheometry measurements at 20-Hz oscillation frequency for G ′ (R 2 = 0.97, p < 0.0001) (Fig.  5(c)) but not for G ″ (R 2 = 0.15, p > 0.1) (Fig. 5(d)). Whereas PF-OCE measured G ″ in the range 40-400 Pa, with a monotonic increasing trend from 3T1C-6T1C, bulk rheometry measured G ″ in the range 10-60 Pa with no apparent trend w.r.t. the polymer concentrations.

Modifications from previous work in order to enable quantitative PF-OCE
We previously demonstrated volumetric PF-OCE measurements of bead mechanical responses in agarose hydrogels, which revealed microscopic heterogeneity in the mechanical properties of the porous medium [16]. In order to be able to utilize PF-OCE as a new tool for microscopic mechanical characterization of soft materials, the technique had to be fine-tuned to improve the accuracy of bead mechanical response measurements and enable quantitative reconstruction of mechanical properties. The first modification was to minimize the magnitude of the photothermal contribution to the measured ΔOPL response. This was achieved by changing the wavelength of the PF beam from 976 nm in our previous study [16] to 789 nm, since the latter has approximately an order of magnitude lower absorption coefficient in water than the former [41]. As a result, A PT measured in this paper with the 789nm PF beam was on the order of 1 nm, a six-fold decrease from the 6 nm with the 976-nm beam in our previous study [16].
The second modification was the redesign of the optical setup to achieve larger bead oscillation amplitude and lower the phase noise floor, in order to maximize A and minimize δA. In a given sample, achievable A mech per unit F rad can be increased by using a smaller bead size, since A mech is inversely proportional to the radius of the bead [16,34,42]. In order to maintain an optimal F rad with the use of smaller beads, the focal spot size of the PF beam must also scale down accordingly [16]. The improved system presented in this paper, with a 2.2-µm PF beam optimized for the 1.7-µm diameter beads, is theoretically expected to demonstrate a 44% increase in A mech per unit PF beam power relative to the system used in our previous study (3.8-µm PF beam and 3-µm diameter beads) [16]. In addition to increasing A mech , the use of smaller beads is also desirable for biological applications with live-cell imaging, where the total volume fraction of exogenous beads in the sample must be kept under the levels that can be tolerated by cells [28,43]. In order to lower the noise floor of the system, mechanical vibration in the optical setup was minimized by lowering the height of entire optical setup on the optical table and designing the collimator arms for both OCT and PF beams to be as compact as possible. Although we did not quantitatively analyze the effects of these factors on the SNR-dependent δA, the improved system demonstrates δA that reaches the theoretical shot-noise limit (Fig. 3(b)) whereas the previous system was not able to [16].
The third modification was to ensure the accuracy of the calculated A mech and φ mech for quantitative PF-OCE by implementing exclusion criteria based on quantifiable noise-related parameters at key processing steps. We found that an OCT SNR ≥ 10 dB was needed to achieve optimal shot-noise limited δA ( Fig. 3(b)). This is particularly important for the photothermal response, which is primarily measured from the low-SNR photothermal reporters (0.1-µm beads). In our previous system using the 976-nm PF beam, an OCT SNR of only 3 dB was used to measure the photothermal response [16]. However, as a direct consequence of the lower photothermal response due to the lower absorption of 789 nm in water, a lower δA is now required to make these measurements. We also found that a stable oscillation amplitude and phase shift could only be calculated from the ΔOPL responses with A/δA > 6 (Figs. 3(c)-3(e)). This corresponds to a relative signal power of approximately 15 dB, which is higher than our previous estimate of 3 dB [16]. Lastly, we also found that the instability of the ΔOPL response over time degraded the accuracy of A mech and φ mech that were calculated from the median values of its amplitude and phase shift, respectively (Figs. 3(d) and 3(e)). This provided the last criteria for obtaining reliable A mech and φ mech values for quantitative PF-OCE.

Reconstruction of complex shear modulus and comparison to shear rheometry
By combining PF-OCE measurements of A mech and φ mech with OCT depth-resolved measurement of F rad (z) [15], the local complex shear modulus of the PAM gel around each measured 1.7-µm bead was reconstructed via GSER (Eq. (8)). We note that under our experimental conditions (e.g., bead size, modulation frequency, sample properties etc [44].), GSER is consistent with the forward prediction of bead dynamics by other more rigorous mathematical models, such as the Oestriecher's model [42] or the solution to the generalized 3D equations of motion for forced harmonic oscillation of a sphere in viscoelastic medium [45,46] (data not shown). Our quantitative PF-OCE results were compared to shear rheometry measurements, under the premise that the PAM gels are microscopically homogeneous at the length scale of the 1.7-µm beads, owing to the mesh size of the polymer networks being on the order of 10-100 nm for chemically cross-linked PAM gels [35,40].
Quantitative PF-OCE was consistent with shear rheometry for the measurements of G ′ , demonstrating a correlation coefficient of 0.97 between the results from the two approaches (Fig. 5(c)). However, results in each PAM gel concentration showed slight discrepancy, with quantitative PF-OCE measuring a larger G ′ (i.e., higher stiffness) in general, suggesting that PF-OCE tends to measure smaller A mech than expected. This can be attributed to any deviation from a perfect lateral alignment between the optical axis of the PF beam and the center of the 1.7-µm beads, either from the initial positioning of the sample or any subsequent displacement of the beads (e.g., via diffusion) during the measurements. Nevertheless, our results show that quantitative PF-OCE is in good agreement with standard shear rheometry for the characterization of the stiffness of the PAM gels.
In contrast, quantitative PF-OCE measurements of G ″ did not show statistically significant correlation to shear rheometry (Fig. 5(d)). In general, quantitative PF-OCE measured larger G ″ (greater relative viscosity) than shear rheometry, and demonstrated an increasing trend as a function of increasing polymer concentration. We note that other microscopic characterizations of PAM gel viscoelasticity using atomic force microscopy (AFM) have yielded G ″ similar to our quantitative PF-OCE results [47,48]. This suggests that the discrepancy between quantitative PF-OCE and shear rheometry may be attributed to different viscoelastic behaviors of the PAM gels at the microscopic versus macroscopic scale. Viscosity in hydrogels arises from the viscous drag due to the flow of fluid and any other frictional interactions in the polymer networks. Based on this physical interpretation, the relative viscosity of the PAM gels is likely to be more affected by the length scale at which the measurements are carried out than the stiffness. Our results suggest that a direct micro-to macroscale connection between different mechanical characterizations of the porous PAM gels may only be valid for elasticity (storage modulus) but not viscosity (loss modulus). For other types of porous hydrogels, particularly those with larger pore sizes such as the agarose hydrogels used in our previous work [16] or those with highly fibrous microstructure such as collagen gels [28], this connection across length scales may not be valid for either elasticity or viscosity if the medium cannot be considered homogeneous at the length scale of the microscopic measurements.

Implications for future applications
By working above the OCT SNR cutoff where the experimental δA obeys the theoretical shotnoise limit, the maximum sample stiffness that can be measured by quantitative PF-OCE (limited by the minimum measurable A mech ) can be estimated based on the A/δA > 6 criterion. For the M-mode acquisition scheme implemented in this paper, the maximum stiffness that can be accurately measured with OCT SNR of 22 dB (minimum threshold for segmenting the 1.7-µm beads) is shear modulus of 5.2 kPa. We note that this is merely an estimate because the actual measured ΔOPL response is the total response, which has the contribution of the photothermal response of the medium. This estimate also applies to 3D mechanical microscopy applications with the previously demonstrated beam-scanning BM-mode acquisition scheme [16].
For the same peak power of the modulated PF beam, several trade-offs exist between maximum measurable stiffness, volumetric throughput, and depth coverage. For a given PF-OCE system, there is a trade-off between volumetric throughput and the maximum measurable stiffness, where the minimum acquisition time per volume is determined by the minimum number of BM-mode frames required per slow-axis position to achieve the desired δA [16]. In order to measure smaller bead oscillation amplitude in a stiffer sample, a longer measurement time is required to achieve lower δA [16,36] and maintain adequate A/δA ratio. For instance, a volumetric measurement at 500-Hz frame rate with lateral pixel size of 1 µm 2 and lateral FOV of 200 µm × 200 µm would require up to 30 minutes of acquisition time in a sample with maximum expected shear modulus of 2 kPa; the same measurement could be made in only 8 minutes if the maximum expected shear modulus were only 1 kPa. For a given desired acquisition time, another trade-off exists between the volumetric depth coverage and the maximum measurable stiffness. On the one hand, the optical system could be designed to have a higher NA (i.e., smaller PF beam focal spot size) so that measurements could be made with smaller beads to increase the achievable A mech in a given sample. This would increase the maximum measurable stiffness, but at the cost of narrowing the depth range in which F rad can be exerted. On the other hand, larger depth coverage could be achieved by lowering the NA of the optical system, but at the cost of decreasing the achievable A mech for a given sample stiffness.
The optical design of the PF-OCE system and the acquisition scheme utilized in each experiment can be adapted to optimize the desired information in an application-specific manner, using the system noise performance and data requirements for quantitative PF-OCE investigated in this paper as a guide. For instance, a study of the spatiotemporal variation ECM mechanical properties induced by the contractile dynamics of a single cell over minuteto-hour timescale may benefit from a relatively high-NA system with limited depth-coverage but large maximum detectable stiffness within 1-10 minutes of acquisition time [49,50]. On the other hand, dynamics over a longer timescale such as the phenotypic change of cells or the migration of cancer cells from cellular spheroids may benefit from a long time-lapsed imaging study with extended depth coverage (on the order of 10 2 µm), but with temporal resolution of 20-40 minutes per volume [51][52][53].

Conclusion
We have demonstrated microrheological quantification of viscoelastic properties of PAM gels with quantitative PF-OCE, combining depth-resolved measurement of radiation-pressure force with OCT [15] and previously introduced PF-OCE technique [16], and compared our results to shear rheometry. Our analysis of measurement noise culminated in data exclusion criteria, based on OCT SNR, oscillation amplitude-to-noise ratio, and instability of the measured response over time, to enable quantitative PF-OCE. These outcomes serve to extend the viability of PF-OCE as a new technique for quantitative mechanical microscopy of viscoelastic media, such as the 3D substrates used in engineered biological systems. The comparison to shear rheometry validates quantitative PF-OCE reconstruction of sample stiffness. However, the comparison between relative viscosities obtained from microscopic measurements by quantitative PF-OCE to bulk measurements by shear rheometry may not be valid. These findings are consistent with previous microscopic measurements with AFM [47,48]. The results presented here can be used to inform the optical system design and data acquisition scheme for specific applications of PF-OCE, such as live-cell imaging studies of the cellular-scale spatiotemporal variation in ECM mechanical properties that are induced by dynamic, biophysical cell-ECM interactions.