Model specification format

Six different xml files must be prepared to run a NeuroRD simulation. The master or model file contains parameters for the simulation run and names of the other five files providing data for the simulation. The other five files include morphology, reaction scheme, initial conditions, output definition and stimulation. The reactions scheme file contains names and characteristics of all molecules and reactions in the simulation, including diffusion constants. The morphology file provides details on compartments, such as dendrites, soma and spines. The initial conditions file specifies the beginning amounts and locations of molecules. The stimulation file provides details of signal inputs during the simulation. The output file specifies which molecules are provided as the results. These results, in the form of the number of particles of each type in each element at each timestep, are stored to files for processing once the simulation has completed.

The independent definition of morphology, reactions and stimulation in separate files makes NeuroRD ideal to investigate the interaction of these characteristics in models. For example, the reaction scheme is independent of morphology to facilitate the investigation of the role of morphology, and to allow different morphologies to be investigated with different signaling pathways and vice versa. It is important to note that not all file are independent of each other: 1. Reactions and morphology _are_ independent of each other. These files define the species involved in the simulations and the morphological compartments in which they can be found. Most of the other files depend on these "core model" files. If these change, it is likely that all subsequent files must change. However, the subsequent files can be changed yet retain the same reaction and morphology files (this is typically how simulation experiments are performed). 2. Initial conditions depends on the reaction file (where molecule species are defined) and and will depend on morphology when some species are confined to subregions of the model. Changes to these files allow exploration of the role of molecule concentration and localization. 3. Similarly, the stimulation file depends on the reaction file (you can only inject a species already defined) and morphology (you can only inject at a defined regions or labels). 4. The output file clearly depends on species defined in the reaction file. Output from specific regions depends also on morphology file.

In all files, lines beginning "<!--" are XML comments and can be left out of real model files. Comments are ended with “-->" and can span multiple lines.

Reactions Scheme File

The reaction file has two parts. The first part is a listing of all molecule species. There is one program statement for each molecule in which each species is given a name, an id, a diffusion constant, and units for the diffusion constant. Species names can have spaces in them, but the id cannot. The id usually does not have to be specified, the name with spaces and slashes replaced with underscores will be used by default. The format is:

.. code-block:: xml

     <Specie name="name" id="id" kdiff="constant" kdiffunit="µm²/s" />

Two examples follow:

1. example for two diffusing molecule specification:

.. code-block:: xml

     <Specie name="IP3"   kdiff="100"   kdiffunit = "µm²/s"/>

2. examples for non-diffusing molecule specification:

.. code-block:: xml

     <Specie name="PIP2"  kdiff="0" />

.. code-block:: xml

    <Reaction name="name" id="id">
        <Reactant specieID="id"/>
        <Product specieID="id"/>
        <forwardRate> r₁ </forwardRate>
        <reverseRate> r₂ </reverseRate>
        <Q10> value </Q10>
    </Reaction>

.. code-block:: xml

    <Reaction>
        <Reactant specieID="PLCaG"                         />
        <Reactant specieID="PIP2"                       />
        <Product  specieID="PLCPIP2"                      />
        <forwardRate> 0.83e-06              </forwardRate>
        <reverseRate> 0.1e-03              </reverseRate>
    </Reaction>

    <Reaction>
        <Reactant specieID="PLCPIP2"                      />
        <Product  specieID="PLCaG"                         />
        <Product  specieID="IP3"                      />
        <forwardRate> 0.58e-03               </forwardRate>
    </Reaction>

.. code-block:: xml

   <Reactant specieID="cAMP" n="2"/>

.. code-block:: xml

   <Reaction name="decay">
      <Reactant specieID="A" power="2" />
      <forwardRate>0.83e-06</forwardRate>
      <Q10>0.2</Q10>
   </Reaction>

Morphology File

The morphology file specifies the morphology of the segments (which may be subdivided into smaller subvolumes for simulation). There are three program statements for each segment. In the first, each segment is given a unique identifier and may be give a region name, which does not have to be unique. Regions are used to group segments with the same initial conditions. The other two statements define the beginning coordinates, ending coordinates, radius, and an optional label. The label can be used as a site at which molecules are injected into the system, but labels cannot be used with "start on" specifications.

In general, segments are specified with a starting x,y,z and radius, and an ending x,y,z and radius:

.. code-block:: xml

     <Segment  id="seg1" region="dendrite">
         <start x="1.0"     y="1.0"   z="0.0"  r="0.75"  label="pointA"  />
         <end   x="1.0"     y="3.0"   z="0.0"  r="0.75"  />
     </Segment>

If a second segment is to be connected to the first it should start relative to the first compartment:

.. code-block:: xml

     <Segment  id="seg2" region="soma">
         <start on="seg1"  at="end"/>
         <end   x="1.0"   y="4.0" z="0.0" r="1.5"/>
     </Segment>

When two segments are connected, they need not have the same radius at the point of connection. In this case, a radius value must be specified when defining a new segment, as in the following example:

.. code-block:: xml

     <Segment id="seg2" region="branch1">
         <start on="seg1" at="end" r="0.3" />
         <end   x="4.0"  y="1.5" z="0.0" r="0.1" />
     </Segment>

If two segments are specified without specifying "start on", then thos two segments will not be connected, and molecules will not diffuse between them.

Branching segments are allowed. Branches are made by creating two segments beginning at the same site but terminating at different points. Where two daughter segments emerge from the same parent, they are offset so as not to overlap. For two dimensional models there is no freedom in the direction of offset. The grid elements are always rectangular, so daughter elements are connected as though they continued in line with the parent even if the morphology file shows them at an angle.

.. code-block:: xml

    <SpineType id="spineA">
        <Section width="0.2" at="0.0" />
        <Section width="0.2" at="0.6" regionClass="neck" />
        <Section width="0.4" at="1.0" />
        <Section width="0.4" at="1.2" />
        <Section width="0.1" at="1.3" label="pointA" />
    </SpineType>

.. code-block:: xml

    <SpineAllocation id="sa1" spineType="spineA" region="region1" lengthDensity="1.1" />

Initial conditions File

The initial conditions file specifies the initial  concentrations or densities of molecules. The file must contain one general concentration set, which applies to everything unless overridden. Each statement names the species and provides a value for its concentration, entered in nanoMoles per litre. An example follows:

.. code-block:: xml

     <ConcentrationSet>
         <NanoMolarity specieID="glu"     value="0"  />
         <NanoMolarity specieID="calcium" value="50"  />
     </ConcentrationSet>

In addition, further sets can be defined with a "region" attributed added after the "ConcentrationSet" . This should correspond to a specified region from the morphology file and indicates the parts of the structure to which the conditions apply. This only makes sense for non-diffusing molecules.

For membrane localized molecules, it is possible to specify initial conditions as a density (picomoles per square meter) which places these molecules only in the submembrane voxels of the morphology. The value attribute for a PicoSD element is the number of picomoles per square metre. For comparison with the volume concentrations, a surface density of 1 picomole/m^2, if spread over a layer 1 micron deep, gives a 1 nanoMolar solution. To average one particle per square micron, you need a PicoSD value of about 1.6. If a region is specified, then that initial condition applies only to that region. If no region is specified, then the initial condition applies to all submembrane voxels. The following example includes the optional region="dendrite" to show its use:

.. code-block:: xml

     <SurfaceDensitySet region="dendrite">
          <PicoSD  specieID="GaGTP"  value="003.729"    />
          <PicoSD  specieID="PLC"    value="2.521e+04"  />
     </SurfaceDensitySet>

Note that the initial conditions file should contain somewhere within it a statement for each species listed in the reactions file. See Purkic.xml for a complete initial condition file.

Stimulation File

The stimulation file specifies the time and location of injection of molecules (which is optional) during a simulation. For example, calcium influx from extracellular space might occur or glutamate might be released from a neighboring terminal. Each program statement must specify the molecule injected and its site of injection. The injectionSite needs to be a labeld point, either a spine (below) or labeled segment. The injection into a segment occurs in the middle mesh element of the labeled end of the segment. There are two ways to specify the timing and rate of injection: in the first a pattern is interpreted by NeuroRD, in the second a simple array of values specifies the pattern.

In the first version, additional statements provide onset (in ms), duration (in ms) and rate (particles/ms). Optional statements can be used to specify a train of input by providing two more parameters: period, and end.

Multiple trains are possible with two more parameters — intertrain interval and number of trains:

.. code-block:: xml

     <InjectionStim specieID="glu"  injectionSite="pointA">
         <onset>              100             </onset>
         <duration>           10              </duration>
         <rate>                50e3            </rate>
         <period>             400              </period>
         <end>                200              </end>
         <interTrainInterval>  10000           </interTrainInterval>
         <numTrains>           2               </numTrains>
     </InjectionStim>

Note that the intertrain interval specifies the interval between repetition of the entire train (excluding the onset time).

In the second version, the rates are given as an array of times and points. The injection rate is assumed to be constant between two time points. Thus a last time point with a rate of 0 should be used to terminate an injection that is bounded in time.

.. code-block:: xml

     <InjectionStim specieID="glu"  injectionSite="pointA">
          <rates> 100 30
                  120  0
                  300 40
                  320  0
                  400 50
                  420  0
          </rates>
     </InjectionStim>

For longer tables of injection values it might be convenient to store the table in an external file. This can be achieved using the standard XML include mechanism:

.. code-block:: xml

     <InjectionStim specieID="glu"  injectionSite="pointA">
          <rates>
               <xi:include href="rates.txt" parse="text" />
          </rates>
     </InjectionStim>

and the rates.txt file:

100 30
120  0
300 40
320  0
400 50
420  0

Since particles can only be injected, and not withdrawn, to produce transient elevations in concentration it may be necessary to inject a "binding partner" and add a reaction between the injected particle and the binding partner to lower the concentration of the unbound injected molecule.

Injection to spines is also possible by specifying a spine as follows:

.. code-block:: xml

    <InjectionStim specieID="a" injectionSite="sa1[3,4,5].pointA">
         <onset>5.0</onset>
         <duration>10.0</duration>
         <rate>210</rate>
    </InjectionStim>

The square brackets can contain:

• a number i - matches just the specified point in the i-th spine
• a comma-separated list of numbers - matches the points on those spines
• a range specified with a colon, such as [1:4]. If the lower or upper limit is missing it is taken to be 0 or the population size respectively.
• an asterisk, [*] to match the whole population.

The rate specified with <rate> is understood to be the total injection rate. In this example, we are injecting into three spines, so the rate will be 70 particles/ms for each of the spines, for a total of 210 particles/ms.

See Purksmlstim.xml, Purkstim.xml, and Purkdifstim.xml for examples.

Output scheme File

The output file specifies the file to which output is written and which molecules from which compartments are output at which dt. Multiple output blocks are allowed, for example if you want some molecules output more frequently than others, or from different regions.

Every OutputSet block must have in its definition one (and only) instance of: * filename Where filename is a string that specifies a sufix appended to the main output file name Additionally, every OutputSet block might have one (and only) instance of: * region or; * outputInterval A separate statement in each outputset block is used to indicate molecules to be included in the output file. For example:

.. code-block:: xml

     <OutputSet filename = "dt01sml"  region="dendrite" outputInterval="1.0">
         <OutputSpecie                 name="glu"  />
     </OutputSet>

If "region" is omitted then the concentrations (or number of particles) for all subvolumes in the system are saved. If "dt" is ommited then the concentrations (or number of particles) will be saved at every timestep (probably not a good idea).

Other examples are in Purkdifio.xml, and Purkio.xml.

Model File

The Model file is the “master file”, and serves to identify the files for other components needed to define the model as well as the type of calculation to be performed, discretization options, simulation seed(s) and various control parameters.

The first part of this files specifies the other files

.. code-block:: xml

     <xi:include href="Purkreactions.xml" />
     <xi:include href="Purkmorph2.xml" />
     <xi:include href="Purksmlstim.xml" />
     <xi:include href="Purkic.xml" />
     <xi:include href="Purkdifio.xml" />

The reactions, morphology, and initial conditions represent the model. The stimulation and output files are not part of the model, but are part of the simulation experiment.

The remainder of the file specifies various run parameters

A geometry statement is used to specify how the morphology is interpreted. 2D implies that there are mutliple voxels in x and y directions, but only a single layer of voxels in the z dimension. Thus, there is a 3 dimensional volume, but diffusion occurs in 2 dimensions only. For 2D, you also specify the depth of the voxel. This parameter is ignored for 3D:

.. code-block:: xml

    <geometry>2D</geometry>
    <depth2D>0.5</depth2D>

