Simulation of stimulated Raman scattering signal generation in scattering tissues excited by Bessel beams

1. Introduction

Label-free microscopy does not require special labeling for in vivo imaging and longitudinal observation of biological systems. As one of them, coherent Raman scattering (CRS) microscopy can image chemical bonds specifically, for example the coherent anti-Stokes Raman scattering (CARS) microscopy and stimulated Raman scattering (SRS) microscopy,^1,2,3,4,5,6 and has become an indispensable tool in biology and medicine. CRS is a third-order nonlinear optical process. The excited resonance signal is much stronger than that generated by spontaneous Raman scattering, so that the imaging speed is significantly faster.³ However, when a CARS signal is generated, it is usually accompanied by a non-resonant signal, which can cause spectral distortion and limit detection sensitivity.⁷ Compared to CARS, SRS has no nonresonance background and can automatically meet the phase matching condition, making it ideal for imaging biological systems.⁸ Traditionally, SRS is based on the principle that two Gaussian beams excite samples at the same time to produce SRS signals. The tightly focused nature of the Gaussian beam enables high axial resolution imaging. However, the focused status of the Gaussian beam is greatly affected by the scattering medium it passes through.⁹ Even after passing through a thin distance of the scattering medium, the focus of the Gaussian beam would be expanded or completely destroyed. This will greatly affect the quality of SRS images, and even make it impossible to image.¹⁰ Therefore, it is difficult for SRS based on Gaussian beams to image the chemicals in the scattering tissues.

Bessel beam, as a kind of nondiffractive beam, has a characteristic that is significantly different from that of ordinary beams. It can be refocused when encountering scattering obstacles.^11,12,13 Such focus self-healing feature has made a breakthrough in the applications of photonics-based biomedical imaging. In fluorescence microscopy, with the help of the Bessel beam, the problem of poor imaging quality of the Gaussian beam in scattering tissues is avoided.^14,15 Here, we demonstrate theoretically the generation of the SRS signal in scattering tissues excited by Bessel beams. The application of Bessel beams can increase and maintain the focal length of SRS imaging in scattering tissues. First, the scattering tissues are simulated by the fractal model of refractive index turbulence which is close to the scattering feature of real tissues.¹⁶ Second, the propagation of Bessel beam is simulated with the fractal propagation method (FPM) which is an implementation of beam propagation method (BPM) in fractal model.^17,18 In the end, the SRS signals of chemicals in scattering tissues stimulated by the Bessel beams are calculated by using the theory of SRS signal generation.² In our validation, a series of simulations is performed to demonstrate the concept of using the Bessel beams to achieve SRS imaging in scattering tissues, and to investigate the influence of fractal parameters on SRS signal generation. The change of the fractal parameters will affect the refractive index turbulence,^{19,20,21,22,23} and therefore affect the scattering characteristics of the simulated tissues. We also prove the superiority of Bessel beam-based SRS in scattering tissues by comparing with the Gaussian beam-based one. Simulation results show that the image quality of SRS imaging worsen with the increase of refractive index turbulence that will enhance the scattering ability of tissues. However, under the excitation of Bessel beams, the SRS signals can be generated in deeper scattering tissues.

2. Theory

The Bessel beam is a type of nondiffracting beam whose electric field generated by a pair of axicons combined with an objective can be expressed as^24,25:

where $E_0$ and $w_0$ are the electric field and the waist radius of the Gaussian beam incident on the first axicon; $z_f$ and $z_B$ are the focal length and the equivalent Rayleigh range of the generated Bessel beam; $r_c$ is the radius of the ring beam incident on the objective; $J_0(g)$ is the Bessel function with zeroth order; $k_0$ is the wave vector of generated Bessel beam; and $\alpha$ is the angle between the propagation direction and the radial direction of Bessel beam.

Commonly, the biological tissues can be modeled as aggregates of discrete spheres.²⁶ However, this model yields a homogeneous set of optical properties, which cannot accurately describe biological tissues with heterogeneous refractive index.¹⁸ In real biological tissues, the scattering of light can be caused by the variation of refractive index due to the irregular structure of organelles in biological tissues.²⁷ The spatial fluctuations of refractive index are fractal in nature, and thus the fractal model can be adopted to characterize the refractive index turbulence.^16,18 The power spectrum of refractive index turbulence in the fractal model can be described as¹⁶ :

