# Pixel simulator

Pixel is an alternative PDE solver which uses the simple FTCS method to solve the PDE using the pixels of the geometry image as the grid.

## Simulation options

The default settings should work well in most cases, but if desired they can be adjusted by going to Advanced->Simulation options

- Integrator
the explicit Runge-Kutta integration scheme used for time integration

default: 2nd order Heun scheme, with embedded 1st order error estimate

a higher order scheme may be more efficient if the maximum allowed error is very small

see the Time integration section for more information on the integrators

- Max relative local error
the maximum relative error allowed on the concentration of any species at any pixel

default: 0.005

local means the estimated error for a single timestep, at a single point

relative means each error estimate is divided by the species concentration

making this number smaller makes the simulation more accurate, but slower

- Max absolute local error
the maximum error allowed on the concentration of any species at any pixel

default: infinite

local means the estimated error for a single timestep, at a single point

absolute means the error estimate is not normalised by the species concentration

making this number smaller makes the simulation more accurate, but slower

- Max timestep
the maximum allowed timestep

default: infinite

- Multithreading
if enabled, multiple CPU threads can be used

default: disabled

enabling this can make simulations of large models run faster

however it can also make small models run slower

- Max CPU threads
limit the maximum number of CPU threads to be used

default: unlimited

- Symbolic optimization
factor out common subexpressions when constructing the reaction terms

default: enabled

- Compiler optimization level
how much optimization is done when compiling the reaction terms

default: 3

## Spatial discretization

Space is discretized using a linear grid with spacing \(\delta x\), \(\delta y\), \(\delta z\). The concentration is defined as a 3d array of values \(c_{i,j,k}\), where the value with index \((i,j,k)\) corresponds to the concentration at the spatial point \((x = i \delta x, y = j \delta y, z = k \delta z)\).

The Laplacian is approximated on this grid using a central difference scheme

which has \(\mathcal{O}(\delta x^2)+\mathcal{O}(\delta y^2)+\mathcal{O}(\delta z^2)\) discretisation errors. Inserting this approximation into the reaction-diffusion equation converts the PDE into a system of coupled ODEs.

## Time integration

Time integration is performed using explicit Runge-Kutta integrators. Compared to implicit integrators, they are easier to implement and offer better performance (for the same timestep). However they become unstable if the timestep \(h\) is made too large, so in practice they can end up being slower than implicit methods for stiff problems, where the timestep is forced to be very small to maintain stability.

Integrators differ in their:

order of truncation error

order of embedded error estimate (if any)

number of stages (i.e. cost of a step)

region of stability (can be increased by adding more stages)

memory requirements

Implemented integrators:

- Euler
1st order solution

no error estimate

1 stage

- Embedded Heun / modified Euler
2nd order solution

1st order error estimate

2 stages

see e.g. eq (2.15) of https://doi.org/10.1016/0021-9991(88)90177-5

- Embedded Shu-Osher
3rd order solution

2nd order error estimate

3 stages

see eq (2.17) of https://doi.org/10.1016/0021-9991(88)90177-5

- RK4(3)5[3S*]
4th order solution

3rd order error estimate

5 stages

see alg.6 & tab.6 of https://doi.org/10.1016/j.jcp.2009.11.006

These integrators have three sources of error:

- Round-off error due to finite precision
mostly only relevant for high order solvers: not relevant here

- Truncation error due to finite order of integration scheme
we are generally forced by the diffusion term to make the timestep small to maintain stability

also no benefit from making the time integration errors significantly smaller than the spatial discretisation errors

so this is also typically not a concern

- Numerical instability of integrator
a problem when ODEs become stiff, e.g. high rate of diffusion, stiff reaction terms

avoiding these instabilities is our main concern

## Adaptive timestep

We use the embedded lower order solution to estimate the error at each timestep, and use this to adapt the stepsize during the integration:

RK gives us a pair of \(u_{n+1}^{(p)} = u_{n} + \mathcal{O}(h^{p+1})\) solutions

difference between \(p, p-1\) solutions gives local error of order \(\mathcal{O}(p)\)

to get the relative error we divide this by \(c = ( |c_{n+1}| + |c_{n}| + \epsilon)/2\)

we use the average of the old and new concentration, plus a small constant, to avoid dividing by zero

we do this for all species, compartments and spatial points, and take the maximum value

if this error is larger than the desired value, the whole step is discarded

the new timestep is given by \(0.98 dt_{old} (err_{desired}/err_{measured})^{1/p}\)

the 0.98 factor is slightly less than 1 to account for the higher order terms that are neglected here

it is better to have a slightly smaller timestep than to have to repeat the whole step

## Maximum timestep

For the Euler method, we don’t have an embedded lower order solution from which we can estimate of the error, so we can’t automatically adjust the stepsize. However, if we ignore the reaction terms, there is an analytic upper bound on the size of timestep that can be used for Euler (see p125 of https://doi.org/10.1017/CBO9780511781438 ), above which the system becomes unstable:

So if the user selects a timestep larger than this, the simulator automatically reduces it to the above value to avoid the system becoming unstable. Note that the system can still become unstable if the reaction terms are stiffer than the diffusion terms.

## Boundary Conditions

The boundary condition for all boundaries is the “zero-flux” Neumann boundary condition. This is implemented in the spatial discretization by setting the concentration in the neighbouring pixel that lies outside the compartment boundary to be equal to the concentration in the boundary pixel value, or equivalently by setting the neighbour of each boundary pixel to itself.

### Compartments

Each compartment is discretized, with the above boundary conditions applied for the diffusion term.

### Membranes

Reactions that take place between two compartments involve a flux across the membrane separating the two compartments. For each neighbouring pair of pixels from the two compartments, whose common boundary constitutes the membrane, the flux term is converted into a reaction term that creates or destroys the appropriate amount of species concentration in each pixel.

### Non-spatial species

A species can be ‘non-spatial’, which means that at each timestep, its time derivative is calculated as normal at each point in the compartment, but is then spatially averaged over the whole compartment. This can be used to approximate a species with a very high diffusion constant without requiring a correspondingly tiny timestep to maintain the stability of the solver.