How to write a program for pattern formation

1. How to write a program for pattern formation Hans Meinhardt Max-Planck-Institut für Entwicklungsbiologie Tübingen / Germany

2. The activator-inhibitor mechanism:how to read the equations Kybernetik 12, 30-39 (1972) Using our standard equation it should be illustrated how to write a computer program that simulates pattern formation. First it should be mentioned that such equations are easy to read: The equations describe the changes of the activator per time unit. The change of the activator concentration depends on a non-linear positive feedback of the activator on its own production; the inhibitor slows down the activator production. The activator displays normal decay (- ra a)and diffusion. The inhibitor production also depends on the activator concentration in a non-linear fashion. The diffusion of inhibitor is assumed to be much more rapid than that of the activator (Db >> Da). The activator-independent activator production bais required to initiate the self-enhancing process, e.g. during regeneration.

3. Flow diagram To program such a reaction, one has to discretize the process; the space into units (such as cells) and the time into small time steps. The computer has to calculate the concentration change for each cell at each time step and add this change to the existing concentrations. In this way we obtain the new distribution after a certain time unit. Many such “iterations” gives the complete time course. Setting parameters Setting initial concentration Calculating changes / new distributions Display ?

4. How to code such a reaction? The new concentrations in the cells have to be computed sequentially from cell i = 1 to cell n (assuming a linear array of cells) : 1 i n This is the concentration change per time unit according to our original equation: This equation has to be re-written as a difference equation: Thus, the new concentration of a cell i ( i = 1…n) at the time t +1 (i.e., after a time unit) is given by the old concentration plus the concentration change according to the equation: ai,t+1 = ai,t + si (ai,t2 + ba) /bi,t – ra ai,t + Da( (ai-1,t - ai,t ) + (ai+1,i – ai,t) )

5. Handling diffusion Concentrations are stored in arrays, e.g. the activator in ax(i) , i = 1…n; the actual local concentrations in the cell i are stored in a, b.., e.g. a =ax(i) [Left cell => al (see below)] Diffusion: Da (aleft cell – aactual cell ) + Da (aright cell – aactual cell ) => Da ((al – a)+ (ax(i+1) – a)) al ax(i) ax(i+1) 1 i n Important: For calculating diffusion use always non-updated concentrations ! In a sequential updating, the left cell will be already updated, the right cell not. This leads to severe errors. To avoid this problem, the non-updated concentration ais saved as al. During calculation for the next cell, al is used as activator concentration in the left cell, which is the correct non-updated concentration. Note that the present cell will be the left cell when the next cell is calculated.

6. Boundary condition Special conditions have to be introduced for the cells at the boundaries as they do not have left or right neighbors. Impermeable boundaries can be introduced by assuming virtual cells to the left and to the right of the terminal cells that have the same concentrations as the terminal cells: Boundary condition: impermeable al = ax(1) ax(n+1) = ax( n) n 1 (no exchange by diffusion takes place if two cells have the same concentrations, i.e., nothing leaks out at the boundaries)

7. Danger: numerical instability Be careful: the numerical value of the diffusion constant should not exceed a certain limit, otherwise numerical instabilities will occur. In other words, the time steps have to be small enough. In the two examples below a central cell has the concentration one, the neighbors zero. With Da = 1 the central cell would become negative (-2), the concentration of the two others would be one, which is nonsense! In linear calculations Da should be not larger than about 0.4. Da = 1 Da = 0.33

8. The code for a BASIC-program The updating has to proceed in a sequential way in a loop, i is the running index, going from 1 to n (number of cells). First, the array values (e.g., ax(i) for the activator) are stored as local concentrations (e.g., a). Then, the new concentration depending only on the decay and diffusion is calculated (e.g., olddecaydiffuseA), followed by the change due to the reaction. FOR i = 1 TO n' i = current cell, n = right cell a = ax(i) 'local activator-concentration b = bx(i) 'local inhibitor-concentration s = sx(i) 'local ability to perform the pattern-forming reaction ‘REM olddecaydiffA contains the new concentration considering diffusion and decay olddecaydiffA = a - RA * a + DA * ((al - a) + (ax(i + 1) - a)) olddecaydiffB = b - RB * b + DB * ((bl - b) + (bx(i + 1) - b)) ‘ adding the concentration change based on the reaction: ax(i) = olddecaydiffA + s * (a * a + ba) / b ’ updating of the activator bx(i) = olddecaydiffB + s * a * a ‘ updating the inhibitor concentration al = a: bl = b 'The not yet changed concentrations are used as left cell 'concentrations when the subsequent cell is calculated NEXT i

