chapter3a

contents software data java models

Implementing reaction diffusion equations - one modelers perspective
Robert C. Cannon

This tutorial covers the development of three small programs sphere.i, wave.i and spark.i for modeling reaction-diffusion systems. To run the examples you need to install Yorick from the software section of the CD.
For those who are familiar with it, mathematics can be a very convenient language for conveying ideas. for those who are not, and even some who are, it can be rather intimidating. What's more, for biological applications it is sometimes unnecessary or even misleading. We should bear in mind that many biological phenomena are discrete on a rather coarse scale compared with the coarseness common in physics. For example a big neuron may have 10^5 ion channels, but contains more than 10^16 water molecules: continuum fluid mechanics is an excellent approximation and worth a big mathematical investment, whereas writing smooth continuous partial differential equations for the membrane potential will simply never be as comprehensive or relevant. So, while maths is undeniably useful, it isn't the whole story. Given that frequently the goal is to model one discrete system, the biology, with another discrete system, a numerical model implemented in software, one can sometimes skip a lot of the associated mathematics. Here we take that approach to that task of modeling calcium signalling in neurons. The focus is firmly on the numerical methods and software. Programs are developed which will allow the user to perform the exercises in chapter 3. The presentation follows the order in which the programs were designed and written.
What language to use?

The first question to answer when approaching a new modeling problem (apart that is from the obvious ones - can I get someone else to do it?; is there something out there that does it already ) is what language to use. The choice is pretty large: perl, fortran, c, lisp, java, matlab, IDL, python, Maple, Mathematica, Yorick and more. All can be used without too much trouble, except perhaps perl, for the type of numerical models considered here.
Here we will use Yorick because it is convenient, powerful and free. You can install it under MacOS, Windows or unix from the software distribution on the CD. In fact the above list roughly follows my personal view of increasing order of suitability for small numerical computations. I hope by the end you'll see why Yorick comes at the top. If you use it under unix, be sure to also install rlterm as described in the distribution. It aids usability enormously.
Numerical methods for ODEs

The next question, for those in the know, is what sort of numerical methods to use, what packages or libraries are available to do the work for you, and generally how not to reinvent the wheel. At the cost of a little reinvention, almost everything will be done from scratch here to show how it all works.
Numerical methods for differential equations fall into two broad camps: finite difference and finite element methods. The problems they handle overlap extensively, but the approaches are very different. In the first case one writes down the continuous equations to be solved, then splits the solution space up into discrete points and derives a new set of difference approximations to the original equations using only those points. Then you find a numerical scheme to solve these difference equations. This is a rather technical procedure full of considerations or truncation errors and order of approximations. It is the preferred choice of many applied mathematicians, who delight in proving that in the limit of small steps their approximate scheme does indeed converge to the solution of the continuous equations. Happily it can be dispensed with in much biological modeling where the continuous equations are only an approximation to a discrete system in the first place.
The finite element method on the other hand proceeds with a leap of intuition by looking for a discrete system which we feel ought to approximate the system represented by the continuous equations which we haven't bothered to write down. One can then often solve the discrete system exactly, so the accuracy issue resides in the question "how close is my discrete picture to the real thing?", a question which one can generally ignore, unless your intuition is particularly bad. The finite element method is much favoured by engineers and physicists, who know enough maths to realise when it can be dispensed with. Biologists are divided but are frequently in awe of the density of incomprehensible equations in finite difference papers and try to emulate it with varying levels of success. Generally, however, the greatest advantage of the finite difference approach is when you trust your maths but not your feeling for the physics. Here we'll go with the presumed majority view, skip the maths and go for the finite element approach.

A first program

