###################################################################################
//GTM_code.zip
###################################################################################

The GTM_code.zip contains:

	GTM_code.c   
	RUN_GTM
	mtt.c
	INPUT/
	OUTPUT/	

###################################################################################
//Global translation model
################################################################################### 
The GTM is a global mathematical model of translation of mRNAs by ribosomes. 
The model is based on the Totally Asymmetric Simple Exclusion Process (TASEP) and it includes competition for ribosomes and tRNAs among the mRNAs, as well as competition among tRNAs for the aminoacyl-tRNA synthetases. 
The model also incorporates ribosome drop-off, where with a given probability the ribosomes detach from the mRNA before reaching a stop codon. 
The mRNAs are represented by one-dimensional lattices, all together representing the yeast genome, with each site of the lattice being a codon. 
Ribosomes are represented by particles with a footprint of 9 codons that enter the lattice at a specific initiation rate and hop from one codon to the next until they reach the stop codon, at which point they leave the lattice, thereby producing a protein. 

Due to the autogeneous feedback loop, the number of ribosomes within the reservoir (=cytoplasm) is proportional to the rate of protein production. In the presence of doxycycline, the number of glutamine synthetases decreases. As a consequence, the ribosomal speed along the mRNA decreases and less proteins are produced per unit time, and so, less ribosomes. This in turn decreases the initiation rate. Hence, the ribosomal current also decreases, thereby decreasing the demand for glutamine synthetases.

The main outputs of the GTM are the average protein production rate (the ribosomal ‘current’ along the mRNA), which can be used as a proxy for the growth
rate of the cell, the density profils along the mRNAs, which can be used to detect queuing, and charging levels of the tRNAs.


###################################################################################
//SETUP
###################################################################################


//--------------------------------------------------------------------------------//
//General notes
//--------------------------------------------------------------------------------//

The default version is set to a factor=280 of the cell, i.e., 102 mRNAs, Gln4 ratio of 1 (no doxy), all switches on (ribosomal reservoir, different initiation rates for each mRNA, drop off, autogeneous feedback mechanism). Note that for the publication, a 5% cell was used.

The default version (single run) should take approx. 6 min. Using the doxycycline loop in line (262) will take approx. 145 min.
A summary of all (chosen) parameters can be found in OUTPUT/other/SETTINGS_j_ID_d_Ww.dat (where j corresponds to the scaling factor, d corresponds to the Gln4 protein ratio and w corresponds to the footprint).

Make sure, the INPUT and OUTPUT folders are within the same location as the code and the program.

The OUTPUT folder should contain the following subfolders: other, out_k, out_J, out_rho and out_Q. Note that most of the files would be replaced/overwritten if you run the program a second time.

If you want to change the name or location of the INPUT or OUTPUT or subfolders, your may find all relevant strings in lines (278-327).
Note that there are TWO (!) strings within the code, in line (2516) and line (2523), which cannot be put to the top, but needs to be changed along, if you change the name or location of the folders.


Due to reproducability, in line (768) the seed for the random number generator has been set arbitrarily default to <seed = 260>.
<seed=time_t> would use the internal clock as a seed



//--------------------------------------------------------------------------------//
//INPUT/
//--------------------------------------------------------------------------------//
The INPUT folder contains: 
(index i = task_id; index j = system size scaling factor)

	input_i.dat 		
	   (the different Gln4 protein ratios)

	tableAlphas_i.dat
	   (the 93 different initiation rate constants (cell size) for each GO-Slim category)

	tableC_i.dat 		
	   (codon indice table, 5 columns: amino acid, tRNA, codon, tRNA index, codon index within tRNA family)



	tableCAT_j_i.dat 		
	   (2 columns: GO-Slim category, number of mRNAs per category)

	tableSeq_j_i.dat 
	   (each row contains one mRNA sequence of the cell fraction, e.g., tableSeq_6098_i.dat contains 2 mRNAs.)

	tableSIZE_j_i.dat 
	   (each row shows the length (ORF) of the corresponding mRNA in tableSeq_j_i.dat)



	tableSeq_PDS_j_i.dat
	   (tableSeq_j_i.dat for the PDS system; PDS= perfect demand and supply (artificial setup); uses tableV_0_i.dat)

	tableSIZE_PDS_j_i.dat 
	   (tableSIZE_j_i.dat for the PDS system)


	
	tableV_0_i.dat 	
	   (charging process, 5 columns: Vmax_PDS (cell size), dummy, gene copy number of the tRNA, amino acid index, number of codons affiliated to the tRNA) 	
	   (Note: tRNAs within amino acid family are sorted corresponding to tableC_i.dat. )

	tableV_1_i.dat 	
	   (tableV_0_i.dat but Vmax_aaE (cell size); aaE=amino acid estimation)

	tableW_k_i.dat 
	   (wobble base pair table, 2 columns: tRNA index, wobble value)
	   (k=0: no wobble; k=1: G-U and I-G wobble; k=2: G-U and I-G wobble and missing tRNAs)



