IGEM:Harvard/2006/Cyanobacteria/Model: Difference between revisions

From OpenWetWare
Jump to navigationJump to search
No edit summary
No edit summary
Line 1: Line 1:
{{IGEM:/Harvard/2006/Cyanobacteria}}
{{IGEM:/Harvard/2006/Cyanobacteria}}


Stuff
<code>
/*
kai.c
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 <stdlib.h>
//#include <unistd.h>
//#include <getopt.h>
#include <stdio.h>
#include <vector>
#include <math.h>
#include <string.h>
#include <algorithm>
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> & kai_vector)
{
int i;
for(i = 0; i < alpha_constant; i++)
kai_vector.push_back(kai());
return 0;
}
int alpha_function_sine(vector<kai> 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> 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> & 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<kai>::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<char *> 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> 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<char *> 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;
}
</code>

Revision as of 11:43, 24 July 2006

<html><style type='text/css'> .tabs { 

 font-size:80%;
 font-weight:none;
 width: 100%;
 color: #FFFFFF;
 background:#FFFFFF url("/images/5/54/DarkgreenTab-bg.gif") repeat-x bottom;

}

.tabs li { 

 background:url("/images/3/36/DarkgeenTab-left.gif") no-repeat left top;

}

.tabs a,.tabs strong { 

 background:url("/images/d/d3/DarkgreenTab-right.gif") no-repeat right top;
 color:#FFFFFF;
 padding: 3px 10px 3px 4px;

}

.tabs strong{ 

 color:#CCFF00;
 background-image:url("/images/b/b1/DarkgreenTab-right_on.gif");

}

.tabs a:hover{ 

 color:#66FF00;

}


</style></html>


/* kai.c 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

  • /


  1. include <stdlib.h>

//#include <unistd.h> //#include <getopt.h>

  1. include <stdio.h>
  2. include <vector>
  3. include <math.h>
  4. include <string.h>
  5. include <algorithm>

using namespace std;

  1. 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> & kai_vector) { int i;

for(i = 0; i < alpha_constant; i++) kai_vector.push_back(kai()); return 0; }

int alpha_function_sine(vector<kai> 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> 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> & 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<kai>::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<char *> 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> 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<char *> 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; }

Retrieved from "https://openwetware.org/mediawiki/index.php?title=IGEM:Harvard/2006/Cyanobacteria/Model&oldid=53243"