Transcranial dipole localization and decoding study based on ultrasonic phased array for acoustoelectric brain imaging

The nonlinear Westervelt equation is used to describe the acoustic field in numerical simulations [12, 13]. It incorporates the effects of nonlinearity, attenuation, reflection and scattering and can be written as:

$$\begin{align}{\nabla ^2}p - \frac{1}{{{c^2}}}\frac{{{\partial ^2}p}}{{\partial {t^2}}} + \frac{\delta }{{{c^4}}}\frac{{{\partial ^3}p}}{{\partial {t^3}}} + \frac{\beta }{{\rho {c^4}}}\frac{{{\partial ^2}{p^2}}}{{\partial {t^2}}} = 0\end{align} \tag{ 1 }$$

where $p$ is the acoustic pressure, $t$ is the time, $\delta = 2{c^3}\alpha /{\omega ^2}$ is the thermoviscous diffusivity, $\beta$ is the nonlinear coefficient, $\omega = 2\pi f$ is the angular frequency, and $f$ is the ultrasonic center frequency. The equation was calculated by the finite difference time domain (FDTD) method in 3D Cartesian coordinates with a spatial step of 0.1 mm and time step of 10 ns. The Mur first-order absorbing boundary condition was used at the boundary to avoid ultrasonic reflections [14]. Here the effect of shear waves is ignored since most of the incidence angles are less than 20° and the maximum incidence angle is less than 25° [15].

Transcranial acoustic field modulation is realized based on the time reversal method. Specifically, a virtual point source ${S_0}\left( t \right) = {P_0}{\text{sin}}\left( {\omega t} \right)$ is first set up at the intracranial target focal point, and the acoustic pressure of the ultrasound emitted by this point source propagating to the center position of each array of the transducer is calculated by numerical simulation. The signal ${P_m}\left( {t - {T_t}} \right)$ received by array $m$ is then inverted in ${T_t}$ in a time sequence to obtain the inverted signal ${P_m}\left( t \right)$. Next, the corresponding fitting parameters are set based on the least squares function fitting method to realize the extraction of the initial phase of the target fitted signal, and the initial phase of array element $m$ is obtained as $\Delta {t_m}$. Finally, under the same acoustic pressure ${P_0}$, each array of the transducer realizes transcranial phase-modulated focusing at the target focal point based on different initial phases $\Delta {t_m}$, so the excitation signal of each array is:

$$\begin{align}{S_{0m}}\left( t \right) = {P_0}{\text{sin}}\left( {\omega \left( {t + \Delta {t_m}} \right)} \right).\end{align} \tag{ 2 }$$

The most basic imaging principle of ABI is the AE effect. The AE effect is a phenomenon in which the current passing through an electrolyte under the coupling of acoustic and electric fields is subjected to the periodic modulation of the focused acoustic pressure field, causing the local conductivity of the medium to vary with the ultrasonic rhythmicity [16]. The phenomenon was first discovered in 1946 by Fox et al. Ultrasound can periodically change the electrical conductivity of salt water by pressure [17]. The next study found that ultrasound affects the dielectric conductivity mainly by changing the physical and chemical properties of the medium such as molar concentration, ion ionization equilibrium and ion mobility [18–20].

The AE effect indicates that the change in medium conductivity $\sigma$ due to FUS is [21, 22]:

$$\begin{align}\frac{1}{{\Delta \sigma }} = {K_I}\frac{1}{{{\sigma _0}}}\Delta P\end{align} \tag{ 3 }$$

where $\Delta \sigma$ is the change in medium conductivity, ${\sigma _0}$ is the intrinsic conductivity, $\Delta P$ is the change in acoustic pressure, and ${K_I}$ is the interaction constant, which is on the order of ${10^{ - 9}}{\text{P}}{{\text{a}}^{ - 1}}$ in a 0.9% NaCl solution.

The voltage ${V_i}$ obtained from the measurement of lead numbered $i$ in the EEG acquisition is generally expressed as:

