From ac2cf488923f0ea709e7e9c5a5f4e8785762884e Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Tue, 11 Feb 2025 11:34:32 -0500 Subject: [PATCH 1/2] fix ode-energy-sample and clean/fix librr stuff - introduce `update_time_step` in `RoadRunnerIntracellular` - set for each instance of `RoadRunnerIntracellular` - this controls the timestep for the model - will be set to `intracellular_dt` if not specified in the xml for the given cell type - add the `intracellular_dt` element as a child of the `intracellular` element in the xml - could be set on a cell-by-cell basis if desired - introduce `previous_update_time` in `RoadRunnerIntracellular` - tracks the previous update time - used to determine how long to run the ODE when updating - set to `current_time` when copying/cloning - fix `next_librr_run` in `RoadRunnerIntracellular` - was not being updated and thus ODE models were running at every time step - check `current_time >= next_librr_run - 0.5 * diffusion_dt` to determine if an update is needed - cleaned up librr files - fixed and cleand ode-energy-sample - it was running the intracellular models in `main.cpp` and in the update cells function - it seemed to try to set intracellular_dt inside the intracellular element of the xml - messy custom.cpp and custom.h --- .../libRoadrunner/src/librr_intracellular.cpp | 242 +++--------------- .../libRoadrunner/src/librr_intracellular.h | 21 +- core/PhysiCell_cell.cpp | 2 - core/PhysiCell_cell_container.cpp | 8 +- .../ode_energy/config/PhysiCell_settings.xml | 5 +- .../ode/ode_energy/custom_modules/custom.cpp | 96 +++---- .../ode/ode_energy/custom_modules/custom.h | 8 +- .../ode/ode_energy/main.cpp | 27 -- 8 files changed, 87 insertions(+), 322 deletions(-) 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..ff080363f 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -3144,13 +3144,11 @@ 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); 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; } From 6d8ab2d1a7ac89ad070efbcceb661e64e14a5269 Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Wed, 12 Feb 2025 08:19:50 -0500 Subject: [PATCH 2/2] make sure to validate every librr model before, only new ones were validated. those "copied" from previously read defs were not Presumably, it was assumed that they fully inherited their parent's model, but as of now, they read the xml to get their model, so they need to be validated. --- core/PhysiCell_cell.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index ff080363f..b4f1a5c8f 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -3151,9 +3151,9 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node { 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