Can ionic concentration changes due to mechanical deformation be responsible for the neurostimulation caused by focused ultrasound? A simulation study

We conducted a set of computational experiments to investigate the feasibility of the ultrasound-induced ionic concentration to cause the generation of AP in neuronal axons (Figure 1). We used HH (Hodgkin and Huxley 1952) and C-fibre (Tigerholm et al 2014) models of an axon, which were electrically coupled to the extracellular space through 3D FEM. We established a simple analytical relationship between ultrasound-induced geometric mechanical strain (Zheng and Wei 2011, Pods et al 2013) and the boundary conditions for purely electric axonal models—the second derivative of the potential with respect to a coordinate along the main axis of the axon (activating function). Finally, we conducted a series of computational experiments varying different parameters of ultrasound stimulation (amplitude, frequency, ultrasound stop phase, pulse width, and inter-stimulus interval), and observed the generation of AP propagating along the axon. We summarized and qualitatively assessed the results and made conclusions about the size of the effect, the feasibility of it being the prime mechanism for AP generation, and how well the observed qualitative effects could match existing experimental data.

Zoom In Zoom Out Reset image size

Ultrasound is a wave of changing density and pressure where changes in density can be described by mechanical deformations (strains). Although we aim to model the focused ultrasound beam, a plane wave with the axis parallel to the axis of an axon was considered in this study (i.e. the wave propagated along the axon). This assumption allowed us to simplify analytical equations without loss of generality. Axons angled to the plane wave could be represented in the same way with addition of a scaling factor obtained using the angle between the plane wave axis to the axonal axis. The only guiding parameter here was the bulk strain, which did not significantly depend on this angle. So, a wave propagating in any direction would influence a section of the axon in a similar way. We assumed that deformation had a sinusoidal form along the axonal axis, with constant amplitude, which is in good agreement with existing literature (Canney et al 2008). We also assumed that the changes in ionic concentration were proportional to the changes in density and that ionic motion was several orders of magnitude lower than the wave speed (Ma et al 2018), thus also having a sinusoidal spatial waveform.

Under the mentioned assumptions, mechanical deformations (strains) are periodical in time and space. Two approaches were chosen to approximate them: a moving wave (1), which modelled the continuous ultrasound, and the standing wave (2), which better described the focused ultrasound:

$$\begin{eqnarray}&&\varepsilon =A\,\sin \left(2\pi f\left(t-\displaystyle \frac{x}{c}\right)\right)\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\varepsilon =A\,\sin \left(2\pi ft\right)\sin \left(\displaystyle \frac{2\pi x}{\displaystyle \frac{c}{f}}\right),\end{eqnarray} \tag{ 2 }$$

where $\varepsilon$—mechanical strain along the axis of ultrasound wave propagation, $t$—time,$x$—axis coordinate of an axon, $f$—frequency, $c$—the speed of sound in this medium, $A$—amplitude of the mechanical strain.

The only unknown parameter here is the strain amplitude $A.$ Since the pressure in the tissue under ultrasound is known (it can reach $1\,{\rm{MPa}}$ and more), the strain amplitude could be found from constitutive equations governing the elastic behaviour of the media. We assume that the surrounding liquid is compressible and cannot pass through the membrane at all. Then, mechanical strain in the membrane in any direction is simply related to mechanical stress:

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{\sigma }{3K},\end{eqnarray} \tag{ 3 }$$