For the impatient, here is the whole program for the full equations or the reaction-diffusion system. To run it, start up yorick, then type
#include "sphere.i"
sphere
Type help to read the documentation to see what parameters you can set and so on.
Lets go through it line by line, in the order it was written. Numbers refer to the numbering in comment lines which begin "//".
1 Declare a function to contain the code. This isn't essential. You could just type a file full of yorick commands, but it is neater. For now, ignore the arguments: they can be put in later as it becomes clear what is necessary.
2, 3 One good thing about yorick is that you don't have to remember the order of the arguments. It can have keyword arguments as here, specified as "name=value". When you do this, they should also have a default value to use if the user hasn't set their own. Again, this block can be filled in later when we know what is needed.
4 Time to think about the problem. Lets focus on Box3.3 which requires a numerical solution for the calcium diffusion and buffering around a point source. So space and time will have to be subdivided into little intervals for the finite elements. Since everything is happening at the center and at the beginning, it seems a good idea to take small intervals near the center and at early times, largest ones further out and later on. Line 4 does this for time, and line 5 for space. In each case the points are taken at uniform intervals of the logarithm. This also introduces the first few parameters we will need: how long to run the simulation for, trun, and what spatial range to cover; from rmin to rmax. These should be added to the keyword list in the function call and given reasonable default values.
Now we have a set of radial points, describing shells around the source, but what we need are finite elements. Line 6 finds the centers of the spaces between shells, and line 7 their volumes. Note the advantages of an array based language (ie, most of the above list except fortran, c, and java) are already clear. It takes one line to handle all the points. These lines also indicate the convenience of Yorick: the index range functions zcen and dif return new arrays giving the zone centers and differences respectively.
8 We have the space and time discretization. Now we have to think about the physics. At each element, there are going to be variables for the concentrations of calcium, the buffers, and the Ca+buffer complexes. Two sorts of process are needed: reactions between the various species, and diffusion from one element (or cell from now on) to its neighbours. The second seems easier. How fast is this diffusion going to be for a particular par of neighbouring cells? One could go back to the equations, but if you ever read a derivation of diffusion equations, it will probably begin "Imagine you split the volume into small cells..." and proceed to write down the equations for finite cells and then take a limit to get the continuous case. Lets skip this little diversion. The diffusion rate is going to be proportional to the are of contact between the cells, inversely proportional to the distance between their centers, and for any particular species we'll have to multiply in its diffusion constant later. That's it. Line 8 gives an array of diffusion rates from cell 1 to and from 2 and from 2, to and from 3 etc.
9 Next come the reaction equations. These are more difficult, so to put them off as far as possible, we can tidy up a few other things that will be necessary. In particular, we can count the buffers and set up the initial conditions - the concentrations of each species in each cell. Buffers haven't figured so far. What do we need? Well, each has an initial concentration, a forward rate for binding calcium, a backward rate for releasing the calcium again, and a diffusion constant. So the buffers will be provided as a two-dimensional array, each row having these four elements. For generality, lets not set the number of buffers but work it out from the array provided. And put the buffers= option in the argument list. Now we can set up the concentration array, conc with 1 + 2 * nbuf elements (Ca, B1, CaB1, B2, CaB2, ...) at each of _tt(np) points. The indgen business fills just the buffer variables with the corresponding bits of he buffers array.
10 This was an afterthought. bc0 and bc1 are flags change the boundary conditions. bc0==1 means the Ca concentration at the center should be fixed at sigma throughout.
11 Now it is time to think about the reactions. We could just plow on and wait till the last minute when evaluating the equations and use the data in the buffers array directly. But here we have a chance to forget that this is a buffering problem and convert it into a generic reaction scheme problem. For one thing the code will be more useful later, and for another it will look cleaner. So what is needed for a general reaction scheme? Just the rates of the reactions, rreac, and a list of which species are involved as reactants or products in each one - ireac. For each row of ireac the first three elements give the reactants and the second three the products, with 0 where there are fewer. Three of each is unnecessary here but lets leave it for now. Each buffer gives rise to two reactions: one for binding and one for unbinding.
Now everything is set up. It is time to think about what numerical scheme to use to solve the equations. But before that, we can put in the outer loop that will take the timesteps (12) and even fill in some little details to do with drawing graphs as the program runs (21). The xlog and ylogkeyword seem like a good idea, to specify whether to plot the distance and concentrations directly or their logs.
Another handy feature (22) would be to make it pause after each step, so we can see how it is going. In yorick this can be done with the read command which just waits for you to hit return. And after that, we can finish the loop "}", and end the function "}" which is always good, even if the guts are still missing.
Numerical Methods