###################################################################################
//SETTINGS and possible VARIATIONs
###################################################################################


//--------------------------------------------------------------------------------//
//System size
//--------------------------------------------------------------------------------//

Adopt code for another fraction of the cell: 

table splitted:

47 /*
48 cell factor	1/20 	1/36	1/50	1/75	1/100 	1/280 		
49 no of mRNAs	1516	837	602	404	300	102		
50 runtime WT	48h	11h	6h	86min 	52 min	6min		
51 */

artifical cell:
47 /*
48 cell factor	1/252		1/279		1/24		1/6098
49 no of mRNAs	93		100		500		2
50 runtime WT	artifical cell	artifical PDS	artifical PDS	artifical PDS
51 */



Change the following lines in the header accordingly:

54 #define factor 280//scaling factor for Go slim cat sequences and sizes, downscaled from real cell;
55 #define Mnull 102 //number of mRNAs which determines the Systemsize
56 #define M 102 // M=M0+Mover; // actual number of mRnAs in the system

Note that
//--------------------------------------------------------------------------------//
//Gln4 protein ratio - doxycycline
//--------------------------------------------------------------------------------//

The input_%d.dat, which are loaded in line (278) define the Gln4 protein ratio

261 	int task_id=1;
262 	// for(task_id=1;task_id<4;task_id++){//codecheck input loop
2926 	// }//codecheck input loop

if you want to enter more Gln4 protein ratio values, you have to adjust the task_id loopin line (262). Don't forget to close it in line (2926)!
if you want to enter more then 13 different values, you have to add more task_ids, input_%d.dat files AND accordingly add copies of task_id dependent files in the INPUT folder, i.e., also e.g. "table_Seq_20_%d.dat"



//--------------------------------------------------------------------------------//
//Footprint, charging
//--------------------------------------------------------------------------------//


26 #define W 9//footprint width of the ribosome, how many codons are covered by the ribosome (literature: 9)


38 #define ValueVmax 1 //0:perfect demand and supply, 1: aa-estimation
39 #define ValueWobble 2 //0: no wobble, 1: normal wobbles 2: all wobbles (incl. 60%)
40
41 #define RealCell_No_codons 1.0e07 // number of codons in a real cell GOSlim (n_c_new)=1.0e07 //PDS (n_c_old)=3.5e06


//--------------------------------------------------------------------------------//
//SWITCHES for  ribosomal reservoir, footprint, dropoff and autogeneous feedback
//--------------------------------------------------------------------------------//
 
In general: A switch is turned on for 'switch'>0 and off for 'switch'=0 

28 #define pool 1 //ribosomal reservoir switch on (1) and off(0)

30 #define switch_alphas 1 //switch to associate different alphas for each mRNA/category; on (1) and off(0)
31			//off means all alphas are the same (differentiation between M0 and M_over mRNAs possible: alp1 and alp2)
32
33 #define DROPOFF_on 4 // 0: no dropoff; 1: normal, 4: stress, (gamma=0.0014*DROPOFF_on) 
34
35 #define feedback1 1 //switch for feedback on the number of ribosomes from doxy;  on (1) and off(0)
36	#define lambda 1.0 //strength of the feedback l=1, i.e., l=doxyscale; l<1 strong, l>1 weak



