We developed an in vitro experimental setup for measuring ultrasound stimulation-induced activity in mouse brain slices. We assumed that our measurement system could record extracellular potentials in brain tissue using an MEA substrate (Figure 1A). For extracellular recording, we conducted ultrasound irradiation using a transducer (V301-SU, Olympus, Waltham, MA, United States) with an MEA probe (MED-P530A, Alpha MED Scientific; glass substrate with a thickness of 0.7 mm and 64 microelectrodes in an area 2.1 × 2.1 mm²; Figure 1Ab). Each of the 64 recording sites (or channels; Chs) in the MEA covered an area of 50 × 50 μm², and the interelectrode interval between the centers of adjacent sites was 300 μm. In the experiment, ultrasound waves were applied to a brain slice located on the glass substrate of an MEA probe. The waves were delivered from the transducer to the brain tissue through the bottom of the substrate (Figure 1Aa) via a custom-designed waveguide (collimator) (Figure 1Ba). Once we determined the size of the brain slices and the parameters of the measurement equipment, which included an ultrasound transducer, MEA probe, and perfusion system, the waveguide was the left as the only adjustable component in the measuring system. Therefore, in this study, we mainly focused on waveguide design by numerically and experimentally evaluating the waveguide properties.

When choosing the parameters for the ultrasound stimulation, we assumed that the focal length would need to be such that the stimulation would pass through glass MEA substrate (0.7 mm) and travel further to penetrate the brain slice (0.4 mm): i.e., 0.7–1.1 mm from the bottom edge of the MEA probe to the brain slice. Otherwise, the standing waves in the chamber containing the MEA probe would be dominant, and the power distribution in the chamber during ultrasound stimulation would be too complex for our analysis.

We calculated the average acoustic intensity (power) of the ultrasound pulse according to the acoustic pressure distribution of the ultrasound wave traveling through the waveguide (Figure 1A). We defined the spatial resolution as the length corresponding to more than half of the maximum ultrasound intensity. To determine the localization of the ultrasound intensity on the MEA substrate amongst the 64 electrodes, we set the spatial resolution (centered at the peak position) to be less than 2.1 mm. This was the size of the area covered by the multielectrode array and the central area of a typical brain slice, i.e., 3.5 × 4.0 mm² (Figure 1Ab).

To design the waveguide used to focus the ultrasound irradiation within a brain slice, we numerically calculated the spatial distribution of ultrasound intensities within the slice. We calculated the focal length from the transducer to the focal point using the waveguide located on top of the transducer. Because of the axial (cylindrical) symmetry around the center axis, we first created a two-dimensional (2D) simulation of the ultrasound propagation using the finite-difference time-domain (FDTD) method (Toda and Tateno, 2020). A three-dimensional (3D) simulation was also performed (as explained in the next subsection) to examine the 3D distribution of the ultrasound intensity.

Here, we explored the desired acoustic pressure field by searching for appropriate structural parameters regarding the length (l) and diameter (d) of the opening end of the waveguide (Figure 1A). Because of a size restriction in our recording system, the height of the waveguide was set at 35 mm (Figure 1Aa). Thus, the length l must be equal to or smaller than the height (i.e., l ≤ 35 mm).

In the 2D FDTD method, the governing equations of medium particle motion are described as:

where ρ_w is water density (ρ_w = 1.0 × 10³ kg/m³), v_x and v_y are the particle velocities of the medium in the x and y directions, respectively, p is pressure, and t is time. The continuous equation concerning pressure p is described as:

where κ is the bulk modulus of the medium (κ = 2.19 × 10⁹ kg/m²s).

The edge of the computational area was conditioned as a free boundary, and the perfectly matched layer (PML) boundary (Chew and Liu, 1996) was used as a boundary condition (Eq. 3a). In addition, the boundary between the waveguide and the water was set as the sound-absorbing boundary. It was mounted using surface normal impedance Z_n, described as:

where c_w is the sound velocity in water (c_w = 1.48 × 10³ m/s) and a_n is the sound absorption coefficient at normal incidence. The minus sign in Eq. 3a indicates the reversal of the particle velocity direction because of reflection. The spatial-discrete width Δh was set to 1/20 of the wavelength (input wave number = 5) of the upper limit of the frequency to be analyzed; Δh = 5.0 × 10^–5 m. The time-discrete width was set to satisfy the Courant–Friedrichs–Lewy (CFL) condition, which was necessary for calculation stability. The coolant C, with the CFL condition as a dimensionless number, is described as:

where Δt is the time-discrete width (Δt = 5.0 × 10^–9 s) because the coolant number is usually smaller than 1/2 in 2D simulations.

The simulation field area is illustrated in Figure 2. The area was surrounded by a rectangle of 50 × 100 mm². The PML boundary was defined as the edges of the rectangle. The input ultrasound (at y = 0 in Figure 2A) was a sinusoidal wave (f = 500 kHz) that passed through the waveguide, which was filled with water. The transducer diameter was 25 mm in our experimental setup (Figures 1Ab,B). The waveguide sidewall was made of polypropylene (PP) with a sound absorption coefficient a_n of 0.876, calculated using the acoustic impedance between water and PP (Lochab et al., 2004). At the top of the waveguide, the bottom of the glass MEA substrate was attached to the surface of the waveguide, which effectively transmitted the power of the ultrasound waves (Figures 1Aa, 2Aa).