where $l$ denotes the correlation length in tissues with a minimum and a maximum cutoff of $l_{min}$ and $l_{max}$; $S_A$ describes the overall fluctuation magnitude of refractive index; $D_f$ represents the fractal dimension of the modeled tissue, which dictates the shape of the power spectral density; $(k_x,k_y,k_z)$ are coordinates in power spectral domain. The power spectrum of refractive index turbulence can be converted into a synthetic refractive index distribution that is randomized in spatial domain. By multiplying Eq. (2) by random phases, the synthetic refractive index distribution can be expressed as : where $R$ is the random phase uniformly distributed at (0, 2$\pi )$.²⁸

FPM that is a combination of BPM and the fractal model, is used to simulate the propagation of Bessel beam in the scattering tissues. The FPM can be regarded as the BPM based on the stepped spectral-domain. In the method, the beam propagation is divided into three steps, and the refraction and diffraction of Bessel beam during propagation are handled separately at each step. The electric field of the Bessel beam after propagating a step $\mathrm{\Delta }z$, regarded as $E_{BP}(x,y,z+\mathrm{\Delta }z)$ can be determined by¹⁸:

where $E_B(x,y,z)$ is the initial electric field of Bessel beam incident on the scattering tissues; $(x,y)$ and $(k_x,k_y)$ are the coordinates in spatial domain and corresponding frequency domain, respectively; $n_0$ is the average refractive index of the surroundings; $\mathrm{\Delta }n$ is the refractive index fluctuation that denotes the difference between the scatters and surroundings and can be calculated by $\mathrm{\Delta }n=|n_0-n(x,y,z)|$; ${ℱ}$ and ${{ℱ}}^{-1}$ represent the operations of Fourier transform and inverse Fourier transform, respectively.

In theory, SRS is a third-order nonlinear optical process, which requires both a pump beam and a Stokes beam for signal stimulation. When the pump and Stokes beams are collinearly overlapped both in space and in time, the SRS signal can be generated. The generated SRS signal is proportional to the intensity of the pump and Stokes beams, and can be expressed as² :

where $E_{Bp}$ and $E_{BS}$ are the electric fields of the Bessel pump and Stokes beams calculated with Eq. (4).

3. Results and Discussion

In the experimental validation, a scattering tissue model was first generated using Eqs. (2) and (3), with the dimension of 300 × 300 × 600 and the grid resolution being 0.5 $\mu$m. The fractal parameters of the refractive index were selected according to the relevant literature,^16,18 and these parameters would have a positive impact on the refractive index of the scattering tissue model. To maintain the variance of the refractive index fluctuation being around 0.0004, the fractal dimension was designed to vary between 2.50 and 4.50, and the maximum cutoff length $l_{max}$ varies between 0.50 and 10.00. The wavelengths of the pump and Stokes beams were set as 800 nm and 1040 nm, which stimulated the SRS signals at the Raman shift around 2885cm${}^{-1}$. The radius of the central lobe was designed to be 0.64$\mu$m for the Bessel pump beam and 0.86$\mu$m for the Bessel Stokes beam, respectively.

We first investigated the generation of SRS signals in a medium with homogeneous refractive index, and compared them with those generated in a scattering tissue with heterogeneous refractive index. Each simulation was run 100 times and the average was taken as the final results. Figures 1(a) and 1(b) show the longitudinal distributions and intensity profiles of the pump, Stokes, and SRS signal intensity generated in the media with homogeneous. The oscillations of the peak intensity of both the light beams and SRS signals can be observed in the longitudinal intensity profiles. This is in line with the results in the existing literature²⁹ and may be due to the exchange of energy between the central lobe and the side rings. Figures 1(c) and 1(d) show the longitudinal distributions of the SRS signal intensity generated in the media with homogeneous and heterogeneous refractive index, respectively. Here and after, the media with homogeneous refractive index was recorded as the homogeneous medium, and the media with heterogeneous refractive index was defined as heterogeneous tissue for short. We find that the SRS signal intensity fluctuated to some extent, especially in the deep layer, in both homogeneous media and heterogeneous tissue. In the heterogeneous tissues, the deepest depth that SRS signal can produce was slightly smaller than that in the homogenous medium, but the strength of the SRS signal intensity was much less. In order to show these results more intuitively, we extracted the cross-sectional images at different z positions and the longitudinal profiles along the $z$ direction of the SRS signal intensity. The cross-sectional images were extracted at positions of z = 35, 70, 105, 140, 175, 210, 245, and 280$\mu$m in the homogenous medium [Fig. 1(e)] and in the heterogeneous tissue [Fig. 1(f)]. The longitudinal profiles were extracted along the central lobe of the Bessel beams, as shown in Fig. 1(g). Almost the same conclusion can be addressed from the cross-sectional images and the longitudinal profiles. From Fig. 1(g), we found that in homogeneous medium, the SRS signal intensity keeps stable excitation in the depth of 100$\mu$m, and then the signal strength will fluctuate. The amplitude of the fluctuation will gradually increase until it increases to the maximum value and then rapidly decreases to disappear, as plotted in black curve. On the contrary, in the heterogeneous tissue, the SRS signal intensity cannot be kept constant. From the generation of SRS signal, its intensity attenuates nearly linearly in a short depth, around 175$\mu$m in this investigation case. This approximate linear attenuation is due to the fact that the heterogeneous tissue is approximately linear in its attenuation of the light beam. The heterogeneous tissue is a three-dimensional volume describing the spatial distribution of the refractive index, and at each point in this volume the refractive index fluctuates somewhat slightly. These differences in refractive index can affect the properties of light transmission. However, if this three-dimensional volume is layered in the direction of light transmission, the spatial distribution of refractive index variations is approximately uniform at each layer, so that the attenuation of light intensity and SRS signals is approximately linear. Subsequently, the SRS signal fluctuates similar until the signal disappears, but the signal strength is much lower than the SRS signal generated in the homogeneous medium. The simulation results show that in the heterogeneous tissue, due to the existence of scatters, the strength of SRS signal is attenuated faster. This is consistent with our previous results of multilayered beads-based investigation.²⁵

