A Simple Basic Program to simulate the Patterns on Oliva Porphyria
This short program provides an example for a program that allows the simulation the a shell pattern. It is written for PowerBasic and for compilers of the Microsoft BASIC family. The program is also compatible with the new BASIC for WINDOWS, FREEBASIC. For an executable file for a PC and the source code in a compressed form click here . The program can be changed and recompiled with compilers that are available on the web. Details are given here . The picture below shows a pattern generated by this program together with two real shells for comparison.
'------------------------------------------------------------- DEFDBL A-G DEFDBL O-Z DEFINT H-N ' Models for the simulations of the color pattern on the shells of mollusks ' see also: Meinhardt,H. and Klingler,M. (1987) J. theor. Biol 126, 63-69 ' see also: H.Meinhardt: "Algorithmic beauty of sea shells" ' (Springer Verlag) (c) H.Meinhardt, Tübingen 'This is a short version of a program for the simulations of the color 'patterns on tropical sea shells, here 'Oliva porphyria'. 'An autocatalytic activator a(i) leads to a burst-like activation 'that is regulated back by the action of an inhibitor b(i). The life 'time of the inhibitor is regulated via a hormone c, that is 'homogeneously distributed along the growing edge. Whenever the number 'of activated cells cells become too small, active cells remain activated 'until backwards waves are triggered 'The program runs with the interpreter QBASIC, but this is very slow. 'Better are the compiler Power Basic, Microsoft QB 4.5, Professional 'Basic 7.1, and Visual Basic für DOS. A compiled version is included ' i = 1...kx < imax = cells at the growing edge imax = 1200: DIM ax(imax), bx(imax), zx(imax) RANDOMIZE TIMER ' By different fluctuations, 'simulation will be slightly different KT = 460 'Number of displays ' KT * KP = number of iterations in total KP = 12 'number of iterations between the displays ( = lines on the screen) kx = 630 'number of cells dx = 1 'width of a cell in pixel; with kp=6 ; kx=315 and dx=2 => ' simulation in a smaller field DA = .015 'Diffusion of the activator RA = .1 'Decay rate of the inhibitor BA = .1 'Basic production of the activator SA = .25 'Saturation of the autocatalysis DB = 0 'Diffusion of the inhibitor RB = .014 'Decay rate of the inhibitor SB = .1 'Michaelis-Menten constant of the inhibition RC = .1 'Decay rate of the hormone start: REM ----------- Initial condition -------------------------- FOR i = 1 TO kx ax(i) = 0 'Activator, general initial concentration bx(i) = .1 'Inhibitor, general initial concentration zx(i) = RA * (.96 + .08 * RND)'Fluctuation of the autocatalysis NEXT i C = .5 'Hormone-concentration, homogeneous in all cells i = 10: FOR j = 1 TO 30 'initially active cells ax(i) = 1: i = i + 100 * RND + 10: IF i > kx THEN EXIT FOR NEXT DAC = 1! - RA - 2! * DA ' These constant factors are used again and again DBC = 1! - RB - 2! * DB ' therefore, they are calculated only once DBCC = DBC' ' at the begin of the calculation SCREEN 12 WINDOW (1, 1)-(640, 480) continuo: v1 = TIMER LINE (1, 1)-(640, 480), 15, BF'background blue x0 = 5: y1 = 475 'Initial position of the first line itot = 0 'Number of iteration calculated so far FOR itot = 0 TO KT FOR iprint% = 1 TO KP ' Begin of the iteration REM ----- --- Boundary impermeable A1 = ax(1) ' a1 is the concentration of the actual cell. since this B1 = bx(1) ' concentration is identical, no diffusion through the border. ax(kx + 1) = ax(kx) ' Concentration in a virtual right cell bx(kx + 1) = bx(kx) BSA = 0! ' This will carry the sum of all activations of all cells REM ---------- Reactions ------ FOR i = 1 TO kx' i = actual cell, kx = right cell AF = ax(i) 'local activator concentration BF = bx(i) 'local inhibitor concentration AQ = zx(i) * AF * AF / (1! + SA * AF * AF) 'Saturating autocatalysis ' Calculation of the new activator and inhibitor concentration in cell i: ax(i) = AF * DAC + DA * (A1 + ax(i + 1)) + AQ / (SB + BF) ' 1/BF => Action of the inhibitor; SB = Michaelis-Menten constant bx(i) = BF * DBCC + DB * (B1 + bx(i + 1)) + AQ 'new inhibitor conc. BSA = BSA + RC * AF 'Hormone production -> Sum of activations A1 = AF ' actual concentration will be concentration in left cell B1 = BF ' in the concentration change of the next cell NEXT i C = C * (1! - RC) + BSA / kx ' New hormone concentration , 1/kx=normalization RBB = RB / C ' on total number of cells 'RBB => Effective life time of the inhibitor DBCC = 1! - 2! * DB - RBB ' Change in a cell by diffusion NEXT iprint% ' and decay. Must be recalculated since ' lifetime of the inhibitor changes REM ----------------Plot ------------- y1 = y1 - 1 'Next plot, one line below LINE (x0, y1)-(x0 + dx * kx, y1), 15 'Background white FOR ix% = 1 TO kx 'Pigment is drawn when a is higher than a threshold IF ax(ix%) > .5 THEN LINE (x0 + dx * (ix% - 1), y1)-(x0 + dx * ix%, y1), 6 NEXT ix% IF INKEY$ > "" THEN EXIT FOR NEXT itot v2 = v1 - TIMER LOCATE 30, 1: PRINT "c = continue; s = a new start, all other keys = End (c) H.Meinhardt "; PRINT USING "##.##"; TIMER - v1; DO: a$ = INKEY$: IF a$ > "" THEN EXIT DO LOOP IF a$ = "c" GOTO continuo IF a$ = "s" GOTO start END '------------------ End of the Program --------------------------------------------