As an index of ultrasound intensity, we used the power integral intensity (spatial peak pulse average, I_sppa), which was calculated from acoustic pressure p(t) and described as:

where the time window T for the integral ranged from 0 to 125 ms (i.e., T = 125 ms). The maximum ultrasound intensity was analyzed using I_sppa (Eq. 5). The area (length in the 2D simulation) in which the intensity was more than the half of the maximum intensity was defined as the half-value area (length), and the longest width in the area (length) was defined as the half width (HW) (Figure 2Ab). When standing waves of pressure were observed, as was the case for some parameter sets, several local peaks appeared. In such cases, the half-value area (length) was obtained after calculating the envelope of the corresponding peak intensities. In vitro, when the amount of artificial cerebral spinal fluid (ACSF) solution in the MEA chamber changes, the height (h) of the fluid level changes, affecting the intensity distribution. Therefore, we varied the amount of ACSF solution such that the fluid level h increased from 1.0 to 10 mm in 1.0-mm steps. For numerical calculations, we used a parallel computer system (PRIMERGY CX400/CX2550, Fujitsu, Japan) at the Hokkaido University Computer Center. In addition, we used GNU Fortran 5.4.0 for the compiler and OpenMP for the parallel programming.

To confirm the results obtained from the 2D FDTD method under more practical conditions, we conducted a simulation with a 3D model using the finite element method (FEM). The simulation was performed using commercially available multiphysics software, COMSOL Multiphysics (ver. 5.5, COMSOL AB, Sweden). The 3D model consisted of a transducer, a waveguide, an MEA probe with a glass chamber, a brain slice, ACSF solution, and the air surrounding these components (Figure 3A). As with the 2D model, the amount of ACSF solution in the MEA chamber was varied in 10 levels (h) ranging from 1.0 to 10 mm in 1.0-mm steps.

The geometry of the model was discretized to the extra fine mesh (the maximum is 1.23 mm and the minimum is 0.123 mm) setting with a swept mesh for the infinite air domain (Figure 3B). The acoustic pressure at the surface of the transducer was set to 100 kPa (I_sppa = 0.34 W/cm²; for a definition of I_sppa, see Eq. 5). The rectangular brain slice (4.0 × 3.5 × 0.4 mm³) was located on the surface of the MEA substrate, which had an acoustic impedance of 1,558 kPa⋅s/m. With respect to the pressure acoustics, the frequency domain-interface and laminar flow-interface in the software were used in the numerical calculations.

The boundary conditions of the sidewalls in the chamber and waveguide were the set using the “sound hard boundary” option in the program. The “sound soft boundary” was applied for the outer boundaries of the air domain and the boundary between the saline solution and air. We performed 3D modeling with the FEM using the same parallel computer system (PRIMERGY CX400/CX2550, Fujitsu, Japan) at the Hokkaido University Computer Center. In the model, the boundary conditions of the surrounding air domain were set to approximate infinity so that they minimally affected the numerical solution.

To examine the irradiation properties of ultrasound waves traveling through our designed waveguide, we measured the acoustic pressure using a calibrated needle hydrophone (HY05N, Toray Engineering Co., Chuo, Japan) with a calibration certificate supplied by the manufacturer. The hydrophone was vertically lowered into the ACSF solution in the chamber of an MEA probe, and the measured data were recorded using an oscilloscope with a digital data acquisition system (SDS 1104X-E, SIGLENT Technologies, Shenzhen, China).

To characterize the waveguide properties, we used a continuous ultrasound waveform (CW) with a fixed window length (duration, 20 μs) to record the hydrophone output. We started recording once the signal had stabilized. Because of the central symmetricity, we measured pressure in one quadrant of the circular chamber filled with ACSF solution. To measure hydrophone output at different planar positions, the device was initially positioned at the approximate center of the MEA plane. Then, the hydrophone was repositioned in 0.5 mm increments to a maximum distance of 4.5 mm from the center on each axis. The distance between the ultrasound transducer and the hydrophone was set to approximately 1.0 mm. The distance was confirmed using sound velocity of the trigger delay time in distilled water or ACSF solution. The sound velocity was calculated as a function of temperature, and the temperature of water (or ACSF solution) was measured using a thermometer (RT-31S, ESPEC MIC CORP., Niwa, Japan). To place a brain slice on an MEA substrate, we used a nylon mesh and an anchor (see Section “Multielectrode-Array for Recording Activity in Brain Slices”). Because the nylon mesh was 1.1 mm mesh opening and 0.12 mm thick, there was little difference between our pressure measurements in the presence and absence of the mesh. Also, a U-shaped anchor was made from a 0.3-mm stainless wire and spacing by the edge distance of 8 mm. The anchor itself did not affect the pressure measurement on the area of multielectrodes (2.1 × 2.1 mm²).

For the CW, the original rectangular wave with a fundamental frequency of 500 kHz was generated by a function generator (WF1947, NF Co., Yokohama, Japan). In terms of the temporal window pattern, the CW stimuli had a duration of 100 or 200 ms. The acoustic pressure was varied from 110 to 410 kPa in 100 kPa steps: i.e., 110, 210, 310, and 410 kPa, and the corresponding power integral intensity (spatial peak pulse average, I_sppa in Eq. 5) was 0.41, 1.49, 3.24, and 5.68 W/cm², respectively. The input voltage signals from the function generator to the transducer were all amplified with a 50-dB radio-frequency amplifier (ZSA5064, RAD Co., Fuji, Japan).