9. Activator-Inhibitor Mechanism: Thus, the concentration change per time unit as used in the original equation… …and the code for the new concentrations after a time unit loo very similar and are thus easy to read: ax(i) = olddecaydiffA + s * (a * a + ba) / b bx(i) = olddecaydiffB + s * a * a

10. To obtain a running program, few additional things have to be done: 1. Setting parameters KT = 100 'Number of displays KP = 20 ' number of iterations between the displays ‘( KT * KP = number of iterations in total) n = 40 ' number of cells DA = .01 ' Diffusion of the activator RA = .02 ' Removal rate of the activator BA = .001 ' Activator-independent activator production rate DB = .4 ' Diffusion of the inhibitor RB = .03 ' Removal rate of the inhibitor ’ 2. ------- Setting the initial conditions -------------------------- KR = 0: AA = 1.2: GA = 1 REM KR = 1: AA = 1: GA = 1 ' removal of the REM{ark} leads ' initiation by random fluctuations FOR i = 1 TO n s(i) = RA * (.99 + KR / 100 * RND)'"Source density" = Production of the 'activator, proportional to the decay rate + KR% fluctuation ax(i) = GA'general initial activator concentration bx(i) = 1'general initial inhibitor concentration NEXT ax(n/3) = AA ‘ one cell may have a different concentration

11. To obtain a running program, few additional things have to be done: ‘3. setting boundary conditions (impermeable) al = ax(1): bl = bx(1): 'al is the concentration in the cell left of the ' actual cell. Left-most cell = virtual cell with the same concentration ax(n + 1) = ax(n)' concentration in a virtual cell to the right of the bx(n + 1) = bx(n)' right-most cell is equal to the concentration in the ‘right terminal cell

12. This is the essential part of the program plott: ' REM ----------------Plot ------------- LINE (20, 45)-(620, 50), 1, BF x1 = 20' Position of the first rectangle FOR i = 1 TO n x2 = x1 + idx afl = 51 + ax(i) * fa: LINE (x1, 51)-(x2, afl), 2, BF'green=activator LINE (x1, afl)-(x2, 450), 15, BF 'remaining part white: erase old plot bfl = 51 + bx(i) * fb: LINE (x1, bfl)-(x2, bfl + 10), 12, BF sfl = 51 + s(i) * fs: LINE (x1, sfl)-(x2, sfl + 3), 1, BF x1 = x2: NEXT IF itot >= KT OR INKEY$ > "" GOTO alldone FOR iprint% = 1 TO KP' Calculations between plots ‘-----------' BOUNDARIES IMPERMEABLE ------------ al = ax(1): bl = bx(1): 'al is the concentration in the cell left of the ' actual cell. Left-most cell = virtual cell with the same concentration ax(n + 1) = ax(n)' concentration in a virtual cell to the right of the bx(n + 1) = bx(n)' right-most cell is equal to the concentration in the REM ---------- Reactions ------ FOR i = 1 TO n' i = current cell, n = right cell a = ax(i) 'local activator-concentration b = bx(i) 'local inhibitor1-concentration REM ----- Calculation of a new activator and inhibitor concentration in cell i olddecaydiffA = a * (1 - RA) + DA * ((al - a) + (ax(i + 1) - a)) olddecaydiffB = b * (1 - RB) + DB * ((bl - b) + (bx(i + 1) - b)) ax(i) = olddecaydiffA + s(i) * (a * a + ba) / b bx(i) = olddecaydiffB + s(i) * a * a al = a: bl = b 'The not yet changed concentrations are used as left cell 'concentrations in the subsequent cell NEXT i NEXT iprint% Itot = itot + 1: GOTO plott