There is a vast literature on numerical methods for differential equations. Happily, most of it is devoted to highly accurate or highly efficient schemes which we can comfortably ignore. This problem is small and quick to solve. Almost any scheme would work. One could even try a forward Euler scheme: evaluate reaction rates, diffusion rates and increment the concentrations accordingly. This could do the trick, but even in small programs there is a limit to the amount of waste that you can put up with! If one day you want to see what happens with a fast buffer, you would have to take very small steps and keep checking that the step size is small enough to get reasonably close to the answer. It is much better to invest in something a bit more sophisticated which will not complain at fast buffers or rapid diffusion. See chapter 1 for a discussion of numerical methods. For now, lets aim at the backward Euler method. It couldn't be more different from the forward Euler: it is unconditionally stable. What you have to do, is take a step not according to the slope at the start, but according to the slope at the (as yet unknown) endpoint. In other words, you must find an endpoint to the step such that when you evaluate the gradients there, they give you the slope of the step you have just taken.
Now, this may seem a bit tricky and it is, (a bit). But the structure is easy to set up. What we must do is make an initial guess of the changes to be made in this step (line 13). Then we need some procedure (15 -- 19) to adjust these changes until the above condition holds. There are many ways of doing this too, but a pretty standard choice is the Newton-Raphson method. What you do is write the system of equations to be solved in the form F(V1,v2,...)=0, evaluate the current deviations from 0 (line 14) and then evaluate their derivatives with respect to the increments on this step, and move your trial increments to where you'd expect this to be zero if the whole thing were linear. If it is indeed linear, then you get there in one step. Otherwise you must repeat the process until the equations are close enough to zero. There are three steps: evaluate the equations(14, 18); evaluate their derivatives (16); and project to zero (17).
Equations

The equations are relatively straightforward now that we have sorted out the buffers into the ireac and rreac arrays, and the diffusion rates into rdif. Each elution states that the change in concentration dconc (32). + the sum of reactions turning this species into something else - the sum of reactions turning other things into this should be zero. Thus, each reaction contributes to various equations as in the block at (32). As mentioned above the diffusion terms need a product of the concentration difference with the neighboring cell, the mobility, the timestep, and the geometrical term rdif. There is also a factor of the volume in the denominator since the concentration change is the flux into the cell divided by its volume. That's it. The yorick syntax with constructions like (-,) adds dimensions to arrays so they be come conformable, allowing all the diffusion terms to be evaluated in two lines. This type if shortcut is what yorick is best at and helps dramatically in eliminating mistakes: if what you want to do can't be done concisely then you are probably trying to do the wrong thing.
Derivatives

