SONIC
SONIC Model — Simulation of Ultrasound Neuronal Interaction Circuits
Overview
SONIC (Lemaire et al. 2019) was developed to solve the fundamental computational intractability of the NICE model. Rather than numerically integrating the full coupled electromechanical system at acoustic timescales, SONIC introduces a coarse-graining procedure based on precomputed lookup tables that separate the fast acoustic dynamics from the slow neuronal dynamics. This reduces simulation time by several orders of magnitude while preserving the essential physics of the NICE framework, making it feasible to simulate realistic sonication protocols across multiple cell types and parameter combinations.
1. The Core Problem: Timescale Separation
As established in the NICE node, the full electromechanical system spans timescales from nanoseconds (acoustic oscillation) to hundreds of milliseconds (neuronal firing patterns). SONIC exploits the fact that these timescales are well-separated: the acoustic oscillation (μs) is much faster than ion channel gating kinetics (ms). This separation justifies replacing the rapidly oscillating quantities with their cycle-averaged effective values, computed once offline and stored in lookup tables.
The fundamental insight is that ion channel gating variables evolve slowly relative to the acoustic cycle. Therefore, at each instant of the slow neuronal simulation, the membrane is effectively experiencing a time-averaged acoustic drive rather than the full instantaneous oscillation.
2. The Coarse-Graining Procedure
For a given set of sonication parameters $(f, A)$ — carrier frequency and pressure amplitude — and a given membrane charge density $Q_m$, SONIC precomputes the cycle-averaged effective quantities by running the fast NICE bubble dynamics for one acoustic cycle:
\[\langle \dot{Q}_m \rangle = \frac{1}{T} \int_0^T \dot{Q}_m(t) \, dt\]
\[\langle \lambda_x(Q_m, A, f) \rangle = \frac{1}{T} \int_0^T \lambda_x(V_m(t)) \, dt\]
Where $\lambda_x$ represents the voltage-dependent rate constants (opening rate $\alpha_x$ and closing rate $\beta_x$) for each gating variable $x$ of each ion channel. These averaged quantities are functions only of the slow variables ($Q_m$ and gating states) and the fixed sonication parameters.
The resulting lookup tables store, for each combination of $(Q_m, A, f)$:
- $\langle \dot{Q}_m \rangle$: effective charge drift rate
- $\langle \alpha_x \rangle$, $\langle \beta_x \rangle$: effective gating rate constants for each channel
- $\langle V_m \rangle$: effective membrane voltage
3. The SONIC Effective System
Once lookup tables are constructed, the slow neuronal dynamics are governed by a modified HH system where all quantities are replaced by their cycle-averaged counterparts:
\[\frac{dQ_m}{dt} = \langle \dot{Q}_m(Q_m, A, f) \rangle\]
\[\frac{dx}{dt} = \langle \alpha_x(Q_m, A, f) \rangle (1 - x) - \langle \beta_x(Q_m, A, f) \rangle x\]
This system evolves on the slow (ms) timescale only, with the fast acoustic physics entirely encoded in the precomputed tables. The acoustic timescale is no longer present in the integration — eliminating the multi-scale stiffness problem of NICE entirely.
For the OFF phase (between sonication pulses), the standard HH system is recovered exactly, since at $A = 0$ the lookup tables reduce to the unperturbed rate constants.
4. Lookup Table Structure
The lookup tables are constructed over a discrete grid of:
| Dimension | Typical range | Resolution |
|---|---|---|
| Charge density $Q_m$ | −80 to +50 nC/cm² | ~1 nC/cm² |
| Pressure amplitude $A$ | 0 to ~600 kPa | ~10 kPa |
| Carrier frequency $f$ | 0.25 to 4 MHz | discrete values |
For each cell type, a separate set of tables is required because the gating rate functions $\alpha_x$, $\beta_x$ differ across neuron models. The table construction (the "offline" phase) requires running the full NICE bubble dynamics simulation for each grid point — computationally expensive once, but done only once per cell type and sonication frequency.
Computational gain: direct NICE integration requires timesteps of ~1 ns; SONIC operates at timesteps of ~0.1 ms — a speedup of approximately $105$, reducing hours-long NICE simulations to seconds.
5. Cell Types Modeled
Lemaire et al. (2019) constructed SONIC lookup tables and validated responses for:
| Cell type | Key distinguishing parameter |
|---|---|
| RS (Regular Spiking) pyramidal | Low $C_{m0}$, adapting firing |
| FS (Fast Spiking) interneuron | Low $C_{m0}$, non-adapting high-frequency firing |
| LTS (Low-Threshold Spiking) interneuron | Low $C_{m0}$, T-type Ca²⁺ current |
| TC (Thalamocortical) | HCN current ($I_h$), burst/tonic modes |
| TRN (Thalamic Reticular Nucleus) | Strong T-type Ca²⁺, burst firing |
| STN (Subthalamic Nucleus) | Autonomous pacemaking, mixed conductances |
A key finding was that $C_{m0}$ is the primary biophysical determinant of cell-type selectivity: cells with lower resting capacitance experience larger voltage excursions per unit charge displacement, making them more sensitive to the SONIC capacitive drive at a given acoustic amplitude.
6. Threshold and Response Maps
A major contribution of SONIC was the systematic mapping of firing threshold curves in the $(A, DC)$ parameter space for each cell type. Key findings:
- Threshold pressure decreases monotonically with increasing DC for all cell types
- At fixed DC, threshold pressure decreases as frequency approaches sonophore resonance (0.5–0.7 MHz)
- Non-monotonic intensity response: above a saturation amplitude, firing rate decreases and eventually ceases (depolarization block), reproducing King et al. (2013) experimental observations
- Cell-type ordering: at equivalent parameters, TC and STN neurons show lower thresholds than cortical RS/FS neurons, consistent with clinical observations of thalamic sensitivity to TUS
7. Validation
SONIC predictions were validated against:
- Direct NICE simulations: cycle-averaged SONIC output matches full NICE integration with < 5% error in spike timing across tested parameter ranges
- Experimental literature: threshold intensities and frequency dependence consistent with published electrophysiology data from Legon et al., King et al., and Deffieux et al.
- PRF and DC dependence: firing probability maps reproduce qualitative trends from in vitro and in vivo ultrasound stimulation studies
8. Strengths and Limitations
Strengths
- Computational efficiency: ~$105$ speedup over NICE enables multi-cell, multi-parameter, and protocol-level simulations
- Physical fidelity preserved: lookup tables encode exact NICE physics — no phenomenological approximations introduced
- Multi cell-type coverage: first framework to systematically compare RS, FS, LTS, TC, TRN, STN responses under identical sonication conditions
- Modular: lookup tables can be regenerated for any HH-compatible neuron model, making SONIC extensible to new cell types
- Reproduces non-linear phenomena: depolarization block, threshold behavior, PRF/DC interactions all emerge naturally
Limitations
- Point neuron: like NICE, SONIC models each neuron as a single isopotential compartment — no spatial integration of acoustic field across dendrites or axon
- Single sonophore assumption: one representative R₀ = 32 nm sonophore per membrane patch; biological heterogeneity in sonophore size distribution is ignored
- Lookup table interpolation error: response accuracy depends on table resolution; coarse grids introduce interpolation artifacts at steep response boundaries
- No MS channel coupling: membrane tension computed implicitly in the bubble dynamics is not exported to mechanosensitive channel gating equations
- No thermal module: sonication-induced heating and Q₁₀ effects on channel kinetics are absent
- Static sonophore coverage fraction: the fraction of membrane occupied by sonophores is fixed; in MorphoSONIC this is relaxed to vary with morphology
- Bilayer sonophore remains unconfirmed experimentally: all limitations of the physical hypothesis carry over from NICE