⇦Chinese Physics Letters, 2016, Vol. 33, No. 11, Article code 114301⇨ A Modified Monte Carlo Model of Speckle Tracking of Shear Wave Induced by Acoustic Radiation Force for Acousto-Optic Elasticity Imaging * Yu-Jiao Li(李玉娇), Wei-Jun Huang(黄伟骏), Feng-Chao Ma(马风超), Rui Wang(王睿), Ming-Zhu Lu(陆明珠)**, Ming-Xi Wan(万明习)** Affiliations Key Laboratory of Biomedical Information Engineering of Ministry of Education, Department of Biomedical Engineering, School of Life Science and Technology, Xi'an Jiaotong University, Xi'an 710049 Received 25 June 2016 ^*Supported by the National Key Scientific Instrument and Equipment Development Projects of China under Grant No 81127901, and the National Natural Science Foundation of China under Grant Nos 61372017 and 30970828.
^**Corresponding author. Email: mzlu@mail.xjtu.edu.cn; mxwan@mail.xjtu.edu.cn Citation Text: Li Y J, Huang W J, Ma F C, Wang R and Lu M Z et al 2016 Chin. Phys. Lett. 33 114301 Abstract A modified Monte Carlo model of speckle tracking of shear wave propagation in scattering media is proposed. The established Monte Carlo model mainly concerns the variations of optical electric field and speckle. The two-dimensional intensity distribution and the time evolution of speckles in different probe locations are obtained. The fluctuation of speckle intensity tracks the acoustic-radiation-force shear wave propagation, and especially the reduction of speckle intensity implies attenuation of shear wave. Then, the shear wave velocity is estimated quantitatively on the basis of the time-to-peak algorithm and linear regression processing. The results reveal that a smaller sampling interval yields higher estimation precision and the shear wave velocity is estimated more efficiently by using speckle intensity difference than by using speckle contrast difference according to the estimation error. Hence, the shear wave velocity is estimated to be 2.25 m/s with relatively high accuracy for the estimation error reaches the minimum (0.071). DOI:10.1088/0256-307X/33/11/114301 PACS:43.25.+y, 42.25.Dd, 42.30.Ms, 42.90.+m © 2016 Chinese Physics Society Article Text Acousto-optic shear-wave elasticity imaging is a novel quantitative and real-time imaging modality, which combines optical contrast with ultrasound spatial resolution for soft tissue functional imaging.[1,2] This technique is of great significance in medicine and has proven to be a promising tool for the early prediction of breast tumors.[3,4] An analytic model of the ultrasonic modulation of multiply scattered coherent light in scattering media as a theoretical reference standard was provided by Wang.[5] Three possible mechanisms have been identified. The first mechanism is based on ultrasound-induced variations of the optical properties of the media. The second mechanism is based on variations of the optical phase variation in response to ultrasound-induced displacements of the scatterers. The third mechanism is based on variations of the optical phase in response to ultrasonic modulation of the index of refraction.[5,6] Due to the strong optical scattering and minimal acoustically induced modulation, the optical modulation signal is weak, thus the acoustic radiation force (ARF) and shear wave are used to increase the ultrasound-induced particle displacement and to strengthen the optical modulation signal.[3] Hence, we are ready to use the Monte Carlo model to perform theoretical simulation of shear-wave-modulated optical detection. Shear wave elasticity imaging (SWEI) was first proposed in 1998 by Sarvazyan et al.[7] Physical and mathematical bases of shear wave elasticity imaging and the propagation of shear wave induced by acoustic radiation force were presented.[7,8] Shear wave propagation induced by acoustic radiation force in a non-dissipative inhomogeneous tissue was investigated, and the shear wave equation of radiation force in inhomogeneous media was solved numerically with a finite-difference time-domain method.[9] A Monte Carlo model which was mainly related to optical phase increments due to the acoustic-radiation-force shear-wave-induced displacements of light scatterers was established by Lu et al.[10] As optical phase increments are just part of the optical electric fields consisting of amplitude and phase, the model of Lu et al. merely tracks the shear wave propagation from the view of phase. Therefore, it is of significance to synthesize the amplitude and phase of optical electric fields which decide the distributions of optical signal intensities in a speckle pattern. A fast Monte Carlo model for simulating ultrasound-modulated light with a graphics processing unit (GPU) which can significantly accelerate calculations was proposed by Leung et al.[11] Nevertheless, they merely focus on the effect of ultrasonic modulation on speckle pattern and do not combine the acoustic-radiation-force-induced shear wave with the light model. Hence we would like to propose a speckle Monte Carlo model which not only synthesizes the amplitude and phase of optical electric field but also introduces the shear wave induced by the acoustic radiation force. Experimentally, a transient optoelastography technique to detect shear wave propagation was proposed, and the detection and discrimination of optical and mechanical properties in phantoms based on the spatial correlation between successive optical speckle patterns were demonstrated.[1,12] Li et al. established a specific UOT system and utilized a parameter of speckle contrast to observe the shear wave propagation in phantoms.[3] In this Letter, we introduce an optical imaging system in an established Monte Carlo model[10] with respect to acoustic-radiation-force shear-wave-induced optical phase increments and put forward a modified Monte Carlo model to investigate the dynamics of the photon electric field and the distribution of optical signal intensity in a speckle pattern. In this study, the amplitude and phase increments, which are influenced by acoustic-radiation-force shear-wave-induced displacements of light scatterers,[13] of an optical electric field were synthesized. The coordinate systems of our model are presented in Fig. 1. As the optical phase increments with a Gaussian distribution play an important role in the distributions of the autocorrelation function of photon electric fields, the optical phase increments induced by the variation of the displacements of the scatterers can be defined as[10] $$\begin{align} \Delta\phi_{j}(t)= n_{0}k_{0}\hat{\boldsymbol k}_{j}({\boldsymbol s}_{j}(t)-{\boldsymbol s}_{j+1}(t)),~~ \tag {1} \end{align} $$ where $n_{0}$ is the refractive index of the tissue, $k_{0}$ is the optical wave number in vacuum, $\hat{\boldsymbol k}_{j}$ is the unit optical wave vector after the $j$th scattering event, ${\boldsymbol s}_{j}(t)-{\boldsymbol s}_{j+1}(t)$ is the variation of displacement vector at time $t$ between two successive scattering events, and ${\boldsymbol s}_{j}(t)$ is calculated in accordance with the previously described methods.[7,10] The electric field of the $m$th transmitted photon is determined as follows:[11] $$\begin{alignat}{1} E_{s_{m}}(t;x_{m},y_{m})=W_{m}\exp\{-i[\sum_{j=1}^{N}\Delta\phi_{j}(t)]\},~~ \tag {2} \end{alignat} $$ where $s_{m}$ is the path length of the $m$th photon, $(x_{m},y_{m})$ are the coordinates of the exiting $m$th photon, $W_{m}$ is the weight of the $m$th photon and represents the probability of the $m$th photon to reach the detector after multiple scattering and absorption events, and $N$ is the total number of scattering events of the $m$th photon along its whole photon path. During photon propagation, the absorption of the photon weight alters the propagation property of the photon and the distribution of the electric field. For the sake of producing a speckle pattern, the electric field is spatially discretized by adding the contributions of each photon exiting within a certain region of the transmission plane[11] $$\begin{align} &E_{\rm s}(t;x_{p},y_{q})\\ =\,&\sum_{\scriptsize{\begin{matrix} x_{p}-\Delta x/2 < x_{m}\leq x_{p}+\Delta x/2\\ y_{p}-\Delta y/2 < y_{m}\leq y_{p}+\Delta y/2\end{matrix}}} \!\!\!E_{s_{m}}(t;x_{m},y_{m}),~~ \tag {3} \end{align} $$ where $x_{p}=p\Delta x$ and $y_{q}=q\Delta y$ with $p$ and $q$ as indices of the $x$ and $y$ coordinates of a region, and $\Delta x$ and $\Delta y$ are the widths of each discrete unit. A speckle pattern is formed by modeling an optical imaging system consisting of a lens with a radius of 0.2 cm and a pupil with a rectangular shape. The rectangular pupil function is expressed as follows: $$ P(x_{p},y_{q})=\begin{cases} \!\! 1, &-\frac{d_{x}}{2}\leq x_{p} \leq \frac{d_{x}}{2},~-\frac{d_{y}}{2} \leq y_{p} \leq \frac{d_{y}}{2},\\\!\! 0, &{\rm else}, \end{cases}~~ \tag {4} $$ where $d_{x}$ and $d_{y}$ are the length and width of the rectangular pupil, respectively. The speckle pattern intensity is calculated by using the following equation[11] $$\begin{alignat}{1} I(t;X_{u},Y_{v})=|{\rm FT}\{P(x_{p},y_{q})E_{\rm s}(t;x_{p},y_{q})\}|^{2},~~ \tag {5} \end{alignat} $$ where $I(t;X_{u},Y_{v})$ is the time-varying speckle pattern intensity with $u$ and $v$ being the indices of the $x$ and $y$ coordinates, respectively, $P(x_{p},y_{q})$ is the pupil function, and ${\rm FT}$ is the two-dimensional Fourier transform. The contrast parameter $C$ of each acquired speckle pattern in the presence of the ultrasound beam is calculated by[14] $$\begin{align} C=\frac{\sigma}{\langle I \rangle},~~ \tag {6} \end{align} $$ where $\sigma$ and $\langle I \rangle$ stand for standard deviation and mean of the speckle intensity, respectively.