$$\begin{align}{V_i} = \iiint {\frac{1}{\sigma }\left( {\widetilde {J_i^L} \cdot {J^I}} \right)}{\text{d}}x{\text{d}}y{\text{d}}z\end{align} \tag{ 4 }$$

where $\tilde J_i^L$ is the current density distribution of lead $i$ and ${J^I}$ is the current density distribution of neurons.

Based on equation (3), the local conductivity of the medium at the FUS modulation is:

$$\begin{align}\frac{1}{\sigma } = \frac{1}{{{\sigma _0}}} - {K_I}\frac{1}{{{\sigma _0}}}\Delta P.\end{align} \tag{ 5 }$$

Substituting (5) into (4) gives

$$\begin{align}{V_i} & = V_i^{{\text{LF}}} + V_i^{{\text{AE}}} = \iiint {\frac{1}{{{\sigma _0}}}\left( {\tilde J_i^L \cdot {J^I}} \right)}{\text{d}}x{\text{d}}y{\text{d}}z \notag\\ & \quad + \iiint {\left( { - {K_I}\frac{1}{{{\sigma _0}}}\Delta P} \right)}\left( {\tilde J_i^L \cdot {J^I}} \right){\text{d}}x{\text{d}}y{\text{d}}z\end{align} \tag{ 6 }$$

where $V_i^{{\text{LF}}}$ is the low-frequency EEG signal of ${V_i}$ and $V_i^{{\text{AE}}}$ is the high-frequency AE signal. The high-frequency AE signal acquired at the point $\left( {{x_1},{y_1}} \right)$ on the plane of FUS irradiation after high-pass filtering is:

$$\begin{align} V_i^{{\text{AE}}}\left( {{x_1},{y_1},t} \right)& = - {K_I}\!\!\iiint\!\! {\frac{1}{{{\sigma _0}}}\left( {\tilde J_i^L \cdot {J^I}} \right)}b\left( {x - {x_1},y - {y_1},z} \right)\notag\\ & \quad \times {P_0}a\left( {t - z/c} \right){\text{d}}x{\text{d}}y{\text{d}}z\end{align} \tag{ 7 }$$

where $b\left( {x,y,z} \right)$ is the ultrasonic beam pattern, ${P_0}$ is the acoustic pressure amplitude, $a\left( t \right)$ is the pulse waveform, and $c$ is the acoustic velocity.

Equation (7) can be simplified as [23]:

$$\begin{align}V_i^{\text{AE}}\left( {{x_1},{y_1},t} \right) \approx {K_a}\left( {\tilde J_i^L\left( {{x_1},{y_1}} \right) \cdot {J^I}\left( {{x_1},{y_1}} \right)} \right)A\left( t \right)\end{align} \tag{ 8 }$$

where ${K_a} = - {K_I}{P_0}\frac{1}{{{\sigma _0}}}b\left( {x,y} \right)$, $A\left( t \right) = \int a\left( {t - z/c} \right){\text{d}}z$. When ${A_0} = \max \left\{ {A\left( t \right)} \right\}$, equation (8) becomes:

$$\begin{align}V_i^{\text{AE}}\left( {{x_1},{y_1}} \right) \approx {K_a}{A_0}\left( {\tilde J_i^L\left( {{x_1},{y_1}} \right) \cdot {J^I}\left( {{x_1},{y_1}} \right)} \right).\end{align} \tag{ 9 }$$

The minimum paradigm estimate of the current density ${J^I}\left( {{x_1},{y_1}} \right)$ at point $\left( {{x_1},{y_1}} \right)$ is:

$$\begin{align}{J^I}\left( {{x_1},{y_1}} \right) = \frac{1}{{{K_a}{A_0}}}{\left( {T\left( {{x_1},{y_1}} \right)} \right)^ + }V_i^{{\text{AE}}}\left( {{x_1},{y_1}} \right)\end{align} \tag{ 10 }$$

where Moore–Penrose pseudoinverse $T{\left( {{x_1},{y_1}} \right)^ + } = {\left( {{T^{\text{T}}}T} \right)^{ - 1}}{T^{\text{T}}}$.

