The simulation scripts described in this chapter is available at STEPS_Example repository.

This chapter introduces how to model and simulate surface diffusion systems. In STEPS, surface diffusion means movement of molecules between triangle elements within a patch and is used to model e.g. mobility of surface receptors within membranes.

In practice, simulating diffusion in surfaces is analogous to simulating diffusion in volumes as introduced in Simulating Diffusion in Volumes, replacing volume system with surface system and compartment with patch. This section will demonstrate the simple case of free diffusion from a source on a large circular surface, which is somewhat analogous to free diffusion in a sphere described in Simulating Diffusion in Volumes. Therefore the code in this chapter has obvious similarities to that in Simulating Diffusion in Volumes and we will not dwell on familiar concepts, but will focus on the key differences instead.

After running the spatial stochastic simulation in familiar solver 'Tetexact', we are then introduced to spatial deterministic solver 'TetODE' (steps.solver.TetODE).

In the deterministic limit, C, the number of diffusing molecules per unit area, at a distance
r from a point source on an infinite plane surface
at time *t* is (Crank, J. (1975) The Mathematics of Diffusion. Oxford: Clarendon Press):

where *N* is the total number of injected molecules and *D* is the diffusion constant (in units $m^{\text{2}}/s$ if length units are metres and time units seconds).

As for previous models we will create a Python script to run the surface diffusion model. Since we're now familiar with the concept of building Python scripts to run STEPS models we will no longer use the Python prompt syntax.

First we import STEPS modules and other modules `pylab`

and `math`

:

```
In [1]:
```import steps.model as smodel
import steps.geom as stetmesh
import steps.utilities.meshio as smeshio
import steps.rng as srng
import steps.solver as solvmod
import pylab
import math

```
In [2]:
```# Number of iterations; plotting dt; sim endtime:
NITER = 100
# The data collection time increment (s)
DT = 1.0
# The simulation endtime (s)
INT = 21.0
# Number of molecules injected in centre
NINJECT = 1000
# The diffusion constant for our diffusing species (m^2/s)
DCST = 0.08e-12

Now we move on to constructing the biochemical model which, as in Simulating Diffusion in Volumes, we will organise into a small function that returns the important steps.model.Model container object. First we'll look at the complete function:

```
In [3]:
```def gen_model():
mdl = smodel.Model()
A = smodel.Spec('A', mdl)
ssys = smodel.Surfsys('ssys', mdl)
diff_A = smodel.Diff('diffA', ssys, A, DCST)
return mdl

`surface system`

. This is the only difference between creating a diffusion rule for a volume (compartment)
or a surface (patch): if the diffusion rule belongs to a volume system (as in Simulating Diffusion in Volumes) the diffusion rule will determine how the specified molecular
species diffuse in any compartments to which that volume system is added, and if the diffusion rule belongs to a surface system (steps.model.Surfsys ),
as in this example, it specifies how the species will diffuse in any patch to which that surface system is added.
There is only one steps.model.Diff object that takes care of both eventualities (volume or surface diffusion), adapting behaviour depending on whether
it belongs to a volume system or surface system. Any given species may of course diffuse in both volumes and surfaces (by different diffusion coefficients) and
therefore may appear in both volume and surface diffusion rules in the same model.
Since this is the same steps.model.Diff object as introduced in Simulating Diffusion in Volumes its construction is similar and the same general behaviour applies.
To reiterate, the arguments to the constructor are, in order: the usual identifier string,
a reference to the parent surface system, a reference to the molecular species to which this diffusion
rule applies, the diffusion constant (which is an optional parameter to the object constructor) given in s.i. units ($m^{2}/s$). This will become the
value for this diffusion rule wherever it appears in the model. Function steps.model.Diff.setDcst is another option for setting the diffusion constant.

The next step, as usual, is to create the geometry for the simulation which, as in the previous chapter Simulating Diffusion in Volumes, will require a tetrahedral mesh because this is a diffusion model. This time we will use a flat disk-shaped mesh (of radius 10 microns) and run the diffusion on one of the circular faces.

Please see previous chapter Simulating Diffusion in Volumes for a more detailed discussion of tetrahedral mesh generation and import. Functions that have already been described in that chapter will not be described here in detail.

The first step, within a function to group all geometry creation code, is to import the tetrahedral mesh shown in the above figure using function steps.utilities.meshio.loadMesh, which returns a tuple with a steps.geom.Tetmesh object as the 0th element.

```
In [4]:
```mesh = smeshio.loadMesh('meshes/coin_10r_1h_13861')[0]

