ARTICLE DOI: 10.1038/s41467-017-00190-7

OPEN

Gradient light interference microscopy for 3D imaging of unlabeled specimens Tan H. Nguyen1, Mikhail E. Kandel1, Marcello Rubessa2, Matthew B. Wheeler2 & Gabriel Popescu1

Multiple scattering limits the contrast in optical imaging of thick specimens. Here, we present gradient light interference microscopy (GLIM) to extract three-dimensional information from both thin and thick unlabeled specimens. GLIM exploits a special case of low-coherence interferometry to extract phase information from the specimen, which in turn can be used to measure cell mass, volume, surface area, and their evolutions in time. Because it combines multiple intensity images that correspond to controlled phase shifts between two interfering waves, gradient light interference microscopy is capable of suppressing the incoherent background due to multiple scattering. GLIM can potentially become a valuable tool for in vitro fertilization, where contrast agents and fluorophores may impact the viability of the embryo. Since GLIM is implemented as an add-on module to an existing inverted microscope, we anticipate that it will be adopted rapidly by the biological community.

1 Department of Electrical and Computer Engineering, University of Illinois at Urbana-Champaign, Champaign, IL 61801, USA. 2 Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign, Champaign, IL 61801, USA. Tan H. Nguyen and Mikhail E. Kandel contributed equally to this work. Correspondence and requests for materials should be addressed to G.P. (email: [email protected])

NATURE COMMUNICATIONS | 8: 210

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

1

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-00190-7

I

microscopy17 and differential interference contrast (DIC) microscopy18. The contrast in these methods is generated by visualizing the modifications of the wavefront when light propagates through the sample. Unfortunately, both PC and DIC are qualitative—that is to say, they do not measure the wavefront deformation quantitatively, and the image recorded on the detector is often substantially different from the scattering potential of the object. This deformation is characterized by a spatially dependent phase shift, defined as ϕ(r) = (2π/λo) h(r)Δn(r), where λo is the central wavelength of the illumination, h(r) and Δn(r) are the sample thickness and refractive index difference, both evaluated at the transverse coordinate r, respectively. Quantitative phase imaging (QPI) is an approach focused precisely on quantifying this phase shift19. Along these lines, Cogswell et al.20 proposed DIC with geometrically-induced phase shifting, applied for two-dimensional (2D) imaging. DIC with two orthogonal shear directions has been used to obtain 2D quantitative phase images21, 22, whereas Mehta et al.23 reported a partially coherent model for DIC 2D imaging. Shribak et al.24 used liquid crystal modulators to change the polarization directions and phase shift modulation. This setup allows them to record the 2D phase-gradient information at two orthogonal directions and reconstruct the optical phase front. QPI has recently gained significant scientific interest, especially in the biomedical field19, thanks to several advancements. For example, common-path interferometry replaced traditional interferometry for better stability and sensitivity25–28. Low temporal coherence illumination methods significantly improve image resolution when suppressing speckles29–32. An interesting direction of study is using QPI to extract scattering information from extremely weakly scattering objects33. This approach is referred to as Fourier transform light scattering (FTLS), a spatial analog to Fourier transform (infrared) spectroscopy34. The idea is that the knowledge of amplitude and phase of an image field allows us to numerically propagate that field to any plane, including the far field, where angular scattering

t has become increasingly clear that understanding morphogenesis and disease requires three-dimensional (3D) tissue cultures and models1. Effective 3D imaging techniques, capable of reporting on subcellular as well as multicellular scales, in a time-resolved manner, are crucial for achieving this goal2. Although the light microscope has been the main tool of investigation in biomedicine for four centuries, the current requirements for 3D imaging pose new, difficult challenges. Owing to their insignificant absorption in the visible spectrum, most living cells exhibit very low contrast under light microscopy. As a result, fluorescence microscopy has become the main tool of investigation in cell biology3. Due to the significant progress in designing fluorescence tags, structures in the cell can now be imaged with high specificity. More recently, super-resolution microscopy methods based on fluorescence have opened new directions of investigation, toward nanoscale subcellular structure4. However, fluorescence imaging is subject to several limitations. Absorption of the excitation light may cause the fluorophore to irreversibly alter its molecular structure and stop fluorescing. This process, known as photobleaching, limits the time interval over which continuous imaging can be performed5. The excitation light is typically toxic to cells, a phenomenon referred to as phototoxicity6. Overcoming these limitations becomes extremely challenging7, 8 when imaging thick objects over an extended period of time9, 10 as acquiring data over the time and the axial dimension increases exposure of the specimen to the excitation light, lowering its viability. Confocal11 and two-photon fluorescence microscopy12 have been the mainstay tools of imaging thick 3D specimens. Although these methods can provide excellent sectioning through tissue, due to the focused, short wavelength excitation, the amount of power required may be harmful. Thus, recent advances in light sheet microscopy were dedicated specifically to reducing phototoxicity and photobleaching13–16. Label-free microscopy provides an alternative solution to overcoming these limitations, albeit at the expense of molecular specificity. Two classical methods are phase contrast (PC)

b

a

GLIM image (r⊥)

0.5 π

0π

GLIM optical setup

c

Intensity images

2.5 1.2 0.0 –1.2 –2.5

10 μm

∇x (red μm–1)

1.5 π

1.0 π

10 μm

P1

L1

d

In (count) 10E+3 4.0 0.5

L2

–1.2

SLM

Camera

e

2.2

 (red)

5 μm

5 4

Phase (rad)

6.5E+3

Integrated Ground truth

3 2 1 0 –1 –2 –15 –10 –5

0 5 x (μm)

10

15

Fig. 1 Optical setup and working principle of GLIM. a GLIM optical setup. The GLIM module is designed as an add-on module connected to the output port of an inverted microscope. This module shifts the phase of one polarization component using an SLM while keeping the other unmodified. Interference patterns generated by these two components are recorded and transferred to a computer for phase-gradient extraction. b Four frames are acquired by the GLIM module, one for each phase shift applied by the SLM. Using these images, we obtain the phase-gradient map and integrate it along the direction of the shift to recover the quantitative phase map. c Extracted quantitative-gradient map of two 3 µm polystyrene beads immersed in oil. d Integrated phase map of a 4.5 µm polystyrene microbeads at NAcon = 0.09. e Cross-sections of the reconstructed phase and the computed ground truth (black dashed curve) taking into account blurring due to diffraction (Eq. (3)) 2

NATURE COMMUNICATIONS | 8: 210

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-00190-7

a

d

b

Norm. mass

t = 0 min.

1.19

Growth rate G(t ) ∝

1.12 DT = 36 h 1.05 0.98 0

100 μm

e c

(rad μm–1) 0.08 0.04 0 –0.04

