DynaSim Getting Started Tutorial

Erik A. Roberts edited this page Aug 1, 2018 · 53 revisions

Overview


This tutorial and code is also available as an executable script.

Setup

First, install DynaSim.

Set where to save outputs:

output_directory = dsGetConfig('demos_path'); % or this can be any other directory
  • Note: this makes a file called 'dsConfig.txt' in the dynasim root directory which can be edited with any text editor

Make and move to the root directory where outputs will be saved:

mkdir(output_directory);
cd(output_directory);
  • ignore a warning like Warning: Directory already exists.

There are two ways to use DynaSim:

Basics of Defining and Simulating Models

Lorenz Attractor

DynaSim makes it easy to simulate arbitrary systems of ordinary differential equations. Our first example is the classic Lorenz system from general dynamical systems theory. We will use DynaSim to simulate and plot the Lorenz attractor, a set of chaotic solutions of the Lorenz system resembling a butterfly or figure eight.

To specify a DynaSim model (i.e., a system of ODEs), simply write out the equations in a string or cell array of string(s). Equations can be placed in different strings of a cell array, in the same string separated by semicolons, or using a mixture of the two techniques as with the example below.

Here are the equations that define the Lorenz attractor, written in DynaSim notation:

eqns={
  's=10; r=27; b=2.666'
  'dx/dt = s*(y-x);       x(0)=1'
  'dy/dt = r*x - y - x*z; y(0)=2'
  'dz/dt = -b*z + x*y;    z(0)=.5'
};

The function dsSimulate is called to simulate the model:

data = dsSimulate(eqns, 'tspan',[0 100], 'solver','rk4');

dsSimulate

The dsSimulate function has the form: data = dsSimulate(model, 'option1',value1, 'option2',value2,...). The options and values are paired. Here, we introduce two 'option' keys:

  • tspan: time limits for the integration in units of milliseconds [ms], given as a numeric vector of length 2.
  • solver: numerical method used to solve (i.e., integrate) the ODEs. The default is rk4, a fixed time-step 4th-order Runge-Kutta method. Other options include:

The choice of solver is a trade-off between accuracy, simulation speed, and memory usage.

When solving a model, DynaSim creates and saves a unique m-file function that will solve the model. This file is stored in a directory named 'solve'. The location of the file that solves the system (i.e., numerically integrates it) is stored in data.simulator_options and can be viewed or rerun afterwards, with the commands edit(data.simulator_options.solve_file) and (data.simulator_options.solve_file).

Every component of the model is assigned to a population (a.k.a., node), and the population name (default: 'pop1') is prepended to all variable and function names.

Simulated data can be easily plotted using the returned data structure:

figure; % make a new figure window
plot(data.pop1_x,data.pop1_z); % <-- Figure 1 in DynaSim paper
title('Lorenz equations'); % asign the figure a title
xlabel('x'); % asign the x-axis label
ylabel('z'); % asign the y-axis label

Here is the resulting plot:

lorenz attractor simulation plot

Izhikevich Neuron

While DynaSim can solve an arbitrary system of ODEs, it was created with neuroscience applications in mind. Hence, our next example is Eugene M. Izhikevich's simple model of spiking neurons, commonly referred to as the Izhikevich neuron. As Izhikevich writes,

"The model combines the biologically plausibility of Hodgkin-Huxley-type dynamics and the computational efficiency of integrate-and-fire neurons."

For more information, see page 274 of Izhikevich's book, Dynamical Systems in Neuroscience.

Simulating the Izhikevich neuron:

eqns={
  'C=100; vr=-60; vt=-40; k=.7; Iapp=70; ton=200; toff=800'
  'a=.03; b=-2; c=-50; d=100; vpeak=35'
  'dv/dt = (k*(v-vr)*(v-vt)-u+I(t))/C; v(0)=vr'
  'du/dt = a*(b*(v-vr)-u); u(0)=0'
  'if(v>vpeak)(v=c; u=u+d)'
  'I(t) = Iapp*(t>ton&t<toff) .*(1+.5*rand(1,Npop))' % define applied input using reserved variables 't' for time and 'dt' for fixed time step of numerical integration
  'monitor I'                              % indicate to store applied input during simulation
  'monitor w(t) = u.*v'                    % defining a function of state variables to monitor (note that monitor expressions follow Matlab's syntax)
};
data = dsSimulate(eqns, 'tspan',[0 1000]);
  • The resulting DynaSim data structure always contains the model state variables, time vector, and a copy of the DynaSim model structure that was simulated.
  • Additionally, functions can be recorded and returned in the DynaSim data structure if indicated using the monitor keyword.
    • Syntax: 'monitor FUNCTION_NAME'
  • t, dt, and T are special variables that can be used in model equations.
    • t represents the current time point of the simulation.
    • dt is the fixed time step used for numeric integration.
    • T is the full simulated time vector defined before simulation begins.

Plotting the simulated voltage and monitored input function:

figure; % <-- Figure 2 in DynaSim paper
subplot(2,1,1);
plot(data.time,data.pop1_v); % plot voltage
xlabel('time (ms)');
ylabel('output, v');
title('Izhikevich neuron with noisy drive')

subplot(2,1,2);
plot(data.time,data.pop1_I); % plot input function
xlabel('time (ms)');
ylabel('INput, I(t)');

izhikevich neuron simulation

Hodgkin-Huxley Neuron

The Izhikevich neuron is a simplified mathematical neuron model whose parameters do not directly map onto physiology. In contrast, the Hodgkin-Huxley neuron model is a more biophysically realistic conductance-based model with an electrical circuit analog. It has the property that each parameter has a physiological correlate that is measurable experimentally.

A Hodgkin-Huxley neuron simulation:

eqns={
  'gNa=120; gK=36; Cm=1'
  'INa(v,m,h) = gNa.*m.^3.*h.*(v-50)' % Inactivating sodium current
  'IK(v,n) = gK.*n.^4.*(v+77)' % Potassium current
  'dv/dt = (10-INa(v,m,h)-IK(v,n))/Cm; v(0)=-65'
  'dm/dt = aM(v).*(1-m)-bM(v).*m; m(0)=.1'
  'dh/dt = aH(v).*(1-h)-bH(v).*h; h(0)=.1'
  'dn/dt = aN(v).*(1-n)-bN(v).*n; n(0)=0'
  'aM(v) = (2.5-.1*(v+65))./(exp(2.5-.1*(v+65))-1)'
  'bM(v) = 4*exp(-(v+65)/18)'
  'aH(v) = .07*exp(-(v+65)/20)'
  'bH(v) = 1./(exp(3-.1*(v+65))+1)'
  'aN(v) = (.1-.01*(v+65))./(exp(1-.1*(v+65))-1)'
  'bN(v) = .125*exp(-(v+65)/80)'
};
data = dsSimulate(eqns);

% plotting
figure;
plot(data.time,data.(data.labels{1}))
xlabel('time (ms)');
ylabel('membrane potential (mV)');
title('Hodgkin-Huxley neuron')

Mechanisms

Mechanisms are predefined reusable sub-models meant to be incorporated into larger complete models. Examples of mechanisms in neuron models include ion currents, pumps, etc. Once defined, they can be easily incorporated into larger models by simply listing the name of the file containing their equations, without the need to re-write any of the mechanism equations. This greatly simplifies large model prototyping and re-configuration.

With mechanisms, one can rewrite the above Hodgkin-Huxley simulation to make it more modular and reusable, with the added benefit of simplifying the code.

Here is the equivalent Hodgkin-Huxley neuron using predefined mechanisms:

eqns = 'dv/dt = 10 + @current/Cm; Cm=1; v(0)=-65; {iNa,iK}'; % the cell {iNa,iK} specifies the mechanisms to use
data = dsSimulate(eqns);

% plotting
figure;
plot(data.time,data.(data.labels{1}))
xlabel('time (ms)');
ylabel('membrane potential (mV)');
title('Hodgkin-Huxley neuron')

% View the mechanism files:
dsEditModelFiles('iNa.mech');
dsEditModelFiles('iK.mech');

The built-in DynaSim mechanism library has models for many common ion currents (see the 'models' folder). Custom mechanisms can also be used.

Using this modular framework, let's add another ion current to our previous model in order to make it burst: the M-current (we will also use slightly different Na and K currents). As described in the spike frequency adaptation chapter of Christoph Börgers' modeling book, An Introduction to Modeling Neuronal Dynamics),

