Give each cell definition its own fixed_duration and death-phase parameters

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" ); 

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; 
 }	
 

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<double> > 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<bool> 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<char> > 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<double> rates; 
 	std::vector<Cycle_Model*> 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<Cycle_Data> model_data; 
 	std::vector<Death_Parameters> 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

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; 
 }