/* --------------------------------------------------------------------- This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. --------------------------------------------------------------------- */ #define VERSION "v7" // Magnification of .bmp files now possible. // Stack enzymes in a different way to avoid // an unintended extra synergy. // Syn, by Allan Crossman // A cellulase model from Edinburgh iGEM 2011 // // A note on how SDL is used: if Syn is compiled from the Makefile, // USE_SDL becomes defined, which changes everything. But if Syn is // just compiled with "gcc syn.c" this doesn't happen. #define SYNERGY_LEFT 0 // These are just booleans for whether #define SYNERGY_RIGHT 1 // to use synergy or not. #define ENZYMES 20 // Enzymes of each type on each side. #define HALFWIDTH 155 #define HEIGHT 155 #define SDL_MAG 3 // Magnification of the graphics in SDL windowed mode. #define BMP_MAG 1 // Magnification of the graphics in BMP files. #define BIOSE_INHIB_LEFT (float) 1.0 // Strength of cellobiose inhibition #define BIOSE_INHIB_RIGHT (float) 1.0 // of exoglucanase (per molecule). #define ENDO_REQ 5 // How many sugars are needed on each side // of a bond before endoglucanase can cut it. // When to end the simulation (0 for never)... #ifndef USE_SDL #define MAX_ITERATIONS 40000 // ... in No-SDL mode. #else #define MAX_ITERATIONS 0 // ... in SDL mode. #endif // How often to draw a .bmp file (0 for never)... #ifndef USE_SDL #define BMPTICK 1000 // ... in No-SDL mode. #else #define BMPTICK 0 // ... in SDL mode. #endif #define REPORTTICK 1000 // How often to print count of glucose & cellobiose molecules. #define DELAY 20 // Wait this many milliseconds between frames. SDL mode only. // Might work OK as 0, but some machines choke if it's too low. #define MARGIN 4 // Don't reduce this below 1, a margin is expected. // (e.g. by the chain-length counting logic) #include #include #include #include #include #include #include #ifdef USE_SDL #include #endif #define EMPTY 0 #define SUGAR 1 #define BOND 2 #define ENDO 3 #define EXO 4 #define GLU 5 // Colours... #define ENDO_R 255 #define ENDO_G 255 #define ENDO_B 0 #define EXO_R 0 #define EXO_G 255 #define EXO_B 255 #define GLU_R 255 #define GLU_G 255 #define GLU_B 255 #define SUGAR_R 0 #define SUGAR_G 200 #define SUGAR_B 0 #define BOND_R 150 #define BOND_G 200 #define BOND_B 150 #ifdef USE_SDL #define setrgb(screen, x, y, red, green, blue) putpixel(screen, x, y, SDL_MapRGB(screen->format, red, green, blue)) #endif #define WIDTH ((((HALFWIDTH) * 2))) struct cellulase_struct { int x; int y; }; #ifdef USE_SDL // Some globals... SDL_Surface * virtue; // Main virtual screen. SDL_Event event; #endif // The main data array... char world[WIDTH][HEIGHT]; // A pixmap to draw to .bmp files if needed... char worldtodraw[WIDTH * BMP_MAG][HEIGHT * BMP_MAG]; // Cellulases... struct cellulase_struct endo[ENZYMES * 2]; struct cellulase_struct exo[ENZYMES * 2]; struct cellulase_struct glu[ENZYMES * 2]; // Exit cleanly, shutting down SDL. void cleanexit (int exitstatus) { #ifdef USE_SDL SDL_Quit(); #endif // If we ever use sound we must also do Mix_CloseAudio(); exit(exitstatus); } #ifdef USE_SDL // Manage relevant events. Should test for key presses here. void manageevents (void) { while (SDL_PollEvent(&event)) { if (event.type == SDL_QUIT) { cleanexit(0); } if (event.type == SDL_KEYDOWN && event.key.keysym.sym == SDLK_ESCAPE) { cleanexit(0); } } return; } #endif // Random number between 0 and 1. float rnd (void) { return (float)rand() / RAND_MAX; } // Return integer between 0 and max inclusive. int intrnd (int max) { int result; result = rnd() * (max + 1); result %= max + 1; return result; } void randomwalk (struct cellulase_struct * c) { c->x += intrnd(2) - 1; c->y += intrnd(2) - 1; // Prevent enzymes falling off the left or right side... if (c->x == 0) c->x = 1; if (c->x == WIDTH - 1) c->x = WIDTH - 2; // Prevent enzymes crossing the centre... if (c->x == HALFWIDTH - 1) c->x = HALFWIDTH - 2; if (c->x == HALFWIDTH) c->x = HALFWIDTH + 1; // Prevent enzymes falling off the top or bottom... if (c->y == 0) c->y = 1; if (c->y == HEIGHT - 1) c->y = HEIGHT - 2; return; } // How many sugars to the left of this spot, plus 1 if there's a sugar on this spot // But zero if this spot is empty. Call with result = 0. int leftchain (int x, int y, int result) { if (world[x][y] == EMPTY) return 0; if (world[x][y] == BOND) result = 0 + leftchain(x - 1, y, result); if (world[x][y] == SUGAR) result = 1 + leftchain(x - 1, y, result); return result; } // How many sugars to the right of this spot, plus 1 if there's a sugar on this spot // But zero if this spot is empty. Call with result = 0. int rightchain (int x, int y, int result) { if (world[x][y] == EMPTY) return 0; if (world[x][y] == BOND) result = 0 + rightchain(x + 1, y, result); if (world[x][y] == SUGAR) result = 1 + rightchain(x + 1, y, result); return result; } #ifdef USE_SDL // Puts a pixel on the SDL surface. Surface should really be locked. // From Bill Kendrick's GPL game Vectoroids (http://www.newbreedsoftware.com/vectoroids) void putpixel(SDL_Surface * surface, int x, int y, Uint32 pixel) { int bpp; Uint8 * p; // Assuming the X/Y values are within the bounds of this surface... if (x >= 0 && y >= 0 && x < surface->w && y < surface->h) { // Determine bytes-per-pixel for the surface in question: */ bpp = surface->format->BytesPerPixel; // Set a pointer to the exact location in memory of the pixel // in question: p = (((Uint8 *) surface->pixels) + // Start at beginning of RAM (y * surface->pitch) + // Go down Y lines (x * bpp)); // Go in X pixels // Set the (correctly-sized) piece of data in the surface's RAM // to the pixel value sent in: if (bpp == 1) { *p = pixel; } else if (bpp == 2) { *(Uint16 *)p = pixel; } else if (bpp == 3) { if (SDL_BYTEORDER == SDL_BIG_ENDIAN) { p[0] = (pixel >> 16) & 0xff; p[1] = (pixel >> 8) & 0xff; p[2] = pixel & 0xff; } else { p[0] = pixel & 0xff; p[1] = (pixel >> 8) & 0xff; p[2] = (pixel >> 16) & 0xff; } } else if (bpp == 4) { *(Uint32 *)p = pixel; } } return; } #endif #ifdef USE_SDL // Pixel putting routine that scales stuff up... void drawmag (SDL_Surface * screen, int x, int y, int magnify, int red, int green, int blue) { int i; int j; for (i = 0; i < magnify; i++) { for (j = 0; j < magnify; j++) { setrgb(screen, x * magnify + i, y * magnify + j, red, green, blue); } } return; } #endif void drawbmp (char * filename, int magnify) { unsigned int headers[13]; FILE * outfile; int extrabytes; int paddedsize; int x; int y; int n; int i; int j; int red, green, blue; // First we need to update the pixmap, so we know what pixels to draw as enzymes... for (x = 0; x < WIDTH; x++) { for (y = 0; y < HEIGHT; y++) { for (i = 0; i < magnify; i++) { for (j = 0; j < magnify; j++) { worldtodraw[x * magnify + i][y * magnify + j] = world[x][y]; } } } } for (n = 0; n < ENZYMES * 2; n++) { for (i = 0; i < magnify; i++) { for (j = 0; j < magnify; j++) { worldtodraw[endo[n].x * magnify + i][endo[n].y * magnify + j] = ENDO; worldtodraw[exo[n].x * magnify + i][exo[n].y * magnify + j] = EXO; worldtodraw[glu[n].x * magnify + i][glu[n].y * magnify + j] = GLU; } } } extrabytes = 4 - ((WIDTH * magnify * 3) % 4); // How many bytes of padding to add to each // horizontal line - the size of which must // be a multiple of 4 bytes. if (extrabytes == 4) extrabytes = 0; paddedsize = ((WIDTH * magnify * 3) + extrabytes) * HEIGHT * magnify; // Headers... // Note that the "BM" identifier in bytes 0 and 1 is NOT included in these "headers". headers[0] = paddedsize + 54; // bfSize (whole file size) headers[1] = 0; // bfReserved (both) headers[2] = 54; // bfOffbits headers[3] = 40; // biSize headers[4] = WIDTH * magnify; // biWidth headers[5] = HEIGHT * magnify; // biHeight // Would have biPlanes and biBitCount in position 6, but they're shorts. // It's easier to write them out separately (see below) than pretend // they're a single int, especially with endian issues... headers[7] = 0; // biCompression headers[8] = paddedsize; // biSizeImage headers[9] = 0; // biXPelsPerMeter headers[10] = 0; // biYPelsPerMeter headers[11] = 0; // biClrUsed headers[12] = 0; // biClrImportant outfile = fopen(filename, "wb"); // // Headers begin... // When printing ints and shorts, we write out 1 character at a time to avoid endian issues. // fprintf(outfile, "BM"); for (n = 0; n <= 5; n++) { fprintf(outfile, "%c", headers[n] & 0x000000FF); fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8); fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16); fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24); } // These next 4 characters are for the biPlanes and biBitCount fields. fprintf(outfile, "%c", 1); fprintf(outfile, "%c", 0); fprintf(outfile, "%c", 24); fprintf(outfile, "%c", 0); for (n = 7; n <= 12; n++) { fprintf(outfile, "%c", headers[n] & 0x000000FF); fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8); fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16); fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24); } // // Headers done, now write the data... // for (y = HEIGHT * magnify - 1; y >= 0; y--) // BMP image format is written from bottom to top... { for (x = 0; x <= WIDTH * magnify - 1; x++) { if (worldtodraw[x][y] == SUGAR) { red = SUGAR_R; green = SUGAR_G; blue = SUGAR_B; } else if (worldtodraw[x][y] == BOND) { red = BOND_R; green = BOND_G; blue = BOND_B; } else if (worldtodraw[x][y] == ENDO) { red = ENDO_R; green = ENDO_G; blue = ENDO_B; } else if (worldtodraw[x][y] == EXO) { red = EXO_R; green = EXO_G; blue = EXO_B; } else if (worldtodraw[x][y] == GLU) { red = GLU_R; green = GLU_G; blue = GLU_B; } else { red = 0; green = 0; blue = 0; } // Also, it's written in (b,g,r) format... fprintf(outfile, "%c", blue); fprintf(outfile, "%c", green); fprintf(outfile, "%c", red); } if (extrabytes) // See above - BMP lines must be of lengths divisible by 4. { for (n = 1; n <= extrabytes; n++) { fprintf(outfile, "%c", 0); } } } fclose(outfile); return; } // Return 1 if the spot is the bond of a cellobiose molecule, return 0 otherwise. // Assumes bonds that aren't flanked by sugars can't exist. int is_cellobiose_bond (int x, int y) { if (x < 0 || x >= WIDTH || y < 0 || y >= HEIGHT) return 0; if (world[x][y] != BOND) { return 0; } else { if (x > 1) { if (world[x - 2][y] == BOND) return 0; } if (x < WIDTH - 2) { if (world[x + 2][y] == BOND) return 0; } } return 1; } // Return 1 if the spot is a glucose sugar, return 0 otherwise. int is_glucose (int x, int y) { if (x < 0 || x >= WIDTH || y < 0 || y >= HEIGHT) return 0; if (world[x][y] != SUGAR) { return 0; } else { if (x > 0) { if (world[x - 1][y] == BOND) return 0; } if (x < WIDTH - 1) { if (world[x + 1][y] == BOND) return 0; } } return 1; } int main (void) { #ifdef USE_SDL char wintitle[200]; #endif char filename[200]; int x; int y; int n; float chance; int tick = 0; int glucoseleft = 0; int glucoseright = 0; int bioseleft = 0; int bioseright = 0; int totalsugars = 0; int foo_bmptick; int seed; printf("\nSyn %s, a cellulose degradation model by Edinburgh iGEM\n\n", VERSION); seed = (int) time(NULL); srand(seed); // Clear the world... for (x = 0; x < WIDTH; x++) { for (y = 0; y < HEIGHT; y++) { world[x][y] = EMPTY; } } // Set up the cellulose... for (x = MARGIN + 1; x < HALFWIDTH - (MARGIN + 1); x += 2) { for (y = MARGIN; y < HEIGHT - MARGIN; y += 2) { // Left side... world[x][y] = BOND; world[x + 1][y] = SUGAR; world[x - 1][y] = SUGAR; // Right side... world[(WIDTH - 1) - x][y] = BOND; world[(WIDTH - 1) - (x + 1)][y] = SUGAR; world[(WIDTH - 1) - (x - 1)][y] = SUGAR; } } // How many sugars are there on each side, anyway? for (x = 0; x < HALFWIDTH; x++) { for (y = 0; y < HEIGHT; y++) { if (world[x][y] == SUGAR) totalsugars++; } } // Place the enzymes in random spots... // Cellulases on left side get assigned to 0 <= n <= ENZYMES - 1 for (n = 0; n < ENZYMES; n++) { endo[n].x = intrnd(HALFWIDTH - 3) + 1; // i.e. don't put things on the very edge. endo[n].y = intrnd(HEIGHT - 3) + 1; exo[n].x = intrnd(HALFWIDTH - 3) + 1; exo[n].y = intrnd(HEIGHT - 3) + 1; glu[n].x = intrnd(HALFWIDTH - 3) + 1; glu[n].y = intrnd(HEIGHT - 3) + 1; } // Cellulases on left side get assigned to ENZYMES <= n <= ENZYMES * 2 - 1 for (n = ENZYMES; n < ENZYMES * 2; n++) { endo[n].x = intrnd(HALFWIDTH - 3) + 1 + HALFWIDTH; endo[n].y = intrnd(HEIGHT - 3) + 1; exo[n].x = intrnd(HALFWIDTH - 3) + 1 + HALFWIDTH; exo[n].y = intrnd(HEIGHT - 3) + 1; glu[n].x = intrnd(HALFWIDTH - 3) + 1 + HALFWIDTH; glu[n].y = intrnd(HEIGHT - 3) + 1; } #ifdef USE_SDL // Start SDL and set up the screen... if (SDL_Init(SDL_INIT_VIDEO) < 0) { fprintf(stderr, "\nFatal error: Unable to initialize SDL!\n\n"); exit(EXIT_FAILURE); } virtue = SDL_SetVideoMode(WIDTH * SDL_MAG, HEIGHT * SDL_MAG, 0, SDL_ANYFORMAT); if (virtue == NULL) { fprintf(stderr, "\nFatal error: Unable to create virtual screen!\n\n"); exit(EXIT_FAILURE); } #endif printf(" Left: Synergy %s biose_inhib = %.3f\nRight: Synergy %s biose_inhib = %.3f\n\n", SYNERGY_LEFT ? "ON, " : "OFF,", BIOSE_INHIB_LEFT, SYNERGY_RIGHT ? "ON, " : "OFF,", BIOSE_INHIB_RIGHT); printf("Enzymes per type present per side: %d\n", ENZYMES); printf("Total sugars present per side: %d\n\n", totalsugars); printf("Iteration <--Glucose <--Cellobiose Glucose--> Cellobiose-->\n"); while (tick <= MAX_ITERATIONS || MAX_ITERATIONS == 0) { if (tick % REPORTTICK == 0) { printf("%9d%12d%15d%12d%15d\n", tick, glucoseleft, bioseleft, glucoseright, bioseright); } if (BMPTICK) { foo_bmptick = BMPTICK; // Use a variable for the test to avoid if (tick % foo_bmptick == 0 && tick) // div_by_zero warning at compile time. { sprintf(filename, "Syn-%05d.bmp", tick); drawbmp(filename, BMP_MAG); } } /* We now move the 3 enzyme types about and degrade appropriate bonds. Note that in the synergistic system, the order of the enzymes matters; in particular, if the Glu is 2 away from the Exo, it will automatically cut any cellobiose that the Exo creates on that side. This is sort of cheating, so it has been avoided by using the order Glu-Exo-Endo. */ // Do exoglucanase stuff... for (n = 0; n < ENZYMES * 2; n++) { randomwalk(&exo[n]); if (world[exo[n].x][exo[n].y] == BOND) { if (leftchain(exo[n].x, exo[n].y, 0) == 2 || rightchain(exo[n].x, exo[n].y, 0) == 2) { // Inhibition if cellobiose present in the right spots... chance = 1.0; if (is_cellobiose_bond(exo[n].x - 6, exo[n].y)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x + 6, exo[n].y)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x - 4, exo[n].y)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x + 4, exo[n].y)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x - 2, exo[n].y - 2)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x + 2, exo[n].y - 2)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x - 2, exo[n].y + 2)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x + 2, exo[n].y + 2)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x, exo[n].y - 2)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (is_cellobiose_bond(exo[n].x, exo[n].y + 2)) chance -= (n < ENZYMES) ? BIOSE_INHIB_LEFT : BIOSE_INHIB_RIGHT; if (rnd() < chance) { world[exo[n].x][exo[n].y] = EMPTY; // Exo can produce both glucose and cellobiose. // Keep track of these... if (n < ENZYMES) { if (is_glucose(exo[n].x - 1, exo[n].y)) glucoseleft++; if (is_glucose(exo[n].x + 1, exo[n].y)) glucoseleft++; if (is_cellobiose_bond(exo[n].x - 2, exo[n].y)) bioseleft++; if (is_cellobiose_bond(exo[n].x + 2, exo[n].y)) bioseleft++; } else { if (is_glucose(exo[n].x - 1, exo[n].y)) glucoseright++; if (is_glucose(exo[n].x + 1, exo[n].y)) glucoseright++; if (is_cellobiose_bond(exo[n].x - 2, exo[n].y)) bioseright++; if (is_cellobiose_bond(exo[n].x + 2, exo[n].y)) bioseright++; } } } } } // Do endoglucanase stuff... for (n = 0; n < ENZYMES * 2; n++) { if (n < ENZYMES) // Left side { if (SYNERGY_LEFT) { endo[n].x = exo[n].x + 1; endo[n].y = exo[n].y; } else { randomwalk(&endo[n]); } } else { // Right side if (SYNERGY_RIGHT) { endo[n].x = exo[n].x + 1; endo[n].y = exo[n].y; } else { randomwalk(&endo[n]); } } if (world[endo[n].x][endo[n].y] == BOND) { if (leftchain(endo[n].x, endo[n].y, 0) >= ENDO_REQ && rightchain(endo[n].x, endo[n].y, 0) >= ENDO_REQ) { world[endo[n].x][endo[n].y] = EMPTY; } } } // Do B-glucosidase stuff... for (n = 0; n < ENZYMES * 2; n++) { if (n < ENZYMES) // Left side { if (SYNERGY_LEFT) { glu[n].x = exo[n].x - 1; glu[n].y = exo[n].y; } else { randomwalk(&glu[n]); } } else { // Right side if (SYNERGY_RIGHT) { glu[n].x = exo[n].x - 1; glu[n].y = exo[n].y; } else { randomwalk(&glu[n]); } } if (world[glu[n].x][glu[n].y] == BOND) { if (is_cellobiose_bond(glu[n].x, glu[n].y)) { world[glu[n].x][glu[n].y] = EMPTY; // Glu destroys a cellobiose to make 2 glucose. // Keep track of these... if (n < ENZYMES) { bioseleft--; glucoseleft += 2; } else { bioseright--; glucoseright += 2; } } } } #ifdef USE_SDL // Draw the world... for (x = 0; x < WIDTH; x++) { for (y = 0; y < HEIGHT; y++) { if (world[x][y] == EMPTY) drawmag(virtue, x, y, SDL_MAG, 0, 0, 0); if (world[x][y] == SUGAR) drawmag(virtue, x, y, SDL_MAG, SUGAR_R, SUGAR_G, SUGAR_B); if (world[x][y] == BOND) drawmag(virtue, x, y, SDL_MAG, BOND_R, BOND_G, BOND_B); } } // Draw the enzymes... for (n = 0; n < ENZYMES * 2; n++) { drawmag(virtue, endo[n].x, endo[n].y, SDL_MAG, ENDO_R, ENDO_G, ENDO_B); drawmag(virtue, exo[n].x, exo[n].y, SDL_MAG, EXO_R, EXO_G, EXO_B); drawmag(virtue, glu[n].x, glu[n].y, SDL_MAG, GLU_R, GLU_G, GLU_B); } // Draw to the monitor... SDL_UpdateRect(virtue, 0, 0, 0, 0); manageevents(); SDL_Delay(DELAY); sprintf(wintitle, "i: %d | <--gluc: %d | <--biose: %d | gluc-->: %d | biose-->: %d", tick, glucoseleft, bioseleft, glucoseright, bioseright); SDL_WM_SetCaption(wintitle, wintitle); #endif tick++; } return 0; }