"M-currents are depolarization-activated, slowly decaying, hyperpolarizing potassium currents. The name comes from the fact that they are down-regulated by activation of the muscarinic acetylcholine receptors (the M stands for muscarinic)."

As Izhikevich mentions in his chapter on bursting, an M-current provides negative feedback to stop spiking temporarily, thus creating a burst.

Here is an example of a bursting neuron model using three predefined current mechanisms:

eqns='dv/dt=Iapp+@current; {iNaF,iKDR,iM}; Iapp=5; gNaF=100; gKDR=5; gM=1.5; v(0)=-70';
data = dsSimulate(eqns, 'tspan',[0 200]);

% plotting
figure;
plot(data.time,data.(data.labels{1})) % <-- Figure 3 in DynaSim paper
xlabel('time (ms)');
ylabel('membrane potential (mV)');
title('Hodgkin-Huxley-type intrinsically bursting neuron')

bursting neuron simulation

Predefined Populations

In addition to predefining the model mechanisms, it is also possible to predefine the population (i.e., a model with a collection of mechanisms). This permits the modular combination of different predefined neuron types for example. Predefined populations are stored in text files (e.g., 'IB.pop', using the preferred pop extension), and simulated by passing the file name to dsSimulate (e.g., dsSimulate('IB')) or by equating population equations to it in the DynaSim specification structure (see below).

Example of the same bursting neuron using a predefined population:

eqns='IB'; % file name of predefined model of intrinsically bursting (IB) neuron
data = dsSimulate(eqns, 'vary',{'IB','Iapp',5}, 'tspan',[0 200]);

% plotting
figure;
plot(data.time,data.(data.labels{1}))
xlabel('time (ms)');
ylabel('membrane potential (mV)');
title('Predefined Intrinsically Bursting neuron');

% View the predefined population file:
dsEditModelFiles('IB.pop');

For more information, see the pages on using predefined populations and making predefined populations.

Defining and Exploring Network Models

The above approach is sufficient for building single-compartment models with arbitrary complexity. However, larger multi-compartment and network models require defining multiple compartments or cell types and connecting them. DynaSim simplifies building these models as well. Behind the scenes, DynaSim always converts the user-supplied equations into a DynaSim "specification" structure that can be used to specify any model at a low or high level of abstraction. Learning to define specification structures can greatly facilitate (1) constructing large network models and (2) parameterizing control of model parameters and mechanism lists in Matlab scripts. The latter advantage provides a higher degree of control that is recommended for conducting computational research without the need to frequently manipulate equation strings.

In the above examples, DynaSim first splits the user-supplied population equations into mechanism_list, parameters, and equations variables; it then stores those separately in a structure (the "specification" structure) and assigns a name and size to the (implicit) population. Practically speaking, this approach requires the user to store the above model information in a specification structure (instead of an equation string or cell array of strings). The specification structure allows multiple populations to be defined (e.g., different compartments or cell types) and connected (e.g., by synapse mechanisms from a 'source' to a 'target').