Switching of the ribosomal reservoir in line (28) corresponds to unlimited supply of ribosomes.
Switching of the the alphas in line (30) means that all mRNAs have the same initiation rate constant of alp=0.1/s (see line (636)) 
For drop off in line (33), we can choose different values:  0=no dropoff; 1=normal (0.0014), 4=stress (4*0.0014=0.0056),  
Feedback can be "tuned" by the factor lambda in line (36). Note that weak feedback is capped, feedback cannot be larger than doxyscale.


//--------------------------------------------------------------------------------//
//LOOPS
//--------------------------------------------------------------------------------//
In gerenal, if a loop is opened in one of the following lines, don't forget to close it at the very end of the code!

1.) Initiation rate: Note that this is only effective, if the switch in line (30) for the different alphas is set to 0 

652   //for (alp=0.0001;alp<10.1;alp=alp*10) {// begin alpha loop

2928  // } //alpha loop



2.) Scaling factors kappa for the Michaelis Menten constant Km (analogeous for the scaling factor nu feo the maximal velocity Vmax):

695 //for (kappa=0.025;kappa<0.21;kappa=kappa+0.02) // begin kappa loop
696 //{  xi=kappa/nu;scaling_K=eta*kappa;scaling_V=eta*nu;//scaling factor for vmax and Km: Vmax/Km=xi;

2929 //} //nu loop



3.) Charging efficiency \xi=Vmax/Km:

698 //for (xi=0.025;xi<1.01;xi=xi+0.025)  // begin xi loop
699 //for (xi=1.1;xi<2.01;xi=xi+0.1) 
700 //{kappa= xi*nu;scaling_K=eta*kappa;scaling_V=eta*nu;//scaling factor for vmax and Km: Vmax/Km=xi;

2930 //}//kappa loop


4.) Any kind of doxycycline scale via the Gln4 protein ratio different from input files.


//--------------------------------------------------------------------------------//
//OVEREXPRESSION of tRNAs
//--------------------------------------------------------------------------------//
We can overexpress one of the Gln tRNAs UUG (1) or CUG (2) with <GLNoverex> and the "strength" of the overexpressing by <GLNoverFactor>

68 #define GLNoverex 0 //switch for overexpression of Gln tRNAs: 0= no overexpression, 1: UUG overex. 2== CUG overexp.
69 #define GLNoverFactor 10 //overexpression factor for gln tRNAs

834 	if(GLNoverex==1)	{TRNA[10].tnum=GLNoverFactor*TRNA[10].tnum;}//UUG overexpressed
835	else if(GLNoverex==2){TRNA[9].tnum=GLNoverFactor*TRNA[9].tnum;}//CUG overexpressed


//--------------------------------------------------------------------------------//
//OVEREXPRESSION of mRNAs
//--------------------------------------------------------------------------------//

additional mRNAs in the system are interpreted as a stress situation. We define the number of extra mRNAs in line (58) and summarise all of these extra mRNAs in a single, additional category cat=93, the maximal number of categories then is 94, see line (63). In line (639) the corresponding initiation rate can be assigned to the extra mRNAs. 

58 #define Mover 0	//number of extra mRNAs, overexpression, stress situation // scan: (string1,"mRNAs_sequences_%d_o%d.dat",M0,Mover)

63 #define catovermax 94//catovermax=catmax+overexpression

639 alp2=alp;//overexpression parts




###################################################################################
//INSTALLATION and RUN
###################################################################################

This code is written in C.

Linux: 

No installation needed.

1.) The following files are required to be in the same directory:

	GTM_code.c   
	    (this is the original code)
	RUN_GTM
	    (this is the executable program)
	
 	mtt.c 	 
	    (get from http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html )
	
	INPUT/
 	    (content as described above)

	OUTPUT/ 
	    (content as described above)


2.) Open a terminal and go to the directory where your files are.


If you changed (parameters or settings within) the code go to step 3.a) otherwise go directly to step 3.b)

3.a) compile the code using 

	gcc -O3 -o RUN_GTM GTM_code.c -lm