As mentioned in the previous subsection, a 500-kHz ultrasound signal was delivered from a planar ultrasound transducer (V301-SU, Olympus) through the custom-designed waveguide attached to the transducer, which was positioned vertically below the MEA probe (at a 90-degree angle) (Figure 1A). The waveguide was optimized to activate a small region of a brain slice (HW ≤ 2.1 mm). This design helped minimize the generation of standing waves and created a favorable interface on the bottom of the MEA dish. Thus, it was possible to use this configuration to produce a controllable ultrasound intensity distribution. The custom-designed waveguide was then produced using a 3D printer (Dreamer, Flashforge, Zhejiang, China) with a polypropylene (PP) filament (1.75-mm diameter; Flashforge, Osaka, Japan). To align the properties of the transducer with those of the waveguide, the waveguide end was sealed with polyethylene, and the inside was filled with degassed water.

The acoustic frequency and intensity characteristics of an ultrasound waveform stimulus are essential elements of its core effects on brain activity. Waveform patterns can also influence the efficiency of neural activation. For our neural stimulation, we used a simple sinusoidal waveform: a CW signal with a fixed window length. For the CW, we generated an original rectangular wave with a fundamental frequency of 500 kHz using a function generator (WF1947). In terms of the temporal window pattern, the duration of the CW stimuli was 100 or 200 ms. The acoustic pressure varied from 110 to 410 kPa in 100 kPa steps: i.e., 110, 210, 310, and 410 kPa, such that the intensity (I_sppa in Eq. 5) corresponded to 0.41, 1.49, 3.24, and 5.68 W/cm², respectively. The signals from the function generator were amplified using a 50-dB radiofrequency amplifier (ZSA5064).

All animal experiments described below were approved by the Institutional Animal Care and Use Committee of Hokkaido University and carried out in accordance with the National Institutes of Health Guidelines for the Care and Use of Laboratory Animals. In the present study, we used five C57BL/6J mice (three male and two female mice, 4–8 weeks old, Japan SLC Inc., Hamamatsu, Japan). Each mouse was deeply anesthetized with isoflurane and decapitated. When sectioning brain slices that included the auditory cortex, the distance from the bregma along the rostral/caudal axis in each coronal section was defined in the following way. We referred to a digitized atlas (Franklin and Paxinos, 2008) for illustrations of coronal sections of the mouse brain. We chose the coronal section that best matched the illustration at 2.92 mm caudal to bregma (Figure 5 in Franklin and Paxinos, 2008). Coronal 400-μm-thick slices that included the auditory cortex were cut using a tissue slicer (Linear Slicer Pro7, D.S.K., Kyoto, Japan) in chilled ACSF. The ACSF contained (in mM) 119 NaCl, 2.5 KCl, 2.5 CaCl₂, 1.3 MgSO₄, 1.0 NaH₂PO₄, and 11.0 D-glucose (pH = 7.4). All slices were left in a submerged-type holding chamber at 28°C in a water bath for at least 2 h before recording.

To record ultrasound-driven local field potentials (LFPs) from mouse brain slices in vitro, we developed a custom in vitro ultrasound stimulation system that we combined with an MEA-based recording system (Figure 1A). All electrophysiological recordings in brain slices were performed with perfusion of ACSF solution with a stable flow rate of 1.26 ± 0.01 ml/min. During the recording, the flow rate was nearly constant, such that we expected the height of the ACSF solution in the MEA chamber to be unchanged. During the perfusion, mixed gas was supplied from the top of the recording chamber (APC-30, Asteck Co., Kasuya, Japan), which was maintained at 28.0°C. During recording, slices were placed on MEA substrates (MED-P515A, Alpha MED Scientific, Ibaraki, Japan) and covered with nylon mesh and a stainless slice anchor. An alternative way to couple a brain slice to an MEA substrate is to use cellulose nitrate as a coating pre-treatment (Egert et al., 2002).

Each of the 64 Chs in the array covered an area that was 50 × 50 μm², and the distance between the centers of the adjacent sites was 300 μm. Each brain slice was carefully located on an MEA substrate, and the arrangement details between an electrode array and the brain slice including the auditory cortex were described in our previous report (Yamamura et al., 2017). The top electrode row was always located at cortical layer 1, so that layer 6 was always located in the fourth or fifth row from the top electrode row.

In each slice on the MEA probe, evoked LFP responses driven by ultrasound stimulation were recorded at a sampling rate of 20 kHz, and the signals were filtered using a range of frequencies between 100 Hz and 10 kHz. The locations of all recording positions on the multielectrode substrate were digitally imaged before and after recording.

Before conducting the main recordings, we recorded current-driven responses to examine the stability of the brain slices. We applied a bipolar square-pulse stimulus (intensity, 50 μA; duration, 200 μs) from one electrode located in layer 2/3 or layer 4: either one of Chs 11–14 (second row in the 64-ch array) or one of Chs 19–22 (third row). To monitor the excitability of brain slices in our main sessions, we performed alternating ultrasound or bipolar current square-pulse stimulation (the same as in the test session: intensity, 25 μA; duration, 200 μs) from one electrode located in layer 2/3 or layer 4 to the slice; one of Chs 11–14 and Chs 19–22. Each stimulation block consisted of 20 trials with the same stimulation conditions and a time interval of 5 s. The pair of blocks was repeated at least three times (total of 30 trials) with the same stimulation parameters. Overall, we obtained 351 ultrasound-evoked responses and 1,069 current-evoked responses for seven brain slices obtained from five animals. The number of current-evoked responses was around three times larger than the ultrasound-evoked responses because the current-evoked responses included those obtained from the test sessions prior to the main sessions, as well as post sessions conducted after the main sessions for the control experiments.