Therefore, the divergence of the neuronal current density distribution ${J^I}\left( {x,y} \right)$ is:

$$\begin{align}\triangledown \cdot {J^I}\left( {x,y} \right) = - {I_{{\text{source}}}}\left( {x,y} \right).\end{align} \tag{ 11 }$$

From the current density distribution of the dipole, it can be seen that the current density at the dipole is the largest, and as the current spreads outward, the current density in the region farther away from the dipole is smaller. From the AE effect, it can be seen that when FUS irradiates the target area at a certain frequency (center frequency/pulse repetition frequency), the AE signals generated in the region of higher current density are stronger, while the AE signals in the region of lower current density are relatively weaker. Therefore, the ABI can localize the dipole position based on the AE signal value at each point of the ultrasound scan, and the larger the AE signal amplitude, the closer the localization position is to the actual position of the dipole.

Figure 1 shows the schematic diagram of the numerical simulation model. The simulation model mainly consists of a 128-array ultrasound transducer, an acrylic simulated skull and water. Among them, the diameter of the transducer opening is 100 mm, the radius of curvature is 60 mm, and the diameter of the array is 4 mm; the thickness of the acrylic plate is 8 mm, the radius of curvature is 100 mm, and it is 15 mm from the surface of the transducer. The numerical simulation area is 100 mm × 100 mm × 100 mm. The numerical simulation parameters are shown in table 1.

Zoom In Zoom Out Reset image size

Acoustic field tests are conducted based on the 128-array ultrasound phased array experimental platform (passed assembly error test) shown in figure 2, in which figure 2(a) shows the overall system diagram, which mainly includes the upper control module (the red circles represent the transmitting arrays and the blue circles represent the receiving arrays), high precision oscilloscope, power amplifier module, data acquisition module, 3D adjustment frame, positioning module, and acoustic field scanning module; figure 2(b) shows a partial enlargement of the positioning module, which utilizes a laser collimator to perform the geometrical focal point and simulated cranial localization; and figure 2(c) shows a partial enlargement of the acoustic field scanning module, which mainly includes the hydrophones, bone plate fixation device, acrylic plate, and transducer. In this case, the thickness of the acrylic plate is 8 mm, the radius of curvature is 100 mm, and the distance of the acrylic plate from the transducer is 15 mm, and its center is aligned with the center of the transducer along the acoustic axis. The acquisition range is the XZ section enclosed by the four points (−4, 0, 55), (4, 0, 55), (4, 0, 65) and (−4, 0, 65), with a scanning length of 8 mm in the X direction and 10 mm in the Z direction, and a scanning step of 0.1 mm.

Zoom In Zoom Out Reset image size

Under the condition that the initial excitation voltage of each array is 5 V, the transcranial self-focusing acoustic field distribution without any modulation of the arrays is firstly investigated, and then the modulated transcranial focusing acoustic field distribution is investigated based on the acoustic field modulation method described in 2.1.2. The acoustic field experiments were all done in a water tank with dimensions of 260 mm × 260 mm × 170 mm.

Figure 3 shows the schematic diagram of the ultrasonic phased array-based ABI experimental platform and ultrasonic irradiation range. Among them, figure 3(a) adds a multi-channel signal source (QDAC-II, QDevil, Denmark) and an EEG signal acquisition module (Neuroscan Curry 9, USA) to figure 2; figure 3(b) shows a partial enlargement of the signal source and ultrasonic irradiation region in figure 3(a); figure 3(c) shows a partial enlargement of the electrode and acrylic plate irradiation area in figure 3(b); figure 3(d) shows the schematic diagram of the imaging chamber and ultrasound incidence direction and electrode insertion direction, in which the imaging chamber was designed by SolidWorks and fabricated by 3D printing (40 mm × 40 mm × 30 mm); figure 3(e) shows the ultrasound phased array grid scanning area and electrode placement position, the gray circle is the ultrasound modulated irradiation site, and the purple and yellow areas are the scanning ranges during single and double dipole experiments, respectively; figure 3(f) shows the phantom physical image.

