Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
261 changes: 136 additions & 125 deletions BioFVM/BioFVM_MultiCellDS.cpp

Large diffs are not rendered by default.

22 changes: 12 additions & 10 deletions BioFVM/BioFVM_MultiCellDS.h
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,8 @@
#ifndef __BioFVM_MultiCellDS_h__
#define __BioFVM_MultiCellDS_h__

#include "BioFVM_microenvironment_interface.h"
#include "BioFVM_basic_agent_interface.h"
#include "pugixml.hpp"

#include <cstdio>
Expand Down Expand Up @@ -160,7 +162,7 @@ class MultiCellDS_Metadata

MultiCellDS_Metadata();
void display_information( std::ostream& os);
void sync_to_microenvironment( Microenvironment& M );
void sync_to_microenvironment( Microenvironment_Interface& M );
void restart_runtime( void );

void add_to_open_xml_pugi( double current_simulation_time, pugi::xml_document& xml_dom );
Expand All @@ -182,29 +184,29 @@ void set_save_biofvm_cell_data_as_custom_matlab( bool newvalue ); // default: tr
/* writing parts of BioFVM to a MultiCellDS file */

void reset_BioFVM_substrates_initialized_in_dom( void );
void add_BioFVM_substrates_to_open_xml_pugi( pugi::xml_document& xml_dom , std::string filename_base , Microenvironment& M );
void add_BioFVM_basic_agent_to_open_xml_pugi( pugi::xml_document& xml_dom, Basic_Agent& BA ); // not implemented -- future edition
void add_BioFVM_agents_to_open_xml_pugi( pugi::xml_document& xml_dom, std::string filename_base, Microenvironment& M );
void add_BioFVM_to_open_xml_pugi( pugi::xml_document& xml_dom , std::string filename_base, double current_simulation_time , Microenvironment& M );
void add_BioFVM_substrates_to_open_xml_pugi( pugi::xml_document& xml_dom , std::string filename_base , Microenvironment_Interface& M );
void add_BioFVM_basic_agent_to_open_xml_pugi( pugi::xml_document& xml_dom, Basic_Agent_Interface& BA ); // not implemented -- future edition
void add_BioFVM_agents_to_open_xml_pugi( pugi::xml_document& xml_dom, std::string filename_base, Microenvironment_Interface& M );
void add_BioFVM_to_open_xml_pugi( pugi::xml_document& xml_dom , std::string filename_base, double current_simulation_time , Microenvironment_Interface& M );

void write_coordinates_node(pugi::xml_node &node, const std::vector<double> &coordinates, std::string name);

void save_BioFVM_to_MultiCellDS_xml_pugi( std::string filename_base , Microenvironment& M , double current_simulation_time);
void save_BioFVM_to_MultiCellDS_xml_pugi( std::string filename_base , Microenvironment_Interface& M , double current_simulation_time);

/* beta in PhysiCell 1.11.0 */

bool read_microenvironment_from_matlab( std::string mat_filename );

/* future / not yet supported */

void read_BioFVM_from_open_xml_pugi( pugi::xml_document& xml_dom , std::string filename_base, double& current_simulation_time , Microenvironment& M );
void read_BioFVM_to_MultiCellDS_xml_pugi( std::string filename_base , Microenvironment& M , double& current_simulation_time );
void read_BioFVM_from_open_xml_pugi( pugi::xml_document& xml_dom , std::string filename_base, double& current_simulation_time , Microenvironment_Interface& M );
void read_BioFVM_to_MultiCellDS_xml_pugi( std::string filename_base , Microenvironment_Interface& M , double& current_simulation_time );

/* partly-implemented code snippets -- not to be used as of March 2016 */

// functions to read multiscale_microenvironment from MultiCellDS file (requires pugixml)
void read_microenvironment_from_MultiCellDS_xml( Microenvironment& M_destination , std::string filename );
void read_microenvironment_from_MultiCellDS_xml( Microenvironment& M_destination , pugi::xml_document& xml_dom );
void read_microenvironment_from_MultiCellDS_xml( Microenvironment_Interface& M_destination , std::string filename );
void read_microenvironment_from_MultiCellDS_xml( Microenvironment_Interface& M_destination , pugi::xml_document& xml_dom );

};

