Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Discrepancy in Synapse-Replay Simulation Outputs between Neurodamus and BlueCelluLab #75

Open
dhuruvapriyan opened this issue Sep 4, 2023 · 11 comments

Comments

@dhuruvapriyan
Copy link

Hi @anilbey,

I would like to report a discrepancy between the two codebases when using the synapse replay feature. I attempted to simulate synapses of a cell using SScx dataset, specifically focusing on the edge population labeled "S1_nonbarrel_neurons__S1_nonbarrel_neurons." The somatic voltage traces obtained from BlueCelluLab do not align with those from Neurodamus. Interestingly, when the release probability is set to 1, the traces do match, suggesting that the issue may lie in my random seed settings. I have attached my simulation_config, spikes files, and results for your review. Please let me know if you need additional details to help resolve this issue. Thank you in advance!

image

Matching Results for Use = 1:
image

My simulation_config.json:

{

  "run": {
    "dt": 0.025,
    "tstop": 2000,
    "random_seed": 12345
  },

  "conditions": {
    "extracellular_calcium": 1.05,
    "v_init": -80.0,
    "spike_location": "AIS",
    "mechanisms": {
        "ProbAMPANMDA_EMS": {
            "init_depleted": false,
            "minis_single_vesicle": false
        },
        "ProbGABAAB_EMS": {
            "init_depleted": false,
            "minis_single_vesicle": false
        }
    }
  },

  "target_simulator": "CORENEURON",
  "network": "./circuit_config.json",
  "node_sets_file": "./node_sets.json",
    "node_set": "my_circuit",
  "output": {
    "output_dir": "./reporting",
    "spikes_file": "spikes.h5"
  },
  
  "inputs":{ 
	"s1_syns_spikes": {
            "node_set": "hex_O1",
            "input_type": "spikes",
            "delay": 0.0,
            "duration": 2000.0,
            "module": "synapse_replay",
            "spike_file": "./input_spikes/s1_syn_spikes.dat",
            "source": "hex_O1"
            }
},
    "connection_overrides": [
        {
            "name": "init",
            "source": "hex_O1",
            "target": "hex_O1",
            "spont_minis":0.0,
            "weight": 1.0
        },
        {
            "name": "VPM_init",
            "source": "proj_Thalamocortical_VPM_Source",
            "target": "hex_O1",
            "spont_minis": 0.0,
            "weight": 1.0
        },
        {
            "name": "VPM-L5E_meanE",
            "source": "proj_Thalamocortical_VPM_Source",
            "target": "Layer5Excitatory",
            "spont_minis": 1.0,
            "weight": 0.733
        },
        {
            "name": "POm_init",
            "source": "proj_Thalamocortical_POM_Source",
            "target": "hex_O1",
            "spont_minis": 0.0,
            "weight": 1.0
        }
    ],
    "reports": {
        "soma_report": {
            "type": "compartment",
            "cells": "my_circuit",
            "variable_name": "v",
            "dt": 0.25,
            "start_time": 0.0,
            "end_time": 11000.0
        }
    }    

    
}

My spikes data:

/scatter
200 375
400 413
800 744
1000 3228
1400 3869
1800 12503

Thanks,
Dhuruva Priyan G M

@anilbey
Copy link
Contributor

anilbey commented Sep 5, 2023

Thanks a lot @dhuruvapriyan for reporting this. May I ask how you are changing the release probability? Do you modify the simulation config for it?

@dhuruvapriyan
Copy link
Author

Hi @anilbey,
Yes. To modify the release probability, I have added this snippet to the connection_overrides block of simulation_config.json.

"connection_overrides":[
       {
              "name": "Make_Use_1",
              "source": "hex_O1",
              "target": "hex_O1",
              "synapse_configure": "%s.Use = 1"
        }
]

@dhuruva1ca
Copy link

Hi @anilbey,

I believe the issue might be related to random number streams. This seems to be connected to #72.

Please let me know if you would like me to conduct validation experiments using both of these code bases.

Thanks,
@dhuruvapriyan

@anilbey
Copy link
Contributor

anilbey commented Sep 26, 2023

Hi again @dhuruvapriyan. It's great that you are advancing on the issue.

Please let me know if you would like me to conduct validation experiments using both of these code bases.

Oh that would be great. I am currently stuck on this below issue with the Zenodo circuit.

The simulation you are using - as well as the simulations provided on Zenodo - are setting the extracellular_calcium parameter in the config file.