where $K$ is compressibility modulus (bulk modulus) of the media which includes both axonal membrane and surrounding liquid, and $\sigma$ is the average normal stress, which for the considered case of volumetric compression/expansion equals to the applied pressure with the opposite sign. However, the actual estimation of $K$ for a brain poses significant controversy in the literature. Some researchers claim that the bulk modulus of the brain tissue must be close to that of water, $K\approx 2\,{\rm{GPa}}$ (Ganpule et al 2018). Similar results can be obtained from the measurements of ultrasound propagation velocity (Etoh et al 1994, Zoric et al 1999, Reyes Hernandez et al 2018), which under the assumption that elastic wave propagation speed $c=\sqrt{K/\rho },$ where $\rho$—density of a medium, allows estimation of $K\approx 2.4\,{\rm{GPa}}.$ A usual acoustic pressure $\sigma \,=\,1\,{\rm{MPa}}$ (Gavrilov et al 1977, Juan et al 2014, Blackmore et al 2019, Feng et al 2019) leads to the strain amplitudes of the order of $0.0001\,(0.01 \% ).$ On the other hand, assuming that the brain tissue is a linear isotropic elastic material, one can use relationships between elastic constants ($K=\tfrac{E}{3\left(1-2\nu \right)}$). This approach gives us the much smaller value of $K$ of the brain tissue: $10\,{\rm{kPa}}$ according to Morin et al (2017) and about $1\,{\rm{MPa}}$ in Miller et al (2000), Miller and Chinzei (2002), Soza et al (2004), Laksari et al (2012). In those cases, the same acoustic pressure $\sigma =1\,{\rm{MPa}}$ leads to significantly higher strains up to $1\,(100 \% ).$ Experimental and analytical estimations of deformations in the peripheral nervous system also show that the tissue can deform up to 1 (100%) during the application of ultrasound (Julian and Goldman 1962, Downs et al 2018).

There are therefore several contradicting pieces of evidence that lead to almost opposite results with respect to mechanical deformations of the membrane and do not allow to definitively state the relationship between the acoustic pressure and the mechanical deformations induced by the ultrasound. Here the influence of the mechanical deformations on the electrical behaviour of the axon was investigated, irrespective of the acoustic pressure. The detailed implications of using different pressure-deformation models are then provided in the discussion.

In this section, the analytical relationship between the ultrasound-induced mechanical deformation and electric cable theory boundary conditions (activating function) is delivered. The voltage across the membrane of the axon could be expressed using the charge and the membrane capacitance through the following equations:

$$\begin{eqnarray}&&{V}_{m0}=\displaystyle \frac{Q}{{C}_{m0}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{V}_{m}=\displaystyle \frac{Q}{{C}_{m}},\end{eqnarray} \tag{ 5 }$$

where ${V}_{m0}$—the resting potential of a membrane, ${V}_{m}$—membrane potential, $Q$—full charge, ${C}_{0m}$ and ${C}_{m}$—resting and changed capacitances.

There are two ways to express the changes to the voltage due to mechanical deformation. We considered an infinitely small cylindrical cut-off of the axon as shown in figure 2. The rate of mechanical deformations considered in this study is much quicker than transmembrane ionic movement (Hille 2001), therefore the total charge on both sides of the axon was assumed to be constant.

Zoom In Zoom Out Reset image size

The voltage across the membrane for the deformed state relative to the undeformed state can be expressed as a ratio of the cylindrical membrane capacitors, where we assumed that the membrane thickness did not change:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{V}_{m}(x)}{{V}_{m0}}=\displaystyle \frac{{C}_{m0}}{{C}_{m}(x)}=\displaystyle \frac{{A}_{m0}}{{A}_{m}(x)}=\displaystyle \frac{l\,2\pi R\,}{2\pi (l+dl)(R+dR)}\\ \,=\displaystyle \frac{1}{{\left(1+\varepsilon \left(x\right)\right)}^{2}}\sim 1-2\varepsilon (x),\end{array}\end{eqnarray} \tag{ 6 }$$

where ${A}_{0m}$ and ${A}_{m}$ are areas of a membrane conductor, $R$—radius of a membrane, $dR$—radius change due to ultrasound, $l$—length of an axon, $dl$—length change due to ultrasound, $\varepsilon$—strain (deformation) ($\varepsilon =\tfrac{dR}{R}=\left.\tfrac{dl}{l}\right).$ We considered the strain in the point to be uniform in an infinitely small volume around the segment. For the small strains ($\varepsilon (x)\ll 1$) we can express $\tfrac{1}{{\left(1+\varepsilon \left(x\right)\right)}^{2}}$ in a Taylor series as $1-2\varepsilon (x).$