Expand Down
1 change: 0 additions & 1 deletion BioFVM/BioFVM_agent_container.h
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,6 @@
namespace BioFVM{

class Basic_Agent;
class Microenvironment;

class Agent_Container
{
Expand Down
237 changes: 195 additions & 42 deletions BioFVM/BioFVM_basic_agent.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -75,16 +75,16 @@ Basic_Agent::Basic_Agent()
velocity.assign( 3 , 0.0 );
previous_velocity.assign( 3 , 0.0 );
// link into the microenvironment, if one is defined
secretion_rates= new std::vector<double>(0);
uptake_rates= new std::vector<double>(0);
saturation_densities= new std::vector<double>(0);
net_export_rates = new std::vector<double>(0);
secretion_rates= std::vector<double>(0);
uptake_rates= std::vector<double>(0);
saturation_densities= std::vector<double>(0);
net_export_rates = 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>(1.0);
internalized_substrates = std::vector<double>(0); //
fraction_released_at_death = std::vector<double>(0);
fraction_transferred_when_ingested = std::vector<double>(1.0);
register_microenvironment( get_default_microenvironment() );

// these are done in register_microenvironment
Expand Down Expand Up @@ -151,20 +151,20 @@ void Basic_Agent::set_internal_uptake_constants( double dt )
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.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 );
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 );

// temp for net export
cell_source_sink_solver_temp_export1 = *net_export_rates;
cell_source_sink_solver_temp_export1 = net_export_rates;
cell_source_sink_solver_temp_export1 *= dt; // amount exported in dt of time

cell_source_sink_solver_temp_export2 = cell_source_sink_solver_temp_export1;
Expand All @@ -176,27 +176,37 @@ void Basic_Agent::set_internal_uptake_constants( double dt )
return;
}

void Basic_Agent::register_microenvironment( Microenvironment_Interface* microenvironment_in )
{
auto me = dynamic_cast<BioFVM::Microenvironment*>(microenvironment_in);
if (!me) {
throw std::invalid_argument("Basic_Agent::register_microenvironment: Provided Microenvironment_Interface is not a BioFVM::Microenvironment.");
}
register_microenvironment(me);
}

void Basic_Agent::register_microenvironment( Microenvironment* microenvironment_in )
{
microenvironment = microenvironment_in;
secretion_rates->resize( microenvironment->density_vector(0).size() , 0.0 );
saturation_densities->resize( microenvironment->density_vector(0).size() , 0.0 );
uptake_rates->resize( microenvironment->density_vector(0).size() , 0.0 );
net_export_rates->resize( microenvironment->density_vector(0).size() , 0.0 );
unsigned int num_densities = microenvironment->number_of_densities();
secretion_rates.resize( num_densities , 0.0 );
saturation_densities.resize( num_densities , 0.0 );
uptake_rates.resize( num_densities , 0.0 );
net_export_rates.resize( num_densities , 0.0 );

// 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 );
cell_source_sink_solver_temp1.resize( num_densities , 0.0 );
cell_source_sink_solver_temp2.resize( num_densities , 1.0 );

cell_source_sink_solver_temp_export1.resize( microenvironment->density_vector(0).size() , 0.0 );
cell_source_sink_solver_temp_export2.resize( microenvironment->density_vector(0).size() , 0.0 );
cell_source_sink_solver_temp_export1.resize( num_densities , 0.0 );
cell_source_sink_solver_temp_export2.resize( num_densities , 0.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 );
internalized_substrates.resize( num_densities , 0.0 );
total_extracellular_substrate_change.resize( num_densities , 1.0 );

fraction_released_at_death->resize( microenvironment->density_vector(0).size() , 0.0 );
fraction_transferred_when_ingested->resize( microenvironment->density_vector(0).size() , 1.0 );
fraction_released_at_death.resize( num_densities , 0.0 );
fraction_transferred_when_ingested.resize( num_densities , 1.0 );

return;
}
Expand All @@ -211,16 +221,16 @@ void Basic_Agent::release_internalized_substrates( void )
// 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?
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;
(*pS)(current_voxel_index) += internalized_substrates;

// zero out the now-removed substrates

internalized_substrates->assign( internalized_substrates->size() , 0.0 );
internalized_substrates.assign( internalized_substrates.size() , 0.0 );

return;
}
Expand Down Expand Up @@ -277,9 +287,9 @@ int Basic_Agent::get_current_voxel_index( void )
return current_voxel_index;
}

std::vector<double>& Basic_Agent::nearest_density_vector( void )
double* Basic_Agent::nearest_density_vector( void )
{
return microenvironment->nearest_density_vector( current_voxel_index );
return (*microenvironment)( current_voxel_index ).data();
}


