Highly parallel, interferometric diffusing wave spectroscopy for monitoring cerebral blood flow dynamics: supplementary material

This document provides supplementary information to “Highly parallel, interferometric diffusing wave spectroscopy for monitoring cerebral blood flow dynamics,” Here we relate the complex optical field autocorrelation to that of the in -phase field component measured by interferometric diffusing wave spectroscopy (iDWS), as well as summarize the multimode fiber solver, statistical simulation, and diffuse correlation spectroscopy fitting model used in the main text. Also , we include complete SANR and speckle number simulation results for variable horizontal pixel binning and vertical slit height. Furthermore, we provide experimental confirmation of optimal horizontal pixel binning and shot noise limited performance of the experimental setup of multimode iDWS. Finally, we investigate repeatability of in vivo cerebral blood flow monitoring, contact iDWS measurements, and effects of ambient light.


S1. AUTOCORRELATION OF IN-PHASE OPTICAL FIELD COMPONENT
From Eq. (7), the time-dependent heterodyne signals for an interference pattern measured by a senor array can be written as which can be expanded as a matrix equation by including delay time, td, as an additional dimension:  T  T  T  T  P  P  t  A  A  t  T  T  T  T  P  P  t   T  T  T  T  P  P  t   T  T  T  T where td is a delay time, δtd is the sampling interval, and Tp, m is given by Eq. (6). The heterodyne signal, AC, ( ), for a given sensor element only includes the in-phase (real) part of ∑ , • , * ( ). We assume that ∑ , • , * ( ) is a complex, zero-mean, and circularly symmetric Gaussian random variable, ( ) , with a known autocorrelation, 1 ( ) = 〈 * ( ) ( + )〉. The corresponding autocorrelation function, 1 ( ), based on the in-phase (real) part of ( ) can thus be expressed as where 〈 ( ) ( + )〉 = 〈 * ( ) * ( + )〉 = 0, due to circular symmetry and independence of the real and imaginary parts of ( ). Then, Eq. (S3) can be further written as Therefore, if 1 ( ) is real, as is the case for many dynamical systems of interest, 1 ( ), calculated from the in-phase (real) field component, provides the same information as 1 ( ). Note that all 1 ( ) mentioned in the following and primary manuscript indicate the in-phase optical field component-based autocorrelation function 1 ( ). To generate MMITMs, vectorial electric fields of all core modes guided in the step-index MMF, with the geometry and numerical aperture (NA) shown in Fig. S1, are determined by a general fiber solver (GFS) (FIMMWAVE, Photon Design) [2]. This GFS can model fibers of arbitrary (real) refractive index (RI), using a rigorous solution to the vectorial wave equation in cylindrical coordinates. The solution assumes isotropic media and a perfect electrical conductor as the outer boundary condition [3]. In Fig. S1, the diameters of the core and cladding layers are directly determined from the specifications of the MMF used in the multimode beamsplitter (FOBS-22P-1111-105/125-MMMM-850-95/5-35-3A3A3A3A-3-1-NA=0.15, OZ Optics), while the RI of the core is set to achieve an NA of 0.15, given an assumed pure silica

S3. MMITM COMPUTATION, STATISTICAL SIMULATION, AND DATA PROCESSING
Here we describe steps in the simulation to determine the SANR and speckle number in Sec. 3.A of the primary manuscript.
1. As mentioned above, the MMITM T can be created from the transverse electric field components ( , ) = ( , ) + ( , ) of vectorial core modes of the MMF (see Sec. S2). Each element of the complex MMITM T can be calculated by Eq. (6) for a given aSlit and NMode. To determine T, the complex amplitudes AR, n are assumed to have identical magnitudes and random phases, corresponding to NMode uncorrelated code modes in the reference arm. The total guided reference power, Mode � , � 2 , is always the same. In this work, four MMITMs with NMode of 1702, 1300, 900, and 500 are created, removing modes in order of increasing mode index (decreasing mode number) to decrease NMode. Because of the fixed input reference light polarization, any pair of orthogonally polarized modes with the same mode index have a fixed phase relationship, and therefore add to form a single mode with a superimposed polarization. Compared to the 1702 excited core modes in the sample arm, the effective maximum number of core modes is 851 (i.e. half of 1702) for the reference arm in our setup [ Fig. 2(a)]. So we set the reference modes with the same mode index have the identical phase in the MMITM computation. By comparison, in the sample arm, the complex amplitudes of any pair of orthogonally polarized modes with the same mode index are uncorrelated due to sample dynamics, leading to a random polarization. Finally, P (the output dimension) is always 512 for preliminary generation of MMITMs. Pixel binning (e.g. Sec. 3.A of the primary manuscript) is accomplished by summing either heterodyne signals or MMITM rows.

