Evaluation and Correction for Optical Scattering Variations in Laser Speckle Rheology of Biological Fluids

Biological fluids fulfill key functionalities such as hydrating, protecting, and nourishing cells and tissues in various organ systems. They are capable of these versatile tasks owing to their distinct structural and viscoelastic properties. Characterizing the viscoelastic properties of bio-fluids is of pivotal importance for monitoring the development of certain pathologies as well as engineering synthetic replacements. Laser Speckle Rheology (LSR) is a novel optical technology that enables mechanical evaluation of tissue. In LSR, a coherent laser beam illuminates the tissue and temporal speckle intensity fluctuations are analyzed to evaluate mechanical properties. The rate of temporal speckle fluctuations is, however, influenced by both optical and mechanical properties of tissue. Therefore, in this paper, we develop and validate an approach to estimate and compensate for the contributions of light scattering to speckle dynamics and demonstrate the capability of LSR for the accurate extraction of viscoelastic moduli in phantom samples and biological fluids of varying optical and mechanical properties.


Introduction
Biological fluids like synovial fluid, vitreous humor, cerebrospinal fluid, blood, lymph, and mucus are biopolymer solutions of water, protein macromolecules, and cells [1,2]. They serve as shock-absorbers, allergen and bacteria trappers, nutrient and oxygen distributers, and lubricants in different organ systems [2][3][4][5][6]. To fulfill these roles, bio-fluids maintain distinct viscoelastic behavior, exhibiting both solid and fluid-like features under different loading conditions and size scales [3][4][5][6][7]. Disease progression in multiple organ systems is frequently accompanied by altered viscoelastic properties of bio-fluids. For instance, in rheumatoid arthritis and osteoarthritis, the reduction of glycosaminoglycan hyaluronate and lubricin contents alters the viscoelastic properties of synovial fluid, compromising its shockabsorbing capability, in turn damaging cartilage under loading condition [8]. The abundance of evidence on the reduced viscosity of synovial fluid in the course of osteoarthritis has led to development of Viscosupplementation, a treatment approach in which diseased synovial fluid is replaced with an elastoviscous hyaluronan solution [9]. In the case of another bio-fluid, vitreous humor, altered viscoelastic properties are believed to be associated with age and numerous ocular pathologies such as retina tear and detachment [6,10]. The significant evidence on the role of biofluid viscoelasticity in disease initiation and progression, therefore, calls for the development of novel technologies for mechanical evaluation of bio-fluids in their native state to advance our understanding of bio-fluid pathologies, improve clinical disease diagnosis and facilitate the development of treatment strategies.
The viscoelastic modulus, G*(v) = G9(v)+i G0(v), defines the mechanical behavior of materials. It is traditionally measured using a mechanical rheometer by evaluating the ratio of a shear oscillatory stress at frequency v, exerted upon the sample to the induced strain. The elastic modulus, G9(v), is the real part of G*(v) and defines the energy stored in the sample. The viscous modulus, G0(v), is the imaginary part and represents viscous dissipation of the material [2,11]. The moduli G*, G9 and G0 are often expressed in the units of Pascal's (Pa). For instance, the typical values of |G*(v)| for soft tissues (at v = 1 Hz) are 60 m Pa for blood, 1 Pa for vitreous humor, 1.6 K Pa for fat, and 4.5 K Pa for muscle [6,12].
A major challenge is that temporal speckle intensity fluctuations, given by g 2 (t) curve, not only depend on the viscoelastic properties, but are also intimately influenced by optical properties of the medium, particularly by optical scattering [33,42,43]. Traditionally, in the limits of single scattering or strong multiple scattering (diffusive regime) media, dynamic light scattering (DLS) and diffusing wave spectroscopy (DWS) formalisms have been used to extract MSD from the measured g 2 (t) curve [17,19,[21][22][23]28,44,45].
Estimating G*(v) of bio-fluids that do not meet the limits of single scattering or diffusive regime, however, is not straightforward. Bio-fluids span a range of optical properties and may be weakly, moderately, or highly scattering. Moreover, pathologies often simultaneously alter both optical and mechanical properties of bio-fluids, further confounding the accurate estimation of viscoelastic properties using LSR [46]. In this paper, we show that for a majority of turbid media, that do not meet the limits of single scattering or diffusion approximations, DLS and DWS formalisms lead to erroneous measurements of MSD, and in turn result in inaccurate moduli estimates. We demonstrate that in order to accurately measure bio-fluid viscoelasticity, the influence of light scattering must be decoupled from that of mechanical properties in interpreting the speckle dynamics and g 2 (t). We therefore introduce a novel approach to improve the accuracy of LSR for the mechanical characterization of bio-fluids of varying optical properties by correcting for the influence of arbitrary optical scattering. In this approach, we measure sample optical properties from time-averaged speckle patterns, and implement a polarization sensitive correlation transfer Monte-Carlo ray tracing (PSCT-MCRT) algorithm to correct for the contribution of optical scattering in MSD evaluation. The close correspondence between LSR measurements and conventional rheology of phantom and bio-fluid samples, presented below, establishes the capability of the new approach in accurately evaluating the viscoelastic modulus, G*(v), for biological fluids of arbitrary optical properties.