Expand All @@ -301,16 +311,159 @@ void Basic_Agent::set_total_volume(double volume)
volume_is_changed = true;
}

double Basic_Agent::get_total_volume()
double& Basic_Agent::get_total_volume()
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a specialized setter for volume that sets the volume_is_changed flag. Are we sure we need to access this by reference?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

void Basic_Agent::set_total_volume(double volume)
{
	this->volume = volume;
	volume_is_changed = true;
}

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah so this is necessary because now that Basic_Agent::get_total_volume() is virtual, it is clashing with Cell::get_total_volume(), which returns the reference.

So we either unify Cell's and Basic_Agents's get_total_volume to return double& or we can actually remove Cell::get_total_volume (since its deprecated) and keep double.

I am open to ideas

{
return volume;
}

const std::vector<double>& Basic_Agent::get_previous_velocity( void ) {
// Implementation of interface methods

double* Basic_Agent::get_position_internal()
{
return position.data();
}

const std::vector<double>& Basic_Agent::get_position() const
{
return position;
}

std::vector<double>& Basic_Agent::get_velocity()
{
return velocity;
}

const std::vector<double>& Basic_Agent::get_velocity() const
{
return velocity;
}

std::vector<double>& Basic_Agent::get_previous_velocity( void )
{
return previous_velocity;
}

void Basic_Agent::simulate_secretion_and_uptake( Microenvironment* pS, double dt )
const std::vector<double>& Basic_Agent::get_previous_velocity( void ) const
{
return previous_velocity;
}

int Basic_Agent::get_ID() const
{
return ID;
}

void Basic_Agent::set_ID(int new_ID)
{
ID = new_ID;
}

int Basic_Agent::get_index() const
{
return index;
}

void Basic_Agent::set_index(int new_index)
{
index = new_index;
}

int Basic_Agent::get_type() const
{
return type;
}

void Basic_Agent::set_type(int new_type)
{
type = new_type;
}

bool Basic_Agent::get_is_active() const
{
return is_active;
}

void Basic_Agent::set_is_active(bool active)
{
is_active = active;
}

double* Basic_Agent::get_secretion_rates()
{
return secretion_rates.data();
}

const double* Basic_Agent::get_secretion_rates() const
{
return secretion_rates.data();
}

double* Basic_Agent::get_saturation_densities()
{
return saturation_densities.data();
}

const double* Basic_Agent::get_saturation_densities() const
{
return saturation_densities.data();
}

double* Basic_Agent::get_uptake_rates()
{
return uptake_rates.data();
}

const double* Basic_Agent::get_uptake_rates() const
{
return uptake_rates.data();
}

double* Basic_Agent::get_net_export_rates()
{
return net_export_rates.data();
}

const double* Basic_Agent::get_net_export_rates() const
{
return net_export_rates.data();
}

double* Basic_Agent::get_internalized_total_substrates()
{
return internalized_substrates.data();
}

const double* Basic_Agent::get_internalized_total_substrates() const
{
return internalized_substrates.data();
}

double* Basic_Agent::get_fraction_released_at_death()
{
return fraction_released_at_death.data();
}

const double* Basic_Agent::get_fraction_released_at_death() const
{
return fraction_released_at_death.data();
}

double* Basic_Agent::get_fraction_transferred_when_ingested()
{
return fraction_transferred_when_ingested.data();
}

const double* Basic_Agent::get_fraction_transferred_when_ingested() const
{
return fraction_transferred_when_ingested.data();
}

Microenvironment_Interface* Basic_Agent::get_microenvironment_interface( void )
{
return microenvironment;
}

void Basic_Agent::simulate_secretion_and_uptake( double dt )
{
if(!is_active)
{ return; }
Expand All @@ -326,22 +479,22 @@ void Basic_Agent::simulate_secretion_and_uptake( Microenvironment* pS, double dt
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 *= (*microenvironment)(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
total_extracellular_substrate_change *= microenvironment->voxels(current_voxel_index).volume; // W*((1-c2)*rho+c1)/c2

*internalized_substrates -= total_extracellular_substrate_change; // opposite of net extracellular change
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;
(*microenvironment)(current_voxel_index) += cell_source_sink_solver_temp1;
(*microenvironment)(current_voxel_index) /= cell_source_sink_solver_temp2;

// now do net export
(*pS)(current_voxel_index) += cell_source_sink_solver_temp_export2;
(*microenvironment)(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;
internalized_substrates -= cell_source_sink_solver_temp_export1;
}

return;
Expand Down
Loading