The Newton Raphson scheme needs to know how each of the equations will vary as you change the trial concentration increments for the step - it needs the derivatives of the result of the equations function with respect to the input dconc. There are two ways to do this - numerically or analytically. To do it numerically, you would just change each variable in turn a tiny amount, reevaluate the equations, and there is your answer. If the equations are very complicated this is probably the only hope of getting near the right answer in a reasonable amount of time (or use a computer algebra package that will write the code for you). In our case, though, the elution routine only has nineteen statements so we can probably write the derivatives without too much trouble.
First of all, how should the results be stored? We could fill one big matrix, with the (i,j) element containing the derivative of the i'th equation with respect to the j'th variable. This would work fine for small systems, but could get awkward to handle if you want very high spatial resolution or a lot of species. The problem is that most of the elements would be zero: the size of the matrix would be proportional to the square of the number of mesh cells, whereas the number of non-zero elements is just linear in the number of cells. This is because each equation only links a cell to its neighbours: derivatives with respect to the rest of the variables are zero. This is a simple case of a very common situation and can be handled with a tri-block-diagonal matrix which is made of sets of three blocks around the diagonal. (see the picture in the rdDerivs routine). Each set of blocks corresponds to one mesh cell, and the rows within the blocks correspond to the different species at that mesh cell. The central block lies directly on the diagonal and contains derivatives with respect to concentrations in that mesh cell. This block contains all the reaction terms. The neighbouring blocks contain derivatives with respect to the concentrations in neighbouring cells so they are diagonal and just contain the diffusion terms.
We won't go through the derivatives routine in detail. Suffice it to say, all you do is look at each term in the equations routine and for each concentration it involves, put the term divided by that concentration in the appropriate slot in the matrix.
Matrix solution

Given the current values of the equations R and the matrix of their derivatives M then the solution X to M X = R is what you should subtract from the variables to get the equations nearer to zero. This isn't necessarily obvious, but what else can you do with M and R? And why subtract it? Well, if you add it and look at the average error as the Newton Raphson iteration proceeds then you'll see it get bigger, not smaller. This is a very useful programming technique : when you've finished but the program does the wrong thing, permute all the signs you're not sure of until it works.
In general, you would go to some linear algebra package to solve a matrix equation but this matrix has a particular structure we imposed. The function tbdMatrixSoln does a form of block-by-block gaussian elimination, and calls a standard LU decomposition routine from the yorick library to invert individual blocks. This last part is important. You may be wondering whether this yorick code is going to be much slower that a low level language like c or fortran. The answer is no, because, for large reaction networks, all the work is done within the LUsolve routine. For small networks we don't care since it will only take a second or two. In fact, yorick's LUsolve is so well optimised, that if one wrote a new compiled routine just for this class of matrix in c it would almost certainly be slower than the present combination of 99% yorick and one function optimised by someone else.
Running the program

The program so far is dimensionless. That avoids a lot of possible mistakes, but you have to make sure any parameters you give it are consistent. Lets measure time in milliseconds and distance in microns. This means we have to specify diffusion constants in micron^2/ms. The concentration unit is still free so we can fix it at microMolar. Then a value of 1 for the sigma parameter for calcium injection is the injection rate which would raise the concentration of a one micron^3 box by one microMolar in one ms. What current is this? The box is 10^-15 litres, so we're getting 10^-21 moles per millisecond, 10^-18 moles per second or 2 * 10^-13 coulombs per second (96500 coulombs in a mole, each ion carries two charges), or 0.2 pA. So, with these units, (microns, milliseconds, microMolar), sigma=5 corresponds to 1pA. This is the default value of sigma.
Now, assuming you've installed yorick, to run the program just type
#include "sphere.i"
sphere
You should see the successive calcium profile as it solves the equations. Ca is in white. The blue and black lines are the dummy buffer that does nothing for now. By default this shows a log-log plot. To see a linear plot run
sphere, xlog=0, ylog=0
Note that you can use the mouse to rescale the graph - see the yorick instructions, or just try clicking at random. You can now answer the questions in Box3.3. It may help to zoom in on a particular range in radius and limit the total time with the trun= keyword.
To explore the buffering, try
sphere, buffers= [1, 1, 1, 1]
and change each of the buffer parameters up or down by a factor of 10. Remember that the order is: initial concentration, binding rate (reactions per ms); unbinding rate; mobility. The blue line shows the free buffer and the black the bound buffer.
You may notice that the concentration near the center is rather high. This is because the equations break down for small radius: for a fixed injection rate the concentration tends to infinity as radius approaches zero. Of course, it is meaningless to talk of concentration in a cube of side 1nm. Such a cube would contain (10^-8)^3 moles, or on average 1 ion of a 1M solution. At 1mM a cube of side 10nM would only contain a single ion so one should be wary of making rmin so small.
A more physical approach might be to say that when it is open the channel quickly ensures that the volume just inside the pore reaches the same concentration as outside the channel. This can be done in the simulation by setting bc0=1 That has the effect of fixing the concentration of the innermost element at sigma. A similar trick can be played with the outermost point. By default it is free, so it is as though there is a fixed wall at rmax - not such a bad approximation since not many calcium channels find themselves in regions without a wall of some sort within a few microns. But you may sometimes want an infinite sink at the outer boundary. To do this, set (bc1=1).
Excess buffer and rapid buffering approximations

