Creating Simulations in Python

The Pype9 package is organised into sub-packages loosely corresponding to each pipeline (e.g. simulate, plot). The simulate package contains the sub-packages, neuron and nest, which provide the simulator-specific calls to their respective backends.

All classes required to design and run simulations in these packages derive from corresponding classes in the common package, which defines a consistent Public API across all backends. Therefore, code designed to run on with one backend can be switched to another by simply changing the package the simulator-specific classes are imported from (like in PyNN).

Note

The neuron and nest packages can be imported separately. Therefore, only the simulator you plan to use needs to be available on your system.

Simulation Control

Simulation parameters such as time step, delay limits and seeds for pseudo random number generators are set within an instance of the Simulation class. Simulator objects (e.g. cells and connections) can only be instantiated within the context of an active Simulation instance, and there can only be one active Simulation instance at any time.

A Simulation is activated with the with keyword

with Simulation(dt=0.1 * un.ms, seed=12345) as sim:
    # Design simulation here

The simulation is advanced using the run method of the Simulation instance

with Simulation(dt=0.1 * un.ms, seed=12345) as sim:
     # Create simulator objects here
     sim.run(100.0 * un.ms)

this can be done in stages if states or parameters need to be updated mid-simulation

with Simulation(dt=0.1 * un.ms, seed=12345) as sim:
     # Create simulator objects here
     sim.run(50.0 * un.ms)
     # Update simulator object parameters/state-variables
     sim.run(50.0 * un.ms)

After the simulation context exits all objects in the simulator backend are destroyed (unless an exception is thrown) and only recordings can be reliably accessed from the “dead” Pype9 objects.

There are pseudo-random number generator (RNG) seeds that can be passed to Simulation instances, seed which is used to seed the RNGs that generate random dynamic processes during the simulation (e.g. poisson processes), and properties_seed which is used to seed the RNG used to generate probabilistic connectivity rules and the random distribution of cell properties over populations. If seed is None (the default) then it is generated from the current timestamp, if properties_seed is None then it is derived from seed.

Cell Simulations

NineML Dynamics classes can be translated into simulator cell objects using the CellMetaClass class. A metaclass is class of classes, i.e. one whose instantiation is itself a class, such as the type class. CellMetaClass instantiations derive from the Cell class and can be used to represent different classes of neural models, such as Izhikevich or Hodgkin-Huxley for example. From these Cell classes as many cell instances (with their corresponding simulator objects) can be created as required e.g:

# Create Izhikevich cell class by instantiating the CellMetaClass with a
# ninml.Dynamics Izhikevich model
Izhikevich = CellMetaClass('./izhikevich.xml#Izhikevich')
# Parameters and states of the cell class must be provided when the cells
# are instantiated.
# either as keyword args
izhi1 = Izhikevich(a=1, b=2, c=3, d=4, v=-65 * un.mV, u=14 * un.mV / un.ms)
# or from a nineml.DynamicsProperties object
izhi3 = Izhikevich('./izhikevich.xml#IzhikevichBurster')

If the specified Dynamics class has not been built before the CellMetaClass will automatically generate the required source code for the model, compile it, and load it into the simulator namespace. This can happen either inside or outside of an active Simulation instance. However, the cells objects themselves must be instantiated within a Simulation instance.

# The cell class can be created outside the simulation context
Izhikevich = CellMetaClass('./izhikevich.xml#Izhikevich')
with Simulation(dt=0.1 * un.ms) as sim:
    # The cell object must be instantiated within the simulation context
    izhi = Izhikevich(a=1, b=2, c=3, d=4, v=-65 * un.mV,
                      u=14 * un.mV / un.ms)
    sim.run(1000.0 * un.ms)

The data can be recorded from every send port and state variable in the NineML Dynamics class using the record method of the Cell class. The recorded data can then be accessed with the recording method. The recordings will be Neo format.

Izhikevich = CellMetaClass('./izhikevich.xml#Izhikevich')
with Simulation(dt=0.1 * un.ms) as sim:
    izhi = Izhikevich(a=1, b=2, c=3, d=4, v=-65 * un.mV,
                      u=14 * un.mV / un.ms)
    # Specify the variables to record
    izhi.record('v')
    sim.run(1000.0 * un.ms)
# Retrieve the recording
v = izhi.recording('v')

Transitions between regimes can be recorded using record_regime and retrieved using regime_epochs

