/* kai.cpp 2006 Harvard iGEM team This program simulates KaiC's phosphorylation cycle in a cell. Arguments: -t time: time to run -o filename: path of output file -b beta: probability that a given protein will decay at a given time step -i initial: initial number of proteins -a alpha_function: transcribe new protein as a function of time. constant: produce protein at a constant amount per time step sine: produce protein in a sinusoidal pattern feedback: produce protein in a self-repressing feedback pattern -c alpha_constant: has different meanings depending on which alpha function is selected constant: the amount of protein to produce per timestep sine: the average amount of protein to produce per timestep (y-offset of the sine wave) feedback: the constant multiplier in the feedback equation [-m amplitude]: amplitude of the sine wave (only used if -a sine is specified) [-k k_m): Michaelis-Menten constant (only used if -a feedback is specified) [-s]: sparse output (comma delimited numbers, suitable for graphing) Examples: ./kai -t 1000 -o output.txt -b 0.03 -i 1000 -a constant -c 30 -s ./kai -t 1000 -o output.txt -b 0.03 -i 1000 -a feedback -c 60 -k 0.5 */ #include //#include //#include #include #include #include #include #include using namespace std; #define PI 3.1415926 /* globals */ int t = 0; /* current time */ int period = 24; /* period of one kai protein's oscillation */ double beta = 0; /* probability that any given protein will decay at a given time step */ enum alpha_functions {constant, sine, feedback} alpha; int alpha_constant = -1; int amplitude = -1; float k_m = -1; int sparse = 0; double last_P = 1.0; /* proportion of phosphorylated protein at the last oscillation */ struct kai { int step; int P; kai():step(0),P(0){} int phosphorylate(){ return (this->P = ((this->step < (period / 2)) ? 1 : 0)); } }; int alpha_function_constant(vector & kai_vector) { int i; for(i = 0; i < alpha_constant; i++) kai_vector.push_back(kai()); return 0; } int alpha_function_sine(vector kai_vector) { int i, amount; double x = (t % period) * PI / (double)period; amount = (int)(amplitude * sin(x) + alpha_constant); for(i = 0; i < amount; i++) kai_vector.push_back(kai()); return 0; } int alpha_function_feedback(vector kai_vector) { int i, amount; amount = (int)(alpha_constant * (k_m / (k_m + last_P))); for(i = 0; i < amount; i++) kai_vector.push_back(kai()); return 0; } /* Called at every timestep to update the array of kai proteins. For each protein, this function does the following: * Randomly decide to kill it or not (representing protein degradation) * If it survives, * Increment its step by 1 * Phosphorylate or unphosphorylate this protein as necessary */ /*double update_kai_array(struct array *kai_array) { int sum = 0; int i, size; struct kai *this_kai; //Add new proteins as necessary switch(alpha) { case constant: if(alpha_function_constant(kai_array)) { fprintf(stderr, "Error: insufficient memory\n"); exit(1); } break; case sine: if(alpha_function_sine(kai_array)) { fprintf(stderr, "Error: insufficient memory\n"); exit(1); } break; case feedback: if(alpha_function_feedback(kai_array)) { fprintf(stderr, "Error: insufficient memory\n"); exit(1); } break; default: fprintf(stderr, "Error: bad alpha function\n"); exit(1); break; } for(i = 0; i < array_getnum(kai_array); i++) { this_kai = array_getguy(kai_array, i); // Ask the RNG to see if this protein dies if((rand() / ((double)RAND_MAX + 1)) < beta) { array_remove(kai_array, i); continue; } // If the protein lives, increment its step and set its phosphorylation state this_kai->step = (this_kai->step + 1) % period; sum += phosphorylate(this_kai); } size = array_getnum(kai_array); if(size > 0) { return (last_P = sum / (double)array_getnum(kai_array)); } else { return 0; } }*/ double update_kai_vector(vector & kai_vector) { int sum = 0; /* Add new proteins as necessary */ switch(alpha) { case constant: if(alpha_function_constant(kai_vector)) { exit(1); } break; case sine: if(alpha_function_sine(kai_vector)) { fprintf(stderr, "Error: insufficient memory\n"); exit(1); } break; case feedback: if(alpha_function_feedback(kai_vector)) { fprintf(stderr, "Error: insufficient memory\n"); exit(1); } break; default: fprintf(stderr, "Error: bad alpha function\n"); exit(1); break; } vector::iterator p = kai_vector.begin(); while (true){ if((rand() / ((double)RAND_MAX + 1)) < beta) p = kai_vector.erase(p); else { (*p).step = (p->step + 1) % period; sum += p->phosphorylate(); p++; } if(p == kai_vector.end()) break; } if(kai_vector.size() > 0) return (last_P = sum / kai_vector.size()); else return 0; } void print_usage(char *error_string) { fprintf(stderr, "%s\n", error_string); fprintf(stderr, "Usage:\n\t-t time: time to run\n\t-o filename: path of output file\n\t-b beta: probability that a given protein will decay at a given time step\n\t-i initial: initial number of proteins\n\t-a alpha_function: transcribe new protein as a function of time (defaults to constant).\n\t\tconstant: produce protein at a constant amount per time step\n\t\tsine: produce protein in a sinusoidal pattern\n\t\tfeedback: produce protein in a self-repressing feedback pattern\n\t-c alpha_constant: has different meanings depending on which alpha function is selected\n\t\tconstant: the amount of protein to produce per timestep\n\t\tsine: the average amount of protein to produce per timestep (y-offset of the sine wave)\n\t\tfeedback: the constant multiplier in the feedback equation\n\t[-m amplitude]: amplitude of the sine wave (only used if -a sine is specified)\n\t[-k k_m): Michaelis-Menten constant (only used if -a feedback is specified)\n\t[-s]: sparse output (comma delimited numbers, suitable for graphing) "); exit(1); } struct opt_arg{ char * name, * val; }; class get_args{ public: get_args(int _argc, char ** _argv):num(0),argc(_argc),argv(_argv){} opt_arg get(const vector value_args){ opt_arg ret; scanf(argv[num],"-%s",&(ret.name)); if (find(value_args.begin(),value_args.end(),ret.name) != value_args.end()){ if (num+1 == argc){ printf("Invalid argument"); exit(1); } ret.val = argv[num+1]; num+=2; } else num++; return ret; } bool done(){ return num < argc; } int num; int argc; char ** argv; }; int main (int argc, char *argv[]) { int i; int t_max = -1; int initial_size = -1; FILE *output; char *output_name = NULL; //struct array *kai_array; vector kai_vector; //struct kai *this_kai; double P_percent; /* parse arguments */ if(argc < 7) { print_usage("Too few arguments given\n"); } get_args args(argc, argv); vector a; a.push_back("t"); a.push_back("o"); a.push_back("b"); a.push_back("i"); a.push_back("a"); a.push_back("c"); a.push_back("m"); a.push_back("k"); while(!args.done()){ opt_arg cur_arg = args.get(a); char * optarg = cur_arg.val; switch(cur_arg.name[0]) { case 't': t_max = atoi(optarg); if(t_max < 0) { print_usage("t_max must be >= 0"); } break; case 'o': output_name = optarg; break; case 'b': beta = atof(optarg); if(beta < 0 || beta > 1) { print_usage("beta must be >= 0 and <= 1)"); } break; case 'i': initial_size = atoi(optarg); if(initial_size < 0){ print_usage("initial_size must be >= 0"); } break; case 'a': if((strcmp(optarg, "constant")) == 0) { alpha = constant; } else if((strcmp(optarg, "sine")) == 0) { alpha = sine; } else if ((strcmp(optarg, "feedback")) == 0) { alpha = feedback; } else { print_usage("Bad alpha function specified\n"); } break; case 'c': alpha_constant = atoi(optarg); if(alpha_constant < 0) { print_usage("alpha_constant must be >= 0"); } break; case 'm': amplitude = atoi(optarg); if(amplitude <0) { print_usage("amplitude must be < 0\n"); } break; case 'k': k_m = atof(optarg); if(k_m < 0) { print_usage("k_m must be < 0\n"); } break; case 's': sparse = 1; break; default: print_usage("Unknown argument\n"); } } /* Verify that reuired arguments are present */ if(t_max < 0) { print_usage("-t time is a required argument"); } if(initial_size < 0) { print_usage("-i initial_size is a required argument"); } if(alpha_constant < 0) { print_usage("-c alpha_constant is a required argument"); } if(alpha == sine && amplitude < 0) { print_usage("-m amplitude is required if -a sine is specified"); } if(alpha == feedback && k_m < 0) { print_usage("-k k_m is required if -a feedback is specified"); } output = fopen(output_name, "w"); if(output == NULL) { print_usage("Failed to open output file"); } if(sparse) { fprintf(output, "Arguments:\nt_max: %d,\noutput file: %s,\ninitial amount of protein: %d,\nbeta: %f,\nalpha function: %i,\nalpha_constant: %i \namplitude: %i \nk_m: %f \n", t_max, output_name, initial_size, beta, alpha, alpha_constant, amplitude, k_m); } /* Initialize array of kai proteins */ //kai_array = array_create(); for(i = 0; i < initial_size; i++) { kai_vector.push_back(kai()); //this_kai = malloc(sizeof(struct kai)); //if(this_kai == NULL) { // fprintf(stderr, "Error: insufficient memory\n"); //} //this_kai->step = 0; //this_kai->P = 0; //array_add(kai_array, this_kai); } /* Run! */ for(t = 0; t < t_max; t++) { if(sparse) { fprintf(output, "%d,%d,%f\n", t, kai_vector.size(), update_kai_vector(kai_vector)); } else { fprintf(output, "t=%d: size=%d, %f\n", t, kai_vector.size(), update_kai_vector(kai_vector)); } } fclose(output); return 0; }