From 48eae7056dd0429fc7763f640477e11e5ec7d210 Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Thu, 6 Aug 2026 20:25:40 -0400 Subject: [PATCH 1/3] Give each cell definition its own fixed_duration and death-phase parameters Fixes #199. Cycle_Model objects are shared by pointer between cell definitions: apoptosis and necrosis are globals, and every definition registers the same address via add_death_model( rate, &apoptosis, ... ). Cell_Definition's copy constructor copies the pointer, not the model. fixed_duration lived in Phase_Link, inside that shared object, while the XML parser wrote it per cell definition. So each definition overwrote the flag for every other one and the last definition parsed decided for all of them. A model asking for deterministic apoptosis on most of its cell types silently got the stochastic branch on all of them, with no warning and nothing in the saved output to reveal it. Migrate the flag from the graph structure to the param structure, as Paul suggested on the issue: Cycle_Data gains fixed_durations, indexed exactly like transition_rates, with fixed_duration() / exit_fixed_duration() accessors mirroring transition_rate() / exit_rate(). Cycle_Model::advance_model now reads the per-cell copy. Phase_Link::fixed_duration is removed; the six standard models declare their defaults through the model's own Cycle_Data instead, and those defaults still reach cells that do not override them in XML. That alone is not enough for the death models. There is no per-definition Cycle_Data for a death model -- Death holds rates, models and parameters, none of which carry cycle parameters -- so the parser wrote durations to the shared models[i]->data, and Cycle::sync_to_cycle_model() copies that same shared data over the cell at start_death(). Relocating the flag would have moved it from one shared object to another. This is also why rheiland observed on the issue that the value is flattened along with the flag. So Death gains model_data, one Cycle_Data per death model, seeded from the model at add_death_model(), written by the parser, and applied at start_death() through a new two-argument Cycle::sync_to_cycle_model( cm, cd ). The single-argument form delegates to it and keeps its exact previous behaviour. Verified on a 7-cell-type model where six definitions ask for and the last asks for a rate. Before: all seven report fixed=0. After: six report 1 and the last reports 0, while all seven still share one Cycle_Model object, so the phase graph stays shared and Phenotype copies stay cheap. A definition with no apoptosis at all still inherits the standard model's default. Note: this changes results for any model that mixes and across its cell definitions. Those models were not getting the death timing their XML asked for; they will now. Co-Authored-By: Claude Opus 5 --- core/PhysiCell_cell.cpp | 27 +++++++++++++++----------- core/PhysiCell_phenotype.cpp | 31 +++++++++++++++++++++++++++--- core/PhysiCell_phenotype.h | 23 ++++++++++++++++++++-- core/PhysiCell_standard_models.cpp | 12 ++++++------ 4 files changed, 71 insertions(+), 22 deletions(-) diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index 1f3e2646b..dc499dfd4 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -500,8 +500,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(); @@ -2252,8 +2254,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 +2284,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 +2513,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 +2538,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..f9b38ae5a 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,6 +125,7 @@ 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 @@ -136,6 +136,7 @@ void Cycle_Data::sync_to_cycle_model( void ) { inverse_index_maps[i][ pCycle_Model->phase_links[i][j].end_phase_index ] = j; transition_rates[i].resize( pCycle_Model->phase_links[i].size() ); + fixed_durations[i].resize( pCycle_Model->phase_links[i].size() , (char) 0 ); } } @@ -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 , 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..e20ed6624 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 , 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; } From 720652bfd0a074c71cab381e634e64d1bbc171ae Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Thu, 6 Aug 2026 20:48:44 -0400 Subject: [PATCH 2/3] Apply the per-definition death parameters on the runtime death path too Cell::advance_bundled_phenotype_functions syncs the cycle to the death model directly when check_for_death fires, rather than going through start_death. That site was still using the single-argument sync_to_cycle_model, so cells dying by death_rate -- the common case -- kept getting parameters from the shared Cycle_Model and #199 was unfixed on that path. Also from review: hoist the transition_rates / fixed_durations resize out of the inner phase-link loop in Cycle_Data::sync_to_cycle_model, and take the Cycle_Data by const reference in the new Cycle overload. Co-Authored-By: Claude Opus 5 --- core/PhysiCell_cell.cpp | 6 ++++-- core/PhysiCell_phenotype.cpp | 6 +++--- core/PhysiCell_phenotype.h | 2 +- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index dc499dfd4..5cd87fb8d 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. diff --git a/core/PhysiCell_phenotype.cpp b/core/PhysiCell_phenotype.cpp index f9b38ae5a..8beb7240f 100644 --- a/core/PhysiCell_phenotype.cpp +++ b/core/PhysiCell_phenotype.cpp @@ -132,11 +132,11 @@ void Cycle_Data::sync_to_cycle_model( void ) 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() ); - fixed_durations[i].resize( pCycle_Model->phase_links[i].size() , (char) 0 ); } } @@ -560,7 +560,7 @@ void Cycle::sync_to_cycle_model( Cycle_Model& cm ) return; } -void Cycle::sync_to_cycle_model( Cycle_Model& cm , Cycle_Data& cd ) +void Cycle::sync_to_cycle_model( Cycle_Model& cm , const Cycle_Data& cd ) { pCycle_Model = &cm; data = cd; diff --git a/core/PhysiCell_phenotype.h b/core/PhysiCell_phenotype.h index e20ed6624..7c9ac21d8 100644 --- a/core/PhysiCell_phenotype.h +++ b/core/PhysiCell_phenotype.h @@ -259,7 +259,7 @@ class Cycle // 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 , Cycle_Data& cd ); // done + void sync_to_cycle_model( Cycle_Model& cm , const Cycle_Data& cd ); // done Asymmetric_Division asymmetric_division; }; From edc522493b285aea2818dd852b022f2182601a72 Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Fri, 7 Aug 2026 07:41:47 -0400 Subject: [PATCH 3/3] Report per-definition death parameters in the cell definition summary display_cell_definitions() read the death phase durations out of death.models[k]->data, the shared Cycle_Model. The parser no longer writes that object, so after the rest of this PR the summary printed the compiled-in standard-model defaults for every definition regardless of the XML -- and that summary is exactly what a user reads to check whether #199 is fixed. Take the parameters from death.model_data[k] instead. Phases and phase links still come from the shared model, which is correct: only the parameters moved. Co-Authored-By: Claude Opus 5 --- core/PhysiCell_cell.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index 5cd87fb8d..5f2376f21 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -1824,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