1 mm

–0.08

10 μm

200 400 Time (min.)

600

Transport relation

1.E–03  (rad s–1)

∇

2t/DT

5.E–04 10

(k⊥) = k⊥  = 58.2 nm min–1

2π / k⊥ (μm)

40

Fig. 2 Time-lapse GLIM imaging of a HeLa cell culture. a The image sequence is composed of 16 × 20 mosaic tiles, each of 280 × 277 μm2 area with an acquisition performed every 16 min. b Image of each mosaic tile of a small region in a. c A zoomed in region of b. The resulting image is well suited to observe cellular morphology, in particular, cells nucleoli are well represented (white arrows in c). The same colorbar applies for all figures from a–c. From the GLIM image, it is possible to measure the mass doubling time, DT. d Measured normalized mass profile (black dotted) and an exponentially fitted one (red). e The slope measurement Γ(k⊥) of the temporal autocorrelation function as a function of k⊥. Red profile: fitting curve of Γ(k⊥) using the expression in Eq. (6). We estimated variance of the advection velocities corresponding to active transport, Δν ≈ 58.2 nm min−1. The contribution of diffusion motion is negligible from the fit

measurements are typically performed. For weakly scattering objects such as live cells, it is much more signal-effective to perform the measurement at the image plane, where all scattering angles overlap at each point, rather than measuring angle-byangle in the far field. As a result, QPI can be used to solve inverse scattering problems and extract the 3D structure of inhomogeneous objects35. Three-dimensional information of the specimen is accessible by measuring the phase across multiple angles of the illumination or axial specimen positions36–39. An equivalent approach is fixing the illumination direction while rotating the sample to obtain phase maps from different viewing angles. For example, see Merola et al.40 Lens-free holography method is also used to obtained tomographic information on a chip by combining with multiple illumination angles41. However, imaging optically thick, multiple scattering specimens is still challenging for any optical method, including QPI. The fundamental obstacle is that multiple scattering generates an incoherent background, which ultimately degrades the image contrast. An imaging method dedicated to imaging these thick specimens must include a mechanism to subdue the multiple scattering backgrounds and exhibit strong sectioning to suppress the out-of-focus light. To overcome these challenges, here we introduce a new QPI method, referred to as gradient light interference microscopy (GLIM). GLIM combines DIC microscopy with low-coherence interferometry and holography. In GLIM, the two interfering fields are identical except for a small transverse spatial shift. This geometry ensures that the two fields suffer equal degradation due to multiple scattering. By accurately controlling the phase shift between the two waves, we acquire multiple intensity images, which have the same incoherent background, but different coherent contributions. As a result, GLIM rejects much of the multiple scattering contributions and yield high contrast of thick objects. Furthermore, the illumination condenser aperture is fully open, which lands GLIM very strong optical sectioning. GLIM can provide tomographic imaging of both thin samples, for example, single cells, and thick specimens, such as multicellular systems. Below, we present the principle of GLIM operation, validation results on test samples, and timeresolved tomography of cells in culture, as well as embryo development. NATURE COMMUNICATIONS | 8: 210

Results GLIM principle. GLIM is an add-on module to a commercial DIC microscope as shown in Fig. 1a. Via a Wollaston prism, a typical DIC microscope generates two replicas of the image field, cross-polarized, shifted transversely by a distance smaller than the diffraction spot. We removed the analyzer that normally renders the two polarizations parallel in DIC and, instead, let the fields enter the GLIM module. These fields are spatially Fourier transformed by the lens L1 at its back focal plane. A spatial light modulator (SLM), placed at this plane with its active axis aligned to the polarization direction of one field, retards its phase by ϕn = nπ/2 with n = 0, 1, 2, 3, and leaves the other field unmodified. Both fields are Fourier-transformed again by lens L2 to generate the image at the camera plane. A linear polarizer, P1, is aligned at 45° with respect to both polarizations to render them parallel. The resulting field at the detector is a coherent superposition of these two fields, namely, Un ðrÞ ¼ U ðrÞ þ U ðr þ δrÞeiϕn ;

ð1Þ

where δr ¼ δx^ x is the spatial offset between the two fields and U is the image field. The intensity at each phase shift, In(r) = |Un(r)|2, can be written as In ðrÞ ¼ I ðrÞ þ I ðr þ δrÞ þ 2jγ ðr; δrÞjcos½ϕðr þ δrÞ  ϕðrÞ þ ϕn ; ð2Þ where I(r) and ϕ(r) are, respectively, the intensity and phase of the image field, and γ is the mutual intensity or the temporal cross-correlation function at zero-delay between these two fields, γ(r,δr) = 〈U*(r)U(r + δr)〉t. The quantity ϕn = nπ/2 is the modulated phase offset between the two fields, externally controlled by the SLM. From the four intensity images, In, with n = 0, 1, 2, 3 (Fig. 1b), we are able to solve for I(r), |γ(r,δr)|, and Δϕ(r) = ϕ(r + δr)−ϕ(r). These data render quantitatively the gradient of the phase along the direction of the shift (Fig. 1c), ∇xϕ(r) ≈ Δϕ(r)/δx. Details on the optical setup, procedures for extracting the phase gradient and estimating δx can be found in the “Methods” section. Before running the experiments, the SLM needs a one-time calibration to ensure proper phase modulation. The calibration

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

3

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-00190-7

∇ (rad μm–1)

GLIM

a

x–y

b

0.50

x–y

weiner,bf

GLIM + filtering

c

(dB) 0

d

0.25

–25

0.00

–50

–0.25 15 μm

25 μm

–0.50

–75

25 μm x–y

f

e

Mass (pg)

81.2 μm

20

h

i

0 0.E+00 2.E+04 4.E+04 6.E+04

μm

j

k

19

0

154 min

242 min

264 min

Surface area (μm2)

0 33 μm

40

110 min

132 min

–100

60

66 min

m

0 min

x–z

80

l

g

15 μm

Volume (μm3) 3.E+04 2.E+04 2.E+04 1.E+04 5.E+03 0.E+00 0.E+00 2.E+04 4.E+04 6.E+04 Volume (μm3)

Fig. 3 GLIM time-lapsed tomography of a HeLa cell culture. a x−y cross-section and b x−z cross-section of the GLIM phase-gradient measurement. Dashed line denotes the locations of the cross-section. c x−y cross-section and d x−z cross-section after spatial filtering. e–k Rendering of cell tomography results at different time, as indicated. Here, the cell nuclei are segmented and displayed in orange. The cell membranes are in green. l Mass vs. volume and m area vs. volume measurements extracted from the GLIM tomography data. Data points with the same color belong to the same cell but were acquired a different time points