Now, suppose you wish to see how close the full solutions are to various approximations that you might want to make in practice to simplify the equations. In particular consider the idea that there may be so much buffer that its depletion is negligible, or that it binds so rapidly that you don't need to consider the binding reaction.
At first it may appear that we need another program to evaluate what happens with the approximate equations. Fortunately not. We can legitimately assume that the excess buffer approximation (EBA) gives the behaviour you'd get if there really was excess buffer, and that the rapid buffer approximation (RBA) gives the solution if the buffering really was rapid. So to work out what they do, all one need do is add lots of buffer or increase the rate constants. At this stage the time invested in writing an implicit solution pays off. Other methods might break down for rapid reaction rates, but the implicit Euler handles them well. However, it may have problems if the change in a single step is very large. To overcome this, the keyword parameter t0 can be used to set the size of the first step.
A suitable set of keywords to start investigating the RBA and EBA is:
sphere, buffers= [50, 0.1, 1.0, 0.0], bc1=1, rmax=3, rmin=0.03, sigma=50, xlog=0, ylog=0, trun=10
This does not actually give you the steady state, but the situation after 10ms. If you really want the steady state, try trun=100 or trun=1000. Note that for these latter cases you really do have to set the outer boundary condition, bc1=1. If you try bc1=0 as though the system was confined to a 3 micron sphere, you see it filling up with calcium and exhausting the buffer entirely. This type of check is useful to ensure that the elegance of your steady state equations is not leading you into unphysical situations. HEAD3(Calcium waves in linear geometry)
So far we have a program for reaction diffusion systems in spherical geometry. How should this be converted to linear geometry? There are only three lines which concern the geometry, those labelled 5, 7 and 8, so it is sufficient to replace the first with simple linear increments and get rid of the pi r^n terms in th others. Since there are more modifications coming, this is done in the file wave.i. A few other modifications can be made at the first pass. Since we're after a continuous event, the logarithmic time intervals are no longer useful, nor is rmin. Indeed, since all cells are the same, we can get rid of the cell by cell volume and diffusion terms entirely.
Wave propagation requires two compartments at each point with different calcium concentration and buffers in each compartment. It also needs reactions for uptake and release of calcium, the latter gated by a calcium dependent process. How much of this can be implemented in the existing reaction scheme? The two compartments are not a problem. We can just duplicate the variables in each cell, and if we do not include reactions between species in different compartments they will not interact. Uptake and release are also straightforward reaction process. The uptake, being enzyme mediated could be modelled with Michaelis Menten kinetics but (Chapter 2) this is just equivalent to two of the normal reactions we have already so there is no need to complicate the code. A potential problem comes from the possibility that the volumes of the two compartments at any given location are not equal. This would mean that when calcium moves between them the concentration of one changes more than that of the other. This does not fit in the framework we have so far: it would either need a reaction of the form Ca_1 <--> f Ca_2 where f != 1, or for the geometry to be included explicitly. But there is so much degeneracy already in the problem that it is probably no great loss to force the compartments to have equal volumes. That is, even if the real system had different volumes, its behaviour could be reproduced in our equal-volume model just by scaling the parameters in one of the compartments.
The last potential problem is the gating of calcium release by ryanodin receptors. Such channels are usually modelled with kinetic schemes in which one or more of the transitions is a function of the calcium concentration. In reality, of course they (like almost everything) are reaction schemes: Ca ions have to bind to the closed state to move it into the open state, so the reaction scheme formalism covers these too. But it is stretching the point a bit to regard the passage of calcium through an open channel as a reaction process. In a sense it is, because the Ca ion lodges in the pore before moving out, so one could treat it the same way as an enzyme reaction with very fast kinetics. This is not very efficient since all we really want to do is scale the rate of the release process according to the fraction in the open state. More importantly, would it work? Can the program handle the reaction rates it would require to get the mechanism working? Is there an alternative? Well, yes. What would the system say to a reaction of the form A + C <--> B + C ? It seems a perfectly legitimate reaction in the terms of the equations, and the rate at which A turns into B is going to be proportional to C - just what we want, where C is the open channel state and A and B are the calcium concentrations on the two sides. So the little calcium domain program appears only to be half a dozen lines different from a Calcium Induced Calcium release (CICR) wave model! All it needs is the right set of reactions, so lets begin by making ireac and rreac arguments. We'll keep the buffers array and also include them in ireac as before.
A wave propagating reaction scheme

