Displacement and strain mapping for osteocytes under fluid shear stress using digital holographic microscopy and digital image correlation

: Osteocytes, as the mechano-sensors in bone, are always subjected to fluid shear stress (FSS) from the surrounding matrix. Quantification of FSS-induced cellular deformation is significant for clarifying the “perceive and transmit” process of cellular mechanotransduction. In this research, a label-free displacement and strain mapping method based on digital holographic microscopy (DHM) and digital image correlation (DIC) is introduced. The method, which is termed DHM-DIC, innovatively utilizes surface features extracted from holographic phase images instead of speckles as the metric for DIC searching. Simulation results on a hemisphere validate the feasibility of DHM-DIC. Displacement and strain maps of living osteocytes under 1.5 Pa FSS are evaluated from DHM-DIC and present good agreement with our previous finite element modeling results. the shear 0.0318 perpendicular to shear direction. results validate the feasibility of DHM-DIC estimating full-field cell displacement and strain. of the advantages of DHM-DIC subregion size ( m ( m the estimation error and processing time.


Introduction
Osteocytes, which act as the mechanosensors in bone [1], are always under fluid shear stress (FSS) from the unmineralized bone matrix [2]. The FSS-induced cell strain is a triggering step for activation of downstream cell mechanotransduction [3]. Quantifying FSS-induced cell strain is meaningful for uncovering the "perceive and transmit" process of osteocytes after mechanical stimuli. In recent years, various efforts have been made to quantify the FSS-induced cell strain. Nicolella and Rath innovatively utilize a machine vision photogrammetry on individual cells for strain quantification [4]. Ueki et al. combine the digital image correlation (DIC) with finite element analysis for subcellular strain calculation [5,6]. Baik and Qiu et al. evaluated cell strain in three-dimension (3D) using a novel 3D fluorescence microscopic setup which obtains cell images from both bottom and side view [7][8][9][10][11].
In the above-mentioned studies, cell images with or without fluorescence labelling are captured before (undeformed) and after (deformed) FSS. DIC-based algorithms are used to obtain the FSS-induced displacements on the cell surface. In common DIC processing, inherent or manually added speckles on sample surface are principally required as the metric to search for displacement vectors [12,13]. However, inherent speckles on the cell surface are often indiscernible, and extrinsic speckles usually require additional manipulation such as fluorescence labeling [4][5][6][7][8][9][10][11], microbeads attaching [14,15] or a specialized substrate [15,16].
Digital holographic microscopy (DHM), as an innovative and fast-developing quantitative phase imaging approach [17][18][19], has been widely applied on non-invasive and label-free studies for living cell [20][21][22][23][24]. In DHM, optical path length, which is reconstructed from hologram, contains abundant morphological features of sample surface. These inherent features can be employed as the metrics for DIC searching process to evaluate cellular displacement fields while under FSS.
In this study, a label-free method based on DHM and DIC is proposed for full-field cell displacement and strain measurement. The method, which is termed as DHM-DIC, innovatively utilizes cell surface features extracted from holographic phase image instead of speckle as the metric for DIC matching. Validation of DHM-DIC is first performed using a numerically simulated hemisphere with certain in-plane displacements. Phase maps of osteocytes under 1.5Pa FSS are then recorded by an off-axis DHM setup for 5 minutes with 1 frame per second and full-field displacement and strain maps are obtained by DHM-DIC. The experimental results of mean and maximum strain along the shear direction, which is 0.042 and 0.060, are in good agreement with our previous finite element modeling study [25,26], in which mean strain and maximum strain are 0.029 and 0.051, indicating that DHM-DIC is an effective and promising approach for cell mechanical study.

Cell preparation
The murine long bone osteocyte Y4 (MLO-Y4) cell line is kindly provided by Dr. Lynda F Bonewald of Indiana University [27]. Cells are maintained in alpha Minimum Essential Medium (α-MEM) (Gibco, USA) supplemented with 2.5% (vol/vol) fetal bovine serum (FBS) (HyClone, South America) and 2.5% (vol/vol) calf serum (CS) (HyClone, New Zealand) on culture flask (Nunc, Roskilde, Denmark) pre-coated with rat tail type I collagen (0.15 mg/ml, Millipore, USA) at 37°C, 5% CO 2 . Cells are seeded on glass slides adaptive for the flow chamber at the density of 5000/cm 2 for 24 hours before shear experiments. Serum-free culture medium are used for clearer visual field during shear experiments and image recording.

DHM setup
An off-axis Mach-Zehnder DHM system coupled with a custom-made flow chamber is employed to record real-time holograms of osteocytes under FSS. The schematic of the integrated system is shown in Fig. 1(a). In the setup, a 532 nm laser (Cobolt, Sweden) with an output power of 50 mW is utilized as the light source. A neutral density filter (NDF) is placed after the laser to control the total input power. The laser beam is then expanded through a beam expander (BE), which consists of an objective, a pinhole, and a lens. The expanded laser beam is then split into an object path and a reference path by a polarizing beam splitter (PBS). A half-wave plate (HFP) is placed between the NDF and BE in order to work with the PBS to adjust the intensity ratio between the object beam and reference beam. The laser beam in object path is focused by a condenser lens (CL) before passing through the imaging part. The microscopic imaging path consists of an objective (UPLFLN, 20x, NA = 0.50, Olympus) and a tube lens (U-MF2, Olympus), which create the magnified image of the sample. In the reference path, another HFP is employed to adjust the polarization direction of the reference beam. The function of the two half-wave plates is to regulate the intensity and polarization directions of both beams to produce the maximum fringe contrast for holograms. A beam splitter (BS) is then placed at the front of the camera to merge two beams so that they interfere at the CCD camera (1920×1200 pixels, 5.86 µm, PointGrey, Canada). In this study, the holograms of living cells are recorded at 1 frame per second with the size of 1024×1024 pixels.

FSS system
The FSS is generated by a custom-made flow chamber system shown in Fig. 1(b) and (c) [28]. Briefly, the system consists of a peristaltic pump (BT600-2J, LongerPump), a parallel plate flow chamber, a thermostatic bath, a pulse dampener, and a reservoir bottle with a 3-pot connector cap. The flow chamber is self-designed and custom-made and fully introduced in [29]. Briefly, the chamber consists of a commercial slide (CITOTEST, China) and a custom-designed silica cover (manufactured by ZHONGXU SHIYING COMPANY, China). The two parts are mounted into a custom-designed metal frame. All devices are connected by rubber hoses to maintain airtight. A pulse dampener is employed to retain a steady flow. The reservoir bottle is placed in a thermostatic bath, maintaining the internal temperature of the system at 37°C during the entire experiment.
The magnitude of shear stress on cells is determined by: here, µ is the viscosity of the medium measured by a rheometer and the value is 0.0107 dyn·s/cm 2 . Q is the volumetric flow rate controlled by a peristaltic pump. The width (b) and height (h) of the flow chamber are 1.5 cm and 0.06 cm, respectively. Thus, FSS in this experiment is 15 dyn/cm 2 (1.5 Pa).

Image acquisition and analysis
In this study, holograms are recorded with the speed of 1 frame per second. For every group, a single hologram is first recorded as the reference cell image before FSS. Then the pump is turned on to generate FSS on cells and holograms are recorded for a total period of 5 minutes. For phase reconstruction, the angular spectrum algorithm [19] and Zernike polynomial fitting method [30] are applied to retrieve phase maps and compensate phase aberrations. Thereafter, an offset procedure is employed to keep the average background phase value of each phase map to zero. The overall procedures are fully described in our former study [24]. All data processing including hologram reconstruction, displacement and strain map estimation are performed by the MATLAB software (MathWorks, US). All the displacement and strain data are presented by mean ± SD.

Workflow of DHM-DIC for displacement map estimation
Distinct from conventional DIC, in which intrinsic or manually added speckles on sample surface are used to search the most correlated subregion, DHM-DIC utilizes surface structural features instead of speckles as the correlating metric. The flow chart of estimating cell displacement and strain fields using DHM-DIC is exhibited in Fig. 2. Briefly, cells under FSS are first recorded by DHM, and reconstructed phase maps of cells before and after FSS are then imported into DHM-DIC as the reference and deformed images. In-plane displacement fields u(x, y) (x direction) and v(x, y) (y direction) are then obtained after DHM-DIC processing. Afterwards, strain fields ε xx , ε yy and ε xy are estimated from u and v using a 2D-SG-filter method, which will be fully described in section 3.2. Step-by-step procedures of DHM-DIC are exhibited in Fig. 2(b). S 0 (x, y) and S ′ 0 (x, y) are respectively the reference and deformed subregions with the size of (2m + 1) × (2m + 1). Vector ⃗ p noted by the red arrow is the target displacement vector to be measured and can be calculated by: where p 0 (x 0 , y 0 ) and p ′ (x ′ 0 , y ′ 0 ) are the centers of S 0 (x, y) and S ′ 0 (x, y). Since p 0 (x 0 , y 0 ) is known, determination of ⃗ p can be achieved by searching for the most correlated p ′ (x ′ 0 , y ′ 0 ) within the search area from reference image. The search areas with the size of (2M + 1) × (2M + 1) are shown in Fig. 2(b) at the top row. S 0 (x, y) is then fitted by a five-order polynomial: Here, S(x, y) is the fitted surface and shown in pseudo-3D at the right bottom in Fig. 2(b); a ij denotes the coefficients of the five-order polynomial; i and j are the order index ranging from 0 to 5. The fitting quality can be evaluated by the determination coefficient R 2 . For S(x, y) shown in Fig. 2, R 2 is 0.9934, indicating a good fitting of S 0 (x, y). R 2 is influenced by the size of subregions and in this study, m and M are chosen to be 10 and 25, resulting in a size of 21×21 for subregions and a size of 51×51 for search areas. Further considerations on window sizes are discussed in section 5.
The fitting process is then performed on every subregion in the deformed search area to find out the best matched subregion S ′ 0 . C ij (S) is defined as the correlation coefficient between each S ′ 0 and S 0 : Here, a ′ ij (S) is the polynomial coefficient of each fitted surface S in deformed search area; i and j range from [0,5]. The minimum C ij (S) indicates the best matched subregion S ′ 0 (x, y).

Strain map estimation
Strain map is estimated from u(x, y) and v(x, y) using the 2D-SG filter method in [31]. In short, considering a subregion in the displacement fields with the size of (2N + 1) × (2N + 1), the displacement fields within the subregion can be written as: where x, y ∈ [−N, N]. Here the u(x, y) and v(x, y) denote the displacement fields in x and y directions, respectively. a i and b i (i = 0, 1, 2) are the unknown polynomial coefficients. Since the relationships between strain and displacement can be described as a numerical differentiation process, strain fields subjected on the subregion can be calculated as: where ε xx and ε yy denote the strain fields in x and y directions, respectively, and ε xy denotes the shear strain. The a 1 and b 2 are then directly estimated by filtering the u(x, y) and v(x, y) with a 2D-SG filter.

Validation of DHM-DIC for displacement estimation using a simulated surface
Validation using a numerical surface h 0 (x, y) fitted from an MLO-Y4 osteocyte is first performed. An in-plane displacement with [u 0 (x, y) = −0.5µm, v 0 (x, y) = 0.3µm] is manually added to h 0 (x, y), the deformed surface h D (x, y) = h 0 (x − u 0 , y − v 0 ). The quasi-3D height map of h 0 (x, y) and the height difference map ∆h(x, y) = h D (x, y) − h 0 (x, y) are displayed in Fig. 3(a) and (b), and displacement fields u(x, y) and v(x, y) estimated by DHM-DIC are shown in Fig. 3(c) and (d).
The averaged displacements areū(x, y) = −0.49µm,v(x, y) = 0.3µm, with an absolute error of 0.005 and 0.001 µm, respectively. The estimated u and v are in good accord with the preset u 0 and v 0 . The simulation results verify that surface feature is reliable as the matching metrics for DIC searching, and also validate the feasibility of the DHM-DIC for displacement estimation.

Evaluation of whole-cell deformation under FSS
To validate the displacement and strain exhibited on cells are induced by FSS, overall morphological cell deformation is first analyzed. Temporal evolution of phase difference (PD), center of mass (COM) and projected area (PA) of osteocytes are calculated. Holographic phase images of a single cell at the 0 min, 1 min, 3 min and 5 min are shown in Fig. 4(a)-(d), and related PD maps are displayed in Fig. 4

(e)-(h). Here PD is calculated by
where φ ∆i (x, y) denotes the PD map at the moment t, φ 0 (x, y) and φ t (x, y) denote phase maps at the beginning and moment t, respectively. The entire process of cell phase and PD variation is exhibited in Visualization 1. The red dash line throughout Fig. 4(a)-(d) is shown to better reveal cell movement tendency. It can be seen that a round-like area of cell surface marked with a white arrow in Fig. 4(a) and (d) is moving towards the FSS direction over time. From the temporal variation of PD maps, it can be easily found that phase values at upper right edge of the cell increase over time while decrease at the bottom right. Variations of phase and PD maps jointly demonstrate that the overall cell is deforming along with FSS.
To quantitatively demonstrate the cell movement, we further calculate the change of cell COM and PA over time, which are shown in Fig. 4(i) and (j). Here the COM X and COM Y are the horizontal and vertical coordinates of COM and both calculated from cell projected shape. As seen in Fig. 4(i), the cell experiences a COM change of ∆COM X,max = −0.43 µm and ∆COM Y,max = 1.05 µm after 5-minute FSS, which shows the same tendency of movement with the PD change. Meanwhile, the cellular PA displayed in Fig. 4(h) nearly remains unchanged during the 5-minute FSS, indicating that the displacements induced by FSS in this work are more related to cell deformation rather than migration.

Displacement field of MLO-Y4 cells under FSS
The time-lapse holographic phase maps of MLO-Y4 cells are then imported into DHM-DIC to estimate the displacement fields. All phase maps are first low pass filtered by a 2D median filter before calculation. A video of time-lapse u(x, y) and v(x, y) during the entire 5-minute FSS is shown in Visualization 2, and Fig. 5(a)-(e) and (f)-(j) display the frames at each minute.
As seen from Fig. 5(a)-(e) and (f)-(j), displacement maps u ( x, y) and v(x, y) change over time under the influence of FSS. The u t and v t denote the average displacement of the cell, and the subscript t indicate time point. Specifically, theū(x, y) over the cell decrease by 0.375 µm while thev(x, y) increased by 0.535 µm after 5-minute FSS.
Further, displacement maps of 10 cells with and without FSS are estimated using DHM-DIC. The cells are randomly selected from three field of views. The mean u FSS , v FSS , u CON and v CON over time are shown in Fig. 5(k). As shown in Fig. 5(i), the mean total changes of displacements after 5-minute FSS are ∆u FSS (x, y) = −0.071 ± 0.363 µm and ∆v FSS (x, y) = 0.412 ± 0.220 µm,

Strain field of MLO-Y4 cells under FSS
A video of time-lapse strain fields ε xx , ε yy and ε xy when cells are under FSS is shown in Visualization 3 and frames at each minute are displayed in Fig. 6(a)-(o). The mean strains of cells with and without FSS over time are displayed in Fig. 6(p) and averaged changes of ε xx , ε yy and ε xy after FSS and without FSS are shown in Fig. 6(q). Note that in Fig. 6(p), in order to independently reveal the temporal change trends of cellular shear strain, the strains are all converted into absolute value so that mean strain over a cell can directly reveal the mechanical deformation level. It can be easily found from Fig. 6(p) that the cell strains in both directions are increasing along with FSS, while without FSS, cell strains keep at a lower level. As shown in Fig. 6(p) and (q), the absolute strain ε xx and ε yy increase by 0.029 ± 0.009 and 0.042 ± 0.018 after 5-minute FSS, respectively, and the shear strain ε xy accordingly increases by 0.026 ± 0.012. For cells without FSS (control group), the absolute strain ε xx and ε yy increase by 0.009 ± 0.012 and 0.005 ± 0.002, and the shear strain ε xy increases by 0.007 ± 0.007. All the strains present significant increase of strain when under FSS. The variations of ε xx and ε yy are in accord with our expectation, in which change level of strain along the shear direction is supposed to be larger than strain perpendicular to shear direction.

Discussion
Full-field displacement and strain map of osteocytes are quantitatively evaluated using DHM-DIC in in this study. The mean strain and maximum strain of MLO-Y4 cells along the shear direction after 5-minute 15dyn/cm 2 FSS are 0.042 and 0.060, respectively. In our previous study [25,26], a finite element model (FEM) based on the real geometrical shape of an osteocyte was built to predict the cell strain under 15dyn/cm 2 , and an average cell strain of 0.029 and a maximum cell strain of 0.051 along the shear direction was calculated. The maximum strain in this study is nearly the same as our FEM result, while the average strain seems a little bit larger. This, from our point of view, is caused by the natural fluctuations when cells are in the experimental environment. In FEM method, cells are simplified as a static model, therefore limited information can be obtained from the model owing to some unavoidable simplification during FE modeling. For example, the linear elastic simplification of the cytoplasm, cytoskeleton, and nuclei, and also the structural simplification during FEM, lost the actual viscoelastic behavior and time-dependent cell strain change under FSS. Quantifying these subtle differences between models and experimental data are our next step, and we believe DHM-DIC is potential to be a powerful tool to help us explore cell mechanics. Besides, our results are also in similar range with strains estimated by Qiu et. al. [8], in which the averaged whole-cell mean strains are 0.0520 along the shear direction and 0.0318 perpendicular to shear direction. These results validate the feasibility of DHM-DIC for estimating full-field cell displacement and strain. One of the advantages of DHM-DIC is that it is especially suitable for phase images which includes the local 3D features of the sample, and this makes DHM-DIC particularly suitable for quantitative phase imaging of living cells. Surface fitting provides correlation metrics for sample images with little natural speckle like cells. This also leads to limitations for DHM-DIC. Unlike conventional DIC method, DHM-DIC performs well for images with little speckle but sufficient 3D surface features, while it underperforms for samples with few surface features.
Three points are worth noting in the fitting procedure of DHM-DIC. The order of fitting is the first key factor to the estimation accuracy. In the present study, a five-order fitting is chosen as a result of the trading-off between accuracy and the processing duration. Secondly, R 2 of fitting is also critical to the estimation accuracy of DHM-DIC. R 2 is greatly sensitive to the size of S 0 (x, y). A larger size of S 0 (x, y) will guarantee a relatively high R 2 but lead to much longer processing time. Here we choose the size of each subregion as 21 × 21 pixel for all searching process, resulting in an overall mean R 2 of 0.9232 ± 0.0703, which validates the accuracy for surface fitting. The size of subregion also influences the error of displacement estimation. To show the relationship between subregion size and the estimation error, different m is applied during the simulation in section 4.1. Figure 7 shows the variations of processing time and the average error of displacement estimation when subregion size changes.
In Fig. 7, the horizontal axis represents the different values of m, while the left vertical axis represents the processing time and the right vertical axis represents the average estimation error.
Here, the average estimation error ε is defined as ε = √︂ (ū − u 0 ) 2 + (v − v 0 ) 2 . According to Fig. 7, it is from m = 9 when the error ε starts to keep at a lower level and the value of 10 will lead to a small error and a relatively fast processing time. Besides, the size of search area is also influential to the search process. A larger search area could include more potential matched subregions while it could also lead to longer processing time. Here the search area size is set as 51 × 51 pixels to reach a good trade-off between estimation accuracy and computing time.  7. Relationship between subregion size (2m + 1) × (2m + 1) and the estimation error and processing time.

Conclusion
In this research, a novel method named DHM-DIC is proposed for evaluation of full-field displacement and strain map for living cells. DHM-DIC innovatively employs cell surface structural parameters obtained from holographic phase image, instead of speckles, as the metrics for correlation matching. Displacement and strain maps of living osteocytes under 5-minute FSS calculated from DHM-DIC are in good agreement with our previous finite element model prediction, indicating DHM-DIC a promising tool for full-field strain measurement of living cells under mechanical loading.