Then, we investigated the influence of fractal parameters on the SRS signal generation, including the fractal dimension and the maximum cutoff length.¹⁸ Two groups of simulations were performed. First, the maximum cutoff length of $l_{max}$was fixed to be 0.9$\mu$m, and the fractal dimension of $D_f$ changed from 3.25 to 4.15. This leads to an increase of the fluctuation range of refractive index from [1.35, 1.38] to [1.34, 1.39],¹⁸ which indicates that the fluctuation magnitude of refractive index become larger. Second, the fractal dimension of $D_f$ kept as 3.25, and the maximum cutoff length of $l_{max}$ was selected as 0.9, 2, 3, 4.5 and 6 $\mu$m. In this case, the fluctuation of refractive index in generated tissue was increased from [1.35, 1.38] to [1.345, 1.39].¹⁸ In these investigations, five spheres with a radius of 4 $\mu$m were located on the central lobe of the Bessel beams with the position of $z=50$, 100, 150, 200, and 250 $\mu$m, respectively to mimic the chemicals that can produce SRS signals under the irradiation of the Bessel beams. Similarly, each group of simulation was run 100 times to investigate the statistical results.

Figure 2 shows the investigation results of the influence of the fractal dimension on the SRS signal generation, presenting the average results of 100 times simulations. The representative results of single simulation were showed in Supplementary Fig. S1. We drew the longitudinal [Fig. 2(a)] and cross-sectional [Fig. 2(b)] distributions of the SRS signal intensity, as well as the longitudinal profiles along the $z$ direction [Fig. 2(c)]. The cross-sectional images in Fig. 2(b) were selected at positions of z = 50, 100, 150, 200, and 250$\mu$m. We find that the results of single time simulation had a great fluctuation (as shown in Supplementary Fig. S1), and the average of hundreds of simulations can eliminate this kind of fluctuation. From the longitudinal and cross-sectional distributions, we find that with the increase of fractal dimension, the strength of SRS signal decreased gradually. However, the degree of such signal strength decrease was different at different depths. For example, at positions of z = 100, 150, and 200$\mu$m, the SRS signal strength decreased much more obviously than the positions of shallow and deeper depths, as shown in Figs. 2(a) and 2(b). In other words, the influence of fractal dimension on SRS signal generation was also related to the depth of the chemicals. The longitudinal profiles in Fig. 2(c) showed this phenomenon more quantitatively. To quantitate this difference, we calculated the attenuation ratios of the scattering tissue to SRS signal, and then drew them as a function of fractal dimension at different depths, as plotted in Fig. 2(d). We can find that the deeper the depth was, the larger the attenuation ratio could be. At all depths, the attenuation ratio increased with the increase of the fractal dimension. At deeper depths, however, the rate of growth of this attenuation ratio decreased. The fractal dimension determines the degree of fluctuation between adjacent voxels in the generated tissue. The larger the fractal dimension is, the greater the difference of refractive index between adjacent voxels in tissue will be, and the more severe the scattering will be. When the beam passes through the tissue with large difference in refractive index distribution, the energy attenuation of the beam will be large, resulting in the shorter propagating distance and finally a great impact on the generation of SRS signal. Our simulation results right confirmed this statement.

