IGEM:Harvard/2006/Cyanobacteria/Model: Difference between revisions
JeffreyLau (talk | contribs) No edit summary |
JeffreyLau (talk | contribs) No edit summary |
||
Line 1: | Line 1: | ||
{{IGEM:/Harvard/2006/Cyanobacteria}} | {{IGEM:/Harvard/2006/Cyanobacteria}} | ||
<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
- /
- 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;
}