DSSRT: A dynamical systems reduction & modeling software tool

DSSRT: Dominant-Scale System Reduction Tool

User documentation

(c) Robert Clewley, Center for Biodynamics, Boston University, 2002-2004.

 

 

Contents

 

1. Version information and bug list

2. Introduction to dominant-scale analysis

3. Getting started

            3.1 Setting up the software on your computer

            3.2 Running the program

            3.3 Display a functional network diagram

            3.4 Interpreting the graphical output

            3.5 Basic navigation

            3.6 Compute a transition sequence

            3.7 Other documentation

4. About the included demos

            4.1 Single excitatory Hodgkin-Huxley cell with no synaptic input

            4.2 Excitatory and Inhibitory two-cell network

            4.3 Single excitatory cell with inhibitory synaptic input

            4.4 A stellate cell

            4.5 Fitzhugh-Nagumo oscillator

5. Overview of DSSRT Matlab functions

            5.1 Functional Networks (FuncNet.m)

            5.2 Transition Sequences

            5.3 Main window interface and commands (StepNet.m)

            5.4 Specification of differential equations

            5.5 Reduced regime determination (RegimeDet.m)

            5.6 Domain of validity / attractor estimation (AttEst.m)

6. Definitions, and main parameters

            6.1 Internal, external variables, and observables

            6.2 Passive and non-passive variables

            6.3 Influence strength (Psi)

            6.4 Dominance scale threshold (sigma)

            6.5 Active vs. potentiated variables

            6.6 Re-sampling the data (nOut)

            6.7 Time scale threshold

7. Understanding the .cfg specification format

            7.1 Overview

            7.2 Correct order of system specifications

            7.3 Full syntax for DSSRT specifications

            7.4 Additional node states and state maps

            7.5 Links

            7.6 User-defined .m functions for use in DE specification

8. Reduced Dynamical Regimes

            8.1 Overview

            8.2 Using the RegimeDet function

            8.3 More information about the algorithm

9. Domain of Validity / Attractor Estimation

            9.1 General guidelines

            9.2 How to calculate an attractor estimate

            9.3 'Bolstering' the original epoch/regime set

            9.4 How DSSRT deems what are qualitatively important changes in Acts

10. Troubleshooting and other practical tips

11. Future applications, contributions

12. Glossary of new terms

 

 

1. Version information and bug list

 

This is DSSRT Version 1.12, May 2004

DSSRT was built on Matlab v6 R13

 

Please send any comments and bug reports to rclewley@bu.edu

 

Known compatibility issue: Pre-release 13 versions of MATLAB do not support the predefined boolean constants used in this program.

Known bug: When 'P' command changes param file, time bar is not correctly restored and subsequent movement in time is not displayed with the vertical marker (nor are the markers).

Known bug: Occasionally, after switching windows, Matlab momentarily redisplays the DSSRT window with the time/position information incorrectly placed, and continues to update it in the correct place.  However, the errant text remains on the window until the system is reset.

Possible bug: Lowering 'low-resolution multiple' parameter for AttEst too much makes computed "attractor" get smaller.

Known bug: If variable data .dat files are renamed, then reloading _FN.dat files referring to the previous names prevents the user from being able to change .dat files using 'O' command.  The command will fail but DSSRT recovers.

Known bug: When functional net re-calculated, most recently calculated transition sequence still thinks it's "fresh" and command 'X' suggests its data is still being displayed in the main window, which it is not.

 

Version histories of individual components are listed in the associated .m files

 

 

 

2. Introduction to dominant-scale analysis

 

Systems of oscillators described by ODEs arise commonly in the natural sciences, often with multiple time-scales.  State-dependent coupling can add to the complexity, introducing new time-scales.  These time-scales may not be explicit in the equations.

 

DSSRT embodies a computational technique known that can be used to reduce the study of such networks, near a known orbit, to a set of low-dimensional approximate models.  This method focuses on dominant scales in the relationships between variables during transient (far from equilibrium) dynamics.  The technique adds rigor to intuitive reduction techniques that are ubiquitous in modeling high-dimensional coupled systems.  In particular, it quantifies the robustness and parametric dependence of coherent temporal activity in more detail than a return-map stability analysis of periodic orbits, because it provides information along the entire length of the trajectory.

 

