/*

3-state model for tRNA charging cycles for all 20 amino acids including synthetase sequestration:

	tRNA is uncharged and empty (state 0) --> synthetase is bound to tRNA (state 1) --> tRNA is charged with amino acid (state 2)
	transitions:

	0 --> 1 binding: depends on the enzymatic concentration E of the synthetase, the doxycycling concentration and the number of uncharged tRNAs of this type
	1 --> 2	charging: depends on the number of bound tRNAs, rate constant is affiliated with kcat, the chargin process does NOT depend on doxycycline (Decoupling)
	2 --> 0 usage: depends on the number of charged tRNAs (corresponds to hopping rate of ribosome in the GTM) and via a feedback factor on doxycycline which reflects that doxycycline reduces the growth rate, i.e., the ribosomal current which is connected to the usage rate for ALL tRNAs (Coupling)

This model does not distinguish between the different tRNAs which are served by the same amino acid and therefore, competition is ruled out here.

*/

//####################################################################################

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <math.h>
#include "mtt.c"	// from  http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html

//####################################################################################
//CONSTANTS
//####################################################################################

#define STATEMAX 3	//number of states per cycle
#define TRANSMAX 4	//number of transitions
#define CYCLEMAX 20	//number of cycles, here: only GLN
#define GLN 5		// index of glutamin 


#define STEPS 	10000000000//Gillespie steps
//#define STEPS 	10000000//Gillespie steps, too small!


#define SCALE 	1	//scales all synthesases to system size; cellsize=1
#define CHARGE 	10490	//scales all tRNA copy numbers to sysytem size;  cellsize=10490
#define HOPP 	0.1	//rate constant for "uncharging" (=usage) rate;  cellsize=0.1
#define allcharged 1.0	// 1.0: all tRNAs initially charged; (1.0-allcharged)=uncharged;

#define k0 	43.0	//degredation rate constant
#define k1phys 	0.007	//physiological binding rate for the whole cell (=1.2*10^8 1/(sM))
#define CHARGEFACTOR	0.2	//scales the charging rate so that we get Q=80% inital charging for no doxy
#define construct	4	//Gln4 strain construct has ca. 4 times more Gln synthetase cp to physiological level


#define kcat_factor 1.0 // tuning the kcat-values: 1.0=physiological value; 0.1=small, 10=large



//####################################################################################
//####################################################################################
//####################################################################################
//MAIN
//####################################################################################
//####################################################################################
//####################################################################################