To distinguish ultrasound-driven responses from stimulation artifacts, tetrodotoxin (TTX), one inhibitor of voltage-gated sodium channels, was applied to the ACFS solution perfused in some experiments. The application of 1-μM TTX in the ACSF solution significantly reduced ultrasound-driven LFP responses or completely blocked the evoked activity (Supplementary Figure 1).

To classify a current-driven LFP response pattern evoked at one stimulation site in one cortical layer (layer 2/3 or layer 4) as one of several clusters, we considered the individual pattern as a matrix A of 8 × 8 (64) elements. Peak amplitudes were recorded at each electrode site during the period 100 ms after the ultrasound stimulus onset and during the period from 3 to 40 ms after the current stimulation onset. Each matrix A was subsequently transformed into a corresponding vector (a) located in a multidimensional space. To achieve this, we first constructed matrix A for each response pattern, as follows. When an electrical stimulation was applied at an electrode site, for example, channel p (p = 1, …, 24, i.e., channel numbers from the top electrode row) in layer q (q = 2/3 or 4), the peak amplitudes of the LFP responses at 63 sites were detected as non-negative values, and the element value of the stimulation site was set to be zero. We denoted the elements of matrix A as A_m,n, where subscripts m and n represent the row and column numbers of the electrode array, respectively, and are both integers (i.e., m, n = 1, …, 8). Then, we denoted a labeled response vector, a_p, as a vector of all the elements of A with the label p (one stimulation channel). This response vector represents the specific layer and channel where the stimulation site is located. Similarly, the ultrasound-driven response pattern was represented as a matrix B with 8 × 8 (64) elements and the peak amplitudes of the LFP response, and the corresponding variable b was obtained as a vector of all the elements of B.

After we obtained a spatiotemporal response pattern as the vector variable a_p (or b), we next applied a pattern classification method. In this study, we used a hierarchical clustering method (Duda et al., 2001). We constructed a set of unlabeled sample vectors from all of the vectors a_p (or b) obtained in our experiments. To perform clustering, we followed the following three steps (Duda et al., 2001). (i) To determine the similarity (dissimilarity) between each vector pair, the Euclidean distance was calculated. (ii) Vector pairs that were in close proximity were linked to each other, and a binary, hierarchical cluster tree (dendrogram) was constructed. (iii) Finally, a branch point was determined in the tree, so that all vectors below the branch point in the direction of singleton clusters were grouped into individually identical clusters. In this study, we classified the vectors a_p and b obtained from the current-driven and ultrasound-driven responses into seven and five clusters, respectively. The centers of the clusters were determined via k-means clustering (Duda et al., 2001). As usual, the Euclidean distance (d_euc) was defined as follows:

In this study, the similarity values (I_sim) were defined as I_sim = (1 - d_norm), where d_norm is the normalized value of d_euc with respect to the maximum for all distances.

In the initial state of clustering, we assumed two clusters and calculated the inter-cluster distance between the two geometric centers of the clusters. Then, while the number of clusters was increased, the inter-cluster distances were repetitively calculated for the corresponding cluster centers in all cluster combinations. If all inter-cluster distances were over a threshold value (0.5), the repetitive calculation was terminated, indicating that the distances between any two clusters were relatively further apart, compared with the initial state. Although the threshold value was determined ad hoc, shifting the threshold value ±0.1 did not profoundly change the terminated number of clusters. Additionally, to visualize the clusters consisting of individual response patterns, we computed a 3D subspace using principle component analysis (PCA) and constructed the corresponding vectors in the reduced dimensional space by projecting all original vectors onto the three principal components in the PCA space (Szymanski et al., 2011). As is usual for PCA, we calculated the contribution ratios, which represent the importance of each component in the data set. In our PCA calculation, the contribution of the first three principal component variables was more than 74%.

To characterize the similarity between any two response patterns, we calculated the correlation (R) of a pair of two centroid vectors c_a and c_b in the PCA space with the first three principal components. Here, the correlation R between two vectors was calculated as

where c_ai and c_bi denote the i components of vector c_a and c_b, respectively.

Furthermore, we analyzed the pattern similarity between an ultrasound pressure distribution obtained from the numerical simulation and the LFP response peak values, under the condition that the centers of the transducer and the MEA recording area were out of alignment. To characterize the similarity between the ultrasound pressure distribution and the LFP response peak values, we constructed a matrix V (x₀, y₀) from the numerically simulated pressure distribution I on an MEA probe surface. The position of (x₀, y₀) represents the center of the recording area over the MEA probe, while the center of the acoustic pressure distribution is always located at the origin of the 2D space (i.e., x = 0 and y = 0 mm). When the center position of the recording area is at (x₀, y₀), the pressure distribution I’ centered at that position was described as I’(x, y) = I (x + x₀, y + y₀), indicating just the shift in the position of the pressure distribution center. Because the centers of the recording electrodes are 0.30 mm apart, an element (U_m,n) in the electrode pressure matrix U is defined as

