Abstract
Diverse translational and research applications could benefit from the noninvasive ability to reversibly modulate (excite or suppress) CNS activity using ultrasound pulses, however, without clarifying the underlying mechanism, advanced design-based ultrasonic neuromodulation remains elusive. Recently, intramembrane cavitation within the bilayer membrane was proposed to underlie both the biomechanics and the biophysics of acoustic bio-effects, potentially explaining cortical stimulation results through a neuronal intramembrane cavitation excitation (NICE) model. Here, NICE theory is shown to provide a detailed predictive explanation for the ability of ultrasonic (US) pulses to also suppress neural circuits through cell-type-selective mechanisms: according to the predicted mechanism T-type calcium channels boost charge accumulation between short US pulses selectively in low threshold spiking interneurons, promoting net cortical network inhibition. The theoretical results fit and clarify a wide array of earlier empirical observations in both the cortex and thalamus regarding the dependence of ultrasonic neuromodulation outcomes (excitation-suppression) on stimulation and network parameters. These results further support a unifying hypothesis for ultrasonic neuromodulation, highlighting the potential of advanced waveform design for obtaining cell-type-selective network control.
Significance Statement
Recent studies have demonstrated that ultrasound waves are capable of stimulating or suppressing neural circuits, thereby opening up an important new route toward targeted noninvasive neuromodulation. However, the underlying mechanism for ultrasonically eliciting specific neuromodulatory effects has not been clarified. Our new theoretical analysis reveals that ultrasound can selectively excite different cortical neuron subtypes simply by changing the stimulation pattern, driven by the response properties of T-type calcium channels. Interestingly, the model’s predictions at the single-neuron and network levels are shown to closely agree with and explain the emerging field’s entire body of experimental results, spanning from rodents to humans, and can thus facilitate the development of new ways of treating or diagnosing brain disorders.
Introduction
Both classical and recently emerging data have demonstrated that the interaction of therapeutic ultrasound (US) and excitable tissues leads to a diverse and complex set of reversible physiological phenomena which are collectively referred to as ultrasonic neuromodulation (Bystritsky et al., 2011; Yoo et al., 2011a; Kim et al., 2015; Naor et al., 2016). US-induced neuromodulation phenomena include the generation of action potentials (APs) in CNS neurons in the brain (Tyler et al., 2008; Tufail et al., 2010, 2011; Kim et al., 2012, 2014; King et al., 2013, 2014) and the retina (Naor et al., 2012; Menz et al., 2013), suppression of CNS activity (Yoo et al., 2011a; Min et al., 2011; Kim et al., 2015) and of nerve conduction (Young and Henneman, 1961; Lele, 1963; Colucci et al., 2009), as well as more subtle changes in excitability (Moore et al., 2000; Tsui et al., 2005). The development of ultrasonic neuromodulation is largely motivated by future therapeutic applications. This is the only neuromodulation technology currently capable of being selective, targeted, reversible and noninvasive with millimeter-scale precision essentially across the entire human and nonhuman primate brain (Deffieux et al., 2013; Coluccia et al., 2014; Legon et al., 2014; Lee et al., 2015) and peripheral nervous system (Tsui et al., 2005; Colucci et al., 2009; Juan et al., 2014). Although ultrasonic neuromodulation phenomena could potentially result from a diverse set of different biophysical interaction modes, including temperature elevation (Lele, 1963; Colucci et al., 2009), acoustic streaming (Tyler, 2011), radiation pressure (Tyler, 2011; Prieto et al., 2013), and stable or inertial cavitation (Tyler, 2011), only the recently proposed intramembrane cavitation hypothesis (Krasovitski et al., 2011) and the related neuronal intramembrane cavitation excitation (NICE) model (Plaksin et al., 2014) have led to a detailed predictive explanation of specific empirical results for CNS stimulation by low-intensity US. In the NICE framework, US-induced oscillatory expansions and contractions of small membrane regions, the bilayer sonophores (BLSs; Krasovitski et al., 2011), cause capacitive displacement currents that indirectly lead to slow membrane charge accumulation (on a time scale of tens of milliseconds), ultimately bringing the neurons to their AP discharge threshold. The NICE model explains the basic features of acoustic cortical neuro-stimulation of pyramidal neurons, including the long durations required and the apparent higher efficiency of continuous wave (CW) versus pulse-mode (discontinuous) stimulation (Tufail et al., 2010, 2011; King et al., 2013), and its predictions have a very good agreement with the recent experimental results of detailed parameter-dependent efficacy for mouse motor cortex direct ultrasonic stimulation (King et al., 2013). In contrast, the biophysical basis of CNS neural suppression (Yoo et al., 2011a; Kim et al., 2015) remains elusive, and the ability to differentially mediate suppression simply by varying the stimulation parameters, did not receive a rigorous, quantitative, and predictive treatment. Fundamentally, it is unclear whether suppression is mediated by directly suppressing active excitatory neurons, by activating inhibitory neurons, or by depleting excitatory synaptic pools, to mention just a few possibilities.
Here, we extend the emerging intramembrane cavitation framework to examine the predicted effects of US on additional cell types and stimulation modes. We analyze the effect of pulsed mode US with parameters ranging from CW to low duty-cycle (previously shown to suppress cortical and thalamic activity; Yoo et al., 2011a; Min et al., 2011; Kim et al., 2015) on NICE-type models of multiple types of cortical neurons and related thalamic neurons. The new analysis reveals that low-threshold spiking (LTS) inhibitory cortical neurons (Pospischil et al., 2008) and major types of thalamic neurons (Destexhe et al., 1998a) are hyper-sensitive during discontinuous pulsed US stimulation. In these neurons, the presence of T-type voltage-gated calcium channels boosts charge accumulation between the short US bursts, leading to their selective activation when low duty-cycle waveforms are applied. These differential US-induced responses of different cell types are shown to explain the emergence of neuromodulation parameter ranges for obtaining network stimulation or suppression that match and clarify empirical parameter choices. Moreover, how these effective parameter ranges are modulated by baseline input into the cortical network provide further empirical predictions regarding the dependence on thalamic input, anesthetic modulation, etc.
Theoretical framework
The biomechanical-biophysical models
We simulated and analyzed the expected effect of US on five types of mammalian cortical and thalamic cell types where detailed BLS-type membrane interactions are coupled with Hodgkin-Huxley (H&H) single-compartment models of three cortical neuron types (Pospischil et al., 2008) and two central types of thalamic neurons (Destexhe et al., 1998a). The cortical neuron types include an excitatory regular spiking (RS) pyramidal neuron and two principal models for cortical inhibitory interneurons: LTS and fast spiking (FS) neurons, whereas the thalamic types included thalamocortical (TC) and thalamic reticular (RE) neurons. The various NICE model parameters were taken “as is” without retuning or post hoc adjustments, and are generally based on previously measured or estimated physical and biophysical quantities (summarized in Table 1 with the respective sources). In the model, a circular, uniform phospholipid bilayer membrane dome is bounded between immobile protein islands (Fig. 1A; Plaksin et al.,2014); structured intramembrane cavities appear to be consistent with the observed clustering of protein distributions in real cells’ membranes (Lillemeier et al., 2006, their Fig. 1, where protein-free patches of 50–100 nm diameters are evident). Electrically, the bilayer membrane has a varying capacitance, and each ion has a Nernst equilibrium potential and a time-dependent conductance (Fig. 1A, right), which generally depends on the product probabilities of multiple voltage-dependent gates (M and H gates for sodium channels; N and P gates for potassium channels and S and U for calcium channels; Pospischil et al., 2008). When US is applied, the membrane leaflets oscillate according to modified Rayleigh–Plesset bubble dynamics (Krasovitski et al., 2011), causing oscillatory variations in the membrane capacitance (Fig. 1B), which lead to an oscillatory displacement current VmdCm/dt that directly modulates the membrane potential (Plaksin et al., 2014). The responses of the nanometer-scale BLS model to US are assumed to be representative of the responses of the whole cell; US waves with sub-millihertz frequencies have wavelengths of several millimeters, orders of magnitude larger than the dimensions of CNS cortical neuron somata, so all BLS elements are subject to essentially the same acoustic effect (Plaksin et al., 2014).
Biomechanical and biophysical parameters for the simulation runs
Cortical and thalamic NICE models. A, Geometrical and biophysical representation structure of the NICE models: top view (left) of the US-induced dome-shaped BLS intermembrane cavities (light gray) in the plasma membrane bare zones (dark grey), bounded by cholesterol-rich protein islands (red areas). The equivalent electrical circuit of this biophysical complex structure (right) includes a potential (Vm), time-varying capacitance (Cm), and Hodgkin–Huxley type ionic conductances (gi) and sources (Vi). Each neuron type channels' composition is summarized in the neocortical and thalamic tables. B, Electrical dynamics during first three cycles of the model membrane exposed to US (f=0.69 MHz, 3.3 W/cm2): acoustic pressure (kPa), membrane capacitance (μF/cm2), and membrane potential (mV). C, A simplified network of RS, FS, and LTS cortical neurons. The filled black circles and open triangles are GABAA and AMPA- type synapses, respectively. The excitatory connections to the two FS and LTS inhibitory neurons are depressing and facilitating, respectively. The synaptic strength is represented by changes of the lines' thickness (logarithmically scaled) and ITh-RS and ITh-FS are the thalamic inputs.
To explore the significance of network interactions and baseline activity on acoustic neuromodulation we explored the effect of US stimuli on a spiking cortical network model with thalamic input (Fig. 1C; Hayut et al., 2011) where the RS, FS, and LTS cell types are coupled through excitatory (AMPA) and inhibitory (GABAA) synaptic connections whose strength and short-term dynamics are based on Vierling-Claassen et al. (2010).
Models' equations
The governing mathematical expression used in our study are H&H-based single compartment model equations, adapted to mammalian RS, FS, and LTS cortical neurons (Pospischil et al., 2008) and thalamic TC and RE neurons (Destexhe et al., 1998a); in addition to the H&H ionic currents, the models include US-induced displacement currents (Plaksin et al., 2014). Excitatory and inhibitory synaptic currents (Vierling-Claassen et al., 2010) and thalamic inputs to the RS and FS neurons (but not to the LTS) were added to the neocortical network simulations (Hayut et al., 2011):
where 
is the cell membrane potential; 
is the cell membrane capacitance; 
represents the capacitive displacement current induced by the US-subjected BLS dynamics, leading to changes in the average membrane capacitance (Plaksin et al., 2014; Fig. 1B); 
, 
, 
, 
, 
, 
and ILeak are the currents of the sodium, delayed-rectifier potassium, slow non-inactivating potassium, low-threshold calcium, hyperpolarization-activated mixed cationic, leak potassium and the non-voltage-dependent nonspecific ion channels, respectively [
, 
, and 
exists in all NICE models and their dynamics were taken from Pospischil et al. (2008); 
only in the LTS and TC neurons, Huguenard and McCormick (1992); 
only in the RE neuron, Destexhe et al. (1996b); 
and 
only in the TC neuron, Destexhe et al. (1996a) and 
only in the RS, FS, and LTS neurons, Pospischil et al., (2008)]; 
is the sum of FS and LTS neurons related GABAA synapses induced inhibitory currents; 
is the RS neuron related AMPA synapses induced excitatory current and 
is the thalamic DC current input to the RS and FS neurons with ITh-FS/ITh-RS
=1.4, following Hayut et al. (2011; see also Cruikshank et al., 2007), and ITh-RS
=0.17 nA, which generates baseline cortical pyramidal neuron average firing rate of about 7Hz, matching the average spontaneous firing rate of ketamine-xylazine anesthetized animals (Destexhe and Paré., 1999).
The network model structure was adopted from Hayut et al. (2011), with synaptic currents modeled as in Vierling-Claassen et al. (2010):
where i and j indexes can be RS, FS, and LTS related to the presynaptic and postsynaptic neurons, respectively; 
relates to the short term synaptic plasticity for the excitatory connections to both the inhibitory neurons as described by Vierling-Claassen et al. (2010) and modeled as by Varela et al. (1997); 
relates to synaptic open probability modeled as bi-exponential function (Kleppe and Robinson, 1999) and 
and 
represents the AMPA or GABAA synaptic currents reversal potentials and strengths, respectively.
All the models' biomechanical, biophysical, and synaptic parameters and their respective sources can be found in Table 1.
LTS and TC neurons T-type calcium channels current dynamics
The LTS inhibitory interneuron's T-type calcium channels current was modeled as by Pospischil et al. (2008):
where VCa is the Nernst potential of the calcium ions; s and u, designate the open probabilities of the activation and inactivation S-type and U-type gates, respectively, and 
, 
, and 
are the open probability, voltage-dependent steady state open probability and voltage-dependent time constant of the S-type or U-type gates, respectively.
The gates' dynamics are based on Huguenard and McCormick (1992):
where the time constants were adapted from of 24°C to 36°C, using a Q10 factor of 3 (Pospischil et al., 2008) and Vx=−7 mV (Vx=−2 mV for Fig. 7) is a voltage-dependence uniform shift factor defined by Pospischil et al. (2008). The TC neuron's T-type calcium channels gates’ dynamics are also described by these expressions, with Vx=0 mV (Huguenard and McCormick, 1992; Destexhe et al., 1998a).
Models’ implementation
The model’s set of equations was numerically solved in MATLAB (using the function ODE113). The time difference between the calculated points was set to 
µs (where f is the US frequency in megahertz). The NICE models, including the network model, were solved in the same manner as explained in our previous study (Plaksin et al., 2014).
The US intensities (I) reported in this study were calculated in the form of spatial peak-pulse average intensities for propagating planar US waves (Hendee and Ritenour, 2002):
where PA
is the pressure amplitude, 
is the surrounding medium density, and c is the speed of sound in the medium.
Results
Prediction I: the responses of cortical neurons to ultrasound are strongly cell-type and waveform-dependent
We first studied the fundamental response of the three NICE-type models of cortical neurons to CW and low duty-cycle (5%) US stimulation, as used by Yoo et al. (2011a) to respectively excite or suppress cortical activity (Fig. 2; 0.69 MHz, acoustic pressure amplitude 320 kPa, intensity 3.3 W/cm2). During CW excitation, the US-frequency-driven oscillations of the NICE model neurons' intramembrane space lead to strongly hyperpolarized oscillations of the membrane potential, and the generation of a delayed AP train (after 10–25 ms; Fig. 2A); this excitation process occurs as a result of a charge accumulation mechanism described previously for the RS-NICE model (Plaksin et al., 2014). In contrast to this fairly uniform response of the three neuron types to CW excitation, their responses to low duty-cycle excitation is highly divergent (Fig. 2B), and only the LTS-type neuron is effectively excited by this excitation mode and tonically fires a volley of APs. Pulsed excitation with a low 5% duty-cycle and with varying US durations and frequencies, has excitation thresholds for LTS neuron that are upward of three orders of magnitude lower than those for FS and RS neurons; this ratio decreased rapidly as the duty-cycle increased, down to comparable thresholds at 50% duty-cycle (Fig. 2C–E). These excitation thresholds for the three cell types varied essentially independently of the excitation duration, and had only a weak dependence on US frequency, varying by <10% between 0.2–1 MHz (data not shown).
Effect of continuous and pulsed US stimuli on the different cortical NICE-neurons (f=0.69 MHz). A, B, Effect of US stimulus (3.3 W/cm2, indicated by bars) on membrane potential and charge (top), sodium and potassium channels kinetics (middle), and on LTS neuron T-type calcium channels kinetics (bottom). Fifty millisecond continuous stimulus, effectively stimulates all neuron types (A), whereas a 300-ms-long pulsed stimulus (pulse repetition frequency (PRF) 100 Hz and duty-cycle 5%) causes only the LTS neuron to tonically fire a volley of APs (B). This selective LTS excitation is mediated through the elevation of the T-type calcium channels' S-gates open probability during the US off times (right), which elevates these channels' conductance and consequently amplifies the charge accumulation process that occurs during US's-on periods. C–E, Threshold intensity versus duration required to generate a single AP using constant duty-cycle (PRF, 100 Hz). The excitation thresholds for the RS and FS neurons at 5% duty-cycle are >3.5 orders of magnitude higher than for the LTS neuron (E), decreasing rapidly to ∼2× at 50% duty-cycle (C).
What is the detailed biophysical basis underlying this selective AP generation by low duty-cycle US in the cortical LTS neurons? During US-on periods, the induced displacement currents cause a rapidly fluctuating hyperpolarized membrane potentials (Fig. 1B) that quickly suppresses the voltage-gated channels, whereas the non-voltage-dependent (leak) ion channels accumulate charge (Fig. 3A, left). Examining this process in high temporal detail (Fig. 3B), the asymmetric hyperpolarization is seen driving negative leak currents and charge accumulation while acting to close voltage-gated activation gates and to open the inactivation gates. The main difference between the cell types occurs during the US-breaks: the LTS neuron's T-type channels are the only VGCs that continue the charge accumulation process during these extended durations (Fig. 3A, left, C, top), due to intrinsic features that enhance their conductance during the US-off periods. Specifically, the contrast between the relatively slow recovery time of U-type gates (τ≈15 ms) versus fast recovery time (τ≈2 ms) of the S-type gates augments the channels conductance by allowing activation to recover relatively quickly (Fig. 3C, left) and net charge accumulation (eg, compare to the Na+ currents, which are negligible prior to AP initiation; Fig. 3C, right).
Detailed US response of LTS neurons (f=0.69 MHz). A, The contribution of each channel type to the accumulated membrane charge during 10 ms of CW versus a short-pulsed US stimulus (5% duty-cycle, PRF=100Hz): leak channels have the biggest contribution during the US-on period, whereas the T-type calcium channels dominate the US-off period. B, Leak and calcium channels' dynamical response to the first few US cycles (1.3 W/cm2); the hyperpolarized phase drives negative leak currents that insert positive charge into the cell, while rapidly suppressing the calcium conductance due to the changes in S- and U-type gates open probability p(t), through dynamical perturbations of the steady state probability (p∞), and the gates' time constants (τ). C, T-type calcium versus sodium channels' dynamical responses during sparse stimulation (5% duty-cycle, 1.3 W/cm2); the comparison highlights the dramatic changes during the US breaks in the calcium currents, open probability p(t) and the steady-state open probability (p∞) of the S- and U-type gates, whereas the Na+ gates are mostly dormant prior to action potential initiation (arrow). D, The pulsed US excitation thresholds of native RS and FS neurons versus following the chimeric addition of T-type calcium channels (RS+ and FS+).
Finally, to test whether the T-type channels’ impact on low duty-cycle excitation is independent of the other specific biophysical properties of the LTS neuron model (eg, different AP thresholds parameters), we added these channels to RS and FS neuron models. This chimeric manipulation was found to create similar orders-of-magnitude differences in the neural thresholds and sensitivity to low duty-cycle excitation (Fig. 3D) through the same mechanisms.
Because the excitation thresholds are nearly independent of US frequency and duration, we next analyzed the joint dependence of the LTS and RS neurons on both duty-cycle and intensity, at a constant frequency and duration (0.69 MHz and 500 ms; Fig. 4). This analysis lead to a two-dimensional phase diagram that naturally separates into two activation domains where either LTS neurons are activated alone (at low duty-cycles) or both neuron types are jointly activated (at high duty-cycles). Our basic “first order” expectation (examined in detail in the next section) is that selectively exciting LTS inhibitory neurons will lead to network suppression, whereas jointly activating all neuron types with generally similar rates (where >75% are excitatory; Markram et al., 2004) will lead to net network stimulation. Next, we examined how these model-based predictions fit empirical data by marking on the parametric phase diagram the empirically-reported stimulation parameter ranges used in nearly all cortical ultrasonic neuromodulation studies (color-coded according to whether the reported effect was excitatory or suppressive). Strikingly, the result of this comparison (Fig. 4) demonstrates an excellent qualitative match between the model-based phase diagram and essentially the entire range of empirical stimulation parameters used for ultrasonic excitation and suppression. The diagram also shows that a two-order-of-magnitude variation in pulse repetition frequency (PRF changes from 10 to 1000 Hz) has a relatively minor influence on the response thresholds; anecdotally, the lowest duty-cycle empirical activation conditions (leftmost borders for King et al., 2013 results, Fig. 4) correspond to an intermediate PRF of about 400 Hz, and are seen to be in a fairly good agreement with the expected thresholds.
Phase plane diagram of single-neuron responses to varying US stimulation duty-cycle and intensity versus experimental cortical neuromodulation parameters. The phase diagram boundaries denote threshold intensities for US-mediated responses (frequency 0.69 MHz, duration 500 ms) from excitatory RS neurons (green dashed lines indicating 10 Hz and 1 kHz PRFs) and inhibitory LTS interneurons (red dashed lines, changes only slightly for different PRFs, not shown). These boundaries separate the phase diagram into regions where either the inhibitory LTS neurons are activated alone (red, “suppression zone”) or the RS and the LTS neurons are jointly activated leading to net network stimulation (green, “activation zone”). The superposed bars indicate the experimental parameter ranges used in seven published cortical ultrasonic neuromodulation studies, color-coded according to the mediated responses: Ref. 1 (King et al., 2013; bars with diagonal lines), Ref. 2 (Yoo et al., 2011a), Ref. 3 (Kim et al., 2015), Ref. 4 (Kim et al., 2012), Ref. 5 (Kim et al., 2014), Ref. 6 (King et al., 2014), and Ref. 7 (Tufail et al., 2011). The excitation parameters reported for King et al. (2013) were those that caused stimulation success rates significantly higher than their noise floor (∼20%), with low-frequency CW intensities corrected for the expected formation of standing waves (Plaksin et al., 2014).
Prediction II: effect of ultrasound on cortical network activity reflects cell-type selective responses and is modulated by thalamic inputs’ strength
To examine how different US stimulation parameters affect cortical activity at the network level, we studied the responses of a simplified but physiologically plausible network population model to US pulses with varying parameters. The network contains the three cortical neuron population types studied above, where the RS and FS neurons receive thalamic inputs and all three NICE-neuron types are subjected to a 1-s-long US stimulus (Fig. 5, black bars). The model was adapted from Hayut et al. (2011) to studying the effect of US by moving from a rate-based network model to a H&H-based model connected by dynamic GABAA and AMPA synapses, (Fig. 1C; see the Theoretical framework section for more information). When US stimulation is applied with a low 5% duty-cycle and intensity (0.1 W/cm2 or 56 kPa pressure amplitude, 100Hz PRF), no significant network response is observed (Fig. 5A), but at a higher intensity of 3.3 W/cm2 [320 kPa, used by Yoo et al. (2011a) and Kim et al. (2015) to obtain visual cortex suppression], the baseline activity of RS and FS populations was inhibited by the US-excitation of the LTS population, which responded at an elevated rate of ∼40 Hz (Fig. 5B). In contrast, using waveform parameters that were used to elicit motor cortical excitation: 50% duty-cycle and 10 Hz PRF, but having the same intensity (Yoo et al., 2011a), the RS and FS responses are dramatically increased by the extended US on times, suppressing the LTS responses and thus leading to overall excitation (Fig. 5C). These predictions from the network analysis are thus seen to be consistent with empirical findings; overall, the analysis predicts that the network's two-dimensional parametric phase diagram (Fig. 5D) naturally separates into a suppression domain, where LTS neurons are activated (at low duty-cycles), an activation domain, where the RS are activated, and a transition zone, where neither effect dominates. Interestingly, examination of experimentally tested neuromodulation parameters (Yoo et al., 2011a; Kim et al., 2015; Lee et al., 2015) with similar US intensities but varying duty-cycles or with the same 50% duty-cycle but varying US intensities, respectively, are found to have an excellent match with these predictions (Fig. 5D, points b–f and the green vertical bar). In these data either predominant suppression or stimulation are observed (points b, c, and the green bar; Yoo et al., 2011a; Lee et al., 2015) or a “transition” where neither clearly dominates (points d–f; Kim et al., 2015; Lee et al., 2015). Importantly, the zone borders are modulated by the strength of external (thalamic) input (Fig. 5D): increasing this input to the RS and FS cells lowers the threshold ultrasonic intensities required for RS stimulation and also shrinks the red suppression zone (Fig. 5D, inset). This shrinking effect is mediated indirectly: the stronger thalamic input drive synergistically compounds the US-induced excitation of the FS activity which acts to counteract the LTS's suppressive activity; this effect explains both the upward trend in the threshold drive required to excite LTS neurons and the (weaker) downward trend in the threshold required to overcome their effect (Fig. 5D, inset).
Simplified cortical NICE-network responses to different US waveforms and intensities. The US stimuli (US frequency and duration: 0.69 MHz and 1 s) are indicated by black bars (A–C). A, For a stimulus duty-cycle of 5% and 0.1 W/cm2 intensity (PRF, 100 Hz) no significant response to US is observed. B, Increasing the intensity to 3.3 W/cm2 causes FS and RS activity suppression due to strong LTS activation (∼40 Hz). C, Increasing the duty-cycle to 50% (PRF, 10 Hz) leads to high frequency activation of the RS and FS neurons, unsuppressed by the weaker LTS firing (only at the beginning of each US pulse). D, Phase plane diagram for the network responses to US with varying duty-cycle and intensity (PRF, 100 Hz). Marks a–c indicate the conditions of the respective simulations (matching the experimental observations of Yoo et al. (2011a) and marks d, e, indicate parameters from Kim et al. (2015) where the experimental responses were no longer suppressive. The vertical green bar represent human primary somatosensory cortex stimulation parameters used to evoke tactile sensations (Lee et al., 2015); f marks the only case where no response was observed. The green and red arrows and the inset depict the effect of increased thalamic input on the activation and suppression thresholds.
Prediction III: thalamic neurons respond to low duty-cycle ultrasound, with exquisite sensitivity to waveform parameters
Finally, we examined whether the characteristics of ultrasonic neuromodulation of thalamic circuits can be predicted by modeling the responses to ultrasound of two major thalamic neural cell types: thalamic RE and TC neurons, to the experimentally tested US pulse stimulation parameters (100 Hz pulses, 5% duty-cycle; Min et al., 2011; Yoo et al., 2011b; Yang et al., 2012). These neuron types also contain T-type calcium channels currents: IT and ITs in the TC and RE neurons, respectively (Destexhe et al., 1998a; see Theoretical framework section). Both TC and RE neuron models were found to respond to low duty-cycle US stimulus (Fig. 6) with a common response mechanism that is related to the calcium channel dynamics (Fig. 6A,B, bottom), similarly to the one described above for the LTS cortical neurons (Figs. 2B, 3). However, whereas the TC neuron fires a tonic volley of APs almost immediately after the onset of the simulated ultrasonic stimulus, the RE neuron fires only one volley of APs before stopping. The TC neurons' calcium channels show stronger responses to both the hyperpolarized oscillations and the depolarized potentials, leading to stronger calcium currents in the US-off periods and therefore to a higher sensitivity to the US-pulse (Fig. 6, compare A, B, bottom left). In addition, the T-type calcium channels' S gates deactivation in the after-hyperpolarized state (following an AP) is relatively slow, and charge accumulation during this period repeatedly brings them to a state that is susceptible to re-excitation via brief bursts of applied US (Fig. 6A, bottom right). To explore the impact of these differential responses across a range of waveform parameters, we studied the models’ responses to waveforms with increasing duty-cycle, finding that even slight increases in the pulses' duty-cycle (from 5% to 7%; Fig. 6B,C) cause a dramatic shift in the US responses of the RE neuron type, eliminating its rapid adaptation to the US stimulus and pushing the neuron to a tonic firing regime; in contrast, the TC responses are not modified by these subtle changes (Fig. 6D). We also note that the T-type calcium channels in both of these thalamic neuron types have a much higher maximal conductance than in the cortical LTS neuron model (Table 1), leading to a generally higher level of excitability.
The response of thalamic NICE-TC and NICE-RE models to low duty-cycle US stimulation waveforms. The US stimuli (US intensity: 5.2 W/cm2; US frequency: 0.69 MHz; PRF, 100 Hz) are indicated by black bars (A–C). A, B, For a 1.5 s, 5% duty-cycle US stimulus, the TC cell fires a tonic 100 Hz volley of APs, whereas the RE cell fires only one volley and stops. Bottom, The currents' profiles of the segments marked in the top, where IU is the sum of Ih, 
, and ILeak currents (see complete channel composition in the Theoretical framework section). C, Increasing the duty-cycle to 6% and 7% brings the RE neurons to fire periodical volleys and a constant volley of APs after two braked volleys, respectively. D, The relation between the TC and RE neurons' spike rates and the US stimulation duty-cycle, calculated for the last 0.5 s period of the 1.5-s-long US stimulation.
Discussion
We theoretically studied the effect of this new mode of biomechanical-biophysical interaction between US and neurons with the goal of understanding the observed suppression of cortical neurons by low duty-cycle US. Such an understanding may guide the development of both therapeutic applications, such as seizure suppression (Min et al., 2011) and a myriad of research applications. The fundamental charge accumulation mechanism resulting from the imbalance of ionic currents during hyperpolarizing oscillations, which we have previously shown to underlie the generation of APs in cortical pyramidal (RS) NICE neurons (Plaksin et al., 2014), was found here to behave fundamentally differently in LTS neurons under low duty-cycle pulsed excitation. In these neurons, the overall suppression of all voltage-gated channels during the US-on periods (Fig. 2B) is accompanied by partial opening of the T-type calcium channels' U-gates; coupled with US-induced charge accumulation, this results in opening the calcium channels' S-gates, whose increased conductance effectively boosts membrane depolarization during the prolonged US-off periods (Fig. 3A, left, C, top) and leads to effective, cell-type selective excitation (Fig. 2B, bottom). Analysis of chimeric models of RS and FS neurons where the T-type channels were added further demonstrates their crucial role in this selective excitation (Fig. 3D).
The model-based predictions were found to provide a post hoc justification and understanding for the results of a significant body of experimental studies (Fig. 4), where parameter regimes naturally separated into suppression versus activation domains according to the pulse train’s duty-cycle, essentially independently of other stimulation parameters. This separation in terms of cell-type-selective interactions is, to the best of our knowledge, completely novel and was also found to carry over to the network level (Fig. 5). Here, additional interesting issues were observed including modulation of the suppression and stimulation thresholds by the thalamic input level (Fig. 5D) and the emergence of a transition zone between the suppression and excitation domains, where joint stimulation of multiple cell types results in a balanced state where none of the effects dominate. The network model's predictions have excellent agreement with very recently reported experimental results in rats (Kim et al., 2015), as well as in human primary somatosensory cortex (Lee et al., 2015), where artificial tactile sensations were reportedly evoked in all but one subject (Fig. 5D). Additionally, we note that these network-level observations could also explain several other effects that were associated with US stimulation, for example, it predicts a dependence of the US excitation parameters on the doses and type of anesthesia used, which is known to directly affect the level of thalamic activity (Alkire and Miller, 2005), moreover, it predicts that an extended ultrasonic suppression regimen will strongly drive LTS neurons and can lead to post-tetanic potentiation of their synapses with plasticity effects that can last for minutes after the stimulation ends (Storozhuk et al., 2005; Yoo et al., 2011a). Finally, we note that the network analysis predicts the existence of an optimum duty-cycle where the activation threshold is minimal (∼70%; Fig. 5D); this convex behavior, first predicted by Plaksin et al. (2014), reflects a tradeoff in RS and FS neurons between charge accumulation during US-on and the ability to discharge during US-off periods, and was recently experimentally validated by Kim et al. (2014). It is important to note that the model apparently has a very strong predictive power despite the multiple simplifying assumptions that our theoretical approach invokes regarding the biomechanics and biophysics of membranes, and ignoring the morphological complexity, which exists at both the single neuron and the network levels. This robustness could potentially indicate that the presented model already captures the essential features of the underlying biophysics, a conclusion supported by early analyses using NEURON models of detailed cortical NICE-networks and detailed analyses showing negligible accessory effects from radiation pressure (Prieto et al., 2013; M. Plaksin, S. Shoham, and E. Kimmel, unpublished observations). Moreover, inspecting the responses of NICE models with different sets of parameters adapted from cortical neuron models for ferrets and cats (Fig. 7), shows an insensitivity of the basic model-based US-induced behavior, and suggests that ultrasonic neuromodulation and its underlying mechanisms may be fairly universal among different mammalian species. In addition, NICE US-induced responses are robust to partial sonophore membrane coverage and the associated reduced membrane potential oscillation range (Fig. 8), nor do they strictly depend on the BLS biomechanics: the same qualitative results can be obtained by an arbitrary source of sinusoidal membrane capacitance variations (Fig. 9).
Effect on cortical NICE-neuron models from different mammalian species (Pospischil et al., 2008) of continuous and 5% duty-cycle pulsed US stimuli (US intensity: 3.3 W/cm2; US frequency: 0.69 MHz; PRF, 100 Hz, indicated by bars). US stimulus effects on membrane potential, charge and channels kinetics for cortical neurons of two different mammals (RS and FS, ferret visual cortex; LTS, cat association cortex). The panel organization and responses were similar to those described in Figure 2 and are explained by the very same underlying mechanisms. For continuous stimuli (A; 100 ms duration) there isn't a major difference between the responses of the different neuron types, except for a delay in the LTS neuron firing due to low leaky channels' conductances that cause slower charge accumulation. For pulsed stimuli (B; 1500 ms duration), only the LTS neuron responded.
Effect of partial sonophore membrane area coverage during continuous and 5% duty-cycle pulsed US stimuli (US intensity: 3.3 W/cm2; US frequency: 0.69 MHz; PRF, 100 Hz, indicated by bars) on cortical RS (A) and LTS-NICE (B) neuron models, respectively. Partial coverage (here 75%) reduces the membrane potential oscillations down to a narrower range (>−150mV). Although the potential oscillations were more limited, the neurons' response to continuous and pulsed stimulation is still evident. Membrane capacitance was calculated as a weighted mean of the resting and dynamic capacitances: 
, where fs is the active area fraction.
Effect of purely sinusoidal capacitive drive on cortical RS and LTS-neuron models in continuous (A) and 5% duty-cycle (B) stimulation modes (
, CAmp
≈0.8 μF/cm2, f=0.69 MHz; PRF, 100 Hz, indicated by bars). 
is the resting membrane capacitance. Although the sinusoidal and the intramembrane cavitation theory-based capacitance variations are fundamentally different, the basic qualitative neural responses remain the same. The CAmp was determined when 80% decline in the membrane capacitance (Fig. 1B; f=0.69 MHz and intensity 3.3 W/cm2) was taken into account.
In comparison to the relatively smooth parameter dependence observed in the cortical simulations, it appears that key thalamic neurons display an even stronger cell-type-selective dependence on stimulation parameters. Interestingly, activity volleys in thalamic RE neurons play a key role in the formation of certain types of epileptic seizures (Slaght et al., 2002; Lee et al., 2014). The observed super-sensitivity of the NICE-RE model to the pulsed duty-cycle parameter could thus potentially underlie a predicted transition from suppression to seizures during very mild corresponding increases in the duty-cycle of US applied to the thalamic structures (from 5% to 7%), and perhaps explain why only 5% has been experimentally reported (Min et al., 2011; Yoo et al., 2011b; Yang et al., 2012).
In summary, we have demonstrated that the NICE framework, based on analyzing the expected effects of intramembrane cavitation on multiple cell types using well characterized, Hodgkin-Huxley type canonical models of membrane biophysics, not only agrees with, but also sheds new light on a large array US neuromodulation experimental results. Importantly, our results show that waveform modulation can potentially be used to selectively drive cell-type-selective responses leading to suppression. Achieving cell-type-selective stimulation is considered one of the crowning achievements of optogenetics, and it is thus fascinating to consider that this capability can also be translated and further refined to this noninvasive neuromodulation modality and potentially easily applied even to human subjects. While this study explored several major neural cell types of the neocortex and the thalamus, future studies could further explore the responses of other neuron types in additional brain regions, to achieve a more complete selective ultrasonic stimulation ability. In fact, the T-type voltage gated calcium channels are expressed not only in the neocortex and the thalamus, but also in hippocampal pyramidal neurons (Nikonenko et al., 2005), putatively explaining the ability to stimulate the hippocampus neurons by low duty-cycles and very low US intensities (Tyler et al., 2008; Tufail et al., 2010). Given the wide design space for acoustic excitation waveforms and strategies, these results highlight the potential of this model-based approach to optimized design for various brain regions and target applications.
Footnotes
↵1 The authors report no conflict of interest.
↵3 This work was supported by the MAGNET program of the Israeli OCS at the ministry of economy, the Russell Berrie Nanotechnology Institute, and the Israel Science Foundation Grant 1725/13. We thank Steve Krupa, Roman Shusterman, Omer Naor, Amit Livneh, and Josh Nehrer for their comments on the paper.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
