Calcium Calculator, Victor Matveev, Version 7.7.2, GPL Jan 11, 2014

Physical units

General syntax
Math expression syntax
List of functions
Operator priorities
Arrays
Character string expressions
The "include" statement
Command-line parameters

Implemented Geometries
Spherically-symmetric geometry
Cylindrical symmetry
Conical symmetry

Geometry definitions:
volume Define an enclosure element
obstacle Define an obstacle
B.obstacle Special buffer obstacle
sphere Define a spherical enclosure
sobstacle Define a spherical obstacle
grid Grid definition
Non-uniform grid:
stretch.factor Define the stretch factor
stretch Stretch the grid along an axis
clamp Clamp a grid node to a location

Calcium and buffer definitions:
bgr Initial background concentration
D Diffusion coefficient
tortuosity Space-dependent diffusion coefficient
Boundary condition definition:
bc Specifying b.c's for each volume surface
bc all Specifying b.c's for entire surface
bc.define User-defined boundary condition
Pump Defining calcium pumps on the boundaries
Old backward-compatible syntax
Buffer definitions:
total Total concentration
kplus Calcium binding rate
kminus Calcium unbinding rate
KD Calcium affinity
Buffers with cooperative calcium binding

Ca.source Defining calcium channels
Uptake Defining volume-distributed calcium uptake rate
Ca[x,y,z] / Buf[x,y,z] Local concentration variable
Ca[] / Buf[] Average concentration variable
_ICa / _Charge Total current and total charge variables
Defining ODEs (calcium "sensors")
equilibrate Find a steady state

Time-dependent variables
table Data table from a file
min / max Tracking variable minimum / maximum
markov Discrete-state continuous-time Markov variables

Simulation specifications:
Run adaptive Use adaptive step size method
current / currents Define current through Ca2+ channels
Run Use fixed time-step scheme
Multiple Run definitions
Numerical accuracy control: Number_Of_Iterations_Per_PDE_Step
method = ODE Run ODEs only

Graphics / data output:
plot.method xmgr Pipe output to xmgr session
plot.method mute Ignore non-mute plot commands
plot.print Save plot data to disc files
plot types:
plot / plot mute Plot a time-dependent variable
plot 1D / plot 1D.mute Plot of concentration along an axis
plot 2D / plot 2D.mute Plot of the concentration in a 2D plane
plot binary Binary plot of entire concentration field
plot dump Save the entire concentration field to disc
Logarithmic plots
plot.string Putting text on the graph

More data input / output:
import Reading calcium/buffer configuration from a file
Import, Export Saving and retrieving the entire simulation state
print / append printing data to a file

Flow control statements (preprocessor commands)
for Running the simulation for multiple parameter values
track Data output when using "for" cycle
if / then / else / endif Conditional statement
exit Program termination command
continue Skip a for loop iteration

Appendix A: controlling verbosity level
Appendix B: debugging the script: "debug" and "echo" commands

Δ Physical Units

The following units are implied for numerical constants: Units for other quantities are derived from the above; for example, units for the diffusion coefficient are [l]2/[t] = µm2/ms
Of course, one is free to use a different choice of units; in that case one should ignore the unit symbols prompted by the program in various text messages.

Δ General syntax

  1. Definition language is Case-Sensitive.
  2. Comment character is "%":     % this is a comment
  3. Semicolon (";") separates statements that follow one another on the same line:
    grid 20 20 20 ; stretch.factor = 1.07
  4. To continue a definition on another line, put a line continuation string "..." at the end of the previous line:
    Ca.bc Noflux Noflux Noflux Noflux ...
          Pump Pump
  5. A parameter name can be used wherever a numerical value is expected:
    grid N N N ; N = 40
  6. A constant numerical expression may be used for defining any parameter, i.e. grid N M M ; N = 2 * M ; M = 10
  7. If a parameter or construct is defined multiple times within the script, only the first definition will be used.
  8. Parameter definitions may appear in arbitrary order throughout the script, unless the conditional if .. then .. else .. endif statement is used, in which case the order of appearance of definitions is arbitrary within each of the then .. else/endif and else .. endif blocks. The only exception are Run and current statements which specify the simulation runs themselves and are executed in the order of their appearance.
  9. Non-trivial program flow can only be arranged via combined use of the for loop statement, the conditional "if" .. "then" .. "else" .. "endif" statement, the "include" statement and the "print/append" statements (see the description of the "if" statement).
  10. Parameter/identifier names can include all kinds of special characters, if the names are enclosed in quotes: "Ca\s0\N\S2+\N" = 0.1 (this creates an xmgr-friendly name with a subscript and a superscript).

Δ Mathematical expression syntax

  1. A space between two numerical values implies multiplication:
    x = 20 exp( y ) is equivalent to x = 20 * exp( y )
  2. The following functions can be used in numerical expressions:
    sinh, cosh, tanh, sin, cos, tan, atan, exp, log, log10, sqr, sqrt, theta, sigma (sigma(x) = ( 1 + tanh(x) )/2 ), int (integer part of float), not (logical negation), rand(0) (random number between 0 and 1)
  3. Time variable is "t" (to be used only in time-dependent variable definitions - see below).
  4. Operator execution priorities:
    Operators are executed in the same order as in the C programming language. Below all operators are listed according to their execution order, from highest to lowest:

    Unary boolean not operator:
    !
    Binary exponentiation operator:
    ^
    Unary sign operator:
    -
    Binary arithmetic operators:
    * / mod
    + -
    Binary logical operators:
    < <= >>=
    == (equal) != (not equal)
    and
    or

  5. Boolean operations: logical operations evaluate to 1 if true, and to 0 if false. An argument of a logical operator is considered true if greater than zero, and false otherwise. Logical operators include the binary operators from the table above, and the not( ) function. The primary use of logical expressions is in the "if" conditional statement, but they may also be used in other expressions:

    F := 2 * (t > 1.0) + (t <= 1.0)

    - defines a function of time equal to 1.0 when simulation time is less than 1.0 ms, and equal to 2.0 otherwise.
  6. Constant arrays may be defined, using a space-separated list:

    Data = 1.2   0.9   0.8   1.1   1.3   1.5       % defines the array
    entry = Data{4}       % getting the 4th element (Data{4} = 1.1)

    The curly-bracket references can also be used to extract arguments of a compound definition, such as the grid definition. Example:

    grid 20 40 30
    gridString = "The grid size is " grid{1} " by " grid{2} " by " grid{3} " points"

    To obtain the number of elements in the array, use array{0}. For instance, in the above examples grid{0} returns 3, and Data{0} returns 6.

Δ String expressions

A string is a sequence of characters in quotes, either "..." or '...'; the closing quote must be the same as the opening quote:

string_variable = "@#$%^!"

To include quote characters in a string, mix-and-match quote types: doubleQuote = '"'
Any string label/name (for instance, a file name) can be defined as a concatenation of several substrings:

string = token1  token2  ...

where token1, token2, ... are either strings or numerical expressions. The latter will be evaluated and their values converted to strings. Example:

