add_drives_from_params vs add_evoked_drive

[mkhalil8]

Hi,

I am trying to use the same parameters to run a simulation with two different methods but not getting the same result.

The first method is where I add the drives from params:

net = hnn_core.calcium_model(params, add_drives_from_params = True)

The second method I add the drives manually (using same parameters from the param file:

net = hnn_core.calcium_model(params, add_drives_from_params = False)
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net1.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal', event_seed = -14, synaptic_delays=synaptic_delays_prox)

I am not getting the same results from both, the dipole is different between both. I tracked the add_drives_from_params function and I found that it substract 18 from the event seed so I adjusted that. The rest appears to be the same, not sure where the problem

[ntolley]

@mkhalil8 do you mind pointing to where 18 is subtracted from event_seed, as well as email the param file to nicholas_tolley@brown.edu ?

I suspect there is a deeper seeding discrepancy, but there may be a more obvious explanation that @jasmainak or @rythorpe remember from when we developed the drives API. I do recall that there were some awkward workarounds to deal with the old param seeding/format.

[mkhalil8]

Hi Nicholas and Happy New Year :)

The 18 is subtracted from event_seed at :
drives.py

def _extract_drive_specs_from_hnn_params(params, cellname_list):

223     drive['event_seed'] = par['prng_seedcore'] - 18

[jasmainak]

The add_drives_from_params should be deprecated. Why do we expose it in the calcium model @ntolley ? We've gone down this rabbit hole many times and it's not worth it maintaining these two versions. Zen of Python says:

There should be one-- and preferably only one --obvious way to do it.

Or is this for IO?

I think @rythorpe at some point figured that subtracting 18 maintained backward compatibility with previous versions in GUI

[ntolley]

It was leftover from when we were originally planning to replace the model in all examples with calcium_model:

# Remove params argument after updating examples
# (only relevant for Jones 2009 model)
def calcium_model(params=None, add_drives_from_params=False):

As of right now the main reason to keep it is compatibility with HNN-GUI. It is still the standard way people are introduced to HNN and is important to show that they produce identical results. Until we have a GUI running hnn-core I think we will be stuck with param files for a bit longer.

[rythorpe]

 I think I came across the same issue @mkhalil8 has encountered when I was putting together the hnn_core gamma tutorial but forgot to follow up. I found that the order in which drives are added is important when trying to match the results of automatically added drives from a param file. The -18 applied to the seed core is important for maintaining continuity with HNN-GUI, but I can’t recall the details off the top of my head. I’m away from my computer but will follow up sometime in the next few days.

…

On Jan 2, 2022, at 10:45 AM, Nick Tolley ***@***.***> wrote:  It was leftover from when we were originally planning to replace the model in all examples with calcium_model: # Remove params argument after updating examples # (only relevant for Jones 2009 model) def calcium_model(params=None, add_drives_from_params=False): As of right now the main reason to keep it is compatibility with HNN-GUI. It is still the standard way people are introduced to HNN and is important to show that they produce identical results. Until we have a GUI running hnn-core I think we will be stuck with param files for a bit longer. — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you were mentioned.

[ntolley]

@mkhalil8 is the parameter file you send me supposed to only add a single proximal drive, or was that just a small example from the code above?

Using the .param you sent me, there are quite a large number of drives created. Try inspecting net.external_drives to see what I mean:

params = read_params('L_Contra.param')
net_params = calcium_model(params, add_drives_from_params = True)
print(net_params.external_drives.keys())

dict_keys(['bursty1', 'bursty2', 'evdist1', 'evprox1', 'evprox2', 'extgauss', 'extpois'])

[mkhalil8]

@ntolley Yes it does add these drives but the weights of all these drives (except the evoked) are Zeros, so I guess they don't end up affecting the network.

[mkhalil8]

@ntolley when comparing the two methods of adding drives I get this result (one trial simulation for each method with the same seed)

[ntolley]

Do you mind sharing the full code you're using for manually adding the evoked drives? It seems like there should be a distal drive as well since they have non-zero weights.

