Auryn simulator

Simulator for spiking neural networks with synaptic plasticity

User Tools

Site Tools

Trace:
examples:sim_isp_orig
no way to compare when less than two revisions

Differences

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


—examples:sim_isp_orig [2014/01/20 10:19] (current) – created zenke
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 10:19 by zenke