param = 4
fileName = "var."   2 * param
plot  mute  var  fileName

The time dependence of a user-defined variable "var" will be saved in file "var.8". This feature is particularly useful when running a for loop over a parameter value.

There is a side-effect of this free-form mixing of string and numerical expressions: if a string expression matches a variable name, it will be interpreted as a variable. In the above example, string defined as fileName = "var." "param" is equal to var.4, not var.param.

Δ The include statement

A script file can import another script file using command include file_name. The script file_name will be inserted at the point where the include command appears. One may also use a string variable as a file name, if it is defined earlier in the file, i.e.:

script = path "/script"
include script


Δ Command-line parameters

Command-line arguments passed to CalC are accessible within the script; the n-th argument is denoted $n (note that the 1st argument is the script file name itself). The total number of arguments appearing on the command line is stored into a constant named $$

Δ Implemented Geometries

Although the description of various definition commands that follow below assume a 3D model in cartesian coordinates (the default geometry), various lower-dimensional geometries have also been implemented. To specify the model geometry, use the command

geometry = geometry_name

where the geometry_name is a geometry label that assumes one of the string values listed in the table below:

Dimensionality Geometry label Axis labelsGeometry description
1D cartesian.1Dx 1D model: diffusion along a line
discr Diffusion inside a disc: rotationally symmetric geometry
sphericalr Fully spherically symmetric geometry
2D cartesian.2Dx y 2D model in cartesian coordinates
cylindricalr z 2D model in cylindrical coordinates
polarr phi 2D model in polar coordinates
conicalr theta Model in spherical coordinates, symmetric with respect to phi azimuth angle
3D cartesian.3D
<< default >>
x y z Full 3D model in cartesian coordinates
spherical.3Dr theta phi Full 3D model in spherical coordinates
cylindrical.3Dr phi z Full 3D model in cylindrical coordinates


When using any of the lower-dimensionality geometries listed above, the argument lists of many of the definition instructions described below should be shortened. Additionally, some commands (e.g., plot and axis stretch instructions) use axis labels as arguments, which assume different string values depending on the geometry, as listed in the above table. Below is an example comparing some of the statements for the full 3D geometry in cartesian coordinates (geometry = cartesian.3D), and the reduced rotationally-symmetric geometry in conical and spherical coordinates:

geometry = cartesian.3D geometry = conical geometry = spherical
volume  0  1  0  1  0  1 volume  0  1  0  1.57 volume  0  1
grid  30  30  30 grid  30  30 grid  30
stretch x  0  0.1 stretch r  0  0.1 stretch r  0  0.1
stretch y  0  0.1 stretch theta  0  0.1  
stretch x  0  0.1    
Ca.source  0  0  0 Ca.source  0  0 Ca.source  0
plot  Ca[0.2,0,0] plot  Ca[0.2,0] plot  Ca[0.2]
plot 1D  Ca  x  0  0 plot 1D  Ca  r  0 plot 1D  Ca
plot 2D  Ca  z  0.2 plot 2D  Ca  

Few special geometry cases (spherical, cylindrical and conical) are described in more detail below, as an example.

Δ Spherically-symmetric geometry

When the geometry is spherically symmetric, the 3D problem reduces to a 1D problem in spherical coordinates. Examples are: a source in the center of a sphere, or a current uniformly distributed over the surface of a spherical cell. In this case it is much more efficient to solve the resulting one-dimensional problem. To do so, use the command

geometry = spherical

and only use two parameters when defining the volume:

volume  R0  R1

where R0 and R1 are the inner and outer radii of the spherical shell that serves as the diffusion space in this case. Only two boundary conditions need to be defined, one for the inner and one for the outer surfaces of the shell.
Similarly, only two parameters specify the Ca current source:

Ca.source  R  dR

where "R" specifies the radius of the spherical surface over which the current is uniformly distributed (R=0 for a problem with a source in the center of a sphere), and "dR" specifies the width of the source current layer in radial direction.
All other commands also remain unchanged; wherever an axis name is needed, use "r", and omit the "y", "z" coordinates. For example, use stretch r rmin rmax to stretch the grid outside of the shell [rmin, rmax], use Ca[r] to access calcium concentration at distance r from the center, and plot 1D Ca to plot the calcium concentration as a function of distance from the center of the sphere.

Δ Cylindrical geometry

If the geometry of the problem is symmetric with respect to the azimuth angle phi in cylindrical coordinates, the 3D problem reduces to a 2D problem in coordinates r and z. To solve the corresponding 2D problem, include the command

geometry = cylindrical

and only use four parameters when defining the volume:

volume  R0  R1  Z0  Z1

where R0 and R1 are the inner and outer radii of the volume, and Z0 and Z1 are the limiting z-values. Only four boundary conditions need to be defined. Similarly, only four parameters specify the Ca current source:

Ca.source  R  Z  [  dR  dZ ]

where R, Theta specify the location of the source, and optional parameters dR, dZ specify the spread of the source.
All other commands are also analogous to their 3D counterparts. Wherever an axis name is needed, use "r" for the radial coordinate, "z" for the 2nd coordinate. and omit the "z" coordinate. For example, use stretch z zmin zmax to stretch the grid outside of the domain [zmin, zmax], use Ca[r,z] to access calcium concentration at location {r, z}, use plot 1D Ca z r0 to plot the calcium concentration as a function of z, at a distance r0 from the center of the sphere.

Δ Conical geometry

If the geometry of the problem is symmetric with respect to the azimuth angle phi in spherical coordinates, the 3D problem reduces to a 2D problem in coordinates r and theta (e.g., see Klingauf and Neher (1997) Biophys J 72: 674). To solve the corresponding 2D problem, include the command

geometry = conical

and only use four parameters when defining the volume:

volume  R0  R1  Theta0  Theta1

where R0 and R1 are the inner and outer radii of the volume, and Theta0 and Theta1 are the limiting angle values, in radian (0 < Theta0, Theta1 < pi). For Theta0=0, the volume is a cone.
Only four boundary conditions need to be defined. Similarly, only four parameters specify the Ca current source:

Ca.source  R  Theta  [  dR  dTheta ]

where R, Theta specify the location of the source, and optional parameters dR, dTheta specify the spread of the source.
All other commands are also analogous to their 3D counterparts. Wherever an axis name is needed, use "r" for the radial coordinate, "theta" for the angle coordinate, and omit the "z" coordinate. For example, use stretch theta thetamin thetamax to stretch the grid outside of the angle domain [thetamin, thetamax], use Ca[r,theta] to access calcium concentration at location {r, theta}, use plot 1D Ca theta r0 to plot the calcium concentration as a function of angle, at a distance r0 from the center of the sphere.

Δ Geometry definitions:

The simulation space volume can be composed of (is a union of) an arbitrary number of "boxes", each defined via keyword "volume", with an arbitrary number of obstacles (barriers to diffusion of calcium and buffers), each defined using keyword "obstacle". Calcium and buffer concentration fields are computed on a mesh of points, specified using keyword "grid", covering the entire diffusion space. The grid may be non-uniform.

volume  xmin  xmax  ymin  ymax  zmin  zmax