A runtime statement is used to specify run time in milliseconds:

.. code-block:: xml

    <runtime>                    2000                 </runtime>

A required simulationSeed statement specifies the seed for the random number generator. If spines are placed randomly, a separate random number generator is used, and a spineSeed statement must be added:

.. code-block:: xml

    <spineSeed>123</spineSeed>
    <simulationSeed>123</simulationSeed>

The discretization statement indicates how to subdivide the segments, i.e. the size of voxels or subvolumes used in running simulations. Smaller sizes require more calculations and result in a longer run time. Three different statements can be used within a discretization statement. The defaultMaxElementSide specifies the largest size (in microns) for each side of the subvolume in each segment. This is the default, and can be overriden by <maxElementSide> for specific regions, which specifies a region and the size of its voxels. The value supplied will be used only for that region. The <spineDeltaX> element specifies the size of subvolumes in spines. spineDeltaX defaults to defaultMaxElementSide. Spines have a one dimensional discretization. Example:

.. code-block:: xml

    <discretization>
        <!-- discretization for spines, microns -->
        <spineDeltaX>0.1</spineDeltaX>

        <!-- default largest size for elements in bulk volumes, microns -->
        <defaultMaxElementSide>0.2</defaultMaxElementSide>
        <maxElementSide region="dendrite">0.4</maxElementSide>
     </discretization>

The actual size of the elements depends on the total radius or length of the compartment, and also the constraint that there are an odd number of voxels across the radius of the structure. The maxElementSize is approximately the maximum that you will achieve. I.e., with a maxElemetSide of 0.4, you may generate subvolumes of size 0.33, depending on the size of the compartment. It calculates the length of compartment, divided by maxElementSide. Then, it determines the number of subvolumes along the length by rounding the results. Then, it divides length by number of compartments to yield the actual element size. Thus, you can end up with a value slightly larger than max element size.

The timestep statement specifies the time step, in milliseconds, used in fixed step calculations:

.. code-block:: xml

    <fixedStepDt>  0.01     </fixedStepDt>

The outputQuantity statement specifies whether quantity of molecules in the output is number of molecules (NUMBER) or concentration (CONCENTRATION):

.. code-block:: xml

    <outputQuantity>NUMBER</outputQuantity>

The outputInterval statement specifies the frequency (time interval) for writing out the complete state of the system. This file can be read by ccviz.

.. code-block:: xml

   <outputInterval>1</outputInterval>

The few other parameters in this file you can ignore for now.

See Purkmodel.xml and Purkdifmodel.xml for complete model files

If you run Purkmodel.xml, the output should match the XPPAUT output from filename.ode. Then, if you use purkmorphsml.xml, and purksmlstim.xml, you should get similar results, but the stochastic fluctuations are significant for some of the molecules. Lastly, to create a two compartment model to see the effect of diffusion, you can use Purkdifmodel.xml

Running a Simulation

To run a simulation from the command line the following command should be issued:

.. code-block:: bash

   java -jar NeuroRD.jar model.xml [output]

where NeuroRD.jar file contains the NeuroRD executable byte-code and model.xml is the model file (“master” file that contains the <SDRun> XML tag). The optional argument output specifies the base name of output files; their names will be created by suffix appending to the main output file name: output.out, output.h5, output.log, ... If the last argument is ommitted, the output name is the input file name without the .xml suffix. If the output name is a path to an existing directory, the output will be created in this directory using the input file name without the .xml suffix.

A number of messages will be printed at execution time. The same set of messages is printed to standard output and to the log file (output.log). Creation of the log file may be disabled, see below.

1. The java executable is in the PATH for the current user in UNIX.
2. NeuroRD.jar and model.xml are located in the same directory from where the command is issued or the full paths for these files have to be included as well.
3. output is in the same directory from where the command is issued or the full path for the output file has to be included as well.

With the default options, output is written as a set of text files and three (or more) output files are generated. One is the model.out file, which contains every molecule in every subvolume at a time interval specified by output interval. A second is the mesh file, which lists four xyz coordinates, depth and volume of every mesh element in the system. This can be used to check the morphology, and to convert from molecule quantity to concentration. The third file (or files) are those specified in the IO.xml file.

