sigma32.cpp

用数值方法模拟了E.coli中sigma32因子对热休克(Heat shock)基因的调控机理

#include <stdio.h>
#include <math.h>
#include <iostream.h>

#define Ve 2.54e-9
#define a0 0.3
#define ad 0.3
#define as 4
#define Kd 20
#define enta(T) ((T>40)?200:40)
#define Kfold 15.0e3
#define K(T) ((T>40)?150:75)
#define a 0.2
#define DT 0.0001
#define Ku (1.0e6*Ve)
#define Ks (1.0e7*Ve)
#define Kf (1.0e6*Ve)

double S=40.0,Sf;
double D=10000.0,Df;
double U=6.0e5;
double P=2.0e6;

int T=37;
long int stime=-1000;

void go();
void print(FILE *fp);

void main()
{
	FILE *fp;
	fp=fopen("out.txt","w");
	int i=0,j;

	fprintf(fp,"stime\tS\tSf\tD\tDf\tPfolded\n");

	//等待达到稳态
	do
	{
		for(j=0;j<100;j++)go();
		i++;
	}
	while(i<2000);

	//初始化曲线
	i=0;
	do
	{
		for(j=0;j<100;j++)go();
		print(fp);
		i++;
		stime++;
	}
	while(i<1000);

	//升温
	i=0;
	T=42;
	do
	{
		for(j=0;j<100;j++)go();
		print(fp);
		i++;
		stime++;
	}
	while(i<4000);

	fclose(fp);
}

void go()
{
	double dS,dD;

	Df=K(T)*K(T)+2*Ku*K(T)*K(T)+2*Kfold*Ku*D*K(T)+Ku*Ku*K(T)*K(T)*P*P-2*Ku*Ku*K(T)*P*D*Kfold+D*D*Kfold*Kfold*Ku*Ku;
	Df=pow(Df,0.5);
	Df+=D*Kfold*Ku-K(T)-Ku*K(T)*P;
	Df/=2*Kfold*Ku;
	Sf=S/(1+Ks*Df*(1+Kf*a*D));
	U=(K(T)*P)/(Kfold*Ku*Df+K(T));

	dS=enta(T)-a0*S-as*Kf*Ks*a*D*Sf*Df;
	dD=Kd*Sf-ad*D;

	S+=DT*dS;
	D+=DT*dD;
}

void print(FILE *fp)
{
	fprintf(fp,"%d\t%e\t%e\t%e\t%e\t%e\n",stime,S,Sf,D,Df,P-U);
}