The same final relationships can be also obtained by considering the infinitely small section of the membrane and calculating the change in concentration of the charge due to the volume changes, assuming the capacitance of the section stays the same. The second spatial derivative of external voltage along the axis of the axon, which is commonly called the activating function (Rattay 1989) can be expressed by the following set of equations according to the cable theory (Hodgkin and Rushton 1946, Rall 1959):

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{V}_{e}}{\partial {x}^{2}}=-{\left(1+\displaystyle \frac{{R}_{i}}{{R}_{e}}\right)}^{-1}\displaystyle \frac{{\partial }^{2}{V}_{m}}{\partial {x}^{2}},\end{eqnarray} \tag{ 7 }$$

where ${R}_{i}$ and ${R}_{e}$ are electrical resistivities of intracellular and extracellular space respectively. Substituting the relationship (6) for the mechanical strains results in:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\partial }^{2}{V}_{e}}{\partial {x}^{2}}=-{\left(1+\displaystyle \frac{{R}_{i}}{{R}_{e}}\right)}^{-1}\displaystyle \frac{{\partial }^{2}\left({V}_{0}\displaystyle \frac{1}{{\left(1+\varepsilon \left(x\right)\right)}^{2}}\right)}{\partial {x}^{2}}\\ \,\approx -{\left(1+\displaystyle \frac{{R}_{i}}{{R}_{e}}\right)}^{-1}\displaystyle \frac{{\partial }^{2}({V}_{0}(1-2\varepsilon (x)))}{\partial {x}^{2}}.\end{array}\end{eqnarray} \tag{ 8 }$$

During the ultrasound propagation, the strain is ${\varepsilon }_{0}\,\sin \left(\tfrac{2\pi x}{\lambda }\right)$ as in (1, 2), where ${\varepsilon }_{0}$ is a strain amplitude, and $\lambda$ is a spatial wavelength of the ultrasonic wave. Even focused ultrasound deformations can be approximated by sinusoidal function (Canney et al 2008). In this case, the equation for the activating function takes the form:

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{V}_{e}}{\partial {x}^{2}}=-2{\varepsilon }_{0}{\left(\displaystyle \frac{2\pi }{\lambda }\right)}^{2}{V}_{0}{\left(1+\displaystyle \frac{{R}_{i}}{{R}_{e}}\right)}^{-1}\,\sin \left(\displaystyle \frac{2\pi x}{\lambda }\right).\end{eqnarray} \tag{ 9 }$$

Assuming that $\tfrac{{\partial }^{2}{V}_{e}}{\partial {x}^{2}}={A}_{V}\,\sin \left(\tfrac{2\pi x}{\lambda }\right),$ we can get an amplitude:

$$\begin{eqnarray}&&{A}_{V}=-2{\varepsilon }_{0}{V}_{0}{\left(\displaystyle \frac{2\pi }{\lambda }\right)}^{2}{\left(1+\displaystyle \frac{{R}_{i}}{{R}_{e}}\right)}^{-1}.\end{eqnarray} \tag{ 10 }$$

For $\lambda =0.1\,{\rm{cm}},$ ${R}_{i}=0.05\,$kΩ cm (HH), ${V}_{0}=60\,{\rm{mV}},$ ${R}_{i}=0.0354\,$kΩ cm (C-Fibre) and ${R}_{e}=0.1\,$kΩ cm we have ${A}_{V}\approx 1.4\cdot {10}^{5}\,{\varepsilon }_{0}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for the HH model and ${A}_{V}\approx 1.6\cdot {10}^{5}\,{\varepsilon }_{0}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for C-Fibre model. The above relationships allow us to connect the amplitude of strain to the amplitude of activating function that we apply as a boundary condition to the simulated axon.

