"

From 2014.igem.org

1. include<stdio.h>
2. include<math.h>
3. include <stdlib.h>
4. include "MT.h"

int main(int argc,char **argv){

   int seed;
   seed = atoi(argv[1]);
   init_genrand(seed);
   int n;
   int x;
   int space = 1000;
   double Asp[space+1];
   double Cd[space+1];
   double CdCd = 10.0;
   double range = 50.0;
   for(x=0;x<space+1;x++){
       Asp[x] = 0.0;
   }
   for(x=0;x<space+1;x++){
       Cd[x] = CdCd/(1.0+((x-1.0*space/2.0)/range)*((x-1.0*space/2.0)/range));
   }
   int producer = 100;
   int P_basyo[producer];
   int P_muki[producer];
   double P_Asp_past[producer];
   double P_Asp_now[producer];
   double P_Cd_past[producer];
   double P_Cd_now[producer];
   for(n=0;n<producer;n++){
       P_basyo[n] = (genrand_real3()-0.5)*space/3 + space/2;
       if(genrand_real3()<0.5){
           P_muki[n] = 1;
       }else{
           P_muki[n] = -1;
       }
       P_Asp_past[n] = 0.0;
       P_Asp_now[n] = 0.0;
       P_Cd_past[n] = P_basyo[n];
       P_Cd_now[n] = P_basyo[n];
   }
   int detecter = 100;
   int D_basyo[detecter];
   int D_muki[detecter];
   double D_Asp_past[detecter];
   double D_Asp_now[detecter];
   double D_Cd_past[detecter];
   double D_Cd_now[detecter];
   for(n=0;n<producer;n++){
       D_basyo[n] = (genrand_real3()-0.5)*space/3 + space/2;
       if(genrand_real3()<0.5){
           D_muki[n] = 1;
       }else{
           D_muki[n] = -1;
       }
       D_Asp_past[n] = 0.0;
       D_Asp_now[n] = 0.0;
       D_Cd_past[n] = Cd[D_basyo[n]];
       D_Cd_now[n] = Cd[D_basyo[n]];
   }
       double a = 0.2;
   double b = 0.004;
   double c = 0.005;
   double bias = 1.0;
   int timescale = 100000;
   FILE *fp;
   fp=fopen("result.txt","w");
   int t;
   for(t=0;t<timescale;t++){
       double new[space+1];
       //cyclic boundary condition
       Asp[0] = Asp[space-1];
       Asp[space] = Asp[1];
       for(n=0;n<producer;n++){
           Asp[P_basyo[n]] += b * Cd[P_basyo[n]];
       }
       for(x=1;x<space;x++){
           new[x] = a*Asp[x-1] + (1.0 - 2.0*a)*Asp[x] + a*Asp[x+1];
       }
       for(x=1;x<space;x++){
           Asp[x] = new[x];
       }
       for(x=1;x<space;x++){
           Asp[x] -= Asp[x]*c;
       }
       for(n=0;n<producer;n++){
           double kakuritu;
           //Chemotaxis for Cd and Asp
           kakuritu = bias*exp(P_Asp_now[n]-P_Asp_past[n])/(bias*exp(P_Asp_now[n]-P_Asp_past[n])+exp(P_Cd_now[n]-P_Cd_past[n]));
           //Randomly move
           //kakuritu = 0.5;
           //Chemotaxis for Asp (in  this case, we regard Cd as Asp)
           //kakuritu = exp(P_Cd_now[n]-P_Cd_past[n])/(1.0+exp(P_Cd_now[n]-P_Cd_past[n]));
           if(genrand_real3()<kakuritu){
               P_basyo[n] += P_muki[n];
               if(P_basyo[n] == 0){
                   P_basyo[n] = space-1;
               }else if(P_basyo[n] == space){
                   P_basyo[n] = 1;
               }
           }else{
               P_muki[n] = -P_muki[n];
           }
       }
       for(n=0;n<detecter;n++){
           double kakuritu;
           kakuritu = bias*exp(D_Asp_now[n]-D_Asp_past[n])/(bias*exp(D_Asp_now[n]-D_Asp_past[n])+exp(D_Cd_now[n]-D_Cd_past[n]));
           //kakuritu = 0.5;
           //kakuritu = exp(D_Cd_now[n]-D_Cd_past[n])/(1.0+exp(D_Cd_now[n]-D_Cd_past[n]));
           if(genrand_real3()<kakuritu){
               D_basyo[n] += D_muki[n];
               if(D_basyo[n] == 0){
                   D_basyo[n] = space-1;
               }else if(D_basyo[n] == space){
                   D_basyo[n] = 1;
               }
           }else{
               D_muki[n] = -D_muki[n];
           }
       }
       for(n=0;n<producer;n++){
           P_Asp_past[n] = P_Asp_now[n];
           P_Asp_now[n] = Asp[P_basyo[n]];
           P_Cd_past[n] = P_Cd_now[n];
           P_Cd_now[n] = Cd[P_basyo[n]];
       }
       for(n=0;n<detecter;n++){
           D_Asp_past[n] = D_Asp_now[n];
           D_Asp_now[n] = Asp[D_basyo[n]];
           D_Cd_past[n] = D_Cd_now[n];
           D_Cd_now[n] = Cd[D_basyo[n]];
       }
   }
   int number[space+1];
   for(x=0;x<space;x++){
       number[x] = 0;
   }
   for(n=0;n<producer;n++){
       number[P_basyo[n]] += 1;
   }
   for(n=0;n<detecter;n++){
       number[D_basyo[n]] += 1;
   }
   int tmp = 0;
   double tmp_Cd = 0.0;
   double tmp_asp = 0.0;
   int length = 10;
   /*
   for(x=0;x<space;x++){
       tmp_Cd += Cd[x];
       tmp_asp += Asp[x];
       tmp += number[x];
       if(x%length == length-1){
           fprintf(fp,"%d %d %f %f\n",x/length,tmp,tmp_Cd,tmp_asp);
           tmp = 0;
           tmp_Cd = 0.0;
           tmp_asp = 0.0;
       }
   }
   */
   for(x=0;x<space;x++){
       fprintf(fp,"%d %d %f %f\n",x,number[x],Cd[x],Asp[x]);
   }
   fclose(fp);
   FILE *gp;
   gp = popen("gnuplot" , "w" );
   fprintf(gp, "set terminal png\n" );
    fprintf(gp, "set output \"/dev/null\"\n" );
    fprintf(gp,"plot \"result.txt\" using 1:2 with points title 'Ecoli'\n" );
    fprintf(gp,"replot \"result.txt\" using 1:3 with points title 'Cd_concentration'\n" );
    fprintf(gp, "set output \"result.png\"\n" );
    fprintf(gp,"replot \"result.txt\" using 1:4 with points title 'Asp_concentration'\n" );
   /*
   fprintf(gp, "set output \"/dev/null\"\n" );
   fprintf(gp,"plot \"result.txt\" using 1:2 with points title 'Ecoli'\n" );
   fprintf(gp, "set output \"result.png\"\n" );
   fprintf(gp,"replot \"result.txt\" using 1:3 with points title 'Asp_concentration'\n" );//In this case, we regard Cd as Asp.
   */
   /*
    fprintf(gp, "set output \"result.png\"\n" );
    fprintf(gp,"plot \"result.txt\" using 1:2 with points title 'Ecoli'\n" );
    */
   /*
    fprintf(gp, "set output \"result.png\"\n" );
    fprintf(gp,"plot \"result.txt\" with lines title 'kazu'\n" );
    */
   fflush(gp);
   fprintf(gp, "exit\n" );
   pclose(gp);
   return 0;

}