Intramembrane Cavitation as a Predictive Bio-Piezoelectric Mechanism for Ultrasonic Brain Stimulation

Low-intensity ultrasonic waves can remotely and nondestructively excite central nervous system (CNS) neurons. While diverse applications for this effect are already emerging, the biophysical transduction mechanism underlying this excitation remains unclear. Recently, we suggested that ultrasound-induced intramembrane cavitation within the bilayer membrane could underlie the biomechanics of a range of observed acoustic bioeffects. In this paper, we show that, in CNS neurons, ultrasound-induced cavitation of these nanometric bilayer sonophores can induce a complex mechanoelectrical interplay leading to excitation, primarily through the effect of currents induced by membrane capacitance changes. Our model explains the basic features of CNS acoustostimulation and predicts how the experimentally observed efficacy of mouse motor cortical ultrasonic stimulation depends on stimulation parameters. These results support the hypothesis that neuronal intramembrane piezoelectricity underlies ultrasound-induced neurostimulation, and suggest that other interactions between the nervous system and pressure waves or perturbations could be explained by this new mode of biological piezoelectric transduction.


INTRODUCTION
Not only is ultrasound (US) widely used for imaging [1]; its interaction with biological tissues is known to induce a wide variety of nonthermal effects ranging from hemorrhage and necrosis [2] to more delicate manipulations of cells and their membranes such as permeability enhancement [3], angiogenesis induction [4][5][6], and increased gene transfection [7]. In particular, both classical and recent studies have demonstrated that US can interact with the physiology of excitable tissues, inducing the generation of action potentials (APs) [8][9][10][11][12][13][14][15][16][17][18], suppression of nerve conduction [19][20][21], as well as more subtle changes in excitability [22][23][24]. While the suppression of nerve conduction is putatively dependent on temperature elevation [20,21], the biophysical basis of neural stimulation is not understood and has not received a rigorous, quantitative, and predictive treatment. Recently, we introduced the bilayer sonophore (BLS) model as a unifying hypothesis for the underlying mechanism of multiple bioacoustic interactions, wherein US preferentially induces intramembrane cavitation (bubble formation) in the intramembrane space between the two lipid leaflets of the cell's membrane [25]. In the BLS model, the negative pressure phase of the US wave pulls the two leaflets away while the positive pressure pushes the leaflets towards each other; dissolved gas accumulates in the hydrophobic zone, creating pockets of gas that expand and contract periodically. BLS formation is predicted to induce various alterations in cells (bioeffects) including the initiation of cellular mechano-transduction processes, the induction of membrane pore formation, and permeability changes. These bioeffects, which are naturally considered as associated with mechanical loading, were shown to systematically and predictably intensify as the US pressure amplitude increases or the frequency decreases, in softer tissues or close to free surfaces, and in the presence of microbubbles [25].
In an attempt to understand the mechanisms underlying the effect of US on excitable tissues, we analyze a neuronal BLS (NBLS) framework where the biomechanics of intramembrane cavitation is coupled to membrane bioelectrical mechanisms in a complex interplay. This mechanoelectrical coupling is shown to induce displacement currents that excite action potentials through an indirect mechanism, whose features explain the requirement for long ultrasonic stimulation pulses [12, 16,18] and predict the experimentally observed efficacy of ultrasonic stimulation in mouse motor cortex [18].