[mkhalil8]

Yes there are p1, d1, p2

#code for manually adding drives
net_ev = hnn_core.calcium_model(params, add_drives_from_params=False)
net1 = hnn_core.calcium_model(params, add_drives_from_params=False)
    # p1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722, 'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247, 'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net1.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal', event_seed =4, synaptic_delays=synaptic_delays_prox)

    # d1 drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.25807}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net1.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal', event_seed = 4, synaptic_delays = synaptic_delays_d1)
  
    # p2 drive
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454, 'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 , 'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net1.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal', event_seed = 4)
dpls_ev = simulate_dipole(net1, tstop=250, n_trials=1)
dpl_ev = dpls_ev[0]
dpl_ev.smooth(30).scale(1500)

[ntolley]

@mkhalil8 I've figured out that manually adding the blank bursty drives that are added by the .param file are necessary. This because the state of the RNG is altered depending on what was previously added. The end result is that that spike times of the evoked drives match exactly. Unfortunately this hasn't completely resolved the issue, as there seems to be some drift between the waveforms later in the simulation:

#code for manually adding drives
net_manual = calcium_model(params, add_drives_from_params=False)

# Add empty bursty drives

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_manual.add_bursty_drive(
    'bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,
    spike_isi=10, n_drive_cells=10, location=location,
    weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_manual.add_bursty_drive(
    'bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,
    spike_isi=10, n_drive_cells=10, location=location,
    weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

    # d1 drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.25807}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_manual.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal', synaptic_delays = synaptic_delays_d1, event_seed=-14)

    # p1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722, 'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247, 'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_manual.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal', synaptic_delays=synaptic_delays_prox, event_seed=-14)
  
    # p2 drive
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454, 'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 , 'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_manual.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal', event_seed=-14)

It's a bit fishy that the waveforms don't match exactly despite identical driving inputs. This may be an artifact of running two simulations back to back, and global variables aren't being cleared properly. I'll discuss with the dev team this week and report back.

[mkhalil8]

@rythorpe Thanks Ryan, in what manner do the order of adding drives matter? shall I in the case of ERP add drives with the following order:
proximal 1, distal 1, proximal 2

[mkhalil8]

@ntolley Thanks a lot Nick, I believe this answers the question for now, will wait for any updates on the waveforms!!

[mkhalil8]

@ntolley So we don't need to add the extgauss or extpois, right?

Do you mind also pointing out where in the code does adding the burst drive updates the RNG? just to make sure I understood :)

[jasmainak]

do we need to expose legacy_mode here? See:

hnn-core/hnn_core/network.py

Lines 842 to 844 in 1c2662c

	# allows tests must match HNN GUI output by preserving original
	# gid assignment convention
	target_populations = list(self.cell_types.keys())

[jasmainak]

Btw @ntolley @mkhalil8 for future, if you can post the param files on github/gist, that would be great. Having as much of the communication in public as possible has the benefit that anyone can help ... plus it's helpful to newcomers in the future.

[jasmainak]

As of right now the main reason to keep it is compatibility with HNN-GUI. It is still the standard way people are introduced to HNN and is important to show that they produce identical results. Until we have a GUI running hnn-core I think we will be stuck with param files for a bit longer.

But is calcium model implemented in the GUI?

[mkhalil8]

@jasmainak Yes sure will share the param file later.
For the GUI, I think yes. There was mention of an L5Pyr.py file that is added two the GUI files to be able to use the calcium model

[rythorpe]

@ntolley have you explicitly verified that the spike times are exactly the same when the drives are automatically vs manually added? I suspect that somewhere there is a very minor difference....

[ntolley]

@rythorpe they actually match exactly using a == check on the entire drive element! It's rather strange :/

[mkhalil8]

@jasmainak the param file I used from Kohl simulation:

https://github.com/kohl-carmen/HNN-AEF/blob/main/HNN_Parameters/L_Contra.param

[mkhalil8]

@ntolley Nick how is the order of adding the evoked drives makes a difference in the dipole? I shuffled their order and got different results.
Also, if I am adjusting one of the three drives, let's say d1 to be the only one synchronized, would that interact with their order?

[ntolley]

Sorry the delay @mkhalil8! I've been swamped the past few days but have some time this week to dig into this a bit deeper.

At the end of the day the difference in the dipole all comes from seeding. When you instantiate a drive, a random seed is used with the random number generator to produce the spike times (distributed according to the type of drive). A different seed produces different spike times.

I didn't write this part of the code base so I need to look into this more carefully before pointing you to the correct sections. In any case, the desired behavior is that every time a drive is added, the seed for the random number generator is incremented. This to avoid weird behavior like two proximal and distal drives having completely coincident spike times.

[ntolley]

Also, if I am adjusting one of the three drives, let's say d1 to be the only one synchronized, would that interact with their order?
I don't believe it should have an effect but will double check.

This brings up a more general question @jasmainak , when we say our code is reproducible should that include the order in which steps are performed? Or should it be reproducible when the calls used to construct the network occur in different order.

[jasmainak]

yes our code should be reproducible no matter what order you add the drives because it's the same network

@mkhalil8 I'd need a small script that demos this to investigate further

[jasmainak]

@mkhalil8 I can't import that param file. I tried:

from hnn_core import read_params

params = read_params('L_contra.param')

Did you make any modifications to it?

[ntolley]

@jasmainak I think the text export messed up, you have to remove the blank lines

[jasmainak]

okay I created a script that tries different permutations of adding the drives. It doesn't seem to make a difference. Can you show me how to break this script @ntolley @mkhalil8 ?

[mkhalil8]

@jasmainak @ntolley I was originally running batch simulations varying the synchronization of the drives but got mixed results, So I used this code below to experiment and figure out where the error is coming from. In the code below I instantiated new network with different name for each simulation(just to make sure this is not the culprit for the error) and changed the order of the drives while exploring the effect of differential synchronization between drives. I will post down two sets of code each with two different orders of drives added to the network. The first one: d1, p1,p2. The second one is : p1,p2,d1.

The first script: adding d1, p1,p2:

## simulating contrl
import os.path as op
from itertools import product
import numpy as np
from joblib import Parallel, delayed
import hnn_core
from hnn_core import simulate_dipole, jones_2009_model, average_dipoles
from hnn_core import read_params
from hnn_core.viz import  plot_dipole
import json
from hnn_core import JoblibBackend


hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14)
# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
 
