// LAST UPDATE: Oct. 18, 2007. #include #include #include #include typedef int TYPE; typedef TYPE *ARRAY; typedef double TYPE1; typedef TYPE1 *ARRAY_d; typedef int TYPE2; typedef TYPE2 **MATRIX; #define pi 3.14159265358979 #define sqrtPI sqrt(pi) #define MIN(x,y) ((xy) ? x:y) struct Mol_bound_Rec { ARRAY GEF, PI3K; ARRAY GEF_Ras, PI3K_Ras; double x,y; }; struct Mol_bound_Rec *Rec_bound; struct Mol_bound_Ras { int Enz; // 0 for FREE; 1 for GEF; 2 for PI3K int state; // 0 for GDP; 1 for GTP int Rec_bound; // 0 for receptor-free; 1 for receptor-bound int Rec_closest; // index of the closest receptor double x,y; // coordinates double d_min; // distance to the closest receptor }; struct Mol_bound_Ras *Ras_bound; struct List_of_2 { int i; int s; }; ARRAY Ras_GDP_free, Ras_GDP_GEF, Ras_GTP_free, Ras_GTP_PI3K, Ras_Inside; ARRAY p_Ras_GDP_free, p_Ras_GDP_GEF, p_Ras_GTP_free, p_Ras_GTP_PI3K,InDist; int N_Ras_GDP_free, N_Ras_GDP_GEF, N_Ras_GTP_free, N_Ras_GTP_PI3K; struct List_of_2 *Rec_GEF_free,*Rec_GEF_bound,*Rec_PI3K_free,*Rec_PI3K_bound; MATRIX p_Rec_GEF_free, p_Rec_GEF_bound, p_Rec_PI3K_free, p_Rec_PI3K_bound; int N_Rec_GEF_free, N_Rec_GEF_bound, N_Rec_PI3K_free, N_Rec_PI3K_bound, N_Ras_Inside; ARRAY_d rate,constant,Loop_I,Loop_II,Loop_IIIA,Loop_IIIB,Loop_IV,Loop_V; double T_Loop_I=0., T_Loop_II=0., T_Loop_IIIA=0., T_Loop_IIIB=0., T_Loop_IV=0., T_Loop_V=0.; int N_Loop_I=0, N_Loop_II=0, N_Loop_IIIA=0, N_Loop_IIIB=0, N_Loop_IV=0, N_Loop_V=0; int Rec_n,Ras_n; double S,dS,eta_R,TSTEP,Diff,phi_RE,tau_RE,tau_SE,kci_E,kappa_RE,phi_act,phi_in,phi_RP, tau_RP,tau_SP,kci_P,kappa_RP,K_R_ES,K_RE_S,K_R_PS,K_RP_S,t_end,t_print,a_tot, tau_GAP,TIME,Dt_next,SIGMA,mu_plus,Lxy,r_bimol; #define GRID_POINTS 20000 // Required for evaluation of Error function typedef struct {double XMIN, XMAX, DX; double Table[GRID_POINTS][2];} LookUpTable; LookUpTable ERFCX, *ptr_erfcx=&ERFCX, ERF, *ptr_erf=&ERF; void Create_Molecules(),Update_Ras_molecules(int,int,int,int), Update_Rec_molecules(int,int,int,int,int),Evaluate_rate_constants(), Add_Ras_site(int *,int *,int,int *),Remove_Ras_site(int *,int *,int,int *), Add_Rec_site(struct List_of_2 *,int *,int,int,int **),Evaluate_other_parameters(), Remove_Rec_site(struct List_of_2 *,int *,int,int **),PrintInDist(), Read_Parameters(char []),Output(),Ras_Rec_binding(), Circle_Point_Picking(double,double *,double *),PeriodicBC(int),Count_PI3K_time(int), Gauss_propagation(double,double,double,double,double,double *,double *,double *,double *), Min_Dist(double,double,int *,double *),GET_DATA(LookUpTable *,LookUpTable *), Update_Ras_molecules_2D_release(int,int,double),Update_Rec_molecules_2D_release(int,int,int,double), Release_Ras(int,int,double),Count_GEF_time(int),Count_GEF_both(int,int),Count_PI3K_both(int,int), Count_GEF_PI3K(int,int),Count_inactive_particles(); void Make_1st_order_reaction_step(int,double),Update_1st_order_rates(),Ras_diffusion(double); int Pick_next_reaction(),Make_2nd_order_reaction_step(int,int); int GEF_PI3K_complex(int); double RAN(),PDF_React(double,double,double),PDF_Reflect(double,double),GAUSS(), EVALUATE(double,LookUpTable *,int),Rotation_angle(double,double), MIN_DIST(double [],double [],double),Ras_diffusion_step(),Diffusion_time_step(); int P_bound[2],E_bound[2]; double Nin=0., Nin_times=0., Nact=0., Nact_rec=0., t_PI3K_0[2], T_P[2], t_GEF_0[2], T_E[2], tt_PI3K_0=0., tt_GEF_0=0., TT_P=0., TT_E=0., tt_EP_0=0., TT_EP=0., n_ES=0., n_PS=0., n_E_tot=0.,n_P_tot=0.,n_RP_S=0.; FILE *OutputFile, *OutputLoop, *OutputInactivation, *alterm; void Test_Ras_list(int,int *,int *), Test_Rec_list(int,struct List_of_2 *,int **), Test_Ras_mol(int),Test_Rec_mol(int); double gtp_free=0., pi3k=0., pi3k_bulk=0., n_count=0., gdp_free=0., gef=0., gef_bulk=0.; void Ras_GTP_and_PI3K(),Ras_GDP_and_GEF(); main(){ char name[50]; int R_type; double t_print_last,Dt_react,Dt_diff,p,t; sprintf(name,"Ras_model.data"); Read_Parameters(name); Evaluate_other_parameters(); Evaluate_rate_constants(); Create_Molecules(); Count_inactive_particles(); // **** INITIALIZE RECEPTOR OCCUPANCY TIMES AND Ras-BINDING EVENTS ****** P_bound[0]=0; P_bound[1]=0; t_PI3K_0[0]=0.; t_PI3K_0[1]=0.; T_P[0]=0.; T_P[1]=0.; E_bound[0]=0; E_bound[1]=0; t_GEF_0[0]=0.; t_GEF_0[1]=0.; T_E[0]=0.; T_E[1]=0.; // ********************************************************************** sprintf(name,"out.data"); OutputFile=fopen(name,"w"); sprintf(name,"loop.data"); OutputLoop=fopen(name,"w"); sprintf(name,"inact.data"); OutputInactivation=fopen(name,"w"); sprintf(name,"alterm.data"); alterm=fopen(name,"w"); // Add_Rec_site(Rec_GEF_bound,&N_Rec_GEF_bound,0,0,p_Rec_GEF_bound); // Update_Rec_molecules(0,1,0,1,-1); // Remove_Rec_site(Rec_GEF_free,&N_Rec_GEF_free,p_Rec_GEF_free[0][0],p_Rec_GEF_free); TIME = 0.; t_print_last = 0.; do{ //Dt_diff=Diffusion_time_step(); Dt_diff=0.01; Update_1st_order_rates(); if(a_tot>0.){ Dt_react=-log(RAN())/a_tot; if(Dt_diff>Dt_react){ R_type=Pick_next_reaction(); Make_1st_order_reaction_step(R_type,Dt_react); Dt_next=Dt_react; }else{ //p=Dt_diff/Dt_react; //if(p > RAN()){ // R_type=Pick_next_reaction(); // Make_1st_order_reaction_step(R_type,Dt_diff); //} Dt_next=Dt_diff; } }else{ Dt_next=Dt_diff; } Ras_diffusion(Dt_next); /*Update_1st_order_rates(); if(a_tot>0.){ Dt_next=-log(RAN())/a_tot; Ras_diffusion(Dt_next); Update_1st_order_rates(); R_type=Pick_next_reaction(); Make_1st_order_reaction_step(R_type,Dt_next); }else{ Dt_next=Diffusion_time_step(); Ras_diffusion(Dt_next); }*/ TIME+=Dt_next; Ras_GTP_and_PI3K(); //Ras_GDP_and_GEF(); if( (TIME-t_print_last)>t_print ){ t_print_last += t_print; Output(); } }while(TIME Ras-GDP-GEF ****** k=(int)(N_Ras_GDP_free*RAN()); i_Ras=Ras_GDP_free[k]; Add_Ras_site(Ras_GDP_GEF,&N_Ras_GDP_GEF,i_Ras,p_Ras_GDP_GEF); Remove_Ras_site(Ras_GDP_free,&N_Ras_GDP_free,k,p_Ras_GDP_free); Update_Ras_molecules(i_Ras,0,1,0); break; case 2: // ************************** Ras-GDP-GEF -> Ras-GDP + GEF ****** k=(int)(N_Ras_GDP_GEF*RAN()); i_Ras=Ras_GDP_GEF[k]; Add_Ras_site(Ras_GDP_free,&N_Ras_GDP_free,i_Ras,p_Ras_GDP_free); Remove_Ras_site(Ras_GDP_GEF,&N_Ras_GDP_GEF,k,p_Ras_GDP_GEF); Update_Ras_molecules_2D_release(i_Ras,1,dt); Update_Ras_molecules(i_Ras,0,0,0); break; case 3: // ************************** Ras-GDP-GEF -> Ras-GTP + GEF ****** k=(int)(N_Ras_GDP_GEF*RAN()); i_Ras=Ras_GDP_GEF[k]; Count_inactive_particles(); Nact++; if( (Ras_bound+i_Ras)->Rec_bound==1 ){ Nact_rec++; // Start Loop I Loop_I[i_Ras]=TIME; // Start Loop II Loop_II[i_Ras]=TIME; }else{ // Start Loop IV Loop_IV[i_Ras]=TIME; // Start Loop V Loop_V[i_Ras]=TIME; } // ************************************************ Add_Ras_site(Ras_GTP_free,&N_Ras_GTP_free,i_Ras,p_Ras_GTP_free); Remove_Ras_site(Ras_GDP_GEF,&N_Ras_GDP_GEF,k,p_Ras_GDP_GEF); Update_Ras_molecules_2D_release(i_Ras,1,dt); Update_Ras_molecules(i_Ras,1,0,0); break; case 4: // ************************** Ras-GTP -> Ras-GDP **************** k=(int)(N_Ras_GTP_free*RAN()); i_Ras=Ras_GTP_free[k]; Add_Ras_site(Ras_GDP_free,&N_Ras_GDP_free,i_Ras,p_Ras_GDP_free); Remove_Ras_site(Ras_GTP_free,&N_Ras_GTP_free,k,p_Ras_GTP_free); Update_Ras_molecules(i_Ras,0,0,0); // Destroy Loop I if(Loop_I[i_Ras]!=0.) Loop_I[i_Ras]=0.; // Destroy Loop II if(Loop_II[i_Ras]!=0.) Loop_II[i_Ras]=0.; // Destroy Loop IIIB if(Loop_IIIB[i_Ras]!=0.) Loop_IIIB[i_Ras]=0.; // Destroy Loop IV if(Loop_IV[i_Ras]!=0.) Loop_IV[i_Ras]=0.; // Destroy Loop V if(Loop_V[i_Ras]!=0.) Loop_V[i_Ras]=0.; // ************************************************ fprintf(OutputInactivation,"%f %e\n",TIME,(Ras_bound+i_Ras)->d_min); InDist[(int)((Ras_bound+i_Ras)->d_min+0.5)]++; break; case 5: // ************************** Ras-GTP + PI3K -> Ras-GTP-PI3K **** k=(int)(N_Ras_GTP_free*RAN()); i_Ras=Ras_GTP_free[k]; Add_Ras_site(Ras_GTP_PI3K,&N_Ras_GTP_PI3K,i_Ras,p_Ras_GTP_PI3K); Remove_Ras_site(Ras_GTP_free,&N_Ras_GTP_free,k,p_Ras_GTP_free); Update_Ras_molecules(i_Ras,1,2,0); // ************************************************ break; case 6: // ************************** Ras-GTP-PI3K -> Ras-GTP + PI3K **** k=(int)(N_Ras_GTP_PI3K*RAN()); i_Ras=Ras_GTP_PI3K[k]; // Destroy Loop IIIA if(Loop_IIIA[i_Ras]!=0.) Loop_IIIA[i_Ras]=0.; // Start Loop IIIB if( (Ras_bound+i_Ras)->Rec_bound==1 ) Loop_IIIB[i_Ras]=TIME; // ************************************************ Add_Ras_site(Ras_GTP_free,&N_Ras_GTP_free,i_Ras,p_Ras_GTP_free); Remove_Ras_site(Ras_GTP_PI3K,&N_Ras_GTP_PI3K,k,p_Ras_GTP_PI3K); Update_Ras_molecules_2D_release(i_Ras,2,dt); Update_Ras_molecules(i_Ras,1,0,0); break; case 7: // ************************** Rec + GEF -> Rec-GEF ************** k=(int)(N_Rec_GEF_free*RAN()); i_Rec=(Rec_GEF_free+k)->i; s_Rec=(Rec_GEF_free+k)->s; m0=(Rec_bound+i_Rec)->GEF[0]+(Rec_bound+i_Rec)->GEF[1]; k0=GEF_PI3K_complex(i_Rec); Add_Rec_site(Rec_GEF_bound,&N_Rec_GEF_bound,i_Rec,s_Rec,p_Rec_GEF_bound); Remove_Rec_site(Rec_GEF_free,&N_Rec_GEF_free,k,p_Rec_GEF_free); Update_Rec_molecules(i_Rec,1,s_Rec,1,-1); // Update receptor occupancy times Count_GEF_time(s_Rec); Count_GEF_both(m0,(Rec_bound+i_Rec)->GEF[0]+(Rec_bound+i_Rec)->GEF[1]); Count_GEF_PI3K(k0,GEF_PI3K_complex(i_Rec)); n_E_tot++; // ******************************** break; case 8: // ************************** Rec-GEF -> Rec + GEF ************** k=(int)(N_Rec_GEF_bound*RAN()); i_Rec=(Rec_GEF_bound+k)->i; s_Rec=(Rec_GEF_bound+k)->s; m0=(Rec_bound+i_Rec)->GEF[0]+(Rec_bound+i_Rec)->GEF[1]; k0=GEF_PI3K_complex(i_Rec); Add_Rec_site(Rec_GEF_free,&N_Rec_GEF_free,i_Rec,s_Rec,p_Rec_GEF_free); Remove_Rec_site(Rec_GEF_bound,&N_Rec_GEF_bound,k,p_Rec_GEF_bound); Update_Rec_molecules_2D_release(i_Rec,s_Rec,1,dt); Update_Rec_molecules(i_Rec,1,s_Rec,0,-1); // Update receptor occupancy times Count_GEF_time(s_Rec); Count_GEF_both(m0,(Rec_bound+i_Rec)->GEF[0]+(Rec_bound+i_Rec)->GEF[1]); Count_GEF_PI3K(k0,GEF_PI3K_complex(i_Rec)); // ******************************** break; case 9: // ************************** Rec + PI3K -> Rec-PI3K ************ k=(int)(N_Rec_PI3K_free*RAN()); i_Rec=(Rec_PI3K_free+k)->i; s_Rec=(Rec_PI3K_free+k)->s; m0=(Rec_bound+i_Rec)->PI3K[0]+(Rec_bound+i_Rec)->PI3K[1]; k0=GEF_PI3K_complex(i_Rec); Add_Rec_site(Rec_PI3K_bound,&N_Rec_PI3K_bound,i_Rec,s_Rec,p_Rec_PI3K_bound); Remove_Rec_site(Rec_PI3K_free,&N_Rec_PI3K_free,k,p_Rec_PI3K_free); Update_Rec_molecules(i_Rec,2,s_Rec,1,-1); // Update receptor occupancy times Count_PI3K_time(s_Rec); Count_PI3K_both(m0,(Rec_bound+i_Rec)->PI3K[0]+(Rec_bound+i_Rec)->PI3K[1]); Count_GEF_PI3K(k0,GEF_PI3K_complex(i_Rec)); n_P_tot++; // ******************************** break; case 10: // ************************** Rec-PI3K -> Rec + PI3K ************ k=(int)(N_Rec_PI3K_bound*RAN()); i_Rec=(Rec_PI3K_bound+k)->i; s_Rec=(Rec_PI3K_bound+k)->s; // Start Loop IIIA i_Ras=(Rec_bound+i_Rec)->PI3K_Ras[s_Rec]; if(i_Ras!=-1) Loop_IIIA[i_Ras]=TIME; m0=(Rec_bound+i_Rec)->PI3K[0]+(Rec_bound+i_Rec)->PI3K[1]; k0=GEF_PI3K_complex(i_Rec); Add_Rec_site(Rec_PI3K_free,&N_Rec_PI3K_free,i_Rec,s_Rec,p_Rec_PI3K_free); Remove_Rec_site(Rec_PI3K_bound,&N_Rec_PI3K_bound,k,p_Rec_PI3K_bound); Update_Rec_molecules_2D_release(i_Rec,s_Rec,2,dt); Update_Rec_molecules(i_Rec,2,s_Rec,0,-1); // Update receptor occupancy times Count_PI3K_time(s_Rec); Count_PI3K_both(m0,(Rec_bound+i_Rec)->PI3K[0]+(Rec_bound+i_Rec)->PI3K[1]); Count_GEF_PI3K(k0,GEF_PI3K_complex(i_Rec)); // ******************************** break; } } void Ras_diffusion(double dt){ double t,dr,dmin,dx,dy,t1,t2,tot,tot_j,tot_j_1,x,Prev_x,Prev_y; int J,i,j,k,m,R; //for(k=0;kRec_bound == 0){ do{ if((Ras_bound+i)->d_minstate == 0 && (Ras_bound+i)->Enz == 0){ if( (Rec_bound+0)->GEF[0]==1 && (Rec_bound+0)->GEF_Ras[0]==-1 ) rate[11]=kappa_RE; if( (Rec_bound+0)->GEF[1]==1 && (Rec_bound+0)->GEF_Ras[1]==-1 ) rate[12]=kappa_RE; } if((Ras_bound+i)->state == 1 && (Ras_bound+i)->Enz == 0){ if( (Rec_bound+0)->PI3K[0]==1 && (Rec_bound+0)->PI3K_Ras[0]==-1 ) rate[13]=kappa_RP; if( (Rec_bound+0)->PI3K[1]==1 && (Rec_bound+0)->PI3K_Ras[1]==-1 ) rate[14]=kappa_RP; } if((Ras_bound+i)->Enz == 1){ if( (Rec_bound+0)->GEF[0]==0 ) rate[15]=kappa_RE; if( (Rec_bound+0)->GEF[1]==0 ) rate[16]=kappa_RE; } if((Ras_bound+i)->Enz == 2){ if( (Rec_bound+0)->PI3K[0]==0 ) rate[17]=kappa_RP; if( (Rec_bound+0)->PI3K[1]==0 ) rate[18]=kappa_RP; } tot=0.; for(j=11;j<=18;j++) tot+=rate[j]; R=0; tot_j=0.; tot_j_1=0.; if(tot>0.){ x=tot*RAN(); for(j=11;j<=18;j++){ tot_j+=rate[j]; if(tot_j_1x; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&dx,&dy); PeriodicBC(i); Min_Dist((Ras_bound+i)->x, (Ras_bound+i)->y, &m, &dr); if(dr<=0.){ dr=RAN()*dS; Circle_Point_Picking(S+dr,&dx,&dy); (Ras_bound+i)->x=(Rec_bound+0)->x+dx; (Ras_bound+i)->y=(Rec_bound+0)->y+dy; } (Ras_bound+i)->d_min=dr; (Ras_bound+i)->Rec_closest=0; } if(R==1) t=1.e10; else t+=TSTEP; }else{ t1=dt-t; t2=((Ras_bound+i)->d_min)*((Ras_bound+i)->d_min)/(4.*Diff); if(t1x = (Ras_bound+i)->x + dx; (Ras_bound+i)->y = (Ras_bound+i)->y + dy; PeriodicBC(i); Min_Dist((Ras_bound+i)->x, (Ras_bound+i)->y, &m, &dmin); if(dmin<=0.){ dmin=RAN()*dS; Circle_Point_Picking(S+dmin,&dx,&dy); (Ras_bound+i)->x=(Rec_bound+0)->x+dx; (Ras_bound+i)->y=(Rec_bound+0)->y+dy; } (Ras_bound+i)->d_min=dmin; (Ras_bound+i)->Rec_closest=0; } }while(tx; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&z0,&z); P=1.-PDF_React(z,z0,K_RE_S)/PDF_Reflect(z,z0); if(P>RAN()){ (Rec_bound+m)->GEF_Ras[0]=i; Add_Ras_site(Ras_GDP_GEF,&N_Ras_GDP_GEF,i,p_Ras_GDP_GEF); k=p_Ras_GDP_free[i]; Remove_Ras_site(Ras_GDP_free,&N_Ras_GDP_free,k,p_Ras_GDP_free); (Ras_bound+i)->Enz=1; (Ras_bound+i)->Rec_bound=1; accept=1; /*k=p_Ras_GDP_GEF[i]; Count_inactive_particles(); Nact++; if( (Ras_bound+i)->Rec_bound==1 ) Nact_rec++; Add_Ras_site(Ras_GTP_free,&N_Ras_GTP_free,i,p_Ras_GTP_free); Remove_Ras_site(Ras_GDP_GEF,&N_Ras_GDP_GEF,k,p_Ras_GDP_GEF); Update_Ras_molecules_2D_release(i,1,TSTEP); Update_Ras_molecules(i,1,0,0);*/ } break; case 12: // ************************** Rec-GEF[1] + Ras ***************** m=0; Prev_x=(Ras_bound+i)->x; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&z0,&z); P=1.-PDF_React(z,z0,K_RE_S)/PDF_Reflect(z,z0); if(P>RAN()){ (Rec_bound+m)->GEF_Ras[1]=i; Add_Ras_site(Ras_GDP_GEF,&N_Ras_GDP_GEF,i,p_Ras_GDP_GEF); k=p_Ras_GDP_free[i]; Remove_Ras_site(Ras_GDP_free,&N_Ras_GDP_free,k,p_Ras_GDP_free); (Ras_bound+i)->Enz=1; (Ras_bound+i)->Rec_bound=1; accept=1; } break; case 13: // ************************** Rec-PI3K[0] + Ras **************** m=0; Prev_x=(Ras_bound+i)->x; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&z0,&z); P=1.-PDF_React(z,z0,K_RP_S)/PDF_Reflect(z,z0); if(P>RAN()){ (Rec_bound+m)->PI3K_Ras[0]=i; Add_Ras_site(Ras_GTP_PI3K,&N_Ras_GTP_PI3K,i,p_Ras_GTP_PI3K); k=p_Ras_GTP_free[i]; Remove_Ras_site(Ras_GTP_free,&N_Ras_GTP_free,k,p_Ras_GTP_free); (Ras_bound+i)->Enz=2; (Ras_bound+i)->Rec_bound=1; n_RP_S++; // Destroy Loop I if(Loop_I[i]!=0.) Loop_I[i]=0.; // Count Loop II if(Loop_II[i]!=0.){ N_Loop_II++; T_Loop_II+=TIME-Loop_II[i]; Loop_II[i]=0.; } // Count Loop IIIB if(Loop_IIIB[i]!=0.){ N_Loop_IIIB++; T_Loop_IIIB=TIME-Loop_IIIB[i]; Loop_IIIB[i]=0.; } // Destroy Loop IV if(Loop_IV[i]!=0.) Loop_IV[i]=0.; // Count Loop V if(Loop_V[i]!=0.){ N_Loop_V++; T_Loop_V=TIME-Loop_V[i]; Loop_V[i]=0.; } accept=1; } break; case 14: // ************************** Rec-PI3K[1] + Ras **************** m=0; Prev_x=(Ras_bound+i)->x; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&z0,&z); P=1.-PDF_React(z,z0,K_RP_S)/PDF_Reflect(z,z0); if(P>RAN()){ (Rec_bound+m)->PI3K_Ras[1]=i; Add_Ras_site(Ras_GTP_PI3K,&N_Ras_GTP_PI3K,i,p_Ras_GTP_PI3K); k=p_Ras_GTP_free[i]; Remove_Ras_site(Ras_GTP_free,&N_Ras_GTP_free,k,p_Ras_GTP_free); (Ras_bound+i)->Enz=2; (Ras_bound+i)->Rec_bound=1; n_RP_S++; // Destroy Loop I if(Loop_I[i]!=0.) Loop_I[i]=0.; // Count Loop II if(Loop_II[i]!=0.){ N_Loop_II++; T_Loop_II+=TIME-Loop_II[i]; Loop_II[i]=0.; } // Count Loop IIIB if(Loop_IIIB[i]!=0.){ N_Loop_IIIB++; T_Loop_IIIB=TIME-Loop_IIIB[i]; Loop_IIIB[i]=0.; } // Destroy Loop IV if(Loop_IV[i]!=0.) Loop_IV[i]=0.; // Count Loop V if(Loop_V[i]!=0.){ N_Loop_V++; T_Loop_V=TIME-Loop_V[i]; Loop_V[i]=0.; } accept=1; } break; case 15: // ************************** Rec[0] + GEF-Ras ***************** m=0; Prev_x=(Ras_bound+i)->x; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&z0,&z); P=1.-PDF_React(z,z0,K_R_ES)/PDF_Reflect(z,z0); if(P>RAN()){ m0=(Rec_bound+m)->GEF[0]+(Rec_bound+m)->GEF[1]; k0=GEF_PI3K_complex(m); (Rec_bound+m)->GEF[0]=1; (Rec_bound+m)->GEF_Ras[0]=i; Add_Rec_site(Rec_GEF_bound,&N_Rec_GEF_bound,m,0,p_Rec_GEF_bound); k=p_Rec_GEF_free[m][0]; Remove_Rec_site(Rec_GEF_free,&N_Rec_GEF_free,k,p_Rec_GEF_free); Count_GEF_time(0); n_ES++; n_E_tot++; (Ras_bound+i)->Rec_bound=1; Count_GEF_both(m0,(Rec_bound+m)->GEF[0]+(Rec_bound+m)->GEF[1]); Count_GEF_PI3K(k0,GEF_PI3K_complex(m)); accept=1; } break; case 16: // ************************** Rec[1] + GEF-Ras ***************** m=0; Prev_x=(Ras_bound+i)->x; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&z0,&z); P=1.-PDF_React(z,z0,K_R_ES)/PDF_Reflect(z,z0); if(P>RAN()){ m0=(Rec_bound+m)->GEF[0]+(Rec_bound+m)->GEF[1]; k0=GEF_PI3K_complex(m); (Rec_bound+m)->GEF[1]=1; (Rec_bound+m)->GEF_Ras[1]=i; Add_Rec_site(Rec_GEF_bound,&N_Rec_GEF_bound,m,1,p_Rec_GEF_bound); k=p_Rec_GEF_free[m][1]; Remove_Rec_site(Rec_GEF_free,&N_Rec_GEF_free,k,p_Rec_GEF_free); Count_GEF_time(1); n_ES++; n_E_tot++; (Ras_bound+i)->Rec_bound=1; Count_GEF_both(m0,(Rec_bound+m)->GEF[0]+(Rec_bound+m)->GEF[1]); Count_GEF_PI3K(k0,GEF_PI3K_complex(m)); accept=1; } break; case 17: // ************************** Rec[0] + PI3K-Ras **************** m=0; Prev_x=(Ras_bound+i)->x; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&z0,&z); P=1.-PDF_React(z,z0,K_R_PS)/PDF_Reflect(z,z0); if(P>RAN()){ m0=(Rec_bound+m)->PI3K[0]+(Rec_bound+m)->PI3K[1]; k0=GEF_PI3K_complex(m); (Rec_bound+m)->PI3K[0]=1; (Rec_bound+m)->PI3K_Ras[0]=i; Add_Rec_site(Rec_PI3K_bound,&N_Rec_PI3K_bound,m,0,p_Rec_PI3K_bound); k=p_Rec_PI3K_free[m][0]; Remove_Rec_site(Rec_PI3K_free,&N_Rec_PI3K_free,k,p_Rec_PI3K_free); Count_PI3K_time(0); n_PS++; n_P_tot++; (Ras_bound+i)->Rec_bound=1; Count_PI3K_both(m0,(Rec_bound+m)->PI3K[0]+(Rec_bound+m)->PI3K[1]); Count_GEF_PI3K(k0,GEF_PI3K_complex(m)); // Count Loop I if(Loop_I[i]!=0.){ N_Loop_I++; T_Loop_I+=TIME-Loop_I[i]; Loop_I[i]=0.; } // Destroy Loop II if(Loop_II[i]!=0.) Loop_II[i]=0.; // Count Loop IIIA if(Loop_IIIA[i]!=0.){ N_Loop_IIIA++; T_Loop_IIIA+=TIME-Loop_IIIA[i]; Loop_IIIA[i]=0.; } // Destroy Loop IIIB if(Loop_IIIB[i]!=0.) Loop_IIIB[i]=0.; // Count Loop IV if(Loop_IV[i]!=0.){ N_Loop_IV++; T_Loop_IV+=TIME-Loop_IV[i]; Loop_IV[i]=0.; } // Destroy Loop V if(Loop_V[i]!=0.) Loop_V[i]=0.; accept=1; } break; case 18: // ************************** Rec[1] + PI3K-Ras **************** m=0; Prev_x=(Ras_bound+i)->x; Prev_y=(Ras_bound+i)->y; Gauss_propagation((Rec_bound+m)->x,(Rec_bound+m)->y,Prev_x,Prev_y, TSTEP,&(Ras_bound+i)->x,&(Ras_bound+i)->y,&z0,&z); P=1.-PDF_React(z,z0,K_R_PS)/PDF_Reflect(z,z0); if(P>RAN()){ m0=(Rec_bound+m)->PI3K[0]+(Rec_bound+m)->PI3K[1]; k0=GEF_PI3K_complex(m); (Rec_bound+m)->PI3K[1]=1; (Rec_bound+m)->PI3K_Ras[1]=i; Add_Rec_site(Rec_PI3K_bound,&N_Rec_PI3K_bound,m,1,p_Rec_PI3K_bound); k=p_Rec_PI3K_free[m][1]; Remove_Rec_site(Rec_PI3K_free,&N_Rec_PI3K_free,k,p_Rec_PI3K_free); Count_PI3K_time(1); n_PS++; n_P_tot++; (Ras_bound+i)->Rec_bound=1; Count_PI3K_both(m0,(Rec_bound+m)->PI3K[0]+(Rec_bound+m)->PI3K[1]); Count_GEF_PI3K(k0,GEF_PI3K_complex(m)); // Count Loop I if(Loop_I[i]!=0.){ N_Loop_I++; T_Loop_I+=TIME-Loop_I[i]; Loop_I[i]=0.; } // Destroy Loop II if(Loop_II[i]!=0.) Loop_II[i]=0.; // Count Loop IIIA if(Loop_IIIA[i]!=0.){ N_Loop_IIIA++; T_Loop_IIIA+=TIME-Loop_IIIA[i]; Loop_IIIA[i]=0.; } // Destroy Loop IIIB if(Loop_IIIB[i]!=0.) Loop_IIIB[i]=0.; // Count Loop IV if(Loop_IV[i]!=0.){ N_Loop_IV++; T_Loop_IV+=TIME-Loop_IV[i]; Loop_IV[i]=0.; } // Destroy Loop V if(Loop_V[i]!=0.) Loop_V[i]=0.; accept=1; } break; } PeriodicBC(i); Min_Dist((Ras_bound+i)->x, (Ras_bound+i)->y, &m, &dr); if(dr<=0.){ dr=RAN()*dS; Circle_Point_Picking(S+dr,&dx,&dy); (Ras_bound+i)->x=(Rec_bound+0)->x+dx; (Ras_bound+i)->y=(Rec_bound+0)->y+dy; } (Ras_bound+i)->d_min=dr; (Ras_bound+i)->Rec_closest=0; return accept; } double Diffusion_time_step(){ double x=1.e10,dt; int i; for(i=0;id_mind_min; } if(x>dS) dt=x*x/(4.*Diff); else dt=TSTEP; return dt; } // ************************************************************************************************** // ************************************************************************************************** // ************************************************************************************************** void Output(){ double T_E_av, T_P_av, m_E, m_P, T_E_2, T_P_2, T_EP_2, a, a_rec,Nin_av; T_E_av=0.5*(T_E[0]+T_E[1])/TIME; T_P_av=0.5*(T_P[0]+T_P[1])/TIME; if(n_E_tot==0) m_E=0.; else m_E=(float)(n_ES)/n_E_tot; if(n_P_tot==0) m_P=0.; else m_P=(float)(n_PS)/n_P_tot; T_E_2=TT_E/TIME; T_P_2=TT_P/TIME; T_EP_2=TT_EP/TIME; Nin_av=(float)(Nin)/Nin_times; a=1.*Nact*Lxy*Lxy/(TIME*Nin_av*Diff); a_rec=1.*Nact_rec*Lxy*Lxy/(TIME*Nin_av*Diff); fprintf(OutputFile,"\n %f %f %f %f %f %f %f %f %f %f", TIME,T_E_av,T_P_av,m_E,m_P,T_E_2,T_P_2,T_EP_2,a,a_rec); fprintf(alterm,"\n %f %f %f",TIME,n_ES/TIME,n_PS/TIME); // Note: n's are average numbers of molecules; r = receptor, b = bulk printf("\n t=%f, T_E=%f, T_P=%f, m_E=%f m_P=%f n_GTP_free=%f n_SP_r=%f n_SP_b=%f a=%f a_r=%f", TIME,T_E_av,T_P_av,m_E,m_P,gtp_free/n_count,(pi3k-pi3k_bulk)/n_count, pi3k_bulk/n_count,a,a_rec); //printf("\n t=%f, T_E=%f, T_P=%f, m_E=%f m_P=%f n_GDP_free=%f n_SE_r=%f n_SE_b=%f a=%f a_r=%f", // TIME,T_E_av,T_P_av,m_E,m_P,gdp_free/n_count,(gef-gef_bulk)/n_count, // gef_bulk/n_count,a,a_rec); //printf("\n t=%f, T_P[0]=%f, T_P[1]=%f, m_P=%f n_GTP_free=%f n_SP_r=%f n_SP_b=%f a=%f a_r=%f", // TIME,T_P[0]/TIME,T_P[1]/TIME,m_P,gtp_free/n_count,(pi3k-pi3k_bulk)/n_count, // pi3k_bulk/n_count,a,a_rec); fprintf(OutputLoop,"\n %f %e %e %e %e %e %e %e %e %e", TIME, n_P_tot/TIME, n_PS/TIME, n_RP_S/TIME, N_Loop_I/TIME, N_Loop_II/TIME, N_Loop_IIIA/TIME, N_Loop_IIIB/TIME, N_Loop_IV/TIME, N_Loop_V/TIME); printf("\n t=%f n_P_tot=%e n_RP_S=%e I=%e II=%e IIIA=%e IIIB=%e IV=%e V=%e", TIME, n_P_tot/TIME, n_PS/TIME, n_RP_S/TIME, N_Loop_I/TIME, N_Loop_II/TIME, N_Loop_IIIA/TIME,N_Loop_IIIB/TIME, N_Loop_IV/TIME, N_Loop_V/TIME); PrintInDist(); } void PrintInDist(){ FILE *hadpasa; int i; hadpasa=fopen("inact.int","w"); for(i=0;i<(int)(Lxy);i++) fprintf(hadpasa,"%d\n",InDist[i]); fclose(hadpasa); } void Count_inactive_particles(){ int i; for(i=0;istate == 0 ) Nin++; } Nin_times++; } void Ras_GTP_and_PI3K(){ int i,n=0; gtp_free += 1.*N_Ras_GTP_free; pi3k += 1.*N_Ras_GTP_PI3K; for(i=0;iEnz==2 && (Ras_bound+i)->Rec_bound==0) n++; } pi3k_bulk += 1.*n; n_count += 1.; } void Ras_GDP_and_GEF(){ int i,n=0; gdp_free += 1.*N_Ras_GDP_free; gef += 1.*N_Ras_GDP_GEF; for(i=0;iEnz==1 && (Ras_bound+i)->Rec_bound==0) n++; } gef_bulk += 1.*n; n_count += 1.; } void Read_Parameters(char datafile[]){ char stam[20]; FILE *InputFile; InputFile=fopen(datafile,"r"); fscanf(InputFile,"%s %d", &stam,&Rec_n); fscanf(InputFile,"%s %d", &stam,&Ras_n); fscanf(InputFile,"%s %lf",&stam,&S); fscanf(InputFile,"%s %lf",&stam,&dS); fscanf(InputFile,"%s %lf",&stam,&eta_R); fscanf(InputFile,"%s %lf",&stam,&TSTEP); fscanf(InputFile,"%s %lf",&stam,&Diff); fscanf(InputFile,"%s %lf",&stam,&phi_RE); fscanf(InputFile,"%s %lf",&stam,&tau_RE); fscanf(InputFile,"%s %lf",&stam,&tau_SE); fscanf(InputFile,"%s %lf",&stam,&kci_E); fscanf(InputFile,"%s %lf",&stam,&kappa_RE); fscanf(InputFile,"%s %lf",&stam,&phi_act); fscanf(InputFile,"%s %lf",&stam,&tau_GAP); fscanf(InputFile,"%s %lf",&stam,&phi_RP); fscanf(InputFile,"%s %lf",&stam,&tau_RP); fscanf(InputFile,"%s %lf",&stam,&tau_SP); fscanf(InputFile,"%s %lf",&stam,&kci_P); fscanf(InputFile,"%s %lf",&stam,&kappa_RP); fscanf(InputFile,"%s %lf",&stam,&t_end); fscanf(InputFile,"%s %lf",&stam,&t_print); fclose(InputFile); GET_DATA(ptr_erf, ptr_erfcx); } void Evaluate_rate_constants(){ double alpha_0=0.325; constant = malloc(20*sizeof(double)); rate = malloc(20*sizeof(double)); constant[0]=0.; rate[0]=0.; K_R_ES = Diff*kappa_RE/(2.*pi*S); K_RE_S = kci_E*Diff*kappa_RE/(2.*pi*S); K_R_PS = Diff*kappa_RP/(2.*pi*S); K_RP_S = kci_P*Diff*kappa_RP/(2.*pi*S); constant[1]=kci_E*Diff*phi_RE/(S*S*tau_RE); // k_on_SE constant[2]=Diff/(1.+phi_act)/(S*S*tau_SE); // k_off_SE constant[3]=Diff*phi_act/(1.+phi_act)/(S*S*tau_SE); // k_act constant[4]=Diff/(S*S*tau_GAP); // k_GAP = k_in constant[5]=kci_P*Diff*phi_RP/(S*S*tau_RP); // k_on_SP constant[6]=Diff/(S*S*tau_SP); // k_off_SP constant[7]=Diff*phi_RE/(S*S*tau_RE); // k_on_RE constant[8]=Diff/(S*S*tau_RE); // k_off_RE constant[9]=Diff*phi_RP/(S*S*tau_RP); // k_on_RP constant[10]=Diff/(S*S*tau_RP); // k_off_RP /*K_R_ES=0.; K_RE_S=0.; constant[1]=(alpha_0*Diff/(Lxy*Lxy))*(constant[2]+constant[3])/constant[3]; constant[7]=0.;*/ /*double Da=0.01; K_R_ES = 0.; K_RE_S = Diff*kappa_RE/(2.*pi*S); K_R_PS = 0.; K_RP_S = 0.; constant[1] =0.; // k_on_SE constant[2] =0.; // k_off_SE constant[3] =0.; // k_act constant[4] =Da*Diff/(S*S); // k_GAP = k_in constant[5] =0.; // k_on_SP constant[6] =0.; // k_off_SP constant[7] =0.; // k_on_RE constant[8] =0.; // k_off_RE constant[9] =0.; // k_on_RP constant[10]=0.; // k_off_RP*/ } void Evaluate_other_parameters(){ SIGMA=sqrt(4.*Diff*TSTEP); mu_plus = 2.*TSTEP/sqrt(4.*Diff*TSTEP); Lxy=S*sqrt(Rec_n/eta_R); r_bimol=1./TSTEP; } void Create_Molecules(){ int i,j,m,count; double dmin,dist,a,x[2],y[2],z[2],ARRAY[10010][2]={0.0}; Rec_bound=malloc( Rec_n*(8*sizeof(int)+2*sizeof(double)) ); for(i=0;iGEF=malloc(2*sizeof(int)); for(j=0;j<2;j++) (Rec_bound+i)->GEF[j]=0; (Rec_bound+i)->PI3K=malloc(2*sizeof(int)); for(j=0;j<2;j++) (Rec_bound+i)->PI3K[j]=0; (Rec_bound+i)->GEF_Ras=malloc(2*sizeof(int)); for(j=0;j<2;j++) (Rec_bound+i)->GEF_Ras[j]=-1; (Rec_bound+i)->PI3K_Ras=malloc(2*sizeof(int)); for(j=0;j<2;j++) (Rec_bound+i)->PI3K_Ras[j]=-1; } if(Rec_n==1){ (Rec_bound+0)->x=0.5*Lxy; (Rec_bound+0)->y=0.5*Lxy; }else{ count=0; a=10.*(S+dS); ARRAY[0][0]=a+(Lxy-2*a)*RAN(); ARRAY[0][1]=a+(Lxy-2*a)*RAN(); z[0]=ARRAY[0][0]; z[1]=ARRAY[0][1]; (Rec_bound+0)->x=z[0]; (Rec_bound+0)->y=z[1]; while(count=(2*S)*(2*S)){ count++; ARRAY[count][0]=y[0]; ARRAY[count][1]=y[1]; (Rec_bound+count)->x=y[0]; (Rec_bound+count)->y=y[1]; } } } Ras_bound=malloc( Ras_n*(4*sizeof(int)+3*sizeof(double)) ); for(i=0;iEnz=0; (Ras_bound+i)->state=0; (Ras_bound+i)->Rec_bound=0; do{ (Ras_bound+i)->x=Lxy*RAN(); (Ras_bound+i)->y=Lxy*RAN(); dmin=1.e10; for(j=0;jx - (Rec_bound+j)->x)*((Ras_bound+i)->x - (Rec_bound+j)->x)+ ((Ras_bound+i)->y - (Rec_bound+j)->y)*((Ras_bound+i)->y - (Rec_bound+j)->y)) - S; if(dmin>dist){ dmin=dist; m=j; } } }while(dmin <= 0.); (Ras_bound+i)->d_min=dmin; (Ras_bound+i)->Rec_closest=m; } Ras_GDP_free = malloc(Ras_n*sizeof(int)); p_Ras_GDP_free = malloc(Ras_n*sizeof(int)); N_Ras_GDP_free = 0; for(i=0;istate = st; (Ras_bound+i)->Enz = enz; (Ras_bound+i)->Rec_bound = rec; } void Update_Ras_molecules_2D_release(int i_ras, int Enz, double DT){ int i_rec; // Rec_closest: index of the closest receptor (the same) // x,y: new coordinates // d_min; new distance to the closest receptor // Rec_bound check before the update // Update receptor-bound enzymes if( (Ras_bound+i_ras)->Rec_bound==1 ){ i_rec=(Ras_bound+i_ras)->Rec_closest; Release_Ras(i_ras,i_rec, DT); switch (Enz){ case 1: if((Rec_bound+i_rec)->GEF_Ras[0] == i_ras ) (Rec_bound+i_rec)->GEF_Ras[0]=-1; else (Rec_bound+i_rec)->GEF_Ras[1]=-1; break; case 2: if((Rec_bound+i_rec)->PI3K_Ras[0] == i_ras ) (Rec_bound+i_rec)->PI3K_Ras[0]=-1; else (Rec_bound+i_rec)->PI3K_Ras[1]=-1; break; } } } void Update_Rec_molecules(int i, int domain, int site, int bound, int ras){ // domain: 1 = for GEF; 2 = for PI3K // ENZ[]: 0 = not bound; 1 = bound // PI3K_Ras[]: -1 = not bound; other = index of Ras bound if(domain==1){ (Rec_bound+i)->GEF[site] = bound; (Rec_bound+i)->GEF_Ras[site] = ras; }else{ (Rec_bound+i)->PI3K[site] = bound; (Rec_bound+i)->PI3K_Ras[site] = ras; } } void Update_Rec_molecules_2D_release(int i_rec, int s_rec, int Enz, double DT){ int i_ras; // x,y: new coordinates for Ras // d_min; new distance to the closest receptor for Ras // Rec_bound update for Ras // Update receptor-bound Ras switch (Enz){ case 1: i_ras=(Rec_bound+i_rec)->GEF_Ras[s_rec]; if(i_ras!=-1){ (Rec_bound+i_rec)->GEF_Ras[s_rec]=-1; (Ras_bound+i_ras)->Rec_bound=0; Release_Ras(i_ras,i_rec,DT); } break; case 2: i_ras=(Rec_bound+i_rec)->PI3K_Ras[s_rec]; if(i_ras!=-1){ (Rec_bound+i_rec)->PI3K_Ras[s_rec]=-1; (Ras_bound+i_ras)->Rec_bound=0; Release_Ras(i_ras,i_rec,DT); } break; } } void Release_Ras(int i_ras, int i_rec, double DT){ double dmin,dx,dy; dmin=RAN()*dS; Circle_Point_Picking(S+dmin,&dx,&dy); (Ras_bound+i_ras)->x=(Rec_bound+i_rec)->x+dx; (Ras_bound+i_ras)->y=(Rec_bound+i_rec)->y+dy; (Ras_bound+i_ras)->d_min=dmin; } // **** COUNT RECEPTOR OCCUPANCY TIMES AND Ras-BINDING EVENTS ****** void Count_PI3K_time(int s){ if( P_bound[s]==0 ){ t_PI3K_0[s]=TIME; P_bound[s]=1; }else{ T_P[s]+=TIME-t_PI3K_0[s]; P_bound[s]=0; } } void Count_GEF_time(int s){ if( E_bound[s]==0 ){ t_GEF_0[s]=TIME; E_bound[s]=1; }else{ T_E[s]+=TIME-t_GEF_0[s]; E_bound[s]=0; } } void Count_GEF_both(int m0, int m1){ if( m0==0 && m1>0 ){ tt_GEF_0=TIME; }else if( m0>0 && m1==0 ){ TT_E+=TIME-tt_GEF_0; } } void Count_PI3K_both(int m0, int m1){ if( m0==0 && m1>0 ){ tt_PI3K_0=TIME; }else if( m0>0 && m1==0 ){ TT_P+=TIME-tt_PI3K_0; } } int GEF_PI3K_complex(int i){ if( (Rec_bound+i)->GEF[0] +(Rec_bound+i)->GEF[1]>0 && (Rec_bound+i)->PI3K[0]+(Rec_bound+i)->PI3K[1]>0 ) return 1; else return 0; } void Count_GEF_PI3K(int m0, int m1){ if( m0==0 && m1>0 ){ tt_EP_0=TIME; }else if( m0>0 && m1==0 ){ TT_EP+=TIME-tt_EP_0; } } void Add_Ras_site(int *list, int *Length, int a, int *pos){ *(list+(*Length)) = a; *(pos+a) = *Length; (*Length)++; } void Remove_Ras_site(int *list, int *Length, int M, int *pos){ if((*Length)==0 || M>(*Length)-1){ printf("\nTrying to remove non-existing element!!!\n"); printf("Failed in Remove_Ras_site()\n"); exit(0); } *(pos+(*(list+M))) = -1; if(M!=(*Length)-1){ *(list+M) = *(list+(*Length)-1); *(pos+(*(list+M))) = M; } (*Length)--; } void Add_Rec_site(struct List_of_2 *list, int *Length, int a, int b, int **pos){ (list+(*Length))->i = a; (list+(*Length))->s = b; *(*(pos+a)+b)=*Length; (*Length)++; } void Remove_Rec_site(struct List_of_2 *list, int *Length, int M, int **pos){ if((*Length)==0 || M>(*Length)-1){ printf("\nTrying to remove non-existing element!!!\n"); printf("Failed in Remove_Rec_site()\n"); exit(0); } *(*(pos+(list+M)->i)+(list+M)->s)=-1; if(M!=(*Length)-1){ (list+M)->i = (list+(*Length)-1)->i; (list+M)->s = (list+(*Length)-1)->s; *(*(pos+(list+M)->i)+(list+M)->s) = M; } (*Length)--; } /*Long period (> 2x10^18) random number generator of L'Ecuyer with Bays-Durham shuffle and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint values). Call with idum a negative integer to initialize; thereafter, do not alter idum between successive deviates in a sequence. RNMX should approximate the largest floating value that is less than 1.*/ #define IM1 2147483563 #define IM2 2147483399 #define AM (1.0/IM1) #define IMM1 (IM1-1) #define IA1 40014 #define IA2 40692 #define IQ1 53668 #define IQ2 52774 #define IR1 12211 #define IR2 3791 #define NTAB 32 #define NDIV (1+IMM1/NTAB) #define EPS 1.2e-7 #define RNMX (1.0-EPS) long idum=-536445; double RAN(){ return drand48(); int j; long k; static long idum2=123456789; static long iy=0; static long iv[NTAB]; float temp; if (idum <= 0){ // Initialize. if (-(idum) < 1) idum=1; // Be sure to prevent idum = 0. else idum = -(idum); idum2=(idum); for(j=NTAB+7;j>=0;j--){ // Load the shuffle table (after 8 warm-ups). k=(idum)/IQ1; idum=IA1*(idum-k*IQ1)-k*IR1; if (idum < 0) idum += IM1; if (j < NTAB) iv[j] = idum; } iy=iv[0]; } k=(idum)/IQ1; // Start here when not initializing. idum=IA1*(idum-k*IQ1)-k*IR1; // Compute idum=(IA1*idum) % IM1 without if (idum < 0) idum += IM1; // overflows by Schrage's method. k=idum2/IQ2; idum2=IA2*(idum2-k*IQ2)-k*IR2; // Compute idum2=(IA2*idum) % IM2 likewise. if (idum2 < 0) idum2 += IM2; j=iy/NDIV; // Will be in the range 0..NTAB-1. iy=iv[j]-idum2; // Here idum is shuffled, idum and idum2 are iv[j] = idum; // combined to generate output. if (iy < 1) iy += IMM1; if ((temp=AM*iy) > RNMX) return RNMX; // Because users don't expect endpoint values. else return temp; } double PDF_React(double x, double x0, double K){ double z,z0,f1,f2; z = x/SIGMA; z0 = x0/SIGMA; f1 = ( exp(-(z-z0)*(z-z0)) + exp(-(z+z0)*(z+z0)) )/(SIGMA*sqrtPI); f2 = (K/Diff)*exp(-(z+z0)*(z+z0))*EVALUATE((z+z0+K*mu_plus),ptr_erfcx,0); return f1 - f2; } double PDF_Reflect(double x, double x0){ double z,z0; z = x/SIGMA; z0 = x0/SIGMA; return ( exp(-(z-z0)*(z-z0)) + exp(-(z+z0)*(z+z0)) )/(SIGMA*sqrtPI); } void Circle_Point_Picking(double r, double *xr, double *yr){ double x1,x2,x12,x22,sum; do{ x1=(1.-2.*RAN()); x12=x1*x1; x2=(1.-2.*RAN()); x22=x2*x2; sum=x12+x22; }while(sum>=1.); *xr=r*(x12-x22)/sum; *yr=r*2.*x1*x2/sum; } void PeriodicBC(int J){ double X,Y; X=(Ras_bound+J)->x; Y=(Ras_bound+J)->y; if(X>0. && XLxy ) (Ras_bound+J)->x = (double) fmod(X,Lxy); else (Ras_bound+J)->x = Lxy + (double) fmod(X,Lxy); } if(Y>0. && YLxy ) (Ras_bound+J)->y = (double) fmod(Y,Lxy); else (Ras_bound+J)->y = Lxy + (double) fmod(Y,Lxy); } } void Gauss_propagation(double Xr, double Yr, double Xi, double Yi, double DT, double *Xnew, double *Ynew, double *DR_curr, double *DR_new){ double ksi,eta,rho,alpha,beta,Ra; Ra=sqrt((Xi-Xr)*(Xi-Xr)+(Yi-Yr)*(Yi-Yr)); *DR_curr = Ra - S; ksi=fabs( *DR_curr+GAUSS()*sqrt(2.*DT*Diff) ); eta=GAUSS()*sqrt(2.*DT*Diff); rho=sqrt(ksi*ksi+eta*eta); alpha=Rotation_angle(ksi,eta); beta=Rotation_angle((Xi-Xr),(Yi-Yr)); *Xnew = Xr + S*cos(beta) + rho*cos(alpha+beta); *Ynew = Yr + S*sin(beta) + rho*sin(alpha+beta); *DR_new=sqrt((*Xnew-Xr)*(*Xnew-Xr)+(*Ynew-Yr)*(*Ynew-Yr)) - S; } double GAUSS(){ static int iset = 0; static double gset; double fac, v1, v2, rsquare=2.0; if(iset == 0) { while(rsquare>1.0 || rsquare==0.0) { v1 = 2*RAN()-1.0; v2 = 2*RAN()-1.0; rsquare = v1*v1 + v2*v2; }; fac=sqrt(-2.0*log(rsquare)/rsquare); gset=v1*fac; iset=1; return v2*fac;} else{ iset=0; return gset; } } double Rotation_angle(double DX, double DY){ double angle; if(fabs(DX)<1.e-10){ if(DY>0.) angle=pi/2.; else angle=pi+pi/2.; } else{ if(DX>0.){ if(DY>=0.) angle=atan(DY/DX); else angle=2.*pi-atan(-DY/DX); }else{ if(DY>0.) angle=pi-atan(DY/(-DX)); else angle=pi+atan(DY/DX); } } return angle; } void Min_Dist(double X, double Y, int *M, double *Dmin){ int i; double d; *Dmin=1.e20; for(i=0;ix)*(X - (Rec_bound+i)->x) + (Y - (Rec_bound+i)->y)*(Y - (Rec_bound+i)->y)) - S; if(*Dmin>d){ *Dmin=d; *M=i; } } } double MIN_DIST(double x[], double y[], double L){ int count=0; double stam, d[2]; for(count=0;count<2;count++){ stam=fabs(x[count]-y[count]); if(stam>L/2) stam=L-stam; d[count]=stam; } return d[0]*d[0]+d[1]*d[1]; } void GET_DATA(LookUpTable *erf, LookUpTable *erfcx){ int I, index; FILE *ifp; double a, b, c; assert( (ifp=fopen("ERFCX.dat","r")) != NULL); fscanf(ifp,"%lf%lf", &(erfcx -> XMIN), &(erfcx -> XMAX)); fscanf(ifp,"%lf%lf", &(erfcx -> DX), &(erfcx -> DX)); a = (erfcx -> XMAX); b = (erfcx -> XMIN); c = (erfcx -> DX); I = (int) ((a-b)/c); for(index=0; index<=I; index++) fscanf(ifp, "%lf%lf", &((*erfcx).Table)[index][0], &((*erfcx).Table)[index][1]); fclose(ifp); assert( (ifp=fopen("ERF.dat","r")) != NULL); fscanf(ifp,"%lf%lf", &(erf -> XMIN), &(erf -> XMAX)); fscanf(ifp,"%lf%lf", &(erf -> DX), &(erf -> DX)); a = (erf -> XMAX); b = (erf -> XMIN); c = (erf -> DX); I = (int) ((a-b)/c); for(index=0; index<=I; index++) fscanf(ifp, "%lf%lf", &((*erf).Table)[index][0], &((*erf).Table)[index][1]); fclose(ifp); return; } double EVALUATE(double x, LookUpTable *POINT, int MODE){ double k, b; int IND; assert( x>=(*POINT).XMIN ); if( x<=(*POINT).XMAX ) { IND = (int) (floor(x/(*POINT).DX)); k = (POINT -> Table)[IND][0]; b = (POINT -> Table)[IND][1]; return (k*x + b); } else { if(MODE) return 1.0; else return 1/(sqrtPI*x); } } void Test_Ras_list(int n, int *list, int *pos){ int i; printf("\n"); for(i=0;ii,(list+i)->s); printf("\n Positions:\n"); for(i=0;iEnz,(Ras_bound+i)->state, (Ras_bound+i)->Rec_bound); printf(" x=%f y=%f Rec_closest=%d d_min=%f\n",(Ras_bound+i)->x,(Ras_bound+i)->y, (Ras_bound+i)->Rec_closest,(Ras_bound+i)->d_min); //} printf("\n"); } void Test_Rec_mol(int i){ printf("\n"); //for(i=0;iGEF[0],(Rec_bound+i)->GEF[1],(Rec_bound+i)->PI3K[0],(Rec_bound+i)->PI3K[1]); printf(" GEF_Ras[0]=%d GEF_Ras[1]=%d PI3K_Ras[0]=%d PI3K_Ras[1]=%d\n", (Rec_bound+i)->GEF_Ras[0],(Rec_bound+i)->GEF_Ras[1], (Rec_bound+i)->PI3K_Ras[0],(Rec_bound+i)->PI3K_Ras[1]); //} printf("\n"); }