procedure and pixel-to-pixel variation of the SLM are described in details in Supplementary Notes 1, 2, and Supplementary Fig. 1. QPI using GLIM. To demonstrate the capability of GLIM to extract quantitatively the phase gradient, we imaged 4.5 ± 5% µm polystyrene micro-beads (Polysciences Inc.), with a refractive index value of 1.59 at the central wavelength. The beads are immersed in immersion oil (Zeiss Inc.) with a refractive index value of 1.518 to generate a total phase shift of 3.87 radians. Figure 1c shows the measured phase gradient at NAcon = 0.09, where the subscript con stands for condenser. Given the phase gradient, ∇xϕ, one can integrate along the gradient direction to get phase value, ϕ(r), using Zx ϕðx; yÞ ¼

½∇x ϕðx0 ; yÞ dx0 þ ϕð0; yÞ;

ð3Þ

0

where ϕ(0,y) is the initial value, which can be obtained with some prior knowledge of the specimen. For example, if (0,y) is a background location, the phase ϕ(0,y) should be set to 0 radians. Figure 1d shows the quantitative phase map, ϕ(r), and Fig. 1e displays a line profile through the center of the bead. Note that our integration result matches well with the expected ground truth, where optical diffraction is taken into account. GLIM imaging of cell cultures. Due to the low phototoxicity, absence of photobleaching, and easy sample preparation, transmitted light modalities appear to be ideal for studying cell growth and proliferation42. Yet, such assays are most frequently conducted with the aid of labels. Although specificity granted by external markers is crucial for certain applications, quantifying cell growth over longer timescales remains a challenge43. It has 4

been known for some time that indicators of cell proliferation do not have equal growth44. More recently, new approaches have been demonstrated using vibrating hollow cantilevers to weigh cell passing through45. This method is limited to non-adherent cells. To expand the principle of this measurement to adherent cells, a method based on vibrating pedestals was demonstrated, at the expense of mass sensitivity46. We show that GLIM is able to quantify the growth and proliferation of large populations of adherent cells over extended periods of time. Specifically, we can characterize the culture by extracting parameters such as single cell mass, volume, surface area, while simultaneously measuring the intracellular transport on timescales associated with the cell cycle. Supplementary Note 3 provides details on 2D image formation in GLIM. Figure 2a shows scanning GLIM data of HeLa cells in culture over a 4.48 × 5.54 mm2 field of view. This image is a mosaic consisting of 16 × 20 individual frames, imaged by a 40×/0.75 NA objective and a condenser aperture adjusted to NAcon = 0.32. Figure 2b and c shows magnified views at different scales for a region denoted by the white box in Fig. 2a. The acquisition took ~3 min for each of the 16 × 20 mosaics, which was assembled into a time-lapse sequence following the procedure outlined in the “Methods” section, Supplementary Note 4, and Supplementary Figs. 3 and 4. The GLIM image (Fig. 2c) provides a quantitative phase-gradient map at the spatial resolution of the objective, clearly showing fine structures such as nucleoli pointed by white arrows. We acquired 38 such large fields of view, over a 10 h time interval (see the Supplementary Movie 1 for a time-lapse sequence). To measure growth rates, the phase values are obtained by integrating the phase gradient at each frame in the time-lapse sequence using Eq. (3). These phase values are used to calculate the cell dry mass, using the linear relationship between the optical

NATURE COMMUNICATIONS | 8: 210

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-00190-7

path-length map of a cell and its dry mass density, ρ, see Barrer47. The black dotted profile in Fig. 2d shows the normalized dry mass to that at time t = 0. Fitting this profile to the function 2t/DT, where DT is the mass doubling time constant, we found that DT ≈ 36 h. Interestingly, this time is approximately 50% longer than the typical cell count-doubling time as cells can divide without doubling in mass48. The fitted profile is shown in red in Fig. 2d. To investigate the physics of cellular mass transport, we use the dispersion phase spectroscopy (DPS) method49. This approach is powerful in extracting spatiotemporal fluctuation information from time-lapse sequence of phase maps, as it requires no manual tracing, making it well suited for fully automated applications50. In DPS, one computes the dispersion relation of the dynamic system, connecting the spatial and temporal frequencies. The behavior of this dispersion curve informs on the nature of the transport (such as, comparisons of diffusion vs. deterministic dominated transport) and numerical fits yield the diffusion coefficient, and the width of the velocity distribution. Let us start by describing the dry mass-density fluctuation, ρ(r⊥,t), via an advection-diffusion equation, namely, D∇2 ρðr? ; t Þ  v:∇ρðr? ; t Þ  ∂ρðr? ; t Þ=∂t ¼ 0;

ð4Þ

where r⊥ = (x,y) is the 2D coordinate vector, v is the advection velocity vector, and D is the average diffusion coefficient. Using this equation, we obtain the temporal autocorrelation function evaluated at each spatial frequency, k⊥, and temporal delay, τ, defined as ρðk? ; t Þ~ ρ ðk ? ; t þ τÞit =h~ρðk? ; t Þ~ρ ðk? ; t Þit . Here, g ðk? ; τÞ ¼ h~ ~ρðk? ; t Þ ¼ Fk? ½ρðr? ; t Þ is the 2D spatial Fourier transform of the dry mass density. We calculated g(k⊥,τ) directly from the phase gradient ∇ϕ instead of the integrated phase ϕ (Supplementary Note 5), namely     g ðk? ; τÞ ¼ expðivo :k? τÞ exp  Δvk? þ Dk2? τ ; ð5Þ where vo is the mean and Δv the standard deviation of the velocity distribution. At each transverse spatial frequency, k⊥, one can fit the measurement of g(k⊥,τ) using Eq. (5) to estimate Δv and the diffusion coefficient, D. We found that vo is negligible for the duration of the experiment or that no dominant velocity vector, vo,was detected in our data. The decay rate of g(k⊥,τ) at each spatial mode k⊥ satisfies Γðk ? Þ ¼ Δvk? þ Dk2? :

ð6Þ

The first term in Eq. (6) dictates the active transport, whereas the remaining term characterizes diffusion. The black dot profile in Fig. 2e is our measured Γ(k⊥). Fitting this profile to Eq. (6), we found that active transport dominated on cellular scales (10–50 μm), with a spread in velocities of Δν ≈ 58 nm min−1. Fitted profiles are shown in red. Tomography of single cells using GLIM. As a result of the high numerical aperture of the illumination, GLIM has strong sectioning capabilities, which yields tomographic imaging of both thin and thick samples. We apply GLIM tomography to a 30% confluence HeLa cell culture over 7.7 h. Seven fields of view (FOVs) were imaged using a 63×/1.4 NA objective with a spatial sampling rate of 10.8 pixels µm−1. Each FOV was scanned every 22 min. For each time point, the sample is scanned over a total depth of 28 µm with a step size of Δz = 0.07 µm. Figures 3a and b show the x−y and x−z cross-sections of the GLIM measurement, namely the quantitative phase gradient ∇ϕ. To remove the background due to weak sectioning at small scattering angles, we perform a spatial high-pass filtering operation, as described in NATURE COMMUNICATIONS | 8: 210