Zoom In Zoom Out Reset image size

The electrodes used in the study are platinum wire electrodes with a length of 15 mm and a diameter of 0.5 mm, which are soldered and connected to the internal signal wires of the BNC, and each pair of dipoles generates sinusoidal signals with the same frequency, equal amplitude, and opposite phase to simulate neuronal discharge process based on a multi-channel signal generator. The electrode tip is located in the y = 0 cross-section and the ultrasound is incident along the z-axis through the open side of the imaging chamber (figure 3(d)). The dipole, acquisition electrode, reference electrode and ground electrode setting positions for the experiment are shown in table 2, where Dipole1+ is set at (57, 0) and Dipole2+ is set at (60, 0).

Each electrode was inserted into the y = 0 section in the phantom along the y-axis direction according to the set position, where the phantom is composed of acrylamide (AM), ammonium persulfate (APS), N, N'-methylenebisacrylamide (MBA), and N, N, N', N'-tetramethylethylenediamine (TEMED), egg white solution, and physiological saline [10, 28].

For the single-source experiments in Dipole1+, the transducer was modulated to perform a $2\,{\text{mm}} \times 2\,{\text{mm}}$ raster scan of the y = 0 cross-section in 0.2 mm steps using the array modulation method in 2.1.2 at an initial excitation voltage of 5 V for each array (purple area in figure 3(e)). For double-source experiments in Dipole1+ and Dipole2+, $4.2\,{\text{mm}} \times 1.2\,{\text{mm}}$ raster scans were performed in 0.2 mm steps (yellow area of figure 3(e)). Pure water single-source phantom experiments, transcranial single-source phantom experiments, and transcranial double-source phantom experiments were performed. The voltage signals generated by the dipole electrodes were acquired, amplified and filtered by the SynAmps2 system (Neuroscan Curry 9, USA) with a sampling rate of 20 kHz, an acquisition time of 20 s for each point, and a bandpass filtering range of 0–3500 Hz.

The voltage signals acquired by SynAmps2 during ultrasonic irradiation of each target point were used as raw data, and the raw data were pre-processed by removing 50 Hz industrial frequency interference and downsampling from 20 kHz to 4 kHz. The pre-processed data were band-pass filtered at ultrasound pulse repetition frequency (PRF) ± 50 Hz to preserve the high-frequency signal at the 500 Hz. The 20 s envelope signal data was selected and divided into 20 segments of data at 1 s intervals, the data was wave-aligned, superimposed and averaged to obtain the 1 s envelope signal, which is the decoded signal of the ultrasonic irradiation point. The peak average value of the decoded signal was taken as the AE signal amplitude at that ultrasound irradiation point, and the AE signal values obtained at each point during the ultrasound raster scan were used to plot the dipole localization map.

Figure 4 shows a schematic of array calibration, transcranial geometric focusing, and offset geometric focusing array delays. Figure 4(a) shows the geometrically focused acoustic pressure of each array under different excitation conditions. Figure 4(b) shows a schematic diagram of the transcranial geometric focusing array delay, where the black line is the delay curve of each array after modulation by the method of 2.1.2, and the red line is the delay curve of each array element of the self-focusing method without any modulation. Figure 4(c) shows a schematic diagram of the delay of the transcranially modulated offset focusing arrays, where the black, blue, green and orange lines are the delay curves of each array at the four focal points of (0,0,60), (0,0,59), (0,0,58), and (0,0,57), respectively, after modulation by the method of 2.1.2.

Zoom In Zoom Out Reset image size

As can be seen from figure 4(a), the average values of acoustic pressure at the geometric focus of each array are 7.8 kPa, 17.95 kPa, and 22.52 kPa for excitation voltages of 2 V, 5 V, and 8 V, respectively, with a good output consistency. As can be seen from figure 4(b), the initial delay of each array after modulation varies, with the earliest excitation of array No. 111, which has an initial delay of 0 ns, and the latest excitation of array No. 52, which has an initial delay of 1590 ns; the initial delays of each array without modulation are all 0 ns. As can be seen from figure 4(c), the initial delays of each array after modulation are all different, when the preset focuses are (0,0,60), (0,0,59), (0,0,58) and (0,0,57), in which the No. 52 array are all the latest excitation, and the initial delays of the No. 52 array at each focus are 1590 ns, 1400 ns, 1210 ns and 990 ns, respectively; the No. 111 array are all the earliest excitation, and its initial delay is 0 ns.