Fig. 1. Schematic geometry of the coordinate systems. Here $(X_{1}, Y_{1}, Z_{1})$ and $(X_{2},Y_{2},Z_{2})$ are the acoustic and optic coordinate systems, respectively. The detection (imaging) plane is the transmission plane with blue gridding.

Fig. 2. Flowchart of the Monte Carlo model of a speckle from shear waves induced by acoustic radiation force.

The process of the speckle Monte Carlo model is constructed and described in the following. The photon number is set to be 1000000, and the simulation is conducted according to the flowchart of Fig. 2. In the simulation, the sample tissue is an infinitely wide slab in the $X_{2}$–$Y_{2}$ plane, with a thickness of tissue $L=1$ cm along the $Z_{2}$ axis (Fig. 1). The main parameters used for the simulation are listed in Table 1.[10] The focused ultrasound beam with a 5 MHz frequency and 6.10 cm focal length is positioned parallel to the $X_{2}$ axis. A light source with 532 nm wavelength and an initial position at (0, 0, 0), is moved along the $Y_{2}$ axis. The duration of the rectangular modulating pulse is 500 μs, which acts on the media. The position of the detected plane is $Z_{2}=1$ cm in the optical coordinate system. In our Monte Carlo simulation model, the effect of the Brownian motion is ignored and the scatterers are regarded as static in the absence of an ultrasound field. We have carried out multiple experiments and a set of typical results are revealed in this study.