the “Methods” section. Figures 3c and d show the corresponding x−y and x−z cross-sections of the filtering, with the yellow arrows pointing to the locations of the nucleus. Clearly, the x−z cross-section of the tomograms shows significant improvements in depth sectioning. Compared to the phase-gradient image, ∇xϕ, this cross-section has no diffraction streaks or shadow artifacts, while preserving clear cell boundaries. Figures 3e–k show the GLIM tomograms obtained via filtering, at seven different time points. The cell nuclei were segmented and shown in orange, whereas the cell membranes displayed in green using isosurface rendering. The rendered images clearly show how the 3D shape of the cell changes over time. It can be further seen that during mitosis (the 110 and 264 min. frames), the cells assumed a spherical shape (pointed by yellow arrows in Fig. 3g). Also, at the 110-min point, while forming a mitotic sphere, the cells appear to leave behind biomass (white arrow in Fig. 3g) that is adherent to the substrate, consistent with previous observations51. Rendered images for the whole time series can be found in Supplementary Movie 2. Thanks to the depth sectioning of GLIM, we can automatically segment different cells in the volume of interest (Supplementary Note 6; Supplementary Fig. 5) to compute several parameters of each cell and study their temporal evolution. An example of GLIM tomograms over a full cell cycle can be found in Supplementary Fig. 6. Figure 3l shows the dry mass (Supplementary Note 7; Supplementary Fig. 7) vs. volume for several different cells during a 21-h window. Each point in these plots corresponds to one cell at one time point. These results (also those from Supplementary Note 8) show that, for the most part, the points align along a straight line we found that the points deviating from this line correspond to cells going through mitosis. Meanwhile, the surface area vs. volume relation shown in Fig. 3m is essentially linear with slightly different slopes for different cells over the whole cell cycle. GLIM investigation of embryo viability. The Centers for Disease Control and Prevention (CDC) report from 2014 shows that 208,768 Assisted Reproduction Technology (ART) cycles were performed with 57,332 live births52. As the numbers indicate, the percentage of live births from these procedures is still rather low. One reason is the lack of objective and accurate evaluation of embryo quality and viability before transfer. Morphological assessment is currently the main method used to determine embryo viability during in vitro fertilization (IVF) cycles. However, studies have shown that the predictive power of the typical day 2 and 3 assessment of morphological parameters has remained low53–55. Various noninvasive analytical tools have recently been used for noninvasive prediction of embryonic potential56–60. One such tool has been the development of quantitative techniques for the non-invasive assessment of embryo metabolism, and its value as a predictor of embryo viability is the subject of ongoing investigations61. But currently, visual observation remains the most used and reliable method. With the improvement of microscopy, it is possible to follow embryo development in real time, and it has been established that morphokinetic parameters can be used to select embryos for higher viability62. One of the most important microscopy techniques is transmission electronic microscopy (TEM), which is considered by many the main tool for intracellular evaluation. The main problem with using TEM for embryo evaluation is that the sample preparation kills the embryo63. Therefore, although this type of microscopy can be considered an important tool for research, it has little value for routine IVF procedures. Other techniques used to evaluate the embryo quality include confocal microscopy64 and two-photon imaging65. Using these

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

5

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-00190-7

a

b

Time-lapse ⏐(r⊥, t )⏐ac 23 h

0h

46 h

69 h

1

0.12

0.9 DIM

30 μm

0.06 0.00

115 h

–0.06 92 h

138 h

162 h

0.8 DIM(t )=ae–t/τ τ = 6.2 h

0.7

–0.12

⏐⏐ac

0.6 115

125

145

135

Time (h) 1 0h

23 h

46 h

69 h

DIM(t )=ae–t/τ

0.9 DIM

0.12 30 μm

0.06

τ = 12.4 h 0.8

0.00 –0.06

92 h

115 h

138 h

162 h

0.7

–0.12

0.6 115

⏐⏐ac

125

135

145

155

165

Time (h)

Fig. 4 GLIM assay of embryo viability. a Bovine embryos were observed under regulated atmospheric conditions, with a 40×/0.75 NA objective and a fully open condenser. As outlined in “Methods”, |γ(r⊥,t)|ac is a high-passed version of |γ|, such that |γ(k⊥,t)|ac = |γ(k⊥,t)|F(k⊥) with F(k⊥) the filter function defined in Eq. (11). b The internal structure dynamics of embryos was quantified by using DIM, based on the time derivative of the stack. DIM shows a continuous decrease over several hours prior to the embryo’s exhaustion. The raw DIM profiles (black dotted) and their fitted profiles (red)

a

∇ (rad μm–1)

b GLIM

4.5 3.0

f

e

1.5 0.0 –1.5

x–y

c

259 μm

2 cells 30 μm

y–z weiner,bf (dB)

d GLIM + filtering

0

g

–25

4 cells

–50

h

–75

5

7.

27

–100

μm

x–y

y–z

275

μm

5 cells

Fig. 5 GLIM imaging for thick embryos. a x−y cross-section and b y−z cross-section of raw GLIM images in a bovine embryo. c x−y and d y−z cross-sections after spatial filtering. e A rendered embryo cut through the center to reveal internal structures. f–h Bovine embryos at different points in their development stages, as indicated

fluorescence methods, the sample must be tagged to be evaluated, which can be detrimental to embryo survival. Due to its sectioning capabilities, GLIM can be used to perform tomography on optically thick specimens such as embryos. For this demonstration, we used bovine embryos, prepared as described in the Supplementary Methods. In a single experiment, we imaged 60 bovine embryos, starting at 12 h after fertilization, sampling every 30 min, over a 7-day period, using a 40×/0.75 NA objective. The embryo thicknesses are within 250–300 µm. Supplementary Movies 3 and 4 illustrate the high contrast that GLIM yields even in these challenging, multiple scattering samples. For example, the lipid droplets, prominent in bovine embryos can be clearly identified. Their contrast switches from dark to bright as they pass through the focus. Our results show that the embryo internal dynamics changes completely when the embryo dies. Specifically, the internal mass transport halts almost 6

