Two adult Sprague–Dawley rats (male: 270–310 g) were purchased from the Laboratory of Animal Center of Zhejiang University. Both rats were anesthetized with chloral hydrate solution (10%, 0.35 ml/100 g). The anesthetized rats were fixed on a stereotaxic apparatus (51603U, RWD, China) with ear bars and a clamping device to keep their heads horizontal. Specifically, the fur covering each rat’s skull was trimmed to expose the cleaner scalp, and the scalp was sterilized with iodine. The bregma and lambda were explored after cutting the scalp and removing the subcutaneous tissue. One hole was drilled into the rats’ meninges to fix reference electrodes. Another hole was drilled at the intersection of the sagittal suture and the interaural line as the ground. Five recorded sites were determined according to the reference (Paxinos and Watson, 2006) (primary sensory cortex, SI: AP = −4.44 mm, ML = 6 mm, DV = 3 mm; primary motor cortex, MI: AP = 0.48, ML = 2.2, DV = 2; lateral habenular nu, LHb: AP = −2.52, ML = 0.6, DV = 4.8; field CA1 of the hippocampus, CA1: AP = −3.96, ML = 1.2, DV = 3.4; field CA2 of the hippocampus, CA2: AP = −5.20, ML = 5.6, DV = 6). Figure 1A shows the anatomy map of the rat. It visualizes the sites of LFP recorded under different ultrasound stimulation parameters. The electrode and the head wound of the rat were sealed with dental cement, and only the interface of the recording electrodes was exposed finally. After the animal surgery, the two rats were housed in cages for a week to recover with free access to food and water. This study’s protocols were approved by the Laboratory Animal Ethical and Welfare Committee of Shandong University Cheeloo College of Medicine.

In order to rule out accidental circumstances, we repeated the experiment three times for each stimulus parameter on two adult Sprague–Dawley rats.

In the TUS system, an arbitrary function generator (YB16051A, Lv yang, China) was used to control the pulsed repetition frequency (PRF), fundamental frequency (FF), number of cycles per pulse (NC/p), and duty cycle (DC) of ultrasound. The pulsed sequence from the generator was amplified by a high voltage amplifier (ATA-2021H, Aigtek, China) and transmitted to an ultrasound transducer (focal length of 45 mm and focus diameter of 20 mm). The ultrasound transducer was connected to the rat skull by a conical acoustic collimator filled with ultrasound coupling gel (Figure 2A). The collimator is a circular area with a radius of 0.4 cm, so the focus area does not exceed 0.51 cm². The transducer connected to the collimator is 8 mm higher than the skull, so the focal length into the brain is less than 37 mm. As illustrated in Figure 1A, in order to reduce ultrasonic attenuation, electrodes were implanted into the nucleus vertically, and the collimator was aimed horizontally from the occipital side. The angle between them was 90°. The electrode and the head skull of the rat were all sealed with dental cement. In order to avoid the obstruction and attenuation of ultrasound by dental cement, we applied ultrasound to the occipital lobe. The stimulation site was chosen as shown in Figure 1A, 6.5 mm lateral of midline and 4 mm posterior of bregma.

Pulsed ultrasound sequences were constructed as illustrated in Figure 2B. The FF and DC of the ultrasound were 500 kHz and 20% respectively. By using the collimator, the ultrasonic focal spot radius was fixed at 0.4 cm. Under the condition of fixed FF and DC, we designed three groups of trials by changing PRF and NC/p. For each trial, we changed the magnification of the amplifier to get the increasing AI. Three PRFs (500 Hz, 1,000 Hz, and 1,500) were used in this study. The range of I_SPTA was from 7.2 to 25.6 mW/cm². The ultrasound parameters, spatial peak temporal average intensity (I_SPPA), spatial peak pulse average (I_SPTA), PRF, and NC/p we used are listed in Table 1.

The ultrasonic power generated by the plane transducer (installed on the collimator) was measured by the radiation force method using the ultrasound power meters model (UPM-DT-1AV, Ohmic, United States). During this process, the transducer was held above the de-gassed water by a positioning clamp. The ultrasonic energy passed through the water to reflect off the conical target and was then absorbed by the rubber lining. The radiant power was proportional to the total downward force in the conical target. This weight was then transferred to the electro-mechanical load cell, which could produce a digital readout in watts of power. In this experiment we assumed that the brain and tissues had the sound, speed, and density of water. The final acoustic attenuation transmitted into the brain was 83.9% because of the absorption and refraction properties of the ultrasound in the rat skull. Since the five nuclei are not far apart, the ultrasound intensity at different nuclei is the same. Thus, we only consider the attenuation of a skull. Details of the AI calculation are described in the Appendix.