E.g.

 "conditions": {
    "extracellular_calcium": 1.05,

The mechanisms provided in Zenodo under O1_mods.xz however do not contain the GluSynapse.mod hence there is no cao_CR_GluSynapse property in NEURON. The simulators need to set that property using the value in the simulation config.

How did you run BlueCellulab? Did you download the GluSynapse.mod elsewhere and compiled all the mechanisms together?

@anilbey
Copy link
Contributor

anilbey commented Sep 26, 2023

Btw. last Friday the docker image of Neurodamus is created. It is easier to install Neurodamus now.

Docker image: https://hub.docker.com/r/bluebrain/neurodamus
Instructions: https://github.com/BlueBrain/neurodamus/tree/main/docker

@dhuruvapriyan
Copy link
Author

Dear @anilbey,

Thank you for your prompt response regarding this issue.

We have recently rebuilt Neurodamus using the newly released Ecker's SSCx plasticity dataset, which can be found here: https://zenodo.org/record/8158472. Subsequently, we recompiled the "mods" present in the SSCx_plasticity dataset and incorporated them into BlueCelluLab.

I successfully simulated VPM fibers with GluSynapse using BlueCelluLab. However, I encountered certain challenges when operating BlueCelluLab with the SSCx_plasticity model. I shall create a new ticket to provide a detailed description of these challenges encountered during the simulation of GluSynapse with BlueCelluLab.

Thanks,
Dhuruva Priyan

@dhuruvapriyan
Copy link
Author

Hey @anilbey,

I'm reaching out to revisit the matter of the plasticity simulations, as they are critical to our current project. Previously, I successfully synchronized the CPre and CPost parameters in both BlueCelluLab and Neurodamus platforms. Yet, I've encountered a discrepancy: the synaptic voltage traces are misaligned when the USE parameter differs from 1. It appears that the 'use_scale_factor' might be improperly set to adjust the USE values. I'd appreciate your input on this observation. I am happy to assist you by running any necessary simulations as needed. Thanks in advance!

Regards,
Dhuruva Priyan

image

When Use = 1:
image

BlueCelluLab Code:

import pickle
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
%matplotlib inline
import bluecellulab
from bluecellulab import SSim

bluecellulab.neuron.load_mechanisms("./circuit_sonata/")
cfg_json = ".../simulation_config.json"
sim = SSim(cfg_json,print_cellstate=True)

post_gid = ('S1nonbarrel_neurons', 170923)
pre_gid = ('S1nonbarrel_neurons', 170982)

bluecellulab.neuron.h.tau_effca_GB_GluSynapse = 278.3177658387
bluecellulab.neuron.h.gamma_d_GB_GluSynapse =  101.5387594661
bluecellulab.neuron.h.gamma_p_GB_GluSynapse =  216.1841700668

pre_spike_times = np.array([ 817,  829,  884,  995,  996, 1000, 1080, 1135, 1335, 1403])

sim.instantiate_gids([post_gid], add_synapses=True,
                          pre_spike_trains={pre_gid: pre_spike_times}, add_projections = False, add_minis = False, intersect_pre_gids=[pre_gid])

t_vec = bluecellulab.neuron.h.Vector()
t_vec.record(bluecellulab.neuron.h._ref_t)
v_synapses = {}

syns = {}

for post_gid in sim.cells:
    if post_gid == post_gid:
        n_of_syns = len(sim.cells[post_gid].synapses)
        for syn_id in sim.cells[post_gid].synapses:
            syns[syn_id] = sim.cells[post_gid].synapses[syn_id].hsynapse
            if ("GluSynapse" in str(syns[syn_id])) and sim.cells[post_gid].synapses[syn_id].hsynapse.gmax0_AMPA:
                v_synapses[syn_id] = neuron.h.Vector()
                _ = v_synapses[syn_id].record(syns[syn_id]._ref_vsyn)

duration = 2500
sim.run(duration, cvode=False)

Simulation Config:

{
    "manifest": {
        "$CURRENT_DIR": "."
    },
    "run": {
        "dt": 0.025,
        "tstop": 2500.0,
        "random_seed": 12345,
        "using_cvode": 0
    },
    "conditions": {
        "extracellular_calcium": 2.0,
        "v_init": -80.0,
        "spike_location": "AIS",
        "mechanisms": {
            "ProbAMPANMDA_EMS": {
                "init_depleted": true,
                "minis_single_vesicle": true
            },
            "ProbGABAAB_EMS": {
                "init_depleted": true,
                "minis_single_vesicle": true
            },
            "GluSynapse": {
                "init_depleted": false,
                "minis_single_vesicle": false,
                "cao_CR": 2.0,
		"tau_effca_GB": 278.3177658387,
		"gamma_d_GB" : 101.5387594661,
		"gamma_p_GB" : 216.1841700668
            }
            }
    },
    "target_simulator": "CORENEURON",
    "network": "/project/ctb-emuller/datasets/SSCx_plasticity/O1_2023a_Ecker/circuit_config.json",
    "node_sets_file": "./node_sets.json",
    "node_set": "test_circuit",
    "output": {
        "output_dir": "./reporting",
        "spikes_file": "out.h5"
    },
    "reports": {
        "soma": {
            "cells": "hex_O1ExcitatoryPlastic",
            "type": "compartment",
            "variable_name": "v",
            "unit": "mV",
            "dt": 1.0,
            "start_time": 0.0,
            "end_time": 122500.0
        },
        "rho": {
            "cells": "hex_O1ExcitatoryPlastic",
            "type": "synapse",
            "sections": "all",
            "variable_name": "GluSynapse.rho_GB",
            "unit": "nd",
            "dt": 0.025,
            "start_time": 0.0,
            "end_time": 122500.0
        },
	"v_syn": {
            "cells": "hex_O1ExcitatoryPlastic",
            "type": "synapse",
            "sections": "all",
            "variable_name": "GluSynapse.vsyn",
            "unit": "mV",
            "dt": 0.025,
            "start_time": 0.0,
            "end_time": 122500.0
        },
        "effcai_GB": {
            "cells": "hex_O1ExcitatoryPlastic",
            "type": "synapse",
            "sections": "all",
            "variable_name": "GluSynapse.effcai_GB",
            "unit": "mM",
            "dt": 0.025,
            "start_time": 0.0,
            "end_time": 122500.0
        }
    },
    "inputs": {
        "s1_syns_spikes": {
            "node_set": "hex_O1ExcitatoryPlastic",
            "input_type": "spikes",
            "delay": 0.0,
            "duration": 2000.0,
            "module": "synapse_replay",
            "spike_file": "./input_spikes/s1_syn_spikes.dat",
            "source": "hex_O1"
            }
    },
    "connection_overrides": [
        {
            "name": "plasticity",
            "source": "hex_O1Excitatory",
            "target": "hex_O1Excitatory",
            "modoverride": "GluSynapse",
            "weight": 1.0
        },
        {
            "name": "init",
            "source": "hex_O1",
            "target": "hex_O1",
	    "spont_minis": 0.0,
            "weight": 1.0
        },
        {
            "name": "VPM_init",
            "source": "proj_Thalamocortical_VPM_Source",
            "target": "hex_O1",
            "spont_minis": 0.0,
            "weight": 1.0
        },
        {
            "name": "POm_init",
            "source": "proj_Thalamocortical_POM_Source",
            "target": "hex_O1",
            "spont_minis": 0.0,
            "weight": 1.0
        }
    ]
}

s1_syn_spikes.dat

/scatter
0 0
817 170983
829 170983
884 170983
995 170983
996 170983
1000 170983
1080 170983
1135 170983
1403 170983

@anilbey
Copy link
Contributor

anilbey commented Nov 9, 2023

Thanks @dhuruvapriyan for the detailed report.

It is good that you narrowed the issue down to the use of USE parameter. I will try debug from here.

@anilbey
Copy link
Contributor

anilbey commented Nov 30, 2023

Hi @dhuruvapriyan sorry for the late reply I returned to the office this week. In the meantime my colleague was looking into this issue.

I believe the issue is raised from these lines.
https://github.com/BlueBrain/BlueCelluLab/blob/main/bluecellulab/synapse/synapse_types.py#L259

        self.hsynapse.Use_d = self.syn_description["Use_d_TM"] * \
            self.syn_description["u_scale_factor"]
        self.hsynapse.Use_p = self.syn_description["Use_p_TM"] * \
            self.syn_description["u_scale_factor"]

In the code above Use_d_TM is multiplied by u_scale_factor which is not done in GluSynapseHelper.hoc shared via Zenodo.
https://github.com/BlueBrain/neurodamus/blob/21a637eba7ab605eae9967a19b3f35e0d47e075e/tests/share/hoc/GluSynapseHelper.hoc#L65

@anilbey
Copy link
Contributor

anilbey commented Nov 30, 2023

Ignore my previous message. That part of the code looks correct.

The issue may be due to the use of these parameters in the simulation_config.

		"tau_effca_GB": 278.3177658387,
		"gamma_d_GB" : 101.5387594661,
		"gamma_p_GB" : 216.1841700668

What are these parameters? Is Neurodamus reading them?
Neither bluecellulab nor GluSynapseHelper.hoc (https://github.com/BlueBrain/neurodamus/blob/21a637eba7ab605eae9967a19b3f35e0d47e075e/tests/share/hoc/GluSynapseHelper.hoc#L65) is accessing those parameters.

@dhuruvapriyan
Copy link
Author

Hi @anilbey,

I'm sorry for not getting back to you sooner.

Those are the parameters fitted for plasticity experiments. All three parameters are computed along with others to determine depression and potentiation thresholds.

https://github.com/BlueBrain/EModelRunner/blob/95afd0e6f293b6a419ab474d40ba14ea43cad71a/examples/synplas_sample_dir/config/fit_params.json#L4

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants