diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index 1f3e2646b..5f2376f21 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -363,8 +363,10 @@ void Cell::advance_bundled_phenotype_functions( double dt_ ) // check for new death events if( phenotype.death.check_for_death( dt_ ) == true ) { - // if so, change the cycle model to the current death model - phenotype.cycle.sync_to_cycle_model( phenotype.death.current_model() ); + // if so, change the cycle model to the current death model, taking the + // parameters from this cell definition rather than from the shared model (#199) + phenotype.cycle.sync_to_cycle_model( phenotype.death.current_model() , + phenotype.death.current_model_data() ); // also, turn off motility. @@ -500,8 +502,10 @@ void Cell::start_death( int death_model_index ) { // set the death data struture to the indicated death model phenotype.death.trigger_death( death_model_index ); - // change the cycle model to the current death model - phenotype.cycle.sync_to_cycle_model( phenotype.death.current_model() ); + // change the cycle model to the current death model, taking the parameters + // from this cell definition rather than from the shared model (#199) + phenotype.cycle.sync_to_cycle_model( phenotype.death.current_model() , + phenotype.death.current_model_data() ); // turn off secretion, and reduce uptake by a factor of 10 phenotype.secretion.set_all_secretion_to_zero(); @@ -1820,7 +1824,9 @@ void display_cell_definitions( std::ostream& os ) << " with rate " << pCD->phenotype.death.rates[k] << " 1/min" << std::endl; Cycle_Model* pCM = (pCD->phenotype.death.models[k] ); - Cycle_Data* pCMD = &(pCD->phenotype.death.models[k]->data ); + // phases and links come from the shared model, but the parameters are + // per-definition now, so take those from model_data (#199) + Cycle_Data* pCMD = &(pCD->phenotype.death.model_data[k] ); os << "\t\tdeath phase transitions: " << std::endl @@ -2252,8 +2258,9 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node // set the transition rate pCD->phenotype.cycle.data.transition_rate(start,end) = value; - // set it to fixed / non-fixed - pCD->phenotype.cycle.model().phase_link(start,end).fixed_duration = fixed; + // set it to fixed / non-fixed -- into this definition's own data, not + // the shared Cycle_Model (#199) + pCD->phenotype.cycle.data.fixed_duration(start,end) = fixed; node = node.next_sibling( "rate" ); } @@ -2281,7 +2288,7 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node // set the transition rate pCD->phenotype.cycle.data.exit_rate(start) = 1.0 / (value+1e-16); // set it to fixed / non-fixed - pCD->phenotype.cycle.model().phase_links[start][0].fixed_duration = fixed; + pCD->phenotype.cycle.data.exit_fixed_duration(start) = fixed; node = node.next_sibling( "duration" ); } @@ -2510,9 +2517,10 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node double value = xml_get_my_double_value( node1 ); // set the transition rate - pCD->phenotype.death.models[death_index]->transition_rate(start,end) = value; - // set it to fixed / non-fixed - pCD->phenotype.death.models[death_index]->phase_link(start,end).fixed_duration = fixed; + pCD->phenotype.death.model_data[death_index].transition_rate(start,end) = value; + // set it to fixed / non-fixed -- into this definition's own copy, not + // the shared Cycle_Model (#199) + pCD->phenotype.death.model_data[death_index].fixed_duration(start,end) = fixed; node1 = node1.next_sibling( "rate" ); } @@ -2534,10 +2542,11 @@ Cell_Definition* initialize_cell_definition_from_pugixml( pugi::xml_node cd_node double value = xml_get_my_double_value( node ); // set the transition rate - pCD->phenotype.death.models[death_index]->data.exit_rate(start) + pCD->phenotype.death.model_data[death_index].exit_rate(start) = 1.0 / (value+1e-16); - // set it to fixed / non-fixed - pCD->phenotype.death.models[death_index]->phase_links[start][0].fixed_duration + // set it to fixed / non-fixed -- into this definition's own copy, not + // the shared Cycle_Model (#199) + pCD->phenotype.death.model_data[death_index].exit_fixed_duration(start) = fixed; node = node.next_sibling( "duration" ); diff --git a/core/PhysiCell_phenotype.cpp b/core/PhysiCell_phenotype.cpp index e725c0170..8beb7240f 100644 --- a/core/PhysiCell_phenotype.cpp +++ b/core/PhysiCell_phenotype.cpp @@ -94,7 +94,6 @@ Phase_Link::Phase_Link() start_phase_index = 0; end_phase_index = 0; - fixed_duration = false; arrest_function = NULL; exit_function = NULL; @@ -126,16 +125,18 @@ void Cycle_Data::sync_to_cycle_model( void ) // querying the phase_links transition_rates.resize( n ); + fixed_durations.resize( n ); // also make sure the transition_rates[] are the right size for( int i=0 ; i < pCycle_Model->phase_links.size() ; i++ ) { inverse_index_maps[i].clear(); + transition_rates[i].resize( pCycle_Model->phase_links[i].size() ); + fixed_durations[i].resize( pCycle_Model->phase_links[i].size() , (char) 0 ); for( int j=0 ; j < pCycle_Model->phase_links[i].size() ; j++ ) { inverse_index_maps[i][ pCycle_Model->phase_links[i][j].end_phase_index ] = j; - transition_rates[i].resize( pCycle_Model->phase_links[i].size() ); } } @@ -151,6 +152,16 @@ double& Cycle_Data::exit_rate(int phase_index ) { return transition_rates[phase_index][0]; } + +char& Cycle_Data::fixed_duration( int start_phase_index , int end_phase_index ) +{ + return fixed_durations[ start_phase_index ][ inverse_index_maps[start_phase_index][end_phase_index] ]; +} + +char& Cycle_Data::exit_fixed_duration( int phase_index ) +{ + return fixed_durations[phase_index][0]; +} Cycle_Model::Cycle_Model() { @@ -301,7 +312,7 @@ void Cycle_Model::advance_model( Cell* pCell, Phenotype& phenotype, double dt ) { // check to see if we should transition bool continue_transition = false; - if( phase_links[i][k].fixed_duration ) + if( phenotype.cycle.data.fixed_durations[i][k] ) { if( phenotype.cycle.data.elapsed_time_in_phase > ((1.0/phenotype.cycle.data.transition_rates[i][k]) - 0.5 * dt) ) { @@ -398,6 +409,7 @@ int Death::add_death_model( double rate , Cycle_Model* pModel ) { rates.push_back( rate ); models.push_back( pModel ); + model_data.push_back( pModel->data ); // seed per-definition params from the model parameters.resize( rates.size() ); @@ -408,6 +420,7 @@ int Death::add_death_model( double rate, Cycle_Model* pModel, Death_Parameters& { rates.push_back( rate ); models.push_back( pModel ); + model_data.push_back( pModel->data ); // seed per-definition params from the model parameters.push_back( death_parameters ); return rates.size() - 1; @@ -496,6 +509,11 @@ Cycle_Model& Death::current_model( void ) return *models[current_death_model_index]; } +Cycle_Data& Death::current_model_data( void ) +{ + return model_data[ current_death_model_index ]; +} + double& Death::apoptosis_rate(void) { static int nApoptosis = find_death_model_index( PhysiCell_constants::apoptosis_death_model ); @@ -537,9 +555,16 @@ int& Cycle::current_phase_index( void ) } void Cycle::sync_to_cycle_model( Cycle_Model& cm ) +{ + sync_to_cycle_model( cm , cm.data ); + return; +} + +void Cycle::sync_to_cycle_model( Cycle_Model& cm , const Cycle_Data& cd ) { pCycle_Model = &cm; - data = cm.data; + data = cd; + data.pCycle_Model = &cm; // cd may have been seeded from a different instance return; } diff --git a/core/PhysiCell_phenotype.h b/core/PhysiCell_phenotype.h index 13edc0c40..7c9ac21d8 100644 --- a/core/PhysiCell_phenotype.h +++ b/core/PhysiCell_phenotype.h @@ -123,8 +123,6 @@ class Phase_Link int start_phase_index; int end_phase_index; - bool fixed_duration; - bool (*arrest_function)( Cell* pCell, Phenotype& phenotype, double dt ); // return true if arrested, false if not @@ -152,6 +150,13 @@ class Cycle_Data std::vector< std::vector > transition_rates; + // Whether each transition is a fixed duration rather than a Poisson rate. + // Indexed exactly like transition_rates. char rather than bool so that + // fixed_duration() can hand back a reference (std::vector is a proxy + // container and cannot). Lives here, not in Phase_Link, because Cycle_Model + // objects are shared by pointer between cell definitions -- see #199. + std::vector< std::vector > fixed_durations; + int current_phase_index; double elapsed_time_in_phase; @@ -169,6 +174,10 @@ class Cycle_Data double& exit_rate(int phase_index ); // This returns the first transition rate out of // phase # phase_index. It is only relevant if the phase has only one phase link // (true for many cycle models). + + // as transition_rate / exit_rate, for the fixed-duration flag + char& fixed_duration(int start_phase_index, int end_phase_index ); // done + char& exit_fixed_duration(int phase_index ); // done }; class Cycle_Model @@ -247,6 +256,10 @@ class Cycle int& current_phase_index( void ); // done void sync_to_cycle_model( Cycle_Model& cm ); // done + // as above, but take the parameters from cd rather than from cm.data. Needed + // because death Cycle_Models are shared between cell definitions, so the + // per-definition parameters live elsewhere -- see Death::model_data (#199). + void sync_to_cycle_model( Cycle_Model& cm , const Cycle_Data& cd ); // done Asymmetric_Division asymmetric_division; }; @@ -275,6 +288,11 @@ class Death public: std::vector rates; std::vector models; + // Per-definition cycle parameters for each death model. models[] are shared + // Cycle_Model objects, so their own .data cannot hold per-cell-definition + // durations or fixed-duration flags -- see #199. Seeded from the model at + // add_death_model(), overridden by the XML, applied at start_death(). + std::vector model_data; std::vector parameters; bool dead; @@ -292,6 +310,7 @@ class Death void trigger_death( int death_model_index ); // done Cycle_Model& current_model( void ); // done + Cycle_Data& current_model_data( void ); // done Death_Parameters& current_parameters( void ); // done ' // ease of access diff --git a/core/PhysiCell_standard_models.cpp b/core/PhysiCell_standard_models.cpp index 4d170a406..91ef9cff7 100644 --- a/core/PhysiCell_standard_models.cpp +++ b/core/PhysiCell_standard_models.cpp @@ -239,7 +239,7 @@ void create_ki67_models( void ) Ki67_basic.transition_rate(0,1) = 1.0/(4.59*60.0); // MCF10A cells are ~4.59 hours in Ki67- state Ki67_basic.transition_rate(1,0) = 1.0/(15.5*60.0); // length of Ki67+ states in advanced model - Ki67_basic.phase_link(1,0).fixed_duration = true; + Ki67_basic.data.fixed_duration(1,0) = true; Ki67_basic.phases[0].entry_function = NULL; // standard_Ki67_negative_phase_entry_function; Ki67_basic.phases[1].entry_function = standard_Ki67_positive_phase_entry_function; @@ -261,8 +261,8 @@ void create_ki67_models( void ) Ki67_advanced.add_phase_link( 1 , 2 , NULL ); // + (pre-mitotic) to + (post-mitotic) Ki67_advanced.add_phase_link( 2 , 0 , NULL ); // + to - - Ki67_advanced.phase_link(1,2).fixed_duration = true; - Ki67_advanced.phase_link(2,0).fixed_duration = true; + Ki67_advanced.data.fixed_duration(1,2) = true; + Ki67_advanced.data.fixed_duration(2,0) = true; Ki67_advanced.transition_rate(0,1) = 1.0/(3.62*60.0); // MCF10A cells ~3.62 hours in Ki67- in this fitted model Ki67_advanced.transition_rate(1,2) = 1.0/(13.0*60.0); @@ -376,7 +376,7 @@ void create_cycling_quiescent_model( void ) cycling_quiescent.transition_rate(0,1) = 1.0/(4.59*60.0); // MCF10A cells are ~4.59 hours in Ki67- state cycling_quiescent.transition_rate(1,0) = 1.0/(15.5*60.0); // length of Ki67+ states in advanced model - cycling_quiescent.phase_link(1,0).fixed_duration = true; + cycling_quiescent.data.fixed_duration(1,0) = true; cycling_quiescent.phases[0].entry_function = NULL; cycling_quiescent.phases[1].entry_function = standard_cycling_entry_function; @@ -440,7 +440,7 @@ void create_standard_apoptosis_model( void ) apoptosis.transition_rate( 0, 1) = 1.0 / (8.6 * 60.0); // Use the deterministic model, where this phase has fixed duration - apoptosis.phase_link(0,1).fixed_duration = true; + apoptosis.data.fixed_duration(0,1) = true; return; } @@ -482,7 +482,7 @@ void create_standard_necrosis_model( void ) necrosis.transition_rate( 1, 2 ) = 1.0 / (60.0 * 24.0 * 60.0 ); // 60 days max // Deterministically remove the necrotic cell if it has been 60 days - necrosis.phase_link(1,2).fixed_duration = true; + necrosis.data.fixed_duration(1,2) = true; return; }