Merge branch 'my-physicell-dev' into my-pc-dev-delete-duplicate-elapsed_time_in_phase

Author: drbergman

.gitignore

@@ -19,5 +19,8 @@ project
 studio_debug.log
 user_projects/*
 !user_projects/empty.txt
+
+project*
 Studio.zip
 /studio
+.vscode

BioFVM/BioFVM_basic_agent.cpp

@@ -46,6 +46,9 @@
 #############################################################################
 */
 
+// #include <mutex>
+// std::mutex mtx;
+
 #include "BioFVM_basic_agent.h"
 #include "BioFVM_agent_container.h"
 #include "BioFVM_vector.h" 
@@ -176,6 +179,36 @@ void Basic_Agent::set_internal_uptake_constants( double dt )
 	return; 
 }
 
+void Basic_Agent::set_transmembrane_diffusion_constants( double dt )
+{
+	// overall form: dI/dt = s*(E-I); dE/dt = s*(I-E) + e
+	// where I is the internal substrate concentration, E is the extracellular substrate concentration (in voxel), s is the diffusion (secretion) rate, and e is the net export rate
+	// use analytical solution to update concentrations (splitting up the diffusion and export components because that's what I solved for)
+	double voxel_coeff = 1 / (volume + (microenvironment->voxels(current_voxel_index)).volume); // 1 / (cell volume + voxel volume)
+	double cell_coeff = (microenvironment->voxels(current_voxel_index)).volume * voxel_coeff; // (voxel volume / (cell volume + voxel volume))
+	double lambda2_base_dt = -dt*(1+volume/(microenvironment->voxels(current_voxel_index)).volume); // -dt * (1 + V_cell/V_voxel)
+
+	double exp_decay_term;
+	for (unsigned int i = 0; i < (*secretion_rates).size(); i++)
+	{
+		exp_decay_term = (1 - exp((*secretion_rates)[i] * lambda2_base_dt));
+		cell_source_sink_solver_temp1[i] = cell_coeff * exp_decay_term;
+		cell_source_sink_solver_temp2[i] = voxel_coeff * exp_decay_term;
+	}
+	
+	// temp for net export 
+	cell_source_sink_solver_temp_export1 = *net_export_rates; 
+	cell_source_sink_solver_temp_export1 *= dt; // amount exported in dt of time 
+		
+	// change in surrounding density 
+	cell_source_sink_solver_temp_export2 = cell_source_sink_solver_temp_export1;
+	cell_source_sink_solver_temp_export2 /= ( (microenvironment->voxels(current_voxel_index)).volume ) ; 
+	
+	volume_is_changed = false; 
+	
+	return; 
+}
+
 void Basic_Agent::register_microenvironment( Microenvironment* microenvironment_in )
 {
 	microenvironment = microenvironment_in; 	
@@ -334,9 +367,9 @@ void Basic_Agent::simulate_secretion_and_uptake( Microenvironment* pS, double dt
 		*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; 
-	
+	(*pS)(current_voxel_index) += cell_source_sink_solver_temp1;
+	(*pS)(current_voxel_index) /= cell_source_sink_solver_temp2;
+
 	// now do net export 
 	(*pS)(current_voxel_index) += cell_source_sink_solver_temp_export2; 
 	if( default_microenvironment_options.track_internalized_substrates_in_each_agent == true ) 
@@ -347,4 +380,44 @@ void Basic_Agent::simulate_secretion_and_uptake( Microenvironment* pS, double dt
 	return; 
 }
 
+void Basic_Agent::simulate_transmembrane_diffusion( Microenvironment* pS, double dt )
+{
+	if(!is_active)
+	{ return; }
+	
+	if( volume_is_changed )
+	{
+		set_transmembrane_diffusion_constants(dt);
+		volume_is_changed = false;
+	}
+	
+	// double conc_diff;
+	double stuff_diff;
+	std::vector<double>& v = nearest_density_vector();
+	for (unsigned int i = 0; i < (*secretion_rates).size(); i++)
+	{
+		// conc_diff = (*internalized_substrates)[i]/volume - nearest_density_vector()[i];
+		stuff_diff = (*internalized_substrates)[i] - volume * v[i]; // this is the concentration difference times the cell volume; done to avoid division by cell volume, which could perhaps be (close to) zero
+		if( default_microenvironment_options.track_internalized_substrates_in_each_agent == true )
+		{
+			// (*internalized_substrates)[i] -= cell_source_sink_solver_temp1[i] * conc_diff;
+			(*internalized_substrates)[i] -= cell_source_sink_solver_temp1[i] * stuff_diff;
+		}
+		// (*pS)(current_voxel_index)[i] += cell_source_sink_solver_temp2[i] * conc_diff;
+		// mtx.lock(); // I think there is the danger of a data race...but I'm not sure how to fix it (or if it's even a problem) (I'm not sure if this would handle it anyway because when v is defined above, it could be changed before it's accessed in defining stuff_diff)
+		v[i] += cell_source_sink_solver_temp2[i] * stuff_diff;
+		v[i] += cell_source_sink_solver_temp_export2[i];
+		// mtx.unlock();
+	}
+
+	// now do net export 
+	// (*pS)(current_voxel_index) += cell_source_sink_solver_temp_export2; 
+	if( default_microenvironment_options.track_internalized_substrates_in_each_agent == true ) 
+	{
+		*internalized_substrates -= cell_source_sink_solver_temp_export1; 
+	}
+
+	return; 
+}
+
 };

BioFVM/BioFVM_basic_agent.h

@@ -88,6 +88,10 @@ class Basic_Agent
 	std::vector<double> * net_export_rates; 
 	double get_total_volume();
 	void set_total_volume(double);
+	void set_volume_is_changed( bool new_value ){
+		volume_is_changed = new_value; 
+		return; 
+	}
 	void update_voxel_index();
 
 	/* new for internalized substrates in 1.5.0 */ 
@@ -97,6 +101,7 @@ class Basic_Agent
 	void release_internalized_substrates( void ); 
 
 	void set_internal_uptake_constants( double dt ); // any time you update the cell volume or rates, should call this function. 
+	void set_transmembrane_diffusion_constants( double dt ); // any time you update the cell volume or rates, should call this function. 
 
 	void register_microenvironment( Microenvironment* );
 	Microenvironment* get_microenvironment( void ); 
@@ -117,6 +122,7 @@ class Basic_Agent
 	// simulate secretion and uptake at the nearest voxel at the indicated microenvironment.
 	// if no microenvironment indicated, use the currently selected microenvironment. 
 	void simulate_secretion_and_uptake( Microenvironment* M, double dt ); 
+	void simulate_transmembrane_diffusion( Microenvironment* M, double dt ); 
 
 	int get_current_voxel_index( void ); 
 	// directly access the substrate vector at the nearest voxel at the indicated microenvironment