II. MODEL & EQUATIONS
The purely mechanical BLS model has been modified to account for the dynamics of membrane charge polarization, capacitance, and voltage-sensitive ion channels in a CNS neuron [NBLS model, Fig. 1(a)]. The NBLS model combines the modified BLS model and the Hodgkin-Huxley (H&H) model adapted for a regular spiking rat cortical pyramidal neuron [26]. Both BLS and H&H model parameters are taken "as is," without re-tuning or post-hoc adjustments, and are based on known or measured physical and biophysical quantities or ranges, wherever attainable (summarized in where G Na , G K , G M and G Leak are the conductances of the sodium, delayed-rectifier potassium, slow noninactivating potassium, and leak channels, respectively, and all voltage-dependent parameters are defined in Ref. [26]. The time-dependent C m is associated with the geometrical shape of the deformed leaflets and can be determined by solving a modified dynamics force (pressure) balance equation that is based on the Rayleigh-Plesset (RP) equation for bubble dynamics [25,29]. The modified BLS equation, is the original dynamics equation [25] with a new electric equivalent pressure term, that are captured by an equivalent "molecular pressure," where a, Z, z(r) and R(Z) are defined in Fig. 1; l r is the membrane surrounding water-like medium density; ε 0 and ε r are the dielectric and relative dielectric constants, respectively, of vacuum and of the transmission medium between the two leaflets (e r chosen here to be equal to 1), and D , * D are the initial gaps between the leaflets if there are or there are no charges on the membrane, respectively. Equation (3) accounts for the effective pressure due to the attraction forces between the electric ion charges on the membrane leaflets. The other pressure terms in Eq. (2) account for the surface tension [25] in the leaflets ( S P ) and for the intramembrane gas (air; N 2 +O 2 ) pressure (P in ). The rest pressure is 0 P and the driving pressure is the US pressure sin( ) The last pressure term on the rhs of Eq. (2) accounts for the viscous loss [25].
The value of P in is determined by the following gas balance equation: ( ) where a D is the diffusion coefficient of the air in the water; and ii) the intra-membrane gaseous compartment, where a k is the Henry coefficient, a n is the mole content of ideal gas, g R is the gas constant, Tem is the temperature, and V a is the intramembrane cavity volume, which is expressed by The gas transport takes place across a boundary layer with thickness x near the leaflet. For simplicity, we assume that x =0.5 nm and as a result the gas reaches almost immediate equilibrium oscillations between the dissolved and the gaseous compartments. Such an assumption saves much time, avoiding the complex calculations of the air concentration field in the surrounding medium, and is justified by the relative short times that the intramembrane space requires to reach stable oscillations [25]. We also assume that [ As the two leaflets of the BLS separate, deform, and curve periodically with the US pressure, the membrane capacitance in Eqs. (1) and (3) is approximately derived by which is based on a parallel-plate capacitor expression per unit area and where the 0 m C is cell membrane capacity at rest.

A. Fundamental response to ultrasound and AP generation
We first studied the NBLS model's fundamental response to US stimulation. When continuous-wave (CW) US stimulates the model neuron, the intramembrane space inflates and deflates at the US frequency [25], and the NBLS membrane potential oscillates strongly between -280 and -60 mV ( does not occur in this particular neuron model [26].)

B. Dependence on ultrasound parameters
Next, we examined the influence of US frequency and duration on the threshold CW US intensity (and energy) required to generate an action potential and the number of

C. Biophysical model predicts in vivo results
Finally, to validate the NBLS model, we compared its predictions to the results of a recent in-vivo study [18] where a wide range of US parameters were used to stimulate mouse primary motor cortex while resulting front limb muscle EMG signals were measured. To compare model predictions with experimental measurements we used Buckingham-Pi dimensional analysis [31] to relate NBLS' measuresthe number of AP spikes (N), the response latency (L), and the overall duration of APs (D), to the experimental success rate (R sr ). Two dimensionless variables, N and the response "effectiveness" D/(D+L), are associated with success rate in our model, but the latter quantity also appears to be solely a function of N (Fig. S1 [27]). Therefore, R sr depends only on N, a dependence that is well approximated by a sigmoidal-shaped logistic function (see the Appendix and Fig. S2 Fig. 4(b)]. When the pulse duration or the amplitude increases, the success rate increases gradually with power and duration until it reaches a plateau of 100 % success rates.
A final simulation study examined excitation by 0.5 MHz US pulse trains for different intensities and duty cycle values. (Pulse-mode US is commonly used in applications in which it is desirable to avoid heating the tissue, including neural stimulation.) We have found that the AP excitation mechanism for pulsed excitation is generally the same as for the CW mode (Fig. 2). Success probability for CW excitation is predictably higher than for 30 % duty-cycle pulsed excitation, a trend which is consistent with the experimental results [18]. Interestingly, however, a certain intermediate duty cycle value (about 70 %) is predicted to have a maximal success rate, higher than for CW excitation [ Fig. 4(c)].

IV. DISCUSSION
We studied a new biomechanical-biophysical mode of interaction between US and neurons with the goal of investigating the non-destructive manipulation of excitable tissues by US [11]. Such an understanding can guide the development of future therapeutic applications of the only technology currently capable of targeted noninvasive brain stimulation. Exposing the NBLS element to simulated US results in rapid oscillatory hyperpolarizing currents and can lead to AP generation through a charge accumulation mechanism that results from the imbalance of ionic currents.
This indirect, energetically inefficient excitation mechanism explains the characteristically very long pulses required for acoustostimulation -tens to hundreds of milliseconds [12, 16,18], compared to submillisecond pulses typically used for electrical stimulation or even for photothermal stimulation [32,33] (recently shown [34]  result lends some support to our underlying assumptions, simplifications, and the natural choice of (unadjusted) parameters: In particular, it putatively supports using the response of a nanometric NBLS as predictive of the compound average behavior of the whole CNS neuron and of a whole neuronal population exposed to US (a realistic assumption given the millimeter-scale beam dimensions and wavelengths).
While the assumption of structured intra-membrane cavities appears perhaps unrealistically simplistic, observed protein distributions in real cells' membranes are in fact clustered and somewhat reminiscent of these idealized circular distributions (see, e.g., Fig. 1 in Ref. [35], where protein-free patches of 50-100nm diameters are evident). The current NBLS model captures only excitatory regular-spiking pyramidal cortical neurons. These are not only the most common neocortical neuron type [26], but they also account for the efferent output onto the motor periphery [36]. It is also worth noting that intramembrane cavitation is not the only possible underlying mechanism for US-induced excitation; for example, the flexoelectricity model is a popular choice for explaining biological mechanoelectrical conversion [37] and is also based on varying the curvature of the bilayer membrane. In Petrov's model, the two membrane leaflets deform together as one sheet with roughly the same radius of curvature, while in our model, the distance between the leaflets varies as they expand (and collapse) away from each other [ Fig. 1(b)]. Our model is thus much more MmPiezo2 [39], presynaptic calcium channels [17,40], or postsynaptic NMDA and AMPA glutamate receptors [41].
The NBLS framework not only agrees with but also sheds new light on experimental results. While the dependence of success rate on frequency [ Fig. 4(a)] can be interpreted as strong [18], Heimburg and Jackson [42] propose an alternative electromechanical basis for action potential propagation that uses these capacitive currents and which is consistent with the action-potential-induced transient swelling observed by Iwasa, Tasaki, and Gibbons [43]. Similar interactions to the ones we have studied may underlie and/or shed new light on other types of biological mechanoelectric coupling, from the sensation of sound vibration in the auditory system (where sub-Pascal-scale thresholdpressure amplitudes [44] could rely on hypersensitized NBLS-type mechanisms) and up to the impact of mega-Pascal pressure shocks that cause concussions and traumatic brain injury [45][46][47]. (1) and (8) where β 0 and β 1 can be estimated from the success rate (R sr ) versus intensity (Table I in Ref. [18]), and the number of AP spikes (N) versus intensity predicted by the NBLS model at US frequency of 0.5 MHz. Two data points were used (R sr in percentage, N) (28.7, 34) and (52.9, 45). The calibration curve can be seen in Fig. S2 [27].

ACKNOWLEDGMENTS
(A1) The intensity (I) in a propagating US CW [48] was calculated by: where P A is the pressure amplitude, l r is the surrounding medium density, and c is the speed of sound in the medium. The spatial-peak-pulse average intensities in the pulsed-mode simulations were calculated using Eq. (A2) as well.  The maximal conductance of the Na + channels.

Mechanism for Ultrasonic Brain Stimulation
The maximal conductance of the delayed-rectifier K + channels.
The maximal conductance of the slow non-inactivating K + channels.
The maximal conductance of the non-voltage-dependent, non-specific ions channels.
The Nernst potential of the Na + Na V mV 50 The Nernst potential of the K + The Nernst potential of the nonvoltage-dependent, non-specific ion channels.
Spike threshold adjustment parameter.
The decay time constant for adaptation at slow non-inactivating K + channels. max t ms 608 The cell membrane capacity at rest. . † The patch's radius a was derived from the average distance between neighboring proteins (53 nm) in native oocytes [S8]. In a rectangular grid geometry, such an average inter-protein distance corresponds to a radius of 31±1 nm (half of the unit cell's diagonal ± a derived error estimate, of which we chose the maximal, least stiff value 32nm).