DSSRT analyzes data that has already been numerically computed using a differential equation solver such as those built into Bard Ermentrout's XPP software, Neuron, GENESIS, or any fixed time-step ODE solver from which you can export variable data to a file.

 

See the PDF documentation and the explanatory sections below, for more information about the theory behind the dominant-scale technique.

 

 

 

3. Getting started

 

3.1 Setting up the software on your computer

 

Make sure you comment / uncomment the appropriate lines in the file DSSRT_useroptions.m depending on whether you are running DSSRT under Unix/Linux or MS Windows, and change the directory settings according to where you've put your main DSSRT folder.  Also, you can make sure that the `START_VERBOSE` option is set to `true` if you're new to DSSRT and wish to receive more system information in the console.

 

In Matlab, change directory to the DSSRT folder you are using.  Type `DSSRT` at the prompt and enter the name of a directory containing network setup information (see Section 4 regarding the included demos).  Then type 'H' or '*' for further help within the main window that appears.

 

Variable data files (.dat) are assumed to have been generated at a fixed time-step in this version of DSSRT.  You can generate these files using any fixed time-step ODE solver of your choice, provided you can output the data to a space separated text file in a column order that corresponds to that specified in a network setup (.cfg) file for your system.  See the later section for details on the specification format, and the included demos for more information.  Files included with the DSSRT system provide examples for use with the XPP software (www.pitt.edu/~phase).

 

Compilation of certain .m files (those known to exhibit speed improvements on compilation) can be done automatically the first time DSSRT is run, if a compiler is available to your Matlab installation.  If you know that your Matlab compile is configured then set FORCE_NOTCOMPILE = true in the DSSRT_useroptions.m file, otherwise DSSRT will continue without compiling certain functions.  Compiling other .m files by hand will most likely not improve the speed of the program (and will make some parts of the system crash on execution).  I have tested the execution speeds of all components to optimize them either for the Matlab performance enhancement mode when interpreted (using the JIT system), or for compilation to MEX files, whenever they perform substantially better in those forms.  This was done on a Windows NT system.  If there is a problem compiling the the code on your system, the compile script may ask you to switch off the compile option in order for you to continue.  The loss of performance would not be great in that situation, as the majority of the DSSRT functions work faster uncompiled when using the interpreter's internal speed enhancement (the Just In Time compiler).

 

 

3.2 Running the program

 

Start the program by changing directory within Matlab to the installation directory of DSSRT.  To study a system that has already been set up (in its own sub-folder of this directory), say 'E_single_cell', type DSSRT E_single_cell at the Matlab prompt. This should initialize the program and load the model-specific configuration.  Alternatively, typing DSSRT at the prompt by itself will bring up a dialog box in which you can specify the name of a system's folder.

 

 

3.3 Display a functional network diagram

 