If strains are not small (i.e. higher than $0.01\,(1 \% )$), we have:

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{V}_{e}}{\partial {x}^{2}}=-{\left(1+\displaystyle \frac{{\rho }_{i}}{{\rho }_{e}}\right)}^{-1}\displaystyle \frac{{\partial }^{2}\left({V}_{0}\displaystyle \frac{1}{{\left(1+{\varepsilon }_{0}\,\sin \left(\displaystyle \frac{2\pi x}{\lambda }\right)\right)}^{2}}\right)}{\partial {x}^{2}}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{V}_{e}}{\partial {x}^{2}}=-8{\left(1+\displaystyle \frac{{\rho }_{i}}{{\rho }_{e}}\right)}^{-1}\displaystyle \frac{{V}_{0}{\varepsilon }_{0}{\pi }^{2}}{{\lambda }^{2}}\left(\displaystyle \frac{3{\varepsilon }_{0}{\left(\cos \left(\displaystyle \frac{2\pi }{\lambda }x\right)\right)}^{2}}{{\left(1+{\varepsilon }_{0}\,\sin \left(\displaystyle \frac{2\pi }{\lambda }x\right)\right)}^{4}}+\displaystyle \frac{\sin \left(\displaystyle \frac{2\pi }{\lambda }x\right)}{{\left(1+{\varepsilon }_{0}\,\sin \left(\displaystyle \frac{2\pi }{\lambda }x\right)\right)}^{3}}\right).\end{eqnarray} \tag{ 12 }$$

This equation allows us to precisely calculate an amplitude of activating function for strains close to $1\,(100 \% )$ and can be solved numerically. For example, when ${\varepsilon }_{0}=0.5,$ $\lambda =0.1\,{\rm{cm}},$ ${R}_{i}=0.05\,$kΩ cm (HH), ${R}_{i}=0.0354\,$kΩ cm (C-Fibre) and ${R}_{e}=0.1\,$kΩ cm we have ${A}_{V}\approx 2\cdot {10}^{6}\,\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for HH model and ${A}_{V}\approx 2.2\cdot {10}^{6}\,\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for the C-Fibre model. Similar, for ${\varepsilon }_{0}=0.9,$ we have ${A}_{V}\approx {10}^{8}\,\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for HH model and ${A}_{V}\approx 1.1\cdot {10}^{8}\,\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for C-Fibre model assuming the same parameters.

The above expressions show that the ultrasound stimulation could be treated based on the same principles as the application of the electric current and its effect could be modelled using the same simulation concepts.

There are two main forms of activating functions which we used in our study, simulating the moving and standing waves of ultrasound respectively, as per (1) and (2), and assuming that there are no nonlinear and transient effects:

$$\begin{eqnarray}&&{f}_{act}\left(x,t\right)=\displaystyle \frac{{\partial }^{2}{V}_{e}}{\partial {x}^{2}}={f}_{act0}\,\sin \left(2\pi f\left(t-\displaystyle \frac{x}{c}\right)\right)\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{f}_{act}\left(x,t\right)=\displaystyle \frac{{\partial }^{2}{V}_{e}}{\partial {x}^{2}}={f}_{act0}\,\sin \left(2\pi ft\right)\sin \left(\displaystyle \frac{2\pi x}{\displaystyle \frac{c}{f}}\right).\end{eqnarray} \tag{ 14 }$$

Two types of unmyelinated nerve fibres with active ion channels were simulated in this study. The first type was based on the HH model of the giant axon of the squid (HH model) (Hodgkin and Huxley 1952), which consisted of three types of ion channels—active sodium and potassium channels as well as passive leakage. The second type of fibre was a realistic mammalian C nociceptor with $10$ active ion channels and voltage-dependent ions' concentrations (Tigerholm et al 2014). This model was required to verify the qualitative results of the HH model and more accurately predict the quantitative values of the considered parameters, as well as give direct relevance to the obtained results. The C-fibre model closely resembles the fibres that can be found in the vagus nerve, stimulation of which is the prime interest of the scientific community for the treatment of a wide variety of drug-resistant disorders (Thompson et al 2019). These two models were used in this study to allow us both to investigate general and very common reliable model (HH) and to check the translation ability of the approach, and the validity of the results on a highly nonlinear detailed system (C-fibre).

