EGFRTraffickingCHOMinimizableDirector.cpp Source File

00001 // EGFRTraffickingCHOMinimizableDirector.cpp: implementation of the CEGFRTraffickingCHOMinimizableDirector class.
00002 //
00004 
00005 #include "EGFRTraffickingCHOMinimizableDirector.h"
00006 
00008 // Construction/Destruction
00010 
00011 CEGFRTraffickingCHOMinimizableDirector::CEGFRTraffickingCHOMinimizableDirector()
00012 {
00013         // create networks,experiments,observers,movers,minimizables
00014         double maxSimTime = 241.0;
00015         double stepSize = 0.0001;
00016         double frequency = 1.0;
00017         int nCells = 1;
00018         double moveTol = 1.0e-06;
00019 
00020         m_dGammaSquared = 1.0e-07;
00021         
00022         // flags for all minimizables
00023         bool rateFlag = true;
00024         bool chemFlag = false;
00025         bool logsFlag = false;
00026         double timeWeight = 1.0;
00027         double rateWeight = 0.0;
00028         double chemWeight = 0.0;
00029         
00030 
00031         // CHO cells
00032         // expt. 0 - CHO cells transfected with ? receptors, treated with EGF at 5 ng/ml
00033         //              concentration
00034         _networkList.push_back(new CEGFRTraffickingNetwork());
00035         _experimentList.push_back(new Experiment(_networkList[0],"EGF at 5 ng/ml, Bao 2000"));
00036         // put data files here
00037         _experimentList[0]->AddDataFileToList("../DataFiles/TimeSeries/InternalizedEGFREGFYarden.dat");
00038         _experimentList[0]->AddForcingFileToList("../DataFiles/TimeSeries/TwoHourEGF.dat");
00039         // mover and plotter
00040         //_moverList.push_back(new CQualityControlCashKarpMover(frequency,stepSize,moveTol));
00041         _moverList.push_back(new CStiffBulirschStoerMover(frequency,stepSize,moveTol));
00042         // minimizable
00043         _minimizableList.push_back(new SingleNetworkMinimizable());
00044         // parameters
00045         _minimizableList[0]->Initialize(_experimentList[0],_moverList[0],nCells,rateFlag,chemFlag,logsFlag,timeWeight,rateWeight,chemWeight);
00046         _minimizableList[0]->SetRateConstantsWeight(0.0);
00047         // attach observer
00048         
00049         _plotterList.push_back(new GnuPlotterTimeSeriesObserver());
00050         _plotterList[0]->SetTitle("EGFR Trafficking Network:EGF at 5 ng/ml (Bao 2000):ODE");
00051         _minimizableList[0]->Attach(_plotterList[0]);
00052         
00053 
00054         // expt. 1 - CHO cells transfected with 8.2(7) X 10^4 receptors, treated 
00055         //              with EGF at 100 nM concentration
00056         _networkList.push_back(new CEGFRTraffickingNetwork());
00057         _experimentList.push_back(new Experiment(_networkList[1],"EGF at 100 nM, Davis 1992"));
00058         // put data files here
00059         _experimentList[1]->AddDataFileToList("../DataFiles/TimeSeries/SurfaceEGFREGFDavis100nM.dat");
00060         _experimentList[1]->AddForcingFileToList("../DataFiles/TimeSeries/TwoHourEGF.dat");
00061         // mover and plotter
00062         //_moverList.push_back(new CQualityControlCashKarpMover(frequency,stepSize,moveTol));
00063         _moverList.push_back(new CStiffBulirschStoerMover(frequency,stepSize,moveTol));
00064         // minimizable
00065         _minimizableList.push_back(new SingleNetworkMinimizable());
00066         // parameters
00067         _minimizableList[1]->Initialize(_experimentList[1],_moverList[1],nCells,rateFlag,chemFlag,logsFlag,timeWeight,rateWeight,chemWeight);
00068         _minimizableList[1]->GetConversionFactor(59)->SetFactorFixed(true);
00069         _minimizableList[1]->GetConversionFactor(59)->SetFactorValue(100.0/82000.0);
00070         _minimizableList[1]->SetRateConstantsWeight(0.0);
00071         // attach observer
00072         
00073         _plotterList.push_back(new GnuPlotterTimeSeriesObserver());
00074         _plotterList[1]->SetTitle("EGFR Trafficking Network:EGF at 100 nM (Davis 1992):ODE");
00075         _minimizableList[1]->Attach(_plotterList[1]);
00076         
00077 
00078         // expt. 2 - CHO cells transfected with 8.2(7) X 10^4 receptors, treated 
00079         //              with EGF at 10 nM concentration
00080         _networkList.push_back(new CEGFRTraffickingNetwork());
00081         _experimentList.push_back(new Experiment(_networkList[2],"EGF at 10 nM, Davis 1992"));
00082         // put data files here
00083         _experimentList[2]->AddDataFileToList("../DataFiles/TimeSeries/SurfaceEGFREGFDavis10nM.dat");
00084         _experimentList[2]->AddForcingFileToList("../DataFiles/TimeSeries/TwoHourEGF.dat");
00085         // mover and plotter
00086         //_moverList.push_back(new CQualityControlCashKarpMover(frequency,stepSize,moveTol));
00087         _moverList.push_back(new CStiffBulirschStoerMover(frequency,stepSize,moveTol));
00088         // minimizable
00089         _minimizableList.push_back(new SingleNetworkMinimizable());
00090         // parameters
00091         _minimizableList[2]->Initialize(_experimentList[2],_moverList[2],nCells,rateFlag,chemFlag,logsFlag,timeWeight,rateWeight,chemWeight);
00092         _minimizableList[2]->GetConversionFactor(59)->SetFactorFixed(true);
00093         _minimizableList[2]->GetConversionFactor(59)->SetFactorValue(100.0/82000.0);
00094         _minimizableList[2]->SetRateConstantsWeight(0.0);
00095         // attach observer
00096         
00097         _plotterList.push_back(new GnuPlotterTimeSeriesObserver());
00098         _plotterList[2]->SetTitle("EGFR Trafficking Network:EGF at 10 nM (Davis 1992):ODE");
00099         _minimizableList[2]->Attach(_plotterList[2]);
00100         
00101 
00102         // expt. 3 - CHO cells transfected with 8.2(7) X 10^4 receptors, treated 
00103         //              with EGF at 1 nM concentration
00104         _networkList.push_back(new CEGFRTraffickingNetwork());
00105         _experimentList.push_back(new Experiment(_networkList[3],"EGF at 1 nM, Davis 1992"));
00106         // put data files here
00107         _experimentList[3]->AddDataFileToList("../DataFiles/TimeSeries/SurfaceEGFREGFDavis1nM.dat");
00108         _experimentList[3]->AddForcingFileToList("../DataFiles/TimeSeries/TwoHourEGF.dat");
00109         // mover and plotter
00110         //_moverList.push_back(new CQualityControlCashKarpMover(frequency,stepSize,moveTol));
00111         _moverList.push_back(new CStiffBulirschStoerMover(frequency,stepSize,moveTol));
00112         // minimizable
00113         _minimizableList.push_back(new SingleNetworkMinimizable());
00114         // parameters
00115         _minimizableList[3]->Initialize(_experimentList[3],_moverList[3],nCells,rateFlag,chemFlag,logsFlag,timeWeight,rateWeight,chemWeight);
00116         _minimizableList[3]->GetConversionFactor(59)->SetFactorFixed(true);
00117         _minimizableList[3]->GetConversionFactor(59)->SetFactorValue(100.0/82000.0);
00118         _minimizableList[3]->SetRateConstantsWeight(0.0);
00119         // attach observer
00120         
00121         _plotterList.push_back(new GnuPlotterTimeSeriesObserver());
00122         _plotterList[3]->SetTitle("EGFR Trafficking Network:EGF at 1 nM (Davis 1992):ODE");
00123         _minimizableList[3]->Attach(_plotterList[3]);
00124         
00125         /*
00126         // expt. 4 - CHO cells transfected with ? EGFR and treated with EGF at 
00127         //              100 ng/ml concentration
00128         _networkList.push_back(new CEGFRTraffickingNetwork());
00129         _experimentList.push_back(new Experiment(_networkList[4],"EGF at 100 ng/ml, Levkowitz 1999"));
00130         // put data files here
00131         _experimentList[4]->AddDataFileToList("../DataFiles/TimeSeries/SurfaceEGFREGFLevkowitz.dat");
00132         _experimentList[4]->AddForcingFileToList("../DataFiles/TimeSeries/TwoHourEGF.dat");
00133         // mover and plotter
00134         //_moverList.push_back(new CQualityControlCashKarpMover(frequency,stepSize,moveTol));
00135         _moverList.push_back(new CStiffBulirschStoerMover(frequency,stepSize,moveTol));
00136         // minimizable
00137         _minimizableList.push_back(new SingleNetworkMinimizable());
00138         // parameters
00139         _minimizableList[4]->Initialize(_experimentList[4],_moverList[4],nCells,rateFlag,chemFlag,logsFlag,timeWeight,rateWeight,chemWeight);
00140         _minimizableList[4]->GetConversionFactor(59)->SetFactorFixed(true);
00141         _minimizableList[4]->GetConversionFactor(59)->SetFactorValue(100.0/80000.0);
00142         _minimizableList[4]->SetRateConstantsWeight(0.0);
00143         // attach observer
00144         
00145         _plotterList.push_back(new GnuPlotterTimeSeriesObserver());
00146         _plotterList[4]->SetTitle("EGFR Trafficking Network:EGF at 100 ng/ml (Levkowitz 1999):ODE");
00147         _minimizableList[4]->Attach(_plotterList[4]);
00148         */
00149 
00150         /*
00151         // expt. 5 - CHO cells transfected with ? EGFR Y1045F and treated with EGF at
00152         //              100 ng/ml concentration
00153         _networkList.push_back(new CEGFRTraffickingNetwork());
00154         _experimentList.push_back(new Experiment(_networkList[5],"EGF at 100 ng/ml, Y1045F EGFR, Levkowitz 1999"));
00155         // put data files here
00156         _experimentList[5]->AddDataFileToList("../DataFiles/TimeSeries/SurfaceEGFREGFY1045FLevkowitz.dat");
00157         // mover and plotter
00158         _moverList.push_back(new CQualityControlCashKarpMover(frequency,stepSize,moveTol));
00159         // minimizable
00160         _minimizableList.push_back(new SingleNetworkMinimizable());
00161         // parameters
00162         _minimizableList[5]->Initialize(_experimentList[5],_moverList[5],nCells,rateFlag,chemFlag,logsFlag,timeWeight,rateWeight,chemWeight);
00163         _minimizableList[5]->GetConversionFactor(47)->SetFactorFixed(true);
00164         _minimizableList[5]->GetConversionFactor(47)->SetFactorValue(100.0/80000.0);
00165         _minimizableList[5]->SetRateConstantsWeight(0.0);
00166         // attach observer
00167         _plotterList.push_back(new GnuPlotterTimeSeriesObserver());
00168         _plotterList[5]->SetTitle("EGFR Trafficking Network:EGF at 100 ng/ml, Y1045F EGFR (Levkowitz 1999):ODE");
00169         _minimizableList[5]->Attach(_plotterList[5]);
00170         
00171         
00172         // expt. 6 - CHO cells transfected with ? EGFR K721A and treated with EGF at
00173         //              100 ng/ml concentration
00174         _networkList.push_back(new CEGFRTraffickingNetwork());
00175         _experimentList.push_back(new Experiment(_networkList[6],"EGF at 100 ng/ml, K721A, Levkowitz 1999"));
00176         // put data files here
00177         _experimentList[6]->AddDataFileToList("../DataFiles/TimeSeries/SurfaceEGFREGFK721ALevkowitz.dat");
00178         // mover and plotter
00179         _moverList.push_back(new CQualityControlCashKarpMover(frequency,stepSize,moveTol));
00180         _plotterList.push_back(new GnuPlotterTimeSeriesObserver());
00181         // title
00182         _plotterList[6]->SetTitle("EGFR Trafficking Network:EGF at 100 ng/ml, K721A EGFR (Levkowitz 1999):ODE");
00183         // minimizable
00184         _minimizableList.push_back(new SingleNetworkMinimizable());
00185         // parameters
00186         _minimizableList[6]->Initialize(_experimentList[6],_moverList[6],nCells,rateFlag,chemFlag,logsFlag,timeWeight,rateWeight,chemWeight);
00187         _minimizableList[6]->GetConversionFactor(7)->SetFactorFixed(true);
00188         _minimizableList[6]->GetConversionFactor(7)->SetFactorValue(100.0/80000.0);
00189         _minimizableList[6]->SetRateConstantsWeight(0.0);
00190         // attach observer
00191         _minimizableList[6]->Attach(_plotterList[6]);
00192         */
00193         
00194         // parameters which need a prior
00195 
00196         // allocate the master residuals list
00197         int masterResidualsSize = 0;
00198         for(int i = 0; i < _minimizableList.size(); i++)
00199         {
00200                 int subRes = _minimizableList[i]->GetNResiduals();
00201                 masterResidualsSize += subRes;
00202                 cout << subRes << " residuals in minimizable " << i << endl;
00203         }
00204         // add residuals for prior parameter costs, if applicable
00205         masterResidualsSize += m_iPriorList.size();
00206         cout << m_iPriorList.size() << " parameters are constrained by a prior " << endl;
00207         if(m_dGammaSquared > 0.0)
00208         {
00209                 masterResidualsSize += 1;
00210         }
00211         Allocate(masterResidualsSize);
00212 
00213         cout << "Master Residuals size : " << this->GetNResiduals() << endl;
00214 
00215         // define the experiments
00216         DefineExperiments();
00217 
00218 }
00219 
00220 CEGFRTraffickingCHOMinimizableDirector::~CEGFRTraffickingCHOMinimizableDirector()
00221 {
00222         for(int i = 0; i < _plotterList.size(); i++)
00223         {
00224                 _minimizableList[i]->Detach(_plotterList[i]);
00225                 delete _plotterList[i];
00226         }
00227 }
00228 
00229 void CEGFRTraffickingCHOMinimizableDirector::DefineExperiments()
00230 {       
00231         
00232         // CHO cells
00233         // 80,000 receptors assumed, EGF at 5 ng/ml
00234         _minimizableList[0]->GetExperiment()->GetReactionNetwork()->RemoveReaction(98);
00235         _minimizableList[0]->GetExperiment()->GetReactionNetwork()->GetChemical(0)->SetInitialAmount(600000.0);
00236         _minimizableList[0]->GetExperiment()->GetReactionNetwork()->GetChemical(3)->SetInitialAmount(10000.0);
00237         _minimizableList[0]->GetExperiment()->GetReactionNetwork()->GetChemical(4)->SetInitialAmount(70000.0);
00238         _minimizableList[0]->GetExperiment()->GetReactionNetwork()->ChemicalReset();
00239 
00240         // 82,000 receptors, EGF at 100 nM
00241         _minimizableList[1]->GetExperiment()->GetReactionNetwork()->RemoveReaction(98);
00242         _minimizableList[1]->GetExperiment()->GetReactionNetwork()->GetChemical(0)->SetInitialAmount(7.2e+007);
00243         _minimizableList[1]->GetExperiment()->GetReactionNetwork()->GetChemical(3)->SetInitialAmount(10000);
00244         _minimizableList[1]->GetExperiment()->GetReactionNetwork()->GetChemical(4)->SetInitialAmount(72000.0);
00245         _minimizableList[1]->GetExperiment()->GetReactionNetwork()->ChemicalReset();
00246 
00247         // 82,000 receptors, EGF at 10 nM
00248         _minimizableList[2]->GetExperiment()->GetReactionNetwork()->RemoveReaction(98);
00249         _minimizableList[2]->GetExperiment()->GetReactionNetwork()->GetChemical(0)->SetInitialAmount(7.2e+006);
00250         _minimizableList[2]->GetExperiment()->GetReactionNetwork()->GetChemical(3)->SetInitialAmount(10000);
00251         _minimizableList[2]->GetExperiment()->GetReactionNetwork()->GetChemical(4)->SetInitialAmount(72000);
00252         _minimizableList[2]->GetExperiment()->GetReactionNetwork()->ChemicalReset();
00253 
00254         // 82,000 receptors, EGF at 1 nM
00255         _minimizableList[3]->GetExperiment()->GetReactionNetwork()->RemoveReaction(98);
00256         _minimizableList[3]->GetExperiment()->GetReactionNetwork()->GetChemical(0)->SetInitialAmount(7.2e+005);
00257         _minimizableList[3]->GetExperiment()->GetReactionNetwork()->GetChemical(3)->SetInitialAmount(10000);
00258         _minimizableList[3]->GetExperiment()->GetReactionNetwork()->GetChemical(4)->SetInitialAmount(72000);
00259         _minimizableList[3]->GetExperiment()->GetReactionNetwork()->ChemicalReset();
00260         
00261         /*
00262         // 80,000 receptors assumed, EGF at 100 ng/ml
00263         _minimizableList[4]->GetExperiment()->GetReactionNetwork()->RemoveReaction(98);
00264         _minimizableList[4]->GetExperiment()->GetReactionNetwork()->GetChemical(0)->SetInitialAmount(1.2e+007);
00265         _minimizableList[4]->GetExperiment()->GetReactionNetwork()->GetChemical(3)->SetInitialAmount(10000);
00266         _minimizableList[4]->GetExperiment()->GetReactionNetwork()->GetChemical(4)->SetInitialAmount(70000);
00267         _minimizableList[4]->GetExperiment()->GetReactionNetwork()->ChemicalReset();
00268         */
00269         /*
00270         // 80,000 receptors assumed, EGF at 100 ng/ml, Y1045F can't bind to Cbl
00271         _minimizableList[5]->GetExperiment()->GetReactionNetwork()->GetChemical(0)->SetInitialAmount(1.2e+007);
00272         _minimizableList[5]->GetExperiment()->GetReactionNetwork()->GetChemical(2)->SetInitialAmount(80000);
00273         _minimizableList[5]->GetExperiment()->GetReactionNetwork()->ChemicalReset();
00274         _minimizableList[5]->GetExperiment()->GetReactionNetwork()->RemoveReaction(40);
00275         _minimizableList[5]->GetExperiment()->GetReactionNetwork()->RemoveReaction(40);
00276         
00277         
00278         // 80,000 receptors assumed, EGF at 100 ng/ml, K721A is kinase dead
00279         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->GetChemical(0)->SetInitialAmount(1.2e+007);
00280         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->GetChemical(2)->SetInitialAmount(80000);
00281         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->ChemicalReset();
00282         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->RemoveReaction(2);
00283         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->RemoveReaction(3);
00284         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->RemoveReaction(4);
00285         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->RemoveReaction(5);
00286         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->RemoveReaction(6);
00287         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->RemoveReaction(7);
00288         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->RemoveReaction(8);
00289         _minimizableList[6]->GetExperiment()->GetReactionNetwork()->RemoveReaction(31);
00290         */
00291 }
00292 
00293 double CEGFRTraffickingCHOMinimizableDirector::GetParameter(int parIndex)
00294 {
00295         // for this director, everyone has the same parameters so any set 
00296         // can be returned
00297         return _minimizableList[0]->GetParameter(parIndex);
00298 }
00299 
00300 int CEGFRTraffickingCHOMinimizableDirector::GetNParameters()
00301 {
00302         // ditto as for GetParameters()
00303         // return _minimizableList[0]->GetNParameters();
00304         return _minimizableList[0]->GetNParameters() + 8;
00305 }
00306 
00307 double CEGFRTraffickingCHOMinimizableDirector::ComputeResiduals(double *parameters)
00308 {
00309         int lastrc = _minimizableList[0]->GetNParameters();
00310         // set all the ICs first
00311         for(int i = 0; i < 4; i++)
00312         {
00313                 _minimizableList[i]->GetExperiment()->GetReactionNetwork()->GetChemical(46)->SetInitialAmount(parameters[lastrc]);
00314                 _minimizableList[i]->GetExperiment()->GetReactionNetwork()->GetChemical(47)->SetInitialAmount(parameters[lastrc+1]);
00315                 _minimizableList[i]->GetExperiment()->GetReactionNetwork()->GetChemical(49)->SetInitialAmount(parameters[lastrc+2]);
00316                 _minimizableList[i]->GetExperiment()->GetReactionNetwork()->GetChemical(53)->SetInitialAmount(parameters[lastrc+3]);
00317                 _minimizableList[i]->GetExperiment()->GetReactionNetwork()->GetChemical(54)->SetInitialAmount(parameters[lastrc+4]);
00318                 _minimizableList[i]->GetExperiment()->GetReactionNetwork()->GetChemical(55)->SetInitialAmount(parameters[lastrc+5]);
00319                 _minimizableList[i]->GetExperiment()->GetReactionNetwork()->GetChemical(71)->SetInitialAmount(parameters[lastrc+6]);    
00320                 _minimizableList[i]->GetExperiment()->GetReactionNetwork()->ChemicalReset();
00321 
00322         }
00323         // now fire away
00324         // just sum up the costs from each director and then add the 
00325         // integration penalty, = (gamma^2/2)*(sum_ti)^2
00326         double summedCost = 0.0;
00327         double timeCost = 0.0;
00328         double pcost = 0.0;
00329         
00330         // compute the residuals/cost from each minimizable
00331         for(i = 0; i < _minimizableList.size(); i++)
00332         {
00333                 summedCost += _minimizableList[i]->ObjectiveFunction(parameters);
00334         }
00335         
00336         // load up the master list with each set of residuals
00337         int nextResidual = 0;
00338         for(i = 0; i < _minimizableList.size(); i++)
00339         {
00340                 int nRes = _minimizableList[i]->GetNResiduals();
00341                 for(int j = 0; j < nRes; j++)
00342                 {
00343                         residuals[nextResidual+j] = (_minimizableList[i]->GetResiduals())[j];
00344                 }
00345                 nextResidual += nRes;
00346         }
00347         
00348         // now compute the cost for deviation from guessed values, if any
00349         for(i = 0; i < m_iPriorList.size(); i++)
00350         {
00351                 int which = m_iPriorList[i];
00352                 //double presid = parameters[which] - _networkList[0]->GetRateConstant(which)->GetInitialValue();
00353                 double presid = _networkList[0]->GetRateConstant(which)->GetValue() - _networkList[0]->GetRateConstant(which)->GetInitialValue();
00354                 presid = presid/(sqrt(2.0)*_networkList[0]->GetRateConstant(which)->GetErrorInInitialValue());
00355                 residuals[nextResidual + i] = presid;
00356                 pcost += presid*presid;
00357         }
00358         nextResidual += m_iPriorList.size();
00359         
00360         // if the prefactor for the integration cost is > 0, compute its residual (otherwise
00361         // that residual shouldn't exist)
00362         if(m_dGammaSquared > 0.0)
00363         {
00364                 double totalSteps = 0.0;
00365                 for(int l = 0; l < _moverList.size(); l++)
00366                 {
00367                         totalSteps += _moverList[l]->GetTotalStepsTaken();
00368                 }
00369                 residuals[nResiduals - 1] = sqrt(m_dGammaSquared/2.0)*totalSteps;
00370                 timeCost = residuals[nResiduals - 1]*residuals[nResiduals - 1];
00371         }
00372         return summedCost + pcost + timeCost;
00373 }