Welcome to dipde’s documentation!

The population density approach in computational neuroscience seeks to understand the statistical evolution of a large population of homogeneous neurons. Beginning with the work of Knight and Sirovich [1] (see also [2]), the approach typically formulates a partial integro-differential equation for the evolution of the voltage probability distribution receiving synaptic activity, and under the influence of neural dynamics. Neuronal dynamics typically follow from the assumption of a leaky integrate-and fire model. We implement a numerical scheme for computing the time evolution of the master equation for populations of leaky integrate-and-fire neurons with shot-noise synapses (for a similar approach, see [3]).

For a Quick Start, you can begin with the Examples (dipde.examples). For a more complicated example, check out the Cortical Column Model, a 16 population example based on the Potjans and Diesmann [4] leaky integrate-and-fire cortical column model.

Details on the classes used to define a dipde simulation can be found in Core functionality (dipde.internals).

Contents:

User Guide

This guide is a resource for installing the dipde package. It is maintained by the Allen Institute for Brain Science.

Quick Start Install Using Pip

1. Pip install for a single user:

   pip install git+https://github.com/AllenInstitute/dipde.git --user
   

2. Uninstalling:

   pip uninstall dipde
   

Required Dependencies

• NumPy
• SciPy
• MatPlotLib
• SymPy

For a Quick Start, you can begin with the Examples (dipde.examples).

Details on the classes used to define a dipde simulation can be found in Core functionality (dipde.internals).

Packages

Examples (dipde.examples)

These quick-start examples demonstrate the configuration, simulation, and plotting of simple dipde networks.

A preamble imports the classes necessary for simulation:

import matplotlib.pyplot as plt
from dipde.internals.internalpopulation import InternalPopulation
from dipde.internals.externalpopulation import ExternalPopulation
from dipde.internals.simulation import Simulation
from dipde.internals.connection import Connection as Connection

Each of the example simulations below use the same general settings:

# Settings:
t0 = 0.
dt = .0001
dv = .0001
tf = .1
update_method = 'approx'
approx_order = 1
tol = 1e-14
verbose = True

Here dt defines the time step of the simulation (in seconds), dv defines the granularity of the voltage domain of the Internal populations (in volts), and tf defines the total duration of the simulation (in seconds). The update method allows the user to change the time-stepping method for the forward evolution of the time domain of each Internal population. The approx_order and tol options fine-tune the time stepping, to provide a tradeoff between simulation time and numerical precision when the time evolution method is ‘approx,’ as opposed to ‘exact.’ Setting verbose to True causes the Simulation to print the current time after the evaluation of each time step.

Singlepop

Download singlepop.py

The singlepop simulation provides a simple feedforward topology that uses every major class in the core library. A single 100 Hz External population population (here specified as a string, although a floating point or integer specification will work also) provides excitatory input. This is connected to an Internal population (modeled as a population density pde) via a delta-distributed synaptic weight distribution, with 5 mV strength. The in-degree (nsyn) of this Connection is set to 1 for this example; in general, this serves as a multiplier of the input firing rate of the source population. The internal population has a linearly binned voltage domain from v_min to v_max. No negative bins (i.e. v_min < 0) are required here, because no negative synaptic inputs (‘weights’ in the Connection object) are defined.

The mean firing rate of the Internal population resulting from this simulation is

plotted below, along with the code used to generate the plot:

Singlepop (recurrent)

Download singlepop_recurrent.py

The next example is identical to the singlepop example, with the exception of an additional recurrent connection from the internal population to itself. This demonstrates that the Connection object is used to connect both Internal and External population sources to an Internal target population. Attempting to create a Connection with an External population as a target will result in an AttributeError. The additional excitatory input resulting from the extra recurrent connection results in a slightly increased steady-state firing rate.

Excitatory/Inhibitory

Download excitatory_inhibitory.py

Inhibitory connections are formed the same as excitatory connections, except the sign of the connection weight distribution is changed. In this example, the background External population is connected to the Internal population with two connections, one excitatory and one inhibitory. Because of this negative input, negative voltage values are possible, requiring that v_min<0.

Singlepop (sine)

Download singlepop_sine.py

By modifying the argument to the ExternalPopulation constructor, we can define a time-varying external input. This input is parsed by Sympy using the lambdify module; any lambdify-able expression with the independent variable ‘t’ is acceptable. Negative values of the firing rate function will throw an Exception.

Singlepop (exponential)

Download singlepop_exponential_distribution.py

Up until now, each connection object was defined with a single synaptic weight. However, connection objects are fundamentally defined by synaptic weight distributions (See Iyer et al. 2013 [1]). Here we consider the effect of an exponentially distributed synaptic weight distribution with mean equal to 5 mV. The steady state firing rate under this stimulus is significantly higher than the delta-distributed example (singlepop, above). In general, any Continuous synaptic weight distribution can be defined using scipy.stats module; however it will be discretized into a finite set of input weights. Discrete distributions can be specified directly, and will be used exactly.

In general, the special case of exponentially distributed synaptic weights admits an analytical steady-state firing rate solution (see Richardson and Swarbrick, 2010 [2], and Iyer et al. 2014 [3]). The figure above includes the mean firing rate computed for the internal population, as well as a ‘*’ indicating a semi-analytic prediction (~8.669 Hz) of the steady-state firing rate for the model. This steady state prediction is computed with the code below, plotted above with an ‘*’.

import scipy.integrate as spint
import numpy as np

R_e = 100
a_e = .005
tau = .02
v_th = .02
f = lambda c: (1./c)*(1-a_e*c)**(tau*R_e)*(np.exp(v_th*c)/(1-a_e*c)-1.)
y, _ = spint.quad(f,0.0,1./a_e)
steady_state_prediction = 1./(tau*y)    # >>> 8.6687760498

[1]	Iyer, R., Menon, V., Buice, M., Koch, C., & Mihalas, S. (2013). The Influence of Synaptic Weight Distribution on Neuronal Population Dynamics. PLoS Computational Biology, 9(10), e1003248. doi:10.1371/journal.pcbi.1003248

[2]	Richardson, Magnus J. E. & Swarbrick, Rupert (2010). Firing-Rate Response of a Neuron Receiving Excitatory and Inhibitory Synaptic Shot Noise. Phys. Rev. Lett., 105, 178102.

[3]	Iyer, R., Cain, N., & Mihalas, S. (2014). Dynamics of excitatory-inhibitory neuronal networks with exponential synaptic weight distributions. Cosyne Abstracts 2014, Salt Lake City USA.

Core functionality (dipde.internals)

External Population

class dipde.internals.externalpopulation.ExternalPopulation(firing_rate, record=False, **kwargs)

Bases: object

External (i.e. background) source for connections to Internal Populations.

This class provides a background drive to internal population. It is used as the source argument to a connection, in order to provide background drive.

firing_rate : numeric, str

Output firing rate of the population. If numeric type, this defines a constant background generator; if str, it is interpreted as a SymPy function with independent variable ‘t’.

record : bool (default=False)

If True, a history of the output firing rate is recorded (firing_rate_record attribute).

**kwargs

Any additional keyword args are stored as metadata (metadata attribute).

self.firing_rate_string : str

String representation of firing_rate input parameter.

self.metadata : dict

Dictonary of metadata, constructed by kwargs not parsed during construction.

self.firing_rate_record : list

List of firing rates recorded during Simulation.

self.t_record : list

List of times that firing rates were recorded during Simulation.

curr_firing_rate

Property that accesses the current firing rate of the population (Hz).

firing_rate(t)

Firing rate of the population at time t (Hz).

initialize()

Initialize the population at the beginning of a simulation.

Calling this method resets the recorder that tracks firing rate during a simulation. This method is called by the Simulation object ( initialization method), but can also be called by a user when defining an alternative time stepping loop.

initialize_firing_rate_recorder()

Initialize recorder at the beginning of a simulation.

This method is typically called by the initialize method rather than on its own. It resets the lists that track the firing rate during execution of the simulation.

update()

Update the population one time step.

This method is called by the Simulation object to update the population one time step. In turn, this method calls the update_firing_rate_recorder method to register the current firing rate with the recorder.

update_firing_rate_recorder()

Record current time and firing rate, if record==True.

This method is called once per time step. If record is True, calling this method will append the current time and firing rate to the firing rate recorder.

Internal Population

class dipde.internals.internalpopulation.InternalPopulation(tau_m=0.02, v_min=-0.1, v_max=0.02, dv=0.0001, record=True, curr_firing_rate=0.0, update_method='approx', approx_order=None, tol=1e-12, norm=inf, **kwargs)

Bases: object

Population density class

This class encapulates all the details necessary to propagate a population density equation driven by a combination of recurrent and background connections. The voltage (spatial) domain discretization is defined by linear binning from v_min to v_max, in steps of dv (All in units of volts). The probability densities on this grid are recorded pv, and must always sum to 1.

tau_m : float (default=.02)