int main(int argc, char *argv[])
{
time_t start,end;
time(&start);



long int t0=0;//t0 counts number of iterations (i.e. reactions)
double r1,r2,tau,a0,runtime,randReaction;// random numbers, dwell time, initial sum of propensities, time, random reaction choosing value;
double suma=0.0;// sum of propensities
int i,j,cycle,state; //counter
int deg=0; //identifies if transition is usage (=0) or degredation (deg=1);
int correctcycle=0;
int correctstate=0;

long int Transient;

double real_t;// realtime t-timetransient rounded
double time_Transient;
int THRESHOLD;

int syn[CYCLEMAX],Etotal[CYCLEMAX],E[CYCLEMAX],copies[CYCLEMAX];//synthetase cell, synthetase small sysytem, GCN of tRNAs
int t[CYCLEMAX],c[CYCLEMAX],b[CYCLEMAX],u[CYCLEMAX],dummy_aa, dummy_tRNA, dummy_codon; //number of tRNAs: total, charged, bound, uncharged, dummies;
double kcat[CYCLEMAX],doxy[CYCLEMAX],doxyCH[CYCLEMAX];//kcat, doxycycling induced downscaling factor for E in binding rate, doxycycline induced downscaling factor for E in charging rate
double BIND[CYCLEMAX],USE[CYCLEMAX];//binding rate constant, usage rate constant
double w[CYCLEMAX][TRANSMAX], trans[CYCLEMAX][TRANSMAX], sum_trans[CYCLEMAX],accum_trans[CYCLEMAX];//propensities
double mean_u[CYCLEMAX],mean_b[CYCLEMAX],mean_c[CYCLEMAX]; //mean values (over time): uncharged, bound, charged levels of tRNAs
double k1[CYCLEMAX],k2[CYCLEMAX],k3;//rate constants

double sum_w=0.0;
double sum_before=0.0;


FILE *f1p, *f2p, *f3p, *f4p;
char string1[100],string2[100],string3[100],string4[100];


sprintf(string3,"mkdir OUTPUT_SynSeq");
  system(string3);


//####################################################################################
//LOOP 
//####################################################################################

double d=1.0;//doxyscale
double FB; //feedback factor, connects the doxy reduction on the growth rate to the usage rate



Transient=STEPS;




int arg_d;
for(arg_d=0;arg_d<27;arg_d++){ //doxy-loop for non-cluster run

if(arg_d==0)		d=	0.01	;
else if (arg_d==1)	d=	0.012	;
else if (arg_d==2)	d=	0.014	;
else if (arg_d==3)	d=	0.016	;
else if (arg_d==4)	d=	0.018	;
else if (arg_d==5)	d=	0.02	;
else if (arg_d==6)	d=	0.03	;
else if (arg_d==7)	d=	0.04	;
else if (arg_d==8)	d=	0.05	;
else if (arg_d==9)	d=	0.06	;
else if (arg_d==10)	d=	0.07	;
else if (arg_d==11)	d=	0.08	;
else if (arg_d==12)	d=	0.09	;
else if (arg_d==13)	d=	0.1	;
else if (arg_d==14)	d=	0.12	;
else if (arg_d==15)	d=	0.14	;
else if (arg_d==16)	d=	0.16	;
else if (arg_d==17)	d=	0.18	;
else if (arg_d==18)	d=	0.2	;
else if (arg_d==19)	d=	0.3	;
else if (arg_d==20)	d=	0.4	;
else if (arg_d==21)	d=	0.5	;
else if (arg_d==22)	d=	0.6	;
else if (arg_d==23)	d=	0.7	;
else if (arg_d==24)	d=	0.8	;
else if (arg_d==25)	d=	0.9	;
else if (arg_d==26)	d=	1.0	;
else printf("ERROR with d");


//####################################################################################
//INITIATION
//####################################################################################
for(j=0;j<CYCLEMAX;j++)	//initiation
    {
	E[j]=0;
	Etotal[j]=0;
	t[j]=0;
	c[j]=0;
	b[j]=0;
	u[j]=0;


	kcat[j]=0.0;
	BIND[j]=0.0;
	USE[j]=0.0;

	k1[j]=0.0;
	k2[j]=0.0;
	k3=0.0;	

	doxy[j]=0.0;
	doxyCH[j]=1.0;

	mean_c[j]=0.0;
	mean_u[j]=0.0;
	mean_b[j]=0.0;
		
        for(i=0;i<TRANSMAX;i++) w[j][i]=0.0;
    }



for(j=0;j<CYCLEMAX;j++)	//initiation
     {
	mean_c[j]=0.0;
	mean_u[j]=0.0;
	mean_b[j]=0.0;


        for(i=0;i<TRANSMAX;i++) {w[j][i]=0.0;trans[j][i]=0.0;}
	sum_trans[j]=0.0;
	accum_trans[j]=0.0;
     } 




	time_t seed; 	// for initialisation of the random generator mtt.c 
	seed = 260;
	init_genrand((unsigned long)seed);
	
	t0=0;
	correctcycle=0;
	correctstate=0;



	suma=0.0;a0=0.0; 
	sum_w=0.0;
	runtime=0.0; tau=0.0;
	time_Transient=0.0;
	real_t=0.0;



//####################################################################################
//FEEDBACK
//####################################################################################
//feedback factor, connects the doxy impact on the growth rate to the usage rate: i.e., if d gets small, the hopping rate gets smaller.

if(d<0.2)	FB=5.0*d;
else 		FB=1.0;


//####################################################################################
//PARAMETERs
//####################################################################################


	sprintf(string1,"VMAX.dat"); 
	f1p=fopen(string1,"r");

  for(j=0; j<CYCLEMAX; j++)
    {
      fscanf(f1p,"%d\t%lf\t%d\t%lf\t%d\t%lf\n",&syn[j],&kcat[j], &dummy_tRNA,&BIND[j],&copies[j],&USE[j]);
	


	Etotal[j]=ceil(syn[j]/SCALE);		//total number of synthetase molecules if type j
	t[j]=copies[j]*CHARGE;			//total number of tRNAs of type j
	c[j]=allcharged*floor(t[j]);		//start config. tRNAs charged
	u[j]=floor((1.0-allcharged)*t[j]);	//start config. tRNAs uncharged
	b[j]=t[j]-c[j]-u[j];			//bound tRNAs

	k1[j]=BIND[j];				//binding rate constant
	k2[j]=CHARGEFACTOR*kcat[j];		//charging rate constant
	k3=FB*HOPP;				//usage rate constant
	
    }

  fclose(f1p);
    

 for(j=0; j<CYCLEMAX; j++) kcat[j]=kcat_factor*kcat[j];	//for the case of artifical kcat-tuning
    

//####################################################################################
//PROPENSITIES 
//####################################################################################

//synthetase
for(j=0;j<CYCLEMAX;j++)
    {
	if(j==GLN)	{
			doxy[j]=d;	
			Etotal[j]=ceil(d*construct*Etotal[j]);	//glutamin synthetase concentration is doxycycline sensitiv
			}	
	else		doxy[j]=1.0;	
	
	if(Etotal[j]-b[j]>=0)	E[j]=Etotal[j]-b[j];	//number of free synthetase 
	else 			{printf("Etotal=%d<0!!!\n",Etotal[j]); E[j]=0;}
	
     }

//transition rates
 for(j=0; j<CYCLEMAX; j++)	
    {

     	
	w[j][0]=k1[j]*E[j]*u[j];	//binding transition from state 0 (empty tRNA) -> 1 (Synthetase bound)		

	w[j][1]=k2[j]*b[j];		//charging transition from state 1 (Synthetase bound) -> 2 (charged tRNA) 

	w[j][2]=k3*c[j];		//usage (=uncharging) transition from state 2 (charged tRNA) -> 0 (empty tRNA)

	w[j][3]=k0*b[j]; 		//dissociation  of the bound complex
		

        for(i=0;i<TRANSMAX;i++) sum_w=sum_w+w[j][i];
   }

//transition probabilities
for(j=0;j<CYCLEMAX;j++)
     {    	

	for(i=0;i<TRANSMAX;i++) 
		{
		 trans[j][i]=w[j][i]/sum_w;	//transition probability within cycle j from state i
		sum_trans[j]=sum_trans[j]+trans[j][i];	//sum of the transition probabilities per cycle
		}

	i=0;
	while(i<=j)
	{   
	accum_trans[j]=accum_trans[j]+sum_trans[i]; i++;	//intervall stacking
	}

     }



//#########################################################################################      
//print PARAMETERS
//#########################################################################################

sprintf(string2,"OUTPUT_SynSeq/Parameters.dat");
f2p=fopen(string2,"a");
	fprintf(f2p,"---------------------------\n");
	fprintf(f2p,"aminoacid, k2=a*kcat, k3=d*hopp*USE, k0,k1\n");
	fprintf(f2p,"T_c,T_e,T_b,E[j] \n");
	fprintf(f2p,"charge, usage, degredation, bind. rate\n");
	fprintf(f2p,"---------------------------\n");
   fclose(f2p);
for(j=0; j<CYCLEMAX; j++)
   	{

		f2p=fopen(string2,"a");
		fprintf(f2p,"%d:\t%.2f \t\t %.2f \t\t %.2f \t\t %.3f\n",j,k2[j],k3,k0,k1[j]);
		fprintf(f2p,"   \t%d \t\t %d \t\t %d \t\t %d\n", c[j],u[j],b[j],	E[j]);
		fprintf(f2p,"   \t%.1e \t %.1e \t %.1e \t %.1e\n",	w[j][1],w[j][2],w[j][3],w[j][0]);
		fprintf(f2p,"\n");
		fclose(f2p);
	}


   f2p=fopen(string2,"a");
	fprintf(f2p,"---------------------------\n");
	fprintf(f2p,"Scaling factor=%d\n",SCALE);
	fprintf(f2p,"a=%f\n",CHARGEFACTOR);
	fprintf(f2p,"Gln4 protein ratio=%f\n",d);
	fprintf(f2p,"Gillespie steps=%d\n",STEPS);
	fprintf(f2p,"---------------------------\n");
	fprintf(f2p,"j\tt[j]\tEtotal[j]\n");
	fprintf(f2p,"---------------------------\n");
	fprintf(f2p,"\n");
   fclose(f2p);

for(j=0; j<CYCLEMAX; j++)
   	{
		f2p=fopen(string2,"a");
		fprintf(f2p,"%d:\t%d\t%d\n",j,t[j],Etotal[j]);
		fclose(f2p);
	}

   f2p=fopen(string2,"a");
	fprintf(f2p,"---------------------------\n");
	fprintf(f2p,"---------------------------\n");
	fprintf(f2p,"\n");
   fclose(f2p);

//printf("START GILLESPIE\n");
//#########################################################################################
//#########################################################################################
//#########################################################################################      
//GILLESPIE START
//#########################################################################################
//#########################################################################################
//#########################################################################################


  THRESHOLD=0;

  while(t0<(Transient+STEPS))//do "2 x steps" iterations in total
    {
	t0++;

	if(t0>Transient && THRESHOLD==0)
        {
	  time_Transient=runtime;//real transient time
	  THRESHOLD=1;
        }

	a0=0.0; //sum up all propensities
  	
     	for(i=0; i<TRANSMAX; i++)	a0=a0+w[j][i];
    	

	state = 0;
        cycle=0;
	suma=0.0;

      
//####################################################################################
//CHOOSE REACTION
//####################################################################################


	r1=genrand_real3();//generate random numbers between 0 and 1
	r2=genrand_real3();
	
	
	tau = 1/a0*(log(1/r1));//dwell time of the next reaction
	runtime = runtime + tau;//update the real time
	

	j=0;


	while(j<CYCLEMAX)	
	   {
		if(j==0) sum_before=0.0;
		else 	 sum_before=accum_trans[j-1];


		if(r2<trans[j][2]+sum_before){state=2;deg=0;cycle=j;break;}
		else if (r2>=trans[j][2]+sum_before && r2 <trans[j][1]+trans[j][2]+sum_before){state=1;cycle=j;break;}
		else if(r2>=trans[j][1]+trans[j][2]+sum_before && r2< trans[j][2]+trans[j][1]+trans[j][0]+sum_before){state=0;cycle=j;break;}
		else if(r2>= trans[j][2]+trans[j][1]+trans[j][0]+sum_before && r2<trans[j][3]+trans[j][2]+trans[j][1]+trans[j][0]+sum_before){state=2;deg=1;cycle=j;break;}
		else j++;

	
		if(r2>=accum_trans[CYCLEMAX-1]) printf("ERROR choosing: r2=%f>%f\n",r2,accum_trans[CYCLEMAX-1]);	

	   }


//######################################################################################
//UPDATE mean levels
//############################################################################################
 if(t0>Transient)
     {//after transient time

	for(j=0;j<CYCLEMAX;j++)
	   {
	      mean_c[j]=mean_c[j]+tau*c[j]/t[j];
	      mean_u[j]=mean_u[j]+tau*u[j]/t[j];
	      mean_b[j]=mean_b[j]+tau*b[j]/t[j];		
           }  
     }

//############################################################################################
//UPDATE numbers:
//############################################################################################
	if(state==0)	//uncharged tRNA gets affiliated with synthetase
	   {
			u[cycle]--;
			b[cycle]++;
			E[cycle]--;
			if(u[cycle]<0 || b[cycle]>t[cycle]) printf("ERROR state 0\n");
	   }


	else if(state==1)	//bound tRNA gets charged
	   {
			b[cycle]--;
			E[cycle]++;
			c[cycle]++;
			if(b[cycle]<0 || c[cycle]>t[cycle]) printf("ERROR state 1\n");
	   }

	else if(state==2)	//charged tRNA used or complex dissociation
	   {
			if(deg==0)	c[cycle]--;			//usage
			else 		{b[cycle]--;	E[cycle]++;}	//dissociated

			u[cycle]++;
			if(c[cycle]<0 || b[cycle]<0 ||  u[cycle]>t[cycle]) printf("ERROR state 2\n");
           }

	else printf("ERROR states\n");


	if(E[j]<0) printf("ERROR:  E<0\n");//debugg



//############################################################################################
//UPDATE propensities:
//############################################################################################


     	w[cycle][0]=k1[cycle]*E[cycle]*u[cycle];	//binding transition 
	    
	w[cycle][1]=k2[cycle]*b[cycle]; //charging transition

	w[cycle][2]=k3*c[cycle]; 	//usage transition

	w[cycle][3]=k0*b[cycle]; 	//dissociation 

//############################################################################################
//UPDATE sum propensities
//############################################################################################
  
//reset sums
	a0=0.0; sum_w=0;	

for(j=0; j<CYCLEMAX; j++)
   	  {
  	for(i=0; i<TRANSMAX; i++)	sum_w=sum_w+w[j][i];//sum of all propensities
	sum_trans[j]=0.0;
	accum_trans[j]=0.0;
	}


//recalculate probabilities
for(j=0;j<CYCLEMAX;j++)
	{    	

	for(i=0;i<TRANSMAX;i++) 
		{
		 trans[j][i]=w[j][i]/sum_w;	//transition probability within cycle j from state i
		sum_trans[j]=sum_trans[j]+trans[j][i];	//sum of the transition probabilities per cycle
		}
	

	i=0;
	while(i<=j)
	{   
	accum_trans[j]=accum_trans[j]+sum_trans[i]; i++;
	}
	
	}


}

//#########################################################################################
//#########################################################################################
//#########################################################################################      
//GILLESPIE END
//#########################################################################################
//#########################################################################################
//#########################################################################################

real_t=floor(runtime-time_Transient+0.5);//real run time of the system as an integer

//############################################################################################
//PRINTs
//############################################################################################


for(j=0; j<CYCLEMAX; j++)
   {
	sprintf(string2,"OUTPUT_SynSeq/Levels_%d.dat",j);
	f2p=fopen(string2,"a");
	fprintf(f2p,"%f\t %f\t%f\t%f\n",d, mean_u[j]/real_t,mean_b[j]/real_t,mean_c[j]/real_t);
	fclose(f2p);
   }


//calculate computational time:
 double diff2;
 char string5[50];
 FILE *f2_p;

double STEPS_number;
STEPS_number=log10((1.0*STEPS));//indication of iteration loops

 time(&end);
  diff2=difftime(end,start);
  sprintf(string5,"OUTPUT_SynSeq/Comp_times_steps_%.0f.dat",STEPS_number);
  f2_p=fopen(string5,"a");
  fprintf(f2_p,"%.2f\t%.2lf s\t%.2lf min\t%.2lf h\t%.2lf d\n",d,diff2, diff2/60,diff2/3600, diff2/(3600*24));
  fclose(f2_p);


}//doxy loop//codecheck

//############################################################################################
//############################################################################################
//PARAMETERs (see VMAX.dat)

/*
synt	kcat	aa index	k1	copies	k3	aa



22920	8.310086	0	0.003	17	1.85	Ala
17540	5.932954	1	0.006	21	0.93	Arg
5950	19.795531	2	0.005	11	5.48	Asn
17660	8.166969	3	0.004	16	1.58	Asp
12970	1.894810	4	0.026	4	0.29	Cys
25830	3.373009	5	0.007	10	1.06	*Gln*
46100	3.721516	6	0.004	17	0.99	Glu
51950	3.023033	7	0.004	21	0.38	Gly
16730	2.872849	8	0.013	8	0.42	His
19790	7.680998	9	0.004	16	1.77	Ile
80830	2.620908	10	0.003	21	0.88	Leu
26230	6.996217	11	0.004	22	1.72	*Lys*
50790	0.974023	12	0.013	11	0.09	Met
10270	9.524099	13	0.007	11	2.15	Phe
13480	7.614081	14	0.006	12	2.18	Pro
18440	9.585765	15	0.004	18	2.53	Ser
30520	4.624517	16	0.005	17	0.82	Thr
11250	2.103773	17	0.027	6	0.19	Trp
9500	7.902359	18	0.009	8	1.87	Tyr
6750	24.349610	19	0.004	19	4.50	Val




*/

 printf("DONE!\n");
  return 0;
}