```
In [5]:
```ntets = mesh.countTets()
comp = stetmesh.TmComp('cyto', mesh, range(ntets))

```
In [6]:
```# Find the indices of the exterior surface triangles
alltris = mesh.getSurfTris()
# Sort patch triangles as those of positive z
patch_tris = []
for t in alltris:
vert0, vert1, vert2 = mesh.getTri(t)
if (mesh.getVertex(vert0)[2] > 0.0 \
and mesh.getVertex(vert1)[2] > 0.0 \
and mesh.getVertex(vert2)[2] > 0.0):
patch_tris.append(t)

`patch_tris`

containing the collection of circular face triangles we are interested in, we can create the
patch for our spatial simulation: a steps.geom.TmPatch object. The required arguments to the constructor are, in order: a string identifier, a
reference to the steps.geom.Tetmesh container object, a sequence of the triangles that comprise this patch, and an 'inner' compartment.

```
In [7]:
```patch = stetmesh.TmPatch('patch', mesh, patch_tris, icomp = comp)
patch.addSurfsys('ssys')

`icomp`

is a required argument to the steps.geom.TmPatch constructor. If the patch is connected to two compartments
one must be labelled the 'inner' compartment (constructor argument `icomp`

) and one the 'outer' compartment (constructor argument `ocomp`

) and the object construction might look something like

```
In [8]:
``````
# ER_membrane = stetmesh.TmPatch('ER_membrane', mesh, Er_memb_tris, icomp = endoplasmic_reticulum, ocomp = cytosol)
```

```
In [9]:
```# Find out how many triangles are in the patch
patch_tris_n = len(patch_tris)
# Create array to store radial distances
trirads = pylab.zeros(patch_tris_n)
# Create array to store triangle areas
triareas = pylab.zeros(patch_tris_n)

```
In [10]:
```# Find the central triangle
# First find the tetrahedron connected to the central triangle
ctetidx = mesh.findTetByPoint([0.0, 0.0, 0.5e-6])
# Find this tetrahedron's neighbours
ctet_trineighbs = mesh.getTetTriNeighb(ctetidx)
# Find which of these (4) neighbours is in the patch
ctri_idx=-1
for t in ctet_trineighbs:
if t in patch_tris:
ctri_idx = t

`patch_tris`

list, the distance of the barycentre of the triangle to the centre
point (in microns) using function steps.geom.Tetmesh.getTriBarycenter, and the triangle area (in square microns) using function steps.geom.Tetmesh.getTriArea.
We will record this information in the arrays `trirads`

and `triareas`

respectively:

```
In [11]:
```# Now find the distance of the centre of each tri to the central tri
cbaryc = mesh.getTriBarycenter(ctri_idx)
for i in range(patch_tris_n):
baryc = mesh.getTriBarycenter(patch_tris[i])
r2 = math.pow((baryc[0]-cbaryc[0]),2) + \
math.pow((baryc[1]-cbaryc[1]),2) + \
math.pow((baryc[2]-cbaryc[2]),2)
r = math.sqrt(r2)
# Convert to microns
trirads[i] = r*1.0e6
triareas[i] = mesh.getTriArea(patch_tris[i])*1.0e12

Finally, we return the important information found within this function.

Our complete geometry function is:

```
In [12]:
```def gen_geom():
mesh = smeshio.loadMesh('meshes/coin_10r_1h_13861')[0]
ntets = mesh.countTets()
comp = stetmesh.TmComp('cyto', mesh, range(ntets))
alltris = mesh.getSurfTris()
# Sort patch triangles as those of positive z
patch_tris = []
for t in alltris:
vert0, vert1, vert2 = mesh.getTri(t)
if (mesh.getVertex(vert0)[2] > 0.0 and mesh.getVertex(vert1)[2] > 0.0 and mesh.getVertex(vert2)[2] > 0.0):
patch_tris.append(t)
# Create the patch
patch = stetmesh.TmPatch('patch', mesh, patch_tris, icomp = comp)
patch.addSurfsys('ssys')
patch_tris_n = len(patch_tris)
trirads = pylab.zeros(patch_tris_n)
triareas = pylab.zeros(patch_tris_n)
# Find the central tri
ctetidx = mesh.findTetByPoint([0.0, 0.0, 0.5e-6])
ctet_trineighbs = mesh.getTetTriNeighb(ctetidx)
ctri_idx=-1
for t in ctet_trineighbs:
if t in patch_tris:
ctri_idx = t
# Now find the distance of the centre of each tri to the central tri
cbaryc = mesh.getTriBarycenter(ctri_idx)
for i in range(patch_tris_n):
baryc = mesh.getTriBarycenter(patch_tris[i])
r2 = math.pow((baryc[0]-cbaryc[0]),2) + math.pow((baryc[1]-cbaryc[1]),2) + math.pow((baryc[2]-cbaryc[2]),2)
r = math.sqrt(r2)
# Convert to microns and square microns
trirads[i] = r*1.0e6
triareas[i] = mesh.getTriArea(patch_tris[i])*1e12
return mesh, patch_tris, patch_tris_n, ctri_idx, trirads, triareas

