Changing a variable while looping for a multicompartment model

Support
1

Description of problem

Hello, I’m new to Brian2 but hope to use it to simulate simple multicompartment neuron models. To get started, I’m working with a somatic compartment where I want to parametrically explore the impact of changing the diameter of the somatic compartment. I’m starting by using a simple loop to change the diameter iteratively. Although the loop appears to be running, the diameter of the spatial neuron does not change during the simulation. I’ve looked at the documentation and the youtube videos on setting up simulations, and optimizing by vectorization when using NeuronGroup but haven’t really been able to solve this. I’m guessing/hoping that this is a fairly simple error and that someone out there with more experience with this simulation environment can help me here.

Minimal code to reproduce problem

defaultclock.dt = 0.025*ms
start_scope()

#Equations for channel definitions removed to make a more concise posting

eqs += eqs_leak + eqs_na + eqs_kht + eqs_klt #add other channels when ready

# Spatial Neuron
diameters = [30,50] 
morpho = Soma(diameter=5*1*um)  #initial definition with very small diameter...to quickly note if the #redefinition during the looping is working.
neuron = SpatialNeuron(morphology=morpho, model=eqs,Cm=1*uF/cm**2, Ri=70*ohm*cm, threshold='v > -40*mV', method='exponential_euler')                      
neuron[0].gl = 0.03e-3*siemens/cm**2 
neuron[0].gnabar = 20e-3*siemens/cm**2
neuron[0].gkhtbar = 2.8e-3*siemens/cm**2
neuron[0].gkltbar = 0*siemens/cm**2
mon_s = StateMonitor(neuron.main, ['v', 'I', 'Im'], record=True)
mon_spike = SpikeMonitor(neuron.main)

store() #this will save the above definition?
for diam in diameters: #start loop to change diameter in the SpatialNeuron object
  restore()
  morpho = Soma(diameter=diam*1*um)
  neuron.morphology=morpho #I thought this would update the diameter.
  print(neuron.morphology)
  #neuron = SpatialNeuron(morphology=morpho, model=eqs,Cm=1*uF/cm**2, Ri=70*ohm*cm, 
 #threshold='v > -40*mV', method='exponential_euler') #this produces a magic error
  ##Begin the simulation
  neuron.I=0*pA
  run(500*ms)
  neuron[0].I = 50*pA
  run(500*ms)
  neuron[0].I = 0*amp
  run(100*ms, report='text') 
  ##Plot the Results
  figure(1)
  subplot(211) 
  plot(mon_s.t/ms, mon_s.v[0]/mV, 'b') 
  axis([0, 20000, -80, 50]) 
  xlabel('Time (ms)') 
  ylabel('v (mV)') 
  title('soma')
  subplot(212)
  axis([0, 20000, 0, 200]) 
  plot(mon_s.t/ms, mon_s.I[0]/pA, 'b') 
  xlabel('Time (ms)')
  ylabel('Membrane potential (mV)')
  legend(['{} um'.format(d) for d in diameters])

What you have aready tried

  1. I tried moving repeating the SpatialNeuron definition after updating neuron.morphology definition within the loop but it produced a magic error
    #neuron = SpatialNeuron(morphology=morpho, model=eqs,Cm=1uF/cm**2, Ri=70ohmcm, threshold='v > -40mV’, method=‘exponential_euler’) #this produces a magic error
  2. I tried defining the SpatialNeuron in functional form with diam as an input variable. Not confident I did that correctly but it also produced a ‘magic error’
  3. I did not try the Neuron Group to vectorize the cell diameter because I’m using this simple case to learn how to run simulations using SpatialNeuron.

Expected output (if relevant)

I expected the current threshold and spiking pattern to scale in response to the two different diameter values.

Actual output (if relevant)

The response is that of the initial setting for diameter. So, it looks like the morphology parameter is not updating as we iterate through the loop.

Full traceback of error (if relevant)

2

Hi @radhak . After creating a SpatialNeuron object, its morphology cannot change anymore – we should indeed make this clearer, e.g. switching out the morphology as you did should raise an error. So the correct way to do this would be to recreate a new SpatialNeuron each time. You shouldn’t use the store and restore methods in this case, they are only meant to be used when you run the same set of objects for another time (e.g. with different inputs), but in your case you are creating new objects, so there is no previous state to restore. Also, you need to recreate the StateMonitor and SpikeMonitor, since they were still attached to the SpatialNeuron object you previously created. So in sum, your loop should look somewhat like this:

diameters = [30,50]

for diam in diameters:
  morpho = Soma(diameter=diam*um)
  neuron = SpatialNeuron(morphology=morpho, model=eqs, Cm=1 * uF / cm ** 2,
                         Ri=70 * ohm * cm, threshold='v > -40*mV',
                         method='exponential_euler')
  neuron[0].gl = 0.03e-3 * siemens / cm ** 2
  neuron.v = El
  mon_s = StateMonitor(neuron.main, ['v', 'I', 'Im'], record=True)
  mon_spike = SpikeMonitor(neuron.main)
  ##Begin the simulation
  neuron.I=0*pA
  ...

Now, this is the most flexible approach, e.g. you can change the morphology completely (add compartments, etc.) between runs. It is quite inefficient, though, vectorizing it with a single NeuronGroup object would be much faster – roughly speaking, running it for 100 different diameters in a single NeuronGroup will take about the same time as running one iteration of your loop. Of course, for more complex morphologies this won’t be straightforward anymore. If you are planning to only simulate simple morphologies (a few compartments), you might want to have a look at the dendrify package which will make this process much easier.