Flexible Computationally Efficient Platform for Simulating Scan Formation in Optical Coherence Tomography with Accounting for Arbitrary Motions of Scatterers

. A computationally efficient and fairly realistic model of OCT-scan formation in spectral-domain optical coherence tomography is described. The model is based on the approximation of discrete scatterers and ballistic character of scattering, these approximations being widely used in literature. An important feature of the model is its ability to easily account for arbitrary scatterer motions and computationally efficiently generate large sequences of OCT scans for gradually varying configurations of scatterers. This makes the proposed simulation platform very convenient for studies related to the development of angiographic processing of OCT scans for visualization of microcirculation of blood, as well as for studies of decorrelation of speckle patterns in OCT scans due to random (Brownian type) motions of scatterers. Examples demonstrating utilization of the proposed model for generation OCT scans imitating perfused vessels in biological tissues, as well as evolution of speckles in OCT scans due to random translational and rotational motions of localized (but not-point-like) scatterers are given. To the best of our knowledge, such numerical simulations of large series of OCT scans in the presence of various types of motion of scatterers have not been demonstrated before. © 2021 Journal of Biomedical Photonics & Engineering.


Introduction
Optical coherence tomography (OCT) has proven to be a powerful tool for characterization biological tissue and diagnostics of various pathologies. In addition to utilization of conventional structural OCT images and solving problems related to enhancement of speed OCT visualization, improvement of resolution, etc., much attention has been paid in recent years to the development of various extensions/modalities of OCT to enable additional types of contrast, in particular, realization of polarization-sensitive OCT imaging [1,2]. In recent years, some other OCT-modalities are emerging in biomedical research and clinical practice and have already proved their high utility for a broad range of various biomedical applications. In particular, one can mention OCT-based angiography [3][4][5][6] (already translated to some clinical applications, e.g., [7,8]), mapping of viscous properties of biological liquids [9,10], visualization of deformations (for elastographic and other applications) [11][12][13][14] including combined application of a few modalities [15,16]. The principles of OCT-based elastography and angiography are essentially based on the analysis of motions of scatterers in acquired sequences of OCT scans.
Realization of such scatterer-motion-based modalities requires application of advanced signalprocessing methods to extract signals related to J of Biomedical Photonics & Eng 7(1) 31 Mar 2021 © J-BPE 010304-2 informative motions of scatterers and, in contrast, to filter out various kinds of distortions related to masking motions. For example, in OCT-angiography it is necessary to separate bulk motions of the "solid" tissue (due to breathing, heart-beat pulses, etc.) from the detected motion of scattering particles in the blood flow, which is not a trivial issue in living tissues [17]. For real biological tissues, performing experimental studies in highly controllable and reproducible conditions often is rather difficult or even impossible. By this reason, for validation and improvement of such OCT-based methods/algorithms for discrimination and quantitative assessment of motions of scatterers, researches often use phantom experiments with better controllable conditions [18,19]. However, even in phantom experiments it is challenging to ensure the required properties of scatterers and parameters of their motions. Furthermore, for real physical experiments (even with phantoms), it is often impossible (or unacceptably complicated and expensive) to flexibly vary in sufficiently broad ranges various parameters/regimes of the used tissue samples, phantoms and OCT setup itself. In view of this, creating models of OCT-scan formation enabling realistic and flexible (in terms of properties of scatterers and OCT-system parameters) attracts high interest of researchers.
Among various approaches to such modeling it is reasonable to point out three main directions. An important group of such methods comprises various realizations of the Monte-Carlo approach, in which pseudo-random generation of trajectories for huge ensembles of photons is made to describe their propagation and scattering (such that multiple scattering can also be accounted for) [20][21][22][23][24][25]. Another group utilizes the concept of point-spread function [26,27] and the third one comprises methods in which optical wave propagation, scattering and backward propagation and reception are consistently described. Such models can reasonably be called full-wave ones, although their realizations differ rather significantly and may comprise elements of analytical and numerical description in various combinations [28][29][30][31][32][33][34][35][36]. Their general feature is that in these methods the point-spread function and the entire speckle pattern are obtained by considering the illuminating-beam propagation first towards scatterers and then the wave scattering and subsequent propagation of the back-scattered field. Its collection over the receiving aperture is also considered. For example, in Ref. [28] wave-field propagation and collection is performed numerically or to a significant degree analytically as in works [32][33][34]. Quite a useful approximation in many models is the concept of discrete scatterers. However, the latter should not necessarily be sub-wavelength in size, so that even larger (but subresolution scatterers with sizes significantly smaller in comparison with the coherence length in the axial direction and the beam cross section laterally) can already be reasonably modeled as localized discrete elements.
For simulations of OCT scans related to the development of such extensions as the above-mentioned angiographic visualization, it is necessary to generate rather large volumes of OCT data scans imitating subsequent acquisition of large series of OCT scans. Therefore, unlike applications for which a single scan or a few scans may be sufficient, for generation of large blocks of OCT data, the ability of the model to perform simulations sufficiently rapidly becomes of key importance. Another requirement to the model is its ability to naturally account for arbitrary motion of scatterers. Such requirements are difficult to satisfy in models based on the Monte-Carlo approach or in computationally demanding essentially numerical variants of full-wave models. In what follows, we present a fast, simple and yet reasonably realistic model of OCTscan formation based on the principles described in Refs. [32][33][34], in the framework of which arbitrary motions of scatterers can easily be introduced. Some examples demonstrating the usefulness of this model for simulations related to angiographic applications OCT will be presented.