Figure 5 shows the results of the transcranial modulation numerical simulation. Figure 5(a) shows the set focal point (0, 0, 57). Figure 5(b) shows the set focal point (0, 0, 58); figure 5(c) shows the set focal point (0, 0, 59). Figure 5(d) shows the set focal point (0, 0, 60). As can be seen from figure 5, the acoustic pressure field after modulation under numerical simulation conditions can be focused at the set position, and the focal acoustic pressure is 0.49 MPa.

Zoom In Zoom Out Reset image size

The results of the comparison of transcranial self-focusing and modulated focusing at the geometric focus are shown in figure 6, where figure 6(a) shows the self-focusing acoustic pressure field, figure 6(b) shows the modulated focusing acoustic pressure field, figure 6(c) shows the self-focusing −6 dB focal area, figure 6(d) shows the modulated focusing −6 dB focal area, figure 6(e) shows the contrast of the focusing field contours, and figure 6(f) shows the contrast of the acoustic axial sound pressure curves. From figures 6(a) and (b), the maximum acoustic pressure is 0.18 MPa for self-focusing and 0.46 MPa for modulated focusing. The modulated acoustic pressure was increased by a factor of 1.56. From figures 6(c)–(e), it can be seen that the area of the self-focusing focal domain is 4.09 mm² and there is a large parafocal flap; after modulation of focusing, the focal domain is 2.64 mm in the long axis and 0.83 mm in the short axis, with an area of 1.65 mm² and basically no parafocal flap. The focal domain area was reduced by 59.66% after modulation. From figure 5(f), it can be seen that the focus of self-focusing exists about 1 mm forward; the focus of modulated focusing is accurately focused at the preset (0, 0, 60), and the effect of modulated focusing is obviously better than self-focusing.

Zoom In Zoom Out Reset image size

Under the conditions of setting the focus as (0,0,60), (0,0,59), (0,0,58), and (0,0,57), the modulated focusing results are shown in figure 7, where figures 7(a)–(d) show the focused acoustic pressure fields when different focuses are set, and figures 7(e)–(h) show the corresponding attenuated −6 dB focusing fields, respectively. From figures 7(a)–(d), it can be seen that the maximum pressure of the acoustic field at the geometrical focal point for focusing and −Z direction shifted forward by 1 mm, 2 mm, and 3 mm are 0.46 MPa, 0.53 MPa, 0.55 MPa, and 0.50 MPa, respectively. From figures 7(e)–(h), the focal domain areas at each point are 1.65 mm², 1.62 mm², 1.25 mm² and 1.51 mm², respectively.

Zoom In Zoom Out Reset image size

The focal domain contour and the acoustic axis pressure curve of the modulated offset focus are shown in figure 8, where figure 8(a) shows the focal domain contour and figure 8(b) shows the acoustic axis pressure curve. As can be seen from figure 8, the regulated focus is basically focused at the preset focus, realizing a regulated offset focus of 3 mm.

Zoom In Zoom Out Reset image size

4.2.1.1. Localization

The voltage signals acquired under the conditions of setting the phased array excitation voltage of 5 V, the Dipole1 discharge frequency of 8 Hz and the set position of (57, 0) are shown in figure 9, in which figure 9(a) shows the voltage time-domain signals with/without ultrasonic irradiation, and blue line is the signal with ultrasonic irradiation, and orange line is the signal without ultrasonic irradiation; figure 9(b) shows the spectrograms of the low-frequency bands of the two parts; and figure 9(c) shows the high-frequency bands of the spectrograms at the PRF of the two parts. As can be seen in figure 9(a), the time-domain signals of the voltages collected with/without ultrasonic irradiation are basically the same, and there is a high-frequency spike 'burr' in the time-domain signals with ultrasonic irradiation (black arrows zoom). As can be seen from figure 9(b), the spectrograms of the low-frequency bands with and without ultrasonic irradiation are basically the same for both conditions, and both are dominated by 8 Hz. As can be seen in figure 9(c), the acquired signal without ultrasonic irradiation has no specific frequency in the high-frequency band, and the acquired signal with ultrasonic irradiation contains the PRF of ultrasound (500 Hz).