13. This is the complete program for copy and paste '------------------------------------------------------------ ' A simple program to simulate biological pattern formation ' (C) Hans Meinhardt, MPI Tuebingen DEFDBL A-G DEFDBL O-Z DEFINT H-N DEFINT H-N ' A simple program to simulate biological pattern formation ' (C) Hans Meinhardt, MPI Tuebingen ' i = 1 ... n= Nsumber of cell imax = 640: DIM ax(imax), bx(imax), cx(imax), s(imax) start: 'Parameter KT = 100 'Number of displays KP = 20 ' 'number of iterations between the displays 'KT*KP = number of iterations in total n = 40' number of cells DA = .01' Diffusion of the activator RA = .02' Removal rate of the activator BA = .001 ' Activator-independent activator production rate DB = .4' Diffusion of the inhibitor RB = .03' Removal rate of the inhibitor fa = 60: fb = 40' Scaling factors for display KR = 0: AA = 1.2: GA = 1 REM KR = 1: AA = 1: GA = 1 ' removal of the REM{ark} leads ' initiation by random fluctuations REM ----------- Initial conditions -------------------------- FOR i = 1 TO n s(i) = RA * (.99 + KR / 100 * RND)'"Source density" = Production of the 'activator, proportional to the decay rate +KR% fluctuation ax(i) = GA ‘ general initial activator concentration bx(i) = 1 ‘ general initial inhibitor concentration NEXT ax(n / 3) = AA ‘ initial perturbation SCREEN 12: t = TIMER: CLS 'Initialization of graphic and timing WINDOW (1, 1)-(640, 480) ' coordinate system LINE (1, 1)-(640, 480), 15, BF 'background white igraph = 1 continuo: ' if calculation is continued... idx = 600 / n: fs = 350 / RA 'Pixel-size of a cell itot = 0 plott: ' REM ----------------Plot ------------- LINE (20, 45)-(620, 50), 1, BF x1 = 20' Position of the first rectangle FOR i = 1 TO n x2 = x1 + idx afl = 51 + ax(i) * fa: LINE (x1, 51)-(x2, afl), 2, BF'green=activator LINE (x1, afl)-(x2, 450), 15, BF 'remaining part white: erase old plot bfl = 51 + bx(i) * fb: LINE (x1, bfl)-(x2, bfl + 10), 12, BF sfl = 51 + s(i) * fs: LINE (x1, sfl)-(x2, sfl + 3), 1, BF x1 = x2: NEXT IF itot = 0 THEN LOCATE 30, 1: PRINT SPACE$(79); LOCATE 30, 1: PRINT " That is the initial situation, green: Activator, red: Inhibitor "; t = TIMER: DO UNTIL TIMER - t > 2: IF INKEY$ > "" THEN EXIT DO: LOOP' END IF DO UNTIL TIMER - t > .1: LOOP: t = TIMER'slow down if computer is too fast IF itot >= KT OR INKEY$ > "" GOTO alldone LOCATE 30, 1: PRINT "green = Activator, red = Inhibitor, blue = Source density; any key: stop"; FOR iprint% = 1 TO KP' Calculations between plots REM ------' boundaries tight ------------ al = ax(1): bl = bx(1): 'al is the concentration in the cell left of the ' actual cell. Left-most cell = virtual cell with the same concentration ax(n + 1) = ax(n)' concentration in a virtual cell to the right of the bx(n + 1) = bx(n)' right-most cell is equal to the concentration in the REM ---------- Reactions ------ FOR i = 1 TO n' i = current cell, n = right cell a = ax(i) 'local activator-concentration b = bx(i) 'local inhibitor1-concentration REM ----- Calculation of a new activator and inhibitor concentration in cell i olddecaydiffA = a * (1 - RA) + DA * ((al - a) + (ax(i + 1) - a)) olddecaydiffB = b * (1 - RB) + DB * ((bl - b) + (bx(i + 1) - b)) ax(i) = olddecaydiffA + s(i) * (a * a + ba) / b bx(i) = olddecaydiffB + s(i) * a * a al = a: bl = b 'The not yet changed concentrations are used as left cell 'concentrations in the subsequent cell NEXT i NEXT iprint% itot = itot + 1: GOTO plott alldone: LOCATE 30, 1: PRINT SPACE$(79); LOCATE 30, 1: PRINT "c = continue; s = start again, other key: END"; a$ = "": DO UNTIL a$ > "": a$ = INKEY$: LOOP' wait for key SELECT CASE a$ CASE "c", "C": GOTO continuo CASE "s", "S": GOTO start END SELECT TheEND: LOCATE 30, 1: t = TIMER PRINT ".... a program for simulation of biol. pattern formation;(c) Hans Meinhardt"; DO UNTIL TIMER - t > 1: LOOP '-------------------End of the Program------------------------------- On our website one can find information how to find the (free) FreeBasic compiler and an appropriate editor http://www.eb.tuebingen.mpg.de/de/forschung/emeriti/hans-meinhardt/biuprom.html

14. …and these are screenshots from the display Green = activator a Red = Inhibitor b Blue = source density of competence s (in this case without random fluctuations, thus a local perturbation is required for initiation.

15. For more information, visit our website http://www.eb.tuebingen.mpg.de/meinhardt Academic Press (1982) Awailable for download