Main features of the model
The approach used here is based on previous works [32][33][34], where the often-used approximation of ballistic single scattering was applied. This assumption is the basic principle of OCT-scan formation, for which effects of multiple scattering the illuminating beam and imperfectly ballistic forward propagation (because of fluctuations of the refractive index) play the role of distorting/degrading factors. Certainly, for some applications, such degrading factors may be also interesting to study (and they may be additionally introduced in modeling, like, e.g., in Ref. [36]). However, the very fact of successful realization of the OCT principle of imaging based on the ballistic approximation convincingly confirms that the ballistic scattering usually gives a strongly dominant contribution to the OCT-scan formation. Furthermore, even remaining within the ballistic approximation many important degrading factors can be efficiently studied, for example, straininduced and displacement-induced decorrelation noises that are critically important for realization of OCT-based elastography [37][38][39]).
The next important note for the present study, is that in the majority of OCT setups the illuminating beams are intentionally weakly focused (with the focal plane located near the middle depth of scan). This allows one to obtain fairly invariable beam diameter enabling nearly uniform lateral resolution within the entire imaged depth (typically ~1-2 mm). Consequently, the illuminating beam in such setups remains close to cylindrical one with uniform diameter. Although the scattered field from localized scatterers, such as nuclei of cells in biological tissues, does not have pronounced directivity and may even be close spherically diverging, the same lens of the OCT setup that forms the highly directed illuminating beam also collects the scattered signal from the medium. As a result, even for isotropically scattering particles, the reception is highly directive, such that ballistically backscattered components are mostly collected by the received aperture. Consequently, the phase of the received signal is almost entirely determined by the forthand-back propagation in the axial direction. For example, fairly small variations in the axial-distance by a quarter of wavelength causes the variation in the received-signal phase by p radians, such that constructive interference may become destructive and vice versa. In contrast, transversal coordinates of scatterers within the illuminating-beam cross section may experience even significantly greater variations with an insignificant effect on the received-signal phase (by this reason the transversal Doppler effect is known to be much smaller in comparison with the effect of axial motions with the same magnitude).
Bearing in mind the above-mentioned arguments, the proposed in Ref. [32] simple model of OCT-scan formation was based on the approximation of a cylindrical illuminating beam for which the transversal profile was supposed uniform for both phase and amplitude (see Fig. 1a). An additional simplification in Ref. [32] was that the amplitudes of the received signals were supposed to be independent of the depth. Despite the fact that the received amplitudes should actually demonstrate depth-dependence, for description of interferometric speckle formation, the absence of this dependence in modeling still is acceptable by the following reasons. Namely, for low-coherence OCT, the efficient mutual interference is possible only for small groups of scatterers located within sample volumes with the axial size limited by the rather small coherence length (typically ~5-10 μm). Consequently, in contrast to phases of scattered signals that may strongly vary for particles located within the sample volume, for the amplitudes of the received backscattered signals, the effect of the scatterer-depth variation within a particular sample volume is negligibly small. Therefore, for any particular depth, the speckle formation is fairly correctly described even neglecting the depth-dependence of amplitudes. However, if necessary, depth-dependent variations in the amplitude in the simulated OCT-scans on scales larger than the speckle size can readily simulated by introducing additional depth-dependent factors to imitate geometrical focusing/divergence of optical waves, their exponential decay due to both absorption and scattering (similarly to introducing these averaged factors in Monte-Carlo approaches). Also, if necessary, roll-off effects in spectral-domain OCT can also be accounted within the same basic assumptions (e.g., Ref. [35]).
Summation of individual spectral harmonics of the illuminating beam in such a model naturally describes the formation of the point-spread function in the axial direction. Such features of model [32] made it possible to efficiently use it for studying in detail the role of variation in the optical-spectrum width on the feasibility of correlational speckle tracking in OCT images of strained materials [37]. The intrinsic exceptional ease in accounting for displacements of scatterers in simple model [32] has also been very efficiently applied for studying important features of strain mapping in phase-resolved OCT [38][39][40][41]. The results of those simulations then have been successfully transferred for mapping strains and elasticity in numerous real situations, in which scatterers experienced predominantly axial interframe displacements [12,[42][43][44]. Sufficient correctness of elastographic simulations based on Ref. [32] has been later verified by comparison with more advanced realizations of this approach in works [33,34], in which rigorous consideration of forward propagation of Gaussian illuminating beams with various degrees of focusing (see Figs. 1b and 1c), backward propagation of scattered signals and accurate description of their collection over the receiving aperture were performed. (a) is a cylindrical beam with invariable width and uniform phase and amplitude distributions as used in model [32]; (b) pronouncedly focused Gaussian beam used in models [33,34]; (c) a weakly focused Gaussian beam rigorously studied in Ref. [33] with accounting for the scattering-field divergence and high directivity of the signal reception; (d) the beam form, which is intermediate between the cases (a) and (c), such that the beam amplitude is approximated by a depth-independent Gaussian profile, whereas the phase is uniform in the cross section. This approximation is used in the following sections to reasonably well imitate weakly-focused beam in angiography-related simulations.
The result obtained in Ref. [33] demonstrated that, for a weakly focused beam, only a slight modification of the illuminating beam form in model [32] already gives OCT scans very similar to the results of rigorous consideration of a weakly diverging Gaussian beam with the same width in the focal region. Namely, instead of uniform amplitude distribution as in Fig. 1a, one may use a beam with a depth-independent Gaussian amplitude profile and uniform distribution of phase as in Fig. 1d. It is interesting that similarly to the assumption adopted in simplified model [32], for a rigorously simulated weakly focused Gaussian beam, the amplitude of the pointspread function demonstrates a very weak dependence on the depth (see examples of such a comparison in Ref. [33], Figs. 2a and 2b). The reason for this nearly depth-independent amplitude is due to interplay between the influence of divergence of scattered waves and their collection over the receiving-lens aperture, see details in Ref. [33]. Therefore, the modification of the beam as shown in Fig. 1d in model [32] gives results very close to those for rigorous consideration of OCT-scan formation for weakly focused Gaussian beams. Bearing in mind that such beams are often used in realizations of OCTangiography to enable nearly depth-independent lateral resolution, the so-modified model [32] can be efficiently applied for simulations of angiographic processing as demonstrated below.
As depicted in Fig. 1d, the illumination-beam intensity has a Gaussian profile : where characteristic radius is assumed depthindependent and determines the lateral resolution. Next, the illuminating-beam spectrum in terms of wavenumbers in the axial direction is described by the amplitude profile having spectral components , summation of which forms an A-scan consisting of pixels with axial coordinates The distance between the spectral components determines the total unambiguously imaged depth . If the characteristic spectral width of the spectral function is , the coherence length determining the axial resolution in the resultant scans can be estimated as or equivalently in terms of the central wavelength and spectral width , (note that more accurate coefficient in the estimate for depends on the exact form of the spectral function ). Then for an A-scan corresponding to lateral coordinates of the illuminating beam, the contribution of a scatterer with coordinates to the received signal at the spectral component with the wavenumber has the following form (2) where is the scattering strength of the scatterer for the spectral component , although often it may be sufficient to consider frequency-independent scattering strength. Now having the complex-valued amplitude of individual spectral components from a single scatterer, the expression for the amplitude of th pixel corresponding to the depth in the A-scan with coordinates can obtained by performing summation of contributions of all scatterers for all spectral components (3) The summation of spectral components on Eq. (3) is equivalent to performing the inverse discrete Fourier transform over the axial wave-vectors , such that the expression for a pixelated A-scan can be represented as (4) Simulations of 2D B-scans and 3D volumes of OCT data is performed by scanning the illumination beam, which corresponds to variation in coordinates of the illuminating beam in the model equations. It is clear from Eqs. (1) and (2) that the model readily describes how OCT-scan varies when scatterers gradually leave and enter the illuminating beam when its coordinates or coordinates of of the scatterers are varied. As mentioned above, if necessary, the depth dependence of the illuminating beam (for example exponential decay with the decay coefficient as in Ref. [39]) or acquisition noise at the receiving array can also be readily imitated by adding random complex-valued terms to spectral amplitudes (2) (as in Ref. [40]).
Expressions (1) -(4) can be evaluated very rapidly, because all signal amplitudes are given in a compact analytical form, so that only summation of contributions of all spectral harmonics and illuminated scatterers should be made numerically. Therefore, because of intrinsic ease of accounting for variations in coordinates of scatterers, the above presented Eqs. (1) -(4) represent a flexible and computationally very fast platform for simulations of evolution of OCT scans due to motion of scattering centers. In simulations related to OCT-based elastography (as in Refs. [38][39][40][41]), for which a couple of simulated B-scans usually was sufficient, the computational efficiency was useful, but not a critical condition. However, for modeling angiographic processing or diffusion phenomena, where large sequences of OCT-scans (in 2D and, moreover, 3D cases) are required, high computational efficiency of the used simulation method becomes critically important. In the following section we demonstrate how the abovepresented model Eqs. (1) -(4) are applied for simulation of angiographic processing using the high-pass filtering principles suggested in Ref. [6,17] and correlation processing discussed in Ref. [45].

Simulation of micro-vasculature visualization using the high-pass filtering principle
In this section we demonstrate the application of the above-formulated model for simulation of OCT-based visualization of microcirculation. There are somewhat different realizations of OCT-based label-free angiography that are based on comparison of sequentially acquired OCT scans. The regions of "solid" tissue in the images correspond to (nearly) invariable parts of OCT scans, whereas the vessels with moving blood particles demonstrate local variability that can be singled out by special processing of the compared images. Different angiographic algorithms proposed during the last two decades use either phase information Doppler methods, or amplitude variability like in the speckle-variancebased methods (see, e.g., review Ref. [46]), or utilize mixed phase-amplitude processing. An important problem in all such method is the influence of masking living tissue motions (related to heart beats, breathing, etc.) that should be somehow eliminated using physical immobilization, special processing of the acquired images or combinations of hardware and software methods. In some methods a series of B-scans are acquired in one position and then a similar series is obtained in another position and so on, such that a 3D stack of OCT data is finally obtained [3].
There are also variants, in which the scanning is performed continuously, but with a significant overlap of neighboring scans. From the viewpoint of controlling the beam position, continuous scanning can be easier realized in comparison with step-wise jumps. This idea can be realized in somewhat different forms. One variant is acquiring dense B-scans with significantly overlapped A-scans (say, the inter-A-scan step is ~50 times smaller than the illuminating-beam diameter) as proposed in Ref. [17]. Then 3D volumetric data can be composed from a series of such dense B-scans. In another variant described in Ref. [6], B-scans have "normal" density of A-scans, but in the 3D-data volumes these B-scans are acquired with an appreciable overlap (say, the inter-Bscan step is 6-7 times smaller than the illuminating-beam diameter). Since the the inter-A-scan interval in the first variants is significantly smaller than inter-B-scan time interval in the second variant, the latter method can easier detect flows with lower velocity. However, by the same reason of higher sensitivity to displacements the inter-Bscan method is less tolerant to masking motions of living tissues. To reduce the magnitude of such motions, usually this method has to be used in contact mode [6].
Schematically this approach based on continuous scanning is shown in Fig. 2. It is clear that in the direction of overlapping, motionless scatterers in the "solid" tissue correspond to essentially elongated speckles with fairly slowly varying amplitude and phase, whereas the speckles corresponding to the places with moving scatterers (vessels with moving blood particles) demonstrate faster variability (and look as "short" speckles).
An example of a real OCT record illustrating these statements is shown in Fig. 3a for a single dense B-scan obtained with a horizontal step between subsequent positions of the illuminating beam ~1/50 of its diameter. It is clear that in the region of the vessel with moving blood particles, the horizontal size of speckles is much smaller than for the regions of "solid" tissue with almost motionless scatterers. As mentioned in Ref. [17], this difference manifests itself in the Fourier spectrum of such OCT-data obtained by applying Fourier transformation along the direction of A-scan overlapping. If the so-found spectrum is subjected to high-pass filtering, then in its inverse Fourier transform the regions of motionless scatterers can be eliminated, whereas the zones of vessels with moving blood particles can be retained (see detail in Ref. [17]). In a similar manner the microvasculature can be imaged for a series of B-scans that are overlapped along the direction orthogonal to the B-scan planes [6].   3 demonstrates that application of high-pass filtering to the initial dense structural scan the regions corresponding to the "solid" tissue with fairly slow horizontal variability of the speckles makes it possible to retain predominantly the zones corresponding to perfused vessels. It is important to point out that besides motions of particles in perfused vessels, for living tissues motions of other nature (related to heart beats, breathing, etc.) are typical, the velocities of particles due to these motions may be comparable or even exceed velocities of scatterers typical for microcirculation of blood. Fortunately, such motions are usually characterized by significantly larger scale of spatial inhomogeneity being closer to bulk translational motions, the axial component of which being strongly dominant in the masking signal. This axial component of masking bulk axial motions can be efficiently suppressed during the angiographic processing by determining and compensating the depthaveraged phase difference between neighboring A-scans caused large-scale bulk motions [17]. Such compensation weakly affects small-scale motions of blood particles in vessels and significantly improves the quality of microvasculature visualization. When obtaining the image of the vessel cross-section in Fig. 3a-2 (by filtering the initial structural scan in Fig. 3а-1), the above described procedures of bulk-motion compensation were applied.
The quality of microvasculature visualization using methods based on particle-motion contrast (not necessarily for high-pass filtering) strongly depends on the parameters used during obtaining the OCT scans, as well as during their processing. Optimization of these parameters is non-trivial and in a rather complex ways depends on the characteristics of the OCT system, regime of scanning, blood-flow velocities, vessel diameters, intensity of masking motions, noises, etc. Sufficiently precise control of such a broad variety of various parameters in real experimental conditions is challenging or may even be impossible. In this context the presented in the previous section method allows one to perform Fig. 3 Elucidation of the high-pass filtering principle of blood-flow visualization. Panel (a-1) is a real example of a dense B-scan with significant overlapping of A-scans along the B-scan plane. In this scan for the motionless tissue, speckles are pronouncedly elongated in the horizontal direction, whereas for the zones of blood flow, speckles are much shorter. This difference also manifests itself in the spatial spectrum found along the horizontal direction for this B-scan. The spectrum characterizing the variability of speckles in the direction of scanning for the middle depth of the scan (a-1) in presented in panel (b-1), where the most intense low-frequency components correspond to elongated speckles within the "solid" tissue. Panel (b-2) shows the spectrum after elimination of the low-frequency components and panel (a-2) shows the structural scan obtained by inverse Fourier transform applied to the spectrum after high-pass filtering. Total time Tscan of the analyzed scan determines the spectral resolution of the Fourier transform. Below we illustrate such possibilities of the described model by generating dense B-scans similar to real scans shown in Fig. 3a. In the simulated region a blood vessel with a Y-like bifurcation was modeled. It was possible to vary the vessel diameter, its orientation, velocity of scatterers inside the vessel, to change the velocity profile in the vessel cross section. Also masking bulk motions of various intensity were simulated by adding vertical displacements of entire vertical columns of scatters between subsequent A-scans (note that reception noises of arbitrary magnitude can readily be introduced as was demonstrated in Ref. [40] in the context of realization of mapping strains). This flexibility of the model allows one to estimate ultimate possibilities of the angiographic visualizations (minimal detectable diameter, minimal velocity) and to estimate maximal level of noise, and degree of A-scan overlapping, for which the visualization is feasible.
In Fig. 4 we illustrate the utilization of the models for assessing the role of masking bulk motions and efficiency of their compensation by correcting the depth-averaged phase shift between neighboring A-scans. In the simulation we imitated the parameters of OCT setup used in Ref. [17]. Namely, the optical wavelength in the tissue was λ0 = 1 µm (corresponding to λ0 ~ 1.3 µm in air), the spectral width of the illuminating beam was about 90 nm, the beam radius was 15 μm, the vertical size of the pixel was 8 µm in air (about 6 µm in the tissue with the refractive index ~1.3). The simulated image depth was 600 µm and horizontal size 150 µm with the horizontal step 1 µm between the subsequent A-scans. The signalto-noise ratio was 23.5 dB in the upper part of the scans and depth decay similar to that for real scans shown in Fig. 3 was introduced. The vessel width before bifurcation was 15 µm and 7.5 µm after bifurcation. The velocity of blood flows was chosen to be 1 mm/s. The simulations were performed for three levels of bulk noises (indicated in the caption). For the upper row, the amplitudes of the inter-A-scan displacements are smaller than λ0/4 for all A-scans, so that the phase shifts induced by the masking bulk motions are unambiguously related to displacements. Therefore, in this case the masking displacement-related noises well-visible for scans in the middle column in Fig. 4 can be nearly perfectly compensated. This is clear from panel (a-3) in which the motion-related noises are very well suppressed. For rows (b) and (c), the inter-A-scan bulk motions mostly have supra-wavelength amplitudes, so that the corresponding phase variations are wrapped and are not unambiguously related to the displacements. Nevertheless, even in these cases, the application of the compensation based on comparison of phases of neighboring A-scans significantly improves the quality of visualizations of the vessels (see the right column in Fig. 4). Fig. 4 Examples of simulations performed in the context of the development of OCT-based angiography. Left column shows initial densely recorded scans, middle column is for filtered scans without correction of the bulk motion, and (right column) with correction of the bulk motion. Rows (a) through (c) correspond to inter-A-scan bulk motions with a Gaussian distribution of the amplitude with various standard deviations (STD). Upper panels (a1-3) correspond to translational displacements with STD = λ0/10, for which the amplitudes of the displacements are under λ0/4 and the inter-A-scam phase difference can be found unambiguously. Panels (b1) -(b3) are for significantly greater displacements with STD = λ0 and (c1) -(c3) are for STD = 1.5λ0.

Model application for simulations of OCT-signal decorrelation in the presence of regular and random scatterer motions
As mentioned above, for perfused vessels, the main contribution to the variability of received OCT signal and, correspondingly, to its temporal decorrelation belongs to motion of erythrocytes. Among other methods applied for detection of perfused microvasculature studies of temporal decorrelation of OCT signals attracts significant attention [47,48]. Although erythrocytes are discrete scatterers that are not yet resolved in conventional OCT images, their size exceeds the optical wavelength. Consequently, even in the absence of translational flows, rotational motions of erythrocytes may give an appreciable contribution to the OCT-signal decorrelation. Real experiments do not allow one to selectively study contributions of transversal and rotational motion of finite-size scatteres, whereas in simulations such studies can readily be implemented. To this end, we made the following modification of the scatterer model as shown in Fig. 5a. Fig. 5a schematically shows the cross section of an erythrocyte with an approximately dumbbell shape. Thus a reasonable simplified approximation of scattering properties of such an object can be a dumbbell made of a pair of localized scattering centers with a distance between the centers corresponding to the erythrocyte diameter. In the below presented simulations we chose µm (which is close to the average diameter of human erythrocytes). The central wavelength is supposed to be µm (which is close to the wavelength in the tissue for widely used OCT systems based on the sources with the wavelength 1.3 µm in vacuum). Using such a model for generation of 2D B-scans we characterize the position of the dumbbell center by coordinate , whereas  Fig. 5b.
Here, we consider three basic types of motions of erythrocytes in vessels: (i) regular translational flow, (ii) translational Brownian motion and (iii) rotational Brownian motion. Let us describe in more detail how these types of motions are characterized in the model. In real experiments it is impossible to realize these motion types in pure form, but can be readily made in simulation.
(i) For the translational flow-type motion, we assume that the erythrocytes move without changing their orientation ( ) with an invariable-in-time velocity. Then for the neighboring A-scans with numbers and the scatterer coordinates are related as , where is the flow velocity and is the time interval between A-scans.
(ii) For the translational Brownian motion, the orientation angle for the scatterers is also constant ( ), whereas the coordinates of the centers vary from one A-scan to another one in random steps, , where the components of the incremental vector have a Gaussian probability distribution for displacements along coordinates and (5) where the standard deviations were supposed equal, . (iii) For the Brownian rotational motion, coordinates of the scattering-particle centers remain constant, whereas the orientation angle randomly varies where the rotation-angle increment has a Gaussian distribution with a standard deviation (6) In all cases the initial position and orientations of the scattering particles are random.
For studying decorrelation of OCT scans due to such motions of scatteres, it was sufficient to form sequences of B-scans (for the below presented examples, scans 255 × 80 pixels in size were used). The scatterers with random initial coordinates were distributed uniformly over the simulated scans and experienced either one of the above-described three types of motions. The crosscorrelation coefficients were calculated for complexvalued amplitudes of pixels in the initial scan for and subsequently simulated B-scans for using the following expression: (7) According to previous studies, exponential dependences of the form could be expected, where parameter depends of the particular character of the scatterer motions [45]. It is clear that if we plot the logarithm against , the assumed exponential dependence should correspond to a straight line, which can be visually clearly seen for the correctly chosen parameter . . The results of the corresponding simulations for the above-mentioned three characteristic types of scatterer motions are presented in Fig. 6. The horizontal axes in Fig. 6 represent rather than so that even if varies nonlinearly with time for , the time is explicitly indicated in the time-axis labels. The dependences found in the performed simulations for regular translational flow demonstrated the exponent (Fig. 6a) and for translational Brownian motion (Fig. 6b), which agrees with the previously known results based on analytical considerations of these motion types [45,49]. Concerning rotational Brownian motion, for which analytical consideration is more difficult, the performed simulations demonstrate different regimes for the cross-correlation decay in the beginning and at larger times. Fig. 6c shows that initially (within the interval below ~25 ms) the decay law corresponds to like for translational Brownian motion, whereas for larger times, the transition to the scaling with is seen within a fairly broad time interval (Fig. 6d). As far as we know, the decay law of the form earlier has been found for speckles when the contribution of multiple scattering was significant. However, our simulations demonstrate that such functional behavior can also occur for speckle decorrelation in the case of single scattering for Brownian rotation.
The fact of initial (for time below ~25 ms) functionally similar behavior of speckle decorrelation due to Brownian rotation and translational Brownian motion may be qualitatively understood as follows. For randomly oriented dumbbell-like scatterers, there is a portion of approximately horizontally-oriented ones that produce in-phase scattered signals. This portion of initially in-phase scattering horizontally-oriented elements gives a dominant contribution to decorrelation. For this portion of dumbbell elements and sufficiently small deviation from the initial approximately horizontal orientation, the difference of optical paths for the signals scatterers by the left and right scatterers is linearly proportional to the rotation angle as , where the factor 2 is due to forward-backward propagation for the incident and backscattered waves. Quantity accumulates in random steps obeying a Gaussian probability law Eq. (6). This is similar to translational Brownian motion, for which the differences in the optical paths of interfering scattered fields also grows with a Gaussian probability, such that (see Eq. (5)).    Fig. 6d agrees with this expected transition to a slower law, although for predicting its functional form, the above-presented arguments are not sufficient. However, they allow one to estimate the characteristic angle of this transition as This simple estimate fairly well agrees with the mean-square rotation angle This value was found in the simulations for the time point ms, around which the transition from to is observed in Fig. 6d. It has also been numerically verified that if parameters , and are varied in such a way that the factor the transition time visible in Fig. 6d also remains invariable. This numerical result confirms the expectation that the characteristic time of the transition is determined by the product Concluding this section it can be mentioned that the time dependence of the correlation function for translational regular flows and Brownian motion of spherical particles was studied analytically in Ref. [49]. Such analysis predicts the decay law for the absolute value of the correlation coefficient of the form (8) where is the characteristic decorrelation time due to Brownian motion, and is the decorrelation time due to translational flow. It is clear that these known cases well agree with the above-presented numerical results, which confirms the correctness of the formulated numerical model. The influence of random rotations of erythrocytes on the decorrelation of speckles was touched in Ref. [50]. Analytical investigation was not performed in that works, and experimentally the decay law of the type was observed. This observation agrees with the decay law for sufficiently small times revealed in our simulations as shown in Fig. 6c.

Conclusions
Here we presented a computationally very efficient and fairly realistic model of OCT-scan formation in spectraldomain optical coherence tomography. The described variant of the model is based on previously proposed variants [32][33][34]. The described variant represents a compromise between the simplest version [32] based on cylindrical illuminating beams with uniform profile and rigorous analysis of OCT systems with fairly focused Gaussian beams [33,34]. The considered compromise variant uses beams with a Gaussian amplitude profile and uniform phase distribution over the cross section, which enables computationally efficient and reasonably realistic simulations of OCT-scan formation for weakly focused Gaussian beams, which is validated by comparison with the rigorously described Gaussian beam in study [33].
Similarly to the previous versions [32][33][34] and many other popular models (e.g. Ref. [30), the present model is based on the approximation of discrete scatterers and ballistic character of scattering. An important advantage of this model demonstrated above is the intrinsic ease in accounting for arbitrary scatterer motions, which allows one to computationally efficiently generate large sequences of OCT scans for gradually varying configurations of scatterers. This makes the described simulation platform very convenient for studies related to the development of angiographic processing of OCT scans for visualization of microcirculation of blood, as well as for studies of decorrelation of speckle patterns in OCT scans due to random (Brownian type) motions of scatterers as was demonstrated in the previous sections. To the best of our knowledge, such numerical simulations of large series of OCT scans in the presence of various types of motion of scatterers have not been demonstrated before. The model is also very useful for simulating images of deformed tissues in the context of the development of Optical Coherence Elastography [37][38][39][40]. If necessary, accounting for non-ideally ballistic propagation of the optical waves can be added as discussed in Ref. [36].

Disclosures
All authors declare that there is no conflict of interests in this paper. supported by RSF grant No 17-72-20249.