Sample preparation
The studies below were performed using glycerol and bio-fluid samples. Various Glycerol-Water mixtures (G-W) were prepared at different proportions (60%G-40%W, 70%G-30%W, 80%G-20%W, 90%G-10%W, and 100%G) over a range of viscosities (0.01-1.4 Pa. s). Frozen bovine synovial fluid and vitreous humor (Animal Technologies, Tyler, TX) were warmed up to 37uC in a water bath prior to LSR testing. The glycerol and bio-fluid samples were chosen for their optical clarity, which allowed us to validate our approach via tuning their scattering properties by adding various concentrations of TiO 2 particles. In all cases, TiO 2 particles (dia. ,400 nm, Anatase, Acros organics, Belgium) were added to glycerol mixtures in multiple concentrations (0.04%-2% volume fractions, corresponding to reduced scattering coefficients, m9 s , : 1.3-84.8 mm 21 , N = 18) and thoroughly mixed in a vortex to ensure even dissemination of scattering particles. For both phantoms and bio-fluids, 1.5 ml of the samples were placed in a clear cuvette (Fischer brand, light path 10 mm, 1.5 ml) for LSR measurements, and 2 ml were used for mechanical testing. The LSR approach described here will be used in the future to evaluate the viscoelastic properties of bio-fluids in their native state, without adding extrinsic scattering particles. However, in this study, extrinsic TiO 2 particles were utilized purely for the purpose of validating our approach over a large range of optical scattering concentrations relevant to tissue.

LSR optical setup
Laser speckle frame series of samples were acquired using the optical setup shown in Fig. 1 [12][13][14]16]. Light from He-Ne Laser (633 nm) was coupled into a single mode fiber (SM600), and focused to a 50 mm spot on the sample. Cross-polarized laserspeckle patterns were acquired at 180u backscattering geometry via a beam-splitter using a high frame rate CMOS camera (PL-761, Pixelink, Ontario, Canada).