where Δx and Δy are inter-electrode intervals (i.e., Δx = Δy = 0.30 mm), and m and n are row and column numbers (both integers, i.e., m, n = 1, …, 8). After normalizing the maximum value of the matrices U (x₀, y₀), we finally obtained the normalized matrix V (x₀, y₀).

The data are based on experiments conducted in seven brain slices that included the auditory cortex from five C57BL/6J mice. Data are given as the mean ± standard error of the mean (SEM). Statistical analysis for comparisons involving multiple groups was performed using an ANOVA for multiple comparisons followed by a post-hoc Tukey–Kramer test. P < 0.05 was taken to be statistically significant. For template matching, we used the Python-wrapped image analysis library (OpenCV).

We first explored the size parameters (d and l in Figure 1Aa) of the waveguide structure to ensure that it satisfied the conditions of the peak intensity (≥10.0 W/cm²) and half width (≤2.1 mm) on the surface of an MEA probe driven by ultrasound stimulation to the backside of the MEA probe. The waveform of the transducer pressure was sinusoidal, with an amplitude of 100 kPa. The typical pressure distribution for the cross section of the waveguide and the MEA chamber is shown in Figure 2A. For a fixed waveguide size of l (l = 15, 20, or 25 mm), the relationships between peak intensity (I_sppa of Eq. 5) and half width as a function of the diameter (d) of the waveguide mouth are illustrated in Figures 2Ba,b, respectively.

When the diameter d of the waveguide mouth was increased, the peak intensity first increased and was almost saturated over 4 mm (Figure 2Ba). Similarly, when the diameter was increased, the half width was increased (Figure 2Bb). The amount of ACSF solution in the MEA chamber was varied with ten height levels (h) ranging from 1.0 to 10 mm in 1.0-mm steps; the error bars in the two plots show the standard deviations.

Among the several candidates for the two size parameters (d and l), we finally selected d = 4 mm and l = 20 mm as the 3D model parameters. These satisfied the two conditions: the peak intensity was 10.1 W/cm, and the half width was 2.03 mm. Furthermore, under the same size parameter conditions, the focal length [peak position in the vertical (z) axis] was 0.9 mm from the bottom of the MEA probe (0.7 mm thickness), which indicates that the focal point was located within the brain slice (0.4 mm thickness).

To examine pressure distribution patterns in a brain slice, we simulated the corresponding 3D model (Figures 3A,C) using the FEM. In the 3D model, a rectangular brain slice of 4.0 × 3.5 × 0.4 mm³ was assumed to be located in the center of the MEA recording area (Figures 1Ab, 3A, 4B). As a waveguide model, the size parameters of d = 4 mm and l = 20 mm were used, which was determined on the basis of the 2D model simulation. The vertical cross section of the pressure distribution in response to sinusoidal ultrasound stimulation (100 kPa) is illustrated in Figure 4A; the top height (h) of the ACSF solution in the MEA chamber was 6 mm. The result obtained from 3D modeling (Figure 4A) was similar to that obtained from 2D modeling (Figure 2Aa), although the colormap plots appear different at first glance.

The horizontal cross sections of the pressure distribution at z = 0 mm (the bottom plane of a brain slice) in Figure 4Aa are illustrated in Figure 4B. For waveguides that are the same size (d = 4 mm and l = 20 mm) with different amounts of ACSF solution in the MEA chamber (h = 3, 6, and 9 mm), three typical examples of pressure distribution patterns for a brain slice are shown in Figure 4B. Figure 4Bb (h = 6 mm) shows several decaying local peaks along the radial axis, with an approximately 2-mm interval on the MEA chamber. This result suggested that concentric standing waves induced by the ultrasound stimulation would be observed around the center of the MED chamber in the actual pressure measurement experiment.

We found that the pressure distributions details depended on the height of the ACSF solution in the chamber (Figures 4B,C). For example, the amount of ACSF solution increased with the number of local peaks around the central peak (Figure 4C). Therefore, the spatial modes (local peak intervals) of acoustic pressure distributions may differ according to the height of the ACSF solution in the MEA chamber. In addition, the relationship between the intensity of the peaks (the largest and second largest peaks) and the height of the ACSF solution was non-monotonic (Figure 4Cc), indicating that this relationship could have a greater impact on the complexity of ultrasound-induced neural activity than the changes in ACSF.

To examine the spatial irradiation properties of ultrasound waves traveling through the custom-designed waveguide attached to an ultrasound transducer, we measured the fluid pressure in the MEA chamber using a needle hydrophone located 1.0 mm from the surface of the MEA. The hydrophone was vertically fixed to the MEA probe in the chamber (Figure 5A). During the measurement, the height of the ACSF solution in the chamber was 3.0 mm.

For one quadrant of the circular chamber, the measured pressure is spatially shown in Figure 5Ba. We compared the measurement with the simulation using the simulated pressure in the corresponding area under the same conditions, as shown in Figure 5Bb. Although the radial symmetry of the pressure distribution was clear in the numerical simulation, the symmetry was not evident in the measurement. While the measurement (in mm) traveled from the origin (0, 0) to (4.5, 4.5) for the sensing area, the time lags among the individual measurements on the 2D plane might have influenced the pressure sensitivity. The difference may have been caused by dewatering between the top of the waveguide and the glass MEA substrate, which was conducted to reduce the impedance mismatch between them.

