Auryn simulator

Simulator for spiking neural networks with synaptic plasticity

User Tools

Site Tools


examples:sim_isp_orig

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

examples:sim_isp_orig [2014/01/20 11:19] (current)
zenke created
Line 1: Line 1:
 +====== sim_isp_orig ======
 +
 +This simulation is the adapted from the network used to generate Figure 4 in [1]. It simulates a balanced network of 10,000 sparsely connected integrate-and-fire neurons[2] and features ​
 +plastic inhibitory-exciatory synapses. ​
 +
 +{{:​examples:​isp_fig4.png?​300|}}
 +
 +[1] Vogels, T.P., Abbott, L.F., 2005. Signal propagation and logic gating in networks of integrate-and-fire neurons. J Neurosci 25, 10786.
 +[[http://​www.ncbi.nlm.nih.gov/​pubmed/​16291952|PubMed]]
 +
 +[2] Vogels, T.P., Sprekeler, H., Zenke, F., Clopath, C., Gerstner, W., 2011. Inhibitory Plasticity Balances Excitation and Inhibition in Sensory Pathways and Memory Networks. Science 334, 1569 –1573.
 +[[http://​www.ncbi.nlm.nih.gov/​pubmed/​22075724|PubMed]]
 +
 +===== Running the program =====
 +
 +To run the program after [[manual:​start#​first steps|compilation]] it suffices to call it as follows
 +<code shell>
 +./​sim_isp_orig
 +</​code>​
 +this will trigger a 1000s simulation which creates the following output files
 +<code shell>
 +   4 -rw-r--r-- 1 zenke lcn1     954 Jan 20 10:52 out.0.log
 +5456 -rw-r--r-- 1 zenke lcn1 5555044 Jan 20 10:52 out.0_e.ras
 +1412 -rw-r--r-- 1 zenke lcn1 1433909 Jan 20 10:52 out.0_i.ras
 + 936 -rw-r--r-- 1 zenke lcn1  951025 Jan 20 10:52 out.0.volt
 + 888 -rw-r--r-- 1 zenke lcn1  901395 Jan 20 10:52 out.0.ampa
 + 928 -rw-r--r-- 1 zenke lcn1  942420 Jan 20 10:52 out.0.gaba
 +</​code>​
 +
 +To see inhibitory plasticity in action the network can be initialized with small inithal I-E inhibitory weights. Try running the network with the following parameter
 +<code shell>
 +./​sim_isp_orig --winh 0.01
 +</​code>​
 +if you then plot the spiking activity recorded in ''​out.0_e.ras''​ after a couple of seconds it will look highly synchronous and pathologic as in the following ​
 +
 +{{:​examples:​isp_orig_out1.png?​300|}}
 +
 +if you wait for roughly a minute this type of activity will give way to asynchronous irregular activity of the following type
 +
 +{{:​examples:​isp_orig_out2.png?​300|}}
 +
 +==== Running in parallel ====
 +
 +This version of the simulation can be run in parallel like any other Auryn simulation. See  [[manual:​parallel_execution|Running simulations in parallel]] for details.
 +
 +
 +===== The important bits =====
 +
 +The following snipped of code first initializes two [[manual:​NeuronGroup|NeuronGroups]] of type [[manual:​TIFGroup]] and then initializes the neuronal membrane potentials to random values that follow a Gaussian distribution centered at -60mV with a standard deviation of 5mV. It furthermore creates an array of Poisson neurons called [[manual:​PoissonGroup]],​ before moving on to set up the inhibitory and excitatory connections. Note that the connections from the inhibitory to the excitatory population is of type [[manual:​SymmetricSTDPConnection]] which implements inhibitory synaptic plasticity as specified in Vogels et al. 2011. 
 +
 +<code c++>
 + ​logger->​msg("​Setting up neuron groups ...",​PROGRESS,​true);​
 +    TIFGroup * neurons_e = new TIFGroup(NE);​
 +    TIFGroup * neurons_i = new TIFGroup(NI);​
 +    neurons_e->​random_mem(-60e-3,​5e-3);​
 +    neurons_i->​random_mem(-60e-3,​5e-3);​
 +    PoissonGroup * poisson = new PoissonGroup(NP,​poisson_rate);​
 +
 +
 +    logger->​msg("​Setting up connections ...",​PROGRESS,​true);​
 +    SparseConnection * con_ei = new SparseConnection(neurons_e,​neurons_i,​wei*w,​sparseness,​GLUT);​
 +    SparseConnection * con_ii = new SparseConnection(neurons_i,​neurons_i,​gamma*w,​sparseness,​GABA);​
 +
 +    SparseConnection * con_exte = new SparseConnection(poisson,​neurons_e,​0,​sparseness_afferents,​GLUT);​
 +
 +    SparseConnection * con_ee = new SparseConnection(neurons_e,​neurons_e,​w,​sparseness,​GLUT);​
 +    SymmetricSTDPConnection * con_ie = new SymmetricSTDPConnection(neurons_i,​neurons_e,​
 +            gamma*w,​sparseness,​
 +            gamma*eta,​kappa,​tau_stdp,​wmax,​
 +            GABA);
 +</​code>​
 +
 +Note also that in the above code the [[manual:​PoissonGroup]] is connected with zero weights only. Certain of the these random sparse zero weights can be changed to non-zero values later in the simulation thus providing external input to the network. By default external input is given by external input currents and the [[manual:​PoissonGroup]] is simply dragged along (it would indeed be more efficient to not define it at all in this case). The actual setting up of a stimulus happens here
 +<code c++>
 +// stimulus
 +    if (!stimfile.empty()) {
 +        char ch;
 +        NeuronID counter = 0;
 +        ifstream fin(stimfile.c_str());​
 +        while (!fin.eof() && counter<​NE) {
 +            ch = fin.get();
 +            if (ch == '​1'​) {
 +                if (poisson_stim==true) {
 +                    for (int i = 0 ; i < NP ; ++i)
 +                        con_exte->​set(i,​counter,​w_ext);​
 +                } else {
 +                    neurons_e->​set_bg_current(counter,​chi*bg_current);​
 +                }
 +            }
 +            counter++;
 +        }
 +        fin.close();​
 +    }
 +</​code>​
 +
 +To run the network and keep activity alive even in the case of highly synchronized activity all cells in the network receive depolarizing external currents. This is set up here
 +<code c++>
 +for ( int j = 0; j < NE ; j++ ) {
 +  neurons_e->​set_bg_current(j,​bg_current);​
 +}
 + 
 +for ( int j = 0; j < NI ; j++ ) {
 +  neurons_i->​set_bg_current(j,​bg_current);​
 +}
 +</​code>​
 +note however that ''​set_bg_current''​ is not a generic function of [[manual:​NeuronGroup]],​ but a pecularity of the neuron model implemented as [[manual:​TIFGroup]]. To give currents to arbitrary [[manual:​NeuronGroup|NeuronGroups]] a [[manual:​CurrentInjector]] has to be used.
 +
 +
 +===== The full program =====
 +<code c++>
 +/* 
 +* Copyright 2014 Friedemann Zenke
 +*
 +* This file is part of Auryn, a simulation package for plastic
 +* spiking neural networks.
 +
 +* Auryn is free software: you can redistribute it and/or modify
 +* it under the terms of the GNU General Public License as published by
 +* the Free Software Foundation, either version 3 of the License, or
 +* (at your option) any later version.
 +
 +* Auryn is distributed in the hope that it will be useful,
 +* but WITHOUT ANY WARRANTY; without even the implied warranty of
 +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. ​ See the
 +* GNU General Public License for more details.
 +
 +* You should have received a copy of the GNU General Public License
 +* along with Auryn. ​ If not, see <​http://​www.gnu.org/​licenses/>​.
 +*/
 +#include <​iomanip>​
 +#include <​stdlib.h>​
 +#include <​string>​
 +#include <​iterator>​
 +
 +#include <​boost/​program_options.hpp>​
 +
 +#include "​auryn_global.h"​
 +#include "​auryn_definitions.h"​
 +#include "​System.h"​
 +#include "​Logger.h"​
 +
 +#include "​NeuronGroup.h"​
 +#include "​TIFGroup.h"​
 +#include "​PoissonGroup.h"​
 +#include "​SparseConnection.h"​
 +#include "​SymmetricSTDPConnection.h"​
 +#include "​StateMonitor.h"​
 +#include "​SpikeMonitor.h"​
 +#include "​RateChecker.h"​
 +
 +#define NE 8000
 +#define NI 2000
 +#define NP 1000
 +#define NSTIM 20
 +
 +using namespace std;
 +namespace po = boost::​program_options;​
 +
 +int main(int ac, char* av[]) 
 +{
 + double w = 0.3 ;
 + double w_ext = w ;
 + double gamma = 10. ;
 + double wmax = 10*gamma*w;
 +
 + double sparseness = 0.02 ;
 +
 + double eta = 1e-4 ;
 + double kappa = 3. ;
 + double tau_stdp = 20e-3 ;
 + bool stdp_active = true;
 + bool poisson_stim = false;
 + double winh = -1;
 + double wei = 1;
 + double chi = 10.;
 + double bg_current = 2e-2;
 +
 + double poisson_rate = 100.;
 + double sparseness_afferents = 0.05;
 +
 + bool quiet = false;
 + double simtime = 1000. ;
 + NeuronID record_neuron = 30;
 + // handle command line options
 +
 + string infilename = "";​
 + string outputfile = "​out";​
 + string stimfile = "";​
 + string strbuf ;
 +
 + int errcode = 0;
 +
 +    try {
 +
 +        po::​options_description desc("​Allowed options"​);​
 +        desc.add_options()
 +            ("​help",​ "​produce help message"​)
 +            ("​quiet",​ "quiet mode")
 +            ("​load",​ po::​value<​string>​(),​ "input weight matrix"​)
 +            ("​out",​ po::​value<​string>​(),​ "​output filename"​)
 +            ("​stimfile",​ po::​value<​string>​(),​ "​stimulus file")
 +            ("​eta",​ po::​value<​double>​(),​ "​learning rate")
 +            ("​kappa",​ po::​value<​double>​(),​ "​target rate")
 +            ("​simtime",​ po::​value<​double>​(),​ "​simulation time")
 +            ("​active",​ po::​value<​bool>​(),​ "​toggle learning"​)
 +            ("​poisson",​ po::​value<​bool>​(),​ "​toggle poisson stimulus"​)
 +            ("​winh",​ po::​value<​double>​(),​ "​inhibitory weight multiplier"​)
 +            ("​wei",​ po::​value<​double>​(),​ "ei weight multiplier"​)
 +            ("​chi",​ po::​value<​double>​(),​ "chi current multiplier"​)
 +        ;
 +
 +        po::​variables_map vm;        ​
 +        po::​store(po::​parse_command_line(ac,​ av, desc), vm);
 +        po::​notify(vm); ​   ​
 +
 +        if (vm.count("​help"​)) {
 +            cout << desc << "​\n";​
 +            return 1;
 +        }
 +
 +        if (vm.count("​quiet"​)) {
 + quiet = true;
 +        } 
 +
 +        if (vm.count("​load"​)) {
 +            cout << "input weight matrix " ​
 +                 <<​ vm["​load"​].as<​string>​() << "​.\n";​
 + infilename = vm["​load"​].as<​string>​();​
 +        } 
 +
 +        if (vm.count("​out"​)) {
 +            cout << "​output filename " ​
 +                 <<​ vm["​out"​].as<​string>​() << "​.\n";​
 + outputfile = vm["​out"​].as<​string>​();​
 +        } 
 +
 +        if (vm.count("​stimfile"​)) {
 +            cout << "​stimfile filename " ​
 +                 <<​ vm["​stimfile"​].as<​string>​() << "​.\n";​
 + stimfile = vm["​stimfile"​].as<​string>​();​
 +        } 
 +
 +        if (vm.count("​eta"​)) {
 +            cout << "eta set to " ​
 +                 <<​ vm["​eta"​].as<​double>​() << "​.\n";​
 + eta = vm["​eta"​].as<​double>​();​
 +        } 
 +
 +        if (vm.count("​kappa"​)) {
 +            cout << "kappa set to " ​
 +                 <<​ vm["​kappa"​].as<​double>​() << "​.\n";​
 + kappa = vm["​kappa"​].as<​double>​();​
 +        } 
 +
 +        if (vm.count("​simtime"​)) {
 +            cout << "​simtime set to " ​
 +                 <<​ vm["​simtime"​].as<​double>​() << "​.\n";​
 + simtime = vm["​simtime"​].as<​double>​();​
 +        } 
 +
 +        if (vm.count("​active"​)) {
 +            cout << "stdp active : " ​
 +                 <<​ vm["​active"​].as<​bool>​() << "​.\n";​
 + stdp_active = vm["​active"​].as<​bool>​();​
 +        } 
 +
 +        if (vm.count("​poisson"​)) {
 +            cout << "​poisson active : " ​
 +                 <<​ vm["​poisson"​].as<​bool>​() << "​.\n";​
 + poisson_stim = vm["​poisson"​].as<​bool>​();​
 +        } 
 +
 +
 +        if (vm.count("​winh"​)) {
 +            cout << "inhib weight multiplier : " ​
 +                 <<​ vm["​winh"​].as<​double>​() << "​.\n";​
 + winh = vm["​winh"​].as<​double>​();​
 +        } 
 +
 +        if (vm.count("​wei"​)) {
 +            cout << "ei weight multiplier : " ​
 +                 <<​ vm["​wei"​].as<​double>​() << "​.\n";​
 + wei = vm["​wei"​].as<​double>​();​
 +        } 
 +
 +        if (vm.count("​chi"​)) {
 +            cout << "chi multiplier : " ​
 +                 <<​ vm["​chi"​].as<​double>​() << "​.\n";​
 + chi = vm["​chi"​].as<​double>​();​
 +        } 
 +
 +    }
 +    catch(exception&​ e) {
 +        cerr << "​error:​ " << e.what() << "​\n";​
 +        return 1;
 +    }
 +    catch(...) {
 +        cerr << "​Exception of unknown type!\n";​
 +    }
 +
 + // BEGIN Global definitions
 + mpi::​environment env(ac, av);
 + mpi::​communicator world;
 + communicator = &world;
 +
 + stringstream oss;
 + oss << outputfile << "​."​ << world.rank();​
 + outputfile = oss.str();
 + oss << "​.log";​
 + string logfile = oss.str();
 + logger = new Logger(logfile,​world.rank());​
 +
 + sys = new System(&​world);​
 + // END Global definitions
 +
 +
 +
 + logger->​msg("​Setting up neuron groups ...",​PROGRESS,​true);​
 + TIFGroup * neurons_e = new TIFGroup(NE);​
 + TIFGroup * neurons_i = new TIFGroup(NI);​
 + neurons_e->​random_mem(-60e-3,​5e-3);​
 + neurons_i->​random_mem(-60e-3,​5e-3);​
 + PoissonGroup * poisson = new PoissonGroup(NP,​poisson_rate);​
 +
 +
 + logger->​msg("​Setting up connections ...",​PROGRESS,​true);​
 + SparseConnection * con_ei = new SparseConnection(neurons_e,​neurons_i,​wei*w,​sparseness,​GLUT);​
 + SparseConnection * con_ii = new SparseConnection(neurons_i,​neurons_i,​gamma*w,​sparseness,​GABA);​
 +
 + SparseConnection * con_exte = new SparseConnection(poisson,​neurons_e,​0,​sparseness_afferents,​GLUT);​
 +
 + SparseConnection * con_ee = new SparseConnection(neurons_e,​neurons_e,​w,​sparseness,​GLUT);​
 + SymmetricSTDPConnection * con_ie = new SymmetricSTDPConnection(neurons_i,​neurons_e,​
 + gamma*w,​sparseness,​
 + gamma*eta,​kappa,​tau_stdp,​wmax,​
 + GABA);
 +
 + if (!infilename.empty()) {
 + sys->​load_network_state(infilename);​
 + }
 +
 + if (winh>​=0)
 + con_ie->​set_all(winh);​
 +
 + logger->​msg("​Setting up monitors ...",​PROGRESS,​true);​
 + strbuf = outputfile;
 + strbuf += "​_e.ras";​
 + SpikeMonitor * smon_e = new SpikeMonitor( neurons_e , strbuf.c_str() );
 +
 + strbuf = outputfile;
 + strbuf += "​_i.ras";​
 + SpikeMonitor * smon_i = new SpikeMonitor( neurons_i, strbuf.c_str() );
 + smon_i->​set_offset(NE);​
 +
 + strbuf = outputfile;
 + strbuf += "​.volt";​
 + StateMonitor * vmon = new StateMonitor( neurons_e, record_neuron,​ "​mem",​ strbuf.c_str() );
 +
 + strbuf = outputfile;
 + strbuf += "​.ampa";​
 + StateMonitor * amon = new StateMonitor( neurons_e, record_neuron,​ "​g_ampa",​ strbuf.c_str() );
 +
 + strbuf = outputfile;
 + strbuf += "​.gaba";​
 + StateMonitor * gmon = new StateMonitor( neurons_e, record_neuron,​ "​g_gaba",​ strbuf.c_str() );
 +
 + RateChecker * chk = new RateChecker( neurons_e , 0.1 , 1000. , 100e-3);
 +
 + for ( int j = 0; j < NE ; j++ ) {
 +   neurons_e->​set_bg_current(j,​bg_current);​
 + }
 +
 + for ( int j = 0; j < NI ; j++ ) {
 +   neurons_i->​set_bg_current(j,​bg_current);​
 + }
 +
 + // stimulus
 + if (!stimfile.empty()) {
 + char ch;
 + NeuronID counter = 0;
 + ifstream fin(stimfile.c_str());​
 + while (!fin.eof() && counter<​NE) { 
 + ch = fin.get(); ​
 + if (ch == '​1'​) {
 + if (poisson_stim==true) {
 + for (int i = 0 ; i < NP ; ++i)
 + con_exte->​set(i,​counter,​w_ext);​
 + } else {
 + neurons_e->​set_bg_current(counter,​chi*bg_current);​
 + }
 + }
 + counter++;​
 +
 + fin.close();​
 + }
 +
 +
 + logger->​msg("​Simulating ...",​PROGRESS,​true);​
 + con_ie->​stdp_active = stdp_active;​
 +
 + sys->​run(simtime);​
 +
 + logger->​msg("​Saving network state ...",​PROGRESS,​true);​
 + sys->​save_network_state(outputfile);​
 +
 + logger->​msg("​Freeing ...",​PROGRESS,​true);​
 + delete sys;
 +
 + logger->​msg("​Exiting ...",​PROGRESS,​true);​
 + return errcode;
 +}
 +</​code>​
  
examples/sim_isp_orig.txt · Last modified: 2014/01/20 11:19 by zenke