2.
The complex amplitudes, , * ( ), for sample core modes are created by the algorithm described in Ref. [5]. * ( ) constitutes 1702 time series of 1702 identically distributed, mutually uncorrelated, complex, circularly symmetric, and zero mean Gaussian random variables with an assumed autocorrelation decay rate of 0.05 μs -1 . This is relevant for in vivo experiments because the solution to the correlation diffusion equation for a semi-infinite homogenous medium [i.e. Eq. (S9)] can be approximated as an exponential decay for small lag times [6,7]. The sampling rate (i.e. 1/δtd) and length of , * ( ) are set as 333 kHz and 33333, respectively, corresponding to an integration time of ~0.1 second for each autocorrelation. The heterodyne signals without any noise [i.e. AC ( )] can then be generated as per Eq. (S1), based on the MMITM and * ( ).

3.
In real measurements, noise is present. For heterodyne detection, if the detected sample photon number is much less than the detected reference photon number, the main noise source is shot noise in the reference photon number, which obeys a Poisson distribution. In the following simulations, the total sample photon number, NS, detected by the 512 camera pixels is set as ~41 per sampling interval, δtd, of 3 μs, where the total reference photon number, NR, is set as 10 5 times higher than NS, estimated from experimental results. Also, camera noise, with uniform statistics across camera pixels, is included. So we include two additional Gaussian random variables of SN, ( ) and CN, ( ) as the shot noise and camera noise for each pixel, respectively, where where te is the exposure time of camera pixel for each sampling period, E is the single photon energy. Thus, the noise-added heterodyne signals, AC, ( ), can be expressed as

4.
Based on AC, ( ), horizontal pixel binning can be applied by summing heterodyne signals over every 512/NPixel pixels to generate NPixel time series that estimate NPixel field autocorrelations. Each field autocorrelation has a different scale due to the nonuniform distribution of the reference and sample light power across binned pixels. Summing the NPixel field autocorrelations achieves speckle averaging, yielding a single autocorrelation, 1 ( ), which can be fitted with the model, where ξ is the decay rate, and A and B are the sum of the squared heterodyne signals and the total white noise variance, respectively. The SANR is estimated by A/B. Based on A and B, the normalized field autocorrelation function, 1 ( ), can be also estimated.

S4. COMPLETE SANR AND SPECKLE NUMBER SIMULATIONS FOR VARIABLE HORIZONTAL PIXEL BINNING AND VERTICAL SLIT HEIGHT
The primary manuscript shows SANR and speckle number results for variable pixel binning only for a slit height aSlit of 1 (Fig. 3), and results for variable slit height only for a binned pixel number NPixel of 64 (Fig. 5). Here, we include the complete SANR and speckle number results in Fig. S2 and S3, respectively, estimated from the squared singular values corresponding to different aSlit, NPixel, and reference mode number NMode. The tradeoff between SANR and speckle number incurred by horizontal pixel binning, shown in Fig.  3, exists for all values of aSlit ( Fig. S2 and S3). Also, the aSlit threshold of 0.2, below which SANR and speckle number decrease in Fig. 5, remains for all other NPixel, except that speckle numbers for 1 ≤ NPixel < 8 in Fig. S3 seem insensitive to vertical slit height. In this regime of weakly correlated or uncorrelated binned pixels, speckle number depends on NPixel instead of sample mode number, leading to speckle numbers that depend less on aSlit. Finally, Fig. S2 and S3 indicate that reference mode number NMode has a minimal effect on SANR and, a minor, but slightly more noticeable effect on speckle number. More comprehensive investigation of reference mode number is needed in the future.

S5. SOLUTION TO CORRELATION DIFFUSION EQUATION
The modified diffuse correlation spectroscopy (DCS) model, for fitting experimental 1 ( , ), is where the fitting coefficients A and B account for the heterodyne signal and additive noise, respectively. Note that 1 DCS ( , ) is the normalized DCS autocorrelation model. Based on the CW Correlation Diffusion Equation, the normalized solution (i.e. normalized field autocorrelation), 1 DCS ( , ), for a semi-infinite homogenous turbid medium with an S-D separation of ρ is given by [8], = −1.44 −2 + 0.71 −1 + 0.668 + 0.064 , n is the ratio of RIs between the medium and air, and k0 is the wavenumber of the light propagating in the medium. For liquid phantoms, DB is the Brownian diffusion coefficient of moving scatters and α = 1; while for biological tissues, the term of αDB is referred to blood flow index (BFI), where the unitless factor α accounts for static scatters in the tissue. Fig. 7 in the primary manuscript only shows experimental phantom and in vivo results for horizontal pixel binning (NPixel of 64), determined to be optimal from simulation. To experimentally confirm the simulation results leading to this choice (Fig. 3), in Fig.  S4, we show experimental estimates of 1 ( ) for 10 different NPixel ranging from 512 to 1, measured from the liquid Intralipid phantom at 3.61 cm S-D separation. As in Fig. 3(a), noise in 1 ( ) estimates appears to be minimized for NPixel ~ 32−64. The minimum 95% confidence intervals for the fitted diffusion coefficient, DB, appear at NPixel of 64 in Fig. S4(k). This NPixel also yields the minimum RMSE in Fig. 3(b). Furthermore, Fig. S4(l) shows SANRs estimated from fitting experimental 1 ( ) (i.e. A/B) for different NPixel. SANRs calculated from the MMITMs with NMode of 1702 and aSlit of 1 [rescaled from Fig. 3(c) to account for the different sample photon number] are included for comparison (black solid curve). The simulated SANRs agree with the experimental results very well. Also, the theoretical SANR for binning fully correlated pixels is shown by the red dashed curve [ Fig. S4(l)]. In agreement with simulations [ Fig. 3(c)], the improvement in experimental SANR achieved by horizontal binning deviates from the theory for fully correlated pixels for NPixel ≲64. Experimental speckle numbers for different NPixel were estimated by substituting mean-subtracted time courses of binned pixels, measured at S-D separations of 0.76 and 1.07 cm, for heterodyne signals in Eq. (12). Since additive noise may cause overestimation of speckle number, relatively small S-D separations with higher SANRs were chosen to estimate speckle number. In Fig. S4(m), the two sets of experimentally obtained speckle numbers agree with simulations very well for NPixel ranging from 1 to 64. However, the SANR decay and consequent additive noise contribution for NPixel ≳ 64 causes overestimation of speckle number. The theoretical speckle number for binning fully uncorrelated pixels is shown by the red dashed curve in Fig. S4(m).

S6. EXPERIMENTAL CONFIRMATION OF OPTIMAL HORIZONTAL PIXEL BINNING
In agreement with simulations [ Fig. 3(d)], the enhancement of experimental speckle number with reduced binning deviates from the theory for uncorrelated pixels for NPixel ≳ 8. Besides confirming optimal horizontal pixel binning parameters experimentally, these results also confirm the predictive capabilities of our statistical iDWS model.
Finally, although the impact of vertical slit height is not directly investigated experimentally, as discussed in Sec. 3.A.2 of the primary manuscript, changes in SANR and speckle number as vertical slit height is varied arise from spatial correlations of measurements. Agreements of Fig. S4(l) and S4(m) with Fig. 3(c)  and 3(d), respectively, suggest that these spatial correlations are indeed well-represented by our simulations.