Five electrodes (nickel-chromium alloy, 35 μm, A-M SYSTEM, United States) were used to record the LFP signals in the designated brain area of rats. Each electrode was made of two nickel-chrome wires that were intertwined with each other. The reference electrode implanted into the rats’ meninges was made in the same way as the recording electrodes. The ground electrode was made of exposed silver wire. Finally, all electrodes were connected to a 16-channel (2 × 8 array) interface.

The raw LFP signals produced in response to TUS were collected by a 32-channel neural signal processor (RHD2132, Intan, United States), which can perform pre-amplification and AD conversion at the same time. These signals were acquired at the sampling frequency of 1 kHz and stored on field programmable gate array (FPGA). Finally, all data were transmitted to the computer and processed using MATLAB 2018b (MathWorks, Inc., United States).

After the final recording session, the animal was perfused with paraformaldehyde (PFA) 4%, and histological analysis was performed to confirm electrode locations. Correct electrode position was confirmed for all animals in this study.

Spectral estimation can transform a signal from the time domain to the frequency domain and provide a description of the signal’s power distribution along with frequency (Hu and Zhang, 2019). The foundation of spectral estimation is the discrete Fourier transform (DFT), which is calculated using the Fast Fourier transform (FFT) algorithm for computational efficiency. Taking discrete-time signalx_n, n = 1, 2…N as an example, the DFT ofx_ncan be written as

Based on DFT, the power spectral density (PSD) can be calculated to describe how the power of signal is distributed over different frequencies. In our study, LFP signals were acquired at the sampling frequency of 1 kHz with software designed Bandpass (1–140 Hz, 1–45 Hz) and notch (50, 100 Hz) filter. And LFP powers were computed for individual channels using the multitaper spectral estimation technique (Babadi and Brown, 2014). These powers were then normalized by mean power for the global spectrum (1–140 Hz) within each time serial.

To observe the time-variant spectral components of LFP signals over whole periods, TFA techniques are needed to study such non-stationary signals in both the time and frequency domains simultaneously. The S-transform expresses the phase information by using the Fourier kernel, and translates the amplitude envelope of the Gaussian window. The S-transform of a time-varying LFP signal x(t) is defined as:

The relative power represented a key power of neuromodulation parameter (Henrie and Shapley, 2005). The upper frequency cutoff of the LFP is often considered to be around 200 Hz. In our study, after spectrum estimation of LFP signals, signals were decomposed into the delta (1–3 Hz), theta (4–8 Hz), alpha (8–12 Hz), beta (12–30 Hz), slow gamma (30–60 Hz), high gamma (60–200 Hz) frequency bands by the discrete wavelet transform (DWT) according to the rhythms (Buzsaki, 2006). Relative power was used to eliminate the difference between experiments. The relative power of each rhythm is equal to the corresponding band-limited power corrected by the global spectrum (1–200 Hz) power. For example, the band-limited power of the alpha rhythm is calculated as the spectral estimation averaged from all frequency points within the range of 8–12 Hz, and the relative power is the ratio between the alpha power and the total power of the LFP signal.

Enhanced neural oscillation (the average amplitude of LFPs were 4–10 times higher than normal) with rhythmic patterns appeared in selected nuclei during VLUTS with PRF of 500 and 1 kHz, when the I_SPTA raised to 12 mW/cm² (Table 1 and Figure 3). Under continuous ultrasonic stimulation, enhanced voltage presented such increase-subside rhythm repeats periodically with very slow frequency, around 0.03–0.05 Hz, which bears some similarity to autorhythmicity. We define these types of oscillations pattern to be auto-rhythm. As shown in Figure 3A, the red line marks one start and end period of CA2 nucleus rhythm. The number of responding nuclei reached maximum when I_SPTA raised to 16.4 mW/cm² for PRF 500 Hz and 20.8 mW/cm² for PRF 1 kHz (Table 1). And the number of responding nuclei started to decrease while I_SPTA continued to increase. All such auto-rhythmic oscillations were highlighted with shadows (Figures 3A,B). The average amplitude of LFP showed in Figure 3 was listed in Table 2. Under PRF 1.5 kHz, neural oscillation was moderately increased during stimulation but without the auto-rhythm.