The investigation results of the influence of the maximum cutoff length on the SRS signal generation is shown in Fig. 3. Similarly, the results of the single time simulation have strong uncertainty (as shown in Supplementary Fig. S2). Although the fluctuation at different depths was not fixed, it generally showed a downward trend. We took the average results of 100 times simulations as the analysis object, and we find that very similar results to the previous simulation were obtained. The strength of the generated SRS signal, or the attenuation ratio of SRS signal, was related to both the maximum cutoff length and the investigated depth. When the maximum cutoff length was less than 3$\mu$m, the SRS signal was linearly attenuated in the depth range of 200$\mu$m. On the contrary, when the maximum cutoff length was larger than 3$\mu$m, the SRS signal had a relationship of approximately inverse square with the depth. On the other hand, with the increase of the maximum cutoff length, the strength of SRS signal decreased, and the attenuation ratio of such decrease was different at different depths, as shown in Fig. 3(d). It should be noted that compared with the previous simulation, in all of the investigated depths, the attenuation ratio of SRS signal generated with the increase of maximum cutoff length was larger than that generated with the increase of fractal dimension. Taking the depth of 200$\mu$m as an example, when the maximum cutoff length was 0.9 $\mu$m, the strength of the SRS signal attenuated to about 70% of that in homogeneous medium. However, when the length increased to 3$\mu$m, the SRS signal strength was only about 20%. This indicated that the maximum cutoff length had much significant impact on SRS signal generation.

Finally, we investigated the superiority of the Bessel beam-based SRS in the scattering tissues by comparing with the Gaussian beam-based one. In the simulation, the grid resolution was increased to 0.2$\mu$m in X and Y directions. The simulation environment was composed of the scattering tissue and target chemicals. The fractal parameters to generate the scattering tissue were set as: $D_f=3.25$, and $l_{\mathrm{\text{max}}}=3$$\mu$m. The target chemicals were designed as a sphere of 4 $\mu$m in radius and placed at the depth of z = 50$\mu$m. It was placed on the central lobe of Bessel and Gaussian beams to stimulate the SRS signals. The energy on the central lobe of Bessel beams were equal to the energy of Gaussian beam at focus, and the same numerical aperture was set to generated these beams. Figure 4 reports the relevant results. First, we found that, compared with the Gaussian beam, the Bessel beam had a longer excitation distance and maintained a more stable excitation over the entire length of the target chemical. Second, the intensity of SRS signal stimulated by Bessel beams was stronger than that generated by Gaussian beams. For example, at the center of the target chemical, the SRS signal intensity excited by Gaussian beams was only 30% of that excited by Bessel beams [Fig. 3(d)]. Third, there is a very interesting phenomenon that under the excitation of Gaussian beams, the highest SRS signal intensity was not in the center of the target chemical, but was slightly shifted to the shallow position, as shown in Figs. 4(b), 4(c) and 4(e). The maximum intensity of the SRS signal excited by Gaussian beams was around 50% of that excited by Bessel beams. Results of this simulation indicated that the scattering tissue had more influence on the focused Gaussian beam. In other words, Bessel beam had stronger penetration ability in the scattering tissue. Collectively, the Bessel beam-based SRS has better application potential for imaging target chemicals in the scattering tissue.

4. Conclusion

In conclusion, we theoretically and numerically proved the potential of the SRS microscopy based on Bessel beams in imaging of deep chemicals in the scattering tissues. The self-healing properties of Bessel beam enable the beam to travel longer distances in the scattering media compared to the Gaussian beam, while allowing the beam focus to be maintained well enough to stimulated the SRS signal in deep tissues. We first generate the scattering tissue using the refractive index turbulence model, and then establish a theoretical model to calculate the generation of SRS signals in it. By changing the fractal parameters of the modeled tissues, which could affect the refractive index turbulence of the scattering tissues, we find that the generation of SRS signal based on Bessel beams are greatly influenced. The larger the refractive index turbulence is, the more serious of the SRS signal attenuation will be. Even though the modeling parameters have great influence on SRS signal generation, the Bessel beams-based SRS can work in depth of the scattering tissues, and provide better signal performance than the Gaussian beam-based one. Future work will concentrate on systematic experiments to demonstrate these simulation results.

Conflicts of Interest

The authors declare that there are no conflicts of interest relevant to this article.

Acknowledgment

This work was supported in part by the National Key R&D Program of China under Grant No. 2018YFC0910600, the National Natural Science Foundation of China under Grant Nos. 81871397, 81627807, 11727813, 91859109, the Shaanxi Science Fund for Distinguished Young Scholars under Grant No. 2020JC-27, the Shaanxi Young Top-notch Talent of “Special Support Program”, the Best Funded Projects for the Scientific and Technological Activities for Excellent Overseas Researchers in Shaanxi Province (2017017).