Fig. 3. Speckles of the two-dimensional intensity distribution of different detection times at a probe location of 0.0 cm, (a) $t=1.8$ and (b) $t=1.9$ ms. The direction of the shear wave propagation is marked by a red arrow, and the fluctuation of the two-dimensional speckle intensity is marked by a red ellipse. The normalized intensity distribution along the $Y$ direction when $X=0$ mm is indicated by the solid line on the bottom of the speckle in the $X$–$Y$ plane, as shown in (c) and (d).

Figure 3 shows the effect of shear wave propagation on speckles. The variation of the brightest zone (marked by the ellipse) in the speckles of two-dimensional intensity distribution follows the propagation of the shear wave very well. The shear wave propagates from the left of the speckle to the right, as indicated by the red arrow. Table 2 summarizes the variation of the normalized peak intensity and full width at half maximum (FWHM) of one-dimensional intensity distribution curves of speckles along the $Y$ direction when $X=0$ mm. The normalized peak intensity initially decreases, then increases and reaches the minimum value at 1.8 ms. By contrast, the FWHM initially increases, then decreases over the detection time and reaches the maximum value at 1.8 ms. This phenomenon is due to the acoustic-radiation-force shear-wave-induced increase and enhancement of scattering events. The increase of scattering events broadens the distribution range of photons in the transmission plane, which results in an increase in FWHM. Simultaneously, stronger scatterings cause the longer photon path and induce (i) a higher probability of absorption by tissue, (ii) a greater loss of energy, and (iii) a lower probability to reach the transmission plane. Therefore, the peak of one-dimensional normalized intensity curves of speckles at $X=0$ mm decreases.