The LFP signal is a continuous stochastic process. It consists of a series of continuously varying voltages in time, which encodes rich information of neural activities and can be quantitatively characterized in the frequency-domain by spectral estimation. The LFP signals of individual channels before and after the ultrasound stimulation were processed by the multitaper spectral estimation technique as illustrated in Figure 3. Compared to neuronal activity prior to ultrasound stimulation, all five nuclei maintained higher excitability with the auto-rhythmic pattern during VLTUS with PRF 500 Hz, I_SPTA = 16.4 mW/cm² and PRF 1 kHz, I_SPTA = 20.8 mW/cm², and increased high-frequency energy (slow gamma and high gamma) was observed. Intriguingly, the CA1 nuclei didn’t show auto-rhythmic pattern during VLTUS with PRF 1 kHz (Figure 3B). Under PRF 1.5 kHz stimulation, the low frequency delta wave (4 Hz) was dramatically increased (Figure 3C).

Table 1 shows numbers of nuclei which showed significant neuromodulations with auto-rhythmic pattern under different UTS parameters. Under PRF 500 Hz, all five nuclei showed the positive response at I_SPTA = 16.4 mW/cm². Under PRF 1 kHz, four nuclei responded at I_SPTA = 20.8 mW/cm². With increase the AI, the number of responding nuclei started dropping. When I_SPTA raised to 25.6 mW/cm² only two nuclei under PRF 500 Hz, three nuclei under PRF 1 kHz showed positive response. Under PRF 1.5 kHz, no auto-rhythmic pattern has been observed.

We measured the average power before and after stimulation (10 s each). Results are shown in Table 2. The average power of normal LFPs was within 0.07 μW, but when PRF was 500 Hz and 1 kHz, the average power of LFPs increased by 13–40 times. It can be seen that the transcranial ultrasound induced neurons to open ion channels and discharge. When the PRF is 1.5 kHz, the power slightly increases with obvious group slow wave, but there is no rhythmic pattern. Definitions of ultrasound parameters please see Table 3.

We utilized the TFA technology to analyze these phenomena further (Figure 4). To avoid high computational complexity, we only analyzed the LFP signal of five nuclei in 1 epoch (2S), starting from 4 s after stimulation to eliminate the noise at the beginning of stimulation. The normal LFP signal energy is mainly concentrated between 2 Hz and 45 Hz. After stimulation, the high frequency component (60–200 Hz) of the LFP is prominent when PRF = 500 Hz and 1 kHz, and the low frequency (1–8 Hz) appear and dominant when PRF = 1.5 kHz. These phenomena are consistent with the results discussed in Figure 3.

The baseline of the stimulation or inhibition is calculated by 6 s normal LFPs (P0: reference power of normal LFPs). And LFP signals within three epochs (6 s) just 4 s after stimulation (to eliminate abnormal interference of brain at the beginning of stimulation) are used to analyze the relative power. To demonstrate the relative stimulation or inhibition effect of different PRF on each rhythm (compared with normal LFP signal without stimulation), the relative power of LFP signals was divided by the relative power of normal LFP signals, and the logarithm of the result was taken.

As we can see from Figure 5 when PRF 500 Hz or PRF 1 kHz the gamma component of the LFP, especially high-gamma component of the LFP (hgLFP), was significantly increased after stimulation, the theta and delta power decreased otherwise. However, when PRF 1.5 kHz high frequency band are suppressed while delta rhythms increased.

Figure 6 analyzed LFP signal of CA2 nucleus in one epoch (2S) just 4 s after stimulation as an example. Figure 6A represents the decomposition of the LFP signal of CA2 nucleus into different rhythms. And we analyzed the power ration of delta and high-gamma rhythms which is changed obviously in Figure 5 under different stimulation parameters. The high-gamma power ratio (Figures 6B,C) was significantly increased after stimulation, while the theta power suppressed. When PRF 1.5 kHz (Figure 6D) the powers ration of the theta and delta frequency bands were increased obvious, however, while the high-gamma power ratio remains the same.