entirely, which suggest either a large increase in viscosity of the material or that the dynamic transport is mostly due to molecular motors, which stop in dead cells. GLIM can be a valuable tool for IVF because it provides an intrinsic marker to predict viability in advance. Toward this goal, we developed a dynamic index marker (DIM), based on the GLIM data. This metric is computed from the phase difference Δϕ(r⊥,t) and the mutual intensity γ(r⊥, δr, t) at each time point t. To measure morphological changes, we compute the time derivative of γ, that is γ′t(r⊥) = dγ(r⊥, δr, t)/dt. On the basis of γ′t, we calculate the spatial cumulative distribution n function (CDF) ofo the time derivative images Ft ðxÞ ¼ P jγ 0 t ðr? Þjr? 2FOV < x ; which is a probability that the amplitude of the time derivative is less than a value x. The cut-off difference (distance in x) at 10 % and 90 % of the CDF for each time point t is defined as Dt ¼ arg minx1 ½F ðx1 Þ>0:9  arg minx2 ½F ðx2 Þ>0:1: Finally, we define the DIM as the relative ratio between Dt and its

NATURE COMMUNICATIONS | 8: 210

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-00190-7

maximum over the imaging duration, maxt (Dt), namely DIMðt Þ ¼ ½Dt =maxt1 ðDt1 Þ:

ð7Þ

Intuitively, during periods of inactivity, the spatial distribution of γ ′i across the embryo is uniform, the histogram is narrow compared to periods of higher activity. We found that between the point of apparent normal, dynamic behavior and the one with a complete lack of dynamic behavior, there exists a continuous process that lasts many hours (Fig. 4a). This process is well captured by the intrinsic DIM quantity, as shown in Fig. 4b. Here, the black dotted profiles are the raw measurement of DIM as a function of time t for different embryos. We found that DIM decreases continuously over several hours until it reaches the point where the embryo motion is suppressed. Further, the DIM (t) profiles are well described by exponentially decay functions (red profiles in Fig. 4b) with time constants of 6.2 h and 12.4 h, respectively. Therefore, we anticipate that this intrinsic dynamic marker can potentially hold valuable viability prediction capability, beyond the current, morphology-based assays. Embryo tomography with GLIM. We obtained 3D GLIM stacks of bovine embryos at different development stages. We used a 63×/1.4 NA oil immersive objective at a transverse sampling rate of 10.8 pixel µm−1. The condenser aperture was fully opened to NAcon = 0.55 to maximize the depth sectioning and spatial resolution. The embryos were scanned in the axial dimension over an interval of (−120 µm, 120 µm) with a step of Δz = 0.05 µm. Figure 5a and b shows the x−y and x−z crosssections of the raw phase gradient, ∇ϕ. As a side note, in comparison with other QPI method, such as the Spatial Light Interference Microscopy (SLIM)29, GLIM is superior when imaging optically thick samples like embryos (Supplementary Note 9; Supplementary Fig. 8). The corresponding cross-sections of the GLIM tomogram after filtering are shown in Fig. 5c and d. More details on this filtering step can be found in the “Methods” section. The GLIM tomography, reveals various structures of the embryos, including their membranes, internal cells, gaps between the membrane of the cells, and their internal content, lipid droplets in each cell, as indicated in Fig. 5e. The x−z crosssections further show the contact between the embryo and the underlying glass substrate (Fig. 5d), along with debris on the substrate. Figure 5f–h shows the rendering results of three different embryos, consisting of two cells, four cells, and five cells, as indicated. One can see clearly how different cells of the embryo stack with respect to each other in 3D. The membranes of the embryos are manually segmented and displayed as transparent surfaces. Rendering videos of these embryos can be found in the Supplementary Movies 5–8. Discussion In summary, we introduced GLIM, as a new QPI method, for 3D imaging of unlabeled specimens. GLIM has all the benefits of common-path white-light methods including nanometer pathlength stability, speckle-free, and diffraction-limited resolution. At the smallest condenser aperture, GLIM gives exact values of the quantitative phase for thin samples. At the largest condenser aperture, GLIM can be used as a tomography method, allowing us to obtain time-lapse 3D information of thick samples. We demonstrated the success of GLIM on various samples, from beads, HeLa cells, to bovine embryos. We believe that this method will set an excellent foundation for other research projects and applications. As a label-free method, GLIM can be applied to imaging live cells and thick samples nondestructively over broad temporal and NATURE COMMUNICATIONS | 8: 210

spatial scales. This technique is not limited by photobleaching and phototoxicity commonly associated with fluorescence microscopy. Also, it provides excellent optical sectioning and obtains 3D information from unlabeled specimens. However, similar to other label-free images, GLIM lacks specificity. Therefore, we envision that GLIM and fluorescence techniques will co-exist and corroborate the advantages of specificity and noninvasiveness. This is completely feasible since GLIM operates on the same optical path as the fluorescence channels, allowing a seamless transition between the two modalities. Methods GLIM optical setup. The GLIM add-on module is mounted to the output camera port of a conventional DIC microscope. The measurements used in Fig. 1b were acquired with an Olympus IX70 microscope equipped with a 20×/0.65 NA objective. Subsequent measurements were conducted using an Axio Observer Z1 microscope with incubation system (Zeiss) under ×40/0.75 NA (420361-9910-000) and 63×/1.4 NA (420781-9910-0000) objectives. The GLIM module contains a polarizer before the camera and the SLM (Meadowlark) is positioned to match the DIC shear angle (45°). To ensure optimal modulation, the SLM was illuminated with filtered white light using a fluorescent emission filter with a central wavelength of λo = 624 nm and a bandwidth of Δλ = 43 nm. Phase-gradient extraction from intensity images. The intensity image at modulation ϕm is given by In ðrÞ ¼ I ðrÞ þ I ðr þ δrÞ þ 2jγ ðr; δrÞj cos½ΔϕðrÞ þ ϕn ;

ð8Þ

where Δϕ = ϕ(r + δr)−ϕ(r) ≈ ∇(ϕ)δr, the phase difference of interest, and ∇ϕ the gradient of the phase in the direction of the shift. The spatial shift δr is the transverse displacement introduced by the DIC prism, estimated experimentally from measurements of the test samples. The quantity γ(r, δr) is the mutual intensity, in other words, the temporal cross-correlation function between these two fields, evaluated at zero delay, γ(r, δr) = 〈U*(r)U(r + δr)〉t. Combining the four intensity frames, we obtain the phase gradient as ∇ϕðrÞ ¼ argf½I4 ðrÞ  I2 ðrÞ; ½I3 ðrÞ  I1 ðrÞg=δr:

ð9Þ