The measured and simulated pressures in the radial axis are superimposed in Figure 5C. In the simulation, the sensor (hydrophone) position was varied with respect to the MEA substrate from 0 to 1.0 mm. Figure 5C indicates that, under the condition in which the sensor position was 1.0 mm above the surface of the MEA substrate, the second largest peak in the measured pressure (black dot) was around 2 mm away from the largest peak. In contrast, under the same conditions, the second largest peak in the simulated pressure (blue curve) was around 3 mm away from the corresponding largest peak. In the simulation, the second peak positions shifted toward the largest peak position as the distance between the sensor and the MEA substrate decreased: e.g., 0.5, 0.2, and 0 mm. This implies that the measurement distance between the sensor and MEA substrate might be smaller than 1.0 mm, although the measured pressures were larger than the simulated pressures. In the pressure measurement, the use of a phantom slice (or brain tissue block) may give a more accurate result consistent of the numerical simulation and an appropriate representation for the biological experiment.

Before and after the ultrasound stimulation experiment, we recorded the activity of brain slices via electric current stimulation, as described in the experimental protocol. As a typical example, when a short current pulse at a low intensity (25 μA) was applied to a brain slice that included the auditory cortex on the MEA probe in a fixed arrangement (Figure 6A), we obtained highly reproducible extracellular LFPs (Figure 6B).

For the recordings, we selected a single stimulation electrode (e.g., Ch 13 in Figure 6A), and the electrode site was located at layer 2/3 or layer 4 (Figure 6A). In layer 2/3, typical LFPs in auditory cortex slices during a 50-ms time window after stimulation onset had three components: (i) an early negative-going potential (component A) from the stimulation onset to approximately 3 ms afterward; (ii) a late positive-going potential (component B) following component A; and (iii) a final decaying potential that approached the initial baseline (component C) (Figure 6C; Yamamura et al., 2017). If the stimulation channel was at layer 5/6, the evoked responses propagated into the upper layers (data not shown; Yamamura et al., 2017). Consequently, the response amplitudes were maximized in layer 2/3 (Yamamura et al., 2017). In this study, we only selected single channels at layer layer 2/3 or layer 4 for each stimulation trial.

To characterize the 1,069 LFP patterns obtained in response to current stimulation, we converted them to pattern vectors a_p and classified the vectors into seven clusters (Figure 7Aa). All inter-cluster distances were over the pre-set threshold value (>0.5), indicating that the distances between any two patterns belonging to different clusters were relatively far apart. Moreover, because we selected stimulation sites at layer 2/3 or layer 4, the maximum amplitude responses were usually obtained from the recording sites in the second, third, or fourth row in the MEA array (Figure 7Ab). A typical LFP response and the pattern matrix A from each cluster are illustrated in Figure 7B. In particular, pattern pairs between clusters 1 and 2 (1, 2), 3 and 4 (3, 4), and 6 and 7 (6, 7) had larger similarity index values (i.e., smaller Euclidean distances). All of these results, including the number of patterns in each cluster, are summarized in Supplementary Table 1.

In short, the obtained response patterns in individual clusters reflected the corresponding stimulation sites and the propagation properties of current-evoked activity in local cortical networks. Identical stimulation sites did not always induce the same response patterns, implying that the propagation properties of individual brain slices varied to some extent.

Next, we examined whether ultrasound stimulation could actually evoke neural activity in a brain slice when our MEA-based recording system was combined with the custom designed waveguide. Ultrasound-driven LFP responses were recorded among seven examined slices. See Figure 8A for examples of typical spatiotemporal ultrasound-driven responses for an ultrasound stimulation with a 100-ms duration, delivered to the same slice shown in Figure 6. For this brain slice, larger responses were locally observed in layer 5/6, although cortical layer-specificity of evoked responses was not always easily characterized. This example demonstrates that the responses were not only obtained from a small local area (e.g., Chs 26, 27, and 34), but from several response subareas (Chs 14 and 22) as well. In addition, with respect to the time courses of LFPs, responses with a very slow transient component were obtained from several electrodes after the ultrasound stimulation onset (Figure 8A). The time scale for typical ultrasound-driven responses was approximately 10 times larger than that for current-driven responses (i.e., Figure 6 vs. Figure 8). As an exception, faster LFP responses were obtained from one brain slice (see below). Furthermore, for ultrasound stimulation with a 100-ms duration, typical LFPs in auditory cortex slices evoked during the 500-ms time window after stimulation onset were also composed of three components (Figure 8B): (i) an early negative-going potential (component A) from the stimulation onset to approximately 100 ms afterward; (ii) a late positive-going potential (component B) following component A; and (iii) a final decaying potential (component C) approaching the baseline level. Here, we characterized the LFP waveforms using two response properties: peak latency and peak amplitude (Figure 8B).

