Artificial Neural Networks architect
Neurocomputational models
Neuro-simulator ANNarchy
Rate-coded networks
Spiking networks
Computational neuroscience is about explaining brain functioning at various levels (neural activity patterns, behavior, etc.) using biologically realistic neuro-computational models.
Different types of neural and synaptic mathematical models are used in the field, abstracting biological complexity at different levels.
There is no “right” level of biological plausibility for a model - you can always add more details -, but you have to find a paradigm that allows you to:
\tau \, \frac{d v(t)}{dt} + v(t) = \sum_{i=1}^d w_{i, j} \, r_i(t) + b
r(t) = f(v(t))
C \, \frac{d v(t)}{dt} = - g_L \, (v(t) - V_L) + I(t)
\text{if} \; v(t) > V_T \; \text{emit a spike and reset.}
\begin{cases} \displaystyle\frac{dv}{dt} = 0.04 \, v^2 + 5 \, v + 140 - u + I \\ \\ \displaystyle\frac{du}{dt} = a \, (b \, v - u) \\ \end{cases}
\begin{cases} \begin{aligned} C \, \frac{dv}{dt} = -g_L \ (v - E_L) + & g_L \, \Delta_T \, \exp(\frac{v - v_T}{\Delta_T}) \\ & + I - w \end{aligned}\\ \\ \tau_w \, \displaystyle\frac{dw}{dt} = a \, (v - E_L) - w\\ \end{cases}
Biological neurons do not all respond the same to an input current: Some fire regularly, some slow down with time., some emit bursts of spikes…
Modern spiking neuron models allow to recreate these dynamics by changing just a few parameters.
Oscillations at the population level.
Excitatory/inhibitory balance.
Spatio-temporal separation of inputs (reservoir computing).
When an axon of cell A is near enough to excite a cell B and repeatedly or persistently takes part in firing it, some growth process or metabolic change takes place in one or both cells such that A’s efficiency, as one of the cells firing B, is increased.
Donald Hebb, 1949
\frac{dw}{dt} = \eta \, r^\text{pre} \, r^\text{post}
\frac{dw}{dt} = \eta \, r^\text{pre} \, r^\text{post} \, (r^\text{post} - \mathbb{E}[(r^\text{post})^2])
\frac{dw}{dt} = \eta \, r^\text{pre} \, (r^\text{post} - \mathbb{E}[r^\text{post}])
\frac{dw}{dt}= \eta \, r^\text{pre} \, r^\text{post} - \alpha \, (r^\text{post})^2 \, w
or virtually anything depending only on the pre- and post-synaptic firing rates, e.g. Vitay and Hamker (2010):
\begin{aligned} \frac{dw}{dt} & = \eta \, ( \text{DA}(t) - \overline{\text{DA}}) \, (r^\text{post} - \mathbb{E}[r^\text{post}] )^+ \, (r^\text{pre} - \mathbb{E}[r^\text{pre}])- \alpha(t) \, ((r^\text{post} - \mathbb{E}[r^\text{post}] )^+ )^2 \, w \end{aligned}
Synaptic efficiencies actually evolve depending on the the causation between the neuron’s firing patterns:
If the pre-synaptic neuron fires before the post-synaptic one, the weight is increased (long-term potentiation). Pre causes Post to fire.
If it fires after, the weight is decreased (long-term depression). Pre does not cause Post to fire.
Bi and Poo (2001) Synaptic Modification by Correlated Activity: Hebb’s Postulate Revisited. Annual Review of Neuroscience 24.
\frac{dw}{dt} = \begin{cases} A^+ \, \exp - \frac{t_\text{pre} - t_\text{post}}{\tau^+} \; \text{if} \; t_\text{post} > t_\text{pre}\\ \\ A^- \, \exp - \frac{t_\text{pre} - t_\text{post}}{\tau^-} \; \text{if} \; t_\text{pre} > t_\text{post}\\ \end{cases}
STDP can be implemented online using traces.
More complex variants of STDP (triplet STDP) exist, but this is the main model of synaptic plasticity in spiking networks.
Bi and Poo (2001) Synaptic Modification by Correlated Activity: Hebb’s Postulate Revisited. Annual Review of Neuroscience 24.
Populations of neurons can be combined in functional neuro-computational models learning to solve various tasks.
Need to implement one (or more) equations per neuron and synapse. There can thousands of neurons and millions of synapses, see Teichmann et al. (2021).
Fixed libraries of models
NEURON
https://neuron.yale.edu/neuron/
NEST
https://nest-initiative.org/
GeNN
https://genn-team.github.io/genn/
Auryn
https://fzenke.net/auryn/doku.php
Code generation
Brian
https://briansimulator.org/
Brian2CUDA
https://github.com/brian-team/brian2cuda
ANNarchy
https://github.com/ANNarchy/ANNarchy
Vitay et al. (2015)
ANNarchy: a code generation approach to neural simulations on parallel hardware.
Frontiers in Neuroinformatics 9. doi:10.3389/fninf.2015.00019
Installation guide: https://annarchy.github.io/Installation/
Using pip:
From source:
Requirements (Linux and MacOS):
Simulation of both rate-coded and spiking neural networks.
Only local biologically realistic mechanisms are possible (no backpropagation).
Equation-oriented description of neural/synaptic dynamics (à la Brian).
Code generation in C++, parallelized using OpenMP on CPU and CUDA on GPU (MPI is coming).
Synaptic, intrinsic and structural plasticity mechanisms.
import ANNarchy as ann
# Create network
net = ann.Network()
# Create neuron and synapse types for transmission and/or plasticity
neuron = ann.Neuron(...)
stdp = ann.Synapse(...)
# Create populations of neurons
pop = net.create(1000, neuron)
# Connect the populations through projections
proj = net.connect(pop, pop, 'exc', stdp)
proj.fixed_probability(weights=ann.Uniform(0.0, 1.0), probability=0.1)
# Generate and compile the code
net.compile()
# Simulate for 1 second
net.simulate(1000.)
\tau \frac{dx(t)}{dt} + x(t) = \sum w^\text{in} \, r^\text{in}(t) + g \, \sum w^\text{rec} \, r(t) + \xi(t)
r(t) = \tanh(x(t))
Parameters can have one value per neuron in the population (local) or be common to all neurons of the population (global).
The parameters are considered global by default. To define a local parameters, use ann.Parameter()
:
int
, bool
).The output variable of a rate-coded neuron must be named r
.
Variables can be updated with assignments (=
, +=
, etc) or by defining first order ODEs.
The math C library symbols can be used (tanh
, cos
, exp
, etc).
Initial values can be specified by using an ann.Variable()
object (default: 0.0).
min
/max
flags:equations = [
ann.Variable('tau * dx/dt + x = sum(exc)', init=-1.0, max=1.0),
ann.Variable('r = x', min=0.0), # ReLU
]
Uniform
, Normal
, LogNormal
, Exponential
, Gamma
…sympy
:dt
) for all neurons:Several numerical methods are available:
Network
is the main object in ANNarchy:pop = net.create(
geometry=1000,
neuron=ESN_Neuron
)
proj = net.connect(
pre=pop,
post=pop,
target='exc'
)
dt
and the seed of the RNG:Network
subclasses:class SimpleNetwork(ann.Network):
def __init__(self, size:int):
self.pop = self.create(
geometry=size,
neuron=ESN_Neuron
)
self.proj = self.connect(
pre=self.pop,
post=self.pop,
target='exc'
)
net = SimpleNetwork(1000)
Network
subclasses allow to run simulations in parallel.PopulationView
and can be addressed separately:Each target ('exc', 'inh', 'AMPA', 'NMDA', 'GABA'
) can be defined as needed and will be treated differently by the post-synaptic neurons.
The weighted sum of inputs for a specific target is accessed in the equations by sum(target)
:
equations = [
ann.Variable('tau * dx/dt + x = sum(exc) / (1 + sum(inh))'),
ann.Variable('r = tanh(x)'),
]
Projections must be populated with a connectivity matrix (who is connected to who), a weight w
and optionally a delay d
(uniform or variable).
Several patterns are predefined:
proj.all_to_all(weights=ann.Normal(0.0, 1.0), delays=2.0, allow_self_connections=False)
proj.fixed_probability(probability=0.2, weights=1.0)
proj.one_to_one(weights=1.0, delays=ann.Uniform(1.0, 10.0))
proj.fixed_number_pre(number=20, weights=1.0)
proj.fixed_number_post(number=20, weights=1.0)
proj.gaussian(amp=1.0, sigma=0.2, limit=0.001)
proj.dog(amp_pos=1.0, sigma_pos=0.2, amp_neg=0.3, sigma_neg=0.7, limit=0.001)
You can now manipulate all parameters/variables from Python.
A simulation is simply run for a fixed duration in milliseconds with:
https://annarchy.github.io/manual/Network.html#early-stopping
By default, a simulation is run in C++ without interaction with Python.
You may want to record some variables (neural or synaptic) during the simulation with a Monitor
:
Calling get()
flushes the underlying arrays.
Recording projections can quickly fill up the RAM…
\tau \frac{dx(t)}{dt} + x(t) = \sum_\text{input} W^\text{IN} \, r^\text{IN}(t) + g \, \sum_\text{rec} W^\text{REC} \, r(t) + \xi(t)
r(t) = \tanh(x(t))
Jaeger (2001) The “echo state” approach to analysing and training recurrent neural networks. Jacobs Universität Bremen.
Synapses can also implement plasticity rules that will be evaluated after each neural update.
Example of the Intrator & Cooper BCM learning rule:
\Delta w = \eta \, r^\text{pre} \, r^\text{post} \, (r^\text{post} - \mathbb{E}[(r^\text{post})^2])
IBCM = ann.Synapse(
parameters = dict(eta = 0.01, tau = 2000.0),
equations = [
ann.Variable(
'tau * dtheta/dt + theta = post.r^2',
locality='semiglobal', method='exponential'),
ann.Variable('dw/dt = eta * post.r * (post.r - theta) * pre.r', min=0.0),
],
psp = "w * pre.r"
)
Each synapse can access pre- and post-synaptic variables with pre.
and post.
.
The semiglobal
locality allows to do computations only once per post-synaptic neuron.
psp
optionally defines what will be summed by the post-synaptic neuron (e.g. psp = "w * log(pre.r)"
).
Network.connect()
:or for a single post-synaptic neuron:
Variant of the BCM (Bienenstock, Cooper, Munro, 1982) learning rule.
The LTP and LTD depend on post-synaptic activity: homeostasis.
\Delta w = \eta \, r^\text{pre} \, r^\text{post} \, (r^\text{post} - \mathbb{E}[(r^\text{post})^2])
Intrator and Cooper (1992) Objective function formulation of the BCM theory of visual cortical plasticity: Statistical connections, stability conditions. Neural Networks, 5(1).
e_{i, j}(t) = e_{i, j}(t-1) + (r_i (t) \, x_j(t))^3
\Delta w_{i, j} = - \eta \, e_{i, j}(T) \, (R - R_\text{mean})
Miconi (2017) Biologically plausible learning in recurrent neural networks reproduces neural dynamics observed during cognitive tasks. eLife. doi:10.7554/eLife.20899
Spiking neurons must also define two additional fields:
spike
: condition for emitting a spike.
reset
: what happens after a spike is emitted (at the start of the refractory period).
A refractory period in ms can also be specified.
C \, \frac{d v(t)}{dt} = - g_L \, (v(t) - V_L) + I(t)
\text{if} \; v(t) > V_T \; \text{emit a spike and reset.}
\tau \cdot \frac{dv (t)}{dt} = E_l - v(t) + g_\text{exc} (t) \, (E_\text{exc} - v(t)) + g_\text{inh} (t) \, (E_\text{inh} - v(t)) + I(t)
Brette and Gerstner (2005) Adaptive Exponential Integrate-and-Fire Model as an Effective Description of Neuronal Activity. Journal of Neurophysiology 94.
g_target
(e.g. g_exc
or g_inh
, depending on the projection).LIF = ann.Neuron(
parameters = {'tau': 20.0},
equations = [
'tau * dv/dt + v = g_exc',
],
spike = "v >= 1",
reset = "v = 0"
)
Each spike increments instantaneously g_target
from the synaptic efficiency w
of the corresponding synapse.
This can be changed through the pre_spike
argument of the synapse model:
LIF = ann.Neuron(
parameters = {'tau'=20.0},
equations = [
'tau * dv/dt + v = g_exc',
'tau_exc * dg_exc/dt = - g_exc',
],
spike = "v >= 1",
reset = "v = 0"
)
equations = [
# Membrane potential
'tau*dv/dt = (E_L- v) + g_a + g_b + alpha_c',
# Exponentially decreasing
ann.Variable(
'tau_b * dg_b/dt = -g_b',
method='exponential'
),
# Alpha-shaped
ann.Variable(
'tau_c * dg_c/dt = -g_c',
method='exponential'
),
ann.Variable('''
tau_c * dalpha_c/dt =
exp((tau_c - dt/2.0)/tau_c) * g_c
- alpha_c''',
method='exponential'
),
],
\tau \, \frac{dv (t)}{dt} = E_l - v(t) + g_\text{exc} (t) \, (E_\text{exc} - v(t)) + g_\text{inh} (t) \, (E_\text{inh} - v(t)) + I(t)
Vogels and Abbott (2005) Signal propagation and logic gating in networks of integrate-and-fire neurons. The Journal of neuroscience 25.
Spiking synapses can define a pre_spike
field, defining what happens when a pre-synaptic spike arrives at the synapse.
g_target
is an alias for the corresponding post-synaptic conductance: it will be replaced by g_exc
or g_inh
depending on how the synapse is used.
By default, a pre-synaptic spike increments the post-synaptic conductance from w
, but this can be changed.
STP = ann.Synapse(
parameters = dict(
tau_rec = 100.0,
tau_facil = 0.01,
U = 0.5
),
equations = [
ann.Variable('dx/dt = (1 - x)/tau_rec', init = 1.0, method='event-driven'),
ann.Variable('du/dt = (U - u)/tau_facil', init = 0.5, method='event-driven'),
],
pre_spike="""
g_target += w * u * x
x *= (1 - u)
u += U * (1 - u)
"""
)
Tsodyks and Markram (1997) The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability. PNAS 94.
post_spike
similarly defines what happens when a post-synaptic spike is emitted.STDP = ann.Synapse(
parameters = dict(
tau_plus = 20.0, tau_minus = 20.0,
A_plus = 0.01, A_minus = 0.01,
w_min = 0.0, w_max = 1.0,
),
equations = [
ann.Variable('tau_plus * dx/dt = -x', method='event-driven'), # pre-synaptic trace
ann.Variable('tau_minus * dy/dt = -y', method='event-driven'), # post-synaptic trace
],
pre_spike="""
g_target += w
x += A_plus * w_max
w = clip(w + y, w_min , w_max)
""",
post_spike="""
y -= A_minus * w_max
w = clip(w + x, w_min , w_max)
""")
\begin{cases} \tau^+ \, \dfrac{d x(t)}{dt} = - x(t) \\ \\ \tau^- \, \dfrac{d y(t)}{dt} = - y(t) \\ \end{cases}
Bi and Poo (2001) Synaptic Modification by Correlated Activity: Hebb’s Postulate Revisited. Annual Review of Neuroscience 24.
Standard populations (SpikeSourceArray
, TimedArray
, PoissonPopulation
, HomogeneousCorrelatedSpikeTrains
), OpenCV bindings.
Standard neurons:
Standard synapses:
Parallel simulations with parallel_run
.
Convolutional and pooling layers.
Hybrid rate-coded / spiking networks.
Structural plasticity.
BOLD monitors.
Tensorboard visualization.
ANN (keras) to SNN conversion tool.