S7. SHOT NOISE LIMITED PERFORMANCE OF MULTIMODE IDWS SYSTEM
To test whether or not the experimental multimode iDWS setup achieves shot noise limited performance, SANR and additive noise variance were estimated at various reference photon numbers. These measurements were performed with the same Intralipid phantom used in Sec. 3.B.1, at an S-D separation of ~1.5 cm. With the experimental setup in Fig. 2(a), dynamic multimode interference time courses were recorded at nine different reference powers, without changing the sample power. Then, SANRs and additive noise variances were estimated by fitting the experimental 1 ( ) with Eq. (S8). A full well capacity of 13,000, specified by the manufacturer, was assumed in converting all quantities from camera gray level (DN) to photon number. SANR and additive noise variance (in units of photon number) for NPixel of 512 and 64 are shown in Fig. S5. With and without horizontal pixel binning, SANRs plateau with increased reference gain [ Fig. S5(a)]. Moreover, by comparing experimental additive noise variances to the shot noise limit (SNL) of the total reference photon number NR, we conclude that the multimode iDWS setup achieves shot noise limited performance for sufficiently large NR [ Fig. S5(b)]. Note that experimental noise variances slightly below the SNL, observed in Fig. S5(b) for NPixel of 512, may be due to an inaccurate assumed full well capacity.

S8. IDWS IS NOT SENSITIVE TO AMBIENT LIGHT
One advantage of the iDWS approach is coherent amplification of weak sample light levels, potentially rendering iDWS less sensitive to ambient light than photon counting methods. The effect of ambient light was tested on the Intralipid phantom used in Sec. 3.B.1 (S-D separation ~4.1 cm, temperature ~18.5°C). Two sets of dynamic interference patterns were continuously recorded for 10 seconds with ambient light levels of <3 lux and >500 lux (measured by a digital illuminance meter, LX1330B, Dr. Meter), respectively. Note that illumination ranges from 300 to 750 lux for medical examination rooms at hospitals. Under each illumination level, 10 field autocorrelations, 1 ( ) , were estimated in sequence with an integration time, tint, of 1 s, from the meansubtracted temporal fluctuations of NPixel = 64 binned pixels. Autocorrelation functions, 1 ( ), and diffusion coefficients, DB, determined by fitting with Eq. (S8), are shown in Fig. S6, both without (a) and with (b) ambient light. Neither the mean nor the standard deviation of DB across 10 trials is affected by ambient light (Fig. S6 insets), suggesting that iDWS is not sensitive to the ambient light.

S9. RELIABILITY AND REPEATABILITY OF CEREBRAL BLOOD FLOW IN VIVO
Repeated measurements of CBF were performed in two subjects to assess repeatability. For each subject, BFI time courses for an S-D separation of 2.5 cm were measured in multiple sessions lasting 30 seconds each, at the same forehead location. Exemplary data from one session for Subjects 1 and 2 are shown in Fig. S7(a) and (b), with BFI across multiple sessions [(0.74 ± 0.04) × 10 -8 and (1.05 ± 0.04) × 10 -8 cm 2 /s] in the upper right insets of each respective panel. Subject 2 consistently exhibited a higher pulsatile BFI than Subject 1, which is reflected in the larger withinsession standard deviation. Intrasubject coefficients of variation were 0.054 for Subject 1 and 0.038 for Subject 2, respectively, indicating repeatability across sessions. Though our results suggest that iDWS measurement errors are small in relation to differences between subjects, future studies designed to quantify reliability are warranted.  Contact measurements may be preferable to non-contact measurements for continuous monitoring in settings such as critical care. Fig. S8 shows a preliminary comparison of contact and non-contact measurements of BFI in the same subject at a similar forehead location, acquired less than 15 minutes apart. The variability between measurements is comparable to the withinsubject variability of the non-contact measurement (Fig. S7). While motion artifacts in non-contact and contact approaches should be rigorously compared, this comparison demonstrates feasibility of contact iDWS probes in the future.