A functional network diagram (defined in Section 6.1) is the most fundamental thing that you need to have in memory in order to use DSSRT to analyze your system.  In most of the included demos (see Section 4 below for details), there is a pre-computed functional network available for a pre-computed orbit of that system.  That orbit is usually a stable limit cycle.  To find out if the required file is present, press 'E' in the main window.  By default, the filename in the dialog box will appear as 'demo'.  Click OK, and if it is present (as the file 'demo_FN.mat' in the system's folder) then it will be loaded.  Such a file is present for the 'E_single_cell' demo, for instance.

 

Alternatively, if the basic parameters were set up correctly in the .par file for the model system, and if the start and end times are set correctly (using commands '1' and '2'), then a new functional network can be computed.  This is done using command 'U', which uses whatever resampling time-step is specified (by command '4'), and the current dominant-scale threshold (set using command '3').  The dominant-scale threshold is explained in Section 6.4, and re-sampling is explained in Section 6.6.

 

 

3.4 Interpreting the graphical output

 

Firstly, you should familiarize yourself with some of the background theory in Section 5, and possibly in the theory document (see Section 3.4) for the more mathematically minded.

 

3.4.1 The ODE model as a network

 

Your ODE model system is represented graphically in the main window as a network of interconnected nodes.  The .cfg file in the system's directory specifies the layout of the network in the window, among other things.  Before a functional network has been loaded or computed, only the nodes of the network will be shown.  These contain the names of the associated variables in the system.  Small vertical "slider bars" will also be shown nearby to each node.  These will show the current value of that variable using a small horizontal bar.  If there is more than one slider beside a node, the additional ones will show the values of hidden ("internal") variables associated with that node, that are not visible to the rest of the network.

 

3.4.2 Color-coded arrows

 

After computing or loading a functional network various red and green arrows should appear on the diagram, which connect colored circles around the variable names shown.  The arrows and the colors give information about which variables are "dominant" in the dynamics in a neighborhood of the computed orbit, at the time shown above the diagram.  See Section 6.5 for details of the meaning of these graphics.  The time-line is shown above the digital time / position readout, with a blue cursor marking the current time.  Raw variable data associated with the current functional network diagram is plotted in the top pane of the main window for any variable of the system, and can be changed using command '6'.

 

But in the simplest terms, at a given time t, a red arrow from a variable 'x' to a variable 'y' in the functional network diagram indicates that x is having a "significant effect" on the time-course of y at that moment.  A green arrow indicates that x does not have a significant effect on y, but that x could do so if it had a different value.  Because DSSRT is told what range x can take (often it is a [0,1]-valued "gating" variable) this is usually telling you that if x was more highly activated then it could become a dominant variable on the dynamics of y.  In the absence of any arrow between variables, there is either no physical connection present or at this time x could not take any value in its known range that would cause it to have a significant effect on y.  The meaning of "significant effect" is described in Sections 5, 6.4, and 6.5, and in the theory documentation.

 

 

3.5 Basic navigation

 

To navigate through the diagram as time changes along the loaded orbit, use the 'N' and 'B' keys for one time-step movements forwards and backwards (at the re-sampled time-step used to compute the diagram data), or ',' and '.' for ten time-step movements at once.  Pressing shift with the latter two commands (i.e. pressing '<' or '>') will jump the cursor to the start and end times of the current diagram data (the very ends of the time-line).

 

3.5.1 Time markers

 

Time markers are useful for several functions within DSSRT.  Think of them as defining a "selection" of time on which you wish to perform an action, as you would do by selecting text in a document editor that you wish to spell-check.  Two markers can be set at any one time: one for a left end-point, and one for a right end-point.  To mark a left time end-point, navigate to that time and press 'J'.  This positions a small red tick on the time line.  Navigate to the end time of your chosen interval and press 'K'.  You can set the right marker before you set the left marker, provided you don't try to set a left marker at a later time position than a right marker that you have already set (and vice versa).  You can change the marker positions by re-pressing the appropriate key at a different time position in the diagram.

 

Once set, you can navigate directly to a left marker by pressing ';', and to a right marker by pressing ':' (i.e. SHIFT plus the semicolon key).  Notice that when you are exactly at an end-point either [L] or [R] will appear at the end of the digital time/position readout.

 

Positions of the time-markers are saved when the current functional network diagram is saved using command 'W', and they are over-written when one is loaded using command 'E' (or when a functional network is recomputed).

 

 

3.6 Compute a transition sequence

 

The next step of analysis using DSSRT is to compute an event transition sequence ('TS') for this functional network.  Usually you will want to restrict this to events involving a particular variable of focus (typically a voltage variable for neural models).  First, check the variable of focus using command '5'.  Mark out a time interval within which to compute the TS, as described above, or you can use one loaded in a demo using command 'E'.  Often you will want the selected interval to be one whole cycle period of a limit cycle, but it can be any sub-interval of the orbit you wish.

 

Now compute the transition sequence by pressing 'C', and press OK to restrict the sequence to the focused variable.  (In very large systems, un-restricted transition sequences are slower to compute and do not highlight features specific to the local dynamics.)  When this is done you have an opportunity to name your TS and add two additional data to be associated with it.  For now just use the default values (by pressing RETURN repeatedly).  Small red dots will appear above the time-line, indicating the positions of events found within the marked time interval.  You can navigate through these event times cyclically by pressing the SPACE bar repeatedly.

 

You are now ready to try out other commands that aid your analysis of the structure inherent in your ODE system.  Most of them require that a TS has been calculated before they can run.  Use the help command 'H' to list them.  The major ones have separate sections in this document explaining their use, e.g. the regime determination and attractor estimation / domain of validity utilities.

 

Please let me know if you have any problems or suggestions.

 

 

3.7 Other documentation

 

The `full` PDF documentation referred to in this document refers to a theory paper, DomScaleTheory.pdf, and a slide presentation that includes information about the attractor estimation algorithm AttEst.m.  Beware that some of the notation in those documents may not yet be entirely consistent with that adopted here.

 

 

 

4. About the included demos

 

4.1 Single excitatory Hodgkin-Huxley cell with no synaptic input

 

The `E_single_cell` demo is the most basic demonstration of the functional network (event, epoch, regime) analysis applied to a neural model.  The cell regularly fires with a period just under 50ms due to a drive current.  The regimes determined from this model agree exactly with those derived by formal methods in a paper by Suckley and Biktashev (see the PDF theory document for details).

 

 

4.2 Excitatory and Inhibitory two-cell network

 

The demo 'EI_demo' consists of two hippocampal neurons E and I (modelling a pyramidal cell and an interneuron) coupled through chemical synapses.  The EI.ode file for XPP contains the set-up for this network.  The principal behaviour of this system for the parameter regime given is a synchronous state in the gamma frequency regime (30-100Hz), in which the I-cell fires twice per cycle (known as a `doublet`) compared to the E-cell, which undergoes a temporary suppression every half cycle due to the incoming strong inhibition.  Various data sets of the variables produced by XPP are listed, as the I->E coupling variable g_ie is varied around its original value.  The saved functional network demo_FN.mat contains one period of data for a functional network for one of those data sets.

 

 

4.3 Single excitatory cell with inhibitory synaptic input

 

The demo 'E_demo' is a stripped-down version of 'EI_demo', for a simpler view of a single neuron.  In particular it more explicitly shows the roles of the action potential currents.  The inhibitory cell has been replaced with a periodic inhibitory synaptic input.  A sample variable data file, and a functional network with markers set for one complete period, is included.

 

 

4.4 A stellate cell

 

'Stellate_single_cell' is a demonstration of a stellate cell incorporating two slow H-type currents, and a persistent potassium current.  The analysis of this system is the subject of an upcoming paper in collaboration with Horacio Rotstein.

 

 

4.5 Fitzhugh-Nagumo oscillator

 

'FHN_demo' takes a look at applying the dominant-scale technique to a non-neural equation.  The FitzHugh-Nagumo equations are set up here in a classic relaxation oscillator regime.  In this case, there is a 'synaptic'-like input added.  In one form of the two FHN systems provided, the dummy `w` variable is dropped from the analysis (see the XPP .ode files for details).  The AttEst function does not provide reliable information about FHN orbits at this time, because this example system involves non-dissipative dynamics, and AttEst does not yet deal with this class of ODE system.

 

 

 

5. Overview of DSSRT Matlab functions

 

Please refer to the full documentation PDF file for more detail about the theoretical concepts discussed briefly below.

 

 

5.1 Functional Networks (FuncNet.m)

 

The matlab function FuncNet.m is run using key command 'U'.  It calculates subsets of the set of observable variables in the system for each time-step over a selected time interval.  At each time step the subset contains the variables are most dominant in affecting the instantaneous target of the system.  This is judged according to a `dominant-scale threshold` i.e. a scale of influence of one variable on another.  The set of all the unique time-ordered subsets can be viewed as a sequence, defining the states of the functional network's diagram.  The dominance of variables on others is sometimes non-intuitive.  Our definition of dominance is given in Section 6.3.

 

Instead, we define a measure of the effective degree of influence that a change in an `input` variable has on the instantaneous `target` of the observable is ranked against all other forces exerted by the `input` variables to the RHS of the observable's equation, and the ranked set (in descending order) of variables having the largest such influence is truncated after their size becomes a sufficiently small fraction (the dominant-scale threshold) of the most influential variable.

 

More precisely:  A type of sensitivity calculation determines the temporal pattern of the .dominant. variables over an orbit (e.g. an attractor), according to the definition of scale of influence.  This partitions the phase space in the temporal direction into epochs, according to a user-defined dominant-scale threshold.  The subsets of dominant variables in each epoch define the transition sequence.  Within each epoch, the local basin of attraction to the invariant manifold is approximated by a lower dimensional manifold corresponding to only the locally dominant variables.

 

 

5.2 Transition Sequences

 

A Transition Sequence (abbrev. TS) is a set of states of a functional net that describe all changes in the network's dominant variables over some time interval, i.e. excluding the contiguous non-changing states that typically occur in the original functional network's state sequence.  Each transition is called an event, and marks the beginning of a new epoch, during which a certain set of variables are `active` (most dominant), and those remaining are one of two types of inactive variable: potentially active, or non-potentially active.  The first set is known as the Actives set, the second is the Potentials set, and the remaining variables are only implicitly defined by exclusion from the first two sets.  (Additionally, the user has the option to distinguish changes in the order of the Active sets as indicating a transition, using command '7'.)  The most recently computed TS is known as the 'fresh' TS.  Any saved TS can be loaded in the current version of DSSRT, regardless of whether it was computed for the same parameters as currently set, or even for the same variable data set.  Until this irritating matter is addressed in a future version of the program, a temporary feature is that many functions requiring a TS specify that the TS be fresh.  This guarantees it to be consistent with the whole ODE system currently set up.

 

A set of related epochs make up a reduced dynamical regime, within which a certain set of equations for dynamical variables is most appropriate for analysis, according to standard multiple-scale analysis techniques (see Section 9 for automated features).

 

Different transition sequences can be calculated for different saved data sets that contain different time-evolutions of the ODE system's variables, e.g. arising from different initial conditions or parameter changes in the original ODE system.  A notion of `distance` from a transition sequence (e.g. a periodic cycle) can be defined in terms of the sum of the errors in comparing a similar sequence of the same length to the original, according to a particular algorithm for making similar sequences the same length.  This leads to speculation of the meaning of considering TS's in the same fashion as discretized dynamics, and thus of the stability and bifurcation of transition sequence cycles.

 

5.2.1 Viewing and comparing transition sequences

 

Command 'V' compares transition sequences of equal length.  If two TS's are not originally of equal length the command 'Z' will attempt to pad the shorter one in an intelligent way so that the positional errors between them are minimized.  Command '#' controls the fussiness of the padding algorithm, although it is still in development.

 

Key to the horizontal `bar` display of viewed transition sequences (using command 'V'):

·             &nbs p;       Contiguous blue lines indicate original portions of the sequence.

·             &nbs p;       Broken lines in black indicate portions that were padded by the sequence matching algorithm in order to minimise the number of errors between sequence members over the length of the sequences, and to make the sequences the same length.

·             &nbs p;       Remaining errors in sequences generated by the matching algorithm (see StepNet) for a sequence are shown by red markers at their sequence positions.

 

 

5.3 Main window interface and commands (StepNet.m)

 

StepNet.m provides the principal user interface to control the derivation and comparison of functional network diagrams and their transition sequences.

 

 

5.4 Specification of differential equations

 

Differential equation specification is achieved using the DEQNS, GAMMA1, GAMMA2, DEPARS specs in the .cfg setup file.

 

Take the example of a Hodgkin-Huxley neuron like that in the included demo `E_demo`.  Typically, a quantity like `mtau_recip` governing the relaxation time of the m variable is not a constant value multiplying the gating variable m in a right-hand side, but a self-contained activation function.  Thus, when they are identified as such functions, argument 3  = `u1` indicates what `mtau_recip` is a function of (note only 1 argument is allowed for those functions).  In such cases the power in argument 4 is ignored (and so typically set to 1).  `mtau_recip`, etc. are known to DSSRT as functions of other variables, and not just constant parameters, because the INPUTS specification for the associated gating variables are non-empty.  (DSSRT takes these names and uses the MATLAB interpreter eval to execute functions with these names, and so .m files having those same names must reside in the ODE system's folder.)

 

For each differential equation, there should be no more than INPUTS(eqn) number of GAM1TERMs + GAM2TERMs.  Also, for fewer than that number, those present in INPUTS that are not in a Gamma set may only be internal variables.  See the examples provided in the demos.

 

 

5.5 Reduced regime determination (RegimeDet.m)

 

This function takes an event transition sequence and consolidates the epochs into a reduced set of regimes, which are more amenable to analysis.  The function does this by looking for commonality between epochs, according to a simple set of rules.  It also explicitly identifies fast and slow variables with respect to the primary focused variable, and predicts appropriate quasi-static bifurcation parameters for each regime.  Section 8 contains more information.

 

 

5.6 Domain of validity / attractor estimation (AttEst.m)

 

Although still in development, AttEst.m uses an efficient algorithm to calculate the domain of validity for a transition sequence or sequence of regimes.  It computes a cross-section of the basin of attraction around a known orbit (under some constraints), so it can be used to study the orbit's robustness to perturbations.  This can be done for any of the external variables.  More information is given in Section 9.

 

 

 

6. Definitions, main parameters

 

Please refer to the full documentation PDF file for more detail about the theoretical concepts introduced below.

 

 

6.1 Internal, external variables, and observables

 

These are subjective definitions which depend on the interpretation of the structural components of the system.  On the RHS of an ODE for a certain variable `V`, if certain terms are considered due to an input `external` to the sub-system that V represents, then the variable contributing to that input is an external variable.  Conversely, if inputs to V arise from dynamics that are considered internal to V's subsystem, then that input's variable is an internal variable.  E.g. for the Hodgkin-Huxley equations, we might consider the sodium and potassium spiking currents to be internal to the cell, and so variables m, n, and h are internal variables, whereas a synaptic input from another cell means that it's variable s is an external variable.  If we consider the equation for the synaptic variable existing in a separate sub-system to the soma then the pre-synaptic membrane potential is an external variable since it inputs to the equation for s.  The synaptic equation has no internal variables in the absence of adaptive processes (like depression, facilitation).  An `observable` is merely a synonym for an external -- it is a variable that can be `observed` by other sub-systems.

 

 

6.2 Passive and non-passive variables

 

We call variables that have an associated differential equation non-passive (or dynamic) variables, and those that do not are passive (non-dynamic).  Passive variables play the same role as non-passive variables in input terms to a differential equation, but they might be truly formal (i.e. actually non-varying) or be an externally-varied signal.  This terminology is used to distinguish this issue from that of dominance, where the term `active` is used (see section 6.5).

 

 

6.3 Influence strength (Psi)

 

Loosely speaking, the influence strength of an input to some variable V's ODE is the amount of control that the input has over the instantaneous target of the variable.  The usefulness of this quantity in DSSRT relies upon certain assumptions about the ODE for V, as discussed in the full documentation (PDF file).  In particular, that trajectories V(t) contract sufficiently quickly towards the instantaneous target, that the motion of the target is a reasonable approximation to the actual trajectory, in some specific sense.

 

 

6.4 Dominance scale threshold (sigma)

 

This is the critical value defining how exclusive or inclusive you wish your sets of active inputs to each observable variable to be.  It is set using command '3' (which also sets the time scale threshold which we discuss later).  Mathematically it defines the Order 1 scale of variables in the system, and is akin to the reciprocal of the ubiquitous small parameter epsilon used in multiple-scale analysis.  Formally, sigma is the reciprocal of the fraction of the largest ranked influence strength (for any observable), such that any other input's influence strength less than this value will not be included in the set of active inputs acting on a the observable variable, at one moment in time (see the full documentation for a mathematical definition).

 

Increasing sigma away from 1.0 will loosen the tolerance on the strength of influence of other candidate-active inputs to the observable, and allow weaker-influencing inputs to become considered `active'.  Decreasing towards 1.0 will tighten the tolerance on influence strength and ignore more and more of the weaker-influencing inputs.  Values of 1 -> 5 give different intuitive perspectives on what is `dominant` at any given time, at least for ODEs describing neural membrane channel dynamics.  Values from 2 -> 3.5 tend to be best for accurate attractor estimates in neural systems.

 

Note that in the full documentatin, the sigma scale threshold is currently referred to only by its reciprocal, and then it has the symbol epsilon.

 

 

6.5 Active vs. Potentiated (a.k.a. potentially active) variables

 

Green-coloured links in the diagram (and certain special colours for neurons' internal variables) show when an input has the potential to be in the set of most dominant variables acting on some observable, if it is not currently in the set of actives.  This can only be defined using knowledge (or estimation) of the maximum absolute value of that input, and is somewhat sensitive to the actual choice of that value (unlike for the actives set itself, which is insensitive to this choice provided it is only erroneous on the conservative side).

 

The procedure for determining this at any time step is to find the maximum value of the influence strength for the variable's gating coefficient s_i over the range 0 -> 1, given the actual current values of all other variables.  This value is then compared to the current set of actives, and if it would belong to that set (using the discrimination of the scale threshold), then that input is considered potentiated.  This information is less useful than knowledge of the actives set, since it is only a guide to the potential influence that variable could be having at some time, i.e. if it were to exactly reach it's point of maximum influence!

 

One helpful trick with making the most from 'potentials' is to do the following.  If you are studying behaviour near a known orbit, observe the range of variable values taken on that orbit and in perturbations from it that you might be interested in.  Then, restrict the BOUNDS for that variable accordingly in the .cfg specification (although leaving a small safety margin is recommended).  That way, the accuracy of the green 'potentiated' arrows is improved, having been made more relevant to this dynamical regime.  Be warned however that significant changes in system parameters may invalidate these bounds and you will have to make sure you adjust them accordingly.

 

 

6.6 Re-sampling the data (nOut)

 

If the sample resolution of the variable data in the .dat file was set too small during the original numerical integration of the ODEs, DSSRT may take a long time to load files and perform major calculations.  Use the nOut setting (command '4') to re-sample data every nOut steps when calculating the functional network to improve speed and efficiency.  However, beware that quantization artifacts can lead to missed or lumped transitions in the state sequence, which can make comparison with other data sets ambiguous.  The visible resolution of the variable displayed in the top pane of the main window is also determined by this parameter, and can be used as an indication of whether it is appropriately set for a particular dataset.

 

 

6.7 Time scale threshold, `fast` and `slow` variables

 

This is also set using command '3', along with the dominance scale threshold.  It is only used by the attractor estimation algorithm (see Section 9) and for generating reduced dynamical regimes (see Section 8).  In the PDF documentation this threshold is referred to only in its reciprocal, using the symbol gamma.

The threshold is used to classify a variable within a particular epoch as either possessing an "order one" speed of response to perturbations, compared to the primary focused variable of an analysis (e.g. the membrane potential of a neuron), or a `fast` or `slow` response compared to the focused variable (order gamma or order 1/gamma, respectively).  These classifications are the basis for further reducing the number of effective dynamical variables in local models.

 

Note that we do not consider variables whose time-derivatives are very small or very large to be necessarily slow or fast.  For instance, along or near an orbit such as a limit cycle, variables with very fast responses can be slaved to a slower variable, so that their time-derivatives remain small, but their response to perturbations would be very fast.  DSSRT uses an additional classification -- `near-constant` -- to describe any dynamic variable that happens to not significantly vary over an epoch (relative to the variable bounds declared in the .cfg file).

 

 

7. Understanding the .cfg specification format

 

7.1 Overview

 

·             &nbs p;       A .cfg file is a list of specifications for an ODE system and its display in the DSSRT window.

·             &nbs p;       Comment lines begin with `#`.  All other non-blank lines will be treated as specification commands.

·             &nbs p;       Command lines must not end with comments or additional text.

·             &nbs p;       Variables not generated by explicit DEs must be declared internal variables (e.g. `u2` variable in `E_demo`), e.g. use for 'external' time-varying inputs in non-autonomous systems.  Strictly, such formal 'internal' variables are still 'observables' (see definitions in Section 6.1).

·             &nbs p;       The order of declared variables must be the same as the order of the variable data file columns.

·             &nbs p;       .dat data files may contain more columns than the total number of system variables, for instance if the software that computed them also stores auxiliary data in those columns (e.g. in XPP).  This is true for some of the data files in the included demos.  Any additional columns are simply ignored by DSSRT.

·             &nbs p;       A DSSRT diagram state is defined to be the states of the set of nodes and the links between them.

·             &nbs p;       Non-regular state maps (see Section 7.4 below) must be used to include activity of internal variables that are not given explicit nodes (e.g. the action-potential gating variables in `EI_demo`).

·             &nbs p;       The DSSRT diagram drawing area has normalized co-ordinates [0,1] x [0,1], for use in the diagram specifications.

 

 

7.2 Correct order of system specifications

 

(1) Variable declarations

            VARSEXT

            VARSINT

The order of the variables listed { VARSEXT, VARSINT } must correspond to the order of numerically integrated variable data in a .dat file loaded into DSSRT.

 

(2) Network connectivity, i.e. variable inter-dependence

            INPUTS

One for each variable, in any order.  A differential equation must exist for any observable having more than one external variable input.

 

(3) Information about variable bounds

            BOUNDS

            UNITBOUNDS

There can be a number of these specifications, and neither is required, provided all variables either have their own BOUNDS spec or are mentioned in the UNITBOUNDS list.  Due to the common occurrence of unit-interval variables, this list obviates the need for individual BOUNDS specs for each of these variables.

 

(4) Differential equation set-up

(4.1)      DEQNS

One of these specifications, listing which variables have an associated differential equation being declared to DSSRT.

(4.2)      GAM1TERM

            GAM2TERM

At least one of either, for every declared differential equation.  These declare all the terms in the right-hand side of the differential equations, according to the classification of Gamma1 and Gamma2 term-types (see full documentation for details).  `Tau-reciprocal` and asymptotic `target` values for each term can either be declared DE parameters (see DEPAR specification), or .m files in the set-up directory that specify the appropriate functions to generate these values (see Section 7.6).  Numerical values for these two parameters cannot be used directly in a GAM?TERM declaration, and no arithmetic combinations of parameters are possible in the declaration (ie. separate DEPARs for every coefficient are needed even if they are algebraic combinations of the fundamental system parameters!).  Inputs in the Gamma1 set are allowed no more than one cross-multiplying internal variable (e.g. the 'h' inactivation variable in the sodium activation term m^3 * h, for the Hodgkin-Huxley equations).  This should be enough for most modelling purposes.  Currently, Gamma2 terms cannot include cross-multiplying variables, but this will be supported in a later version.

(4.3)      DEPARS

As many of these as there are parameters used in the GAM1TERM and GAM2TERM specs

 

(5) Functional network diagram set-up

(5.1)      NODE

            NODELABEL

            NODEADDSTATES

            NODESTATEMAP_ACTS

            NODESTATEMAP_POTS

One set of these specifications, in this order, for each external variable.  If no additional states for the node are needed (the usual case) then these specifications must still be present, using `NODEADDSTATES <Var> 0` and then `NODESTATEMAP_ACTS <Var>', `NODESTATEMAP_POTS <Var>', where <Var> is the relevant observable variable name.  See `EI_demo` for a typical use of the additional node states, in order to avoid the clutter of too many nodes for the action-potential generating variables in the Hodgkin-Huxley equations.  Note that a node label does not have to correspond to the variable name used internally to DSSRT.

(5.2)      LINK

One of these for each variable interdependency that involves variables with their own diagram NODEs (i.e. for every coupling between observable variables described by terms in a differential equation).  This specification displays an arrow between the connected variables in the main graphical display, which will be colour-coded (or absent) depending on the 'functional network' state (see Section 7.5).  The display of unwanted arrows can be suppressed by specifying `1` in the optional 7th argument.  These specifications can be in any order.  (A future version will alleviate the need to specify explicit arrow co-ordinates.)

(5.3)      VBAR

One of these for every declared variable, in any order.  The final argument is the switch for allowing log-scaling: 0 = `off`; 1,2 = `on` ... but 1 means the low-value end is magnified, and 2 means the high-value end is magnified (e.g. for variables that usually rest near their maximum values rather than their minimum, such as the sodium inactivation variable `h` in the Hodgkin-Huxley equations).

 

 

7.3 Full syntax for DSSRT specifications

 

VARSINT <varname> [ <varname> ... ]

VARSEXT <varname> [ <varname> ...