In the models, fibres were represented as one-dimensional cables placed in a cylindrical extracellular space simulated using the FEM approach in COMSOL Multiphysics software, as in (Tarotin et al 2019). In this study, bi-directional coupling of the fibres and extracellular space was omitted since only the effect of the external electrical and acoustic fields on the fibres were investigated. The effect of the external electric field on the fibres was simulated using the concept of activating function (Rattay 1989). All the parameters of the models and equations representing their electrical behaviour were derived from (Tarotin et al 2019).

We have used a $25\,{\rm{cm}}$ and $7\,{\rm{cm}}$ long axons for the HH and C-fibre models respectively. The speed of sound was chosen to be equal to $1500\,{\rm{m}}\,{{\rm{s}}}^{-1}.$ The focus size depends on the wavelength (frequency). If the wavelength of the ultrasound was higher than $1\,{\rm{cm}},$ it was spaced on the area of a single wavelength; otherwise, it was spaced on the length of not less than $1\,{\rm{cm}}.$ This assumption was made to allow one-to-one comparison between different parameters, which would not be possible using the focused approximation, as the focal region was blurred and changed with frequency. Prior calculations show no relationship between the number of waves spaced along an axon and the activation threshold. The possible reason for this is that all the effects investigated in our study occur at the edge of a stimulated zone.

The zone of stimulation was located $5\,{\rm{cm}}$ and $1.5\,{\rm{cm}}$ from the left end of the axon for the HH and C-fibre models respectively. Element size was selected in such a way that there were no less than $32$ and $43$ elements within one ultrasound wavelength for the HH and C-fibre models respectively, to satisfy mesh convergence criteria. The adaptive time step was used for simulations: there were at least $20$ temporal steps for one period of ultrasound stimulation. When ultrasound stimulation was stopped there were at least $1000$ temporal steps per $1\,{\rm{ms}}.$

There were two different simulation cases in our investigation: continuous and pulsed (figure 3(A)). In the continuous case, the ultrasound was applied without any break continuously during the whole simulation. Different durations were utilised, and both standing and moving waves were considered. The simulations lasted at least $1\,{\rm{ms}}$ after the ultrasound modulation so that the AP generation and propagation could be detected. Two main parameters were under investigation: ultrasound amplitude and frequency.

Zoom In Zoom Out Reset image size

In the pulsed case, a short pulse of ultrasound was applied (figure 3(B)). The pulse duration was chosen to be $5* T,$ where $T$ is a period of the ultrasound wave. The frequency was chosen to be $1\,{\rm{MHz}}$ as in the majority of the existing experimental studies.

We have investigated the influence of the stop phase and ultrasound pulse durations on the AP generation. Stop phase is a phase at which the ultrasound pulse is stopped, ranged from 0 to $2\pi$ as a phase of the sinusoidal function. In the finite element (FE) model it is a phase of an activating function which was applied at the last moment of each pulse for a pulsed ultrasound and at the last moment of a continuous ultrasound (which is equivalent to the physical meaning as far as the effect of ultrasound is approximated by the activating function).

A wide range of phases and amplitudes (table 1) were studied to understand the relationship between parameters and the generation and propagation of the AP, measured at $1\,{\rm{mm}}$ left from the ultrasound stimulated area for HH model and at $0.1\,{\rm{mm}}$ for C-Fibre model.

Table 1. Parameter ranges for all studies.

Study	Model	Parameters
		Frequency, ${\rm{MHz}}$	Amplitude, $\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$	Pulse width, $\mathrm{number}\,$ ${\rm{of}}\,{\rm{periods}}$	Stop phase, $\mathrm{part}\,$ ${\rm{of}}\,{\rm{a}}\,{\rm{period}}$	Inter-stimulus interval, ${\rm{ms}}$	Number of pulses
Continuous	HH	$0.01\,-\,10$	${10}^{2}-{10}^{10}$
	C-Fibre	$0.01\,-\,1$	${10}^{2}-{10}^{10}$
One pulse (amplitude and stop phase)	HH	$1$	${10}^{2}-{10}^{7}$	$5$	$0-2\pi$
	C-Fibre	$1$	${10}^{7}-{10}^{9}$	$5$	$0-2\pi$