**Table 1.** The main parameters used for the simulation.[10]
Parameters	Value
Optical absorption coefficient $\mu_{\rm a}$	0.1 cm$^{-1}$
Optical scattering coefficient $\mu_{\rm s}$	10 cm$^{-1}$
Refractive index $n_{0}$	$1.33$
Tissue density $\rho $	10$^{3}$ kg/m$^{3}$
Acoustic wave velocity $c$	1540 m/s
Shear wave velocity $c_{t}$	2.1 m/s
Transducer aperture $a$	2.86 cm
Initial acoustical power $I_{0}$	10 W/cm$^{2}$
Anisotropy coefficient $g$	0
Viscosity coefficient $v$	0 Pa$\cdot$s
Elasto-optical coefficient $\eta $	0.321
Acoustic absorption coefficient $\alpha$	0.5 dB/(cm$\cdot$MHz)

**Table 2.** Peak and FWHM of normalized speckle intensity when $X=0$ mm at different detection times.
Detection time (ms)	1.0	1.8	2.6
Normalized peak intensity	0.767	0.031	0.951
FWHM	27	60	38

Fig. 4. Estimation of the shear wave velocity at a sampling interval of 0.2 ms: (a) the time evolution of speckle contrast difference at different probe locations, (b) the time evolution of speckle intensity difference (logarithmic form) at different probe locations, (c) linear regression processing of (a), and (d) linear regression processing of (b). Here $R$ is the correlation coefficient, and $y$ is the linear regression equation.

To estimate the shear wave velocity which is set to 2.1 m/s in theory, we calculate the speckle contrast difference $\Delta C$ and the speckle intensity difference in logarithmic form $\log(\Delta I)$ in a block with $60\times 60$ pixels. The center of the block is set according to the probe location. Figures 4 and 5 illustrate the time evolution of the speckle contrast difference and the speckle intensity difference (logarithmic form) at each probe location and the corresponding linear regression processing at sampling intervals of 0.2 ms and 0.5 ms, respectively. The peaks of the speckle contrast difference curves at different probe locations can track the peaks of the shear wave front, and a linear relationship is presented between them; a similar relationship exists for speckle intensity difference curves. The speckle contrast difference curve is sharper at a small sampling interval. At a sampling interval of 0.2 ms (Figs. 4(a) and 4(b)), the shear wave velocities are estimated to be 2.32 m/s (Fig. 4(c)) and 2.25 m/s (Fig. 4(d)). The estimation errors are 0.105 and 0.071, respectively. At a sampling interval of 0.5 ms, the shear wave velocities are estimated to be 2.42 m/s (Fig. 5(c)) and 2.30 m/s (Fig. 5(d)). The estimation errors are 0.152 and 0.095, respectively. At a sampling interval of 0.5 ms, the peak of the speckle intensity difference attenuates at 25.2% when the shear wave propagates to a probe location of 2 cm (the logarithmic peak intensities are 7.93 and 5.93 at 0 cm and 2 cm, respectively). These results demonstrate that a smaller sampling interval results in a higher estimation precision, and the shear wave velocity is more efficiently estimated by using speckle intensity difference than by using speckle contrast difference according to the estimation error. Consequently, the shear wave velocity is estimated to be 2.25 m/s with relatively high accuracy for the estimation error reaches the minimum (0.071) under the above described conditions. These results are consistent with those of tracking the shear-wave displacement by using optical phase increments[10] and speckle contrast difference.[14] Our Monte Carlo simulation can be employed to estimate the elastic characteristics of tissues for which there is a quantitative relationship between Young's modulus and the shear wave velocity, $E=3\rho c_{t}^{2}$.[3]

Fig. 5. Estimation of the shear wave velocity at a sampling interval of 0.5 ms: (a) the time evolution of speckle contrast difference at different probe locations, (b) the time evolution of speckle intensity difference (logarithmic form) at different probe locations, (c) linear regression processing of (a), and (d) linear regression processing of (b). Here $R$ is the correlation coefficient, and $y$ is the linear regression equation.