Measurement of speckle intensity temporal autocorrelation curves from time-varying speckle patterns
For glycerol samples, time-varying speckle images were captured for 2 second at 490 frames per second (fps) (Fig. 2). Speckle size was adjusted to at least twice the pixel size (12 mm) to maintain sufficient spatial sampling and contrast. For bio-fluids, a higher frame rate (840 fps) was used due to relatively faster speckle dynamics. The speckle intensity temporal autocorrelation curve, g 2 exp (t), was obtained by measuring the correlation coefficient of pixel intensities in the first speckle image (time point t 0 ) with subsequent images (time points t 0 +t, 0#t#2) in the image series, over a 2 s duration (Fig. 2, Block 1). Spatial averaging was performed over 40640 pixels, and several g 2 (t) curves evolving in time were averaged to enhance the accuracy of temporal statistics as follows [12][13][14][15][16]19,22,23,35,47]: where I(t 0 ) and I(t+t 0 ) referred to the speckle images at times t 0 and t+t 0 , , . pixels and , . t0 indicated spatial and temporal averaging over all the pixels in the images and for entire imaging duration (2 s), respectively.

MSD evaluation using DLS and DWS Formalisms
For single or strong multiple scattering media, DLS and DWS theories, respectively, have expressed the measured g 2 (t) (eqn. (1)) as a function of MSD as below [17,21,32,44,48,49]: Here k 0 is the wave number, n is the refractive index of the sample, and c is an experimental parameter that is generally assumed to be equal to 5/3 (,1.7) [17][18][19][20][21][22][23][24]27,28,31,44,45]. Since use of the DLS formalism is not valid for the moderate to strongly scattering samples used in this paper, the accuracy of our new LSR approach for measuring sample viscoelasticity is compared with that obtained using the DWS formalism (eqn. (3)) as described below.
Algorithm to derive corrected MSD and measure G*(v) using LSR To correct for the influence of arbitrary optical scattering in samples that do not meet the criteria of single scattering or light diffusion regimes, we have developed an algorithm that utilizes experimentally measured optical properties in a PSCT-MCRT model to establish a modified expression for g 2 (t) which corrects for influence of optical scattering in MSD calculations. The steps involved in measuring and correcting for optical scattering, the computational methods employed in the algorithm, and estimation of MSD and G*(v) in LSR are shown in the flowchart (Fig. 2) and detailed below.
i) Experimental evaluation of optical properties from time-averaged speckle images. Optical properties were derived experimentally by measuring the radial remittance or photon flux profile of samples from time-averaged speckle images using previously published methods (Fig. 2, Box 2) [14,[50][51][52][53][54][55]. Briefly, the speckle image series were temporally averaged over an ROI of 2966296 pixels at the CMOS sensor (Field of View (FOV) of 2 mm). The average pixel intensity values were converted to photon flux based on camera responsivity (28.1 DN/(n J/cm2)), gain (12.04 dB), exposure time (1 ms), and the (solid) angle of view. The radial photon flux profile was then fitted to the model predicted by the steady-state diffusion theory to evaluate the absorption and reduced scattering coefficients (m a , m9 s ) [14,51]. . Detailed flow chart of the compensation algorithm. Block 1: Speckle acquisition and g 2 (t) calculation: Speckle frame series are acquired with sufficient frame rate, ROI, and pixel to speckle size ratio. Speckle intensity temporal autocorrelation curves, g 2 (t), are evaluated for phantom and tissue samples using sufficient temporal and spatial averaging. Block 2: Measurement of optical properties: The radial remittance profile is evaluated from temporally averaged speckle intensities and is converted to the photon flux, y(r). Optical properties of the sample (m a and m9 s ) are derived experimentally by fitting the photon flux profile to the model obtained from steady-state diffusion theory. Block 3: PSCT-MCRT for simulating g 2 (t)-MSD expression: Experimentally evaluated optical properties, LSR configuration, and sample geometry are used in the PSCT-MCRT simulation to derive an expression for g 2 (t) as a function of MSD. Block 4: Evaluating MSD and |G*(v)|: Following the measurement of MSD using the modified expression, logarithmic slope of MSD, a(t) = h log ,Dr 2 (t)./h log t, is calculated and replaced in the simplified GSER to evaluate the viscoelastic modulus [18][19][20]22,23,25,26,[36][37][38][39][40][41]. Here K B is the Boltzman constant (1.38610 223 ), T is temperature (degrees kelvin), a is the scattering particle size, ,Dr 2 (1/v). corresponds to ,Dr 2 (t)., evaluated at t = 1/v, v = 1/t is the frequency, and C represents the gamma function. doi:10.1371/journal.pone.0065014.g002 Since the asymmetry parameter, g, is not trivially related to photon flux, the value of g and the corresponding scattering phase function were calculated from Mie theory predictions, which resulted in g = 0.6 for TiO 2 particles suspended in glycerol solutions (see Discussion) [56]. For predominantly scattering samples in this study, the optical absorption coefficients, m a , were negligible. The accuracy of this approach in estimating m9 s for glycerol and biofluid samples with varying TiO 2 concentrations was confirmed via comparison with Mie theory estimates [56].
ii) PSCT-MCRT simulation to establish the modified g 2 (t) and MSD relationship. The PSCT-MCRT algorithm below was employed to simulate g 2 MCRT (t) curves (Fig. 2, Box 3) and derive a modified relationship between MSD and g 2 (t), for samples with arbitrary optical properties. The PSCT-MCRT model incorporated all experimental LSR parameters, for a focused Gaussian beam (50 mm) illuminating the sample placed in a cuvette (10 mm light path, 1.5 ml) with m a and m9 s measured as above. A total of 10 5 photons were tracked from the source to the receiver (FOV of 2 mm). The temporal speckle fluctuations were modified by the polarization state of detected light. Therefore, PSCT-MCRT algorithm incorporated attributes of the polarization state by tracking the Stokes' vector, [I Q U V], with respect to the corresponding reference frame. Euler equations were used to modify the Stokes' vector upon scattering and transport within the medium [57,58]. At the receiver site, a final rotation was applied to redefine the Stokes' vector in the receiver coordinates system and since LSR setup captured the rapidly evolving speckle pattern of the cross-polarized channel ( Fig. 1), only the cross-polarized component of intensity was retained. To account for momentum transfer (Doppler shift) at each scattering event, the scattering wave vector, q = 2k 0 sin(h/2), was tracked, as well [17,35]. The total momentum transfer, defined as Y = gq 2 /(2k 0 2 ), with the summation over all scattering events involved in that path, represented the reduction of speckle intensity temporal autocorrelation due to all scattering events involved in each path [35,59]. Consequently, g 2 MCRT (t) was obtained by integrating the field temporal autocorrelation curves of all received rays, weighted by the corresponding momentum transfer distribution, P(Y), as [43,60]: From eqn. (4), it is noted that the term in brackets is simply the Laplace transform of P(Y), L{P(Y)} evaluated at 1/3k 0 2 n 2 ,Dr 2 (t). [43,61]. L{P(Y)} is equivalent to speckle field autocorrelation, g 1 MCRT (t), in turn related to speckle intensity autocorrelation curve, g 2 MCRT (t), through the Siegert relation as: g 2 (t) = 1+|g 1 (t)| 2 [27,28,35,44,62]. PSCT-MCRT only provided the statistical histogram of photons' P(Y) and generated a numerical solution for L{P(Y)} ( = g 1 MCRT (t)), and consequently g 2 MCRT (t). To simplify eqn. (4), a parametric function was fitted to L{P(Y)} as below: where S was the argument of the transform (complex frequency).
The parameters c and f were derived from PSCT-MCRT simulation by numerical calculation of the total momentum transfer distribution P(Y), based on the experimentally evaluated values for m a and m9 s and Mie theory calculations of g for each individual sample. Consequently, the following expression was derived for g 2 MCRT (t) as a function of MSD (see Discussion): The parametric functions of eqns. (5) and (6) were cross-checked for the limits of single scattering and diffusion approximations (eqns. (2) and (3)). For instance, for TiO 2 concentration of 0.04%, and m9 s = 1.3 mm 21 , the results of PSCT-MCRT gave rise to empirical parameters of c = 2/3 and f = 1, in which case eqn. (6) above converged to the DLS expression of eqn. (2). On the other hand, for TiO 2 concentration of 2%, corresponding to m9 s value of 84.8 mm 21 , MCRT results gave rise to c = 5/3 and f = 0.5, and the eqn. (6) converged to the DWS formalism (eqn. (3)). Since c and f are directly related to the sample optical properties, the need for repetitive execution of MCRT simulations can be eliminated by calculating the c and f parameters for a wide range of m a and m9 s values relevant to tissue, beforehand, and preserve the trends in a look up table (LUT). Thus, in the future by measuring the m a and m9 s of the sample, the corresponding c and f parameters can be simply obtained from the LUT without the need for executing PSCT-MCRT simulations.

