This document describes how to run simulations using xolotl.

# Running simulations

`xolotl`

is designed to solve for state variables of conductance-based neuronal and network models.
The voltage across the membrane is given by the conservation of current equation.

where is the specific membrane capacitance and is a specific transmembrane current. Each current takes the form

where is the instantaneous conductance and the ionic reversal potential. In many models, the conductance is given by

where is the maximal conductance in Siemens per unit area and and are gating variables . The gating variables themselves are defined by differential equations which depend on the membrane potential. These equations are nonlinear and usually quite stiff. For these reasons, numerical simulations are required to adequately solve them.

## Changing the time step and duration

In the case where the simulation time step (`sim_dt`

) and the output time step (`dt`

) are identical, the `dt`

property only needs to be set.

```
% set the simulation time to 10 seconds
x.t_end = 10e3; % ms
% set the time-step to 1 microsecond
x.dt = 1e-3; % ms
```

The `sim_dt`

property determines the number of time-steps actually computed. The `dt`

property determines the number of time-steps output. In the following example, the computer would perform 1 million (`x.t_end/x.sim_dt`

) iterations but the output vector (or matrix) `V`

is only `1000 x nComps`

where `nComps`

is the number of compartments in the `xolotl`

object tree.

```
% perform 1e6 iterations, interpolate at a ratio of 1/1000
x.t_end = 1000; % 1000 milliseconds
x.sim_dt = 1e-3; % 0.001 milliseconds
x.dt = 1; % 1 millisecond
V = x.integrate;
```

## Closed loop vs. open loop

The `closed_loop`

flag (false or true) determines whether initial conditions should be reset before a new simulation. If `closed_loop`

is true, successive simulations will use the current state of the `xolotl`

object (e.g. the end state of the previous simulation if you run `integrate`

twice in a row).

```
% use current state of model as initial conditions
x.closed_loop = true
V = x.integrate;
```

## The outputs of `x.integrate`

### (1) Voltage or Injected current

When the `x.V_clamp`

property is not set, the first output of `x.integrate`

is the voltage trace in the form of a `nSteps x nComps`

matrix where `nSteps`

is the number of time steps and `nComps`

is the number of compartments in the model. The number of time steps in the output is determined by the simulation time and the output time step `x.t_end`

and `x.dt`

.

### (2) Calcium

The calcium trace is in the form of a `nSteps x 2*nComps`

matrix where `nSteps`

is the number of time steps and `nComps`

is the number of compartments in the model. The first `nComps`

columns are the intracellular calcium concentration (in M) for each compartment in the serialized `xolotl`

object tree. The next set of `nComps`

columns are the calcium reversal potential (in mV).

### (3) Currents

### (4) Mechanism variables

## Inject current into compartments

Injected current is mediated by the `I_ext`

property of the `xolotl`

object.

If `I_ext`

is a scalar, that amount of current in is injected into every compartment at all time.

```
% inject 0.2 nA into each compartment
x.I_ext = 0.2;
```

If `I_ext`

is a vector the length of the number of compartments, constant current will be added to each compartment in order (based on the serialized `xolotl`

object).

```
% add current to only the first compartment
x.I_ext = [0.2, 0];
```

If `I_ext`

is a matrix, it should be of size `nSteps x nComps`

where `nSteps`

is the number of time steps and `nComps`

is the number of compartments in the model. Current is added at each time step to each compartment in order (based on the serialized `xolotl`

object).

```
% add a variable current into one of two compartments
nSteps = x.t_end / x.dt;
I_ext = zeros(nSteps, 2);
I_ext(:, 1) = 0.2 * rand(nSteps, 1);
x.I_ext = I_ext;
```

## Voltage clamping compartments

## Switching solvers

```
x.solver_order = 4; % uses Runge Kutta 4
x.solver_order = 0; % default, uses exponential Euler
```

### The exponential Euler method

The exponential Euler method is a time-discrete solution to differential equations of the form:

where is the state variable and and are functions of . This equation can be solved as follows:

Rearrange equation

Divide by , multiply by

Integrate, using the relation

Given a constant term , the equation can be written

For a time step , the voltage at time can be approximated from the voltage at time , where and

This approximation is more accurate than a first order Euler method approximation, and is significantly faster than higher order (*viz.* Runge-Kutta) methods.

See Ch. 5, Dayan, P. and Abbott, LF (2001) Theoretical Neuroscience, MIT Press for more information

### The Runge-Kutta fourth-order method

The Runge-Kutta methods are extensions of forward Euler to higher derivative orders. Given a differential equation in the form

with some initial condition . The first order (Euler) approximation for given is

Euler's method is accurate per step. The Runge-Kutta fourth-order method uses four coefficients to extend this method to accuracy of per step at the cost of speed. The coefficients are, for ,

The time-evolution formula for the Runge-Kutta fourth order method is

The method is more accurate because slope approximations at fractions of are being taken and averaged. The method is slower because the four coefficients must be computed during each integration step.