The variables we will always need are: Ca concentrations in each compartment, variables 1 and 2; the two states of the channel between them, 3 and 4; and perhaps some enzyme concentration to control the uptake of Ca by one compartment, but lets try without first. Any buffers can be expressed as before except that now we need two separate arrays for the two compartments. The file wave.i includes all these modifications. To run this set of reactions, you should #include "wave.i" in yorick.
As a first attempt, lets set the initial concentrations to [0., 1000, 1, 0] - all Ca in the store, all receptors closed. Now, the uptake reaction just takes species 1 to species 2, so its line in ireac would be [1, 0, 2, 0].
wave, c0 = [1., 1000, 1., 0.], ireac = [[1,0,2,0]], rreac = [0.1]
The receptor transition involves external Ca binding: [3, 1, 4, 0].
wave, c0 = [1., 1000, 1., 0.], ireac = [[1,0,2,0], [3,1,4,0], [4,0,3,1]], rreac = [0.1, 1., 0.1]
The receptor gates the free diffusion of calcium between compartments: [1, 4, 2, 4] and its inverse.
wave, c0 = [0., 100, 1., 0.], ireac = [[1,0,2,0], [3,1,4,0], [4,0,3,1], [1,4,2,4], [2,4,1,4]], rreac = [0.1, 1., 0.1, 0.1, 0.1]

This set of parameters does indeed show a propagating wave of release, but the channels remain open so the system does not recover. For that another state is required in the channel for inactivation. We need another species in the simulation for the inactivated channels, and another reaction taking species 4, the open state of the channel, to species 5, the inactivated state.
wave, c0 = [0., 40, 1., 0., 0.], ireac = [[1,0,2,0], [3,1,4,0], [4,0,3,1], [1,4,2,4], [2,4,1,4], [4,0, 5, 0]], rreac = [0.1, 1., 0.1, 0.3, 0.3, 1.1]

By setting the bufa keyword, you can explore the effects of fixed or mobile buffers on the propagation speed of the wave. So far we have concentrated on setting up and solving the reaction diffusion problem, with only minimal use of graphics. The next task is to improve the display of results so it is clearer what is going on. This will be combined with a few minor modifications to allow the software to handle non-uniform calcium release.
Spark mediated propagation

