A computationally highly efficient full-wave model of OCT-scan formation by focused beams in spectral-domain optical coherence tomography (OCT) is presented. Similarly to some previous models, it is based on summation of fields scattered by discrete sub-resolution scatterers and enables accounting for axial and lateral inhomogeneity of the illuminating Gaussian beam. The main feature ensuring high computational efficiency of the described model is that instead of numerical integration of the scattered signal over the receiving aperture we apply analytical description of both the illuminating-beam focusing and collection of the scattered signals over the receiving aperture. Elimination of numerical integration over the receiving aperture has increased the computation speed by a factor of ~310. This is of key importance for practical feasibility of simulations of 3D OCT data volumes for large amounts (~751010) of scatterers corresponding to realistic densities of cells in biological tissues. We demonstrate the model possibilities by simulating digital refocusing of strongly focused OCT beams in the presence moving scatterers. A novel principle of contrast-agent-free visualization of scatterer flows with velocities typical of blood microcirculation is demonstrated. Keywords: OCT image formation; simulation of OCT scans; digital refocusing; lateral resolution improvement; flow imaging 1. Introduction Development of analytical and numerical models of scan formation in optical coherence tomography (OCT) has been attracting much attention over the last two decades since the first works [1,2describing formation of OCT scans including their peculiar speckle structure. The latter sometimes is considered as a negative factor degrading the image quality, so that various methods of speckle suppression has been discussed in literature(e.g., [3,4] and refs therein). However, analysis of speckle structure can also be used as an informative source in the context of the development of multimodal OCT, in particular, for elastographic [5,6, 7] and angiographic applications [8, 9]. A popular trend is the development of OCT-scan modeling is the Monte-Carlo approach [10,11,12]. The Monte-Carlo methods are usually rather computationally demanding (since they intrinsically need to simulate pseudo-random trajectories of huge amount of photons in turbid media). Besides, such an approach is not well-suited for generation of numerous data volumes corresponding to evolution of OCT scans due to motion of scatterers, although the latter is indispensable for testing signal processing in OCT elastography and angiography, where the locations of scatters are random, but their motions may correspond to either deterministic or random laws. There are also simulation methods in which the propagation/scattering is described as a regular process (although the medium is turbid). Such methods may consider fluctuations of refractive index with full-wave numerical simulation of the beam propagation and scattering [13], which is also very demanding computationally. There are also models based on utilization of discrete scatterers either using explicit representation of point-spread function for sub-resolution scatterers [1415] or obtaining the point-spread function in the axial direction by summing the spectral components similarly to what is made in spectral-domain OCT devices 16, 17]. The approach [16] based on summation of signal from discrete sub-resolution scatterers is very computationally efficient and is very convenient for generating series of OCT scans formed by arbitrarily moving scatterers, which has been efficiently used for simulating elastographic processing in OCT (e.g., [18, 19, 20, 21]). However, only weakly focused OCT beams with depth-independent radius could be described by model [16], so that its extension for describing the influence of pronounced beam focusing/defocusing was proposed in [22]. In [22] the influence of the divergence of the Gaussian illuminating beam was rigorously described analytically. The spherically diverged backscattered signal was convoluted with the receiving aperture having the same phase- amplitude mask that initially forms the incident illuminating beam. This collection of the baskscattered signal over the receiving aperture in [22] was made using numerical integration which is rather computationally demanding. Nevertheless, in comparison with Monte-Carlo methods or finite-element full-wave approaches, model [22] is significantly less demanding computationally and makes it possible to simulate a series of 2D B-scans on intervals of several minutes to several tens of minutes (depending on the number of scatterers) on CPUs of Intel-I7 class 2 without the necessity of multi-core GPU-based calculations. However, even for model [22] that is computationally fairly low-demanding, the required computational times increase up to several hours and even days for generation of a set of 3D OCT scans. This is indispensable, for example, for modeling digital refocusing of OCT data acquired using strongly focused illuminating beams. In what follows we present a computationally much faster modification of the semi-analytical full-wave model [22], in which the same basic principles are used to describe the illuminating beams and the signals scattered by localized sub-resolution scatterers. We emphasize that similarly to the case real tissues the scatterers may move along either regular or random trajectories, may demonstrate dependence of the sсattering properties on the wavelength, etc. The key difference from [22] is that the integration of the backscatterered spectral harmonics over the receiving aperture in the below described model is made analytically rather than numerically and only summation of the so-found complex-valued spectral amplitudes should be performed numerically for all scatterers. Due to elimination of numerical integration the computational efficiently of the model described in this paper is strongly increased (by a factor of 3210...10 depending on the number of scatterers) in comparison with [22]. In this paper the new computationally efficient form of the model is applied to simulate 3D volumes of OCT data consisting of hundreds of B-scans in the presence of both static and moving scatterers, including flows imitating perfused blood vessels. The realized acceleration of simulations opens previously inaccessible possibilities to perform highly controllable and realistic massive numerical experiments on formation of OCT scan in various conditions. In particular, we demonstrate the application of the developed model to studying refocusing procedures that attract much attention in recent years due to possibilities to obtain increased lateral resolution in OCT over the entire imaged depth [23]. Based on the simulations of digital refocusing of strongly focused OCT beams in the presence moving scatterers we demonstrate a novel principle of contrast-agent-free visualization of scatterer flows with velocities typical of blood microcirculation. 2. Optimized model of OCT-signal formation by a strongly focused beam In model [22] that is taken here as a basis, the focused illuminating field is described as a Gaussian beam incident onto the imaged tissue through an immersion layer as shown in Fig. 1. Fig. 1. Schematically shown focused illuminating Gaussian beam and a sub-resolution scatterer with coordinates sr producing a spherically-diverging scattered wave. Here, 0W is the beam radius in the focal plane and )(zW is its current radius, dz in the distance of the receiving aperture from the tissue surface. In spectral-domain OCT an elementary step consists of formation of a single one-dimensional A-scan in the direction of z-axis aligned with the probing-beam axis. The complex-valued amplitude )(qzAof qth pixel in an A- scan consisting of N pixels with the coordinate of the centers equal to zq is given by summation of all spectral harmonics with wavenumbers nk received by the array of photodetectors after propagation back from all scatterers with coordinates sr illuminated by the OCT beam:   sq ns nq zHnikrB zA )2exp(, )(. (1) Here the visualized depth H is determined by the distance between the neighboring spectral components | |/1nnkkH   . The complex amplitudes nskrB, of the received spectral components are jointly determined by the forward propagation of the illuminating beam, backward propagation of the spherically diverging fields scattered by the sub-resolution scatterers, as well as the collection of those back-propagated scattered waves over the receiving aperture (shown in Fig. 1 by the dash-dotted line). Equation (1) for the complex-valued pixel amplitudes can also be represented as 3        sns k q krB FFTzA n,1 , (2) where 1 nkFFT denotes the inverse Fourier transform applied to the received complex-valued amplitudes . Let us recall integral expressions obtained in [22] for the illuminating beam amplitude nskru, and the signal nskrB, collected at the receiving aperture. For a Gaussian beam incident at a localized (sub-wavelength) scatterer with coordinates ssss zyxr ,, , the complex-valued amplitude krus, of a wave with wavenumber k can be written in the form:                             02 0 0 0 022 0 00 ,2, exp ,exp,,zzRrrk zzzzki zzWrr zzWWkSkrus s s ss ss . (3) Here, the focus position is  000 0 ,,zyxr , where 0z is the axial coordinate and the lateral coordinates correspond to radius-vector 00 0,yxr . Quantity 0,zzW determines the current amplitude at the axis of illuminating beam and has the meaning of its current radius    2/12 2 00 00 ] 1[ ,     WzzWzzW , (4) where 0W is the beam radius in the focal zone; k/2 is the wavelength. Quantity   )/( ,2 0 0 0 Wzz arctgzz  is the phase difference between the Gaussian beam and an ideally plane wave front; the current radius of wave-front curvature in the beam is given by     ] 1[ ,2 02 0 00     zzWzzzzR , (5) where 0,zzR tends to infinity in the focal region corresponding to 0zz. Factor )(kS in Eq. (3) describes the spectral form of the probing beam. In principle, function )(kS may be considered to be jointly determined by the spectral properties of the source and the receiving system. The field ),,(krrus scat scattered by a sub-resolution scatterer with coordinates sr and scattering strength kKs corresponds to a spherically diverging wave of the form     rrki rrkrukKkrrus ss s s scat exp1, ),,( (6) Collecting over the receiving aperture the scattered field after its backward propagation and passing through the optical system can be written as the following integral expression for the sought complex-valued amplitude krBs, of the received signal for each spectral component nkk and each scatterer:      ddzzRrrk i zzWrr zrrki rrkrukKkrB dd dd dds dss Ds s                  02 0 022 0 ,,2exp ,exp exp1, , , (7) where    0 0,yxrd is the radius-vector of the integration point on the receiving aperture (shown by the dashed line in Fig. 1), where the integration coordinates , are counted from the coordinates ),(00yx of the illuminating-beam axis; vector dr has the following meaning,  d d zyxr ,,0 0  , where dz is the axial coordinate of the receiving aperture. The last two exponential factors in the integral (5) describe transformations of the amplitude and phase of the scattered signal during its backward passing through the optical system that initially has formed the illuminating Gaussian beam. In study [22] integral expression (7) was evaluated numerically, this evaluation being the most time-consuming procedure is the model formulated in [22]. To accelerate evaluation of integral (6) without appreciable loss in the accuracy, the following observations can be made. Taking into account that ( 0dz, 0sz) we consider the scatterers satisfying the condition,  2 0Wzzds , so that it is possible to approximately represent the argument in the first exponential function in Eq. (7): 4   ds ds s dds zz rrkkzzrrk       22 (8) Then integral (7) can be represented in a simplified form                 ddzzRki zzWzki yy xxzzkizzkrukSKkrB d dDs s s dsdss ss s             022 0222 ,2 02 0,2exp ,exp2exp1, , (9) For sufficiently large receiving aperture in comparison with the illuminating-beam cross section, simplified integral (8) can analytically evaluated as:                      2 02 0 222 02 0 4exp2exp1, , yyxx zzPkzkiyyxxzzki PzzkrukSKkrB s s dss s s ds dss ss s , (10) where the factor P is given by:      ds d dzzki zzRki zzWP 2,2,1 0 02 (11) Certainly, the exponential decay of the signals propagating though the tissue is not yet accounted in Eqs.(3)-(11) and if necessary should be introduced additionally. Overall, in comparison with the initial form proposed in [22] (where integral (7) was numerically evaluated), the analytical representation (10) of this integral accelerates computations required to simulate 3D volumes of OCT data by a factor ~102-103 (depending on the number of scatterers). This significant acceleration opens previously inaccessible possibilities for generation of large amounts of 3D data corresponding to the densities of scatterers typical of densities of cells in real biological tissues. Indeed, for other models (including [22] that was faster in comparison with the alternative variants mentioned in the introduction) simulations similar to the presented below would require either practically unacceptable time intervals (several days or greater), or the necessity of utilization of very powerful and expensive computing means to reduce the computation time. 3. Examples of refocusing of simulated OCT scans including the presence of moving scatterers Bearing in mind the interest to digital refocusing of OCT scans obtained with a highly-focused illuminating beam, an important application of the developed model can be its utilization for development of refocusing procedures and studying the effect of motion of scatterers on the results of digital refocusing. By analogy with [22] in what follows we will use the following spectral transformations corresponding to the beam refocusing [23]:                   kvuzi zyxAFFT FFT FFTikz zyxAkzvuyx yxvu kz4exp),,( )exp(,,22 0 ,,1 ,, (12) where zyxAz,, is the optical field corresponding to the numerical shift of the focal plane by a distance of z; yx, are the lateral coordinates and zis the axial coordinate (depth); like above FFTand 1FFT denote the forward and inverse Fourier transforms for the respective pairs of variables indicated in Eq. (13); zyxA,,0 is the complex-valued signal amplitude before the numerical shift of the focal plane; ),(vu are the lateral coordinates in the Fourier domain for the 2D-Fourier-transformed function zyxA,,0. For the below discussion, it will be sufficient to use Eq. (12) based on subsequent forward and inverse Fourier transforms, although other refocusing algorithms can also be used (e.g., convolution of the signal with a specially constructed function [23]). In what follows we consider discrete sub-resolution scatterers with frequency-independent scattering coefficient, which is a reasonable approximation for many real situations. We assume that the illuminating beam has a central wavelength of 1310 nm and spectral width ~90nm, which correspond to parameters widely used in real OCT applications. We assume that the lens plane, over which integration is performed, is located at 100dzm above the tissue. The beam focus is located at a depth of 128ozm inside the tissue. The 1/e radius of the beam in the focal waist is 30Wm for the beam with pronounced focusing like shown in Fig. 1. In the following examples we will show relatively small-size stacks comprising 128 B-scans along y-axis covering 192 m. Along the x-axis B-scan also covers 192 m and having 128 pixels, so that the horizontal inter-pixel step in 5 the (x,y) plane is 1.5 m in both directions. Each A-scan has an in-depth length of 128 pixels covering 512 m in air, so that the vertical pixel size is close to 4 m, so that the imaged volume is 192x192x512 microns in size. To demonstrate that the numerical estimation of integral (7) used in [22] and it analytical estimate given by Eq. (10) coincide with a high accuracy, first we consider images of seven well separated localized scatterers located along the vertical axis at depths zs=32, 64, 128, 200, 300, and 480 µm. The images of these scatterers before digital refocusing are shown in Fig. 2(a-1), where the left half of the images corresponds to the direct numerical estimation of integral (7) and the right half is for the analytical estimate (10) of this integral. The vertical profiles through the center of the scatterer images are shown in Fig. 2(a-2); they clearly demonstrate that the meaningful parts of these images coincide with a high accuracy and only in the dark region between the scatterers (~50-60 dB below the maximal intensities) there is some discrepancy which should be below the noise level in real images. Therefore, in the following examples we will demonstrate only the results based on the utilization of the computationally much more efficient analytical estimate (10). In Fig. 2(b-1) and 2(b-2) the result of refocusing via Eq. (13) over the entire depth is shown, from which it is clear that refocusing enables uniform lateral resolution similar to that in Fig. 2a for the scatterer located at the initial focal depth zs=z0=128 m. Figure 2b demonstrates, however, that the amplitudes of the refocused images of identical localized scatterers outside the focus have lower amplitudes than in the focus. To explain this fact we have to return to the structure of factor P that enters Eq. (10) for the received spectral amplitude krBs,. It can be noted that in Eq. (11) for factor P the first real-part term is related to purely geometrical variations in the amplitude of the illuminating beam. If this geometrical effect solely determined the received signal, then, performing refocusing via Eq. (12) by appropriately correcting the phase fronts described by the exponential factors with imaginary arguments in Eq. (10), it would be possible to collect the distance-dependent energy of the scattered signal and perform its in-phase summation over the receiving aperture. Consequently, for identical scatterers, the amplitudes of their refocused images would be identical. However, Eq.(11) for factor P also contains the other two imaginary terms. These terms are related to inconsistent structure of the phase in the received spherical-waves from the scatterers (with the curvature radius ||0zzs) and the curvature 0,zzRd corresponding to the phase transformation of a plane wave front passing through the optical system that collects the scattered signals over the receiving aperture. This optical system fully compensates the received-wavefrom curvature only for scatterers located in the focal plane, i.e. for 0zzs, so that the imaginary part in Eq. (11) disappears. For a scatterer with 0zzs, the spherical divergence of the scattered wave is not fully compendsated by the lens, so that summation is not perfectly in-phase over the receiving aperture and the resultant amplitude of refocused images is lower than for 0zzs. Therefore, in order to obtain identical amplitudes for refocused images of identical scatterers with various depths sz, the following correction factor corrK should be additionally applied to Eq.(10):     2 2 00 0 10      WzzPzzPKs s corr (13) The additional application of the correcting factor (13) equalizes the amplitudes of the refocused images of the localized scatterers as shown in Figs. 2(c-1) and 2(c-2). Note that in Figs. 2a-2c the noise is absent. It is clear that application of the correcting factor (13) also increases the level of the noise together with the useful signal. However, refocusing collects the noise incoherently, whereas the useful scattered signal is collected coherently, so that the signal-to-noise ratio after refocusing is increased, so that even the scatteres hidden by noises before refocusing may become visible after refocusing. This is demonstrated in Fig. 2(d-1), where the left part of the image is similar to the image in Fig. 2(a-1), but with additionally added noise (simulating the noise at the receiving array) with intensity ~-50 dB with respect to the intensity of the scatterer image near the physical focus at z0=128 m. It is clear that in the lower part of the left half of the simulated scan in Fig. 2(d-1) the images of scatterers are almost hidden below the noise. In contrast, the right half of Fig. 2(d-1) shows that after refocusing and application of the correcting factor (13) the scatterers become clearly visible against the noise and their amplitudes are nearly equal as expected for identical scatterers. After such verification that the results of simulations of OCT scans and subsequent application of refocusing procedures well correspond to the results expected from the physical considerations for motionless scatterers, we consider the influence of motion of scatterers. In the following examples instead of rarified scatterers we will increase their concentrations to make them comparable with the density of cells in real biological tissues, where we assume that the cells have characteristic diameters ~5-6 µm. This corresponds to effective densities of scatterers on the order of (4-6) x106 mm-3. For generation of 3D volumes of OCT data with such concentrations of scatterers, numerical evaluation of integral (7) over the receiving aperture may require unacceptably long calculations (up to many hours and days using a CPU of a desktop or would require specialized and expensive high-performance computational means to accelerate calculations). However, the acceleration due to analytical estimate (10) of the signal collected over the receiving aperture opens the possibility to perform such simulation using a “typical” desktop computer (the described below 6 results were obtained using 16-cores of a CPU Intel I7 6950X and required several hours to simulate 3D volumes of 192x192x512 microns in size containing 5102 scatterers). Figure 2. Examples of simulated OCT images (upper row) and the corresponding vertical profiles (bottom row) through a chain of seven identical localized scatterers placed at the depths zs=32, 64, 128, 300, 420 and 480 m. The focus depth is z0=128 m. The left half of panel (a-1) is simulated using numerical evaluation of integral (7) as in model [22] and the right half of panel (a- 1) is simulated using the analytical estimate (10) of Eq. (7). Panel (a-2) shows the vertical profiles corresponding to the numerical (circles) and analytical (solid line) estimates of integral (7) corresponding to the left and right halves of panel (a-1). Panel (b-1) corresponds to refocused image (a-1) using Eq. (12) without correction of amplitudes of the refocused images of scatterers; the corresponding intensity profile through the scatterers is shown in (b-2). Panel (c-1) is the same refocused image as in (b-1), but with additional correction factor (12) due to application of which the intensities of the identical scatterers have equal amplitudes as shown in the corresponding vertical profile in (c-2). The left half of panel (d-1) is similar to the unfocused image (a-1), but with the amplitude-correction factor (13) and added noise -20 dB with respect to the scatterer image at the focal depth, so that the out-of-focus scatterers are significantly masked. The right half of (d-1) is obtained by digital refocusing of the image in the left half of (d-1) also with the amplitude-correction factor (13). The vertical profiles corresponding to image (d-1) before (solid line) and after refocusing (dotted line) are shown in panel (d-2). From the above discussions it can be expected that refocusing of moving scatterers by collecting defocused images of scatterers via digital correction of the phase fronts (which is represented in Eq.(13)) should be differently affected by lateral and axial motions of scatterers. It is clear that by purely geometrical reasons lateral motions cause weaker phase distortions in comparison with axial motions with the same velocity. Note that by the very same reason the Doppler frequency shift of the OCT signal for axial motions of scatterers is much greater than for lateral ones. In OCT-based angiography usually the main interest represents visualization of microvasculature with diameters of vessels below ~50 microns, whereas the smallest capillaries may have the diameter comparable with sizes of red blood cells ~5-6 microns. The velocities of scatterers in such fairly thin vessels typically range from a few mm/s to fractions of mm/s. Since in real situations orientations of the vessels are not ideally orthogonal to the OCT-beam axis, usually there is an axial projection of the flow velocity (at least ~several per cent of the total velocity). In view of this in what follows we consider simulated examples, in which a vertical flow of scatterers is introduced to imitate a vessel with 10 µm radius. Figure 3 shows an example of such a simulation for the same conditions as for individual scatters in Figs. 2(a-c) (in the absence of additive noise), but for randomly distributed 5102 scatterers, the density of which is close to concentration of cells in biological tissues. In the simulation, the axially moved particles displaced by 2.25 x10-4 m between the consequent A-scans, which corresponds to the velocity 2.25 m/s for the acquisition rate of the OCT scanner described in [9,24] and applied for angiographic mapping in [25,26,27,28] using the high-pass filtering principle. Figure 3a shows the initial B-scan through the axis of the vertical flow and Fig. 3b shows the intensity of profile along the depth indicated by the dashed line in Fig. 3a, in which the of flow is not visible at all. Figures 3c and 3d show the same B-scan and intensity profile, but after application of digital refocusing over the entire visualization depth, so that the presence of the vessel becomes clearly visible. Indeed, it can be expected that, when collecting the scattered fields over the cross-section of the beam to perform refocusing, the phase distortions caused by motion of scatterers should perturb the inter-scan phasing and reduce the amplitude of regions with moving scatterers in the resultant refocused image. This effect should be stronger in out-of-focus regions, where the beam is significantly wider that in the focus. Consequently, the number of the overlapped B-scans, from which the refocused image of a scatterer is formed, is greater than in the focal plane. Unlike focus, where only the nearest B- 7 scans are partially overlapped, in out-of-focus regions the significantly larger stacks of overlapped B-scans are acquired during larger time intervals, so that the motion of scatterers causes stronger mutual de-phasing of those overlapped B-scans. During refocusing this de-phasing impedes in-phase summation of signals from moving scatterers, so that the flow region has significantly lower intensity. However, closer to the focus, where the illuminating beam is insufficiently wide, the motion-induced de-phasing is not yet able to visualize the flow in the refocused image. Figure 3. Examples of focusing/refocusing influence on OCT scans of 3D volumes, in which a cylindrical flow of vertically moving scatterers is introduced to imitate a vessel with 10 m radius. In this example the vertical displacement of the particles is 2.25x10-4 m between A-scans. For these examples, the physical focus of the illuminating beam with radius 30Wm is placed at a depth of 1280z m. Panel (a) shows the B-scan passing through the “vessel” axis before refocusing and panel (b) shows the intensity profile through the vessel along the dashed line in (a). Panel (c) is the refocused B-scan at the same position, where the vessel becomes clearly visible as a region of strongly reduced amplitude because of signal-phase perturbations caused by the motion of scatterers. Panel (d) shows the intensity profile along the dashed line in (c), where the intensity is strongly reduced below the focus (compare the encircled regions in panels (b) and (d)). Taking this observation in mind, in order to enable higher contrast of the so-visualized “vessel” over the desired depth range 5000z m it makes sense to place the focus somewhat deeper. The following Fig. 4 corresponds to such a case, where similarly to Figs. 2 and 3 the depths z<500 µm represent the main interest for imaging, but the location of the physical focus is intentionally chosen deeper (at z 0=628 µm), so that within the imaged region the beam is significantly wider than in the focal region. Such a configuration ensures sufficient sensitivity of refocusing to scatterer motions and also the illuminating-beam convergence may be helpful to compensate the signal attenuation due to absorption and scattering in real conditions. Figure 4. Examples of focusing/refocusing influence on OCT scans of 3D volumes, in which a cylindrical flow of vertically moving scatterers is introduce to imitate a vessel. In this example the vertical displacement of the particles is 2.5 x10-4 m between A-scans. For these examples, the physical focus of the illuminating beam with radius 30Wm is intentionally placed below the depth of interest at 6280z m. Panel (a-1) shows the B-scan passing through the “vessel” axis before refocusing and panel (b-1) is for the refocused B-scan at the same position, where the vessel becomes clearly visible as a region of strongly reduced amplitude because of signal-phase perturbations caused by the motion of scatterers. Panel (c-1) shows cross-sectional C- scan through the “vessel” at depth z=300 m (shown by the horizontal dashed lines in (a-1) and (b-1)) before refocusing and (d- 1) the same cross-section after refocusing. Panels (a-2) and (b-2) show the vertical profiles along the “vessel” axis shown by dashed lines in B-scans (a-1) and (b-1). Panels (c-2) and (d-2) show the horizontal profiles through the “vessel” along the dashed lines in C-scans (c-1) and (d-1). 8 Figure 4(a-1) shows the initial B-scan in which the “vessel” is not visible against the background material and Fig. 4(b-1) shows the refocused B-scan for the same position, where the “vessel” becomes clearly seen as a stripe with strongly reduced intensity by the above-formulated reasons. Figures 4(a-2) and 4(b-2) show the profiles of Figs. 4(a-1) and 4(b-1) along the vertical dashed lines, where the reduced intensity within the vessel’s cross section after refocusing is also clearly seen. Figures 4(c-1) and 4(d-1) are the horizontal C-scans plotted through the dashed lines shown in Figs. 4(a-1) and 4(b-1). In the C-scan before refocusing the vessels is not seen against the background similarly to Fig. 4(a-1), whereas after refocusing the vessel cross-section is well seen in the C-scan in Fig. 4(d-1). The corresponding intensity profiles through the center of the vessel’s cross sections in the initial and refocused C-scans are shown in Figs. 4(c-2) and 4(d-2), where after refocusing the vessel becomes clearly visible. It has been verified that although in refocused cross-sectional images (like in Fig. 4(d-1)) it was difficult to recognize an ultimately thin vessel with diameter of 6 m close to the diameter of individual erythrocytes, in the axial sections (like in Fig. 4(b-1)) trajectories of such thin vessels were still decently well visible, which confirms that the refocusing principle opens the possibility to resolve the ultimately thin capillaries. After refocusing and singling out the vessels as “dark channels” the so-constructed image can be inverted to represent the microvasculature network in the conventionally used way as bright channels against the dark background. Concerning the sensitivity of the refocusing-based method it has been verified that for visual contrasting of the flow in the refocused image it is sufficient that during the acquisition time of the entire B-scan-stack covering the cross-section of the defocused beam the vertical positions of the scatters should change by ~ 2/...4/ , so that the constructive interference of the summed B-scans during the refocusing procedure becomes significantly perturbed. For the discussed “typical” OCT system, this corresponds to the axial component of particle velocities ~ several m/s, which should be usually sufficient to visualize blood microcirculation. The additional advantage is that the lateral resolution due to focusing can be made significantly better than for systems with weakly focused beams having several times larger diameter. 4. Conclusions This paper describes an optimized version the full-wave model of OCT-scan formation [16] , in which the medium is represented as an ensemble of discrete scatterers. Similarly to the earlier variant [16] the present version of the model uses a Gaussian illuminating beam with allowance for pronounced focusing. The OCT scans are formed by collecting over the receiving aperture the field of every spectral component scattered by discrete sub- resolution particles. However, unlike [16] where the integral (7) over the receiving aperture was evaluated numerically, in this work we present an analytical expression (10) for the result of this integration, so that only summation of the analytically found expressions for every spectral component and every scatter has to be performed numerically. This elimination of the numerical integration allowed us to strongly increase the computation speed by ~102-103 times for the same computer configuration. Therefore, now simulations of 3D OCT data sets for regions containing ~105-106 scatterers can be performed on intervals on the order of several hours using a CPU Intel I7 6950X. For the typical density of cells in real biological tissues this corresponds to tissue volumes ~108-109 µm3, which is quite sufficient for a broad range of problems related to testing and perfection of signal-processing methods in OCT. The model is very convenient for modeling situations with arbitrarily moving scatterers allowing for flexible control of all main parameters characterizing the OCT setup and scattering tissue. This enables previously inaccessible ease in studying influences of various factors affecting the OCT-based imaging of biological tissues. Realization of similar highly controllable physical experiments would be very challenging – expensive, extremely laborious or even impossible in practice. In particular, in the above sections we demonstrated results of such highly controllable numerical experiments related to feasibility of numerical refocusing of OCT scans formed by highly-focused beams in the presence of motion of scatterers. The results demonstrated interesting prospects for visualization of particle flows with parameters similar to the motion of erythrocytes in blood microvasculature. Although the suggested method of vessel visualization similarly to other contrast-agent-free approaches is based on the use own motion of the scattering particles, the principles of the above described visualization is significantly different. The proposed computationally efficient modeling of large series of realistic OCT scans can be used to study in detail the proposed new angiographic approach, as well as can be used for a broad range of other applications where realization of physical experiments is challenging or even impossible in sufficiently controllable conditions. Acknowledgements The development of the full-wave model of OCT-signal formation was supported by the RFBR grant No 19-02-00645 and simulations of full-depth refocusing in the presence of moving scatterers were supported by the RSF grant No 17-72-20249.