(or set your own programm name intstead of 'RUN_GTM'.)

 
3.b) Run the programm with

	./RUN_GTM

(or use your own name you chose in 3.a) instead of 'RUN_GTM'.)


4.) The program creates files within the 'OUTPUT/' folder (details below)




###################################################################################
//OUTPUT
###################################################################################


After the successfull run of the simulation, the OUTPUT folder will contain the following files:

(
index i = task_id; 
index j = system size scaling factor; 
index d =Gln4 protein ratio from input_i.dat
index w = footprint, ribosomal codon coverage (default w=9)
index y = total number of mRNAs in cell fraction for scale factor j
index x = specific mRNA out of y
index a =specific GO-Slim category out of 93
)


out_J/
	gene_rho_J_cell_j_ID_d_Ww.dat:
	   (number of the mRNA, initiation rate alpha[i], mean density on mRNA i mean_rho[i], mean current J[i], number of proteins[i])

	meanJ_CAT_cell_j_ID_d_Ww.dat:
	   (Gln4protein ratio, category, mean current per cat, stdev)

	meanJ_cell_j_ID_d_Ww.dat: 
	   (Gln4 protein ratio, mean current, stdev)
	   (WITH mRNA overexpression: Gln4 protein ratio, mean current, stdev,
		 mean current through overexpression parts, stdev, mean current through all mRNAs,stdev, no. of drop_off events)


out_k/
	av_k_cellj_ID_d_Ww.dat:
	   (Gln4protein ratio, average hopping rate, stdev, smallest hopping rate, largest hopping rate)
	mean_k_cellj_ID_d_Ww.dat:
	   (codon index, mean hopping rate for codon i)


out_Q/
	average_Q_level_cell_j_ID_d_Ww.dat:
	   (Gln4protein ratio, av. charging level Q, stdev Q, inhomogenity marker, mean charging level CUG, mean charging level UUG, time, no. of dropoff events)
		
	mean_Q_cell_j_ID_d_Ww.dat:
	   (tRNA index, mean charging level of tRNA i, inhomogenity marker)

		

out_rho/
	density_profile_gene_x_of_y_cat_a_cell_j_ID_d_Ww.dat:
	   (lattice space i of the mRNA, density along the mRNA on lattice space i)



other/	

	No_DO_total_j_ID_d_Ww.dat
	   (Gln4protein ratio, total number of dropoff events, number of dropoff events from stopcodon)

	No_DO_events_j_ID_d_Ww.dat
	   (category index, mRNA index, lattice site i, number of drop off events on lattice site i)

	No_DO_eventsperCat_j_ID_d_Ww.dat
	   (category i, summ of dropoff events in category i)

	No_DO_eventsperCodon_j_ID_d_Ww.dat
	   (codon index, sum of dropoff events per codon)

	No_DO_eventsperMRNA_j_ID_d_Ww.dat
	   (category index of mRNA i, mRNA i, sum of dropoff events on mRNA i, double check of sum of dropoff events on mRNA i)

	No_ribos_cell_j_ID_d_Ww.dat:
	   (mRNA index, category index, number of ribosomes on mRNA)
	   (last line: total number of ribosomes N_t, average number of free ribosomes N_f, difference N_t-N_f, sum of ribosomes in transition)
		
	No_ribos_on_mRNA_j_ID_d_Ww.dat
	   (mRNA index , number of ribosomes on mRNA i, length of the ribosome array rib[j].len, length of the propensity array a[j].len)

	No_ribos_perCat_j_ID_d_Ww.dat
	   (category index, number of mRNAs in cat i, sum of all ribosomes on mRNAs in cat i, av. number N_mRNA if ribosomes per mRNA in cat i, stdev N_mRNA)

	Pos_ribos_on_mRNA_j_ID_d_Ww.dat
	   (snapshot at the end of the simulation)
	   (mRNA index i, ribosome index j, position (=lattice) of the ribsome j on mRNA i)

	time_j_ID_d_Ww.dat
	   (computational time: task_id, time in seconds, time in minutes, time in hours, time in days)

	SETTINGS_j_ID_d_Ww.dat: Example