Comparison of LSR with reference-standard mechanical testing
To validate our approach above, we compared LSR measurements of the magnitude of viscoelastic moduli, |G*(v)|, with those obtained using a reference standard mechanical rheometer (AR-G2, TA Instruments, MA). To conduct mechanical rheometry, the rotating rod and top plate (40 mm dia. stainless steel) exerted a shear oscillatory torque (stress) upon the sample over the frequency range of 0.1-100 Hz and measured the induced strain. For aqueous glycerol mixtures, tests were carried out at 25uC, and 2% strain was applied on the samples. Synovial fluid and vitreous humor samples were evaluated at 37uC with a 1% strain. The accuracy of our new modified approach in measuring the viscoelastic moduli of glycerol and bio-fluid samples was further compared with that obtained using conventional DWS (eqn. (3)).

Results
Influence of optical scattering on LSR measurements Fig. 3 shows the speckle intensity temporal autocorrelation curves, g 2 exp (t) for glycerol-water (G-W) mixtures of varying viscosities with 0.1% of TiO 2 scattering particles. As expected, for liquids of higher viscosity, g 2 exp (t) curves decayed slower due to reduced Brownian motion of TiO 2 particles compared to the lower viscosity samples. Here, given the identical optical properties of samples (0.1% TiO 2 ), a direct comparison of g 2 exp (t) curves enabled an accurate assessment of relative differences in sample viscosities. However, variation in scattering concentrations also modified speckle dynamics even in samples with identical viscosity, as demonstrated in Fig. 4, which displays measured g 2 exp (t) curves for 90% G -10% W (volume fraction) with viscosity ,0.25 Pa. S. While these samples had identical mechanical properties, optical scattering properties were tuned by mixing with different TiO 2 particle concentrations (0.04%-2%). To ensure that the addition of TiO 2 particles did not affect the sample mechanical properties, the viscoelastic moduli of glycerol suspensions with different TiO 2 particle concentrations (up to 2%) were measured. No detectable differences in G* were measured even at the highest TiO 2 concentration. Fig. 4 also displays the g 2 (t) curves obtained from eqns. (2) and (3) based on DLS and DWS formalisms (dashed and dot-dashed lines). From Fig. 4 it was evident that for most scattering concentrations, g 2 exp (t) curves fell somewhere between the dotted curves. In other words, by changing the TiO 2 concentration from 0.04% to 2%, g 2 exp (t) curves swept the gap between theoretical limits of single and multiple scattering, demonstrating a dramatic change in temporal speckle intensity fluctuations for samples of identical viscosities. Further, the results that follow established that the direct usage of DWS approximation for extracting the MSD of scattering particles caused erroneous estimation of |G*(v)|. The results of exploiting DLS formalism were not presented here, since most of our samples were visibly not dilute enough to be considered single scattering. Fig. 5(a) shows the emerging errors by displaying the measured MSD, from g 2 exp (t) curves of Fig. 4, using the DWS approximation for backscattering geometry (eqn. (3)) [31,35]. As expected, a large variation was observed between the curves and scattering dependence in the g 2 exp (t) plots were directly conveyed to MSD curves. Fig. 5(b) displays the resultant |G*(v)| curves, calculated by substituting the raw MSD curves of Fig. 5(a) in the GSER (eqn. (7)) and the |G*(v)| measured using the rheometer (dashed line) [18][19][20]22,23,25,26,[36][37][38][39][40][41]. Results of Fig. 5(b) were clearly biased by variations in scattering concentrations and failed to correspond with conventional rheology results (dashed line), for most curves. In particular, |G*(v)| was over estimated using DWS formalisms, especially at lower concentrations, due to slower decay of the g 2 (t) curve, influenced by lower optical scattering independent of viscoelastic properties. It was only at TiO 2 concentration of 2%, that strong multiple scattering dominated, the diffusion approximation was valid, and |G*(v)| approached the results of mechanical rheometry.
Results of LSR using the new optical scattering correction algorithm i) LSR results for glycerol suspensions. Fig. 6 demonstrates the validity of our methods for evaluating the optical properties of phantom glycerol samples from time-averaged speckle images. Fig. 6(a) shows the radial profile of photon flux measured as a function of distance from the illumination center for the glycerol suspensions of 90%G-10%W with TiO 2 scattering concentrations ranging from 0.04%-2%. In Fig. 6(a) the number of remitted photons per unit area (photon flux) intensified at higher concentration and the inset of curves increased while slope of the photon flux profile became steeper. m9 s and m a were derived by fitting the photon flux profile (Fig. 6(a)) to theoretical models of the steady-state diffusion theory [14,51]. Fig. 6(b) shows the experimentally evaluated m9 s values plotted against corresponding Mie theory predictions [56]. Good agreement was observed (R = 0.96, p,0.0001), demonstrating the validity of the experimental approach in assessing sample optical properties. The results were more accurate for low to moderately scattering samples but started to diverge at higher concentrations as discussed below. Fig. 7(a), plots the MSD of particle dynamics in glycerol suspensions of 90%G-10%W, measured by employing the PSCT-MCRT based optical scattering correction algorithm. As noted, the variability between MSD curves over the range of scattering concentrations (0.04%-2%), was significantly reduced compared to Fig. 5(a) (which employed the DWS formalism). The impact of corrections was more pronounced in the intermediate times and residual small deviations were still observed at very early or long times, corresponding to initial decay and final plateau of g 2 exp (t). These mismatches were most likely due to certain experimental factors, as discussed later. Fig. 7(b) showed the LSR evaluation of |G*(v)| for the 90%G-10%W samples measured by employing optical scattering compensation compared with the corresponding rheometer measurements (dashed line). Compared to Fig. 5(b), the optical scattering dependence of |G*(v)| curves was significantly reduced by employing the compensation algorithm. Moreover, the scattering compensated moduli corresponding to all scattering concentrations closely corresponded with the measurements of mechanical rheometer. Our results showed that while differences in optical properties dramatically modulated the g 2 (t) curves, a significant improvement was achieved in the LSR evaluation of viscoelastic moduli, by compensating for optical scattering variations versus the direct application of DWS formalism in the estimation of MSD, and calculation of the |G*(v)|.
ii) LSR results of biological fluids. Figs. 8 (a) and (b) show g 2 exp (t) curves measured from time-varying speckle images of synovial fluid and vitreous humor, respectively. Similar to the glycerol samples, the g 2 exp (t) decay accelerated with increased scattering in both cases. Since g 2 exp (t) decayed slower for synovial fluid compared to vitreous humor, it was expected that synovial fluid would have a relatively higher modulus. However, it was necessary to correct for the contribution of optical scattering prior to comparing absolute mechanical moduli. Fig. 9 shows the LSR results of |G*(v)| for synovial fluid (Fig. 9a) and vitreous humor (Fig. 9b) measured with and without optical scattering correction. The red diamonds represent average |G*(v)| values of synovial fluid and vitreous humor samples of Fig. 8, estimated using LSR based on the DWS expression (eqn. (3)) which did not take into account optical scattering variations. The purple squares correspond to the moduli resulted from corrected MSD values, using the modified expression of eqn. (6), derived from the compensation algorithm. The red and purple error bars stand for the standard error. Also depicted in this figure are the |G*(v)| values measured using a conventional rheometer (black solid line, round markers). It was evident that in the case of LSR with optical scattering correction, |G*(v)| exhibited a close correspondence with conventional mechanical testing. Moreover, |G*(v)| measured using DWS approximation resulted in an offset of about one decade relative to conventional rheometric testing results. This was due to slower decay of speckle intensity temporal autocorrelation curve, caused by relatively low concentration of TiO 2 particles as discussed later. From the results of Fig. 9, it was clear that synovial fluid had a slightly higher viscoelastic modulus, which was consistent with our initial observation of speckle fluctuations and with standard reference mechanical rheometry. Moreover, the non-Newtonian behavior of these bio-fluids, reflected in smaller slope of |G*(v)| and lower frequency dependence compared to viscous glycerol solutions, pointed to the complex viscoelastic