The option -Dneurord.writers=h5 can be used to specify [HDF5] output. In that case one output file is created, plus the log file. The format of the HDF5 output file uses the .h5 extension.

It is also possible to use an HDF5 file as a source of the model and/or exact population. The .h5 file contains the serialized XML which was used to create the model. If a file with the .h5 extension is given as the first argument (instead of an .xml file), model description will be extracted from it.

The random seed is also extracted from the .h5 file, which means that by using the HDF5 file as input, we can rerun the exact same simulation. In case the source .h5 file contains multiple trials, -Dneurord.source_trial=N can be used to specify which trial's seed should be used. It is also possible to override the seed as usual with -Dneurord.sdrun.simulationSeed=X. If the current simulation includes multiple trials, the seed is ignored; this is the case always when multiple trials are simulated. In effect, if an .h5 file is used as model source, the seed may be the seed specified by the original model, the seed used for the original simulation, the seed specifed through Java properties, or the randomized value if multiple trials are executed.

When an .h5 file is used as model input, the initial population can be taken from the results of that simulation. This is achieved by specifying --ic-time=T, where T is a timestamp of some saved state (T >= 0). It is also possible to specify --ic-time=-1, which uses the last T found in the file. Option --ic-trial=N can be used to pick a specific trial. It is possible but slightly risky to load initial conditions from an h5 file that it also the output file: the original file is overwritten at the beginning of the simulation after loading the model and initial conditions.

Please note: the list of species in the file containing the input state must match exactly. This is checked by the program. The list of reactions doesn't have to match, but be aware that if the list of reactions (or their constants) is changed, if the input state was a steady state, it might not be anymore with the new set of reactions.

.. code-block:: shell

    java -Dneurord.sdrun.<element>=<value> ...

.. code-block:: shell

    java -Dneurord.sdrun.discretization.maxElementSide=3 ...

.. code-block:: shell

    java neurord.StochDiff -t 2000 ...

.. code-block:: shell

    java neurord.StochDiff --help

    usage: java neurord.StochDiff [option...] <model> [<output>]

    where the <model> is an XML specification of the model to run. The optional
    <output> specifies where the results should be stored (w/o extension). When not
    supplied, <output> defaults to <model> without the extension.

     -D <property=value>     set java property
     -i,--ic <arg>           output file to take the initial conditions from
        --ic-time <arg>      time to take the ICs from (default: none)
        --ic-trial <arg>     trial to take the seed from (default: 0)
        --log <arg>          log file name ("no" to disable)
     -s,--statistics <arg>   override statistics gathering
                             ({none|injections|by-channel|by-event}[:interval])
     -t,--runtime <arg>      override simulation time
     -v,--verbose            increase log level
        --version            print version string and exit

    Recognized properties (use as -D<property>=<value>):
    neurord.mmap_appender (default: true)       User MemoryMappedFileAppender
    neurord.writers (default: [h5])             Write output in those formats
    neurord.trials (default: 1)                 How many trials to run
    neurord.threads (default: 0)                How many threads to use (0 == #CPUs)
    neurord.log_events (default: false)         Log detailed information about every event
    neurord.curtail_leaps (default: false)      Do not allow leaps to extend past reporting time
    neurord.neq.update_times (default: true)    Update putative times using Gibson&Bruck
    neurord.neq.only_init (default: false)      Terminate after initialization is finished
    neurord.neq.check_updates (default: false)  Perform extensive checks of propensity changes
    neurord.neq.log_queue (default: false)      Log debug info about queue operations
    neurord.neq.log_propensity (default: false) Log debug info about propensity calculations
    neurord.neq.log_reposition (default: false) Log debug info about movements in the queue
    neurord.neq.log_debug_start (default: NaN)  Turn on debugging at this time in simulation
    neurord.neq.log_start_events (default: 99)  Print information about this many events at startup
    neurord.compression (default: 1)            Compression level in HDF5 output
    neurord.sdrun.* (default: none)             Override values in the XML tree

.. code-block:: shell

    java -Dlog.<logger-name>=<level> ...

.. code-block:: shell

    java -Dlog.neurord.numeric.grid.NextEventQueue=debug ...

.. code-block:: shell

    java --log=no ...

License

This software is licensed under the GNU Publice License, version 2 or any later, at your option.

[HDF5]: https://www.hdfgroup.org/HDF5/