The peak amplitudes increased with the acoustic pressure levels of the ultrasound stimulation (Figures 9Aa–d) for the two stimulation durations (100 and 200 ms; Figure 9A). In response to the ultrasound stimulation with a 200-ms stimulus duration, the differences in the peak amplitude of the ultrasound-driven LFP responses were significant between the smallest pressure (110 kPa) and the largest pressure (410 kPa) levels (paired t-test, P = 0.042, Figure 9Ba), indicating that the pressure level of the ultrasound substantially impacted the LFP response. In addition, the pressure levels of the ultrasound stimulation with a 100-ms duration did not influence the peak latencies in the early phase of the evoked response, and the latencies were around 100 ms (Figure 9Bb). In contrast, as the pressure levels of the ultrasound increased with the 200-ms stimulation duration, the latencies of the large negative peak shifted to around 100 ms (Figure 9Bb). This result also confirmed that the pressure levels of the ultrasound stimulation closely reflected the LFP responses, and that they more strongly influenced the peak amplitudes than the peak latencies. Furthermore, the duration of the ultrasound stimulation (100 vs. 200 ms) influenced the duration of the LFP response, and the response durations were such that the temporal lengths of the LFP response amplitude >50% of the corresponding peak amplitude. For the pressure of 310 kPa, for example, in response to the ultrasound stimulation with the durations of 100 ms and 200 ms, the LFP response durations were 69.0 ± 20.4 ms and 128.0 ± 33.5 ms, respectively. This result indicates that the LFP response durations, which were significantly different, were dependent on the ultrasound durations (paired t-test, P = 0.028). Therefore, for the examined ultrasound range, the temporal integration effect of the ultrasound stimulation was proportional to the stimulation duration.

To characterize the 351 LFP patterns obtained in response to the ultrasound stimulation, we converted them to pattern vectors b and classified the vectors into five clusters (Figure 10Aa). All inter-cluster distances were over the pre-set threshold value (>0.5), although some clusters seemed to be close to one another (e.g., clusters 1 and 2). For each response pattern, the number of sites for which the maximum amplitude was on a negative peak was counted and accumulated on a 2D electrode array, and the results indicated that the sites of the peaks driven by ultrasound stimulation were more dispersed (Figure 10Ab) than those with current stimulation (Figure 7Ab). A typical LFP response and a pattern matrix B from each cluster are given in Figure 10B. In cluster 3, the time courses of the LFP responses obtained from one brain slice were faster than the response time-courses from the other clusters, and the response latencies were longer after the stimulation onset (Figure 10Bf). In cluster 5, negative-going responses were obtained for almost all sites, implying that the ultrasound stimulation was effective for evoking neural activity. Therefore, the obtained response patterns in the individual clusters appeared to reflect the corresponding ultrasound pressure distributions and propagation properties of ultrasound-evoked activity in particular local networks. All of these results, including the number of patterns in each cluster, are summarized in Supplementary Table 2.

The correlation coefficients between two response clusters obtained from current- and ultrasound-driven stimulations are illustrated in Figure 11. We created dendrograms of the clusters obtained from the individual stimulation methods by characterizing the similarities among the cluster pairs for each stimulation method (Figure 11). Although the response times differed in size by approximately 10-fold between the two stimulation methods, some clusters had similar response properties in the 2D space across the two stimulation methods. In particular, three cluster patterns (clusters 1–3 in both stimulations) were closely matched (similarity index ≥ 0.794 in Figure 11). However, the patterns in cluster 5 showed dissimilarity (similarity index < 0.3) with all response patterns obtained by current stimulation.

During the experiments, we were careful to ensure that the center of the transducer coincided with the center of the recording area on the MEA substrate. However, the pressure distribution measurements suggested that the center of the pressure distribution did not always match the exact center of the MEA substrate (Figure 5Ba). Therefore, we examined how discrepancies in the center position could affect ultrasound pressure distributions and pattern similarities. Here, we assumed that the one-site neural response driven by ultrasound stimulation would be proportional to the ultrasound pressure distribution at the corresponding electrode site. In our numerical calculation, we artificially simulated 49 (7 × 7) pressure distribution patterns, in which the center position (x₀, y₀) of the ultrasound transducer was simply shifted by a step of 0.5 mm relative to the center of the MEA substrate (Figure 12A). We found that ultrasound-driven clusters 1, 2, and 3 were most similar to the pressure distribution patterns (D 23 and D 11 in Figure 12A). Their transducer centers were 1.0 mm leftward [i.e., transducer-center position, (x₀, y₀) = (−1.0, 0) in mm] and 1.0 mm upward [(x₀, y₀) = (0, +1.0) in mm]. The correlation coefficients (R) of D 23 vs. Cl 1 and D 11 vs. Cl 3 were 0.546 and 0.504, respectively. The correlation coefficients between all of the ultrasound clusters (Cl 1 to Cl 5) and all artificial pressure distributions (D 1 to D 49) are given in Figure 12Ba. These results of the numerical calculation and the correlation analysis indicate that the alignment between the centers may influence the spatial patterns of ultrasound-driven neural responses. Therefore, a measurement system that is capable of fine alignment is desirable for future research.

A variety of temporal patterns have been used for ultrasound stimulation of the brain (Rinaldi et al., 1991; Bachtold et al., 1998; Khraiche et al., 2008; Tyler et al., 2008; Choi et al., 2013; Menz et al., 2013; Kim et al., 2017; Kubanek et al., 2018; Prieto et al., 2018; Menz et al., 2019; Supplementary Table 3). In this study, we selected a simple sinusoidal oscillation (frequency, 500 kHz; duration, 100 or 200 ms) as the stimulation waveform. Many indices of intensity have been used to characterize the multiple time scales of ultrasound stimulation patterns (Supplementary Table 3). Here, we used the spatial-peak pulse-averaged intensity (I_sppa) as one such index, simply defined as the average intensity of a sinusoidal wave pulse.