Discussion
In this study, we have developed a new approach to significantly improve the accuracy of LSR for evaluating the viscoelastic properties of bio-fluids by correcting for the influence of optical scattering on laser speckle intensity fluctuations. We anticipate that this work will potentially open new avenues for the application of LSR in clinical diagnosis, treatment monitoring, tissue engineering, and drug development.
In LSR it is critical to correctly deduce the MSD of scattering particles from speckle fluctuations, given by g 2 (t), to derive |G*(v)| (eqn. (7)) [18][19][20]22,23,25,26,[36][37][38][39][40][41]. The major difficulty, however, is that g 2 (t) curve not only depends on the sample viscoelasticity but also on its optical properties, which define light transport in the medium. Thus, LSR evaluation of viscoelastic properties is particularly challenging in bio-fluids which span a  It is observed that for samples of higher scattering particles' concentration, the net backscattered signal level increases. At the same time, the curves decay faster as a function of radial distance. Transport albedo (m9 s /(m9 s +m a )) and effective attenuation coefficient (!m a (m9 s +m a )) are derived by fitting the photon flux to theoretical models of the steady-state diffusion theory to further extract m9 s and m a [51]. In panel (b) Mie theory estimates of m9 s are shown, which are derived based on TiO 2 particle size, source wavelength, and the ratio of refractive indices of TiO 2 particles and glycerol solutions(refractive index mismatch). A close correspondence is observed between experimental and theoretical measurements of the m9 s (R = 0.96, P,0.0001), especially at lower scattering concentrations. For higher scattering concentrations, potential sedimentation of scattering particles, and particle interactions lead to distortion of photon flux curves and saturation of evaluated parameters. doi:10.1371/journal.pone.0065014.g006 range of optical properties. For instance, aqueous humor and vitreous humor are almost transparent, cerebrospinal fluid, synovial fluid, plasma, and lymph are moderately scattering, and mucus, sputum, milk, blood and bile are highly scattering and absorbing at different wavelengths. Furthermore, the optical (and mechanical) properties of bio-fluids are often altered in diseased states, which are of particular interest.
Previously, DLS and DWS approximations have been used to derive MSD from g 2 (t). In DLS, the well-defined scattering geometry and single scattering assumption simplify the g 2 (t) to an exponentially decaying function of MSD [32]. In the other extreme, using DWS, light diffusion theory approximates the path length distribution and provides an expression for g 2 (t) [32]. However, in moderately scattering samples and for small sourcedetector distances, DLS and DWS cannot be directly applied to the analysis of MSD (Figs. 5(a) & (b)). The g 2 (t) curves of glycerol samples, displayed in Fig. 4 demonstrate that for samples of identical mechanical properties, g 2 (t) can still be modulated by tuning optical scattering properties and in a moderately scattering material lower numbers of scattering counts can result in a g 2 (t) curve that decays slower compared to a rich scattering medium of similar viscosity (Fig. 4). We have shown that the direct use of DWS formalism leads to an underestimation of the MSD (Fig. 5a), and in turn results in inaccurate |G*(v)| values biased by the  Fig. 4 using eqn. (6). The modified expression of eqn. (6) resulted from PSCT-MCRT simulation of photon propagation and correlation transfer in LSR experimental setup considering the exact sample geometry and optical properties. Compared to Fig. 5(a), variability of MSD curves is significantly reduced, especially at intermediate times. Residual small deviations, still observed at very early or long times, are most likely due to electronic noise and speckle blurring, respectively. In panel (b) Generalized Stokes'-Einstein Relation is used to calculate |G*(v)| from corrected MSD. It is observed that the variability between measured |G*(v)| for different concentrations is considerably reduced, compared to Fig. 5(b). Moreover, a high correspondence is observed between LSR results for |G*(v)| and mechanical rheometry. doi:10.1371/journal.pone.0065014.g007 optical scattering concentration (Fig. 5(b)). Similarly, if DLS approximation is used, it may overestimate the MSD because for moderately scattering media, speckle intensity temporal autocorrelation curve decays faster compared to a single scattering scenario [48]. Therefore, a more precise g 2 (t) -MSD expression is required to account for variations in optical properties in order to estimate accurate |G*(v)| values using LSR.
Leveraging our earlier work, absorption and reduced scattering coefficients (m a , m9 s = m s 6(1-g)) were similarly evaluated in this study from the photon flux profile (Fig. 6(a)), obtained from temporally averaged speckle image series [14]. The reduced scattering coefficient, m9 s , measured using this method demonstrated close agreement with theoretical calculations particularly at low and intermediate scattering concentrations. Deviations at higher concentrations were likely caused by clumping (larger particles, larger g) and sedimentation of TiO 2 particles (given the density of 4.23 g/cm 3 for TiO 2 relative to ,1.26 g/cm 3 for glycerol mixture) which resulted in lower m9 s compared to Mie predictions. Moreover, at higher TiO 2 concentrations close proximity of adjacent particles could also lead to interactions of near-field radiation and reduce the backscattering efficiency, which influenced the measured m9 s values [32,64,65]. For predominantly scattering samples, used here, results were solely focused on the influence of m9 s variations on the speckle dynamics and the role of absorption was not studied. Nonetheless, optical absorption is expected to eliminate rays with longest optical paths, corresponding to a large number of scattering events, and decelerate g 2 (t) curves [32]. In the received back-scattered signal, attributes of scattering angular distribution were extensively washed off by multiple scattering. As a result, experimental evaluation of phase function and g was not trivial and instead theoretical Mie calculations were used to predict these parameters which resulted in g = 0.6 for TiO 2 particles suspended in glycerol suspensions. Thus, in the current study, the effect of scattering anisotropy was not addressed in experiments. In the future application of LSR in bio-fluids, in their native states, the typical values of scattering asymmetry parameter for tissues (g = 0.7-0.9) can be used [53].
Theoretical analyses by others indicate that for fixed m9 s and particle size values, the changes in g have minimal effect on the g 2 (t) trend [43,61]. Our PSCT-MCRT simulations of g 2 MCRT (t) could independently confirm these findings. In reality, however, the influence of the anisotropy parameter cannot be isolated from the MSD. This is because g is directly related with scattering particle size, a, and in turn with the MSD. For instance, larger particles have a larger g value, and also exhibit slower Brownian displacements, i.e. small MSD, and slow down the g 2 (t) trend.
As shown in Fig. 6(b), experimentally measured m9 s values were validated via comparison with corresponding values obtained from Mie theory predictions that assumed near identical size spheres [56]. A DLS-based particle sizer (Zeta Sizer, Malvern Instruments, USA) was used to determine particle size, which provided an average hydrodynamic diameter of 400 nm and a polydispersity index (PdI, i.e. normalized distribution width) of 0.3 for the TiO 2 particles used in our experiments. Given the low PDI of 0.3 (below 0.5), the particles could be considered sufficiently mono-dispersed for Mie theory approximation. The TiO 2 particles were nonspherical pyramidal crystals [66]. However, due to their small sizeparameter, (radius times wave number, kr,,5) and random orientations in the suspension, the use of Mie theory predictions could be further justified [67,68].
The PSCT-MCRT model used here incorporated the experimental configurations of the LSR setup (Fig. 1), the finite sample geometry, and measured optical properties to derive a modified expression that related g 2 (t) with MSD in samples that do not meet the criteria of single scattering or diffusive regime. To this end, at first speckle field temporal autocorrelation function, g 1 (t) was derived in terms of MSD, as shown by the term in brackets in eqn. (4). Next, the Siegert relation [27,28,35,44,62] was used to express g 2 (t) in terms of g 1 (t) and MSD, as shown in eqns. (4) and (6) and the flowchart of Fig. 2. The modified expression of eqn. (6) converged to the classical DLS and DWS expressions (eqns. (2 and 3)) in the limits of single and rich multiple scattering, respectively, as explained in the materials and methods section. The parameter, c, in the above equations is an scattering dependent variable, governed by the sample optical properties, concentration of light scattering particles, the polarization state of collected backscattered light, reflectivity of samples walls (i.e. boundary conditions), and other experimental parameters [69][70][71]. In the limit of strongly scattering medium, the c is traditionally assumed to be 5/3 [17,19,[21][22][23]28,44,45]. Consistent with this assumption, our MCRT results exhibited the asymptotic value of c = 5/3 in the limit of rich, multiple scattering samples. The MSD values ( Fig. 7(a)), estimated using the PSCT-MCRT-derived expression, demonstrated significant improvement compared with the DWSbased approach in compensating for scattering variations. The minimal residual MSD variability at early times was likely due to CMOS sensor noise, pronounced in the few initial points of g 2 (t) curve particularly in low scattering samples with lower speckle intensity. Deviations at long times of Fig. 7(a) were caused by rapid speckle dynamics and insufficient camera frame rate which led to spatio-temporal speckle blurring and an artificial plateau in g 2 (t) curve which was transferred to long time MSD slope. This minor hardware limitation could be avoided by exploiting higher frame rates and maintaining sufficient signal to noise ratio. The |G*(v)| curves of Fig. 7(b), estimated from the corrected MSD values, showed high correspondence with the results of mechanical rheometry. The variability between measured |G*(v)| for different concentrations was considerably reduced and the absolute value of the |G*(v)| curves converged to that obtained using reference-standard mechanical testing. The minimal scattering-dependent discrepancy in the |G*(v)|, at low and high frequencies, resulted from initial and long times variations of MSD as explained above.
Similarly, in the low viscosity bio-fluids, studied here, LSR measurements of |G*(v)| at low frequencies (,0.1 Hz) were influenced by speckle blurring caused by insufficient camera frame rate. These effects were minimized by subtracting the background signal (caused by speckle blurring) to improve speckle contrast prior to evaluation of g 2 (t) (Fig. 8). With these corrections, LSR could measure viscoelastic behavior over a frequency range that spanned three decades (0.1-100 Hz). High frame rate speckle acquisition could potentially extend this frequency range through probing speckle dynamics and sample mechanics at finer temporal resolution and consequently higher frequencies. Also, it would improve speckle contrast and extend the validity of LSR results to lower oscillation frequencies and enable probing structural and rheological behavior of bio-fluids in more details. For these biofluids, constraints of rotational rheometer for evaluating the |G*(v)| at high oscillation frequencies limited the comparison of LSR with mechanical rheometry to below 10 Hz (Figs. 9(a) and (b)). At higher frequencies, the torque (stress) applied by the rheometer motor shaft to rotate the rod and top plate rapidly increased with frequency (due to inertia) and exceeded the torque needed to strain the low viscosity bio-fluids. This is an inherent limitation of mechanical rheometer that cause unreliable readouts of |G*(v)| at higher frequencies.
In the current LSR setup, focused beam illumination and fullfield collection led to a broad distribution of optical path lengths for received rays, with each length probing a different time-scale (frequency) of the sample dynamics [32]. The multi-speckle collection enabled shorter acquisition and better statistical accuracy by exploiting both ensemble and temporal averaging (eqn.(1)) [49,72]. To permit depth-resolved mapping of bio-fluids' mechanical properties, spatio-temporal processing of speckle patterns can be employed as previously described [14]. Low coherence interferometery (e.g. M-mode OCT) techniques can also probe particle dynamics in specific depths within the medium with superior resolution and potentially evaluate the viscosity [29,30]. However, there exists a tradeoff between higher depthresolution capabilities of OCT versus higher measurement sensitivity of LSR to particles' displacements (MSD). This is because in LSR the detected light has scattered multiple times over a volume of interest. Therefore, even minute motions (fraction of a wavelength) of particles give rise to cumulative phase changes over an ensemble of light paths within the measured volume and induce detectable reduction in speckle intensity temporal autocorrelation. By measuring g 2 (t) over multiple speckles as in LSR, particle displacements as small as a few angstroms can be detected [32,37,73]. In OCT however, single scattered light is detected. Therefore, a substantially larger displacement of particles is required to induce a sufficient path length change and noticeable reduction of speckle intensity temporal autocorrelation particularly in highly viscous media with smaller particle MSD [32,37,72,73].
As described earlier, in the current work TiO 2 particles were used to sufficiently validate our approach by tuning the reduced scattering coefficient over a range of scattering concentrations. Adding the extrinsic scattering particles, at various arbitrary concentrations enabled complete evaluation of scattering variations effects, and validation of the proposed compensation algorithm. In future, we anticipate that LSR will be used to measure bio-fluid viscoelasticity in the native state without addition of extrinsic light scattering particles. This will require additional information about particle size parameter, a, to be evaluated experimentally. To this end, the future LSR configuration may be coupled with polarization dependent analysis of diffused back-scattered light or angle-resolved detection of low coherence radiation to enable particle sizing and permit the quantification of G*(v) of bio-fluid in their native states [74][75][76][77]. In addition, LSR can potentially be conducted via needle-based probes or endoscopes to enable future in vivo use. By limiting the illumination and collection volume at the interrogation site residual scattering from surrounding structures could be potentially restricted. We anticipate that the demonstrated capability of LSR for the non-contact and accurate evaluation of viscoelastic properties and the potential of this technology in probing rheological properties of bio-fluids will open multiple new avenues for clinical applications of LSR in the future.