DNA.cpp Source File

00001 // DNA.cpp: implementation of the CDNA class.
00002 //
00004 
00005 #include "Cell.h"
00006 #include "Display.h"
00007 #include "Geometry.h"
00008 #include "Components.h"
00009 #include "DNA.h"
00010 #include <math.h>
00011 
00013 // Construction/Destruction
00015 
00016 CDNA::CDNA(CCell* pCell)
00017 {
00018         m_pCell=pCell;
00019         GD=12.0;                                //Gene dosage
00020         NPOS=1;                 //Number of fork positions
00021         DnaABoxes[0] = 9;
00022         DnaABoxes[1] = 154;
00023         DnaABoxes[2] = 804132;          //number of DnaA boxes
00024         NumForks[0] = 1;
00025         NumForks[1]=0;
00026         NumForks[2]=0;                  //Number of replicating forks in DNA (formerly NF)
00027         ForkPos[0] = 0.56;
00028         ForkPos[1] = 0.0;
00029         ForkPos[2] = 0.0;       //Fork position (% total DNA length, formerly FPOS)
00030         TotalGD = NumForks[0]*4405*(1+ForkPos[0]+2*ForkPos[1]+4*ForkPos[2]);
00031         NumChrome=1;                    //Number of chromosomes
00032         DnaA = 1000/6.023E23*52500;             //DnaA weight
00033         Kd[0] = 0.79;
00034         Kd[1] = 33.0;
00035         Kd[2] = 2.0E6;  //dissociation constants (nM)
00036         MidBound = 0.88;                                        //initial bound fraction
00037         NTOT=1;
00038         NORG=2;
00039         NFTERM=0;
00040 }
00041 
00042 CDNA::CDNA()
00043 {
00044 
00045 }
00046 
00047 CDNA::~CDNA()
00048 {
00049 
00050 }
00051 
00052 
00053 CDNA::CDNA(const CDNA &rDNA, CCell *pCell)
00054 {
00055         //*this = *rDNA;
00056         m_pCell=pCell;
00057 }
00058 
00059 //INITIATION
00060 void CDNA::EcoInit()
00061 {
00062         double CPDnaA=20*0.995*1.5;             //Rate of DnaA production
00063         double DT=m_pCell->GetComponents()->DT;
00064         int BoxNum = 5;                 //Boxes needed to bind for initiation
00065 
00066         
00067 
00068         double Criterion=0.901398215;                   //Criterion for initiation
00069         double MethTime=6;                                              //Time for methylation
00070         double FreeDnaA = 0;
00071         double BoundFrac[3] = {0.0,0.0,0.0};    //Bound fractions of boxes
00072         static bool R1R2signal = 0;                             //Signal for complex II
00073         static double MethState = 1.1;                          //Signal for methylation                                                
00074                                                                                 // (greater than 1 if fully methylated)
00075         int NFINT=0;
00076         R1R2signal = 0;                 
00077         //Compute DnaA gene dosage
00078         double DnaAGD=NumChrome;
00079         for(int i=0;i<3;++i)
00080         {
00081                 if(ForkPos[i]>=0.0181)          //DnaA gene is at this position
00082                         DnaAGD*=2;
00083         }
00084         //Compute total gene dosage
00085         //Assume 4405 genes / chromosome
00086         TotalGD = NumChrome*4405*(1+ForkPos[0]+2*ForkPos[1]+4*ForkPos[2]);
00087                         
00088         ComputeBoxes(); //Compute DnaABoxes
00089         FreeDnaA = Equilibrium();               //Compute Equilibrium concentration
00090 
00091         for(i=0;i<3;++i)
00092                 BoundFrac[i] = FreeDnaA/(FreeDnaA + Kd[i]);                     
00093                         
00094         double FracB1 = FreeDnaA/(FreeDnaA + Kd[0]);
00095         double Prob = gamma(FracB1*DnaABoxes[0]+1)*gamma(DnaABoxes[0]-BoxNum+1)
00096                 /gamma(FracB1*DnaABoxes[0]-BoxNum+1)/gamma(DnaABoxes[0]+1);
00097 
00098         MidBound = BoundFrac[1];
00099 
00100         MethState += MethTime*DT;
00101 
00102         if(Prob>0.9&&MethState>1)
00103         {
00104                 R1R2signal = 1;                 //signal complexII complete
00105                 if(BoundFrac[1]>Criterion)              //Initiation!
00106                 {
00107                         NFINT=NORG;
00108                         NTOT+=NFINT;            //Add the new forks
00109                         NORG=NORG+NFINT;        //Twice the origins
00110                         NPOS += 1;                      //One more fork
00111                         NumForks[NPOS-1]=NFINT; //Count the forks on newly divided strand
00112                         MethState = 0;          //Methylated state is now 0
00113                 }
00114         } 
00115         
00116         double DY6S=m_pCell->GetComponents()->GetDY6S();
00117         double Y1=m_pCell->GetComponents()->GetY(1);
00118         double V=m_pCell->GetGeometry()->GetVolume();
00119         //Update DnaA
00120         if(MethState>=1)
00121         {
00122                 DnaA += DT*(1-BoundFrac[1])*CPDnaA*DY6S*(DnaAGD/TotalGD); //Compute DnaA synthesis
00123         }
00124         DnaA -= DT*(0.025*DnaA+0.025*(1.1E-5*V/(1.1E-5*V+Y1))*DnaA); //DnaA degradation
00125 
00126 
00127 }
00128 
00129 void CDNA::ComputeBoxes()
00130 {
00131         for(int i=0; i<3; ++i)
00132                 DnaABoxes[i]=0;                         //Clear old data
00133         double marker=0;
00134 
00135         double BoxDist=29.596;                          //Box[1] distribution factor
00136         
00137         double FP[3];                           //FP is a modified fork position (distance from term)
00138         for(i=0; i<3; ++i)
00139                 FP[i]=1-ForkPos[i];
00140         
00141         //box distributions
00142         DnaABoxes[0]=4.368*(pow(FP[0],2)+2*(pow(FP[1],2)-pow(FP[0],2))
00143                         +4*(pow(FP[2],2)-pow(FP[1],2))+8*(1-pow(FP[2],2)))
00144                         +5.1201*(FP[0]+2*(FP[1]-FP[0])+4*(FP[2]-FP[1])+8*(1-FP[2]));
00145         DnaABoxes[1]=BoxDist*(pow(FP[0],2)+2*(pow(FP[1],2)-pow(FP[0],2))
00146                         +4*(pow(FP[2],2)-pow(FP[1],2))+8*(1-pow(FP[2],2)))
00147                         +63.694*(FP[0]+2*(FP[1]-FP[0])+4*(FP[2]-FP[1])+8*(1-FP[2]));
00148         DnaABoxes[2]=(4639221/9)*(FP[0]+2*(FP[1]-FP[0])+4*(FP[2]-FP[1])+8*(1-FP[2]));
00149 
00150         for(i=0;i<3;++i)
00151                 DnaABoxes[i]*=NumChrome;
00152 
00153         return;
00154 }
00155 
00156 
00157 void CDNA::Termination()
00158 {
00159         //TERMINATION--Only terminate when fork reaches end
00160         if(ForkPos[0]>=1.0)
00161         {
00162                 NFTERM=NumForks[0];     //Number terminated = forks on 1st strand
00163                 NTOT-=NFTERM;
00164                 NumChrome *=2;  //Twice the separate chromosomes
00165                 NPOS -= 1;              //Lose a fork
00166                 for(int i=0;i<(NPOS+1);++i)
00167                 {
00168                         NumForks[i]=NumForks[i+1];
00169                         ForkPos[i]=ForkPos[i+1];
00170                 }       
00171         }
00172 }
00173 
00174 void CDNA::OldInit()
00175 {
00176         double RP = 3.82E-17;   //Repressor protein
00177         double ARP = 3.82E-17;  //Anti-repressor protein
00178 }
00179 
00180 void CDNA::AdvanceFork()
00181 {
00182         if(NPOS>=1)             //If there's a fork...
00183         {
00184                 double DY8S=m_pCell->GetComponents()->GetDY8S();
00185                 double DT=m_pCell->GetComponents()->DT;
00186                 for(int i=0;i<NPOS;++i)                                 //...increase fork position
00187                         ForkPos[i] += DY8S*DT*2.49E14/NTOT;     //based on DNA growth
00188         }
00189 }
00190 
00191 double CDNA::Equilibrium()
00192 {       
00193         double V=m_pCell->GetGeometry()->GetVolume();
00194         //convert to nanomolar concentrations
00195         double DnaAConc = DnaA/52500/V*1E9*1000;
00196         double BoxConc[3];
00197         for(int i=0;i<3;++i)
00198                 BoxConc[i]=DnaABoxes[i]/V/6.023E23*1E9*1000;
00199 
00200         //FirstGuess for free DnaA concentration        
00201         double FreeDnaA=DnaAConc/(1+BoxConc[0]/Kd[0]+BoxConc[1]/Kd[1]+BoxConc[2]/Kd[2]);
00202 
00203         //Newton's Method for freeDnaA concentration    
00204         double NewFree=1;
00205         double error=0;
00206         for(i=0;error<(0.0001*NewFree);++i)
00207         {
00208                 NewFree = DnaAConc/(1+BoxConc[0]/(FreeDnaA+Kd[0])+BoxConc[1]
00209                         /(FreeDnaA+Kd[1])+BoxConc[2]/(FreeDnaA+Kd[2]));
00210                 error=fabs(NewFree-FreeDnaA);
00211                 FreeDnaA = NewFree;
00212                 if(i==500)
00213                         break;
00214         }
00215         
00216         return FreeDnaA;
00217 }
00218 
00219 double CDNA::gamma(double Z)
00220 {
00221         if(Z<0||Z>100)
00222                 m_pCell->GetDisplay()->ErrorMessage(2);
00223         const double G = 5;
00224         const double PI = 3.14159265358979324E0;
00225         if(Z<1){return gamma(Z+1)/Z;}
00226         return pow((2*PI),0.5) * pow((Z + G - 0.5),(Z - 0.5)) / exp(Z + G - 0.5) *
00227          (1.000000000178E0 + 76.180091729406E0 / Z - 86.505320327112E0/(Z+1) +
00228          24.014098222230E0 / (Z + 2) - 1.231739516140E0 / (Z + 3) +
00229           0.001208580030E0 / (Z + 4) - 0.000005363820E0 / (Z + 5) );
00230 }
00231 
00232 double CDNA::OddGD()    //this makes absolutely no sense but was in the original model
00233 {
00234         if(NTOT>=1)                     //only relevent if there are forks
00235         {
00236                 if(ForkPos[0]>=0.02)
00237                         GD = 8.0*NTOT;
00238                 if(ForkPos[0]>=0.06)
00239                         GD = 9.0*NTOT;
00240                 if(ForkPos[0]>=0.125)
00241                         GD = 11.0*NTOT;
00242                 if(ForkPos[0]>=0.23)
00243                         GD = 12.0*NTOT;
00244                 if(ForkPos[0]>=0.44)
00245                         GD = 13.0*NTOT;
00246                 if(ForkPos[0]>=0.52)
00247                         GD = 14.0*NTOT;
00248         }
00249         return GD;
00250 }
00251 
00252 void CDNA::Divide()
00253 {
00254         NumChrome /= 2;
00255         if(NPOS>=1)                                     //if there were forks in the DNA...
00256         {
00257                 for(int i=0;i<NPOS;++i)
00258                         NumForks[i]=NumForks[i]/2;      //...now there are 1/2 as many
00259         }
00260         NTOT=NTOT/2;            //1/2 the total forks
00261         GD=GD/2.0;                      //1/2 the gene dosage
00262         DnaA=0.5*DnaA;          //1/2 the DnaA
00263         NORG=NORG/2;            //half the origins
00264 }
00265 
00266 int CDNA::CheckTerm()
00267 {
00268         if(NFTERM>=1)
00269         {
00270                 NFTERM=0;
00271                 return 1;
00272         }
00273         else
00274                 return 0;
00275 }