Components.cpp Source File

00001 // Components.cpp: implementation of the CComponents class.
00002 //
00004 // #include <iostream.h>
00005 #include <math.h>
00006 #include "Cell.h"
00007 #include "Display.h"
00008 #include "DNA.h"
00009 #include "Geometry.h"
00010 #include "Components.h"
00011 
00013 // Construction/Destruction
00015 
00016 
00017 //this is the constructor we make to call the parent cell
00018 CComponents::CComponents(CCell* pCell)
00019 {       
00020         m_pCell=pCell;                                                  //Set pointer to parent
00021         Initialize();
00022 }
00023 
00024 //this is the default constructor  
00025 CComponents::CComponents()
00026 {
00027                                 
00028 }
00029 
00030 //this is the default destructor 
00031 CComponents::~CComponents()
00032 {
00033 
00034 }
00035 
00036 CComponents& CComponents::operator =(const CComponents & rhs)   //assignment operator
00037 {
00038         return *this;
00039 }
00040 
00041 CComponents::CComponents(const CComponents& rComp)                      //copy constructor
00042 {
00043         m_pCell=rComp.m_pCell;                                                                  //In this case, same parent                                                                                                             //Same dimensions
00044         for(int i=0;i<NBIOVARS;++i)
00045                 InitialValues[i] = rComp.InitialValues[i];      
00046 }
00047 
00048 //Overloaded copy constructor
00049 CComponents::CComponents(const CComponents &rComp, CCell *pCell)
00050 {
00051         m_pCell=pCell;                                                                                  //New parent cell
00052         test=rComp.test;                                                                                //just for testing out code
00053 
00054         
00055 }
00056 
00057 void CComponents::Integrate()
00058 {
00059         double DY[NBIOVARS]={0};        
00060         double DYB[NBIOVARS]={0.0};             //For computing averages
00061 
00062         if(NT2)
00063         {
00064                 FirstInt();
00065                 NT2=0;
00066         }
00067         PredictorCorrector();
00068         m_pCell->AdvanceTime(DT);
00069         Flux(Y2,DYB);                           //New derivatives
00070         for(int i=0;i<NBIOVARS;++i)
00071                 Averages[i]=Averages[i]+0.5*(2*Y2[i]+DYB[i]*DT)*DT;
00072         DecideStepSize();       
00073         TruncatedError();
00074         m_pCell->GetGeometry()->ComputeVolume();
00075         Flux(Y,DY);                             //we're taking derivatives again
00076         CheckNegConc();
00077 }
00078 
00079 void CComponents::Initialize()
00080 {
00081         InitialValues[0] = 1.00E-4;
00082         InitialValues[1] = 1.54E-3;
00083         InitialValues[2] = 7.00E-4;
00084         InitialValues[3] = 6.60E-4/3;
00085         InitialValues[4] = 2.97E-4*2;
00086         InitialValues[5] = 1.08E-4*2;
00087         InitialValues[6] = 1.76E-1; 
00088         InitialValues[7] = 6.36E-2;
00089         InitialValues[8] = 1.10E-2;
00090         InitialValues[9] = 2.2E-2;
00091         InitialValues[10] = 10.7E-3;
00092         InitialValues[11] = 7.34E-4;
00093         InitialValues[12] = 5.62E-4;
00094         InitialValues[13] = 10.35E-4;
00095         InitialValues[14] = 3.632E-3;
00096         InitialValues[15] = 5.632E-2;
00097         InitialValues[16] = 0.035/10;
00098         InitialValues[17] = 6.67E-8*2;
00099         InitialValues[18] = 0.22*0.035/10; 
00100         InitialValues[19] = 2.0E-3*10;
00101 
00102         double V = 5.65E-13;                    //Volume;
00103         for(int i=0; i<NBIOVARS; i++)
00104                 Y0[i]=InitialValues[i]*V;               //Calculates grams of components
00105         for(i=0;i<NBIOVARS;++i)
00106                 Y1[i]=0;
00107         for(i=0;i<NBIOVARS;++i)
00108                 Y2[i]=0;
00109         for(i=0;i<NBIOVARS;++i)
00110                 Y2EST[i]=0;
00111         for(i=0;i<27;++i)
00112                 KR[i]=1;
00113 
00114 
00115         double P[NBIOVARS]={0.0};               //concentrations
00116         for(i=0; i<NBIOVARS; i++)
00117                 Y[i]=Y0[i];                     //Set components
00118         NT2 = 0;                //Number of times step size decreased
00119         ITN=0;
00120         DT = 0.0013;            //Time step
00121         DEV = 0.001;            //Allowable deviation in corrector
00122         E3Rate=1.6E-2;                          //Rate constant for E3
00123         CA1=1E-3;
00124         CA2=1E-3;
00125         CA3=1E-3;
00126         DY6S=0;
00127         DY8S=0;
00128         DY9=0;
00129 
00130 }
00131 
00132 void CComponents::StartAverages()
00133 {
00134         for(int i=0;i<NBIOVARS;i++)
00135                 Averages[i]=0.0;                //Remove any previous values for BY[i]
00136         TotalTime=0.0;  //Keeps track of total time of cell life
00137 }
00138 
00139 double CComponents::GetY(int i)
00140 {
00141         return Y[i];
00142 }
00143 
00144 void CComponents::FirstInt()
00145 {       
00146         double Y1EST[NBIOVARS]={0};
00147         double DY0[NBIOVARS]={0};
00148         double DY1EST[NBIOVARS]={0};
00149         Flux(Y0,DY0);   //First Euler derivative
00150         for(int i=0;i<NBIOVARS;++i)
00151                 Y1EST[i]=Y0[i]+DY0[i]*DT;
00152         Flux(Y1EST,DY1EST);
00153         for(i=0;i<NBIOVARS;++i)
00154                 Y1[i]=Y0[i]+0.5*DT*(DY0[i]+DY1EST[i]);  //Euler result
00155         NT2=0;
00156 }
00157 
00158 void CComponents::PredictorCorrector()
00159 {
00160         //Predictor--Guess Y2 based on Y0 and Y1 slope
00161         double DY1[NBIOVARS]={0};       
00162         Flux(Y1,DY1);           //Predictor derivative
00163         for(int i=0;i<NBIOVARS;++i)
00164                 Y2EST[i]=Y0[i]+2.0*DT*DY1[i];
00165         ITN=0;
00166         double DARP=0;  //cell envelope slope
00167         double DY2EST[NBIOVARS]={0};
00168         //Correcting                            
00169         for(ITN=0;ITN<10;++ITN)
00170         {
00171                 Flux(Y2EST,DY2EST);
00172                 //Corrector
00173                 int ITN1=0;                             //Counts bad estimates in corrector                             
00174                 for(int i=0;i<NBIOVARS;++i)
00175                         Y2[i]=Y1[i]+0.5*DT*(DY1[i]+DY2EST[i]);  //Corrector result
00176                 for(i=0;i<NBIOVARS;++i)
00177                 {
00178                         if(fabs((Y2[i]-Y2EST[i])/Y2[i])<DEV)    //Acceptable deviation
00179                                 continue;
00180                         Y2EST[i]=Y2[i];                 //Set Y2 estimate to corrector guess
00181                         ++ITN1;                                 //Count unacceptable guesses
00182                 }
00183                 DARP = 0.5*(DY2EST[9]+DY1[9])*DT;       //ARP change=slope of cell envelope *DT
00184                 if(ITN1==0)
00185                         break;                          //Skip ahead if all guesses are OK
00186         }
00187         if(ITN >=10)                    //Loop broke after 10 iterations
00188         {
00189                 NT2++;
00190                 DT*=0.1;
00191                 Integrate();                            //Repeat Euler with new smaller step size
00192         }
00193 }
00194 
00195 void CComponents::DecideStepSize()
00196 {
00197         //Decide step size
00198         if(ITN==0)
00199         {
00200                 DT*=1.01;               //Increase step size by 1%
00201                 NT2=1;                  //signal DT has been modified
00202         }
00203         else
00204         {
00205                 if(ITN>=4)
00206                 {
00207                         DT*=0.5;        //If more than 3 tries, use 1/2 time step size
00208                         NT2=1;          //signal DT has been modified
00209                 }
00210         }
00211 }
00212 
00213 void CComponents::TruncatedError()
00214 {
00215         double DIF[NBIOVARS];
00216         for(int i=0;i<NBIOVARS;++i)
00217         {
00218                 DIF[i]=0.2*(Y2EST[i]-Y2[i]);    //truncated error
00219                 Y[i]=Y2[i]+DIF[i];      //adds truncated error
00220         }
00221 }
00222 
00223 void CComponents::CheckNegConc()
00224 {
00225         //Assuring no negative concentration
00226         for(int i=0;i<NBIOVARS;++i)
00227         {
00228                 if(Y[i]<0.0)
00229                         m_pCell->GetDisplay()->ErrorMessage(1);
00230         }
00231 }
00232 
00233 
00234 
00235 void CComponents::Flux(double* Y, double* DY)
00236 {
00237         const double R=1.5;                             //P/O ratio
00238         double V=m_pCell->GetGeometry()->GetVolume();
00239         double S=m_pCell->GetGeometry()->GetSurface();
00240         //Differential equations
00241         //idling ribosomes for ppGpp balance
00242         //(those not involved in protein synthesis
00243         double RI=(1.0-(Y[1]/(Y[1]+V*7.2E-5))*(Y[2]/(Y[2]+V*22.0E-5)))*Y[15];
00244         
00245         //ppGpp
00246         double DY17S=KR[1]*6.4E-3*RI*(12.65E-5/(12.65E-5+Y[2]/V));      //ppGpp synthesis
00247         double DY17D=KR[2]*125.0*Y[17]*(Y[1]/(Y[1]+V*7.2E-5));          //ppGpp degradation
00248         DY[17]=DY17S-DY17D;                                                                     //ppGpp balance in moles
00249 
00250         //R-protein mRNA synthesis.  It is believed this is growth modulated
00251         //thru ppGpp and varies according to gene dosage based on how Y[16]
00252         //varies.  This is a fraction of total RNA but is already included in
00253         //the total mass thru mRNA(Y[16])
00254 
00255         //total mRNA synthesis
00256         double DY16S=KR[3]*2.0E-14*(Y[3]/(Y[3]+V*7.2E-5))*(Y[8]/4.0E-15);
00257 
00258         //RIBOSOMAL CONTENT
00259         double RIBOSOME=Y[15]/(4566.0*5.4E-22); //ribosomal RNA in molecules
00260         double RIBOGRAM=Y[15];                                  //..and in grams
00261 
00262         //Calculating FRACRT, the fraction of ribosomes translating mRNA
00263         double AFFRIBO=8.75E5;          //affinity of ribosomes for mRNA (Perretti 1987)
00264         double RNAMW=3.9E5;             //molecular weight of average mRNA
00265         double RNAMOL=Y[16]/(RNAMW*V)*1000.0;   //molar concentration of ribosomal
00266                                                                         //binding sites
00267         double FRACRT=(AFFRIBO*RNAMOL)/(AFFRIBO*RNAMOL+1.0);    //fraction of translating
00268                                                                                                         //ribosomes
00269         //Protein synthesis based on mRNAL content
00270         //CP is the rate of protein elongation, 23aa/sec, RIBOSOME
00271         double CP=KR[4]*23.0*(1.8E-22*3600.);
00272 
00273         //Protein Balance (M1)
00274         DY6S=RIBOSOME*FRACRT*CP*(Y[1]/(7.2E-5*V+Y[1]))
00275                 *(Y[2]/(2.2E-4*V+Y[2]));                //Protein synthesis
00276         double DY6D=KR[5]*(0.025*Y[6]+0.025*(1.1E-5*V/(1.1E-5*V+Y[1]))*Y[6]);
00277                                                                                 //Protein degradation
00278         DY[6]=DY6S-DY6D;                                        //Protein balance
00279 
00280         //Degradation of mRNA depends on translational efficiency
00281         double DY16D=KR[6]*21.0*Y[16];                  //rate of mRNA degradation
00282         DY[16]=DY16S-DY16D;                             //overall mRNA balance
00283 
00284         //r-protein translation Y[19]
00285         //Inhibition of translation by free R-proteins
00286         //Guess free R-proteins (M)
00287         double RRNAT=(Y[14]/1.9E6)/V*1000.0;    //molar quantity of precursor rRNA
00288         double RNAMT=(Y[18]/8.4E6)/V*1000.0;    //molar quantity of r-protein mRNA
00289         double RIBOMOL=RIBOSOME/(6.02E23);              //moles of ribosomes
00290         double RPROTEINT = Y[19]/(7336.0*1.8E-22*6.02E23);      //total mol ribosomal protein
00291         double RPROTEIN=RPROTEINT-RIBOMOL;      //ribosomal protein not in ribosomes
00292         if(RPROTEIN<0.0)
00293                 RPROTEIN=0.0;           //don't let free r-protein be negative
00294         double FREEPRO=RPROTEIN*0.2/V*1000.0;   //first guess at amount of free r-proteins
00295         double CRRNA=1.0E6;                             //binding constant to rRNA
00296         double CMRNA=8.0E4;                             //binding constant to r-protein mRNA
00297         int KCOUNTP=0;                                  //counter to zero for each FCT
00298         double RRNAFREE=0;
00299         double RNAMFREE=0;
00300         for(KCOUNTP=0;KCOUNTP<500;KCOUNTP++)    //for 500 iterations
00301         {
00302                 RRNAFREE=(RRNAT/(1.0+FREEPRO*CRRNA));   //equilibrium
00303                 RNAMFREE=(RNAMT/(1.0+FREEPRO*CMRNA));   //equilibrium
00304                 double FREEPROT=(RPROTEIN/V)*1000.0;            
00305                 double FREEPRO1=(FREEPROT/(1.0+(RRNAFREE*CRRNA)+(RNAMFREE*CMRNA)));
00306                                                                                 //new guess via mass balance
00307                 double  ERRPRO=fabs(FREEPRO-FREEPRO1);          //error of guess
00308                 FREEPRO=FREEPRO1;                                       //new guess
00309                 if(ERRPRO<=(0.0001*FREEPRO1))   
00310                         break;                                  //jump the loop if error is small
00311         }
00312         
00313         if(KCOUNTP>499)
00314                 m_pCell->GetDisplay()->WarningMessage(3);
00315 
00316         double RRNA=(1.0-RRNAFREE/RRNAT)*RRNAT*(V/1000.0);
00317                                                                         //calculate bound rRNA precursors
00318         //R-protein balance
00319         double DY19S=DY6S*(RNAMFREE*(V/1000.0)*8.4E6)/Y[16];    //synthesis proportional
00320                                                                                                         //to mRNA fraction
00321         double DY19D=(FREEPRO*(V/1000.0)*(7.9E5))/Y[6]*DY6D;    //only free r-proteins
00322                                                                                                         //degrade
00323         DY[19]=DY19S-DY19D;                                             //balance
00324         
00325         //R-PROTEIN mRNA balance
00326         double DY18S=0.22*(2.0E-7/(2.0E-7+Y[17]/V))*DY16S;      //synthesis
00327         double DY18D=KR[7]*21.0*Y[18];                                          //degradation
00328         DY[18]=DY18S-DY18D;                                                     //balance
00329 
00330         //mature rRNA balance
00331         double DY15S=KR[8]*85.0*RRNA;           //synthesis
00332         double DY15D=KR[9]*0.0;                 //ignore degradation
00333         DY[15]=DY15S*1.5E6-DY15D;       //balance
00334 
00335         //immature sRNA balance
00336         //max rate of sRNA synthesis from Bremer and Dennis
00337         double DY14S=KR[10]*(1.15E-14)*(Y[3]/(Y[3]+V*7.2E-5))*
00338                 (2.0E-7/(2.0E-7+Y[17]/V))*
00339                 (0.009/(0.009+pow((2.1*Y[6]/V-DY6S/V),3)))*m_pCell->GetDNA()->OddGD();  //synthesis
00340         double DY14D=KR[11]*21.0*Y[14];                                                 //degradation
00341         DY[14]=DY14S-DY14D-DY15S*1.9E6; //last term in balance for maturation
00342         
00343         //glycogen balance
00344         double DY10S=KR[12]*2.0E-2*(Y[1]/(1.4E-3*V+Y[1]))*V;            //synthesis
00345         double DY10D=KR[13]*0.14*(1.0E-3*V/(1.0E-3*V+Y[1]))*Y[10];      //degradation
00346         DY[10]=DY10S-DY10D;                                                             //balance
00347 
00348         //cell envelope balance
00349         double DY9S=KR[14]*100.60*(Y[5]/(5.0E-4*V+Y[5]))*
00350                 (Y[1]/(1.80E-4*V+Y[1]))*m_pCell->GetGeometry()->E23;    //synthesis
00351         double DY9D=KR[15]*0.23*Y[9]*(Y[1]/(Y[1]+V*1.8E-4));            //degradation
00352         DY[9]=DY9=DY9S-DY9D;                                                            //balance
00353         
00354         //DNA balance, M3
00355         DY[8]=DY8S=KR[16]*6.0E-15*(Y[4]/(4.0E-6*V+Y[4]))*
00356                 (Y[1]/(2.0E-5*V+Y[1]))*m_pCell->GetDNA()->NTOT;                         
00357 
00358         //total RNA balance, M2
00359         //total RNA synthesized:  26% in rrn operons, 15% as total sRNA(mature)
00360         //      let rate of synthesis be DY14S*0.04/0.26 since about 4% of rrn is
00361         //      tRNA, and degradation is negligible.
00362         DY[7]=DY[14]+DY[15]+DY[16]+DY14S*0.04/0.26;
00363 
00364         //cell envelope precursors balance, P4
00365         double DY5S=KR[17]*0.06*(Y[2]/(9.9E-4*V+Y[2]))
00366                 *(Y[1]/(1.2E-4*V+Y[1]))                         //synthesis
00367                 *(5.0E-3*V/(5.0E-3*V+Y[5]))*V;          
00368         DY[5]=DY5S-1.1*DY[9];                                   //balance
00369 
00370         //deoxyribonucleotides balance, P3
00371         double DY4S=KR[18]*40.0*(Y[3]/(1.44E-5*V+Y[3]))
00372                 *(Y[1]/(7.2E-5*V+Y[1]))*                        //synthesis
00373                 (2.2E-4*V/(2.2E-4*V+Y[4]))*Y[11];       
00374         DY[4]=DY4S-1.053*DY[8];                                 //balance
00375 
00376 
00377         //ribonucleotides balance, P2
00378         double DY3S=KR[19]*1.9E-1*(Y[2]/(9.9E-4*V+Y[2]))*(Y[1]/(2.5E-4*V+Y[1]))
00379                 *(9.0E-3*V/(9.0E-3*V+Y[3]))*V;          //synthesis
00380         double DY3D=KR[20]*0.03*(1.1E-5*V/(1.1E-5*V+Y[1]))*Y[3];        //degradation
00381         DY[3]=DY3S-DY3D-1.057*(DY14S-DY14D-DY15D+DY[16])
00382                 -1.049*DY4S+1.057*(0.2)*DY15S*1.9E6;    //balance, last term for
00383 
00384         //rRNA processing
00385 
00386         //amino acids balance, P1
00387         double DY2S=KR[21]*0.39*(5.0E-2*V/(5.0E-2*V+Y[2]))*(Y[0]/(2.2E-5*V+Y[0]))
00388                 *(Y[1]/(2.5E-4*V+Y[1]))*V;                      //synthesis
00389         double DY2D=KR[22]*0.025*(1.1E-5*V/(1.1E-5*V+Y[1]))*Y[2];       //degradation
00390         DY[2]=DY2S-DY2D-1.167*DY[6]-1.149*(DY3S-DY3D)-0.128*DY5S;       //balance
00391                 //(amino acids become proteins, nulcleotides, and envelope precursors)
00392 
00393         //ammonium ion uptake
00394         double RA1=KR[23]*(2.0E-7*(CA1/(1.08E-7+CA1)-Y[0]/(1.0E-5*V+Y[0]))
00395                 +8.2E-7*(CA1/(3.6E-6+CA1))*(2.3E-3*V/(2.3E-3*V+Y[0]))
00396                 *(Y[1]/(7.2E-5*V+Y[1])));
00397         //internal ammonium ion balance, A1
00398         DY[0]=RA1*S-0.179*(DY2S-DY2D);
00399 
00400         //parameters for energy balance
00401         double Z=0.9*(DY9S/(DY9S+1.0E-14));
00402         double DMT=0;
00403         for(int i=0;i<11;++i)
00404                 DMT += DY[i];
00405         double DVT=DMT/0.29;                    //estimated change in volume
00406         
00407         //glucose anabolism
00408         double DY1A=(DY2S-DY2D)*1.128-(DY3S-DY3D)*0.456+DY5S*1.26+DY[10]*1.11;
00409 
00410         //energy balance
00411         double DATP=DY[1]*0.0056+DY2S*0.0025+DY3S*0.022+DY5S*0.0016+DY10S*0.0124
00412                 +DY6S*0.039+DY[8]*0.0071+DY9S*0.0081+DY17S*3.0+DVT*0.007
00413                 +(DY14S+DY16S)*0.0067;
00414         double DE=RA1*S*0.028+DMT*2.0E-3+S*5.8E-8
00415                 +DMT*1.39E-2+DMT*1.21E-3+DY4S*0.0031;   //dX/dt
00416 
00417         //glucose catabolism
00418         double DY1C=((DATP+DE*R)/(4.0+(12.0-4.0*Z)*R))*180.0;
00419         //glucose uptake
00420         double RA2=KR[24]*1.7E-5*(CA2/(1.2E-5+CA2))
00421                 *(7.78/(7.78+pow(((Y[1]/V)*1.0E3),4)));
00422         //internal glucose balance (A2)
00423         DY[1]=RA2*S-DY1A-DY1C;
00424 
00425         //enzyme balances, E1, E2, and E3
00426         DY[11]=1.3E-3*m_pCell->GetDNA()->NORG*DY6S*m_pCell->GetDNA()->MidBound/0.88-
00427                 KR[25]*(0.025*Y[11]+0.025*(1.1E-5*V/(1.1E-5*V+Y[1]))*Y[11]);
00428         DY[12]=3.2E-3*DY[6];
00429         DY[13]=KR[26]*E3Rate*Y[6];
00430 }
00431 
00432 
00433 double CComponents::GetDY6S()
00434 {
00435         return DY6S;
00436 }
00437 
00438 double CComponents::GetDY8S()
00439 {
00440         return DY8S;
00441 }
00442 
00443 double CComponents::GetDY9()
00444 {
00445         return DY9;
00446 }
00447 
00448 void CComponents::ResetY(int i)
00449 {
00450         Y[i]=0.0;
00451         Y1[i]=0.0;
00452         Y2[i]=0.0;
00453 }
00454 
00455 void CComponents::Divide()
00456 {
00457         for(int i=0;i<NBIOVARS;++i)
00458                 Y0[i]=0.5*Y[i];
00459 }
00460 
00461 void CComponents::NextY()
00462 {
00463         for(int i=0;i<NBIOVARS;++i)
00464         {                                       //set values for the next iteration
00465                 Y0[i]=Y1[i];
00466                 Y1[i]=Y2[i];
00467         }
00468 }