One pulse (pulse width)	HH	$1$	${10}^{6}$	$1-10000$	$0,\tfrac{2\pi }{5},\,\pi ,\tfrac{8\pi }{5}$
	C-Fibre	$1$	$5\cdot {10}^{7}$	$10-5000$	$0,\,\pi$
Pulsed	HH	$1$	${10}^{6}$	$1-1000$	$0,\tfrac{3\pi }{5},\pi$	$0.005-0.05$	$2-20$
	C-Fibre	$1$	${10}^{6}-{10}^{8}$	$10,\,100$	$0,\,\pi$	$0.01,\,0.05$	$2,\,5,\,10$
Pulsed (random phase)	HH	$1$	${10}^{6}$	$10,\,50,\,100$	$random$	$0.01$	$5,\,10,\,20$
	C-Fibre

The influence of the pulse duration on the AP generation was investigated at frequencies of $1\,{\rm{MHz}}$ and amplitudes of the activating function of ${10}^{6}$ $\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for HH model and $5\cdot {10}^{7}$ $\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for C-Fibre model (the amplitudes were determined as a threshold for the AP generation for the respective models).

A series of pulsed ultrasound was applied for close resemblance to experiments and real applications (Gavrilov et al 1977, King et al 2013, Kim et al 2014, Blackmore et al 2019). Two new parameters were under investigation: number of pulses and duration of an inter-stimulus interval between pulses. Stop phase and duration of a pulse were also varied. Frequency and amplitude were chosen to be $1\,{\rm{MHz}}$ and ${10}^{6}$ $\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ for HH model and $1\,{\rm{MHz}}$ with different amplitudes for C-Fibre model.

Finally, the probabilistic behaviour of ultrasound stimulation was investigated. Stop phases were randomly distributed from 0 to $2\pi ,$ with the pulse length of $100$ periods, and inter-stimulus interval of $0.01\,{\rm{ms}}.$ Frequency and amplitude were 1 ${\rm{MHz}}$ and ${10}^{6}$ $\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ respectively. The study was conducted only for the HH model as it was not possible to compute for the C-fibre due to computational difficulties.

Both HH and C-fibre models demonstrated that the proposed approach could reliably generate initiation and propagation of the AP. The ultrasound-induced APs were $43.196\,{\rm{mV}}$ and $26.6\,{\rm{mV}}$ in membrane potential for HH and C-fibre model respectively, and their corresponding velocities were $1.5$ and $0.08\,{\rm{cm}}\,{{\rm{ms}}}^{-1},$ which matches the literature for the typical slow-fibre APs (Hodgkin and Huxley 1952, Rosenthal and Bezanilla 2000, Tarotin et al 2019). The example studies demonstrated the stable AP initiation process and consequent propagation along the axon (figures 4–6) at various ultrasound modulation parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The comparison between standing and moving waves at frequencies of $100\,{\rm{kHz}}$ with continuous ultrasound stimulation revealed that the standing wave has a lower threshold of AP generation by approximately half. At lower frequencies, this effect decreased so that below $10\,{\rm{kHz}}$ there was no significant difference between standing and moving waves. Significantly higher amplitude was required for AP generation at higher frequencies in the case of moving wave, doubling the one for the standing wave at frequencies above $1\,{\rm{MHz}}.$ We were using a standing wave for the subsequent investigations as it was closer to the actual focused ultrasound used for neuromodulation.