`Tetexact`

We're now ready to run the simulation and collect data. First we call the two functions `gen_model`

and `gen_geom`

and store references to the important objects
they return, the steps.model.Model and steps.geom.Tetmesh container objects, and the patch triangles information (being consistent with their
names in the `gen_geom`

function to minimise confusion):

```
In [13]:
```model = gen_model()
tmgeom, patch_tris, patch_tris_n, ctri_idx, trirads, triareas = gen_geom()

We then need to create our random number generator object, as usual for stochastic simulations:

```
In [14]:
```rng = srng.create('mt19937', 512)
rng.initialize(234)

And then we can create a `Tetexact`

spatial stochastic solver object:

```
In [15]:
```# Create solver object
sim = solvmod.Tetexact(model, tmgeom, rng)

```
In [16]:
```tpnts = pylab.arange(0.0, INT, DT)
ntpnts = tpnts.shape[0]
# Create the data structure: iterations x time points x tri samples
res = pylab.zeros((NITER, ntpnts, patch_tris_n))

`triareas`

data
to record, in the `res`

results array, molecule number per square micron. The complete simulation loop is:

```
In [17]:
```from __future__ import print_function # for backward compatibility with Py2
# Run NITER number of iterations:
for j in range(NITER):
if not j%10: print("Running iteration ", j)
sim.reset()
sim.setTriCount(ctri_idx, 'A', NINJECT)
for i in range(ntpnts):
sim.run(tpnts[i])
for k in range(patch_tris_n):
res[j, i, k] = sim.getTriCount(patch_tris[k], 'A')/triareas[k]

```
```

`numpy.mean`

function as in
previous chapters:

```
In [18]:
```res_mean = pylab.mean(res, axis = 0)

As we are quite familiar with plotting STEPS output now we won't dwell on the details. As in Simulating Diffusion in Volumes we wish to also plot the analytical solution alongside our STEPS output for comparison, where the solution is of a similar form to that for unbounded volume diffusion. We create two functions to take care of all our plotting needs. The first function plots the results:

```
In [19]:
```def plotres(res_mean, tidx):
if (tidx >= INT/DT):
print("Time index is out of range.")
return
pylab.scatter(trirads, res_mean[tidx], s=10)
pylab.xlabel('Radial distance ($\mu$m)')
pylab.ylabel('Concentration (/$\mu$m$^2$)')
t = tpnts[tidx]
pylab.title('Unbounded surface diffusion. Time: ' + str(t) + 's')
plotanlyt(t)
pylab.xlim(0,10)
pylab.ylim(0)
pylab.show()

`plotres`

function) plots the analytical solution, this time as number of molecules per square micron:

```
In [20]:
```def plotanlyt(t):
segs = 100
anlytconc = pylab.zeros((segs))
radialds = pylab.zeros((segs))
maxrad = 0.0
for i in trirads:
if (i > maxrad): maxrad = i
maxrad *= 1e-6
intervals = maxrad/segs
rad = 0.0
for i in range((segs)):
anlytconc[i]=(NINJECT/(4*math.pi*DCST*t))* \
(math.exp((-1.0*(rad*rad))/(4*DCST*t)))*1e-12
radialds[i] = rad*1e6
rad += intervals
pylab.plot(radialds, anlytconc, color = 'red')

```
In [21]:
```pylab.figure(figsize=(10,7))
plotres(res_mean, 20)

```
```

`A`

in individualtriangles in STEPS (black dots) is plotted with the analytical solution from the above equation (red).There is a small discrepancy near the centre of the surface due to the injection of molecules into a finite area inSTEPS whereas a point source is assumed for the analytical solution. STEPS output shows some noise due to the effects of stochastic surface diffusion.

`TetODE`