In conclusion, we have studied the effects of shear waves induced by acoustic radiation force on speckles based on an optical phase increment Monte Carlo model. In our simulation, the speckles of two-dimensional intensity distribution at different detection times are presented. The variation of the brightest zone in the speckles follows shear wave propagation well. The reduction of the peak on the speckle intensity difference curve implies attenuation during the shear wave propagation. Shear wave velocities are estimated quantitatively by extracting the time evolution of the speckle contrast difference and speckle intensity difference (logarithmic form). The results demonstrate that a smaller sampling interval yields higher estimation precision and the shear wave velocity is estimated more efficiently by using speckle intensity difference than by using speckle contrast difference according to the estimation error. Consequently, the shear wave velocity is estimated to be 2.25 m/s with relatively high accuracy for the estimation error reaches the minimum (0.071). The estimated shear wave velocity can be used to calculate Young's modulus quantitatively and to serve the estimation of tissue elastic characteristics. The parameter Monte Carlo speckle simulation provides a theoretical foundation for acousto-optic shear wave elasticity imaging and acts as a guidance for further experiments. References Detection and discrimination of optical absorption and shear stiffness at depth in tissue-mimicking phantoms by transient optoelastographyUltrasound-mediated optical tomography: a review of current methodsEffects of acoustic radiation force and shear waves for absorption and stiffness sensing in ultrasound modulated optical tomographyPharmacokinetic Monitoring of Indocyanine Green for Tumor Detection Using Photoacoustic ImagingMechanisms of Ultrasonic Modulation of Multiply Scattered Coherent Light: An Analytic ModelUltrasound-Mediated Biophotonic Imaging: A Review of Acousto-Optical Tomography and Photo-Acoustic TomographyShear wave elasticity imaging: a new ultrasonic technology of medical diagnosticsAcoustic radiation force and streaming induced by focused nonlinear ultrasound in a dissipative mediumPropagation of Shear Waves Generated by Acoustic Radiation Force in Nondissipative Inhomogeneous MediaMonte Carlo Simulation of Scattered Light with Shear Waves Generated by Acoustic Radiation Force for Acousto-Optic ImagingFast Monte Carlo simulations of ultrasound-modulated light using a graphics processing unitTransient optoelastography in optically diffusive mediaShear Wave Elasticity Imaging Based on Acoustic Radiation Force and Optical Detection

[1]	Daoudi K, Boccara A C and Bossy E 2009 Appl. Phys. Lett. 94 154103
[2]	Elson D S, Li R, Dunsby C, Eckersley R and Tang M X 2011 Interface Focus 1 632
[3]	Li R, Elson D S, Dunsby C, Eckersley R and Tang M X 2011 Opt. Express 19 7299
[4]	Yang S H, Yin G Z and Xing D 2010 Chin. Phys. Lett. 27 094302
[5]	Wang L V 2001 Phys. Rev. Lett. 87 043903
[6]	Wang L V 2004 Dis. Markers 19 123
[7]	Sarvazyan A P, Rudenko O V, Swanson S D, Fowlkes J B and Emelianov S Y 1998 Ultrasound Med. Biol. 24 1419
[8]	Rudenko O V, Sarvazyan A P and Emelianov S Y 1996 J. Acoust. Soc. Am. 99 2791
[9]	Lu M Z, Liu X J, Shi Y, Kang Y N, Guan Y B and Wan M X 2012 Chin. Phys. Lett. 29 014301
[10]	Lu M Z, Wu Y P, Shi Y, Guan Y B, Guo X L and Wan M X 2012 Chin. Phys. Lett. 29 124302
[11]	Leung T S and Powell S 2010 J. Biomed. Opt. 15 055007
[12]	Bossy E, Funke A R, Daoudi K, Boccara A C, Tanter M and Fink M 2007 Appl. Phys. Lett. 90 174111
[13]	Ju Z P, Cao W F and Liu X W 2009 Acta Phys. Sin. 58 174 (in Chinese)
[14]	Cheng Y, Li R, Li S, Dunsby C, Eckersley R J, Elson D S and Tang M X 2012 Ultrasound Med. Biol. 38 1637