dpl = simulate_dipole(net, tstop=250, n_trials= 1)
dpl_cntrl= dpl[0]
dpl_cntrl.smooth(30).scale(1500)

##  sync d1

hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net_d = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_d.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_d.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_d.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific)
# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_d.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_d.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
 
dpls_sync_d1 = simulate_dipole(net_d, tstop=250, n_trials= 1)
dpl_sync_d1= dpls_sync_d1[0]
dpl_sync_d1.smooth(30).scale(1500)



# n_drive_cells = n_drive_cells, cell_specific = cell_specific

## sync p1, d1


hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net_p = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_p.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_p.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_p.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific)
# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_p.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_p.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
 
dpls_sync_p1_d1 = simulate_dipole(net_p, tstop=250, n_trials= 1)
dpl_sync_p1_d1= dpls_sync_p1_d1[0]
dpl_sync_p1_d1.smooth(30).scale(1500)

##all_drives synchronized


hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net_all = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_all.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_all.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_all.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific)
# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_all.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_all.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14,n_drive_cells = n_drive_cells, cell_specific = cell_specific)
 
dpls_all = simulate_dipole(net_all, tstop = 250, n_trials=1)
dpl_sync_all = dpls_all[0]
dpl_sync_all.smooth(30).scale(1500)

# plotting after adding them

