Description

Context

Low-Intensity Focused Ultrasound Stimulation (LIFUS) is emerging as a promising technology for noninvasive brain stimulation, but many questions remain regarding its mechanisms of action as well as optimal stimulation parameters.

One candidate mechanism is that of intramembrane cavitation, in which LIFUS induces high-frequency mechanical deflections of the neural membrane to generate depolarizing currents and trigger action potentials. This phernomenon is decribed mathematically by a Neuronal Intramembrane Cavitation Excitation (NICE) model [1-3], whose predictions of cell-type-specific LIFUS responses corroborate a wide range of experimental results on LIFUS-evoked brain activity. However, the NICE model is conceptually complex and tedious to simulate, thereby limiting its adoption by the community.

To adress these limitations, we recently developed the multiScale Optimized Intramembrane Cavitation (SONIC) model [4], a computationally efficient variant of the NICE model that also defines a more intuititve frame of interpretation for LIFUS-evoked responses.

The PySONIC package is a Python implementation of the SONIC model, allowing to simulate the responses of various neuron types to ultrasonic (and electrical) stimuli.

Content of repository

Models

The package contains four model classes:

Model defines the generic interface of a model, including mandatory attributes and methods for simulating it.
BilayerSonophore defines the underlying biomechanical model of intramembrane cavitation.
PointNeuron defines an abstract generic interface to conductance-based point-neuron electrical models. It is inherited by classes defining the different neuron types with specific membrane dynamics.
NeuronalBilayerSonophore defines the full electromechanical model for any given neuron type. To do so, it inherits from BilayerSonophore and receives a specific PointNeuron object at initialization.

All model classes contain a simulate method to simulate the underlying model's behavior for a given set of stimulation and physiological parameters. The NeuronalBilayerSonophore.simulate method contains an additional method argument defining whether to perform a detailed (full), coarse-grained (sonic) or hybrid (hybrid) integration of the differential system.

The default (and fastest) simulation method is the sonic method, and makes use of pre-computed tables that are stored in lookup files within the package architecture (in the lookups subfolder). These large binary files are handled by the git-lfs utility, which sets up dynamic links to these files without storing them physically in the repository, thereby avoiding to store their entire history. Hence, you must install git-lfs in order to download the lookup files together with the repo.

Solvers

Numerical integration routines are implemented outside the models, in a separate solvers module:

PeriodicSolver integrates a differential system periodically until a stable periodic behavior is detected.
EventDrivenSolver integrates a differential system across a specific set of "events" that modulate the stimuluation drive over time.
HybridSolver inherits from both PeriodicSolverand EventDrivenSolver. It integrates a differential system using a hybrid scheme inside each "ON" or "OFF" period. First, the full ODE system is integrated for a few cycles with a dense time granularity until a periodic stabilization detection. Then, the profiles of all variables over the last cycle are resampled to a far lower (i.e. sparse) sampling rate. Next, a subset of the ODE system is integrated with a sparse time granularity, while the remaining variables are periodically expanded from their last cycle profile, until the end of the period or that of an predefined update interval. This process is then repeated throughout the simulation.

Neurons

Several conductance-based point-neuron models are implemented that inherit from the PointNeuron generic interface.

Among them, several CNS models:

CorticalRS: cortical regular spiking (RS) neuron
CorticalFS: cortical fast spiking (FS) neuron
CorticalLTS: cortical low-threshold spiking (LTS) neuron
CorticalIB: cortical intrinsically bursting (IB) neuron
ThalamicRE: thalamic reticular (RE) neuron
ThalamoCortical: thalamo-cortical (TC) neuron
OstukaSTN: subthalamic nucleus (STN) neuron

But also some valuable models used in peripheral axon models:

FrankenhaeuserHuxleyNode: Amphibian (Xenopus) myelinated fiber node (FHnode)
SweeneyNode: Mammalian (rabbit) myelinated motor fiber node (SWnode)
MRGNode: Mammalian (human / cat) myelinated fiber node (MRGnode)
SundtSegment: Mammalian unmyelinated C-fiber segment (SUseg)

Other modules

batches: a generic interface to run simulation batches with or without multiprocessing
parsers: command line parsing utilities
plt: graphing utilities
postpro: post-processing utilities (mostly signal features detection)
constants: algorithmic constants used across modules and classes
utils: generic utilities