Izhikevich = CellMetaClass('./izhikevich.xml#Izhikevich')
with Simulation(dt=0.1 * un.ms) as sim:
    izhi = Izhikevich(a=1, b=2, c=3, d=4, v=-65 * un.mV,
                      u=14 * un.mV / un.ms)
    # Specify the variables to record
    izhi.record_regime()
    sim.run(1000.0 * un.ms)
# Retrieve the regime changes
epochs = izhi.regime_epochs()

Data in Neo format can be “played” into receive ports of the Cell

neo_data = neo.PickleIO('./data/my_recording.neo.pkl').read()
Izhikevich = CellMetaClass('./izhikevich.xml#Izhikevich')
with Simulation(dt=0.1 * un.ms) as sim:
    izhi = Izhikevich(a=1, b=2, c=3, d=4, v=-65 * un.mV,
                      u=14 * un.mV / un.ms)
    # Play analog signal (must be of current dimension) into 'i_syn'
    # analog-receive port.
    izhi.play('i_syn', neo_data.analogsignals[0])
    sim.run(1000.0 * un.ms)

States and parameters can be accessed and set using the attributes of the Cell objects

Izhikevich = CellMetaClass('./izhikevich.xml#Izhikevich')
with Simulation(dt=0.1 * un.ms) as sim:
    izhi = Izhikevich(a=1, b=2, c=3, d=4)
    sim.run(500.0 * un.ms)
    # Update the membrane voltage after 500 ms to 20 mV
    izhi.v = 20 * un.mV
    sim.run(500.0 * un.ms)

Event ports can be connected between individual cells

Poisson = CellMetaClass('./poisson.xml#Poisson')
LIFAlphSyn = CellMetaClass('./liaf_alpha_syn.xml#LIFAlphaSyn')
with Simulation(dt=0.1 * un.ms) as sim:
    poisson = Poisson(rate=10 * un.Hz, t_next=0.5 * un.ms)
    lif = LIFAlphaSyn('./liaf_alpha_syn.xml#LIFAlphaSynProps')
    # Connect 'spike_out' event-send port of the poisson cell to
    # the 'spike_in' event-receive port on the leaky-integrate-and-fire
    # cell
    lif.connect(poisson, 'spike_out', 'spike_in')
    sim.run(1000.0 * un.ms)

Network Simulations

Network simulations are specified in much the same way as Cell Simulations, with the exception that there is no metaclass for Networks (Network metaclasses will be added when the “Structure Layer” is introduced in NineML v2). Therefore, Network objects need to be instantiated within the simulation context.

with Simulation(dt=0.1 * un.ms) as sim:
    # Network objects need to be instantiated within the simulation context
    network = Network('./brunel/AI.xml#AI')
    sim.run(1000.0 * un.ms)

During construction of the network, the NineML Populations and Projections are flattened into Component Array and Connection Group objects such that the synapse dynamics in the projection are included in the dynamics of the Component Array and each port connection is converted into a separate Connection Group of static connections.

To record data, the relevant component array needs to be accessed using the component_array or component_arrays accessors of the network class. Then as in the Cell Simulations case the record method is used to specify which variables to record and the recording method is used to access the recording after the simulation.

with Simulation(dt=0.1 * un.ms) as sim:
    network = Network('./brunel/AI.xml#AI')
    # 'spike_out' is explicitly connected in the connection so it is
    # mapped to the global namespace of the flattened cell + synapses model
    network.component_array('Exc').record('spike_out')
    # State-variables of the cell dynamics are suffixed with '__cell'
    network.component_array('Inh').record('v__cell')
    # State-variables of synapses, in this case synapses from the
    # 'Inhibition' projection, are prefixed with '__<projection-name>'
    network.component_array('Exc').record('a__Inhibition')
    sim.run(1000.0 * un.ms)
exc_spikes = network.component_array('Exc').recording('spike_out')
inh_v = network.component_array('Inh').recording('v__cell')
exc_inh_a = network.component_array('Exc').recording('a__Inhibition')

Note

During the cell and synapse flattening process the names of state variables and unconnected ports will be suffixed with __cell if they belong to the population dynamics or __<my-projection> if they belong to the synapse of the a projection

Network models are simulated via integration with PyNN and therefore will run on multiple processes using Open MPI (and Open MP_ for NEST) if the calling Python script is run with mpirun/mpiexec.