Merge pull request #17 from MathCancer/development

Development

BioFVM/BioFVM_basic_agent.cpp

@@ -74,7 +74,15 @@ Basic_Agent::Basic_Agent()
 	saturation_densities= new std::vector<double>(0);
 	// extern Microenvironment* default_microenvironment;
 	// register_microenvironment( default_microenvironment ); 
+
+	internalized_substrates = new std::vector<double>(0); // 
+	fraction_released_at_death = new std::vector<double>(0); 
+	fraction_transferred_when_ingested = new std::vector<double>(0); 
 	register_microenvironment( get_default_microenvironment() );
+
+	// these are done in register_microenvironment
+	// internalized_substrates.assign( get_default_microenvironment()->number_of_densities() , 0.0 ); 
+
 	return;	
 }
 
@@ -88,7 +96,6 @@ bool Basic_Agent::assign_position(std::vector<double> new_position)
 
 bool Basic_Agent::assign_position(double x, double y, double z)
 {
-	// std::cout << __FILE__ << " " << __LINE__ << std::endl; 
 	if( !get_microenvironment()->mesh.is_position_valid(x,y,z))
 	{	
 		// std::cout<<"Error: the new position for agent "<< ID << " is invalid: "<<x<<","<<y<<","<<"z"<<std::endl;
@@ -122,20 +129,34 @@ void Basic_Agent::set_internal_uptake_constants( double dt )
 	//   p(n+1)*temp2 =  p(n) + temp1
 	//   p(n+1) = (  p(n) + temp1 )/temp2
 	//int nearest_voxel= current_voxel_index;
+
+/*	
+	// new for tracking internal densities
+	if( use_internal_densities_as_targets == true )
+	{
+		*saturation_densities = *internalized_substrates;
+		*saturation_densities /= ( 1e-15 + volume ); 
+	}
+*/
+
 	double internal_constant_to_discretize_the_delta_approximation = dt * volume / ( (microenvironment->voxels(current_voxel_index)).volume ) ; // needs a fix 
-
+	
 	// temp1 = dt*(V_cell/V_voxel)*S*T 
 	cell_source_sink_solver_temp1.assign( (*secretion_rates).size() , 0.0 ); 
 	cell_source_sink_solver_temp1 += *secretion_rates; 
 	cell_source_sink_solver_temp1 *= *saturation_densities; 
 	cell_source_sink_solver_temp1 *= internal_constant_to_discretize_the_delta_approximation; 
+
+//	total_extracellular_substrate_change.assign( (*secretion_rates).size() , 1.0 ); 
 
 	// temp2 = 1 + dt*(V_cell/V_voxel)*( S + U )
 	cell_source_sink_solver_temp2.assign( (*secretion_rates).size() , 1.0 ); 
 	axpy( &(cell_source_sink_solver_temp2) , internal_constant_to_discretize_the_delta_approximation , *secretion_rates );
 	axpy( &(cell_source_sink_solver_temp2) , internal_constant_to_discretize_the_delta_approximation , *uptake_rates );	
 
 	volume_is_changed = false; 
+
+	return; 
 }
 
 void Basic_Agent::register_microenvironment( Microenvironment* microenvironment_in )
@@ -148,6 +169,14 @@ void Basic_Agent::register_microenvironment( Microenvironment* microenvironment_
 	// some solver temporary variables 
 	cell_source_sink_solver_temp1.resize( microenvironment->density_vector(0).size() , 0.0 );
 	cell_source_sink_solver_temp2.resize( microenvironment->density_vector(0).size() , 1.0 );
+
+	// new for internalized substrate tracking 
+	internalized_substrates->resize( microenvironment->density_vector(0).size() , 0.0 );
+	total_extracellular_substrate_change.resize( microenvironment->density_vector(0).size() , 1.0 );
+
+	fraction_released_at_death->resize( microenvironment->density_vector(0).size() , 0.0 ); 
+	fraction_transferred_when_ingested->resize( microenvironment->density_vector(0).size() , 0.0 ); 
+
 	return; 
 }
 
@@ -156,7 +185,22 @@ Microenvironment* Basic_Agent::get_microenvironment( void )
 
 Basic_Agent::~Basic_Agent()
 {
- return; 
+	Microenvironment* pS = get_default_microenvironment(); 
+
+	// change in total in voxel: 
+	// total_ext = total_ext + fraction*total_internal 
+	// total_ext / vol_voxel = total_ext / vol_voxel + fraction*total_internal / vol_voxel 
+	// density_ext += fraction * total_internal / vol_volume 
+
+	// std::cout << "\t\t\t" << (*pS)(current_voxel_index) << "\t\t\t" << std::endl; 
+	*internalized_substrates /=  pS->voxels(current_voxel_index).volume; // turn to density 
+	*internalized_substrates *= *fraction_released_at_death;  // what fraction is released? 
+
+	// release this amount into the environment 
+
+	(*pS)(current_voxel_index) += *internalized_substrates; 
+
+	return; 
 }
 
 Basic_Agent* create_basic_agent( void )
@@ -237,6 +281,7 @@ double Basic_Agent::get_total_volume()
 {
 	return volume;
 }
+
 void Basic_Agent::simulate_secretion_and_uptake( Microenvironment* pS, double dt )
 {
 	if(!is_active)
@@ -247,6 +292,20 @@ void Basic_Agent::simulate_secretion_and_uptake( Microenvironment* pS, double dt
 		set_internal_uptake_constants(dt);
 		volume_is_changed = false;
 	}
+
+	if( default_microenvironment_options.track_internalized_substrates_in_each_agent == true )
+	{
+		total_extracellular_substrate_change.assign( total_extracellular_substrate_change.size() , 1.0 ); // 1
+
+		total_extracellular_substrate_change -= cell_source_sink_solver_temp2; // 1-c2
+		total_extracellular_substrate_change *= (*pS)(current_voxel_index); // (1-c2)*rho 
+		total_extracellular_substrate_change += cell_source_sink_solver_temp1; // (1-c2)*rho+c1 
+		total_extracellular_substrate_change /= cell_source_sink_solver_temp2; // ((1-c2)*rho+c1)/c2
+		total_extracellular_substrate_change *= pS->voxels(current_voxel_index).volume; // W*((1-c2)*rho+c1)/c2 
+
+		*internalized_substrates -= total_extracellular_substrate_change; // opposite of net extracellular change 	
+	}
+
 	(*pS)(current_voxel_index) += cell_source_sink_solver_temp1; 
 	(*pS)(current_voxel_index) /= cell_source_sink_solver_temp2; 
 

BioFVM/BioFVM_basic_agent.h

@@ -74,14 +74,21 @@ class Basic_Agent
 	std::vector<double> previous_velocity; 
 	bool is_active;
 
+	std::vector<double> total_extracellular_substrate_change; 
+
  public:
 	std::vector<double> * secretion_rates; 
 	std::vector<double> * saturation_densities; 
 	std::vector<double> * uptake_rates;  
 	double get_total_volume();
 	void set_total_volume(double);
 	void update_voxel_index();
-
+
+	/* new for internalized substrates in 1.5.0 */ 
+	std::vector<double> * internalized_substrates; 
+	std::vector<double> * fraction_released_at_death; 
+	std::vector<double> * fraction_transferred_when_ingested; 
+
 	void set_internal_uptake_constants( double dt ); // any time you update the cell volume or rates, should call this function. 
 
 	void register_microenvironment( Microenvironment* );

BioFVM/BioFVM_microenvironment.cpp

@@ -217,6 +217,13 @@ void Microenvironment::update_dirichlet_node( int voxel_index , std::vector<doub
 	return; 
 }
 
+void Microenvironment::update_dirichlet_node( int voxel_index , int substrate_index , double new_value )
+{
+	mesh.voxels[voxel_index].is_Dirichlet = true; 
+	dirichlet_value_vectors[voxel_index][substrate_index] = new_value; 
+	return; 
+}
+
 void Microenvironment::remove_dirichlet_node( int voxel_index )
 {
 	mesh.voxels[voxel_index].is_Dirichlet = false; 
@@ -431,6 +438,8 @@ void Microenvironment::resize_densities( int new_size )
 	default_microenvironment_options.Dirichlet_condition_vector.assign( new_size , 1.0 );  
 	default_microenvironment_options.Dirichlet_activation_vector.assign( new_size, true ); 
 
+	default_microenvironment_options.initial_condition_vector.assign( new_size , 1.0 ); 
+
 	return; 
 }
 
@@ -483,6 +492,8 @@ void Microenvironment::add_density( void )
 	default_microenvironment_options.Dirichlet_condition_vector.push_back( 1.0 ); //  = one; 
 	default_microenvironment_options.Dirichlet_activation_vector.push_back( true ); // assign( number_of_densities(), true ); 
 
+	default_microenvironment_options.initial_condition_vector.push_back( 1.0 ); 
+
 	return; 
 }
 
@@ -534,6 +545,8 @@ void Microenvironment::add_density( std::string name , std::string units )
 	default_microenvironment_options.Dirichlet_condition_vector.push_back( 1.0 ); //  = one; 
 	default_microenvironment_options.Dirichlet_activation_vector.push_back( true ); // assign( number_of_densities(), true ); 
 
+	default_microenvironment_options.initial_condition_vector.push_back( 1.0 ); 
+
 	return; 
 }
 