--------------------SETTINGS-----------------------------
Number of reactions/Gillespie steps:
T_single=5.0e+05 	 T=5.1e+07 	 Transient=5.1e+07
 --------------------------------------------------------
Sys. No. codons all= 	 1.00e+07
 Scaling factor eta= 	 306.21
--------------------------------------------------------
Tiny system:
 	 No. mRNAs=	 102 (all=102) 
	 No. codons= 	 32657
	 No. tRNAs= 	 9797 
	 No. ribos= 	 714
--------------------------------------------------------
Switches off/on (0/>0): 
	 drop off= 	 4 
	 ribopool=	 1
 	 feedback=	 1	 strength (lamda)=1.000000
--------------------------------------------------------
	Vmax value=	 1 	 (0:PDS; 1: aaE)
	Wobble value= 	 2 	 (0:no; 1: two wobbles, 2: all wobbles)
	Footprint= 	 9
Transcription rate Tt=	 3.40e+01
Hopping rate factor r=	 5.71e-02
Rates: 
 	 alphas vary 
 	 beta= 		 38.81 
	 gamma=		 5.6e-03
--------------------------------------------------------
Charging rate scaling factors: 
	 xi=kappa/nu=	 5.000
	 kappa=		 5.000
	 nu=		 1.000
--------------------------------------------------------
Gln4 protein ratio= 	 0.050	 Doxycycline [mug/ml]=(0.640000)^(4.05/3)
--------------------------------------------------------
--------------------------------------------------------





OUTPUT on screen (in terminal): shows progress (iteration steps, transient marker)

1.00e+00	0
1.00e+01	0
1.00e+02	0
1.00e+03	0
1.00e+04	0
1.00e+05	0
1.00e+06	0
1.00e+07	0
1.00e+08	1

DONE!





###################################################################################
//ATTRIBUTION PARTIES
###################################################################################

The GTM model as described here has been developed by Corina Keller and M.Carmen Romano.
It is based on the 'Multiple TASEP sharing tRNAs'-model (2nd October 2015) coded by Adélaïde Raguin and M.Carmen Romano.
Input data provided by Ian Stansfield.
The work has been carried out at University of Aberdeen, UK, and has been funded by BBSRC.


-----------------------------------------------------------------------------------------------
//AUTHORSHIP
-----------------------------------------------------------------------------------------------
The authors of the GTM code are Corina Keller, Adélaïde Raguin and M.Carmen Romano.


-----------------------------------------------------------------------------------------------
//PUBLICATIONS
-----------------------------------------------------------------------------------------------
The GTM model was used to support the following publication. If you use the code, please refer to this publication:

- to be added - 
(*
submitted to NUCLEIC ACIDS RESEARCH

The molecular aetiology of tRNA synthetase depletion: induction of a GCN4 amino acid starvation response despite homeostatic maintenance of charged tRNA levels.
Matthew R. McFarland, Corina Keller, Brandon Childers, Holly Corrigall, Adélaïde Raguin, M. Carmen Romano, Ian Stansfield.
*)


-----------------------------------------------------------------------------------------------
//CORRESPONDENCE
-----------------------------------------------------------------------------------------------
i.stansfield@abdn.ac.uk
m.romano@abdn.ac.uk
corina.keller@abdn.ac.uk




###################################################################################
//LICENSE
###################################################################################

Creative Commoms BY-NC find here: https://creativecommons.org/licenses/by-nc/4.0/

You are free to:

    Share — copy and redistribute the material in any medium or format
    Adapt — remix, transform, and build upon the material

    The licensor cannot revoke these freedoms as long as you follow the license terms.

Under the following terms:

    Attribution — You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.

    NonCommercial — You may not use the material for commercial purposes.

    No additional restrictions — You may not apply legal terms or technological measures that legally restrict others from doing anything the license permits.

Notices:

    You do not have to comply with the license for elements of the material in the public domain or where your use is permitted by an applicable exception or limitation.
    No warranties are given. The license may not give you all of the permissions necessary for your intended use. For example, other rights such as publicity, privacy, or moral rights may limit how you use the material.