Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Brian 2: repeated execution of a standalone simulation with changed parameters without recompilation
import os
import time
################################################################################
# Helper functions that expose some internal Brian machinery
def initialize_parameter(variableview, value):
variable = variableview.variable
array_name = device.get_array_name(variable)
static_array_name = device.static_array(array_name, value)
device.main_queue.append(('set_by_array', (array_name,
static_array_name,
False)))
return static_array_name
def set_parameter_value(identifier, value):
np.atleast_1d(value).tofile(os.path.join(device.project_dir,
'static_arrays',
identifier))
def run_again():
device.run(device.project_dir, with_output=False, run_args=[])
from brian2 import *
set_device('cpp_standalone', directory='parallel', clean=False)
num_neurons = 500
duration = 2*second
# Parameters
area = 20000*umetre**2
Cm = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
# The model
eqs = Equations('''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
(exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
(exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
(exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
''')
# Threshold and refractoriness are only used for spike counting
group = NeuronGroup(num_neurons, eqs,
threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler', name='neurongroup')
group.v = El
param_I = initialize_parameter(group.I, np.zeros(num_neurons))
monitor = SpikeMonitor(group, name='spikemonitor')
run(duration)
# Run it again
start = time.time()
set_parameter_value(param_I, 0.7*nA *np.arange(num_neurons) / num_neurons)
run_again()
print('Took: {:.3f}s'.format(time.time() - start))
import os
from brian2 import *
################################################################################
# Helper functions that expose some internal Brian machinery
def initialize_parameter(variableview, value):
variable = variableview.variable
array_name = device.get_array_name(variable)
static_array_name = device.static_array(array_name, value)
device.main_queue.append(('set_by_array', (array_name,
static_array_name,
False)))
return static_array_name
def set_parameter_value(identifier, value):
np.atleast_1d(value).tofile(os.path.join(device.project_dir,
'static_arrays',
identifier))
def run_again():
device.run(device.project_dir, with_output=False, run_args=[])
################################################################################
# Begin of Brian script
# Adapted from: https://brian2.readthedocs.io/en/stable/examples/IF_curve_Hodgkin_Huxley.html
################################################################################
set_device('cpp_standalone', directory='sequential')
n_neurons = 500
duration = 2*second
# Parameters
area = 20000*umetre**2
Cm = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
# The model
eqs = Equations('''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
(exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
(exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
(exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
''')
# Threshold and refractoriness are only used for spike counting
group = NeuronGroup(1, eqs,
threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler')
group.v = El
I_param = initialize_parameter(group.I, 0*nA) # <-- set up I as changeable variable
monitor = SpikeMonitor(group)
import time
run(duration)
# Store results from first run
counts = np.full(n_neurons, np.nan)
counts[0] = monitor.count[0]
start = time.time()
# Sequential run
I_values = 0.7*nA *np.arange(n_neurons) / n_neurons
for idx, I in enumerate(I_values[1:]):
set_parameter_value(I_param, I) # <-- change the value of I
run_again()
counts[idx + 1] = monitor.count[0]
print('Took: {:.3f}s'.format(time.time() - start))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.