#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>




/*********************READ THIS FIRST*******************/
//This program does both single graphs of mRNA or protein and multiple graphs of mRNA or protein. IT ALSO DOES MONTE CARLO SIMULATIONS.
//go to line 202 and set parameter1 and 2 and 3 to the three parameters you want to analyse, and their corresponding ranges.

//run... and have fun



























const int PoissonRandomNumber(const double lambda);




main()
{     
remove("stochastic.dat");     

srand(time(0));  //initialize random generator



 

 
double parameter1, parameter2, parameter3;
int erase,i,m,j,k=0, sign;////////
int good=0;/////////////////
int bad=0;//////////////////
float y[19];////////////////leave these as they are
double a;///////////////////
double d1,d2,d3,d4,d5, d6;//
double lambda, t=0;/////////
double xpercent;////////////


/**********************************CHANGE THESE VARIABLES AT WILL*****************************************/
int simulation_time=21600; //how long the simulation runs for before breaking and "failing" (auto breaks when lysis happens if "LYSIS=1")
int IPTG_time=3600;   //when IPTG is added
int HSL_time=10800;  //when HSL is added
double   tau=1; // the time step (seconds) (this is an important parameter but can be changed between around 0.01 -> 2 seconds)  
int printing;   //set as 1 to print the levels of proteins and mrnas to "stochastic.dat".... or 0 not to print anything
int printwhen=1;   //sets how often the program prints the protein and mrna levels (can be anything from 1 to several thousand)
int LYSIS; //shall prog end when lysis occurs? (1 for yes, 0 for no)
int runs;


//////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
printf ("Press 0 for Monte Carlo \nPress 1 for plotting graphs\n");
scanf ("%d",&printing);

if (printing==0)
{
printing=0;
LYSIS=1;
printf("delete old parameters? (0 for no,1 for yes)\n");
scanf ("%d",&erase);
      if(erase==1)
      {
      remove("montecarlo_parameters.dat");  
      }
      
printf("How many runs?\n");
scanf("%d",&runs);
}

else
{
printing=1;
printf("How many runs would you like?\n");
scanf("%d",&runs);
printf("How often do you want to print?\nType 1 to print every %lf seconds (this is limited by the timestep, tau)\nPress a large integer for quicker plotting (type 10 to print every 10*%lf seconds)\n",tau, tau);
scanf("%d",&printwhen);
printf("Finish simulations when lysis occurs? (0 for no, 1 for yes)\n");
scanf("%d",&LYSIS);
}
////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////

  
FILE *outfile1;
FILE *outfile2;
outfile1=fopen("stochastic.dat","w");   //for printing the protein and mrna levels at each timestep
outfile2=fopen("montecarlo_parameters.dat","a"); //for printing the chosen monte carlo parameters


for (m=1;m<=runs;m++)   //montecarlo loop opens   //set from 1 to 1 to plot one graph (do one simulation) set it from 1 to several thousand to run a montecarlo sim
{
    //////////////////////////////////////////////////////////////////////////////////////////////
    
int lysis_on=0;///////////
int lysis_count=0;//////// leave as they are
t=0;//////////////////////
   
  
  /***********************************************************/
      /********************** PARAMETERS *************************/
      /***********************************************************/
    
    double HSL_outside               =10000; //sets the outside levels of HSL (when it is added)
    double IPTG_outside              =10000; //sets the outside levels of IPTG (when it is added)
    
    
    
    
    double mu                        = 0.45;    //rate at which LuxI produces HSL
    double theta_HSL                 = 0.24;    // Rate at which HSL diffuses into or out of cell
    double theta_IPTG                = 0.0137;   //Rate at whih IPTG diffuses into or out of cell
    
    double rho_production            = 10;  // plasmid copy numbers
    double rho_latch                 = 10;
    double rho_lysis                 = 10;
    double rho_antilysis             = 10;
    double rho_qs                    = 10;
    double rho_luxi                  = 10;
    
    
    double beta_production_max       = 0.44;   //promoter strengths in pops (molecules produced per second)
    double beta_production_min       = 0.013;
    double beta_production_max2      = 0.002;
    double beta_production_min2      = 0.00015; 
    double beta_latch                = 0.28;
    double beta_lysis                = 0.0426;
    double beta_antilysis            = 0.0066;
    double beta_luxr                 = 0.01;
    
    
                                    // min + (max - min) //sets the montecarlo range for these three parameters
    double K_laci                    = 700;  
    double K_ci                      = 7000; 
    double K_tetr                    = 7000; 
    double K_P                       = 700;
    double K_IPTG_laci               = 400;
    
    double n_laci                    = 2;
    double n_ci                      = 2;
    double n_tetr                    = 3; //hill coefficents
    double n_P                       = 2; 
    
    double gamma_protein             = 0.1; //rate that mRNA produces protein molecules
    
    
    double k_on_P                    = 0.00001; //rate for HSL and luxR producing complex P
    double k_off_P                   = 0.00333; // off rate for the above
    
    
    
    //degradation constants
    double alpha_mrna                = 0.00288;                    // ~2 minutes
    double alpha_X                   = 0.0000080225;               // 24 hours
    double alpha_Y                   = 0.0000080225;               // 24 hours
    double alpha_ci                  = 0.002888; // taged          // 4 minutes
    double alpha_laci                = 0.00115;                    // 10 minutes
    double alpha_tetr                = 0.00115;                    // 10 minutes
    double alpha_holin               = 0.0002 ;                    // 1 hour
    double alpha_antiholin           = 0.0002 ;                    // 1 hour
    double alpha_endolysin           = 0.0002;                     // 1 hour
    double alpha_luxi                = 0.00288811; //taged         //4 minutes
    double alpha_luxr                = 0.002888;   //taged         //4 minutes
    
    double alpha_hsl                 = 0.00017;                    // ~1 hour
    
    double alpha_p                   = 0.0002;                     // 1 hour
    
   if (printing==0)
   { 
    /***************MONTE CARLO PARAMETERS************************/
    
                   // min + (max - min) //sets the montecarlo range for  three parameters
     parameter1 =100 + (2000-100)*rand()/((RAND_MAX));             
     parameter2 =700 + (7000-700)*rand()/((RAND_MAX));
     parameter3 =700 + (7000-700)*rand()/((RAND_MAX));
    
      K_laci=parameter1; 
      K_ci=parameter2;///////////////set what the parameters are
      K_tetr=parameter3;
}
    
    
    
//initial conditions 

y[0]=/*PRODUCTION*/0;
y[1]=/*LATCH     */(beta_latch/(alpha_mrna+(log(2)/1800)))*rho_latch; 
y[2]=/*LYSIS     */0;
y[3]=/*ANTILYSIS */(beta_antilysis/(alpha_mrna+(log(2)/1800)))*rho_lysis; 
y[4]=/*QS        */(beta_production_min2/(alpha_mrna+(log(2)/1800)))*rho_qs;                  
y[5]=/*X         */0;
y[6]=/*Y         */0;
y[7]=/*cI        */0;
y[8]=/*lacI      */(gamma_protein/(alpha_laci+(log(2)/1800)))*y[1];
y[9]=/*tetR      */(gamma_protein/(alpha_tetr+(log(2)/1800)))*y[1];
y[10]=/*Holin    */0;
y[11]=/*Endolysin*/0;
y[12]=/*Antiholin*/(gamma_protein/(alpha_antiholin+(log(2)/1800)))*y[3];
y[13]=/*luxI     */(gamma_protein/(alpha_luxi+(log(2)/1800)))*y[4];
y[14]=/*LuxR     */0;//(gamma_protein/(alpha_luxr+(log(2)/1800)))*y[4];
y[15]=/*HSL      */0;
y[16]=/*HSL/LuXR */0;
y[17]=/*IPGT     */0;
y[18]=/*IPGToutside*/0;
y[19]=/*HSLoutside*/0;
y[20]=/*mRNA for luxR*/(beta_luxr/(alpha_mrna+(log(2)/1800)))*rho_qs; 



double antiholinthreshold=0;  //keep as is
double lysis_time=(1/tau)*300;  //300 is the number of seconds at which holin must be above antiholin + threshold (set to 1000 molecules) before the cell lyses


 

 
   
  
while (t<39600) // 11 hours (make arbitarily high)
{
    
   t=t+tau; //t is the count of real time (seconds)
     k++; // a count to determine when printing occurs
    
    
//INPUT SIGNALS
    
// IPTG outside  (after 1 hour add IPTG to the surrounding solution)
if (t>IPTG_time) // 1 hours
{
y[18]=IPTG_outside;//  y[18]++;
}
  
 
  
//HSL outside (afrer 3 hours add HSL to the surrounding solution)
if(t>HSL_time) //3 hours
{
y[19]=HSL_outside    ;  //    y[19]++; 
}
             
              
/***********************************************************/
/******************** MRNA EQUATIONS ***********************/
/***********************************************************/
      
      
     
//  mRNA_production - first term of equation
lambda=rho_production*(1/(1+pow((y[8]/(K_laci*(1+pow(y[17]/K_IPTG_laci,n_laci)))),n_laci)))*((beta_production_max*pow(y[16],n_P)/(pow(y[16],n_P)+pow(K_P,n_P)))+beta_production_min)*tau;
d1=PoissonRandomNumber(lambda);

// mRNA_production - second term of equation
lambda=(alpha_mrna+(log(2)/1800))*y[0]*tau;
d2=PoissonRandomNumber(lambda);
    
//mRNA_production - equation
y[0]=y[0]+d1-d2;
     
      
      
//mRNA_latch - first term of the equation
lambda=rho_latch*beta_latch*(1/(1+pow((y[7]/K_ci),n_ci)))*tau;
d1=PoissonRandomNumber(lambda);
      
//mRNA_latch - second term of the euqation
lambda=(alpha_mrna+(log(2)/1800))*y[1]*tau;
d2=PoissonRandomNumber(lambda);
      
//mRNA_latch - equation
y[1]=y[1]+d1-d2;


      
//mRNA_lysis - first term of the equation
lambda=rho_lysis*beta_lysis*(1/(1+pow((y[9]/K_tetr),n_tetr)))*tau;
d1=PoissonRandomNumber(lambda);
      
//mRNA_lysis - second term of the euqation
lambda=(alpha_mrna+(log(2)/1800))*y[2]*tau;
d2=PoissonRandomNumber(lambda);
      
//mRNA_lysis - equation
y[2]=y[2]+d1-d2;
      
      
      
//mRNA_antilysis - first term of the equation
lambda=rho_antilysis*beta_antilysis*tau;
d1=PoissonRandomNumber(lambda);
      
//mRNA_antilysis - second term of the equation
lambda=(alpha_mrna+(log(2)/1800))*y[3]*tau;
d2=PoissonRandomNumber(lambda);
      
//mRNA_antilysis - equation
y[3]=y[3]+d1-d2;
      
      
//  mRNA LuxI first term of the equation
lambda=tau*rho_qs*beta_production_max2*(pow(y[16],n_P)/(pow(y[16],n_P)+pow(K_P,n_P)));
d1=PoissonRandomNumber(lambda);
      
// mRNA LuxI second term of the equation      
lambda=rho_qs*beta_production_min2*tau;
d2=PoissonRandomNumber(lambda);
      
//mRNA_luxi - third term of the equation
lambda=(alpha_mrna+(log(2)/1800))*y[4]*tau;
d3=PoissonRandomNumber(lambda);
      
//  mRNA_luxi - equation
y[4]=y[4]+d1+d2-d3;
      
      
// mRNA luxR equation first term
lambda=rho_qs*beta_luxr*tau;
d1=PoissonRandomNumber(lambda);
 

//mRNA LuxI equation second term      
lambda=(alpha_mrna+(log(2)/1800))*y[20]*tau;
d2=PoissonRandomNumber(lambda);


//mRNA LuxI - equation 
y[20]=y[20]+d1-d2;
      /***********************************************************/
      /****************** PROTEIN EQUATIONS **********************/
      /***********************************************************/
      
         
      //[X] - first term of the equation
      lambda=gamma_protein*y[0]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[X] - second term of the equation
      lambda=(alpha_X+(log(2)/1800))*y[5]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[X] - equation
      y[5]=y[5]+d1-d2;
      
      
    
      //[Y] - first term of the equation
      lambda=gamma_protein*y[0]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[Y] - second term of the equation
      lambda=(alpha_X+(log(2)/1800))*y[6]*tau; //Y is not important if we have X!
      d2=PoissonRandomNumber(lambda);
      
      //[Y] - equation
      y[6]=y[6]+d1-d2;
      
      
   
       //[cI] - first term of the equation
      lambda=gamma_protein*y[0]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[cI] - second term of the equation
      lambda=(alpha_ci+(log(2)/1800))*y[7]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[cI] - equation
      y[7]=y[7]+d1-d2;
      
      
      
       //[LacI] - first term of the equation
      lambda=gamma_protein*y[1]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[LacI] - second term of the equation
      lambda=(alpha_laci+(log(2)/1800))*y[8]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[LacI] - equation
      y[8]=y[8]+d1-d2;
      
      
      
      //[TetR] - first term of the equation
      lambda=gamma_protein*y[1]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[TetR] - second term of the equation
      lambda=(alpha_tetr+(log(2)/1800))*y[9]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[TetR] - equation
      y[9]=y[9]+d1-d2;
      
      
      
      //[Holin] - first term of the equation
      lambda=gamma_protein*y[2]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[Holin] - second term of the equation
      lambda=(alpha_holin+(log(2)/1800))*y[10]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[Holin] - equation
      y[10]=y[10]+d1-d2;
      
      
      
       //[Endolysin] - first term of the equation
      lambda=gamma_protein*y[2]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[Endolysin] - second term of the equation
      lambda=(alpha_endolysin+(log(2)/1800))*y[11]*tau; //endolysin is not important!
      d2=PoissonRandomNumber(lambda);
      
      //[Endolysin] - equation
      y[11]=y[11]+d1-d2;
      
      
      
      //[Antiholin] - first term of the equation
      lambda=gamma_protein*y[3]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[Antiholin] - second term of the equation
      lambda=(alpha_antiholin+(log(2)/1800))*y[12]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[Antiholin] - equation
      y[12]=y[12]+d1-d2;
      
      
      
      //[LuxI] - first term of the equation
      lambda=gamma_protein*y[4]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[LuxI] - second term of the equation
      lambda=(alpha_luxi+(log(2)/1800))*y[13]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[LuxI] - equation
      y[13]=y[13]+d1-d2;
      
      
      
      //[LuxR] - first term of the equation
      lambda=gamma_protein*y[20]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[LuxR] - second term of the equation
      lambda=k_off_P*y[16]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[LuxR] - third term of the equation
      lambda=k_on_P*(y[15]*y[14])*tau;
      d3=PoissonRandomNumber(lambda);
      
      //[LuxR] - fourth term of the equation
      lambda=(alpha_luxr+(log(2)/1800))*y[14]*tau;
      d4=PoissonRandomNumber(lambda);
      
      //[LuxR] - equation
      y[14]=y[14]+d1+d2-d3-d4;
      
      /***********************************************************/
      /****************** COMPLEX EQUATIONS **********************/
      /***********************************************************/
      
      //[HSL] - first term of the equation
      lambda=mu*y[13]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[HSL] - second term of the equation
      lambda=(k_off_P+alpha_luxr)*y[16]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[HSL] - third term of the equation
      lambda=k_on_P*y[15]*y[14]*tau;
      d3=PoissonRandomNumber(lambda);
      
      //[HSL] - fourth term of the equation (adjusted becuase the poisson random generator does not output negative integers
      lambda=theta_HSL*(y[19]-y[15])*tau;
      
      if (lambda<0)
      {
      sign=-1;
      lambda=lambda*(-1);
      d4=PoissonRandomNumber(lambda);
      }
      else
      { 
      sign= 1;
      d4=PoissonRandomNumber(lambda);
      }
      
      
     
      //[HSL] - fifth term of the equation
      lambda=(alpha_hsl+(log(2)/1800))*y[15]*tau;
      d5=PoissonRandomNumber(lambda);
      
      //[HSL] - equation
      y[15]=y[15]+d1+d2-d3+(sign*d4)-d5; 
      
      
      
      
      //[P] - first term of the equation
      lambda=k_on_P*y[15]*y[14]*tau;
      d1=PoissonRandomNumber(lambda);
      
      //[P] - second term of the equation
      lambda=k_off_P*y[16]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[P] - third term of the equation
      lambda=(alpha_p+(log(2)/1800))*y[16]*tau;
      d3=PoissonRandomNumber(lambda);
     
      //[P] - equation
      y[16]=y[16]+d1-d2-d3;
      
      
      
      
      
      //[IPTG inside cell] - first term of the equation
      lambda=theta_IPTG*y[18]*tau;
      d1=PoissonRandomNumber(lambda);
    
    //[IPTG inside cell] - second term of the equation
      lambda=theta_IPTG*y[17]*tau;
      d2=PoissonRandomNumber(lambda);
      
      //[IPTG inside cell] - equation
      y[17]=y[17]+d1-d2;
      
 
    
 
  
  
 if (k==printwhen &&printing==1) //prints the protein and mrna concentrations
    {            
   fprintf(outfile1,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",t,y[0],y[1],y[2],y[3],y[4],y[5],y[6],y[7],y[8],y[9],y[10],y[11],y[12],y[13],y[14],y[15],y[16],y[17],y[18],y[19],y[20]);    
    k=0;
    }
    
    
    antiholinthreshold=y[12]+1000; //sets the lysis command so that lysis begins when holin is 1000 molecules above antiholin
    
    
    
if  (y[10] > antiholinthreshold &&LYSIS==1) // counts up
{
    lysis_count++;
    lysis_on=1;
}
 
if (lysis_on==1 && y[10]< antiholinthreshold) //counts down
{
    lysis_count=lysis_count-1;
}

if (lysis_count>lysis_time && t<HSL_time)   //IF LYSIS HAPPENS AND BOTH SIGNALS ARE NOT ON.... THEN WE HAVE FAILED!
    { 
                     bad++;
                     printf("Good=%d\tBad=%d\t%lf\t\tLysis too soon!\n", good, bad, xpercent);
                     fprintf(outfile2,"%lf\t %lf\t %lf\t 0\n",parameter1,parameter2,parameter3);
                         
                         break;
                           }



if (lysis_count>lysis_time && t>HSL_time)   //IF LYSIS HAPPENS AND BOTH SIGNALS ARE ON THEN WE HAVE SUCCEEDED
    {   
         xpercent=(y[5] /((gamma_protein/(alpha_X+(log(2)/1800)))*((beta_production_max/(alpha_mrna+(log(2)/1800))*rho_production))))*100;    
          
                     good++;
                     printf("Good=%d\tBad=%d\t%lf\t\tSuccess!\n", good, bad, xpercent);
                     fprintf(outfile2,"%lf\t %lf\t %lf\t %lf\n",parameter1,parameter2,parameter3,xpercent);                                             
                     break;  
    }
      
      
    if (t>simulation_time) //  6 hours  //IF LYSIS DOES NOT HAPPEN IN 6 HOURS THEN WE HAVE FAILED
      {
             bad++;
             xpercent=(y[5] /((gamma_protein/(alpha_X+(log(2)/1800)))*((beta_production_max/(alpha_mrna+(log(2)/1800))*rho_production))))*100;      
             printf ("Good=%d\tBad=%d\t%lf\t\tFailure to lyse!\n",good,bad,xpercent);
               fprintf(outfile2,"%lf\t %lf\t %lf\t 0\n",parameter1,parameter2,parameter3);
             break  ;
    }  
               
} //while loop closes
 fprintf(outfile1,"\n");
} //montecarlo closes

 fclose(outfile1);
 
  return 0;
} //main closes




const int PoissonRandomNumber(const double lambda)
{
  int k=0;                          //Counter
  const int max_k = 3*lambda;           //k upper limit
  double p = 1.*rand()/((RAND_MAX)); //uniform random number
  double P = exp(-lambda);          //probability
  double sum=P;                     //cumulant
  if (sum>=p) return 0;             //done allready
  for (k=1; k<max_k; ++k) {         //Loop over all k:s
    P*=lambda/(double)k;            //Calc next prob
    sum+=P;                         //Increase cumulant
    if (sum>=p) break;              //Leave loop
  }

  return k;                         //return random number
}