To find the activation threshold for the continuous ultrasound, different frequencies and amplitudes of activating function were applied for the whole duration of the simulation. The activating threshold varied with frequency (figure 7) and was ${10}^{3}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{{\rm{2}}}}$ for $10\,{\rm{kHz}},$ and increased to ${10}^{8}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{{\rm{2}}}}$ for $10\,{\rm{MHz}}$ for the HH model, which corresponded to the mechanical strains of the axon by $0.68\,(68 \% )$ and $0.41\,(41 \% )$ respectively. In figure 7 one can also find strain for each simulation case written in each cell. These strains show that an activation threshold in terms of strain is similar at all frequencies. However, there is a slight decrease of strains (from 0.68 to 0.41) for frequencies higher than $1\,{\rm{MHz}}.$ The same analysis for the C-fibre resulted in the threshold ${10}^{7}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{{\rm{2}}}}$ for $10\,{\rm{kHz}}$ and ${10}^{10}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{{\rm{2}}}}$ for $1\,{\rm{MHz}}$ respectively, which corresponded to $0.98$ and $0.96$ mechanical strain. The activation threshold in terms of strain is similar at all frequencies for the C-fibre model.

Zoom In Zoom Out Reset image size

3.3.1. Phase–amplitude relationship for the single pulse stimulation

The activation threshold for the pulsed ultrasound was found for different stop phases and the length of the pulse of $5$ periods ($5$ $\mu {\rm{s}}$) at the frequency of 1 MHz (figure 8(A)). AP was generated when the amplitudes of the activating function were larger than $5\cdot {10}^{6}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.61 strain). The stop phases between $\tfrac{2\pi }{5}$ and $\tfrac{8\pi }{5}$ (phases are expressed as a phase of sinusoidal function, see figure 3) caused the highest membrane voltage change and subsequent AP generation. This suggested that the generation of AP was more sensitive to the total energy translated to the system rather than the maximum membrane potential (figure 8, red dots). The highest energy value translated to the system was at the phase of $\pi$ (the wave stopped just after the positive phase), which corresponded to the highest membrane voltage for all amplitudes. The amplitudes of the activating function below ${10}^{4}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.05 strain) did not significantly affect the membrane potential. For the C-fibre model the pattern was very similar (figure 8). Phases between $\tfrac{3\pi }{5}$ and $\tfrac{7\pi }{5}$ were the most effective ones. The amplitudes below ${10}^{7}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.8 strain) did not significantly change the membrane voltage.

Zoom In Zoom Out Reset image size

3.3.2. Pulse width analysis

Pulse widths and stop phases were varied at $1\,{\rm{MHz}}$ with the activating function amplitudes ${10}^{6}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.41 strain) and $5\cdot {10}^{7}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.8 strain) for the HH and C-Fibre models respectively (figure 9). The membrane potential in the HH model decreased with the increase of pulse width. As in the previous case, the stop-phase significantly affected the results with $\pi$ being optimal for the generation of the AP.

Zoom In Zoom Out Reset image size

C-Fibre model behaved differently (figure 9(B)). Membrane potential increased with the increase of a pulse width when the ultrasound was stopped at the zero stop phase. However, if the wave was stopped at $\pi ,$ the membrane potential decreased with the increase of pulse width. When pulse width reached the length of $1000$ periods there was no significant difference in the output between the two stop phases.

3.3.3. Analysis of the ultrasound modulation with pulse trains.

Here the activating thresholds and membrane potentials were investigated with respect to the number of pulses, pulse widths, stop phases, and interstimulus intervals (table 1) for the pulsed ultrasound delivered using the series of pulses, or a train (figures 10–12). A range of pulses of activating function with different widths was applied at $1\,{\rm{MHz}}.$ The amplitudes were ${10}^{6}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.41 strain) for the HH model and ${10}^{6}-{10}^{8}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.40–0.83 strain) for C-Fibre model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the HH model, inter-stimulus interval duration did not significantly affect the membrane voltage changes, however shorter inter-stimulus intervals resulted in generally higher membrane potentials (figure 10). The increase in the number of pulses was generally more effective in generating the AP (figure 11). Pulses for which the stop phase was $0$ did not generate APs at all for our study parameters. The increase of a pulse width resulted in a significant decrease in the AP generation ability (figure 11).

C-Fibre model behaved in a similar way (figure 12). Amplitude had a major influence on the qualitative behaviour of the system. The number of pulses significantly changed the activating threshold at the amplitude of $5\cdot {10}^{7}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.8 strain) (figures 12(E), (F)).