Repetitive bursts of short ultrasound oscillation pulses are often selected as another common stimulation pattern (Tyler et al., 2008; Choi et al., 2013; Kim et al., 2017). For these, the spatial-peak temporal-averaged intensity (I_spta) is commonly used as an index, defined as the average intensity over the total duration of each stimulus. Burst parameters, including the burst duration and interval used in I_spta, are not relevant to the present study. However, the effect of repetitive ultrasound bursts on neural activity in brain slices is of interest for future work because these bursts elicit stronger neural activation (Bystritsky et al., 2011; Blackmore et al., 2019). The stimulation intensities and carrier frequencies used in previous in vitro preparations are summarized in Supplementary Table 3.

To the best of our knowledge, this study is the first to report ultrasound stimulation-evoked responses in acute cortical slices from mice. In terms of diagnostic ultrasound imaging devices for human applications, the Food and Drug Administration (FDA) guidelines provide the following safety criterion regarding stimulation intensities: the I_spta must not exceed 720 mW/cm² and the I_sppa must not exceed 190 W/cm² (Duck, 2007). These intensities are much larger than those used in our experiments (I_sppa ≤ 5.68 mW/cm²), although in vitro and in vivo situations may vary in terms of acoustically reflective surfaces and the surrounding environment (i.e., in vivo vs. glass MEA probe, recording chamber, and brain slices in experimental ACSF solution). Understanding the limitations associated with in vitro systems is important when interpreting our experimental data and considering how our findings will be applicable to future work. In particular, the standing ultrasound waves that we detected in the MEA chamber could cause complex distributions of ultrasound intensity within slices, and this could reflect the observed response activity (e.g., Figures 4B,C).

Although neural impulses have historically been considered to be electrical signals, the super-threshold depolarization mechanisms of neural membranes involve a combination of electrical, mechanical, chemical, and conformational changes in the excited cells. Thus, describing a neural impulse via a mechanical pathway could support a physical basis for ultrasound-driven neuromodulation. Several recent review papers (Blackmore et al., 2019; Pasquinelli et al., 2019; Cardenas-Rojas et al., 2020) summarized four potential mechanisms by which ultrasound signals could trigger action potentials: (i) the generation of capacitive currents as a result of membrane displacements, (ii) the activation of mechanosensitive channels, (iii) sonoporation in the lipid bilayer, and (iv) coupling with membrane waves along the axon. Effective ultrasonic neurostimulation may be possible by combining these mechanisms. In particular, to enable sufficient radiation force for larger ultrasound amplitudes, the conformational change could directly polarize the membrane potentials via mechanocapacitive coupling (Zecchi et al., 2017), resulting in rapid excitation without a time delay. In contrast, at lower ultrasound amplitudes such as those used in this study, the radiation force is weak, which may depolarize the membrane in the presence of very slowly depleting ion gradients. Therefore, the duration of a spike will need to be sufficiently long for polarization to cross the threshold required to produce an action potential.

As reported in our previous studies (Namikawa et al., 2017; Yamamura et al., 2017), we reliably observed cortical layer-specificity of evoked responses among brain slices that included the mouse auditory cortex. In this study, we recorded spatiotemporal LFP patterns in response to current stimulation via a single electrode at layer 2/3 or layer 4 of mouse auditory cortical slices. A short current-pulse activated LFP responses around the stimulation site, and the activity propagations were locally restricted (Figure 6). In addition, the LFP patterns could be classified into seven clusters (Figure 7A) using a pattern classification method known as hierarchical clustering. Thus, the clustered patterns were closely dependent on the stimulation sites (Figure 7B) and the laminar properties of the cortical slices.

We recorded and analyzed activity patterns for MEA sites at which ultrasound-driven LFP responses were evoked (Figure 8), and found that the activity patterns could be classified into five clusters (Figure 10A). That ultrasound-driven LFP responses were more complex than the response patterns for current stimulation (Figure 10B vs. Figure 7B), implying that the ultrasound pressure distribution itself was complex (e.g., Figure 4B). Furthermore, among the two response cluster sets, we compared the response patterns driven by current and ultrasound stimulations. We found that three cluster patterns (clusters 1–3 in both stimulations) were closely matched (similarity index ≥ 0.794 in Figure 11). In addition, we compared the response patterns with the patterns of ultrasound pressure distributions, and found that some pressure patterns were matched with ultrasound-driven response patterns (Figure 12). Thus, the response similarity appears to be influenced by the laminar properties of cortical slices and/or spatial patterns of ultrasound stimulations. We plan to investigate this experimentally in our future work.

One factor over which we had limited control was the volume of the ACSF solution flowing in the MEA chamber. We are planning to construct a feedback control system to regulate this in our future work. Using such a control system, we expect to find that ultrasound-driven cortical responses have more similar activity patterns among different cortical slices, reflecting the relationship between ultrasound pressure patterns and neural response patterns.

In addition, in our future MEA-based recording system, we hope to improve the accuracy with which the center axis of the transducer can be aligned with the center of the electrode recording area. We found that a short shift (≤1 mm) in the distance between the two centers led to different ultrasound pressure distributions (Figure 12A). Therefore, pressure measurements will be conducted before ultrasound stimulation in future procedures. Furthermore, we hope to significantly reduce the measurement time with respect to hydrophone sensitivity in a recording area over a 2D plain.