Here we convincingly prove that the hgLFP is active when PRF 500 Hz and 1 kHz. And the awakening of delta and theta rhythm is a particular phenomenon when PRF 1.5 kHz.

To the best of our knowledge, this is the first report of neurons sensitive to the specific frequency and AI during ultrasound stimulation. And we could find that neurons responding selectively to inputs of a specific frequency.

When PRF 500 Hz or PRF 1 kHz, neurons show prefer at gamma frequencies. The increased power of hgLFP during VLTUS is a kind of neural excitement. The physiological reasons for the elevated power of hgLFP phenomenon are not clear, but previous studies can provide some insights for its source. One potential explanation is that high-frequency (>60 Hz) components of signals represent activation of neuronal populations in the underlying cortex (Crone et al., 2006). Contamination of spike waveforms from nearby cells, either too small or too variable to be detected as single unit spikes (Berens et al., 2008). TUS under PRF 500 Hz or 1 kHz could wake up and activate hgLFP. Recent data recorded from animals and humans have suggested that gamma-frequency activity plays an important role in both attention and long-term memory (Jensen et al., 2007). Such hgLFP with auto-rhythm performance could aid in understanding neuronal processing in complex cognitive functions in future studies.

The slow waves (theta and delta frequency band) appear and occupy when PRF is 1.5 kHz. Theta rhythms are produced by local interactions between hippocampal interneurons and pyramidal cells. There is evidence that mammalian brain display intrinsic resonance with frequency selectivity for electrical stimulation inputs within the theta-range (4–10 Hz) (Hutcheon and Yarom, 2000; Vera et al., 2017). Seidenbecher et al. (2003) have found that synchronization of theta activities in the amygdalohippocampal network represents a neuronal correlate of conditioned fear, apt to improve neuronal communication during memory retrieval. This may be related to the fact that ultrasound causes abnormal neuronal activity, so neuronal resonance of theta and delta rhythms might be waked up, relevant to cognitive operations and learning in the brain.

The mechanism of auto-rhythmic pattern of neural oscillation during VLTUS remains unknown but could relate to the refractory period observed in neurons. Refractory periods are caused by the inactivation property of voltage-gated sodium channels and the lag of potassium channels in closing. Voltage-gated sodium channels have two gating mechanisms, the activation mechanism that opens the channel with depolarization, and the inactivation mechanism that closes the channel with repolarization. While the channel is in the inactive state, it will not open in response to depolarizing stimuli. LFP reflects the neural population activity. The synchronization of refractory periods amongst groups of neurons changes population neuronal activity between depolarized and hyperpolarized phases. This may cause that rhythmic pattern. Interestingly, the rhythmic patterns in five nuclei are not identical (Figure 3), which may be due to physiological features of the neuron and the characteristics of the extracellular field in each nucleus.

All five nuclei showed the auto-rhythmic neuromodulation with high frequency energy during the VLTUS (PRF 500 Hz, I_SPTA 16.4 mW/cm²). All nuclei are in distinct brain networks and located distally, while focus ultrasound stimulation is localized. It’s hard to conclude that all electrical signals were from post synaptic action. Previous studies showed that spatiotemporal fluctuations of the extracellular field may influence via purely electrostatic, so-called ephaptic coupling (Jefferys, 1995; Weiss and Faber, 2010). Network oscillations across and within brain areas are critical for brain function (Roach et al., 2018). Coherent oscillations could provide such a mechanism in coordinating the timing of activity and enhancing the functional link between distant areas. Ephaptic coupling may play a big role in this procedure. Our data support those findings and indicate that the ephaptic coupling could widely exist in brain networks, which strongly entrain action potentials (Anastassiou et al., 2011). Ephaptic coupling has been shown to affect population activity during hypersynchronous epileptic discharges (Jefferys, 1995; McCormick and Contreras, 2001). However, how ephaptic coupling alters the functioning of neurons under physiological conditions remains unclear. Our data suggest that VLTUS could be a valid tool in this field.

A larger sample study will help us better evaluate the neuromodulation effect of VLTUS; Only selective PRF and AI have been tested in this pilot study. Need further experiments to prove the mapping relationship between simulation position with the five recording nuclei; The hydrophone can be used to characterize the acoustic field with the focal length and focal length of the transducer measured accurately.