@@ -585,6 +598,8 @@ void Microenvironment::add_density( std::string name , std::string units, double
 	default_microenvironment_options.Dirichlet_condition_vector.push_back( 1.0 ); // = one; 
 	default_microenvironment_options.Dirichlet_activation_vector.push_back( true ); // assign( number_of_densities(), true ); 
 
+	default_microenvironment_options.initial_condition_vector.push_back( 1.0 ); 
+
 	return; 
 }
 
@@ -1043,6 +1058,8 @@ Microenvironment_Options::Microenvironment_Options()
 	Dirichlet_condition_vector.assign( pMicroenvironment->number_of_densities() , 1.0 ); 
 	Dirichlet_activation_vector.assign( pMicroenvironment->number_of_densities() , true ); 
 
+	initial_condition_vector.resize(0); //  = Dirichlet_condition_vector; 
+
 	// set a far-field value for oxygen (assumed to be in the first field)
 	Dirichlet_condition_vector[0] = 38.0; 
 
@@ -1059,20 +1076,15 @@ Microenvironment_Options::Microenvironment_Options()
 
 	calculate_gradients = false; 
 
+	track_internalized_substrates_in_each_agent = false; 
+
 	return; 
 }
 
 Microenvironment_Options default_microenvironment_options; 
 
 void initialize_microenvironment( void )
 {
-	/*
-	if( default_microenvironment_options.use_oxygen_as_first_field == false )
-	{
-		std::cout << "wha???" << std::endl; 
-	}
-	*/
-
 	// create and name a microenvironment; 
 	microenvironment.name = default_microenvironment_options.name;
 	// register the diffusion solver 
@@ -1109,9 +1121,21 @@ void initialize_microenvironment( void )
 	microenvironment.time_units = default_microenvironment_options.time_units;
 	microenvironment.mesh.units = default_microenvironment_options.spatial_units;
 
-	// set the initial densities to the values set in the Dirichlet_condition_vector
+	// set the initial densities to the values set in the initial_condition_vector
+
+	// if the initial condition vector has not been set, use the Dirichlet condition vector 
+	if( default_microenvironment_options.initial_condition_vector.size() != 
+		microenvironment.number_of_densities() )
+	{
+		std::cout << "BioFVM Warning: Initial conditions not set. " << std::endl 
+				  << "                Using Dirichlet condition vector to set initial substrate values!" << std::endl 
+				  << "                In the future, set default_microenvironment_options.initial_condition_vector." 
+				  << std::endl << std::endl;  
+		default_microenvironment_options.initial_condition_vector = default_microenvironment_options.Dirichlet_condition_vector; 
+	}
+
 	for( unsigned int n=0; n < microenvironment.number_of_voxels() ; n++ )
-	{ microenvironment.density_vector(n) = default_microenvironment_options.Dirichlet_condition_vector; }
+	{ microenvironment.density_vector(n) = default_microenvironment_options.initial_condition_vector; }
 
 	if( default_microenvironment_options.outer_Dirichlet_conditions == true ) 
 	{

BioFVM/BioFVM_microenvironment.h

@@ -234,6 +234,7 @@ class Microenvironment
 
 	void add_dirichlet_node( int voxel_index, std::vector<double>& value ); 
 	void update_dirichlet_node( int voxel_index , std::vector<double>& new_value ); 
+	void update_dirichlet_node( int voxel_index , int substrate_index , double new_value );
 	void remove_dirichlet_node( int voxel_index ); 
 	void apply_dirichlet_conditions( void ); 
 
@@ -300,6 +301,8 @@ class Microenvironment_Options
 	std::vector<double> Dirichlet_condition_vector; 
 	std::vector<bool> Dirichlet_activation_vector; 
 
+	std::vector<double> initial_condition_vector; 
+
 	bool simulate_2D; 
 	std::vector<double> X_range; 
 	std::vector<double> Y_range; 
@@ -310,6 +313,8 @@ class Microenvironment_Options
 	bool calculate_gradients; 
 
 	bool use_oxygen_as_first_field;
+
+	bool track_internalized_substrates_in_each_agent; 	
 };
 
 extern Microenvironment_Options default_microenvironment_options; 

BioFVM/BioFVM_vector.cpp

@@ -172,11 +172,13 @@ void operator/=( std::vector<double>& v1, const double& a )
 
 std::ostream& operator<<(std::ostream& os, const std::vector<double>& v )
 {
+ /*
  if( v.size() == 3 )
  {
   os << "x=\"" << v[0] << "\" y=\"" << v[1] << "\" z=\"" << v[2] << "\"" ; 
   return os; 
  }
+ */
 
  for( unsigned int i=0; i < v.size(); i++ )
  { os << v[i] << " " ; }

CITATION.txt

@@ -1,7 +1,7 @@
 If you use PhysiCell in your project, please cite PhysiCell and the version 
 number, such as below:                                                      
 
-We implemented and solved the model using PhysiCell (Version 1.4.1) [1].    
+We implemented and solved the model using PhysiCell (Version 1.5.0) [1].    
 
 [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, 
     PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu-  
@@ -11,7 +11,7 @@ We implemented and solved the model using PhysiCell (Version 1.4.1) [1].
 Because PhysiCell extensively uses BioFVM, we suggest you also cite BioFVM  
     as below:                                                               
 
-We implemented and solved the model using PhysiCell (Version 1.4.1) [1],    
+We implemented and solved the model using PhysiCell (Version 1.5.0) [1],    
 with BioFVM [2] to solve the transport equations.                           
 
 [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, 

Makefile

@@ -118,6 +118,15 @@ cancer-immune-sample:
 	cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml 
 	cp ./sample_projects/cancer_immune/config/* ./config/
 
+virus-macrophage-sample:
+	cp ./sample_projects/virus_macrophage/custom_modules/* ./custom_modules/
+	touch main.cpp && cp main.cpp main-backup.cpp
+	cp ./sample_projects/virus_macrophage/main-virus_macrophage.cpp ./main.cpp 
+	cp Makefile Makefile-backup
+	cp ./sample_projects/virus_macrophage/Makefile .
+	cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml 
+	cp ./sample_projects/virus_macrophage/config/* ./config/
+
 beta-testing:
 	cp ./sample_projects/beta_testing/custom_modules/* ./custom_modules/
 	touch main.cpp && cp main.cpp main-backup.cpp