A comprehensive workflow and its validation for simulating diffuse speckle statistics for optical blood flow measurements

Diffuse optical methods including speckle contrast optical spectroscopy and tomography (SCOS and SCOT), use speckle contrast (κ) to measure deep blood flow. In order to design practical systems, parameters such as signal-to-noise ratio (SNR) and the effects of limited sampling of statistical quantities, should be considered. To that end, we have developed a method for simulating speckle contrast signals including effects of detector noise. The method was validated experimentally, and the simulations were used to study the effects of physical and experimental parameters on the accuracy and precision of κ. These results revealed that systematic detector effects resulted in decreased accuracy and precision of κ in the regime of low detected signals. The method can provide guidelines for the design and usage of SCOS and/or SCOT instruments.


Introduction
An accurate and often continuous assessment of microvascular, regional blood flow has many 27 implications for diagnosis and treatment of diseases and for the study of healthy physiology. 28 Despite continued efforts to establish practical means for measuring microvascular, regional 29 blood flow in a non-invasive manner, this remains an important unmet need. One potential 30 approach uses near-infrared, coherent light and the arising speckles after its diffusion [1][2][3][4]. 31 Coherent laser light can be used to non-invasively measure local microvascular blood flow in 32 tissue by detecting the fluctuating speckle patterns after light interaction with the tissue [5-9]. 33 For the purposes of this manuscript, we will focus on deep-tissue, i.e. those that utilize light that 34 penetrates up to several centimeters, measurements using photon diffusion. This is possible since Depending on the experimental setup, the imaged field of view will differ. In this example, source and the detector fibers are placed a certain distance ( ) from each other and are coupled to the laser and detector. The imaged field-of-view (imaged over × pixels includes the fiber core which in later steps will be used to calculate 2 over a specified region of interest ( ). Sub-figure b illustrates Step 1 of the simulations. In this step, the rate at which the speckles decorrelate, , is determined from the correlation diffusion equation (CDE). Using this value of , consecutive frames of correlated speckles are simulated so that their electric-field autocorrelation decays with . The intensity of these simulations are in arbitrary units, and independent of exposure time, . Instead they represent speckles measured during a finite time-bin width, , on the 1 curve. In order to simulate several values of , the process illustrated in b can be repeated several times to simulate the dependent change in . In Step 2 (sub-figure c), the arbitrary units of the simulated frames is scaled to represent realistic values of photon current rate, Φ, in units of photons/second. In Step 3 (sub-figure d), an exposure time is introduced to the simulations by summing over frames. This process additionally converts the units of the simulations from photons/s to photons. Various values of can be simulated from the same set of simulated frames of Step 1. In this case, the simulation of two values of exposure time, and , is shown. Multiplying the summed frames in units of photons by the quantum efficiency (QE) of the camera converts the units of the simulations to electrons (e − ). In Step 4 (sub-figure e), the detector effects are simulated by altering the simulated intensity statistics according to the specifications of real detectors. In Step 5 (sub-figures f and g), speckles are sampled over an area, or over pixels of several repetitions of simulations to estimate a value of 2 . The yellow dots represent 2 simulated for the and therefore simulated in Step 1. The two values of simulated in Step 3 are also shown. In the final step (Step 6, sub-figure h), the discrepancies in the exact form of the speckle autocorrelation decay between the solution for the CDE for a semi-infinite medium and the developed model is corrected for.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint and is also a function of these parameters. 144 2.2. Speckle statistics detected by a two dimensional detector array 145 We have simulated 2 for tissue with specific optical properties and blood flow by simulating 146 consecutive frames of correlated speckles which simulate their electric field autocorrelation with 147 a decorrelation time, , defined by the solution of the CDE for a semi-infinite medium [10]. The 148 methodology presented is independent of this solution and other solutions (layered, heterogeneous, 149 numerical) of the CDE could be utilized. For clarity, electric-field autocorrelation curves following 150 the solution of the CDE will be referred to asˆ1, while the simulated electric-field autocorrelation 151 curves are referred to as 1 . While the two are similar, there are slight differences which are 152 discussed below. Furthermore, the theoretical value of 2 derived from the CDE will be referred 153 to asˆ2 while the simulated values will be referred to as 2 .

154
In the first step of the simulation pipeline (Figure 1b), is derived fromˆ1. The derived 155 value of was used to simulate frames of individual speckles by modifying the copula method 156 developed in Ref. [42]. This method simulates consecutive two dimensional matrices of numbers individual matrix can be considered as a camera frame acquired in a speckle contrast experiment.

160
These matrices are referred to as "frames" ( ) simulating pixel coordinates , while imaging 161 speckles with diameter, Ø. Ø behaves as a scaling factor to put physical units for the pixel 162 size since the speckle diameter is approximately equal to the wavelength of light being used. 163 Therefore, choosing Ø to be equal to three pixels for a system modeling = 785 nm will scale 164 the width of each pixel to be equal to approximately 262 nm.

165
The autocorrelation, 1 , of the first frame, = 1 to the k th frame, = is given by where is the frame number and is a parameter related to the decorrelation of the frames.

167
In our adaptation we have defined to be a function of . Since has been defined as . (2) Each of the individual simulations of 1 consisting of = frames of speckles patterns 170 constitute an experiment, defined by . This process together with notation is illustrated in Figure   171 2. The simulations are simulated in arbitrary copula units. In addition, the frames are only 175 dependent on and every simulated frame represents a point on the 1 curve with a finite time-bin 176 width, . Since each frame has a defined and is simulated over an array × , the complete 177 notation is, ∼ ( ) . In this notation, the pre-superscript indicates the units of the simulated 178 frame. In this case, refers to the arbitrary "copula" units. The pre-subscript, ∼, indicates that 179 no effect of detector noise has been included in the simulated frame. The indices , and refer 180 to the pixel and frame. In order to convert ∼ ( ) to physical units, the arbitrary copula units must be scaled to a Figure 2. Illustration of how frames with a defined are simulated. First individual speckles are simulated on a grid of × pixels. These individual frames, , are correlated to each other and their electric-field autocorrelation, 1 , decay according to defined from semi-infinite theory ( Figure 1). One full simulation of a theoretical 1 curve ( 1 ) consisting of frames corresponds to one experiment, . This process is repeated several times resulting in several simulations of 1 . or experimentally. According to the photon diffusion theory, in a semi-infinite geometry, the 185 measured photon current rate, Φ( ), in units of photons/second, decreases with as:

202
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint This is done by adding = / consecutive frames: Note that with the introduction of exposure time, the simulated frames drop their indexing of 204 .

205
Finally, the simulated frames are converted from photons to electrons: Where is the quantum efficiency of the camera.
207 Table 1 summarizes the introduced notation to refer to the simulated frames.  The final step before using the simulations to calculate 2 is to simulate the effects of the main 210 types of detector noise on the simulated frames previously described, namely: photon shot noise,

211
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint dark signal non-uniformity (DSNU), dark current shot noise, and read-out noise [44,45]. This 212 step is illustrated in Fig. 1e. To simulate detector noise, the distribution of each of the types of 213 noise is considered, and random numbers are generated following the distribution. The notation 214 used to describe the generation of random numbers and their distributions is shown in Eq. 7 is the random number generated representing a certain intensity (in e − ) at pixel , .

217
Photon shot noise is a Poisson distributed noise source [44, 46]. Using the notation in Eq. 7, 218 the contribution of photon shot noise at each pixel , is described as: Where we have applied the definition of a Poisson distribution, ( ) = 2 ( ). In this case 220 ( ) = ∼ ( , ) (i.e. the measured intensity in e − (Eq. 6)). We have also included a new Where ( ) and 2 ( ) are the mean and variance of the DSNU specific to each detector.

228
Their values can typically be found in camera specification sheets. The variance of a logistic 229 distribution is given by 2 ( ) = ( 2 2 )/3 where is the shape parameter of the logistic 230 distribution.

231
The dark shot noise, similar to the photon shot noise (Eq. 8) is simulated by applying Poisson 232 distributed random numbers [44] to each pixel simulated in Eq 9: Finally, read out noise is simulated by assuming that it is a normally distributed noise 234 source [47]. Read out noise in CMOS cameras is added at each pixel and is independent of the 235 dark noise and the detected signal. Therefore, the contribution of the read out signal at each pixel,

236
, is simulated: where the mean and variance of the read-out signal ( ( ) and 2 ( )) are specific to each 238 detector and can be found in specification sheets or estimated from online camera simulators.

239
The total dark frame, , is then given by Putting everything together, the frames with shot noise, DSNU, dark shot noise, and read-out • ′ ′ ′ : shot noise and dark frame added, dark frame and shot noise corrected.

251
The definitions and notation for simulating detector noise is summarized in Table 2: Table 2. Table of definitions of the noise sources that are included in the simulations along with their corresponding distributions. The notation ( ; , 2 ) is used to define random numbers, , originating from a distribution, , with a mean value of, , and variance, 2 . † denotes parameters that can be found in camera specification sheets.

253
The final steps of the simulation pipeline require the calculation of 2 using the frames that 254 have been simulated. In the first step, 2 is directly calculated using the simulated frames. The will be affected by the model used for 1 , the noise is well described using the simplified single 263 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint exponential model. The definitions and notation related to 2 are summarized in Table 3. The 264 following sections will describe their calculations.
corrected for semi-infinite theory (Eq. 20) Table 3. Table of definitions for 2 . Three different variations of 2 are calculated: first 2 calculated directly from the integration of the double exponential 1 from CDE. This isˆ2. Secondly, 2 calculated directly from the simulated frames whose 1 ( 1 ) follows a single exponential form. This is 2 and outlined in Section 2.7. Thirdly, the model differences due to the differences in 1 is corrected. This is 2 ′ and is outlined in Section 2.8. Moreover, 2 and 2 ′ can be calculated either spatially or temporally.

Model uncorrected speckle contrast 266
So far the process for simulating the detection of speckle statistics on a 2D detector array and the 267 detector properties (Fig. 1 b to e) has been described. These steps can be repeated in order to 268 simulate several experiments ( , Fig. 2) for several different values of and therefore , for 269 calculating 2 in the temporal domain over , or for determining ( 2 ).

270
The next step in the pipeline is to use these frames to calculate values of 2 ( Fig. 1 f and g). 271 As mentioned previously, 2 can be measured spatially or temporally i.e. speckle statistics can 272 be determined spatially by using an area, , of pixels or temporally over the pixels in a set of 273 experiments, .
274 Spatial 2 is given by: . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint After the dark frame offset is corrected, the additional variance due to shot ( 2 ℎ ) and the 287 dark frame (dark and read out noise, 2 ) is corrected by subtracting their respective variances 288 from the signal variance, 2 = 2 ( ′ ′ ( , ) ) .

292
Variations in the noise correction can also be simulated. For example, the shot noise only 293 added frames, 2 , can be corrected in the following way: Where in this case, 2 = 2 ( ′ ( , ) ) and 2 ℎ = ( ′ ( , ) ) . that the values of 2 derived from the two will differ. In particular,ˆ1 describes a measurement 300 in a semi-infinite medium and a multi-scattering (diffuse) regime. Sinceˆ1 is a more realistic 301 solution to the CDE, rather than working with 2 derived from 1 , we introduce another variable, semi-infinite geometry.

307
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint Finally 2 ′ values are generated by generating normally distributed random numbers, , with 309 mean equal toˆ2 + and variance equal to 2 ( 2 ) : 2.9. Using the simulations to evaluate system performance The speckle contrast noise model was further used to design a speckle contrast system and define 339 the required detected electron count rate (e − /pixel/second) in order to accurately measure blood 340 flow in the adult human brain. An sCMOS camera by Basler (daA1920-160um, Basler AG, 341 Ahrensburg, Germany) was considered and simulated due its lightweight (15 g), compact size 342 (19.9 mm x 29.3 mm x 29 mm) and cheap price (<300e). Measurements were chosen to be 343 taken at of 2.5 cm and of 5 ms.

344
The required detected electron count rate to accurately measure 2 was determined by (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint using an 800 m core multi-mode fiber (0.22 NA). The imaged speckles had a size of Ø = 5 348 pixels. The value of of the system was previously determined to be approximately 0.2. Speckle 349 contrast data was acquired over 600 frames, and data was analyzed using an ROI of approximately 350 1100 pixels.

351
As in the setup (A) to validate the simulations, of the simulations was obtained from 2 352 recorded using a standard DCS device. In order to approximate the required detected electron  Table 4. The simulations used obtained from the 1 curve recorded using DCS 362 (Figure 3 a). was simulated to be 0.2 and Ø was set to 4 pixels to agree with the values of and 363 Ø of the experimental data. Both experimental and simulation results were obtained for exposure 364 times ranging between 0.1 ms and 5 ms in order to cover a range of detected electron intensities.

365
It was ensured that the average value of the simulated detected electron intensity matched the 366 experimental data (Figure 3 b). The resulting experimental and simulated standard deviation of 367 2 is shown in Figure 3 c. The calculated signal to noise ratio of 2 in Figure 3 Table 5 were simulated. These values were chosen as they are roughly 376 the expected values when measuring in human tissue. 1 was simulated for ranging from 0.5 to 377 4.5 cm for = 5 ms. Ø was chosen to equal three pixels in order to meet the requirements of 378 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint    (Figure 5 a), at =0.1 ms, the majority of the detected electron signal 420 after =2 cm originate from the detector rather than from speckles. Therefore, without applying 421 corrections, any value of 2 in this regime is not a reflection of speckle contrast, rather reflects a 422 "detector signal" contrast.

423
The bias term, (Eq. 19), is shown in Fig. 5 c and f Fig. 6 a and d.

428
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint is shown in Fig. 6 b and e, reflected as the percent error. The percent error 432 increases (accuracy decreases) with increasing reaching 5% at approximately 1.8 cm for short 433 (Fig. 6 b) and 2.5 cm for long ( (Fig. 6 e). Similarly, the precision of 2 ′ , represented as  in order to be able to sample at fast enough acquisition rates while also maximizing the number 448 of detected photons (Figure 5 d). 449 In speckle contrast optical tomography (SCOT) or speckle contrast diffuse correlation tomog-450 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint corresponding to the sampling of 1100 independent speckles. In Figure 7, was changed to 459 simulate the effects of the number of independently sampled speckles on the CV of 2 ′ .

460
As expected in Fig. 7 2048×2048 pixels by choosing a larger region of pixels.

467
As observed in Fig. 6 b and e, accuracy was seen to be higher at shorter and longer ,  We have further demonstrated in detail (without experimental comparison) the entire simulation 482 pipeline. Finally, in the following section we will demonstrate how these simulations can be used 483 to design and optimize a speckle contrast system.

484
Speckles were simulated using the parameters specified in Table 6. These parameters were 485 derived from the experimental results ( and Ø), properties of the camera defined by the 486 manufacturer, as well as data analysis ( ). The resulting experimental and simulated percent 487 error in 2 for varying detected electron count rates is shown in Fig. 9.

488
The experimental and simulated results are in good agreement with each other and suggest 489 that for the chosen detector, a minimum detected count rate on the order between 4 to 5×10 4 490 e − /pixel/second allows us to calculate 2 with approximately 5% error.

491
Using the derived acceptable minimum detected count rate as a guide in determining the 492 accuracy of raw data signal, the same device was placed on a human subject's forehead using 493 a of 2.53 cm and of 5 ms. Data was acquired at a frame rate of 100 fps. A summary of  . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  Table 6. Parameters that were used to simulate synthetic speckles based on experimental data taken using a Basler (daA1920-160um) CMOS camera on a liquid phantom. contrast can account for parameters such as the speckle to pixel size and . 515 We have shown that the simulations accurately represent experimentally observed behavior

583
In the present work we have introduced a method for simulating the formation and detection 584 of dynamic speckle patterns. The main application that we have focused on was the design 585 and characterization of a speckle a contrast system capable of measuring human adult cerebral 586 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint blood flow non-invasively. To this end, the simulation method was validated on a dynamic 587 liquid phantom, the details of speckle contrast signal as a function of and were studied, and 588 finally a system designed for human cerebral blood flow was characterized and validated on an 589 adult human subject. The simulation method has been shown to be useful when identifying 590 the lower bounds of detected electron count-rate to achieve the desired accuracy and precision 591 of speckle contrast signal. As speckle contrast signal is sensitive to detector noise effects at 592 low detected electron count-rates, characterizing these limits is advisable when developing any 593 speckle contrast system. 25, 097004 (2020).

641
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 4, 2023. ; https://doi.org/10.1101/2023.08.03.551830 doi: bioRxiv preprint