DIC shear estimation. The lateral offset between the two DIC beams, δr, relates the phase difference image, Δϕ, linearly to a quantitative phase gradient, ∇ϕ, via the following relation Δϕ = ∇ϕδr. Although this parameter is known to the microscope manufacturer (Zeiss, Olympus, etc.), to the best of our knowledge, it is not publicly listed. To estimate the spacing between the two beams, we matched the associated phase shift to a known calibration sample. In our procedure, we acquired a fine tomographic stack of a small object (300 nm bead), and performed line integration in the direction of the DIC gradient. As expected, the phase was always maximized at the plane of best focus. The peak of the integrated phase corresponds to the theoretical maximum phase shift due to the control structure, namely φ* = (2πd Δn/λo), where Δn is the refractive index mismatch between the sample and the background, λo is the central wavelength, and d is the diameter of the bead. However, dimension of R x given a Rpixel Pp (µm), we also have x another relationship φ ¼ xo ∇x φdx  xo Δφdx=δr  p kk¼ko Δφ½k=δr; where xo, x* are locations of the background and somewhere inside the control structure. These locations correspond to pixel indices of ko and k*, respectively. Combining these two relations yields pλo δr ¼

k P

Δφ½k

k¼ko

2πdΔn

:

ð10Þ

Using this formula, we estimated that the DIC prism shifts are 175 nm for 63×/ 1.4 NA and 345 nm for 40×/ 0.75 NA sliders. These values are in good agreement with reported literature66. 2D real-time interferometric reconstruction. The 2D image formation model in GLIM is shown in Supplementary Note 3, where we relate the measured phase difference Δϕ with the sample transmission, T. To fully automate the data acquisition for Δϕ, we developed a software platform capable of mechanical automation and real-time phase retrieval. Our image acquisition platform is designed to overlap the GLIM computation with the operation of the camera, SLM, and microscope. The software is developed in C++ using the Qt framework. The real-time reconstruction runs on three threads with the first thread responsible for triggering new camera frames and modulating the SLM. The second thread receives incoming images and transfers them to the graphics card. The third thread is used to display the GUI and render the resulting phase maps (Supplementary Fig. 2a).

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

7

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-00190-7

As we decouple the triggering thread from the data transfer and computation threads, variability in camera transfer rates does not slow down the acquisition. To remove the phase slant found in commercial DIC microscopes (see ref. 42), we perform Fourier bandpass filtering. Here, we construct a filter function that does bandpass filtering as         F ðk? Þ ¼ 1  exp k2? = 2ρ2hi exp k2? = 2ρ2lo ; ð11Þ where ρhi and ρhi define the bandwidth of the high-pass and low-pass filtering operations, respectively. Their values also depend on the system magnification. This operation eliminates the slow-varying oscillation in the GLIM images. The results of this operation are shown in the Supplementary Fig. 2b. To improve the SNR of the reconstruction, we match the mean intensities of the four frames by proportionally adjusting the exposure time. For example, the exposure time is eight times longer for the extinction frame (π modulation) compared to maximum brightness frame (0π modulation). The longer exposure is later compensated numerically in extracting the phase (Supplementary Fig. 2c). This operation results in an increased signal-to-noise ratio (SNR), thanks to a reduction in the phase noise deviation (Supplementary Fig. 2c). Our GLIM system operates at 10 phase images per second with a rendering rate at 40 frames per second. As the computational portion is overlapped with acquisition, the rate-limiting factor in our system is the exposure time. Thus, longer exposure can be replaced by illumination with a brighter source. After acquiring individual images from different FOVs, we combine them together, forming a large mosaic to study large-scale dynamics. See Supplementary Note 4 and Supplementary Fig. 3 for implementation details on the image alignment and registration algorithm. Filtering method to improve depth sectioning of GLIM. To improve the optical sectioning, and push GLIM into a 3D imaging method, we removed the low-frequency (out-of-focus) components from our data using a high-pass filter. Steps of our methods are summarized in the Supplementary Fig. 9. First, we removed the DIC shading artifact using Wiener deconvolution67. The 3D point spread function of the system is given as hðrÞ ¼ Im½ðμi g  ÞðrÞ  ðμi g  Þðr  ^ xδxÞ; (Supplementary Note 10; Supplementary Fig. 10), the transfer function is ~ðkÞ ¼ 2i sinðkx δx=2ÞF fImðμ g  ÞgðkÞ. As a side note, a measurement of the point h i spread function using a microbead can be found in the Supplementary Fig. 11. The Wiener deconvolution result of the susceptibility can be obtained in the frequency domain as ~χweiner ðkÞ ¼

2i sinðkx δx=2Þδx F ½∇x ϕðkÞ; β2o ½4sin2 ðkx δx=2Þ þ ε

ð12Þ

where ε is a small number, set to be 10−4 to avoid amplifying frequency components with small SNRs. To further improve the axial resolution, it is necessary to significantly suppress the low-frequency components in χweiner(r). We achieve this by applying high-pass filtering in the x−y domain for each recorded z-image. In each dimension (x and y), a convolution with a finite-length impulse response (FIR), chosen as hhp (x) = (0.25, −0.25, 0, −0.25, 0.25), is applied. The result of this high-pass filtering, χweiner, hp(r), (Supplementary Fig. 9b) has most of the small transverse frequencies suppressed and, as a result, yields very good depth sectioning. Note that this high-pass filtering step can be combined with the Wiener deconvolution step since both are linear operators. Also, there is no need to perform any z-processing in our proposed method. This allows the processing to interlace with image acquisition. After filtering, we applied a log-compression transform to increase the contrast of the retained high-frequency components in the output image and normalized so that the largest signal is 0 dB. To reject the background noise, we keep signals with amplitude larger than −100 dB. Finally, to smooth the image and remove high-frequency oscillations in the image, we further apply bilateral filtering68 on the transformed results to obtain χweiner, bf(r⊥,z) (Supplementary Fig. 9c). There are two passes of 2D bilateral filtering. In the first pass, bilateral filtering is applied to each 2D Wiener deconvolved image χweiner (r⊥,z) for each value of z. The second pass applies bilateral filtering on the stacked result of the first pass for each 2D image of the same lateral coordinate x. Owing to the similarity, we describe here the first one. For each value z, the bilateral results χweiner,bf (r⊥,z) is obtained from χweiner (r⊥,z) using χ weiner;bf ðr? ; z Þ ¼   RR 2 0 d r ? cr ðr? ; r0 ? Þcs χ weiner;bf ðr? ; z Þ; χ weiner;bf ðr0 ? ; z Þ χ weiner;bf ðr0 ? ; z Þ   RR 2 ; d r00 ? cr ðr? ; r00 ? Þcs χ weiner;bf ðr? ; z Þ; χ weiner;bf ðr00 ? ; zÞ

ð13Þ