import matplotlib.pyplot as plt
plt.ion()
fig, axes = plt.subplots(1, 1, sharex=True, figsize=(6, 6),constrained_layout=True)
axes.plot(dpl_cntrl.times,dpl_cntrl.data['agg'], label ='cntrl')
axes.plot(dpl_sync_d1.times,dpl_sync_d1.data['agg'], label ='d1_sync')
axes.plot(dpl_sync_p1_d1.times,dpl_sync_p1_d1.data['agg'], label ='p1_d1_sync')
axes.plot(dpl_sync_all.times, dpl_sync_all.data['agg'], label= 'all_sync')
plt.legend()
axes.legend()
axes.grid()
axes.set_xlabel("time (ms)")
axes.set_ylabel("nAm x scaling factor")
fig.savefig('/storage/work/mzk5905/Batch hnn_core/data/experimenting_kohl/all_manual_shuffle2_d1p1p2_cntrlunsync_d1sync_p1d1sync_allsync.png')

[mkhalil8]

The second script: adding p1,p2,d1 (while varying synchronization the same way as above):

import os.path as op
from itertools import product
import numpy as np
from joblib import Parallel, delayed
import hnn_core
from hnn_core import simulate_dipole, jones_2009_model, average_dipoles
from hnn_core import read_params
from hnn_core.viz import  plot_dipole
import json
from hnn_core import JoblibBackend


hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14) 

dpl = simulate_dipole(net, tstop=250, n_trials= 1)
dpl_cntrl= dpl[0]
dpl_cntrl.smooth(30).scale(1500)

##  sync d1

hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net_d = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_d.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_d.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_d.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_d.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_d.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific) 

dpls_sync_d1 = simulate_dipole(net_d, tstop=250, n_trials= 1)
dpl_sync_d1= dpls_sync_d1[0]
dpl_sync_d1.smooth(30).scale(1500)



# n_drive_cells = n_drive_cells, cell_specific = cell_specific

## sync p1, d1


hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net_p = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_p.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_p.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_p.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_p.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_p.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific) 

dpls_sync_p1_d1 = simulate_dipole(net_p, tstop=250, n_trials= 1)
dpl_sync_p1_d1= dpls_sync_p1_d1[0]
dpl_sync_p1_d1.smooth(30).scale(1500)

##all_drives synchronized


hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net_all = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_all.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_all.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_all.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_all.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14,n_drive_cells = n_drive_cells, cell_specific = cell_specific)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_all.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific) 

dpls_all = simulate_dipole(net_all, tstop = 250, n_trials=1)
dpl_sync_all = dpls_all[0]
dpl_sync_all.smooth(30).scale(1500)

# plotting after adding them

import matplotlib.pyplot as plt
plt.ion()
fig, axes = plt.subplots(1, 1, sharex=True, figsize=(6, 6),constrained_layout=True)
axes.plot(dpl_cntrl.times,dpl_cntrl.data['agg'], label ='cntrl')
axes.plot(dpl_sync_d1.times,dpl_sync_d1.data['agg'], label ='d1_sync')
axes.plot(dpl_sync_p1_d1.times,dpl_sync_p1_d1.data['agg'], label ='p1_d1_sync')
axes.plot(dpl_sync_all.times, dpl_sync_all.data['agg'], label= 'all_sync')
plt.legend()
axes.legend()
axes.grid()
axes.set_xlabel("time (ms)")
axes.set_ylabel("nAm x scaling factor")
fig.savefig('/storage/work/mzk5905/Batch hnn_core/data/experimenting_kohl/all_manual_shuffle2_p1p2d1_cntrlunsync_d1sync_p1d1sync_allsync.png')

[mkhalil8]

@jasmainak @ntolley I have noticed also that syncrhonizing d1 affect p50 making it stronger, although d1 takes effect around 80 msec in time !!!. I tried narrowing the sigma but it still gave the same effect

[jasmainak]

@mkhalil8 would you mind simplifying these scripts a bit? It's a little too much to look for the problem in this. Ideally a script that plots two dipoles comparing the original and permuted addition order. You can also use loops/functions or just show for one case. Scripts shorter than 30 lines is preferable but I can live with something less than 50 lines ;-)