- defines a box spanning area {xmin < x < xmax, ymin < y < ymax, zmin < z < zmax}. All bounding coordinates should be in micrometers. If the diffusion space is composed of several boxes, several volume definitions should be included.

obstacle  xmin  xmax  ymin  ymax  zmin  zmax

- defines a diffusion obstacle spanning an area {xmin < x < xmax, ymin < y < ymax, zmin < z < zmax}
An arbitrary number of obstacle definitions may be used.

To define areas that are impermeable to buffer molecules, one can also define special buffer obstacles. This is a simple alternative to defining non-uniform diffusivity through the tortuosity definition. The format for a special obstacle is very similar to the definition of a usual obstacle; one only needs to add the buffer name prefix:

buffer_name.obstacle  xmin  xmax  ymin  ymax  zmin  zmax


Δ Spherical elements:

In CalC versions 6.7 and above and 7.7 and above, in the default Cartesian 3D geometry the diffusion volume can be defined as a superposition of spherical enclosures, defined using

sphere  xCenter  yCenter  zCenter  radius

Analogously, one can define spherical dffusion obstacles (e.g. vesicles):

sobstacle  xCenter  yCenter  zCenter  radius

See a simple example here.

Note: This feature is still under testing; one should not place distributed Ca2+ sources near a curved boundary, as any overlap will result in a reduced current. Future versions 6.8 and 7.8 will improve handling of distributed sources near curved boundaries.


Δ Grid specification:

CalC uses a simple cubic grid, which is a direct product of spatial grids along the individual axes. To define a grid of size nx by ny by nz, use

grid  nx  ny  nz

The number of points in each direction should be an odd integer; an even number will be increased by 1.

When using geometries with spatial dimensionality of 1 or 2, use respectively 1 or 2 integer arguments only, corresponding to the respective grid dimensions.

Δ Non-uniform grid:

For the same number of grid points, the accuracy of the computation can be improved by increasing the density of grid point in the vicinity of calcium channel(s), where the concentration gradients are the greatest. If the channel array is enclosed within a region {xmin < x < xmax, ymin < y < ymax, zmin < z < zmax}, one should use:

stretch  x  xmin  xmax
stretch  y  ymin  ymax
stretch  z  zmin  zmax

- each successive grid interval further away from the region {xmin < x < xmax, ymin < y < ymax, zmin < z < zmax} will be multiplied by a factor specified by stretch.factor, which should be slightly greater than 1.0: The grid may also be stretched along a single direction (x, y, or z) only.

stretch.factor = value

- sets the grid non-uniformity factor used by stretch. Recommended values for the factor value are between 1.05 and 1.2. Using larger values, while making the results converge faster with increasing grid dimensions (i.e. spatial resolution), will cause the results to converge to a slightly wrong answer (stretching the grid introduces an additional systematic computation error).

In cases when the grid is relatively coarse, it may be desirable to have a grid node pass through some special coordinate in space:

clamp  x  y  z

- forces the closest grid node to assume coordinates [x,y,z]. Since this will change the neighboring grid spacings, this will distort the grid, so it's better to do without clamping, if possible.

The stretch command syntax is the same for spherical, disc, cylindrical and conical geometries, except that one may use, respectively, the direction labels "r", "r" and "z", or "r" and "theta", i.e.

stretch   r  rMin  rMax
stretch  theta  angleMin  angleMax


Δ Calcium and buffer definitions:

Definition of each Ca2+ property starts with a label Ca followed by a period, then followed by the property name:

Ca.bgr = expression

- defines the initial background calcium concentration (in units of µM)

Ca.D = expression

- defines the value of the calcium diffusion coefficient (in units of µm2/ms)

Buffer properties definitions have the same syntax. For example, definition of a diffusion coefficient of a buffer named Bf is Bf.D = expression

In order to define a diffusion coefficient which is non-uniform in space, use the command

Ca.tortuosity = f(x,y,z)

The diffusion coefficient of Ca2+ will be a product of a space-independent term defined by Ca.D and the space-dependent term defined by a user-specified function f(x,y,z) Example:

Ca.tortuosity = 1 - 0.5 exp( - 20 ( x*x + y*y + z*z ) )

- in this case the Ca2+ diffusion will be reduced in the vicinity of point x = y = z = 0.
Diffusion coefficient of a buffer can be modified the same way, for instance:

buffer Fura
Fura.tortuosity = 1 - 0.5 sigma( 50 (z - 0.2) )
Fura.D = 0.118
............

To define the same tortuosity function for Ca2+ and all the buffers, simply use

tortuosity = f(x,y,z)     or      all.tortuosity = f(x,y,z)

Note that the diffusion coefficient should have (approximately) zero derivative at the boundary in the direction normal to the boundary.

Non-uniform diffusion coefficient is implemented for all cartesian and curvilinear geometries (1D, 2D and 3D).

Δ Boundary condition specification:

By default, both calcium and buffers satisfy zero flux (reflective) boundary conditions at all surfaces of all geometry elements. However, boundary conditions may be specified by the user for each volume and each obstacle (including special buffer obstacles), by providing the following statement for each geometry element:

Ca.bc   bcXmin  bcXmax  bcYmin  bcYmax  bcZmin  bcZmax

This command specifies boundary conditions for each of the six surfaces of a volume element. Each of the labels bcXmin ... bcZmax is either defined by the user via a bc.define instruction described below, or is one of the two predefined types, Noflux or Dirichlet (Bgr), corresponding respectively to the zero flux (reflective) boundary condition, or to the boundary value clamped to background concentration. In 2D geometry, only 4 boundary condition types should be specified, and only 2 are needed for 1D geometry. If a concentration field satisfies the same boundary condition "bcName" at all boundaries of a corresponding volume element, use

Ca.bc  all  bcName

If several volumes and/or obstacles are specified, the boundary condition definitions appearing earlier in the script will correspond to volume definitions (in the order those were defined), and the following boundary condition definitions will correspond to the obstacles, including any special buffer obstacles.


Δ User-defined boundary condition:

The user definition of a boundary condition type is over-loaded for convenience, allowing the following alternative forms:

bc.define  bcName  A  B  C  n  K 

bc.define  bcName  A  B 

bc.define  bcName  R

where bcName is a user-chosen label, and parameters A ... R correspond to the parameters in the respective boundary conditions:
Boundary condition flux formula








Here [Ca] is the Ca2+ concentration at the boundary, and [Ca]0 is the background concentration. In the first two cases Flux is the flux per unit area at the boundary, with default units (in 3D) of µM µm / ms = 10-6 mol / (m2 s), and therefore the units of B are μm/ms = 10-3 m/sec, and units of C are the same as units of flux. Default units of pump affinity K are μM. Non-linearity power n can be any real number greater than 1.

For instance, to simulate a linear saturable calcium pump on the bottom and top ("zmin" and "zmax") surfaces of a volume, with maximal pump rate of 0.01 µM µm / ms, and dissociation rate of K = 0.2 µM, use the following command:

bc.define Pump   1    0   0.01   1    0.2    
Ca.bc Noflux Noflux Noflux Noflux Pump Pump

The Noflux (zero flux) and Dirichlet / Bgr boundary condition types are defined internally using

bc.define  Noflux  1  0 
bc.define  Dirichlet  0 
bc.define  Bgr  0 

Note that if the calcium or a buffer field is to assume a fixed concentration value at the boundary which is not equal to its resting background concentration, user has to provide his own redefinition for the Dirichlet boundary condition (using a different type label/name, since the bc's already defined may not be redefined); for instance:

bc.define  myDirichlet   0.1
Ca.bc  all  myDirichlet   % the Ca concentration value is the background concentration plus 0.1 µM on all boundaries


Finally, note that buffers are assumed to satisfy either no-flux or Dirichlet boundary conditions only. If non-zero flux boundary conditions are desired for a buffer, two additional parameters are required, the background concentration and mobility of the buffer:

bc.define  bcName  A  B  C  n  K  B_bgr   DB


Δ Backward-compatible boundary condition definition:

For compatibility with versions of CalC preceding 7.6, the following alternative definition may be used:

bc.define  bcName  A  B  C  P

corresponding to boundary condition
Boundary condition flux formula




where d[Ca]/dn is the derivative normal to the boundary facing out. The value of P is zero by default. For the pump example above, with maximal pump rate of Vm=0.01 µM µm / ms, and dissociation rate of K = 0.2 µM, use A = 1, B = −Vm / (K DCa) = -0.227 (µm-1), where DCa=0.220 µm2/ms is the calcium diffusion coefficient, and set P = 1 / K = 5 µM-1:

bc.define Pump  1.0  -0.227  0.0   5.0   
Ca.bc Noflux Noflux Noflux Noflux Pump Pump


Δ Buffer definitions:

An arbitrary number of buffers with one-to-one Ca2+ binding stoichiometry can be included in the simulation, each defined using

buffer  buffer_name

Values specified for the calcium concentration field should also be specified for each of the buffers. However, instead of defining the background concentration, it is necessary to define its total concentration, using

buffer_name.total = constant_expression

For a mobile buffer, this total concentration is assumed to be homogeneous throughout the diffusion space. Example:

buffer  myBuffer
myBuffer.D = 0.2
myBuffer.total = 2000
myBuffer.bc  all  Noflux         % (default setting)

Calcium binding kinetics should be specified for each of the buffers, by defining two of the following three parameters:

myBuffer.kplus = expression         % defines calcium binding rate, in 1/(ms uM)
myBuffer.kminus = expression       % defines calcium unbinding rate, in 1/ms
myBuffer.KD = expression         % defines calcium affinity of buffer (in uM): KD = kminus/kplus

Before the simulation is initiated, the program automatically calculates the initial buffer concentration from its kinetic parameters and total concentration value, assuming an equilibrium with the defined background [Ca2+] (unless the concentration field is initialized from a disk file using the import instruction).


A fixed buffer can be localized to a given segment of the simulation volume, by using a spatially-dependent function in the total concentration definition, for instance:

buffer  myBuffer
myBuffer.D = 0        % buffer is fixed
myBuffer.total = 2000 * (x*x + y*y + z*z <= 0.01)    % buffer is confined to a spherical region around the origin, with a radius of 0.1 µM

For non-cartesian geometries, remember to use appropriate space variables (axis labels) in the above expression for the total concentration (see table)
Further, one can define a fixed membrane-bound buffer, localized to a layer of a given depth over the entire boundary surface, using instruction

buffer  buffer_name   membrane [ depth ]

where depth is an optional numerical constant expression determining the thickness of the buffer layer surrounding the boundary.
If the depth is zero (which is its default value), the buffer will be localized to one grid layer around the membrane, and the total concentration should be specified in units of molecules per µm2. [ Note that in this case the buffer concentration will be higher in the corners of an enclosure, since the corner voxels will include contributions from the buffer bound to all surfaces meeting in that corner. ]

Δ Buffers with cooperative Ca2+ binding

To define a buffer with 1-to-2 calcium binding stoichiometry named "B", use keyword buffer cooperative:

buffer cooperative B

This will create three different concentration fields, with names B, Ca.B and Ca2.B, corresponding to the calcium binding states of the reaction

B + 2 Ca2+ CaB + Ca2+ Ca2B

The rates of the two reactions in the above diagram should be declared as properties of buffers Ca and Ca.B, respectively, even though the first unbinding is a property of Ca.B, not B: this simplifies the syntax. For example, the statements below define the properties of calmodulin's N-lobe as reported recently by Faas et al. (2011):

buffer cooperative N
N.kplus = 2 * 0.77     % first reaction rate is 0.77 uM-1ms-1
N.KD = 193 / 2     % affinity of first reaction is 193 uM
Ca.N.kplus = 32     % second reaction rate is 32 uM-1ms-1
Ca.N.KD = 2 * 0.788     % affinity of second reaction is 788 nM

Note that factors of 2 in above definitions are added to simulate equal probability of binding of the two Ca2+ binding sites of the N-lobe of apo-calmodulin, and to simulate equal probability of unbinding of each of the two Ca2+ ions from fully calcium-saturated N-lobe.

To simulate Ca2+ diffusion in the presence of calmodulin with both lobes intact, a second cooperative buffer of equal concentration has to be defined, with reaction rates corresponding to the C-lobe of calmodulin.

Note that diffusion coefficient of cooperative buffers are allowed to change upon Ca2+ biding, so three diffusivities have to be defined for a cooperative buffer named "B": namely, all three parameters B.D, Ca.B.D, and Ca2.B.D have to be assigned.

If the total concentration is defined for the unbound buffer only (e.g. N.total = 100), then the fully and partially bound states of the buffer are assumed to be in equilibrium with the background Ca2+. Alternatively, total concentrations of all buffer states can be given, in which case the simulation will start from the defined non-equilibrium initial condition, e.g.:

buffer cooperative B
B.total = 100
Ca.B.total = 50
Ca2.B.total = 20

These definitions may or may not represent an equilibrium concentration state, depending on the values of the binding rates and background Ca2+, defined by Ca.bgr.

Finally, note that the total concentration of any immobile buffer state can be declared to be a space-dependent function, as in the case of a non-cooperative buffer.

All CalC instructions straightforwadly generalize to the different states of the cooperative buffer: for example, to plot the concentration of the doubly-bound (saturated) state of cooperative buffer named "B" at a point with coordinates (a; b; c), use the standard syntax plot Ca2.B[a,b,c].

Δ Calcium channels:

An arbitrary number of calcium channels can be defined. For each, channel location and width of the spatial distribution of the current (which is of gaussian shape) need to be defined, using

Ca.source  x  y  z  [  dx  dy  dz ]

where the optional parameters dx, dy and dz specify the width of the source (square brackets indicate that parameters are optional; do not include brackets in the argument list!). The width parameters dx, dy, dz equal the Gaussian parameter σ multiplied by √2.

The default value for the width parameters is zero, which corresponds to a point-like source. However, the point-like nature of calcium channels introduces a spatial singularity into the problem, leading to questionable convergence of the numerical method with increasing spatial resolution. Therefore, unless the calcium concentration in the microdomain of a single channel is studied, one should set the width parameters so that the source covers at least 1-2 grid intervals. This will improve the numerical accuracy and the convergence properties of the method. If only the "dx" parameter is defined, the source spread in all three directions is assumed to be the same.

The shape of the spatial current distribution function can be changed from a gaussian function to a step function, using instruction

current.shape    square

In this case the current strength is assumed to be spread uniformly over a rectangular area defined by [x - dx, x + dx] * [y - dy, y + dy] * [z - dz, z + dz].

Channels can be located anywhere, either inside the simulation volume or on the boundary.

The Ca2+ current magnitude through each channel is defined using the current / currents commands.

Δ Volume-distributed Ca2+ uptake rate:

There are two ways to define intracellular Ca2+ uptake into internal stores such as ER or mitochondria. First, one can define diffusion obstacles with an appropriate obstacles boundary condition on all sides of the obstacle. Alternatively, one can simulate a volume-distributed ER syncytium with a simple linear calcium uptake, using a single rate constant uptake:

uptake =   rate

where rate is the rate of uptake in units of ms-1. Versions 6.7.1+ and 7.7.1+ allow the rate to be defined as a function of the spatial coordinates, implementing space-dependent calcium uptake.

This linear uptake corresponds to an additional term in the reaction-diffusion equation for Ca2+:

d[Ca]/dt = D ∇2[Ca] + R([Ca], Buffers) - uptake * ([Ca] - [Ca]0)

where [Ca]0 is the background calcium concentration, and R([Ca], Buffers) contains terms describing buffer-calcium binding/unbinding reactions.
Note that as of now CalC does not allow to define calcium release from internal stores.

Δ Local concentration variable:

Local concentration of Ca2+ at some location {x,y,z} is accessed using the following time-dependent construct: Ca[x,y,z].

Buffer concentration at a given location can be accessed the same way, via variable buffer_name[x,y,z]. If the coordinate {x,y,z} does not coincide with any grid node locations, linear nearest-neighbour interpolation is used.

Δ Average concentration:

To monitor the average concentration of Ca2+, use the built-in variable Ca[]; analogously, the variable tracking average concentration of unbound buffer named buffer_name is buffer_name[].

Δ Current and Charge variables

Built-in time dependent variable _ICa gives the current value of the total Ca2+ current entering the volume through all defined channels. Another predefined variable, _Charge, gives the total Ca2+ charge that entered since the simulation start (it is defined via _Charge = d _ICa / dt). The discrepancy between the value of the _Charge variable and the actual integral of Ca2+ over the simulation volume is output at the end of each simulation, and can be accessed through the built-in variable Charge.loss. In the absence of Ca2+ uptake, Dirichlet or pump boundary conditions (which would cause the charge to escape through the surface), this variable can be used as a simulation accuracy measure.

Δ Defining ODEs

To simulate the binding of calcium to putative neurosecretory/facilitation sensors, a corresponding system of differential equations (ODEs) has to be defined. A differential equation is defined using dvariable / dt = expression_for_derivative, for example

dR / dt = - koff R + kon Ca[0.0,0.0,0.1] ( 1 - R )

Alternatively, one may denote differentiation using the backquote "`" instead of the d/dt construct, i.e.:

R`= - koff R + kon Ca[0.0,0.0,0.1] ( 1 - R )

Note that this differs from the xpp syntax, where the symbol for defining an ODE is a normal quote. In the above, the construct Ca[x,y,z] denotes calcium concentration at point {x,y,z} (see above).

For each equation, initial boundary condition can be defined:

R(0) = 0.0

The default initial condition is zero.

Note that all initial conditions are ignored if an Import data input statement is included.

Δ Finding the steady state of the ODE system

In order to integrate (see Run definition) the declared system of ODEs starting from a steady state, rather than from the initial conditions, include the command

equilibrate  step_num  accuracy  initial_time_step

where step_num is the maximal number of equilibration steps, accuracy is the desired accuracy (for the zero derivative condition), and time_step is the initial equilibration time. Default values are step_num = 40000, accuracy = 1.0e-6 and initial_time_step = 1 ms.

Δ Time-dependent variables:

It may be necessary to introduce additional time-dependent variables that are functions of ODE variables; this is done using a command of the form "variable := expression", i.e.

dR / dt = - koff R + kon Ca[0.0,0.0,0.1] U        % same equation as before
U := 1 - R        % variable denoting the unbound receptor

Variables (and expressions for derivatives) may depend on time explicitly, through time variable "t"

There are two more user-defined time-dependent variable types, in addition to variables defined via "d /dt", ":=", "_ICa"/"_Charge", and "Ca[x,y,z]"/"Ca[]" (or "buffer_name[x,y,z]"/"buffer_name[]"). These are table and min (max) variables described below:

Δ Data table from a file:

A time-dependent variable can be defined using a two-column data file (satisfying the format "tn f(tn)", without a comma or other delimiter in between):

variable_name  table  file_name

Time points don't have to be regularly spaced; linear interpolation will be used to fill in for other time values.

This command may be useful in cases when one wants to change parameters of the calcium binding kinetics without having to re-calculate the time-course of calcium concentration. For instance, in the example script demo.par, one may save [Ca] evolution at the site x=y=z=30 nm, assigned to a variable named ca, using plot mute ca "ca.dat", and then run the ODE-only version of the script after including the command ca table "ca.dat". Since at this second step only the ODEs corresponding to the calcium binding scheme have to be re-calculated, one should also include the mode = ODE command.

Δ Tracking extremum of a variable:

To obtain the maximum of a variable var on an interval [T1,T2], use

[  time_var_name   ]  max_var_name  max  var  T1  T2

where max_var_name is the name of the variable tracking the maximum of var, and time_var_name (optional) specifies the name of the variable that will contain the time value at which the maximum of var is achieved.
Although one is only interested in the value of max_var_name once the time point T2 is reached, CalC treats max_var_name as a time-dependent variable, so usual plotting routines may be used to monitor max_var_name.

This object is particularly useful in modeling synaptic facilitation, where one is interested in the peak value of [Ca2+], or a peak of a variable simulating calcium binding site occupancy (such as "S" in example file demo.par).

The syntax is analogous for tracking the minimum of a variable:

[  time_var_name   ]  min_var_name  min  var  T1  T2


Δ Markov process variable:

A discrete-state continuous-time stochastic variable can be defined using

variableName  markov  stateNumber   initialState

where the integer parameter stateNumber specifies the number of markovian states, and the parameter initialState specifies the state of the variable at simulation start.

Only non-zero off-diagonal transition probabilities should be defined, using constant parameters or time-dependent variables with names of form variableName.i.j. The following example specifies a markov gating variable with two states:

gate markov  2  0
gate.0.1  =  0.2
gate.1.0  :=  0.01 * Ca[0,0,0]

In this example the transition from state 0 to state 1 is constant at 0.2 per millisecond, while the backward transition probability 1 -> 0 depends on the concentration of Ca2+ at the origin.
Note that special care should be taken if the Ca2+ current depends on the state of a markov variable, since the Ca2+ and buffer PDEs are integrated deterministically, assuming no stochasticity. In particular, only a fixed time-step PDE integration method should be used in this situation, with a small time step (see an example). Using markov variables in ODE-only mode is not recommended.

To change the seed of the random number generator, use the command

seed  =   longInteger

where longInteger is a 32-bit integer value that will be used to initialize the high-order 32 bits of the random number generator state. On Unix platforms, a standard 48-bit random number generator drand48() is used.

Δ Simulation specifications:

The simulation run itself is defined using instruction Run, and requires specifying the integration time and the value of the calcium channel current during the run. There are two main variation for the Run command, corresponding to the adaptive (variable time-step) and the non-adaptive (fixed time-step) methods.
In both cases, it is also possible to change a particular finite-difference scheme used by the integrator: see Appendix A.

Δ Adaptive time-step method

When using the adaptive time-step method, the size of the integration time step is constantly modified to maintain a given level of accuracy (but see below). The simulation definition syntax in this case is

Run  adaptive  Time  [   accuracy    dtMax   dt0   dtStretch  ODEaccuracy   ]

where Time is the total integration time (in ms). The values in brackets are optional and specify the parameters for the adaptive algorithm. The initial time-step is dt0, which is reduced if the error tolerance accuracy is not met. Each successive time-step is increased by a multiplicative factor dtStretch, slightly greater than one; the time-step cannot exceed the value of optional parameter dtMax. Accuracy checks are performed every n steps: if accuracy check is successful, n is increased by one, otherwise the time step is halved, and n is reset depending on the number of steps since last time step reduction. If the error is greater than five times the accuracy tolerance, the integrator takes a reverse step back to the previous successful accuracy check. Finally, the value ODEaccuracy specifies the error tolerance for integration of the system of ODEs (e.g., corresponding to the user-defined calcium binding receptor kinetics). Default values are: accuracy=1e-5, dt0=1e-3, dtMax=0.1, dtStretch=1.03, ODEaccuracy=1.0e-4. Values of the adaptive scheme parameters can also be set individually by assigning values to corresponding pre-defined identifiers (which are the properties of the adaptive method), i.e.: adaptive.dtStretch = 1.1; adaptive.accuracy = 1e-8; adaptive.dtMax = 0.4. Finally, there are two parameters that control the minimal and the maximal number of integrator steps between successive accuracy checks: adaptive.steps = 3; adaptive.maxSteps = 20. There is little need to change these default values under most circumstances.

Note: parameter adaptive.accuracy controls the accuracy of a single time-step of the finite difference method. The only rigorous way to check the solution accuracy is to examine its convergence using the non-adaptive intergator, by gradually decreasing the (non-adaptive) time-step size. That said, under most conditions the adaptive integrator provides a much more efficient way of obtaining a fairly accurate result, as compared to the non-adaptive method. If the solution is not sensitive to a change in adaptive.accuracy parameter, chances are it is pretty accurate. An additional indirect accuracy check is provided by the total charge conservation check echoed after each simulation (see here).

Each Run statement should be accompanied by its own current definition statement, which accepts one of the following two forms:

current  = expression
  or
currents  { blank-separated variable list}

where expression is a numerical expression, variable name or a constant specifying the calcium current per channel. The first form of the current statement is used when the same current flows through each of the defined Ca++ channels. The second form of the current statement is used when several channels are present, with distinct Ca++ current magnitudes. In this case the keyword currents should be followed by a list of variable names, with the number of variables in the list equal to the number of channels defined by the Ca.source commands.

Example 1:

Run  adaptive  1.0    0.001  1e-6 
current  = 0.2 pA

- Runs the adaptive scheme for 1 ms, with Ca++ current of 0.2 pA per channel, with the initial time step of 1 µsec, and error tolerance of 1.0e-6.

Example 2:

Run  adaptive  10
currents  I1  I2
I1 = 0.1 pA
I2 := t * 0.2 pA

- Runs the adaptive scheme for 10 ms; each of the two variables following the currents keyword specifies the current through each of the two Ca++ sources, defined elsewhere in the script. The Ca++ current through the first channel is 0.1 pA, while the current through the second channel is time-dependent (linearly increasing with time).


Δ Fixed time-step method

The syntax for the non-adaptive simulation run is:

Run  Time  dt  [  ODEaccuracy  ]

where "Time" is the total simulation time (in ms), "dt" is the simulation time-step, ODEaccuracy is an optional parameter setting the accuracy for the integration of ODEs (if any are defined).

As in the case of the adaptive simulation, each Run statement should be accompanied by a current definition statement, described above.

Example 1:

Run  1.0    0.001  1e-5 
current  :=   0.2  t^2  pA

- Runs the fixed time-step scheme for 1 ms, with a Ca++ current that is increasing with time quadratically, with a fixed time step of 1 µsec, and Runge-Kutta ODE integration error tolerance of 1.0e-5.

Example 2:

Run   3.0    0.02
currents   I_Ca1   I_Ca2

- Runs the fixed time-step scheme for 3 ms, with a fixed time step of 20 µsec, and the Ca++ currents of I_Ca1 and I_Ca2 through the two channels (two Ca.source commands should be present in the script), where I_Ca1 and I_Ca2 are either constants or time-dependent variables defined by the user.


Δ Multiple "Run" definitions

The simulation may be broken down into an arbitrary number of shorter simulations; in this case, there should be a current definition for each Run definition. Different algorithms may be used for successive runs.

Breaking up the simulation into multiple "runs" is desirable when the calcium channel currents change abruptly during the run, as when simulating channel opening and closing events. If an adaptive method is used, it is important that the channel current is continuous during each of the individual "runs". This will limit the error arising from large current discontinuity. Thus, it is desirable to use

Run  adaptive  1.0    % 1 ms-long channel opening
current  = 0.2 pA
Run adaptive 3.0    % 3 ms-long channel closed period
current  = 0

rather than

Run  adaptive  4.0    % 4 ms-long channel opening
current  = 0.2 pA theta( 3 - t )

although formally the above two definitions are equivalent.

Δ Numerical accuracy control: Number_Of_Iterations_Per_PDE_Step

The most important hidden parameter controlling the accuracy of the numerical PDE solver is called Number_Of_Iterations_Per_PDE_Step, which specifies how many iterations the numerical engine performs to integrate the non-linear buffering reaction in each time step. The default value is equal to 3, but in the case of very fast Ca2+ buffering reactions, its value should be increased to insure accuracy, e.g.

Number_Of_Iterations_Per_PDE_Step = 5

To determine whether the value of this parameter should be changed, run the simulation with no-flux boundary conditions and in the absense of Ca2+ extrusion and clearance mechanisms, and monitor the Charge.loss parameter, which is prompted after each simulation run.

Δ Integrating ODEs only

It is also possible to run the ODEs only (see also the table command), using instruction

mode = ODE

In this case the current definitions are ignored, and the Run command syntax becomes

Run  Time  [   accuracy  dt0  ]

The adaptive 5th-order Cash-Karp Runge-Kutta scheme is used to integrate the ODEs. Here Time is the total integration time (in ms), dt0 is the initial time step, and accuracy is the accuracy setting. Default values are: dt0 = 1 ms, accuracy = 1.0e-6. These values can be also set through assignment to the pre-defined constant names ODE.dt0 and ODE.accuracy:

model = ODE
Run 2.0
ODE.dt0 = 0.1 ; ODE.accuracy = 1e-5


Δ Graphics / data output

To plot or save the simulation results, plot instructions should be included. The type of output produced by plot instructions is controlled by the plot.method command.

To view the results in real time using the xmgr or the xmgrace programs, include the following command:

plot.method   xmgr [ rows cols ]

This instructs the program to produce an ASCII stream of commands interpretable by xmgrace/xmgr. Xmgrace / xmgr is a powerful GNU data graphing application. Optional parameters "rows" and "cols" specify the arrangement of graphs within xmgr session window. If this plot method is used, the program output should be piped into an xmgr/xmgrace session, using "calc file.par | xmgr -pipe". Only simple (point) and 1D plot types will work in this regime; other plot types will be ignored. xmgr is the recommended plot method, since in this case data goes directly into a dedicated graphics application that can be used to produce publication-quality graphs. Note that to run xmgr/xmgrace under Windows environment, an X-server must be installed first.

Note that there is a slight difference in viewport coordinates used by xmgrace vs xmgr. To correct for these differences, include one of the following two commands in the script file if you are using xmgrace:

xmgrace.style  portrait

or

xmgrace.style  landscape

In addition, you can control the margins and graph spacing through the use of the following built-in parameters (default values are shown for each parameter):

xmgr.bottom = 0.1
xmgr.top = 0.1
xmgr.left = 0.1
xmgr.right = 0.06
xmgr.xspace = 0.12
xmgr.yspace = 0.12



The default plot method is

plot.method  mute

- instructs the program to ignore all plot commands with the exception of mute, 1D.mute and 2D.mute file plots.




Data for each of the (non-mute) plots will be written to a file, if one includes the following command:

plot.print  file_prefix

File names will be composed of the file_prefix and the name of the variable to be plotted. This will work even in the plot.method mute regime. The format of the data files produced in this case is the same as the format used by the mute, 1D.mute and 2D.mute plot types (see below).

Below follows a list of implemented plot types.

Δ Plot command (simple variable plot)

Plot of a time-dependent variable (defined via "d /dt=", ":=", "Ca[x,y,z]" (or "buffer_name[x,y,z]"), "peak" or "table").

Syntax in the plot.method xmgr mode is

plot  variable1  [  variable2  variable3  ...  ]

where variable1, variable2, variable3, ... is a sequence of variable names to be plotted on the same panel of an xmgr graph.

To produce a data file instead of plotting, use

plot  mute  variable  file_name     

where variable is the name of the variable to be saved to file file_name, in the two-column format "time variable(time)". One can also use the plot.print command to write all defined simple plots to two-column data files.

The number of times each plot is refreshed is controlled by two parameters specifying respectively the resolution along the absicca (time) and ordinate (concentration value) axes; their names and default values are:

plot.steps.point  =   600      % minimal number of update steps along x-axis
update.accuracy  =   0.001      % relative resolution along y-axis


Δ "1D" plot

Plot of a concentration field (calcium or buffer) along an axis.
Syntax in the plot.method xmgr mode is

plot  1D  field  axis  coord1  coord2  [  sets  ]

where field is the concentration field label (either Ca or the buffer name defined via buffer command), axis is the axis label ("x", "y" or "z" for 3D geometry, "r" or "theta" for the conical geometry, etc.), {coord1,coord2} is the coordinate of the point in the plane perpendicular to the chosen axis, at the intersection with that axis, and the optional parameter sets denotes the number of data sets (corresponding to several previous time steps) that will be shown (default is 5).

To produce a data file instead of plotting, use

plot  1D.mute  field  axis  coord1  coord2  filename

Data will be saved as a 3-column ASCII file filename, for several different time values; in each row the 1st number is the time, the 2nd value is the space coordinate along the chosen direction, and the third value is the corresponding value of the concentration field.

The number of refresh (plot update) steps per simulation is set to 200 by default; it may be changed by assigning a value to the parameter plot.steps.1D. If this number is equal to "1", the concentration field values are only saved at the end of the simulation, as a two-column ASCII file, with the first column specifying the coordinate along the chosen axis, and the second number specifying the concentration value.

Note that for all 2D geometries, the argument coord2 should be omitted. For all 1D geometries, one should omit all three arguments axis, coord1 and coord2.

One can also use the plot.print command to write 1D non-mute plot data to a file.

Δ "2D" plot

To produce a file plot of the concentration profile in a 2D plane, use the command

plot  2D.mute  field  axis  coord  Time  filename

where field is the concentration field label (either Ca or the buffer name defined via the buffer command), axis is the name of the axis perpendicular to the plane of the 2D plot ("x", "y" or "z"), and coord is the coordinate of a point along this axis at which it intersects the chosen 2D plane. This command will produce a 3-column ASCII file; in each row the first two numbers specify the 2D coordinate in the chosen plane, and the 3rd number is the value of the field at that point, at time T. The file format is chosen for compatibility with the gnuplot plotting program, available for both Unix and Windows platforms.

Note that for all 2D geometries, the two arguments axis and coord should be omitted. The plot 2D.mute statement is equivalent to a non-mute plot 2D statement combined with the plot.print command instructing CalC to save such 2D plot to a file.

Δ "Binary" plot: saving a concentration field to a file as binary data

It may be desirable to store the entire 3D concentration field to disc in binary format, in order to plot different "slices" of concentration field later using a MATLAB script, for instance. This is accomplished by a command

plot  binary  field  file_name

where field is the concentration field (either Ca or the buffer name defined via buffer command), to be saved in file file_name. The default number of saved time points per simulation is 40; this may be adjusted by assigning a different value to internal control parameter plot.steps.binary.

The binary data is written in the following sequence: (1) an integer specifying the geometry (for Cartesian geometry, this is simply the number of dimensions); (2) one to three integers specifying the number of nodes along each axis; (3) one to three double arrays specifying the coordinates of the axis nodes; (4) finally, for each time point, a time value (double precision) followed by the field concentration as a double precision array over all nodes. A simple MATLAB script (see example readBinaryPlot.m) can be used to plot this binary data once the simulation is complete; consult this script for information about the binary file format.

To save the concentration fields at a single specified time point only, see dump plot below.

Δ "Dump" plot: saving a concentration field to a file at specific time point

Similar to "binary" plot above, but saves the binary concentration field at one specific time-point only, and without any other simulation-specific information:

plot  dump  field  T  file_name

where field is the concentration field (either Ca or the buffer name defined via buffer command), to be saved in file file_name at simulation time T. This dump plot is particularly useful when combined with the import command to read it back into another simulation. To save all of the concentration fields and ODE variables to a file, use command Export instead.

Δ Logarithmic plots

To plot data on a linear-log scale, simply append the suffix .log to a corresponding plot type, i.e.

plot  mute.log  Ca[0.1,0.1,0.1]  "Ca.dat"
plot  2D.log Ca  z  0.2


Δ Putting text on the graphs

This feature is currently implemented for the xmgr plot mode only. To print a string anywhere within the xmgr plot, use the command

plot.string  string  [  x  y  charSize  font  color  ]

where string is the string to be plotted, {x, y} is the string position within the graph in xmgr viewport coordinates (between 0 and 1), charSize is the character size (default is 0.8), font it the xmgr font code (default is 4), and color it the xmgr color code (default is 4 - blue).

If several plot.string commands are included, and the optional coordinate arguments are omitted, each successive string will be printed directly below the preceding string.

Δ More data input / output

Δ Importing 3D concentration field from a file:

Instead of initializing the concentration of a field (where field is either Ca or the buffer name defined via buffer command) with the defined resting background value or the value derived from the total value, one can read the 3D concentration field from a file (generated via the plot dump command), using

field.import  file_name

where file_name is the name of the file to be read.
To import all of the defined concentration fields (i.e., calcium plus buffers) at once, use command Import described below.

Δ Saving and retrieving data from a file: Import/Export commands

Apart from the import command and the plot dump command used to read and save the data corresponding to a single concentration field, one can also store and retrieve the entire simulation state (calcium and buffer fields plus the values of ODE variables) via a single command:

Export  T  file_name

- instructs the program to store the entire simulation state to file file_name at simulation time T.

Import  filename

- instructs the program to import the simulation state stored using command Export.

In ODE mode, the Import statement overrides all initial condition definitions in this case.

Δ Printing strings and numbers to a file: print / append commands

To write the values of parameters and variables to a file at the end of the simulation, use

print  file_name  token1  token2  token3  ...

where file_name is the name of the file to which the output is directed, and token1, token2, ... are either strings or numerical expressions. To print the data to the terminal, use the names of the standard i/o streams stderr or stdout for the file_name.

To write several lines of output to a file, after the first print command use commands

append  file_name  token1  token2  token3  ...

- same as the print command, except that file_name has to exist when this command is invoked.

Any string that can be interpreted as a numerical expression will be converted to a numerical value, so the statement

var = 4
print  stderr  "var" "=" var

will print 4=4. However, the command print stderr "var=" var will produce the expected result var=4.

The first file_name argument can be omitted, and instead a single output file can be specified for all print / append commands, using the statement

print.file  =  file_name

A carriage-return (new line) control character is added to the file file_name for each print / append command.

Note that the print statements are invoked only after the end of the simulation runs, so all expressions in the argument list will be computed based on the final values of time-dependent variables

Δ Flow control statements (preprocessor commands)

The most flexible flow control can be achieved by calling CalC from a MATLAB or Mathematica script (or any other wrapper script), and using the include statement to load the part of the script that is changed by the wrapper program. However, CalC does have some useful built-in flow-control features, described below.

The flow-control instructions are special in the sense that they are not definitions, like other commands described above, but represent tools controlling the parsing of the script itself. To this extent they are analogous to the C language pre-processor statements. Also, unlike all other CalC commands (definitions), the effect of some of these pre-processor statements may depend on their exact position in the script file.

Δ Repeating the simulation for multiple parameter values

If it is necessary to run the entire simulation many times as one or several simulation parameters are changed, one can use the "for" command. The syntax of the for cycle is:

for  parameter_name = initial_value_expression  to  end_value_expression  step  increment_expression

For instance, one may want to explore the convergence of the results as spatial resolution is increased, using:

grid N N N
for N = 10 to 40 step 2

Nested for loops are realized by using several "for" commands. The top-most for loop is the outer loop, and each successive for command implies a cycle nested within the preceding cycle.

If the "for" loop command is present, then all non-mute plot commands are ignored, and one can use the track command to monitor the results:

track  var1  var2  ....

where var1, var2, ... are names of variable to be monitored. The results will be automatically printed to a file, and graph plots will be produced also, depending on the specified plot.method. The data is output in the form of pairs of values "parameterValue  varValueAtEndOfSimulation". The file names will be composed of the file prefix specified by the plot.print command, if it is present, and followed by the name of the variable. When using xmgr plot method, arguments of the same track command will be plotted in a single graph panel. To plot different variables in separate graph panels, use several track commands.

Do not use concentration construct variables in the track statements; instead, assign it to a variable and track the variable, i.e.
X := Ca[0.5,0.5,0.5]; track X

NOTE: Although the for statement can appear anywhere within the script, like all parameter definitions, in the presence of the conditional if statement one has to make sure that it is in the accessible part of the if .. then .. else block.

Δ The conditional statement: if / then / else / endif

The conditional if statement allows conditional processing of the script (this is analogous to C preprocessor commands):

if  expression  then
  ...........
  statements
  ...........
else
  ...........
  statements
  ...........
endif

where the else clause is optional and the expression is a logical expression. Note that the if statement is really a "preprocessor" command, meaning that this logical instruction controls the parsing of the script. Therefore, the logical expression may only depend on parameters defined earlier in the script.

Although the expression only depends on the predefined time-independent parameters, the if statement allows implementing non-trivial algorithms when used in conjunction with the for loop statement, the include command, and the print/append commands. Here's how it works: at each for loop iteration the results of the simulation can be printed to a file using the print/append commands, in the format of parameter definitions (i.e. print file "x0=" x). If the include statement appearing at the beginning of the script file imports the resulting file, then at each new for loop iteration the results of the previous iteration may influence the execution. For instance, an if conditional statement may contain a clause depending on the value of the saved parameter (i.e., if x0 > 1.0 then ...).

Δ Program termination command: exit

When both "if" and "for" statements are used, it may be desirable to terminate the program once a certain condition is satisfied. For this purpose, include the statement

exit

If this statement appears in an accessible part of the script file (as determined by conditional "if" statements), the program will be terminated before any useful action is performed.

Δ Skip a for-loop iteration: continue

When a "for" statement is present, it may desirable to skip a loop iteration, for instance when a certain condition is satisfied. To skip an iteration, simply include the command

continue

The action of the continue statement is analogous to its use in other high-level programming languages.

Δ Appendix B: controlling verbosity level

To control the program's commentary, use

verbose = level

To stop the program from commenting on everything it does, use verbose = 0. The default level value is 2. Values higher than 5 are only useful for debugging.

Δ Appendix B: debugging the script

When using the if/then/else/endif statement for conditional processing of the parameter file, it may be necessary to track the parsing of the file to ensure the desired order of processing of definitions. There are two tools for such debugging:
  • Instruction

    debug

    will cause the parser to print out to the terminal every token following the debug statement, as the file is parsed.

  • The command

    echo string

    will cause the parser to print out the value of the string at parse time. If string matches a name of a parameter (either string or numerical) defined earlier in the file, the value of that parameter will be printed instead.



[back to top]
Victor Matveev

Last modified: Jan 11, 2014