where the radially symmetric Gaussian  functions were used  for the closeness function cr ðr? ; r0 ? Þ ¼ exp 0:5kr?  r0 ?k22 =σ 2r and cs ðχ; χ 0 Þ ¼  2 exp 0:5kχ  χ 0 k2 =σ 2s : Here, k:k2 denotes the l2− norm. The coefficients, σr, σs are chosen to determine the amount of filtering. In our case, they are set to σr = 2.2μm and σs = 3% of the maximum value of the input data χweiner,bf (r⊥,z), respectively. Clearly, the output of the post-processing has better depth sectioning compared to the input image. Different structures and materials, which are not visible in the raw input, appear nicely in the output. 8

Data availability. The data that support the findings in this paper are available upon request.

Received: 27 April 2017 Accepted: 8 June 2017

References 1. Yamada, K. M. & Cukierman, E. Modeling tissue morphogenesis and cancer in 3D. Cell 130, 601–610 (2007). 2. Pampaloni, F., Reynaud, E. G. & Stelzer, E. H. The third dimension bridges the gap between cell culture and live tissue. Nat. Rev. Mol. Cell Biol. 8, 839–845 (2007). 3. Diaspro A. Optical Fluorescence Microscopy. Optical Fluorescence Microscopy by Alberto Diaspro Berlin (Springer, 2011). 4. Hell, S. W. Far-field optical nanoscopy. Science 316, 1153–1158 (2007). 5. Lippincott-Schwartz, J., Altan-Bonnet, N. & Patterson, G. H. Photobleaching and photoactivation: following protein dynamics in living cells. Nat. Cell Biol. 5, Suppl S7-S14 (2003). 6. Hoebe, R. A. et al. Controlled light-exposure microscopy reduces photobleaching and phototoxicity in fluorescence live-cell imaging. Nat. Biotechnol. 25, 249–253 (2007). 7. Tsien, R. Y. The green fluorescent protein. Annu. Rev. Biochem. 67, 509–544 (1998). 8. Crivat, G. & Taraska, J. W. Imaging proteins inside cells with fluorescent tags. Trends Biotechnol. 30, 8–16 (2012). 9. North, A. J. Seeing is believing? A beginners’ guide to practical pitfalls in image acquisition. J. Cell Biol. 172, 9–18 (2006). 10. Graf, B. W. & Boppart, S. A. Imaging and analysis of three-dimensional cell culture models. Methods Mol. Biol. 591, 211–227 (2010). 11. Pawley, J. & Masters, B. R. Handbook of biological confocal microscopy. Opt. Eng. 35, 2765–2766 (1996). 12. Masters B. R. & So P. T. C. Handbook of Biomedical Nonlinear Optical Microscopy (Oxford University Press, 2008). 13. Keller, P. J., Schmidt, A. D., Wittbrodt, J. & Stelzer, E. H. Reconstruction of zebrafish early embryonic development by scanned light sheet microscopy. Science 322, 1065–1069 (2008). 14. Planchon, T. A. et al. Rapid three-dimensional isotropic imaging of living cells using Bessel beam plane illumination. Nat. Methods 8, 417–423 (2011). 15. Mertz, J. Optical sectioning microscopy with planar or structured illumination. Nat. Methods 8, 811–819 (2011). 16. Chen, B. C. et al. Lattice light-sheet microscopy: imaging molecules to embryos at high spatiotemporal resolution. Science 346, 1257998 (2014). 17. Zernike, F. How I discovered phase contrast. Science 121, 345–349 (1955). 18. Way, D. Principles and applications of differential interference contrast light microscopy. Micros. Anal. Light Microsc. Suppl. 20, 9–11 (2006). 19. Popescu, G. Quantitative Phase Imaging of Cells and Tissues (McGraw-Hill, 2011). 20. Cogswell, C. J., Smith, N. I., Larkin, K. G. & Hariharan, P. Quantitative DIC microscopy using a geometric phase shifter. in Part of Photonics West BiOS’97 (International Society for Optics and Photonics, 1997). 21. Shribak, M., LaFountain, J., Biggs, D. & Inoue, S. Orientation-independent differential interference contrast microscopy and its combination with an orientation-independent polarization system. J. Biomed. Opt. 13, 014011 (2008). 22. Preza C., Snyder D. L., Rosenberger F. U., Markham J. & Conchello J.-A. in Optical Science, Engineering and Instrumentation'97 (International Society for Optics and Photonics, 1997). 23. Mehta, S. B. & Sheppard, C. J. R. Equivalent of the point spread function for partially coherent imaging. Optica 2, 736–739 (2015). 24. Shribak, M., Larkin, K. G. & Biggs, D. Mapping optical path length and image enhancement using quantitative orientation-independent differential interference contrast microscopy. J. Biomed. Opt. 22, 16006 (2017). 25. Popescu, G. et al. Fourier phase microscopy for investigation of biological structures and dynamics. Opt. Lett. 29, 2503–2505 (2004). 26. Popescu, G., Ikeda, T., Dasari, R. R. & Feld, M. S. Diffraction phase microscopy for quantifying cell structure and dynamics. Opt. Lett. 31, 775–777 (2006). 27. Bhaduri, B. et al. Diffraction phase microscopy: principles and applications in materials and life sciences. Adv. Opt. Photon. 6, 57–119 (2014). 28. Kandel, M. E., Teng, K. W., Selvin, P. R. & Popescu, G. Label-free imaging of single microtubule dynamics using spatial light interference microscopy. ACS Nano 11, 647–655 (2017). 29. Wang, Z. et al. Spatial light interference microscopy (SLIM). Opt. Express 19, 1016–1026 (2011). 30. Wang, Z. et al. Spatial light interference tomography (SLIT). Opt. Express 19, 19907–19918 (2011).

NATURE COMMUNICATIONS | 8: 210

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-00190-7