The final version of the program, spark.i, includes a few more keyword parameters allowing the initial profile of calcium release channels to be specified. It has a parameter dx for the distance between release sites, and the rmax and nr keywords are displaced in favour of nsite and npps (numper of points per site), the latter being the number of cells separating release sites. The more points there are, the more accurate the calculation, but the slower it runs. Frequently it is most convenient when working on a model to use low resolution and relatively poor accuracy, only running slower, finer grained calculations at the end. It is rare that you learn much more from a vary accurate calculation, given the uncertainties in many of the parameters and the gross approximations to the true physical system that have to be made.
Given the new parameters, the modifications necessary to implement the spark mediated propagation amount to only a few lines: at (6) a list of the indices of the release sites are generated, and these are used at (10) to set the initial concentrations for these points. At the same time we add a new array (11a) to accumulate the calcium concentration history, add each new line of data as it arrives (21a) and plotting it as an image.
The following set of parameters which made a continuous wave before now shows spark mediated propagation
spark, c0 = [0., 40, 1., 0., 0.], ireac = [[1,0,2,0], [3,1,4,0], [4,0,3,1], [1,4,2,4], [2,4,1,4], [4,0, 5, 0]], rreac = [0.1, 1., 0.1, 0.3, 0.3, 1.1]
In order to eliminate flickering in the image window, spark.i calls yorick's animate function, which causes the full contents of the window to be prepared off-screen before it is shown.
Since the new display is rather pretty, it seems a shame not to be able to use it with the previous continuous wave simulation. How best to do this? One could duplicate the extra lines in the previous program but this would mean that any new enhancement would also have to be duplicated for it to work in both. Better, since they are already 98% the same, is to modify spark.i to cover all the functionality in wave.i and throw the earlier wave.i away. In fact this has already been done: if instead of specifying dx and nsite you specify nr as before, it will do a continuous wave simulation.
Fire-diffuse-fire models

One last option that may be of interest is to cut out some of the complexity of the CICR. A complex but incomplete model of a complex system is frequently worse (less informative or more misleading) than a simple model of the system. About the simplest model one can make of CICR is to set a threshold at which the stores instantaneously dump all their calcium which then diffuses and may trigger the next site. Can this also be done with the same program? Certainly. If you were starting from scratch to make a fire-diffuse-fire you would not bother with all the machinery for handling reaction schemes, but since we have that already, and we are not particularly concerned with speed, lets use the existing code and add a few blocks of instructions for the fire-diffuse-fire case.
First of all there is a keyword argument fdf to set the threshold. If this is null or zero then we will do the previous spike or wave calculation. Next, (11b) we need an array to record which sites have already fired, or, slightly more conveniently, which are still able to fire. Finally (20a) we check if any firable sites have exceeded threshold this step,and equilibrate their calcium if they have.
A fire-diffuse-fire simulation can be run with the command:
spark, c0=[0.,40], ireac=[1,0,2,0], rreac=[0.1], fdf=0.01

Conclusion

This tutorial has been concerned with computational and numerical side of building models. We have ended up with short (~170 lines) but versatile program capable of handling a wide variety of reaction diffusion systems by specifying command line parameters only. With the addition of a few lines here and there it could be extended to encompass special features often added to models in the interest of realism, such as more complex expressions for reaction rates or initial conditions.
For many applications it is far from the most efficient solution, but it works. In the last few models, the only diffusing species has been calcium, whereas the calculation included the possibility that any of them might have had non-zero diffusion constants. Likewise the behaviour of the rapid or excess buffer approximations was explored simply by making the buffering very rapid or providing an excess of buffer while still solving the full equations. In either case, the computational task could have been reduced, but at the expense of making the software more complicated and error prone. Unless the point is reached where the simple general purpose program really is too slow, it is much better to use a few more minutes or even hours of computer time than risking a few more days of your own!
We have only touched briefly on the issue of parameterisation. The models implemented are dimensionless, and the parameters used in examples are only very loosely based on those which are appropriate to neuronal calcium dynamics. The behaviour,however, is qualitatively similar even with these parameters, and the general trends are the same. The programming is done: in order to do the exercises in the book you must find the relevant regions of parameter space for the distances, concentrations, times, reaction rates and diffusion constants.