Schema for DynaSim specification structure (select fields shown):

  • .populations(i): contains info for defining independent population models (required)
    • .name : name of population (default: 'pop1')
    • .size : number of elements in population (i.e., # cells) (default: 1)
    • .equations : string listing equations (required)
    • .mechanism_list : cell array listing mechanisms (default: [])
    • .parameters : parameters to assign across all equations in the population. provide as cell array list of key/value pairs (default: []) {'param1',value1,'param2',value2,...}
  • .connections(j): contains info for linking population models (default: [])
    • .direction : connection direction (source -> target) where 'source' and 'target' are population names
    • .mechanism_list : list of mechanisms that link two populations (required)
    • .parameters : parameters to assign across all equations in mechanisms in this connection's mechanism_list. (default: [])

See dsCheckSpecification in the Function Reference for more details.

Pyramidal Interneuron Network Gamma Rhythm (PING)

One of the simplest biophysical neural network models using Hodgkin-Huxley neurons is the Pyramidal Interneuron Network Gamma Rhythm (PING). As described in the 'The PING Model of Gamma Rhythms' chapter of the Börgers modeling book,

"When populations of excitatory and inhibitory neurons are synaptically connected, oscillations often emerge. The reason is apparent: Activity of the excitatory neurons...generates activity of the inhibitory neurons (I-cells). The activity of the I-cells causes the activity of the E-cells to cease transiently, and when it resumes, the E-cell population is closer to synchrony...In particular, brain oscillations with frequencies of about 30–80Hz are thought to arise in this way in many instances. Oscillations in the 30–80 Hz frequency range are called gamma oscillations or gamma rhythms. They are seen in EEG traces and in local field potentials, and are correlated with perception, attention, and memory."

For more information, checkout the Scholarpedia article on fast oscillations, one of the early publications about PING, or a more recent review paper about gamma oscillations.

The PING network model is incredibly easy to simulate in DynaSim, especially when using predefined mechanisms as described previously. Using predefined populations simplifies this further.

How to setup Sparse Pyramidal-Interneuron-Network-Gamma (sPING):

ping-architecture

% define equations of cell model (same for E and I populations)
eqns={
  'dv/dt=Iapp+@current+noise*randn(1,N_pop); Iapp=0; noise=0'
  'monitor iGABAa.functions, iAMPA.functions'
};
% Tip: monitor all functions of a mechanism using: monitor MECHANISM.functions

% create DynaSim specification structure
s=[];
s.populations(1).name='E';
s.populations(1).size=80;
s.populations(1).equations=eqns;
s.populations(1).mechanism_list={'iNa','iK'};
s.populations(1).parameters={'Iapp',5,'gNa',120,'gK',36,'noise',40};
s.populations(2).name='I';
s.populations(2).size=20;
s.populations(2).equations=eqns;
s.populations(2).mechanism_list={'iNa','iK'};
s.populations(2).parameters={'Iapp',0,'gNa',120,'gK',36,'noise',10};
s.connections(1).direction='I->E';
s.connections(1).mechanism_list={'iGABAa'};
s.connections(1).parameters={'tauD',10,'gGABAa',.1,'netcon','ones(N_pre,N_post)'}; % connectivity matrix defined using a string that evalutes to a numeric matrix
s.connections(2).direction='E->I';
s.connections(2).mechanism_list={'iAMPA'};
s.connections(2).parameters={'tauD',2,'gAMPA',.1,'netcon',ones(80,20)}; % connectivity set using a numeric matrix defined in script
% Simulate Sparse Pyramidal-Interneuron-Network-Gamma (sPING)
data=dsSimulate(s);

dsPlot(data); % <-- Figure 4 in DynaSim paper

ping-dynamics

To view the connection mechanism file, run:

dsEditModelFiles('iAMPA.mech');

Mechanisms iAMPA and iGABAa contain a fixed variable netcon storing a default connectivity matrix (all 1s for all-to-all connectivity). The user can specify the connectivity matrix from a source to target by setting netcon in the .parameters cell array. The dimensions of netcon should be [N_pre x N_post] where N_pre (reserved variable) is the size of the (presynaptic) source population and N_post (reserved variable) is the size of the (postsynaptic) target population. netcon can be set as a numeric matrix or a string that evaluates to a numeric matrix of the correct dimensions.

See Specifying connectivity matrices for more details.

Running sets of simulations by varying model parameters

One of the most fun and useful ways to interact with a DynaSim model is to explore how its behavior changes as the various parameters (and even mechanisms) making up the model specification change. DynaSim makes this easy to do, with the vary argument to dsSimulate. The value of vary is a cell array, indicating which model variables (e.g. parameters or mechanisms) to vary, the values these model variables should take, and the object (population or connection) whose variable should be varied. There are two possible ways to specify the vary array.

Setting specific, single values: vary={{object, variable, value1},{object, variable, value2},...}

  • this is useful for simulating an arbitrary set of parameter values.

Setting ranges of values: vary={object(s), variable(s), values; ...}

  • this is useful for varying parameters systematically, in combination and in parallel (described later).

For instance, to vary parameter gNa, taking on values 100 and 120, in population E, set vary to {{'E','gNa',100},{'E','gNa',120}} (syntax 1), or set vary to {'E','gNa',[100 120]} (syntax 2). To additionally vary gAMPA in the connection mechanism from E to I, set vary to {'E','gNa',[100 120];'E->I','gAMPA',[0 1]}. Mechanism lists and equations can also be varied. (see dsVary2Modifications in the Function Reference for more details and examples).

Below is an Izhikevich study of neuro-computational properties (using Syntax 1), based on http://www.izhikevich.org/publications/izhikevich.m:

eqns={
  'a=.02; b=.2; c=-65; d=6; I=14'
  'dv/dt=.04*v^2+5*v+140-u+I; v(0)=-70'
  'du/dt=a*(b*v-u); u(0)=-20'
  'if(v>=30)(v=c;u=u+d)'
  };
P='pop1'; % name of population
vary={
  {P,'a',.02; P,'b',.2 ; P,'c',-50; P,'d',2;  P,'I',15} % tonic bursting
  {P,'a',.01; P,'b',.2 ; P,'c',-65; P,'d',8;  P,'I',30} % spike frequency adaptation
  {P,'a',.02; P,'b',.2 ; P,'c',-65; P,'d',6;  P,'I',7}  % spike latency
  {P,'a',.03; P,'b',.25; P,'c',-52; P,'d',0;  P,'I',0}  % rebound burst
  {P,'a',1;   P,'b',1.5; P,'c',-60; P,'d',0;  P,'I',-65}% bistability
  {P,'a',.02; P,'b',1  ; P,'c',-55; P,'d',4;  P,'I',1}  % accomodation
  {P,'a',-.02;P,'b',-1 ; P,'c',-60; P,'d',8;  P,'I',80} % inhibition-induced spiking
  {P,'a',-.026;P,'b',-1; P,'c',-45; P,'d',0;  P,'I',70} % inhibition-induced bursting
  };
data=dsSimulate(eqns, 'tspan',[0 250], 'vary',vary);
dsPlot(data);

To vary two parameters (run a simulation for all combinations of values) using Syntax 2:

vary={
  'E'   ,'Iapp',[0 10 20];      % amplitude of tonic input to E-cells
  'I->E','tauD',[5 10 15]       % inhibition decay time constant from I to E
  };
data=dsSimulate(s, 'vary', vary, 'parfor_flag',1);
dsPlot(data,'plot_type','rastergram'); % <-- Figure 5 in DynaSim paper
dsPlotFR(data); % examine how mean firing rate changes with Iapp and tauD

ping-vary

Saving simulated data

A set of simulations, analyses, and plots deriving from a common base model are collectively called a DynaSim "study". All results for a given study are saved in a directory corresponding to the study_dir variable. This directory contains subdirectories 'data' (simulated data and results derived from analyzing the data), 'plots', and 'solve' (containing m- and mex-files that were used for all simulations of the study).

How to: set save_data_flag to 1 and (optionally) study_dir to /path/to/outputs

Saving data from a single simulation:

Example using the previous sPING model:

data = dsSimulate(s, 'save_data_flag',1, 'study_dir','demo_sPING_1');

Loading the saved data from disk:

data = dsImport('demo_sPING_1');

Save data for a set of simulations:

% Specify what to vary
vary = {'E','Iapp',[0 10 20]}; % vary the amplitude of tonic input to E-cells
data = dsSimulate(s, 'save_data_flag',1, 'study_dir','demo_sPING_2', 'vary',vary);

Load and plot the saved data:

data_from_disk = dsImport('demo_sPING_2');
dsPlot(data_from_disk);
dsPlot(data_from_disk,'variable','E_v');

Tip for saving large data files:

  • By default data is saved in compatible mode between Matlab and Octave (matCompatibility_flag set to 1, i.e., data is saved in '-v7' mat format).
    • Unfortunately, '-v7' mat format is not able to save variables > 2GB.
    • If compatible saving fails, data is saved in '-v7.3' format (Matlab), or in '-hdf5' format (Octave). This allows that data can be stored in all its integrity.
  • However, if any var > 2GB is anticipated, we strongly recommend to set matCompatibility_flag manually to 0 to speed up data saving. Example using matCompatibility_flag set to 0: data=dsSimulate(s,'matCompatibility_flag',0,'save_data_flag',1,'study_dir','demo_sPING_1_varlt2GB');

For more information about the DynaSim data structure, see the Data Structure Reference


→ Dynasim Graphical User Interface (GUI) Tutorial

You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.
Press h to open a hovercard with more details.