Requirements

Python 3.6+
Git Large File Storage (LFS)

Installation

Install git-lfs if not already done. This step is crucial as it ensures that lookup tables will be downloaded when you clone the repository.
Download and install Anaconda if not already done.
Open a terminal. For Windows users, you might need to use the Anaconda Prompt if you did not add Anaconda to the path.
Create a new Anaconda environment called sonic with an appropriate Python version (i.e., 3.6 or more recent), and activate it. For instance with Python 3.11:

conda create -n sonic python=3.11
conda activate sonic

Clone the PySONIC repository and install it as a python package:

git clone https://github.com/tjjlemaire/PySONIC.git
cd PySONIC
pip install -e .

All package dependencies (numpy, scipy, ...) should be installed automatically.

Usage

Python code

You can easily run simulations of any implemented point-neuron model under both electrical and ultrasonic stimuli with various temporal protocols, and visualize the simulation results, in just a few lines of code. Here's an example:

import logging
import matplotlib.pyplot as plt
from PySONIC.utils import logger
from PySONIC.neurons import getPointNeuron
from PySONIC.core import NeuronalBilayerSonophore, AcousticDrive
from PySONIC.core.protocols import *
from PySONIC.plt import GroupedTimeSeries
logger.setLevel(logging.INFO)

# Define point-neuron model and corresponding neuronal bilayer sonophore object
nbls = NeuronalBilayerSonophore(32e-9, getPointNeuron('RS'))

# Define stimulus drive and time protocol
drive = AcousticDrive(500e3, 100e3)
pp = PulsedProtocol(100e-3, 100e-3, tstart=10e-3)

# Run simulation and plot results
data, meta = nbls.simulate(drive, pp)
fig = GroupedTimeSeries([(data, meta)]).render()[0]
plt.show()

The output data is formatted as a pandas dataframe, hence you can easily export your simulation results to Excel/CSV for further processing with:

data.dropna(axis=1).to_excel('mysim.xlsx')

Similarly, you can export the output figure to PDF for further editing with:

fig.savefig('myfig.pdf')

For more code examples, have a look at the demo notebook.

From the command line

Running simulations and batches

You can easily run simulations of all 3 model types using the dedicated command line scripts. To do so, open a terminal in the scripts directory.

Use run_mech.py for simulations of the mechanical model upon ultrasonic stimulation. For instance, for a 32 nm radius bilayer sonophore sonicated at 500 kHz and 100 kPa:

python run_mech.py -a 32 -f 500 -A 100 -p Z

Use run_estim.py for simulations of point-neuron models upon intracellular electrical stimulation. For instance, a regular-spiking (RS) neuron injected with 10 mA/m2 intracellular current for 30 ms:

python run_estim.py -n RS -A 10 --tstim 30 -p Vm

Use run_astim.py for simulations of point-neuron models upon ultrasonic stimulation. For instance, for a coarse-grained simulation of a 32 nm radius bilayer sonophore within a regular-spiking (RS) neuron membrane, sonicated at 500 kHz and 100 kPa for 150 ms:

python run_astim.py -n RS -a 32 -f 500 -A 100 --tstim 150 --method sonic -p Qm

Additionally, you can run batches of simulations by specifying more than one value for any given stimulation parameter (e.g. -A 100 200 for sonication with 100 and 200 kPa respectively). These batches can be parallelized using multiprocessing to optimize performance, with the extra argument --mpi.

Saving and visualizing results

By default, simulation results are neither shown, nor saved.

To view results directly upon simulation completion, you can use the -p [xxx] option, where [xxx] can be all (to plot all resulting variables) or a given variable name (e.g. Z for membrane deflection, Vm for membrane potential, Qm for membrane charge density).

To save simulation results in binary .pkl files, you can use the -s option. You will be prompted to choose an output directory, unless you also specify it with the -o <output_directory> option. Output files are automatically named from model and simulation parameters to avoid ambiguity.

When running simulation batches, it is highly advised to specify the -s option in order to save results of each simulation. You can then visualize results at a later stage.

To visualize results, use the plot_timeseries.py script. You will be prompted to select the output files containing the simulation(s) results. By default, separate figures will be created for each simulation, showing the time

PySONIC

Install / Use

README