diff --git a/addons/libRoadrunner/src/librr_intracellular.cpp b/addons/libRoadrunner/src/librr_intracellular.cpp index 1765c9532..c1bec4c9b 100644 --- a/addons/libRoadrunner/src/librr_intracellular.cpp +++ b/addons/libRoadrunner/src/librr_intracellular.cpp @@ -8,59 +8,44 @@ RoadRunnerIntracellular::RoadRunnerIntracellular() : Intracellular() intracellular_type = "sbml"; std::cout << "====== " << __FUNCTION__ << "() intracellular_type=" << intracellular_type << std::endl; std::cout << "====== " << __FUNCTION__ << "() sbml_filename = " << sbml_filename << std::endl; - // initial_values.clear(); - // mutations.clear(); parameters.clear(); } + // constructor using XML node RoadRunnerIntracellular::RoadRunnerIntracellular(pugi::xml_node& node) { - // std::cout << "======rwh " << __FUNCTION__ << ": node.name() =" << node.name() << std::endl; intracellular_type = "roadrunner"; initialize_intracellular_from_pugixml(node); - // std::cout << "======rwh " << __FUNCTION__ << "(node) intracellular_type=" << intracellular_type << std::endl; - // std::cout << "======rwh " << __FUNCTION__ << "(node) sbml_filename = " << sbml_filename << std::endl; - // std::cout << "======rwh " << __FUNCTION__ << "(node) this=" << this << std::endl; - // std::cout << "======rwh " << __FUNCTION__ << "(node) this->sbml_filename=" << this->sbml_filename << std::endl; } -// Intracellular* RoadRunnerIntracellular::clone() // --> 'Intracellular' does not name a type -// { -// return static_cast(new RoadRunnerIntracellular(this)); -// } - -// rwh: review this RoadRunnerIntracellular::RoadRunnerIntracellular(RoadRunnerIntracellular* copy) { + update_time_step = copy->update_time_step; + previous_update_time = PhysiCell::PhysiCell_globals.current_time; + next_librr_run = PhysiCell::PhysiCell_globals.current_time + update_time_step; intracellular_type = copy->intracellular_type; sbml_filename = copy->sbml_filename; - // cfg_filename = copy->cfg_filename; - // time_step = copy->time_step; - // discrete_time = copy->discrete_time; - // time_tick = copy->time_tick; - // scaling = copy->scaling; - // initial_values = copy->initial_values; - // mutations = copy->mutations; parameters = copy->parameters; - } -// Parse the info in the .xml for (possibly) each , e.g. -// -// ./config/Toy_SBML_Model_2.xml -// 1 -// Oxy -// Glc -// Energy void RoadRunnerIntracellular::initialize_intracellular_from_pugixml(pugi::xml_node& node) { + update_time_step = PhysiCell::intracellular_dt; // default to this, but overwrite below if defined in XML + pugi::xml_node node_sbml = node.child( "sbml_filename" ); if ( node_sbml ) { sbml_filename = PhysiCell::xml_get_my_string_value (node_sbml); std::cout << "\n------------- " << __FUNCTION__ << ": sbml_filename = " << sbml_filename << std::endl; } + + pugi::xml_node node_update_time_step = node.child( "intracellular_dt" ); + if ( node_update_time_step ) + { + update_time_step = PhysiCell::xml_get_my_double_value (node_update_time_step); + std::cout << "\n------------- " << __FUNCTION__ << ": intracellular_dt = " << update_time_step << std::endl; + } pugi::xml_node node_species = node.child( "map" ); while( node_species ) @@ -70,7 +55,6 @@ void RoadRunnerIntracellular::initialize_intracellular_from_pugixml(pugi::xml_no std::string substrate_name = node_species.attribute( "PC_substrate" ).value(); if( substrate_name != "" ) { - //std::cout << "-----------" << node_species.attribute( "sbml_species" ).value() << std::endl; std::string species_name = node_species.attribute( "sbml_species" ).value(); substrate_species[substrate_name] = species_name; std::cout << "\n------------- " << __FUNCTION__ << ": species_name= " << species_name << std::endl; @@ -81,7 +65,6 @@ void RoadRunnerIntracellular::initialize_intracellular_from_pugixml(pugi::xml_no { std::string species_name = node_species.attribute( "sbml_species" ).value(); custom_data_species[custom_data_name] = species_name; - // std::cout << "\n------------- " << __FUNCTION__ << ": species_name= " << species_name << std::endl; } @@ -92,7 +75,6 @@ void RoadRunnerIntracellular::initialize_intracellular_from_pugixml(pugi::xml_no { std::string species_name = node_species.attribute( "sbml_species" ).value(); phenotype_species[phenotype_name] = species_name; - // std::cout << "\n------------- " << __FUNCTION__ << ": species_name= " << species_name << std::endl; } node_species = node_species.next_sibling( "map" ); @@ -125,170 +107,61 @@ void RoadRunnerIntracellular::start() // called when a new cell is created; creates the unique 'rrHandle' rrc::RRVectorPtr vptr; - //std::cout << "\n------------ " << __FUNCTION__ << ": librr_intracellular.cpp: start() called\n"; - // this->enabled = true; - - //std::cout << "\n------------ " << __FUNCTION__ << ": doing: rrHandle = createRRInstance()\n"; - rrHandle = createRRInstance(); - //std::cout << "\n------------ " << __FUNCTION__ << ": rrHandle = " << rrHandle << std::endl; - - // if (!rrc::loadSBML (rrHandle, get_cell_definition("lung epithelium").sbml_filename.c_str())) - //std::cout << " sbml_filename = " << sbml_filename << std::endl; - - // TODO: don't hard-code name if ( !rrc::loadSBML(rrHandle, (sbml_filename).c_str() ) ) - // std::cout << " for now, hard-coding sbml_file = ./config/Toy_SBML_Model_1.xml" << std::endl; - // if (!rrc::loadSBML(rrHandle, "./config/Toy_SBML_Model_1.xml") ) { std::cerr << "------------->>>>> Error while loading SBML file <-------------\n\n"; - // return -1; - // printf ("Error message: %s\n", getLastError()); - exit (0); + exit(-1); } - // std::cout << " rrHandle=" << rrHandle << std::endl; - int r = rrc::getNumberOfReactions(rrHandle); int m = rrc::getNumberOfFloatingSpecies(rrHandle); int b = rrc::getNumberOfBoundarySpecies(rrHandle); int p = rrc::getNumberOfGlobalParameters(rrHandle); int c = rrc::getNumberOfCompartments(rrHandle); - - //std::cerr << "Number of reactions = " << r << std::endl; - //std::cerr << "Number of floating species = " << m << std::endl; // 4 - //std::cerr << "Number of boundary species = " << b << std::endl; // 0 - //std::cerr << "Number of compartments = " << c << std::endl; // 1 - - //std::cerr << "Floating species names:\n"; - //std::cerr << "-----------------------\n"; std::string species_names_str = stringArrayToString(rrc::getFloatingSpeciesIds(rrHandle)); - //std::cerr << species_names_str <<"\n"<< std::endl; std::stringstream iss(species_names_str); std::string species_name; int idx = 0; while (iss >> species_name) { species_result_column_index[species_name] = idx; - //std::cout << species_name << " -> " << idx << std::endl; idx++; } vptr = rrc::getFloatingSpeciesConcentrations(rrHandle); - //std::cerr << vptr->Count << std::endl; -/* for (int kdx=0; kdxCount; kdx++) - { - std::cerr << kdx << ") " << vptr->Data[kdx] << std::endl; - } */ - //std::cerr << "---------- end start() -------------\n"; rrc::freeVector(vptr); - // return 0; } bool RoadRunnerIntracellular::need_update() { - return PhysiCell::PhysiCell_globals.current_time >= this->next_librr_run; + return PhysiCell::PhysiCell_globals.current_time >= this->next_librr_run - 0.5 * PhysiCell::diffusion_dt; } // solve the intracellular model void RoadRunnerIntracellular::update() { static double start_time = 0.0; - static double end_time = 0.01; - // static double end_time = 10.0; - // static int num_vals = 1; - // static int num_vals = 10; - static int num_vals = 2; - - // result = rrc::simulateEx (pCell->phenotype.molecular.model_rr, 0, 10, 10); // start time, end time, and number of points - //std::cout << __FUNCTION__ << " ----- update(): this=" << this << std::endl; - //std::cout << __FUNCTION__ << " ----- update(): rrHandle=" << this->rrHandle << std::endl; + static int num_vals = 2; // start time and end time - // if (this->result != 0) // apparently not necessary (done in freeRRCData hopefully) rrc::freeRRCData (this->result); - this->result = rrc::simulateEx (this->rrHandle, start_time, end_time, num_vals); // start time, end time, and number of points - - - // this->next_librr_run += this->rrHandle.get_time_to_update(); - // std::cout << "----- update(): result=" << result << std::endl; - //std::cout << "----- update(): result->CSize=" << this->result->CSize << std::endl; - //std::cout << "----- update(): result->RSize=" << this->result->RSize << std::endl; // should be = num_vals - // std::cout << "----- update(): result->ColumnHeaders[0]=" << result->ColumnHeaders[0] << std::endl; // = "time" - - // debug - does it generate expected data? - int index = 0; - // Print out column headers... typically time and species. - for (int col = 0; col < this->result->CSize; col++) - { - // std::cout << result->ColumnHeaders[index++]; - // std::cout << std::left << std::setw(15) << result->ColumnHeaders[index++]; - //std::cout << std::left << this->result->ColumnHeaders[index++]; - // if (col < result->CSize - 1) - // { - // // std::cout << "\t"; - // std::cout << " "; - // } - } - //std::cout << "\n"; - - index = 0; - // Print out the data - for (int row = 0; row < this->result->RSize; row++) - { - for (int col = 0; col < this->result->CSize; col++) - { - // std::cout << result->Data[index++]; - //std::cout << std::left << std::setw(15) << this->result->Data[index++]; - // if (col < result->CSize -1) - // { - // // std::cout << "\t"; - // std::cout << " "; - // } - } - // std::cout << "\n"; - } - // int idx = (result->RSize - 1) * result->CSize + 1; - // std::cout << "Saving last energy value (cell custom var) = " << result->Data[idx] << std::endl; - // pCell->custom_data[energy_cell_idx] = result->Data[idx]; - - // return 0; + this->result = rrc::simulateEx (this->rrHandle, start_time, PhysiCell::PhysiCell_globals.current_time - previous_update_time, num_vals); // start time, end time, and number of points + this->previous_update_time = PhysiCell::PhysiCell_globals.current_time; + this->next_librr_run = PhysiCell::PhysiCell_globals.current_time + update_time_step; } double RoadRunnerIntracellular::get_parameter_value(std::string param_name) { rrc::RRVectorPtr vptr; - //std::cout << "-----------" << __FUNCTION__ << "-----------" << std::endl; - // std::cout << " substrate_name = " << substrate_name << std::endl; - //std::cout << " param_name = " << param_name << std::endl; - - // TODO: optimize this eventually - // std::map species_result_column_index; - // int num_columns = result->CSize; - // int offset = (num_rows_result_table-1) * result->CSize - 1; - // int offset = (num_rows_result_table-1) * result->CSize; - // offset += (num_rows_result_table-1) * result->CSize - 1; - - // int offset = species_result_column_index[name]; - // std::string species_name = this->substrate_species[substrate_name]; - // std::cout << " species_name = " << species_name << std::endl; - vptr = rrc::getFloatingSpeciesConcentrations(this->rrHandle); - //std::cerr << vptr->Count << std::endl; - for (int kdx=0; kdxCount; kdx++) - { - //std::cerr << kdx << ") " << vptr->Data[kdx] << std::endl; - } int offset = species_result_column_index[param_name]; - //std::cout << " result offset = "<< offset << std::endl; - // double res = this->result->Data[offset]; double res = vptr->Data[offset]; - //std::cout << " res = " << res << std::endl; rrc::freeVector(vptr); return res; } @@ -301,10 +174,8 @@ void RoadRunnerIntracellular::set_parameter_value(std::string species_name, doub vptr = rrc::getFloatingSpeciesConcentrations(this->rrHandle); int idx = species_result_column_index[species_name]; vptr->Data[idx] = value; - // rrc::setFloatingSpeciesConcentrations(pCell->phenotype.molecular.model_rr, vptr); rrc::setFloatingSpeciesConcentrations(this->rrHandle, vptr); rrc::freeVector(vptr); - // return 0; } RoadRunnerIntracellular* getRoadRunnerModel(PhysiCell::Phenotype& phenotype) { @@ -314,7 +185,6 @@ RoadRunnerIntracellular* getRoadRunnerModel(PhysiCell::Phenotype& phenotype) { void RoadRunnerIntracellular::save_libRR(std::string path, std::string index) { std::string state_file_name = path + "/states_" + index + ".dat"; -// std::string state_file_name = path + "/states_" + index + ".csv"; std::ofstream state_file( state_file_name ); state_file << "--------- dummy output from save_libRR ---------" << std::endl; state_file << "ID,state" << std::endl; @@ -386,10 +256,7 @@ int RoadRunnerIntracellular::update_phenotype_parameters(PhysiCell::Phenotype& p //uptake rate if (elm.first.substr(0,3) == "sur") { - //std::cout << sub_index << std::endl; - //std::cout << "Before sur1 : " << phenotype.secretion.uptake_rates[sub_index] << std::endl; phenotype.secretion.uptake_rates[1] = phenotype.intracellular->get_parameter_value(elm.second); - //std::cout << "After sur1 : " << phenotype.secretion.uptake_rates[sub_index] << std::endl; } //secretion rate else if (elm.first.substr(0,3) == "ssr") @@ -426,7 +293,6 @@ int RoadRunnerIntracellular::update_phenotype_parameters(PhysiCell::Phenotype& p int start_index; while ((pos = s.find(delimiter)) != std::string::npos) { token = s.substr(0, pos); - //std::cout << counter << " : "<< token << std::endl; if (counter == 1) { start_index = atoi( token.c_str() ); @@ -435,8 +301,6 @@ int RoadRunnerIntracellular::update_phenotype_parameters(PhysiCell::Phenotype& p counter += 1; } int end_index = atoi( s.c_str() ); - //std::cout << "START INDEX : " << start_index << std::endl; - //std::cout << "END INDEX : " << end_index << std::endl; phenotype.cycle.data.transition_rate(start_index,end_index) = phenotype.intracellular->get_parameter_value(elm.second); } else @@ -468,7 +332,6 @@ int RoadRunnerIntracellular::update_phenotype_parameters(PhysiCell::Phenotype& p } } - //std::cout << std::endl; return 0; } @@ -477,9 +340,6 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe { for(auto elm : phenotype_species) { - //std::cout << "PhysiCell_token_validation" << std::endl; - //std::cout << elm.first << " : " << elm.second << std::endl; - // motility params if (elm.first[0] == 'm') { @@ -498,8 +358,7 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe std::cout << "ERROR: There is no specified token parameters in the name of \"" << elm.first << "\" at motility parameters. Please take a look token specifications." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } } // death params @@ -517,8 +376,7 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe std::cout << "ERROR: There is no specified token parameters in the name of \"" << elm.first << "\" at death parameters. Please take a look token specifications." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } } // secretion params @@ -534,15 +392,13 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe s.erase(0, pos + delimiter.length()); } int sub_index = microenvironment.find_density_index(s); - //std::cout << "SUBSTRATE_INDEX = : " << sub_index << std::endl; if ( sub_index < 0 ) { std::cout<< std::endl; std::cout << "ERROR: There is no substrate named in the name of \"" << s << "\" at microenvironment. Please take a look token specifications." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } if (elm.first.substr(0,3) == "sur") @@ -563,8 +419,7 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe std::cout << "ERROR: There is no specified token parameters in the name of \"" << elm.first << "\" at secretion parameters. Please take a look token specifications." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } } else if (elm.first[0] == 'c') @@ -573,7 +428,6 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe { // getting num of phases int num_of_phases = (&(phenotype.cycle.model()))->phases.size(); - //std::cout << num_of_phases << std::endl; // getting start and end indices std::string s = elm.first; @@ -600,8 +454,7 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe std::cout << "ERROR: Given transition start index is beyond cycle indices. Please double check it." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } if ( end_index > num_of_phases - 1) { @@ -609,8 +462,7 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe std::cout << "ERROR: Given transition end index is beyond cycle indices. Please double check it." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } } else @@ -619,8 +471,7 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe std::cout << "ERROR: There is no specified token parameters in the name of \"" << elm.first << "\" at cycle parameters. Please take a look token specifications." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } } @@ -641,8 +492,7 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe std::cout << "ERROR: There is no specified token parameters in the name of \"" << elm.first << "\" at volume parameters. Please take a look token specifications." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } } else @@ -651,8 +501,7 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe std::cout << "ERROR: There is no specified token parameters in the name of \"" << elm.first << "\" at phenotypic parameters. Please take a look token specifications." << std::endl; std::cout<< std::endl; std::cout<< std::endl; - exit (-1); - return -1; + exit(-1); } } @@ -663,14 +512,12 @@ int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phe int RoadRunnerIntracellular::validate_SBML_species() { - //std::cout << "---------VALIDATING_SBML_SPECIES START-------" << std::endl; - // reading SBML rrHandle = createRRInstance(); if ( !rrc::loadSBML(rrHandle, (sbml_filename).c_str() ) ) { std::cerr << "------------->>>>> Error while loading SBML file <-------------\n\n"; - return -1; + exit(-1); } // getting Species Names std::string species_names_str = stringArrayToString(rrc::getFloatingSpeciesIds(rrHandle)); @@ -684,7 +531,6 @@ int RoadRunnerIntracellular::validate_SBML_species() { species_result_column_index[species_name] = idx; all_species.push_back(species_name); - //std::cout << species_name << " -> " << idx << std::endl; idx++; } @@ -692,14 +538,10 @@ int RoadRunnerIntracellular::validate_SBML_species() for (auto elm : phenotype_species) { bool exist = 0; - // std::cout << species_name.size() << std::endl; for (int i=0; i < all_species.size(); i++) { - //std::cout << all_species[i] << std::endl;; - //std::cout << "Comparing " << all_species[i] << " with " << elm.second << std::endl; if ( all_species[i] == elm.second ) { - //std::cout << "And they are the same..... " < initial_values; std::map parameters; std::map substrate_species; std::map custom_data_species; std::map phenotype_species; std::map species_result_column_index; - // rrc::RRHandle rrHandle = createRRInstance(); rrc::RRHandle rrHandle; - // rrc::RRHandle rrHandle; - // rrc::RRVectorPtr vptr; rrc::RRCDataPtr result = 0; // start time, end time, and number of points - double next_librr_run = 0; + double update_time_step = 0.01; + double previous_update_time = 0.0; + double next_librr_run = 0.0; RoadRunnerIntracellular(); @@ -61,12 +49,9 @@ class RoadRunnerIntracellular : public PhysiCell::Intracellular RoadRunnerIntracellular(RoadRunnerIntracellular* copy); - // rwh: review this Intracellular* clone() { - // return static_cast(new RoadRunnerIntracellular(this)); RoadRunnerIntracellular* clone = new RoadRunnerIntracellular(this); - clone->sbml_filename = this->sbml_filename; clone->substrate_species = this->substrate_species; clone->phenotype_species = this->phenotype_species; clone->custom_data_species = this->custom_data_species; diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index a14c70b99..b4f1a5c8f 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -3144,18 +3144,16 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node // If it has already be copied if (pParent != NULL && pParent->phenotype.intracellular != NULL) { - // std::cout << "------ " << __FUNCTION__ << ": copying another\n"; pCD->phenotype.intracellular->initialize_intracellular_from_pugixml(node); } // Otherwise we need to create a new one else { - std::cout << "\n------ " << __FUNCTION__ << ": creating new RoadRunnerIntracellular\n"; RoadRunnerIntracellular* pIntra = new RoadRunnerIntracellular(node); pCD->phenotype.intracellular = pIntra->getIntracellularModel(); - pCD->phenotype.intracellular->validate_PhysiCell_tokens(pCD->phenotype); - pCD->phenotype.intracellular->validate_SBML_species(); } + pCD->phenotype.intracellular->validate_PhysiCell_tokens(pCD->phenotype); + pCD->phenotype.intracellular->validate_SBML_species(); } #endif diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index 0c5850b1d..6d6cfcd69 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -149,12 +149,16 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me if( (*all_cells)[i]->phenotype.intracellular != NULL && (*all_cells)[i]->phenotype.intracellular->need_update()) { if ((*all_cells)[i]->functions.pre_update_intracellular != NULL) - (*all_cells)[i]->functions.pre_update_intracellular( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); + { + (*all_cells)[i]->functions.pre_update_intracellular((*all_cells)[i], (*all_cells)[i]->phenotype, diffusion_dt_); + } (*all_cells)[i]->phenotype.intracellular->update( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); if ((*all_cells)[i]->functions.post_update_intracellular != NULL) - (*all_cells)[i]->functions.post_update_intracellular( (*all_cells)[i], (*all_cells)[i]->phenotype , diffusion_dt_ ); + { + (*all_cells)[i]->functions.post_update_intracellular((*all_cells)[i], (*all_cells)[i]->phenotype, diffusion_dt_); + } } } } diff --git a/sample_projects_intracellular/ode/ode_energy/config/PhysiCell_settings.xml b/sample_projects_intracellular/ode/ode_energy/config/PhysiCell_settings.xml index f0a40f2f4..ebee4b92a 100644 --- a/sample_projects_intracellular/ode/ode_energy/config/PhysiCell_settings.xml +++ b/sample_projects_intracellular/ode/ode_energy/config/PhysiCell_settings.xml @@ -90,7 +90,8 @@ 1440 min micron - + + 0.01 0.01 0.1 6 @@ -307,7 +308,7 @@ ./config/Toy_Metabolic_Model.xml - 0.01 + 0.1 diff --git a/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.cpp b/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.cpp index acbd054a6..3267a9e51 100644 --- a/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.cpp +++ b/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.cpp @@ -179,80 +179,66 @@ void setup_tissue( void ) set_single_behavior( pCell , "custom:intra_lac" , parameters.doubles("initial_internal_lactate")); set_single_behavior( pCell , "custom:intra_energy" , parameters.doubles("initial_energy")); -/* pCell->custom_data[i_Oxy_i] = parameters.doubles("initial_internal_oxygen"); - pCell->custom_data[i_Glu_i] = parameters.doubles("initial_internal_glucose"); - pCell->custom_data[i_Lac_i] = parameters.doubles("initial_internal_lactate"); - pCell->custom_data[energy_vi] = parameters.doubles("initial_energy"); */ - double cell_volume = pCell->phenotype.volume.total; - //std::cout << "oxygen custom data : " << pCell->custom_data[i_Oxy_i] << std::endl; - //std::cout << "oxygen custom data : SIGNAL" << get_single_signal( pCell, "custom:intra_oxy") << std::endl; - - - set_single_behavior( pCell , "custom:intra_oxy" , parameters.doubles("initial_internal_oxygen")); - - pCell->phenotype.molecular.internalized_total_substrates[oxygen_substrate_index]= get_single_signal( pCell, "custom:intra_oxy") * cell_volume; pCell->phenotype.molecular.internalized_total_substrates[glucose_substrate_index]= get_single_signal( pCell, "custom:intra_glu") * cell_volume; pCell->phenotype.molecular.internalized_total_substrates[lactate_substrate_index]= get_single_signal( pCell, "custom:intra_lac") * cell_volume; - pCell->phenotype.intracellular->start(); - (*all_cells)[i]->phenotype.intracellular->set_parameter_value("Energy",get_single_signal( pCell, "custom:intra_energy")); - + pCell->phenotype.intracellular->set_parameter_value("Energy",get_single_signal( pCell, "custom:intra_energy")); + + pCell->functions.pre_update_intracellular = pre_update_intracellular; + pCell->functions.post_update_intracellular = post_update_intracellular; } return; } -void update_intracellular() +void pre_update_intracellular( Cell* pCell, Phenotype& phenotype, double dt ) { // BioFVM Indices static int oxygen_substrate_index = microenvironment.find_density_index( "oxygen" ); static int glucose_substrate_index = microenvironment.find_density_index( "glucose" ); static int lactate_substrate_index = microenvironment.find_density_index( "lactate"); - #pragma omp parallel for - for( int i=0; i < (*all_cells).size(); i++ ) - { - if( (*all_cells)[i]->is_out_of_domain == false ) - { - // Cell Volume - double cell_volume = (*all_cells)[i]->phenotype.volume.total; - - // Get Intracellular Concentrations - double oxy_val_int = get_single_signal((*all_cells)[i], "intracellular oxygen"); - double glu_val_int = get_single_signal((*all_cells)[i], "intracellular glucose"); - double lac_val_int = get_single_signal((*all_cells)[i], "intracellular lactate"); - - // Update SBML - (*all_cells)[i]->phenotype.intracellular->set_parameter_value("Oxygen",oxy_val_int); - (*all_cells)[i]->phenotype.intracellular->set_parameter_value("Glucose",glu_val_int); - (*all_cells)[i]->phenotype.intracellular->set_parameter_value("Lactate",lac_val_int); - - // SBML Simulation - (*all_cells)[i]->phenotype.intracellular->update(); - - // Phenotype Simulation - (*all_cells)[i]->phenotype.intracellular->update_phenotype_parameters((*all_cells)[i]->phenotype); - - // Internalized Chemical Update After SBML Simulation - (*all_cells)[i]->phenotype.molecular.internalized_total_substrates[oxygen_substrate_index] = (*all_cells)[i]->phenotype.intracellular->get_parameter_value("Oxygen") * cell_volume; - (*all_cells)[i]->phenotype.molecular.internalized_total_substrates[glucose_substrate_index] = (*all_cells)[i]->phenotype.intracellular->get_parameter_value("Glucose") * cell_volume; - (*all_cells)[i]->phenotype.molecular.internalized_total_substrates[lactate_substrate_index] = (*all_cells)[i]->phenotype.intracellular->get_parameter_value("Lactate") * cell_volume; - - - //Save custom data - set_single_behavior( (*all_cells)[i] , "custom:intra_oxy" , (*all_cells)[i]->phenotype.intracellular->get_parameter_value("Oxygen") ); - set_single_behavior( (*all_cells)[i] , "custom:intra_glu" , (*all_cells)[i]->phenotype.intracellular->get_parameter_value("Glucose") ); - set_single_behavior( (*all_cells)[i] , "custom:intra_lac" , (*all_cells)[i]->phenotype.intracellular->get_parameter_value("Lactate") ); - set_single_behavior( (*all_cells)[i] , "custom:intra_energy" , (*all_cells)[i]->phenotype.intracellular->get_parameter_value("Energy") ); - - } - } - + // Cell Volume + double cell_volume = pCell->phenotype.volume.total; + + // Get Intracellular Concentrations + double oxy_val_int = get_single_signal(pCell, "intracellular oxygen"); + double glu_val_int = get_single_signal(pCell, "intracellular glucose"); + double lac_val_int = get_single_signal(pCell, "intracellular lactate"); + + // Update SBML + pCell->phenotype.intracellular->set_parameter_value("Oxygen",oxy_val_int); + pCell->phenotype.intracellular->set_parameter_value("Glucose",glu_val_int); + pCell->phenotype.intracellular->set_parameter_value("Lactate",lac_val_int); } +void post_update_intracellular( Cell* pCell, Phenotype& phenotype, double dt ) +{ + // BioFVM Indices + static int oxygen_substrate_index = microenvironment.find_density_index( "oxygen" ); + static int glucose_substrate_index = microenvironment.find_density_index( "glucose" ); + static int lactate_substrate_index = microenvironment.find_density_index( "lactate"); + + // Cell Volume + double cell_volume = pCell->phenotype.volume.total; + // Phenotype Simulation + pCell->phenotype.intracellular->update_phenotype_parameters(pCell->phenotype); + + // Internalized Chemical Update After SBML Simulation + pCell->phenotype.molecular.internalized_total_substrates[oxygen_substrate_index] = pCell->phenotype.intracellular->get_parameter_value("Oxygen") * cell_volume; + pCell->phenotype.molecular.internalized_total_substrates[glucose_substrate_index] = pCell->phenotype.intracellular->get_parameter_value("Glucose") * cell_volume; + pCell->phenotype.molecular.internalized_total_substrates[lactate_substrate_index] = pCell->phenotype.intracellular->get_parameter_value("Lactate") * cell_volume; + + + //Save custom data + set_single_behavior( pCell , "custom:intra_oxy" , pCell->phenotype.intracellular->get_parameter_value("Oxygen") ); + set_single_behavior( pCell , "custom:intra_glu" , pCell->phenotype.intracellular->get_parameter_value("Glucose") ); + set_single_behavior( pCell , "custom:intra_lac" , pCell->phenotype.intracellular->get_parameter_value("Lactate") ); + set_single_behavior( pCell , "custom:intra_energy" , pCell->phenotype.intracellular->get_parameter_value("Energy") ); +} std::vector my_coloring_function( Cell* pCell ) { diff --git a/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.h b/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.h index 8062d09b0..712054b6f 100644 --- a/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.h +++ b/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.h @@ -82,12 +82,8 @@ void setup_microenvironment( void ); // custom pathology coloring function std::vector my_coloring_function( Cell* ); -void update_intracellular(); +void pre_update_intracellular( Cell* pCell, Phenotype& phenotype, double dt ); +void post_update_intracellular( Cell* pCell, Phenotype& phenotype, double dt ); // custom functions can go here - -void predator_hunting_function( Cell* pCell, Phenotype& phenotype, double dt ); -void predator_cycling_function( Cell* pCell, Phenotype& phenotype, double dt ); - -void prey_cycling_function( Cell* pCell , Phenotype& phenotype, double dt ); std::vector> create_cell_circle_positions(double cell_radius, double sphere_radius); \ No newline at end of file diff --git a/sample_projects_intracellular/ode/ode_energy/main.cpp b/sample_projects_intracellular/ode/ode_energy/main.cpp index 828ce4dc9..d791331c4 100644 --- a/sample_projects_intracellular/ode/ode_energy/main.cpp +++ b/sample_projects_intracellular/ode/ode_energy/main.cpp @@ -161,14 +161,6 @@ int main( int argc, char* argv[] ) report_file<<"simulated time\tnum cells\tnum division\tnum death\twall time"<update_all_cells( PhysiCell_globals.current_time ); @@ -226,17 +210,6 @@ int main( int argc, char* argv[] ) Custom add-ons could potentially go here. */ - double time_since_last_intracellular = PhysiCell_globals.current_time - last_intracellular_time; - - //update_intracellular(); - - if( PhysiCell_globals.current_time >= next_intracellular_update ) - { - update_intracellular(); - - next_intracellular_update += intracellular_dt; - } - PhysiCell_globals.current_time += diffusion_dt; }