[mkhalil8]

@jasmainak sorry Mainak for the long script. This is a shorter one of only two simulations. The only difference between both is the order of the drives (d1 is synchronized in both). it is 86 lines long (sorry the best I could do)
The first one is d1,p1,p2
The second one is p1,p2,d1

from itertools import product
import numpy as np
from joblib import Parallel, delayed
import hnn_core
from hnn_core import simulate_dipole, jones_2009_model, average_dipoles
from hnn_core import read_params
from hnn_core.viz import  plot_dipole
import json
from hnn_core import JoblibBackend
### order of drives d1,p1,p2
hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net_d = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_d.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_d.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_d.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific)
# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_d.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_d.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
 
dpls_sync_d1 = simulate_dipole(net_d, tstop=250, n_trials= 1)
dpl_sync_d1= dpls_sync_d1[0]
dpl_sync_d1.smooth(30).scale(1500)

##################################################################
##  p1p2d1
net_s = hnn_core.calcium_model(params,add_drives_from_params= False)

net_s.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_s.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_s.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)

# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_s.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_s.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific) 

dpls_sync_d1_s = simulate_dipole(net_s, tstop=250, n_trials= 1)
dpl_sync_d1_s= dpls_sync_d1_s[0]
dpl_sync_d1_s.smooth(30).scale(1500)

import matplotlib.pyplot as plt
plt.ion()
fig, axes = plt.subplots(1, 1, sharex=True, figsize=(6, 6),constrained_layout=True)
axes.plot(dpl_sync_d1.times,dpl_sync_d1.data['agg'], label ='d1_p1p2')
axes.plot(dpl_sync_d1_s.times,dpl_sync_d1_s.data['agg'], label ='p1p2d1')
plt.legend()
axes.legend()
axes.grid()
axes.set_xlabel("time (ms)")
axes.set_ylabel("nAm x scaling factor")
fig.savefig('/storage/work/mzk5905/Batch hnn_core/data/experimenting_kohl/all_manual_simple__d1syncp1p2__p1p2d1sync.png')

[jasmainak]

Thanks! Okay I confirm I can reproduce and it's very strange. Will investigate and get back

[jasmainak]

sorry for the delay in responding here @mkhalil8 . I have been "vacation busy".

So I figured out the issue. Let me give you a bit of historical context and we can decide what's the best solution moving forward. The spiking events of the external drives were made such that every cell ID has a separate seed for the events. For example, if you provided event_seed=2 and net.gid_ranges['evprox1'] = range(35, 70) then you will get seeds from 35 + 2 = 37 to 70 + 2 = 72. I was told in a group meeting that the reason for this is to prevent surprises when debugging or when extending the tstop in the GUI. You want the same spike pattern of external drives to show even if you changed the tstop.

A simple solution that comes to mind is to sort net.external_drives alphabetically so the gids are assigned similarly. But this won't solve the problem if you name the drives differently in the two networks. Ultimately the gid assignment affects the events and this causes the issue. Another simple solution is to call event_seed as event_seedcore to make it explicit what it does and prevent such confusions in the future. I think we should discuss this with @stephanie-r-jones and collectively brain storm what could be the solution. Please feel free to share further thoughts or ideas to solve this @rythorpe @ntolley @mkhalil8

[jasmainak]

another idea is to use gid_offset instead of gid for the seeds ... where gid_offset is gid - gid_range[0]

[mkhalil8]

@jasmainak I hope you had nice vacation and sorry for my delay in response also :)

So for now, should keeping the order of the drives similar among different networks mitigate the problem? till we work on the solution?
On the other side, my understanding for synchronizing a drive is that it makes the drive one cell and it fires to the network cells that the drive is connected to. Actually my primary question that got me into the order of drives was while working with synchronization of one drive vs keeping others unsynchronized, which got me some unexpected results.
What I did specifically was to make one of the drives (the distal one) synchronized. The distal drive fires at approximately 82msec and it affected the strength of the N100 as expected (and showed in one of your tutorials), but some how it is strengthening also the P50 (as you can see in the figure below).
So my second question is: Would synchronizing one of the drives affects one of the other two based on order? can that be the explanation?

