diff --git a/addons/libRoadrunner/src/librr_intracellular.cpp b/addons/libRoadrunner/src/librr_intracellular.cpp index 1765c9532..3422a4575 100644 --- a/addons/libRoadrunner/src/librr_intracellular.cpp +++ b/addons/libRoadrunner/src/librr_intracellular.cpp @@ -3,56 +3,528 @@ #include #include +void RoadRunnerMapping::initialize_mapping( void ) +{ + bool use_for_input = io_type == "input"; + + int behavior_ind = PhysiCell::find_behavior_index(physicell_name); + int signal_ind = PhysiCell::find_signal_index(physicell_name); + if (behavior_ind != -1) + { + physicell_dictionary_name = "behaviors"; + index = behavior_ind; + if (use_for_input) + { + value_map = [this](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(this->sbml_species, PhysiCell::get_single_behavior(pCell, this->index)); }; + } + else + { + value_map = [this](PhysiCell::Cell *pCell) + { PhysiCell::set_single_behavior(pCell, this->index, pCell->phenotype.intracellular->get_parameter_value(this->sbml_species)); }; + } + } + else if (signal_ind != -1) + { + physicell_dictionary_name = "signals"; + index = signal_ind; + if (use_for_input) + { + value_map = [this](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(this->sbml_species, PhysiCell::get_single_signal(pCell, this->index)); }; + } + else + { + value_map = select_signal_setter(physicell_name, sbml_species); + } + } + else if (is_physicell_phenotype_token(physicell_name)) + { + physicell_dictionary_name = "tokens"; + if (use_for_input) + { + value_map = select_phenotype_by_token_inputter(physicell_name, sbml_species); + } + else + { + value_map = select_phenotype_by_token_outputter(physicell_name, sbml_species); + } + } + else + { + std::cerr << "ERROR: " << physicell_name << " is not a valid entry for the libRoadRunner mapping." << std::endl + << " The available entries are: " << std::endl + << " - signals (see dictionaries.txt in output)" << std::endl + << " - behaviors (see dictionaries.txt in output)" << std::endl + << " - tokens: " << std::endl + << " - mms, mpt, mmb, da, dn, vtsc, vtsn, vff" << std::endl + << " - ctr__" << std::endl + << " - _ where prefix is one of: sur, ssr, ssd, ser" << std::endl << std::endl + << "You can also defined your own pre- and post-update functions in the custom.cpp file." << std::endl + << " - Set them using `pCD->functions.pre_update_intracellular = foo;` and `pCD->functions.post_update_intracellular = bar;`" << std::endl + << " - These functions should have the signature `void foo(Cell* pCell, Phenotype& phenotype, double dt)`" << std::endl; + exit(-1); + } + mapping_initialized = true; + return; +} + +// PhysiCell does not have an API for setting signals, but libRoadRunner does in limited cases +MappingFunction select_signal_setter(const std::string& name, const std::string& sbml_species) +{ + // if "intracellular " or "internalized ", set the internalized substrate amount + if (name.find("intracellular") == 0 || name.find("internalized") == 0) + { + size_t space_ind = name.find(" "); + if (space_ind != std::string::npos) + { + std::string substrate_name = name.substr(space_ind + 1, std::string::npos); + int density_index = microenvironment.find_density_index(substrate_name); + if (density_index != -1) + { + return [density_index, sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.molecular.internalized_total_substrates[density_index] = pCell->phenotype.intracellular->get_parameter_value(sbml_species) * pCell->phenotype.volume.total; }; + } + } + std::cerr << "ERROR: \"" << name << "\" is not a valid signal that can be set using libRoadRunner." << std::endl + << " Somehow, this signal was recognized by the signals dictionary but now it seems malformed." << std::endl + << " I honestly don't know what to say...or how to help you :/" << std::endl; + exit(-1); + } + + // if "volume", set the cell volume + else if (name == "volume") + { + std::cout << "WARNING: setting the volume using libRoadRunner will do so by rescaling ALL cell volumes, not just setting the total volume." + << " To only set the total volume (or to set other components of the volume), use the `pCell->functions.post_update_intracellular`." << std::endl; + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.volume.multiply_by_ratio(pCell->phenotype.intracellular->get_parameter_value(sbml_species) / pCell->phenotype.volume.total); }; + } + + // if "damage", set the cell damage + else if (name == "damage") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.cell_integrity.damage = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + + // if begins with "custom", check that it is custom data and set that + else if (name.find("custom") == 0) + { + std::cerr << "ERROR: custom data should be handled using the behavior dictionary, not the signal dictionary." << std::endl + << " It is simpler to use the behavior setters." << std::endl + << " How did you even end up here?" << std::endl; + exit(-1); + } + else + { + std::cerr << "ERROR: \"" << name << "\" is not a signal that can be set using libRoadRunner." << std::endl + << " The available signals are: " << std::endl + << " - intracellular " << std::endl + << " - internalized " << std::endl + << " - volume" << std::endl + << " - damage" << std::endl + << " - custom:" << std::endl; + exit(-1); + } +} + +bool is_physicell_phenotype_token(const std::string& name) +{ + if (name[0] == 'm') + { + std::cout << "WARNING: The token \"" << name << "\" can be replaced with the relevant behavior name. The mapping is as follows:" << std::endl + << " - mms -> migration speed" << std::endl + << " - mpt -> persistence time" << std::endl + << " - mmb -> migration bias" << std::endl << std::endl + << "Doing so will allow the libRoadRunner addon to validate that there are not competing mappings." << std::endl; + return name == "mms" || name == "mpt" || name == "mmb"; + } + else if (name[0] == 'd') + { + std::cout << "WARNING: The token \"" << name << "\" can be replaced with the relevant death model name. The mapping is as follows:" << std::endl + << " - da -> apoptosis" << std::endl + << " - dn -> necrosis" << std::endl << std::endl + << "Doing so will allow the libRoadRunner addon to validate that there are not competing mappings." << std::endl; + return name == "da" || name == "dn"; + } + else if (name[0] == 's') + { + std::cout << "WARNING: The token \"" << name << "\" can be replaced with the relevant substrate name. The mapping is as follows:" << std::endl + << " - sur -> uptake" << std::endl + << " - ssr -> secretion" << std::endl + << " - ssd -> secretion target" << std::endl + << " - ser -> export" << std::endl + << "Doing so will allow the libRoadRunner addon to validate that there are not competing mappings." << std::endl; + if (name.substr(0, 3) != "sur" && name.substr(0, 3) != "ssr" && name.substr(0, 3) != "ssd" && name.substr(0, 3) != "ser") + { + return false; + } + if (name[3] != '_') + { + return false; + } + return microenvironment.find_density_index(name.substr(4, std::string::npos)) != -1; + } + else if (name[0] == 'c') + { + std::cout << "WARNING: The token \"" << name << "\" can be replaced with the relevant cycle transition rate name. The mapping is as follows:" << std::endl + << " - ctr__ -> exit from cycle phase " << std::endl + << "Doing so will allow the libRoadRunner addon to validate that there are not competing mappings." << std::endl; + std::vector indices = parse_ctr_token(name); + return indices[0] >= 0 && indices[1] >= 0; + } + else if (name[0] == 'v') + { + return name == "vtsc" || name == "vtsn" || name == "vff"; + } + return false; +} + +std::vector parse_ctr_token(const std::string &name) +{ + if (name.substr(0, 3) != "ctr" || name[3] != '_') + { + throw_invalid_ctr_token(name); + } + size_t pos = name.find("_", 4); + if (pos == std::string::npos) + { + throw_invalid_ctr_token(name); + } + int start_index; + int end_index; + try + { + start_index = atoi(name.substr(4, pos - 1).c_str()); + end_index = atoi(name.substr(pos + 1, std::string::npos).c_str()); + } + catch(const std::exception& e) + { + throw_invalid_ctr_token(name); + } + return {start_index, end_index}; +} + +void throw_invalid_ctr_token(const std::string& name) +{ + std::cerr << "ERROR: \"" << name << "\" is not a valid token format. The available cycle tranisition rate (ctr) tokens are \"ctr__\"." + << " For example: \"ctr_0_1\" or \"ctr_2_3\"." << std::endl; + exit(-1); +} + +MappingFunction select_phenotype_by_token_inputter(const std::string& name, const std::string& sbml_species) +{ + if (name[0] == 'm') + { + if (name == "mms") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.motility.migration_speed); }; + } + else if (name == "mpt") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.motility.persistence_time); }; + } + else if (name == "mmb") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.motility.migration_bias); }; + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized motility token. The available motility tokens are \"mms\", \"mpt\", and \"mmb\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + } + else if (name[0] == 'd') + { + if (name == "da") + { + int death_model_index = PhysiCell::cell_defaults.phenotype.death.find_death_model_index(PhysiCell::PhysiCell_constants::apoptosis_death_model); + return [sbml_species, death_model_index](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.death.rates[death_model_index]); }; + } + else if (name == "dn") + { + int death_model_index = PhysiCell::cell_defaults.phenotype.death.find_death_model_index(PhysiCell::PhysiCell_constants::necrosis_death_model); + return [sbml_species, death_model_index](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.death.rates[death_model_index]); }; + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized death token. The available death tokens are \"da\" and \"dn\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + } + else if (name[0] == 's') + { + size_t pos = name.find("_"); + std::string substrate_name = name.substr(pos + 1, std::string::npos); + int substrate_index = microenvironment.find_density_index(substrate_name); + + if (substrate_index == -1) + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << substrate_name << "\" is not a recognized substrate in this model." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + + std::string token_prefix = name.substr(0, 3); + + if (token_prefix == "sur") + { + return [sbml_species, substrate_index](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.secretion.uptake_rates[substrate_index]); }; + } + else if (token_prefix == "ssr") + { + return [sbml_species, substrate_index](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.secretion.secretion_rates[substrate_index]); }; + } + else if (token_prefix == "ssd") + { + return [sbml_species, substrate_index](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.secretion.saturation_densities[substrate_index]); }; + } + else if (token_prefix == "ser") + { + return [sbml_species, substrate_index](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.secretion.net_export_rates[substrate_index]); }; + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized secretion token. The available secretion tokens are \"sur\", \"ssr\", \"ssd\", and \"ser\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + } + else if (name[0] == 'c') + { + std::vector indices = parse_ctr_token(name); + int start_index = indices[0]; + int end_index = indices[1]; + + return [sbml_species, start_index, end_index](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.cycle.data.transition_rate(start_index, end_index)); }; + } + else if (name[0] == 'v') + { + if (name == "vtsc") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.volume.target_solid_cytoplasmic); }; + } + else if (name == "vtsn") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.volume.target_solid_nuclear); }; + } + else if (name == "vff") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.intracellular->set_parameter_value(sbml_species, pCell->phenotype.volume.target_fluid_fraction); }; + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized volume token. The available volume tokens are \"vtsc\", \"vtsn\", and \"vff\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized token. It must start with \"m\", \"d\", \"s\", \"c\", or \"v\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } +} + +MappingFunction select_phenotype_by_token_outputter(const std::string& name, const std::string& sbml_species) +{ + if (name[0] == 'm') + { + if (name == "mms") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.motility.migration_speed = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (name == "mpt") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.motility.persistence_time = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (name == "mmb") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.motility.migration_bias = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized motility token. The available motility tokens are \"mms\", \"mpt\", and \"mmb\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + } + else if (name[0] == 'd') + { + if (name == "da") + { + int death_model_index = PhysiCell::cell_defaults.phenotype.death.find_death_model_index(PhysiCell::PhysiCell_constants::apoptosis_death_model); + return [sbml_species, death_model_index](PhysiCell::Cell *pCell) + { pCell->phenotype.death.rates[death_model_index] = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (name == "dn") + { + int death_model_index = PhysiCell::cell_defaults.phenotype.death.find_death_model_index(PhysiCell::PhysiCell_constants::necrosis_death_model); + return [sbml_species, death_model_index](PhysiCell::Cell *pCell) + { pCell->phenotype.death.rates[death_model_index] = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized death token. The available death tokens are \"da\" and \"dn\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + } + else if (name[0] == 's') + { + size_t pos = name.find("_"); + std::string substrate_name = name.substr(pos + 1, std::string::npos); + int substrate_index = microenvironment.find_density_index(substrate_name); + + if (substrate_index == -1) + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << substrate_name << "\" is not a substrate name in this model." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + + std::string token_prefix = name.substr(0, 3); + + if (token_prefix == "sur") + { + return [sbml_species, substrate_index](PhysiCell::Cell *pCell) + { pCell->phenotype.secretion.uptake_rates[substrate_index] = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (token_prefix == "ssr") + { + return [sbml_species, substrate_index](PhysiCell::Cell *pCell) + { pCell->phenotype.secretion.secretion_rates[substrate_index] = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (token_prefix == "ssd") + { + return [sbml_species, substrate_index](PhysiCell::Cell *pCell) + { pCell->phenotype.secretion.saturation_densities[substrate_index] = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (token_prefix == "ser") + { + return [sbml_species, substrate_index](PhysiCell::Cell *pCell) + { pCell->phenotype.secretion.net_export_rates[substrate_index] = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized secretion token. The available secretion tokens are \"sur\", \"ssr\", \"ssd\", and \"ser\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + } + else if (name[0] == 'c') + { + std::vector indices = parse_ctr_token(name); + int start_index = indices[0]; + int end_index = indices[1]; + + return [sbml_species, start_index, end_index](PhysiCell::Cell *pCell) + { pCell->phenotype.cycle.data.transition_rate(start_index, end_index) = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (name[0] == 'v') + { + if (name == "vtsc") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.volume.target_solid_cytoplasmic = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (name == "vtsn") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.volume.target_solid_nuclear = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else if (name == "vff") + { + return [sbml_species](PhysiCell::Cell *pCell) + { pCell->phenotype.volume.target_fluid_fraction = pCell->phenotype.intracellular->get_parameter_value(sbml_species); }; + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized volume token. The available volume tokens are \"vtsc\", \"vtsn\", and \"vff\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } + } + else + { + std::cerr<< std::endl; + std::cerr << "ERROR: \"" << name << "\" is not a recognized token. It must start with \"m\", \"d\", \"s\", \"c\", or \"v\"." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); + } +} + 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; - + input_mappings = copy->input_mappings; + output_mappings = copy->output_mappings; + mappings_initialized = copy->mappings_initialized; } -// 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) { pugi::xml_node node_sbml = node.child( "sbml_filename" ); @@ -61,234 +533,189 @@ void RoadRunnerIntracellular::initialize_intracellular_from_pugixml(pugi::xml_no sbml_filename = PhysiCell::xml_get_my_string_value (node_sbml); std::cout << "\n------------- " << __FUNCTION__ << ": sbml_filename = " << sbml_filename << std::endl; } - - pugi::xml_node node_species = node.child( "map" ); - while( node_species ) - { - // --------- substrate - - 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; - } - // --------- custom_data - std::string custom_data_name = node_species.attribute( "PC_custom_data" ).value(); - if( custom_data_name != "" ) - { - 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; - } - - - // --------- phenotype_data - std::string phenotype_name = node_species.attribute( "PC_phenotype" ).value(); - - if( phenotype_name != "" ) - { - 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" ); - } - - std::cout << " ------- substrate_species map:" << std::endl; - for(auto elm : substrate_species) - { - std::cout << " " << elm.first << " -> " << elm.second << std::endl; - } - std::cout << " ------- custom_data_species map:" << std::endl; - for(auto elm : custom_data_species) - { - std::cout << " " << elm.first << " -> " << elm.second << std::endl; + + update_time_step = PhysiCell::intracellular_dt; // default to this, but overwrite below if defined in XML + 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; } - std::cout << std::endl; - std::cout << " ------- phenotype_species map:" << std::endl; - for(auto elm : phenotype_species) + std::vector new_input_mappings; + std::vector new_output_mappings; + + pugi::xml_node node_map = node.child( "map" ); + while ( node_map ) { - std::cout << " " << elm.first << " -> " << elm.second << std::endl; + std::string io_type = node_map.attribute( "type" ).value(); + if (io_type != "io" && io_type != "input" && io_type != "output") + { + std::cout << "\n------------- " << __FUNCTION__ << ": ERROR: type must be io, input, or output" << std::endl; + exit(-1); + } + + std::string physicell_name = node_map.attribute( "physicell_name" ).value(); + std::string sbml_species = node_map.attribute( "sbml_species" ).value(); + + if (io_type == "input" || io_type == "io") + { + new_input_mappings.push_back(new RoadRunnerMapping(physicell_name, sbml_species, "input")); + } + if (io_type == "output" || io_type == "io") + { + new_output_mappings.push_back(new RoadRunnerMapping(physicell_name, sbml_species, "output")); + } + + node_map = node_map.next_sibling("map"); } - std::cout << std::endl; + validate_mappings(new_input_mappings, true); + validate_mappings(new_output_mappings, false); + input_mappings = std::move(new_input_mappings); + output_mappings = std::move(new_output_mappings); } +void validate_mappings(std::vector mappings, bool is_inputs) +{ + if (mappings.empty()) + { return; } + std::vector values_already_set; + std::string name_to_add; + for (auto mapping : mappings) + { + name_to_add = is_inputs ? mapping->sbml_species : mapping->physicell_name; + if (std::find(values_already_set.begin(), values_already_set.end(), name_to_add) != values_already_set.end()) + { + std::cout << "ERROR: the " << (is_inputs ? "SBML species " : "PhysiCell name ") << name_to_add << " is set by multiple " << (is_inputs ? "inputs" : "outputs") << std::endl; + exit(-1); + } + values_already_set.push_back(name_to_add); + } +} 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; + + initialize_mappings(); } 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() +void RoadRunnerIntracellular::update(PhysiCell::Cell* pCell, PhysiCell::Phenotype& phenotype, double dt) { + pre_update(pCell); + // update the intracellular model 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; + static int num_vals = 2; // start time and end time - // 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; - - // 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->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; + post_update(pCell); +} - // 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" +void RoadRunnerIntracellular::initialize_mappings() +{ + for (auto mapping : input_mappings) + { + if (mapping->mapping_initialized) + { continue; } + mapping->initialize_mapping(); + } + for (auto mapping : output_mappings) + { + if (mapping->mapping_initialized) + { continue; } + mapping->initialize_mapping(); + } + mappings_initialized = true; +} - // 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++) +void RoadRunnerIntracellular::pre_update(PhysiCell::Cell* pCell) +{ + for (auto mapping : input_mappings) { - // 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 << " "; - // } + mapping->value_map(pCell); } - //std::cout << "\n"; +} - index = 0; - // Print out the data - for (int row = 0; row < this->result->RSize; row++) +void RoadRunnerIntracellular::post_update(PhysiCell::Cell* pCell) +{ + for (auto mapping : output_mappings) { - for (int col = 0; col < this->result->CSize; col++) + mapping->value_map(pCell); + } +} + +RoadRunnerMapping* RoadRunnerIntracellular::find_input_mapping(std::string sbml_species) +{ + for (auto mapping : input_mappings) + { + if (mapping->sbml_species == sbml_species) { - // 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 << " "; - // } + return mapping; } - // 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 nullptr; +} - // return 0; +RoadRunnerMapping* RoadRunnerIntracellular::find_output_mapping(std::string physicell_name) +{ + for (auto mapping : output_mappings) + { + if (mapping->physicell_name == physicell_name) + { + return mapping; + } + } + return nullptr; } 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,20 +728,21 @@ 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) { return static_cast(phenotype.intracellular); } +RoadRunnerIntracellular* getRoadRunnerModel(PhysiCell::Cell* pCell) { + return getRoadRunnerModel(pCell->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; @@ -328,349 +756,41 @@ std::string RoadRunnerIntracellular::get_state() return sbml_filename; } - -int RoadRunnerIntracellular::update_phenotype_parameters(PhysiCell::Phenotype& phenotype) +int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phenotype) { - for(auto elm : phenotype_species) - { - // motility params - if (elm.first[0] == 'm') - { - if (elm.first == "mms") - { - phenotype.motility.migration_speed = phenotype.intracellular->get_parameter_value(elm.second); - } - else if (elm.first == "mpt") - { - phenotype.motility.persistence_time = phenotype.intracellular->get_parameter_value(elm.second); - } - else if (elm.first == "mmb") - { - phenotype.motility.migration_bias = phenotype.intracellular->get_parameter_value(elm.second); - } - else - { - } - } - // death params - else if (elm.first[0] == 'd') - { - if (elm.first == "da") - { - phenotype.death.rates[0] = phenotype.intracellular->get_parameter_value(elm.second); - } - else if (elm.first == "dn") - { - phenotype.death.rates[1] = phenotype.intracellular->get_parameter_value(elm.second); - } - else - { - } - } - // secretion params - else if (elm.first[0] == 's') - { - // parsing attribute and getting substrate name - std::string s = elm.first; - std::string delimiter = "_"; - - size_t pos = 0; - std::string token; - while ((pos = s.find(delimiter)) != std::string::npos) { - token = s.substr(0, pos); - s.erase(0, pos + delimiter.length()); - } - int sub_index = microenvironment.find_density_index(s); - - //transport types - //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") - { - phenotype.secretion.secretion_rates[sub_index] = phenotype.intracellular->get_parameter_value(elm.second); - } - //secretion density - else if (elm.first.substr(0,3) == "ssd") - { - phenotype.secretion.saturation_densities[sub_index] = phenotype.intracellular->get_parameter_value(elm.second); - } - //net export rate - else if (elm.first.substr(0,3) == "ser") - { - phenotype.secretion.net_export_rates[sub_index] = phenotype.intracellular->get_parameter_value(elm.second); - } - else - { - } - } - - // cycle params - else if (elm.first[0] == 'c') - { - if (elm.first.substr(0,3) == "ctr") - { - // parsing attribute and getting substrate name - std::string s = elm.first; - std::string delimiter = "_"; - - size_t pos = 0; - std::string token; - int counter = 0; - 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() ); - } - s.erase(0, pos + delimiter.length()); - 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 - { - } - } - - // volume params - else if (elm.first[0] == 'v') - { - if (elm.first == "vtsc") - { - phenotype.volume.target_solid_cytoplasmic = phenotype.intracellular->get_parameter_value(elm.second); - } - else if (elm.first == "vtsn") - { - phenotype.volume.target_solid_nuclear = phenotype.intracellular->get_parameter_value(elm.second); - } - else if (elm.first == "vff") - { - phenotype.volume.target_fluid_fraction = phenotype.intracellular->get_parameter_value(elm.second); - } - else - { - } - } - else - { - } - - } - //std::cout << std::endl; + // the mappings are not yet intialized by this time because the signals/behaviors dictionaries are not yet initialized + int num_of_phases = (&(phenotype.cycle.model()))->phases.size(); + validate_cycle_mappings(input_mappings, num_of_phases); + validate_cycle_mappings(output_mappings, num_of_phases); return 0; } - -int RoadRunnerIntracellular::validate_PhysiCell_tokens(PhysiCell::Phenotype& phenotype) +void validate_cycle_mappings(std::vector mappings, int num_of_phases) { - for(auto elm : phenotype_species) + for (auto mapping : mappings) { - //std::cout << "PhysiCell_token_validation" << std::endl; - //std::cout << elm.first << " : " << elm.second << std::endl; - - // motility params - if (elm.first[0] == 'm') - { - if (elm.first == "mms") - { - } - else if (elm.first == "mpt") - { - } - else if (elm.first == "mmb") - { - } - else - { - std::cout<< std::endl; - 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; - } - } - // death params - else if (elm.first[0] == 'd') - { - if (elm.first == "da") - { - } - else if (elm.first == "dn") - { - } - else - { - std::cout<< std::endl; - 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; - } - } - // secretion params - else if (elm.first[0] == 's') - { - // parsing attribute and getting substrate name - std::string s = elm.first; - std::string delimiter = "_"; - size_t pos = 0; - std::string token; - while ((pos = s.find(delimiter)) != std::string::npos) { - token = s.substr(0, pos); - 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; - } - - if (elm.first.substr(0,3) == "sur") - { - } - else if (elm.first.substr(0,3) == "ssr") - { - } - else if (elm.first.substr(0,3) == "ssd") - { - } - else if (elm.first.substr(0,3) == "ser") - { - } - else - { - std::cout<< std::endl; - 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; - } - } - else if (elm.first[0] == 'c') - { - if (elm.first.substr(0,3) == "ctr") - { - // 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; - std::string delimiter = "_"; - size_t pos = 0; - std::string token; - int counter = 0; - int start_index; - while ((pos = s.find(delimiter)) != std::string::npos) { - token = s.substr(0, pos); - if (counter == 1) - { - start_index = atoi( token.c_str() ); - } - s.erase(0, pos + delimiter.length()); - counter += 1; - } - int end_index = atoi( s.c_str() ); - - // validating the indices - if ( start_index > num_of_phases - 1) - { - std::cout<< std::endl; - 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; - } - if ( end_index > num_of_phases - 1) - { - std::cout<< std::endl; - 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; - } - } - else - { - std::cout<< std::endl; - 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; - } - } - - else if (elm.first[0] == 'v') - { - if (elm.first == "vtsc") - { - } - else if (elm.first == "vtsn") - { - } - else if (elm.first == "vff") - { - } - else - { - std::cout<< std::endl; - 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; - } - } - else + // check that the mapping uses a cycle transition rate, i.e. starts with "ctr_" + if (mapping->physicell_name.find("ctr_") != 0) + { continue; } + + std::vector indices = parse_ctr_token(mapping->physicell_name); + if (indices[0] > num_of_phases - 1 || indices[1] > num_of_phases - 1) { - std::cout<< std::endl; - 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; + std::cerr << "ERROR: The token \"" << mapping->physicell_name << "\" is invalid for this cell type. The indices are out of range." << std::endl + << " The phase indices must be between 0 and " << num_of_phases - 1 << "." << std::endl; + exit(-1); } - } - std::cout << "---- Specified PhysiCell tokens at config file are validated. ----- " << std::endl; - - return 0; } 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,105 +804,32 @@ 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++; } - // Phenotype 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..... " < all_species, std::vector mappings) +{ + for (auto mapping : mappings) { - bool exist = 0; - // std::cout << species_name.size() << std::endl; - for (int i=0; i < all_species.size(); i++) + if (std::find(all_species.begin(), all_species.end(), mapping->sbml_species) == all_species.end()) { - //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..... " <sbml_species << "\" at " << mapping->io_type << " mapping. Please take a look SBML species specifications." << std::endl; + std::cerr<< std::endl; + std::cerr<< std::endl; + exit(-1); } - //std::cout << "existence check : " << elm.second <<": " << exist << std::endl; - } - - - //std::cout << "---------VALIDATING_SBML_SPECIES END-------" << std::endl; - - - std::cout << "---- Specified SBML species at config file are validated. ----- " << std::endl; - return 0; + } } int RoadRunnerIntracellular::create_custom_data_for_SBML(PhysiCell::Phenotype& phenotype) { - //std::cout << "Test" << std::endl; - return 0; } diff --git a/addons/libRoadrunner/src/librr_intracellular.h b/addons/libRoadrunner/src/librr_intracellular.h index 4ca8b041a..855a5c266 100644 --- a/addons/libRoadrunner/src/librr_intracellular.h +++ b/addons/libRoadrunner/src/librr_intracellular.h @@ -18,42 +18,66 @@ // #include "rrc_types.h" #include "../roadrunner/include/rr/C/rrc_api.h" #include "../roadrunner/include/rr/C/rrc_types.h" +#include +#include + // #include "rrc_utilities.h" extern "C" rrc::RRHandle createRRInstance(); // #endif +typedef std::function MappingFunction; + +class RoadRunnerMapping +{ +public: + std::string physicell_name; + std::string sbml_species; + std::string io_type; + std::string physicell_dictionary_name; + int index; + MappingFunction value_map = [] (PhysiCell::Cell *pCell) {}; // default to a function that does nothing + bool mapping_initialized = false; + + RoadRunnerMapping() {}; + RoadRunnerMapping(std::string physicell_name, std::string sbml_species, std::string io_type) + : physicell_name(physicell_name), sbml_species(sbml_species), io_type(io_type) {}; + + void initialize_mapping(void); +}; + +MappingFunction select_signal_setter(const std::string& name, const std::string& sbml_species); + +bool is_physicell_phenotype_token(const std::string& name); +MappingFunction select_phenotype_by_token_inputter(const std::string& name, const std::string& sbml_species); +MappingFunction select_phenotype_by_token_outputter(const std::string& name, const std::string& sbml_species); + +void validate_mappings(std::vector mappings, bool is_inputs); + +std::vector parse_ctr_token(const std::string &name); +void throw_invalid_ctr_token(const std::string& name); +void validate_cycle_mappings(std::vector mappings, int num_of_phases); + class RoadRunnerIntracellular : public PhysiCell::Intracellular { private: public: - - // static long counter; - std::string sbml_filename; - // bool enabled = false; - int num_rows_result_table = 1; - - // double time_step = 12; - // bool discrete_time = false; - // double time_tick = 0.5; - // double scaling = 1.0; - - // std::map initial_values; + std::map parameters; - std::map substrate_species; - std::map custom_data_species; - std::map phenotype_species; + bool mappings_initialized = false; + void initialize_mappings(); + std::vector input_mappings; + std::vector output_mappings; 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,15 +85,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; return static_cast(clone); } @@ -78,26 +96,35 @@ class RoadRunnerIntracellular : public PhysiCell::Intracellular std::cout << "------ librr_intracellular: getIntracellularModel called\n"; return static_cast(this); } - - void initialize_intracellular_from_pugixml(pugi::xml_node& node); - - // Need 'int' return type to avoid bizarre compile errors? But 'void' to match MaBoSS. + + void initialize_intracellular_from_pugixml(pugi::xml_node &node); + + // Need 'int' return type to avoid bizarre compile errors? But 'void' to match MaBoSS. void start(); bool need_update(); - // Need 'int' return type to avoid bizarre compile errors. - void update(); - void update(PhysiCell::Cell* cell, PhysiCell::Phenotype& phenotype, double dt) { - update(); - update_phenotype_parameters(phenotype); - } + void update() {}; // needed because the base class has this function + void update(PhysiCell::Cell* cell, PhysiCell::Phenotype& phenotype, double dt); + + void pre_update(PhysiCell::Cell* cell); + void post_update(PhysiCell::Cell* cell); void inherit(PhysiCell::Cell * cell) {} - - int update_phenotype_parameters(PhysiCell::Phenotype& phenotype); + + // These find__mapping functions are not currently used, but since I made them, we'll keep them around. + RoadRunnerMapping *find_input_mapping(std::string sbml_species); // sbml_species is unique for inputs (below is for convenience) + RoadRunnerMapping *find_input_mapping(std::string physicell_name, std::string sbml_species) + { return find_input_mapping(sbml_species); } // sbml_species is unique for inputs + + RoadRunnerMapping *find_output_mapping(std::string physicell_name); // physicell_name is unique for outputs (below is for convenience) + RoadRunnerMapping *find_output_mapping(std::string physicell_name, std::string sbml_species) + { return find_output_mapping(physicell_name); } // physicell_name is unique for outputs + + int update_phenotype_parameters(PhysiCell::Phenotype& phenotype) {return 0;}; // all handled within update int validate_PhysiCell_tokens(PhysiCell::Phenotype& phenotype); int validate_SBML_species(); + void validate_SBML_species(std::vector all_species, std::vector mappings); int create_custom_data_for_SBML(PhysiCell::Phenotype& phenotype); double get_parameter_value(std::string name); @@ -114,4 +141,6 @@ class RoadRunnerIntracellular : public PhysiCell::Intracellular static void save_libRR(std::string path, std::string index); }; +RoadRunnerIntracellular* getRoadRunnerModel(PhysiCell::Phenotype& phenotype); +RoadRunnerIntracellular* getRoadRunnerModel(PhysiCell::Cell* pCell); #endif \ No newline at end of file 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..09bd79e26 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -149,16 +149,20 @@ 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_); + } } } } - + if( time_since_last_cycle > phenotype_dt_ - 0.5 * diffusion_dt_ || !initialzed ) { // Reset the max_radius in each voxel. It will be filled in set_total_volume diff --git a/core/PhysiCell_signal_behavior.cpp b/core/PhysiCell_signal_behavior.cpp index 8780e8d7c..daad14bfa 100644 --- a/core/PhysiCell_signal_behavior.cpp +++ b/core/PhysiCell_signal_behavior.cpp @@ -725,10 +725,8 @@ int find_signal_index( std::string signal_name ) auto search = signal_to_int.find( signal_name ); // safety first! if( search != signal_to_int.end() ) - { return search->second; } - - std::cout << "having trouble finding " << signal_name << std::endl; - + { return search->second; } + return -1; } 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..589c27b19 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,23 +308,27 @@ ./config/Toy_Metabolic_Model.xml - 0.01 - - - - - - - + 0.1 + + + + + + + + + + + - 0.0 - 0.0 + 0.8 + 15.0 0.0 - 0.0 + 450.0 @@ -347,10 +352,6 @@ 0 - 0.8 - 15 - 0 - 450 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..9ca79abec 100644 --- a/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.cpp +++ b/sample_projects_intracellular/ode/ode_energy/custom_modules/custom.cpp @@ -69,12 +69,6 @@ #include "../BioFVM/BioFVM.h" using namespace BioFVM; - -#include "rrc_api.h" -#include "rrc_types.h" -// #include "rrc_utilities.h" -extern "C" rrc::RRHandle createRRInstance(); - void create_cell_types( void ) { // set the random seed @@ -133,38 +127,16 @@ void setup_microenvironment( void ) void setup_tissue( void ) { - 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"); - - double Xmin = microenvironment.mesh.bounding_box[0]; - double Ymin = microenvironment.mesh.bounding_box[1]; - double Zmin = microenvironment.mesh.bounding_box[2]; - - double Xmax = microenvironment.mesh.bounding_box[3]; - double Ymax = microenvironment.mesh.bounding_box[4]; - double Zmax = microenvironment.mesh.bounding_box[5]; - - if( default_microenvironment_options.simulate_2D == true ) - { - Zmin = 0.0; - Zmax = 0.0; - } - - double Xrange = Xmax - Xmin; - double Yrange = Ymax - Ymin; - double Zrange = Zmax - Zmin; - // create cells Cell* pCell; double cell_radius = cell_defaults.phenotype.geometry.radius; - double cell_spacing = 0.8 * 2.0 * cell_radius; double initial_tumor_radius = 100; - double retval; std::vector> positions = create_cell_circle_positions(cell_radius,initial_tumor_radius); @@ -174,86 +146,17 @@ void setup_tissue( void ) pCell = create_cell(get_cell_definition("default")); pCell->assign_position( positions[i] ); - set_single_behavior( pCell , "custom:intra_oxy" , parameters.doubles("initial_internal_oxygen")); - set_single_behavior( pCell , "custom:intra_glu" , parameters.doubles("initial_internal_glucose")); - 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")); } return; } -void update_intracellular() -{ - // 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") ); - - } - } - -} - - - std::vector my_coloring_function( Cell* pCell ) { @@ -288,8 +191,6 @@ std::vector my_coloring_function( Cell* pCell ) return output; } - - std::vector> create_cell_circle_positions(double cell_radius, double sphere_radius) { std::vector> cells; 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..9a3c6f552 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,6 @@ void setup_microenvironment( void ); // custom pathology coloring function std::vector my_coloring_function( Cell* ); -void update_intracellular(); // 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; }