Time constant (unit: 1/sec) of neuronal population.

v_min : float (default=-.1)

Minimum of voltage domain (unit: volt).

v_max : float (default=.02)

Maximum of voltage domain (Absorbing boundary), i.e spiking threshold (unit: volt).

dv : float (default=.0001)

Voltage domain discritization size (unit: volt).

record : bool (default=False)

If True, a history of the output firing rate is recorded (firing_rate_record attribute).

curr_firing_rate : float (default=0.0

Initial/Current firing rate of the population (unit: Hz).

update_method : str ‘approx’ or ‘exact’ (default=’approx’)

Method to update pv (exact can be quite slow).

approx_order : int or None (default=None)

Maximum Taylor series expansion order to use when computing update to pv.

tol : float (default=1e-12)

Error tolerance used when computing update to pv.

norm : non-zero int, np.inf, -np.inf, or ‘fro’ (default=np.inf)

Vector norm used in computation of tol.

**kwargs

Any additional keyword args are stored as metadata (metadata attribute).

self.edges : np.array

Vector defining the boundaries of voltage bins.

self.pv : np.array

Vector defining the probability mass in each voltage bin (self.pv.sum() = 1).

self.firing_rate_record : list

List of firing rates recorded during Simulation.

self.t_record : list

List of times that firing rates were recorded during Simulation.

self.leak_flux_matrix : np.array

Matrix that defines the flux between voltage bins.

get_total_flux_matrix()

Create a total flux matrix by summing presynaptic inputs and the leak matrix.

initialize()

Initialize the population at the beginning of a simulation.

In turn, this method:

1. Initializes the voltage edges (self.edges) and probability mass in each bin (self.pv),
2. Creates an initial dictionary of inputs into the population, and
3. Resets the recorder that tracks firing rate during a simulation.

This method is called by the Simulation object (initialization method), but can also be called by a user when defining an alternative time stepping loop.

initialize_edges()

Initialize self.edges and self.leak_flux_matrix attributes.

This method initializes the self.edges attribute based on the v_min, v_max, and dv settings, and creates a corresponding leak flux matrix based on this voltage discretization.

initialize_firing_rate_recorder()

Initialize recorder at the beginning of a simulation.

This method is typically called by the initialize method rather than on its own. It resets the lists that track the firing rate during execution of the simulation.

initialize_probability()

Initialize self.pv to delta-distribution at v=0.

initialize_total_input_dict()

Initialize dictionary of presynaptic inputs at beginning of simulation

This method is typically called by the initialize method rather than on its own. It creates a dictionary of synaptic inputs to the population.

n_bins

Number of probability mass bins.

n_edges

Number of probability mass bin edges.

plot_probability_distribution(ax=None)

Convenience method to plot voltage distribution.

ax : None or matplotlib.pyplot.axes object (default=None)

Axes object to plot distribution on. If None specified, a figure and axes object will be created.

source_connection_list

List of all connections that are a source for this population.

update()

Update the population one time step.

This method is called by the Simulation object to update the population one time step. In turn, this method:

1. Calls the update_total_input_dict method to gather the current strengths of presynaptic input populations,
2. Calls the update_propability_mass method to propagate self.pv one time-step,
3. Calls the update_firing_rate method to compute the firing rate of the population based on flux over threshold, and
4. Calls the update_firing_rate_recorder method to register the current firing rate with the recorder.

update_firing_rate()

Update curr_firing_rate attribute based on the total flux of probability mass across threshold.

update_firing_rate_recorder()

Record current time and firing rate, if record==True.

This method is called once per time step. If record is True, calling this method will append the current time and firing rate to the firing rate recorder.

update_propability_mass()

Create a total flux matrix, and propogate self.pv one time-step.

update_total_input_dict()

Update the input dictionary based on the current firing rates of presynaptic populations.

Connection

Module containing Connection class, connections between source and target populations.

class dipde.internals.connection.Connection(source, target, nsyn, **kwargs)

Bases: object

Class that connects dipde source population to dipde target population.

The Connection class handles all of the details of propogating connection information between a source population (dipde ExtermalPopulation or InternalPopulation) and a target population (dipde InternalPopulation). Handles delay information via a delay queue that is rolled on each timestep, and reuses connection information by using a ConnectionDistribution object with a specific signature to avoid duplication among identical connections.

source : InternalPopulation or ExternalPopulation

Source population for connection.

target : InternalPopulation

Target population for connection.

nsyn : int

In-degree of connectivity from source to target.

weights : list

Weights defining synaptic distribution (np.ndarray).

probs : list (same length as weights, and sums to 1)

Probabilities corresponding to weights.

delay: float (default=0)

Transmission delay (units: sec).

metadata: Connection metadata, all other kwargs

curr_delayed_firing_rate

Current firing rate of the source (float).

Property that accesses the firing rate at the top of the delay queue, from the source population.

initialize()

Initialize the connection at the beginning of a simulation.

Calling this method:

1. Initializes a delay queue used to store values of inputs in a last-in-first-out rolling queue.
2. Creates a connection_distribution object for the connection, if a suitable object is not already registered with the simulation-level connection distribution collection.

This method is called by the Simulation object (initialization method), but can also be called by a user when defining an alternative time stepping loop.

initialize_connection_distribution()

Create connection distribution, if necessary.

If the signature of this connection is already registered to the simulation-level connection distribution collection, it is associated with self. If not, it adds the connection distribution to the collection, and associates it with self.

initialize_delay_queue()

Initialiaze a delay queue for the connection.

The delay attribute of the connection defines the transmission delay of the signal from the souce to the target. Firing rate values from the source population are held in a queue, discretized by the timestep, that is rolled over once per timestep. if source is an ExternalPopulation, the queue is initialized to the firing rate at t=0; if the source is an InternalPopulation, the queue is initialized to zero.

update()

Update Connection, called once per timestep.

Connection Distribution

Module containing ConnectionDistribution class, organizing connection details.

class dipde.internals.connectiondistribution.ConnectionDistribution(edges, weights, probs, sparse=True)

Bases: object

Container for matrix that computes element flux, and associated methods.

A ConnectionDistribution object contains a matrix used to compute the time propogation of the voltage probability distribution, and a vector used to compute the flux over threshold at the end of the timestep. It is unique up to the descretization of the target voltage distribution, and the specific synaptic weight distribution of the connection. Therefore, whenever possible, a ConnectionDistribution object is reused for multiple connections with identical signatures, to reduce memory consumption.

edges: np.ndarray

Voltage bin discretization of target population.

weights: np.ndarray

Weights defining the discrete synaptic distribution.

probs : np.ndarray

Probabilities corresponding to weights.

self.threshold_flux_vector : np.ndarray

Vector used to compute over-threshold flux.

self.flux_matrix : np.ndarray

Matrix used to propagate voltage distribution.

initialize()

Initialize connection distribution.

Initialization creates the flux_matrix and threshold_flux_vector. Implemented lazily via a try catch when flux matrix is requested, that does not exist.

signature

A unique key used to organize ConnectionDistributions.

Connection Distribution Collection

class dipde.internals.connectiondistributioncollection.ConnectionDistributionCollection

Bases: dict

Container that organizes connection components for a simulation, to reduce redundancy.

In a simulation, connections that share the same weights and probabilities, as well as the same target bin edges, can make use of the same flux_matrix and threshold_flux_vector. This can significantly improve the overall memory efficiency of the simulation. To facilitate this, each simulation creates a ConnectionDistributionCollection object that indexes the ConnectionDistribution objects according to their signature, and re-uses the for multiple connections.

add_connection_distribution(cd)

Try and add a ConnectionDistribution object, if signature not already used.

Simulation

class dipde.internals.simulation.Simulation(population_list, connection_list, verbose=True)

Bases: object

Initialize and run a dipde simulation.

The Simulation class handles the initialization of population and connection objects, and provides a convenience time stepping loop to drive a network simulation. Typical usage involves the create of populations and connections, construction of a simulation object, and then call to simulation.run()

population_list : list of ExternalPopulation, InternalPopulation objects

List of populations to include in simulation.

connection_list : list of Connection objects

List of connections to include in simulation.

verbose : bool (default=True)

Setting True prints current time-step at each update evaluation.

initialize(t0=0.0)

Initialize simulation, populations, and connections.

This function is typically called by the self.run() method, however can be called independently if defining a new time stepping loop.

t0 : float (default=0.)

Simulation start time (unit=seconds).

run(t0=0.0, dt=0.001, tf=0.1)

Main iteration control loop for simulation

The time step selection must be approximately of the same order as dv for the internal populations, if the ‘approx’ time stepping method is selected.

t0 : float (default=0.)

Simulation start time (unit=seconds), passed to initialize call.

tf : float (default=.1)

Simulation end time (unit=seconds).

dt : float (default=0.001)

Time step (unit=seconds).

Utilities

dipde.internals.utilities.approx_update_method_order(J, pv, dt=0.0001, approx_order=2)

Approximate the effect of a matrix exponential, truncating Taylor series at order ‘approx_order’.

dipde.internals.utilities.approx_update_method_tol(J, pv, tol=2.2e-16, dt=0.0001, norm='inf')

Approximate the effect of a matrix exponential, with residual smaller than tol.

dipde.internals.utilities.assert_probability_mass_conserved(pv)

Assert that probability mass in control nodes sums to 1.

dipde.internals.utilities.descretize(distribution, N, loc=0, scale=1, shape=[])

Compute a discrete approzimation to a scipy.stats continuous distribution.

dipde.internals.utilities.exact_update_method(J, pv, dt=0.0001)

Given a flux matrix, pdate probabilty vector by computing the matrix exponential.

dipde.internals.utilities.flux_matrix(v, w, lam, p=1)

Compute a flux matrix for voltage bins v, weight w, firing rate lam, and probability p.

dipde.internals.utilities.fraction_overlap(a1, a2, b1, b2)

Calculate the fractional overlap between range (a1,a2) and (b1,b2).

Used to compute a reallocation of probability mass from one set of bins to another, assuming linear interpolation.

dipde.internals.utilities.get_v_edges(v_min, v_max, dv)

Construct voltage edegs, linear discretization.

dipde.internals.utilities.get_zero_bin_list(v)

Return a list of indices that map flux back to the zero bin.

dipde.internals.utilities.leak_matrix(v, tau)

Given a list of edges, construct a leak matrix with time constant tau.

dipde.internals.utilities.redistribute_probability_mass(A, B)

Takes two ‘edge’ vectors and returns a 2D matrix mapping each ‘bin’ in B to overlapping bins in A. Assumes that A and B contain monotonically increasing edge values.

Cortical Column Model

Download potjans_diesmann_cortical_column.py

Potjans and Diesmann [1] parameterized a four-layer, two cell-type (i.e excitatory and inhibitory) model of a cortical column. Beginning with their detailed model description, we implemented a version of their model using dipde. The population statistic approach lends itself to quickly analyzing the mean response properties of population-scale firing-rate dynamics of neural models.

Model Output at steady-state:

[1]	Potjans T.C., & Diesmann, M. (2014) The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model. Cerebral Cortex 24: 785–806.

dipde Team

DiPDE is developed by the Modeling, Analysis and Theory group at the Allen Institute, with contributions from Nicholas Cain, Ram Iyer, Vilas Menon, Michael Buice, Tim Fliss, Keith Godfrey, David Feng, and Stefan Mihalas.

Indices and tables

• Index
• Module Index
• Search Page

Citations

[1]	Knight, N.W., Manin, D., & Sirovich, L. (1996) Dynamical models of interacting neuron populations in visual cortex. Symposium on Robotics and Cybernetics; Computational Engineering in Systems Application: 1–5.

[2]	Omurtag, A., Knight, B.W., & Sirovich, L. (2000) On the Simulation of Large Populations of Neurons. Journal of Computational Neuroscience 8: 51–63.

[3]	de Kamps M. (2003) A simple and stable numerical solution for the population density equation. Neural Computation 15: 2129–2146.

[4]	Potjans T.C., & Diesmann, M. (2014) The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model. Cerebral Cortex 24: 785–806.

dipde Manuscripts and Posters

[5]	Cain, N., Iyer, R., Koch, C., & Mihalas, S. (2015) The computational properties of a simplified cortical column model. In Preparation.

[6]	Cain, N., Fliss, T., Menon, V., Iyer, R., Koch, C., & Mihalas, S. (2014) Simulations of the statistics of firing activity of neuronal populations in the entire mouse brain. Program No. 160.02/GG10. 2013 Neuroscience Meeting Planner. Washington, DC: Society for Neuroscience, 2014. Online.

[7]	Iyer, R., Menon, V., Buice, M., Koch, C., & Mihalas, S. (2013). The Influence of Synaptic Weight Distribution on Neuronal Population Dynamics. PLoS Computational Biology, 9(10), e1003248. doi:10.1371/journal.pcbi.1003248

[8]	Iyer, R., Cain, N., & Mihalas, S. (2014). Dynamics of excitatory-inhibitory neuronal networks with exponential synaptic weight . Cosyne Abstracts 2014, Salt Lake City USA.

Citing dipde

Website: © 2015 Allen Institute for Brain Science. DiPDE Simulator [Internet]. Available from: https://github.com/AllenInstitute/dipde.