Zoom In Zoom Out Reset image size

Figure 10 shows the results of dipole localization in the phantom via pure water and simulated skull under the above conditions, respectively, where figures 10(a) and (b) shows the results of dipole localization under pure water; figures 10(c) and (d) shows the results of dipole localization under simulated skull. As can be seen from figures 10(a) and (b), the maximum value of the AE signal in the pure water experiment is 725.31 μV, which locates the Dipole1+ position at (0, 56.8), and the z-direction error is 0.2 mm compared with the set position. As can be seen from figures 10(c) and (d), the maximum value of the AE signal in the simulated skull experiment is 134.22 μV, which localizes the Dipole1+ position at (0, 56.6), and the z-direction error is 0.4 mm compared to the set position. As expressed in equation (7) in section 2.2.2, the AE signal is directly related to the focal acoustic pressure. The AE signal was reduced by 81.49% due to the presence of the simulated skull.

Zoom In Zoom Out Reset image size

4.2.1.2. Decoding

The acquired single dipole source voltage signal is processed by the decoding algorithm of 3.2.3 and the decoded signal is obtained as shown in figure 11, in which figure 11(a) shows the comparison between the decoded signal and the source signal under pure water, with the blue curve being the decoded signal and the orange curve being the source signal; figure 11(b) is the spectrogram of the experimental decoded signal under pure water experiment; figure 11(c) is the comparison of the decoded signal and the source signal of the under simulated skull; and figure 11(d) is the spectrogram of the experimental decoded signal of the under simulated skull experiment. Calculations from the data in figures 11(a) and (b) show that the Pearson correlation coefficients of the decoded signal and source signal of Dipole1+ under pure water and simulated skull experiments are 0.94 and 0.93, respectively.

Zoom In Zoom Out Reset image size

4.2.2.1. Localization

Under the conditions of setting the phased array excitation voltage of 5 V, Dipole1+/Dipole2+ positions of (0, 57)/(0, 60) and discharge frequencies of 8 Hz/13 Hz, the localization results of the simulated skull experiments are shown in figure 12, in which figure 12(a) is the localization imaging map of the AE signal; and figure 12(b) is the localization imaging attenuation thermogram. From figure 12(a), the maximum value of AE signal near Dipole1+ is 186.66 μV; the maximum value of AE signal near Dipole2+ is 236.18 μV. From figure 12(b), it can be seen that the position of localized Dipole1+ is at (0, 57.2), and the z-direction error in the simulated skull experiment is 0.2 mm compared with the set position; the position of localized Dipole2+ is at (−0.2, 60), and the x-direction error in the simulated skull experiment is 0.2 mm compared with the set position.

Zoom In Zoom Out Reset image size

4.2.2.2. Decoding

The acquired double dipole source voltage signals are processed by the 3.2.3 decoding algorithm and the decoded signals are shown in figure 13, in which figures 13(a) and (c) are the comparison maps of the decoded signals and source signals of Dipole1+ and Dipole2+ by the simulated skull experiments, the blue curve is the decoded signals, and the orange curve is the source signals; and figures 13(b) and (d) are the spectral maps of the decoded signals of Dipole1+ and Dipole2+ by the simulated skull experiments, respectively. The Pearson correlation coefficients between the Dipole1+/Dipole2+ decoded signal and the source signal in the simulated skull experiment can be calculated from the data in figures 13(a)–(d) to be 0.90 and 0.89, respectively.

Zoom In Zoom Out Reset image size