import os.path as op
from itertools import product
import numpy as np
from joblib import Parallel, delayed
import hnn_core
from hnn_core import simulate_dipole, jones_2009_model, average_dipoles
from hnn_core import read_params
from hnn_core.viz import  plot_dipole
import json
from hnn_core import JoblibBackend
### network 1 unsync
hnn_core_root = op.dirname(hnn_core.__file__)
params_fname = op.join(hnn_core_root, 'param', 'L_Contra.param')
params = read_params(params_fname)
net_d = hnn_core.calcium_model(params,add_drives_from_params= False)

n_drive_cells = 1
cell_specific = False

location = 'proximal'
burst_std = 20
weights_ampa_p = {'L2_pyramidal': 0, 'L5_pyramidal': 0}
syn_delays_p = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.}

net_d.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_d.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_d.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_d.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14)
# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_d.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
 
dpls_sync_d1 = simulate_dipole(net_d, tstop=250, n_trials= 1)
dpl_sync_d1= dpls_sync_d1[0]
dpl_sync_d1.smooth(30).scale(1500)

##################################################################
##  network 2 d1 sync
net_s = hnn_core.calcium_model(params,add_drives_from_params= False)

net_s.add_bursty_drive('bursty1', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2, spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

net_s.add_bursty_drive('bursty2', tstart=50., burst_rate=10, burst_std=burst_std, numspikes=2,spike_isi=10, n_drive_cells=10, location=location,weights_ampa=weights_ampa_p, synaptic_delays=syn_delays_p, event_seed=-14)

# Proximal 1 drive
weights_ampa_p1 = {'L2_basket': 0.997291, 'L2_pyramidal':0.990722,'L5_basket':0.614932, 'L5_pyramidal': 0.004153}
weights_nmda_p1 = {'L2_basket': 0.984337, 'L2_pyramidal':1.714247,'L5_basket':0.061868, 'L5_pyramidal': 0.010042}
synaptic_delays_prox = {'L2_basket': 0.1, 'L2_pyramidal': 0.1,'L5_basket': 1., 'L5_pyramidal': 1.}
net_s.add_evoked_drive('evprox1', mu=54.897936, sigma=5.401034, numspikes=1, weights_ampa=weights_ampa_p1, weights_nmda=weights_nmda_p1, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)
# Distal drive
weights_ampa_d1 = {'L2_basket': 0.624131, 'L2_pyramidal': 0.606619, 'L5_pyramidal':0.258}
weights_nmda_d1 = {'L2_basket': 0.95291, 'L2_pyramidal': 0.242383, 'L5_pyramidal': 0.156725}
synaptic_delays_d1 = {'L2_basket': 0.1, 'L2_pyramidal': 0.1, 'L5_pyramidal': 0.1}
net_s.add_evoked_drive('evdist1', mu=82.9915, sigma=13.208408, numspikes=1, weights_ampa=weights_ampa_d1,weights_nmda=weights_nmda_d1, location='distal',synaptic_delays=synaptic_delays_d1, event_seed=-14, n_drive_cells = n_drive_cells, cell_specific = cell_specific) 
# Second proximal evoked drive.
weights_ampa_p2 = {'L2_basket': 0.758537, 'L2_pyramidal': 0.854454,'L5_basket': 0.979846, 'L5_pyramidal': 0.012483}
weights_nmda_p2 = {'L2_basket': 0.851444, 'L2_pyramidal':0.067491 ,'L5_basket': 0.901834, 'L5_pyramidal': 0.003818}
net_s.add_evoked_drive('evprox2', mu=161.306837, sigma=19.843986, numspikes=1, weights_ampa=weights_ampa_p2,weights_nmda= weights_nmda_p2, location='proximal',synaptic_delays=synaptic_delays_prox, event_seed=-14)


dpls_sync_d1_s = simulate_dipole(net_s, tstop=250, n_trials= 1)
dpl_sync_d1_s= dpls_sync_d1_s[0]
dpl_sync_d1_s.smooth(30).scale(1500)

import matplotlib.pyplot as plt
plt.ion()
fig, axes = plt.subplots(1, 1, sharex=True, figsize=(6, 6),constrained_layout=True)
axes.plot(dpl_sync_d1.times,dpl_sync_d1.data['agg'], label ='no_drives_synchronized')
axes.plot(dpl_sync_d1_s.times,dpl_sync_d1_s.data['agg'], label ='d1_only_sync')
plt.legend()
axes.legend()
axes.grid()
axes.set_xlabel("time (ms)")
axes.set_ylabel("nAm x scaling factor")
fig.savefig(".........")

[ntolley]

@mkhalil8 I think the best way to ensure your results are not an artifact of drive order is to run multiple trials of each simulation `simulate_dipole(..., n_trials=>1). This will change the random seed for the spike times and allow you to see what parts of the simulation are consistent.

In this specific example, I think the reason you're seeing the P50 peak change is because its final amplitude is determined by the balance of proximal and distal drives. Excitatory distal inputs do indeed produce a more negative dipole (due to the direction of current flow in the apical dendrite). However, synchronizing distal drive means that this effect is exclusively observed at the time the single distal artificial cell fires.

The end result is that that early proximal drive has nothing to oppose it at t=50, producing the larger peak. When the distal drives are not synchronized, then distal inputs occurring earlier will pull the P50 peak down (this is speculation, definitely confirm this by looking at the spike times/histograms for the drives).

[ntolley]

For reference, here is what happens when I run your script for 5 trials:

The P50 peak is consistently higher (save for one simulation)when distal drives are synchronized so this is not an artifact of seeding. Additionally you'll notice the time of the negative peak is much more variable since the single time the distal drive fires is random:

You'll notice there is one simulation with an earlier distal spike, I'd put money on that being the one with a small P50 peak!

[mkhalil8]

@ntolley. Thanks Nick, yes it makes sense.I will run multiple trials as you suggested and look at the spiking pattern.

[jasmainak]

@ntolley @mkhalil8 would it make sense to document some "best practices" for running HNN simulations? I can see two points coming out of this discussion:

1. Always look at dipoles in conjunction with raster plots etc to avoid misinterpretation
2. Run multiple trials before drawing conclusions

Could one of you take the initiative in adding this?

[mkhalil8]

@jasmainak yes sure Mainak I would like to do that. which part of the documentation will this be part of? I can make it under title of adding drives automatically vs manually for example, shall I make a draft and send to you for editing?

[jasmainak]

@mkhalil8 try making a pull request following the contribution guidelines. Are you familiar with git workflow? Otherwise I'll update the contributing guideline with that information.

I would add this information in all the evoked example. You can use the "..warning" directive in sphinx if you want to be more emphatic.

[mkhalil8]

@jasmainak yes sure I am familiar with git workflow, will start working on that.

[mkhalil8]

@jasmainak which branch shall I push the modified plot_simulate_evoked.rst file to?

[jasmainak]

@mkhalil8 you should push to a branch on your fork (not main or master) and make a "pull request" from that branch

[mkhalil8]

@jasmainak what I did was making a branch called bestpractice then tried the following:
$ git push origin bestpractice
remote: Permission to jonescompneurolab/hnn-core.git denied to mkhalil8.
fatal: unable to access 'https://github.com/jonescompneurolab/hnn-core/': The requested URL returned error: 403

[mkhalil8]

Am I doing it the wrong way?

[jasmainak]

what you call origin is typically called upstream. You should rename your current origin to be upstream and add a new origin:

$ git remote rename origin upstream
$ git remote add origin https://github.com/mkhalil8/hnn-core/