Another option for spatial simulation is to use the deterministic solver 'TetODE' (steps.solver.TetODE). TetODE shares many similarities with Tetexact in terms of model and geometry construction operating on the same tetrahedral meshes, but solutions are deterministic. TetODE uses CVODE (http://computation.llnl.gov/casc/sundials/description/description.html) for solutions. Although solutions are therefore very different between solver Tetexact and TetODE, in terms of simulation construction there are only a few implementation differences. Therefore, we can use almost the exact same code as already introduced to run the deterministic unbounded surface diffusion model, with a few changes highlighted below.

Firstly, since TetODE solutions are deterministic, we do not need to run more than one iteration so we set

```
In [22]:
```NITER = 1

Finally, the `reset`

solver function is not available for TetODE (in part because only one simulation iteration need be run per model), so we remove the call to `sim.reset()`

.

from our simulation loop. Theses changes are all we minimally need to do to switch this simulation to deterministic mode using solver TetODE. However, there are two important additions to this solver, the functions steps.solver.TetODE.setTolerances and steps.solver.TetODE.setMaxNumSteps. To understand what these functions do requires a little background of how CVODE works. Although there will only be a brief explanation here, thorough descriptions are available in CVODE documentation (http://computation.llnl.gov/casc/sundials/documentation/cv_guide.pdf).

Solving STEPS models in CVODE
requires supplying information of all the variables in a STEPS simulation at any time as a state vector to the CVODE solver. The variables in STEPS
are the molecular species, which have unique populations in individual mesh elements (tetrahedrons and triangles) meaning that
the state vector can be rather large (number_volume_specs*number_tetrahedrons + number_surface_specs*number_triangles). STEPS
must also supply a function that describes the rate of change of each of these variables with time depending on other variables in the system. CVODE then finds
approximate solutions (here STEPS choses the recommended Adams-Moulton formulas with functional iteration) when the system
advances in time.

To do this it takes a number of 'steps', each time estimating the local error and comparing to tolerance conditions: if the test fails, step size is reduced, and this is repeated until tolerance conditions are met. This means that there is a tradeoff between accuracy and simulation speed- with a high tolerance, steps sizes will be large and few steps will have to be taken to advance the simulation some amount of time though accuracy will be low, with low tolerance, steps sizes will be small so a large number of steps will be taken to advance the simulation although accuracy will be high. Therefore, the tolerance is an important consideration both for accuracy and efficiency.

STEPS users can control the tolerances with function steps.solver.TetODE.setTolerances. Two different types of tolerance are specified: relative tolerance and absolute tolerance, and in STEPS both are scalars. Relative tolerance controls relative errors so that e.g. $10^{-3}$ means that errors are controlled to 0.1% (and it is not recommended to go any higher than that). Absolute tolerances can be useful when any components of the vector approach very small numbers when relative error control becomes meaningless. The absolute values in the internal state vectors within TetODE are the (fractional) number of molecules per tetrahedron or triangle, so if a user specifies an absolute tolerance of $10^{-3}$ it means that populations within tetrahedrons and triangles will be accurate to within 1/1000th of a molecule! In TetODE only one value each for absolute tolerance and relative tolerance can be specified, and will be applied to all species in all locations in the system. The default value for both absolute tolerance and relative tolerance is $10^{-3}$ .

We set tolerances with a call to function steps.solver.TetODE.setTolerances: the first function argument is
absolute tolerance, the second is relative tolerance. In this example we set an absolute tolerance of $10^{-3}$ and relative tolerance of $10^{-4}$ with `sim.setTolerances(1e-3, 1e-4)`

.

`sim.run(1)`

.

`sim.setMaxNumSteps(50)`

.

We can now run the simulation and plot the results.

```
In [23]:
```# Create solver object
sim = solvmod.TetODE(model, tmgeom)
sim.setTolerances(1e-3, 1e-4)
sim.setMaxNumSteps(50)
tpnts = pylab.arange(0.0, INT, DT)
ntpnts = tpnts.shape[0]
# Create the data structure: iterations x time points x tri samples
res = pylab.zeros((NITER, ntpnts, patch_tris_n))
# Run NITER number of iterations:
for j in range(NITER):
sim.setTriCount(ctri_idx, 'A', NINJECT)
for i in range(ntpnts):
sim.run(tpnts[i])
for k in range(patch_tris_n):
res[j, i, k] = sim.getTriCount(patch_tris[k], 'A')/ \
triareas[k]
res_mean = pylab.mean(res, axis = 0)
pylab.figure(figsize=(10,7))
plotres(res_mean, 20)

```
```