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

/******************READ THIS FIRST***********************/
// this is a variation on the Monte Carlo simulation program. it varies 2 parameters over a user specified range, with user specified steps.
//simply go to line 110 and 132 and change the ranges and names of parameters you want to study.



const int PoissonRandomNumber(const double lambda);

main()
{
      
remove("parameteranalysis.dat");
FILE *outfile1;
outfile1=fopen("parameteranalysis.dat","a");      
srand(time(0));  //initialize random generator

double a;
int i,m,j,k=0, sign;
float y[20];
double d1,d2,d3,d4,d5, d6;
double lambda, t=0;
double   tau=1; 
double xpercent;
double parameter1, parameter2, parameter1low, parameter1high, parameter2low, parameter2high, parameter1step, parameter2step;
int simulation_time=21600;
int HSL_time=10800;
int IPTG_time=3600;
  


 /***********************************************************/
 /********************** 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
    
    
    
////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////SET THESE PARAMETERS///////////////////////////////////////////////////////////////////               


parameter1=parameter1low=0;
          parameter1high=7000;
          parameter1step=200;

parameter2=parameter2low=0;
          parameter2high=7000;
          parameter2step=200;

               
while(parameter1 >=parameter1low && parameter1<parameter1high)///////
{         ///////////////////////////////////////////////////////////
parameter2=parameter2low; //resets second parameter back to zero/////
parameter1=parameter1+parameter1step;//////////////////////////////// leave this part alone
                                       //////////////////////////////
while (parameter2>=parameter2low && parameter2<parameter2high)///////
{         ///////////////////////////////////////////////////////////
parameter2=parameter2+parameter2step;     ///////////////////////////




//////////////////////////////////SET HERE ASWELL///////////////////////////////////////////////////////////
K_tetr=parameter1;
K_laci=parameter2;
////////////////////////////////////////////////////////////////////////////////////////////////////////////
     
      
      
int lysis_on=0;
int lysis_count=0;
t=0;
   
  

//input 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)

    
    
//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]++; 
}


                  
//define timesteps, t, tau



    
/***********************************************************/
/******************** 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;
      
 
    

  
    antiholinthreshold=y[12]+1000;
if  (y[10] > antiholinthreshold)
{
    lysis_count++;
    lysis_on=1;
}
 
if (lysis_on==1 && y[10]< antiholinthreshold)
{
    lysis_count=lysis_count-1;
}

if (lysis_count>lysis_time && t<HSL_time)  
    { 
    fprintf(outfile1,"%lf\t %lf\t %lf\t 0\n",parameter1,parameter2,t); //LYSIS TOO SOON
    printf (     "%lf\t %lf\t %lf\t 0\n",parameter1,parameter2,t);                                
    break;
    }



if (lysis_count>lysis_time && t>HSL_time)  
    {   
         xpercent=(y[5] /((gamma_protein/(alpha_X+(log(2)/1800)))*((beta_production_max/(alpha_mrna+(log(2)/1800))*rho_production))))*100;    
        
         fprintf(outfile1,"%lf\t %lf\t %lf\t %lf\n" ,parameter1,parameter2, t, xpercent); //Lysis on time
         printf (     "%lf\t %lf\t %lf\t %lf\n" ,parameter1,parameter2, t, xpercent);
           break;
     } 
      
    if (t>simulation_time) //  6 hours
      {
           
             fprintf(outfile1,"%lf\t %lf\t %lf\t 0\n",parameter1,parameter2,t); // no lysis at all
             printf (     "%lf\t %lf\t %lf\t 0\n",parameter1,parameter2,t);   
             break  ;
      }  
      
      
                       
                                        
      
               
} //while loop closes


                                              

} //close inside while loop 
                                               
                                                 fprintf(outfile1,"\n" ); //always on
                                                 printf (     "\n" );
}// close outside while loop


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




const int PoissonRandomNumber(const double lambda)
{
  int k=0;                          //Counter
  const int max_k = 2*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
}