Additional experiment with $0$ stop phase showed that axon did not generate the AP at amplitude 1$e8\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}},$ $10$ pulses, $0.01$ ${\rm{ms}}$ interstimulus interval and $100$ periods ($100\,\mu {\rm{s}}$) pulse width, whilst the AP was generated with a $\pi$ stop phase using the exact same parameters.

3.3.4. Phase–probability relationship: prediction of the experimental outcomes given random phase distribution per pulse.

In real experiments, the stop phase of the ultrasound sometimes is not controlled precisely and is therefore randomly distributed between the ultrasound pulses. This was modelled in the current subsection (figure 13) with the ultrasound pulses at a frequency of $1\,{\rm{MHz}}$ with the amplitude of ${10}^{6}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ (0.43 strain), and with the $0.01$ ms interstimulus interval. The number of pulses and pulse width were varied respectively.

Zoom In Zoom Out Reset image size

The stimulation of an axon had probabilistic nature, and the probability of AP generation significantly increased with a higher number of pulses ($10$ and $20$) and smaller pulse widths (mostly $10$ and $50$ periods, figure 13).

The model of ionic concentration changes as a means of AP generation during ultrasound stimulation showed that this mechanism could indeed theoretically be responsible for ultrasound neuromodulation in a variety of cases. The continuous ultrasound case showed that the activating function amplitudes required for AP generation were ${10}^{5}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ and ${10}^{8}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ at $100\,{\rm{kHz}},$ and ${10}^{7}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ and ${10}^{10}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ at $1\,{\rm{MHz}}$ for HH and C-fibre models respectively. This corresponded to the strain values of $0.68$ and $0.96$ at $100\,{\rm{kHz}}$ and $0.68$ and $0.96$ at $1\,{\rm{MHz}}.$ The strains were equal at both frequencies, that indirectly supported the deformative nature of the effect of ultrasound stimulation on axons.

The pulsed ultrasound stimulation cases showed that the activation threshold varied with the stop phase, the number of pulses in the train, pulse width, and inter-stimulus interval. The stop phase had the most significant effect on the overall probability of AP generation, with the highest probability at the stop phase of $\pi$ for both HH and C-fibre models.

The lowest thresholds for the AP generation with the pulsed ultrasound were ${10}^{6}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ and $5\cdot {10}^{7}\tfrac{{\rm{mV}}}{{{\rm{cm}}}^{2}}$ at $1\,\,{\rm{MHz}}$ for HH and C-fibre models respectively, corresponded to the strain values of $0.41$ and $0.8,$ using the 5 pulses with 10 periods each, which were all stopped at $\pi$ stop phase. The analysis of the energy transfer in the system indicated that the first half of the sine wave generated the positive activating function which had a cumulative effect on the generation of the AP.

4.2.1. Can this mechanism be responsible for AP generation by ultrasound stimulation?

4.2.2. What are the implications for the experimental predictions?

Here we only considered unmyelinated fibres. We also greatly simplified the transition between mechanical strain and boundary conditions. Although we are confident about the qualitative results, we need to develop a more complicated model in order to study the quantitative effects thoroughly. The one-to-one comparison to other hypotheses would be ideal, however, this is challenging given the fact that most parameters for the other hypotheses are usually obtained experimentally.

In future, we are planning to extend the model for the simultaneous solution of Nernst–Plank equations and HH (and C-fibre) equations. We are also planning to create a myelinated coupled electro-acoustic model of the axon. We then are planning to confirm the qualitative results obtained in this work, obtain predictive quantitative results for particular nerves, and conduct an experimental study which would definitively answer the question on the exact influence of the ionic concentration mechanism on neurostimulation by the ultrasound.

One of the potential applications of ultrasound stimulation could be the possibility of affecting the direction of the information flow within the vagus nerve. For example, the threshold of stimulation could be different for A-delta efferent fibres and the rest of the myelinated fibres. Existing models of myelinated fibres (McIntyre et al 2002) offer the possibility of such investigations. We therefore will consider the possibility of assessing this and comparing the myelinated and unmyelinated fibre behaviour.