31. Bhaduri, B., Pham, H., Mir, M. & Popescu, G. Diffraction phase microscopy with white light. Opt. Lett. 37, 1094–1096 (2012). 32. Bhaduri, B., Tangella, K. & Popescu, G. Fourier phase microscopy with white light. Biomed. Opt. Express 4, 1434–1441 (2013). 33. Ding, H., Wang, Z., Nguyen, F., Boppart, S. A. & Popescu, G. Fourier transform light scattering of inhomogeneous and dynamic structures. Phys. Rev. Lett. 101, 238102 (2008). 34. Davis S. P., Abrams M. C. & Brault J. W. Fourier Transform Spectrometry (Academic Press, 2001). 35. Kim, T., Zhou, R. J., Goddard, L. L. & Popescu, G. Solving inverse scattering problems in biological samples by quantitative phase imaging. Laser Photonics Rev. 10, 13–39 (2016). 36. Choi, W. et al. Tomographic phase microscopy. Nat. Methods 4, 717–719 (2007). 37. Charriere, F. et al. Living specimen tomography by digital holographic microscopy: morphometry of testate amoeba. Opt. Express 14, 7005–7013 (2006). 38. Kim, T. et al. White-light diffraction tomography of unlabelled live cells. Nat Photonics 8, 256–263 (2014). 39. Ralston, T. S., Marks, D. L., Carney, P. S. & Boppart, S. A. Interferometric synthetic aperture microscopy. Nat. Phys. 3, 129–134 (2007). 40. Merola, F. et al. Tomographic flow cytometry by digital holography. Light-Sci. Appl. 6, e16241 (2017). 41. Isikman, S. O., Bishara, W., Zhu, H. & Ozcan, A. Optofluidic tomography on a chip. Appl. Phys. Lett. 98, 161109 (2011). 42. Frigault, M. M., Lacoste, J., Swift, J. L. & Brown, C. M. Live-cell microscopy— tips and tools. J. Cell Sci. 122, 753–767 (2009). 43. Tzur, A., Kafri, R., LeBleu, V. S., Lahav, G. & Kirschner, M. W. Cell growth and size homeostasis in proliferating animal cells. Science 325, 167–171 (2009). 44. Su, T. T. & O’Farrell, P. H. Size control: cell proliferation does not equal growth. Curr. Biol. 8, R687–R689 (1998). 45. Godin, M. et al. Using buoyant mass to measure the growth of single cells. Nat. Methods 7, 387–390 (2010). 46. Park, K. et al. Measurement of adherent cell mass and growth. Proc. Natl Acad. Sci. USA 107, 20691–20696 (2010). 47. Barer, R. Interference microscopy and mass determination. Nature 169, 366–367 (1952). 48. Su, T. T. & O’Farrell, P. H. Size control: cell proliferation does not equal growth. Curr. Biol. 8, R687–R689 (1998). 49. Wang, R. et al. Dispersion-relation phase spectroscopy of intracellular transport. Opt. Express 19, 20571–20579 (2011). 50. Kandel, M. E. et al. Three‐dimensional intracellular transport in neuron bodies and neurites investigated by label‐free dispersion‐relation phase spectroscopy. Cytometry Part A 91, 519–526 (2017). 51. Heng, Y. W. & Koh, C. G. Actin cytoskeleton dynamics and the cell division cycle. Int. J. Biochem. Cell. Biol. 42, 1622–1633 (2010). 52. Gnoth, C. et al. Final ART success rates: a 10 years survey. Hum. Reprod. 26, 2239–2246 (2011). 53. Graham, J. et al. Day 3 morphology is a poor predictor of blastocyst quality in extended culture. Fertil. Steril. 74, 495–497 (2000). 54. Guerif, F. et al. Limited value of morphological assessment at days 1 and 2 to predict blastocyst development potential: a prospective study based on 4042 embryos. Hum. Reprod. 22, 1973–1981 (2007). 55. Rijnders, P. M. & Jansen, C. A. The predictive value of day 3 embryo morphology regarding blastocyst formation, pregnancy and implantation rate after day 5 transfer following in-vitro fertilization or intracytoplasmic sperm injection. Hum. Reprod. 13, 2869–2873 (1998). 56. Marhuenda-Egea, F. C. et al. A crucial step in assisted reproduction technology: human embryo selection using metabolomic evaluation. Fertil. Steril. 94, 772–774 (2010). 57. Sakkas, D. & Gardner, D. K. Noninvasive methods to assess embryo quality. Curr. Opin. Obstet. Gynecol. 17, 283–288 (2005). 58. Scott, R. et al. Noninvasive metabolomic profiling of human embryo culture media using Raman spectroscopy predicts embryonic reproductive potential: a prospective blinded pilot study. Fertil. Steril. 90, 77–83 (2008). 59. Seli, E., Botros, L., Sakkas, D. & Burns, D. H. Noninvasive metabolomic profiling of embryo culture media using proton nuclear magnetic resonance correlates with reproductive potential of embryos in women undergoing in vitro fertilization. Fertil. Steril. 90, 2183–2189 (2008).

NATURE COMMUNICATIONS | 8: 210

60. Vergouw, C. G. et al. Metabolomic profiling by near-infrared spectroscopy as a tool to assess embryo viability: a novel, non-invasive method for embryo selection. Hum. Reprod. 23, 1499–1504 (2008). 61. Katz-Jaffe, M. G., McReynolds, S., Gardner, D. K. & Schoolcraft, W. B. The role of proteomics in defining the human embryonic secretome. Mol. Hum. Reprod. 15, 271–277 (2009). 62. Weinerman, R. et al. Morphokinetic evaluation of embryo development in a mouse model: functional and molecular correlates. Biol. Reprod. 94, 84 (2016). 63. Halvaei, I., Khalili, M. A. & Nottola, S. A. A novel method for transmission electron microscopy study of cytoplasmic fragments from preimplantation human embryos. Microsc. Res. Tech. 79, 459–462 (2016). 64. Zucker, R. M. Whole insect and mammalian embryo imaging with confocal microscopy: morphology and apoptosis. Cytometry A. 69, 1143–1152 (2006). 65. Squirrell, J. M., Wokosin, D. L., White, J. G. & Bavister, B. D. Long-term twophoton fluorescence imaging of mammalian embryos without compromising viability. Nat. Biotechnol. 17, 763–767 (1999). 66. Shribak, M. Quantitative orientation-independent DIC microscope with fast switching shear direction and bias modulation. J. Opt. Soc. Am. A. Opt. Image. Sci. Vis. 30, 769–782 (2013). 67. Gonzalez, R. C., Woods, R. E. & Eddins, S. L. Digital Image Processing Using MATLAB. (Pearson, 2004). 68. Tomasi, C. & Manduchi, R. in Sixth International Conference on Computer Vision, 1998 (IEEE, 1998).

Acknowledgements This work was supported by the National Science Foundation (CBET-0939511 STC, DBI 14-50962 EAGER, IIP-1353368). T.H.N. would like to acknowledge the Computational Science Engineering fellowship and Yuen T. Lo award for partially funding this research.

Author contributions G.P. and T.H.N. proposed the idea. G.P. and T.H.N. developed the theoretical description of the method. T.H.N. and M.E.K. developed the algorithm and code to solve the inverse problem. M.E.K. developed the code to control the imaging instrument. T.H.N. and M.E.K. performed imaging. M.R. and M.B.W. performed embryo preparation. All authors contributed to writing the manuscript.

Additional information Supplementary Information accompanies this paper at doi:10.1038/s41467-017-00190-7. Competing interests: G.P. has financial interest in Phi Optics, Inc., a company developing quantitative phase imaging technology for materials and life science applications. The remaining authors declare no competing financial interests. Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/ Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2017

| DOI: 10.1038/s41467-017-00190-7 | www.nature.com/naturecommunications

9