Commit ce105e26 authored by gsaxena's avatar gsaxena
Browse files

(1) Cleaned up TNF MPI example

(2) Have NOT touched any function in any core file - need to address functions like display_simulation_status() so that printing is done only by a single process.
(3) When printing from within any custom function, have put a note that multiple processes will print, to prevent make another custom function which uses the "world" object to control printing.
parent e9af88d5
......@@ -65,6 +65,18 @@
###############################################################################
*/
/*================================================================================
+ If you use PhysiCell-X in your project, we would really appreciate if you can +
+ +
+ [1] Cite the PhysiCell-X repository by giving its URL +
+ +
+ [2] Cite BioFVM-X: +
+ Saxena, Gaurav, Miguel Ponce-de-Leon, Arnau Montagud, David Vicente Dorca, +
+ and Alfonso Valencia. "BioFVM-X: An MPI+ OpenMP 3-D Simulator for Biological +
+ Systems." In International Conference on Computational Methods in Systems +
+ Biology, pp. 266-279. Springer, Cham, 2021. +
=================================================================================*/
#include "./custom.h"
#include "../DistPhy/DistPhy_Utils.h"
#include "../DistPhy/DistPhy_Collective.h"
......@@ -82,7 +94,7 @@ void create_cell_types(void)
SeedRandom(parameters.ints("random_seed")); // or specify a seed here
// housekeeping
std::cout << cell_defaults.name << std::endl;
// std::cout << cell_defaults.name << std::endl; <---- will print using all processes
cell_defaults.functions.volume_update_function = standard_volume_update_function;
cell_defaults.functions.update_velocity = standard_update_cell_velocity;
......@@ -100,7 +112,7 @@ void create_cell_types(void)
// No custom rule needed
cell_defaults.functions.custom_cell_rule = NULL;
/*
/*--------------------------------------------------------------------------------------
update_pc_parameters_O2_based flag controls how is the update_phenotypes
if update_pc_parameters_O2_based is set to true the model will call
tumor_cell_phenotype_with_signaling which just call two functions
......@@ -108,31 +120,30 @@ void create_cell_types(void)
2) tnf_bm_interface_main; the former updates growth and death rates
based on oxygen while the second is the function that update the boolean model.
If the flag is false then only the tnf_bm_interface_main is invoked
*/
---------------------------------------------------------------------------------------*/
if (parameters.bools("update_pc_parameters_O2_based"))
{
cell_defaults.functions.update_phenotype = tumor_cell_phenotype_with_signaling; //No parallelization needed
cell_defaults.functions.update_phenotype = tumor_cell_phenotype_with_signaling;
}
else
{
cell_defaults.functions.update_phenotype = tnf_bm_interface_main; //No parallelization needed
cell_defaults.functions.update_phenotype = tnf_bm_interface_main;
}
cell_defaults.functions.add_cell_basement_membrane_interactions = NULL;
cell_defaults.functions.calculate_distance_to_membrane = NULL;
cell_defaults.functions.set_orientation = NULL;
/*
This parses the cell definitions in the XML config file.
*/
/* This parses the cell definitions in the XML config file. */
initialize_cell_definitions_from_pugixml(); //Only serial version exists and is sufficient.
// initialize tnf
// the following should be changed
//commented today
tnf_receptor_model_setup(); //I think no need to parallelize
tnf_boolean_model_interface_setup(); //I think no need to parallelize
submodel_registry.display(std::cout);
tnf_receptor_model_setup();
tnf_boolean_model_interface_setup();
//submodel_registry.display(std::cout); <------ will print using all processes
// Needs to initialize one of the receptor state to the total receptor value
cell_defaults.custom_data["unbound_external_TNFR"] = cell_defaults.custom_data["TNFR_receptors_per_cell"];
......@@ -140,7 +151,15 @@ void create_cell_types(void)
cell_defaults.custom_data["bound_internal_TNFR"] = 0;
build_cell_definitions_maps();
display_cell_definitions(std::cout);
//display_cell_definitions(std::cout); <------ will print using all processes
/*------------------------------------------------------------------------------------*/
/* Multiple print statements have been disabled above as the printing will be done by */
/* ALL processes. If you want to print (for checking, testing etc.) then create a new */
/* version of void create_cell_types(mpi_Environment &world, mpi_Cartesian &cart_topo)*/
/* and then use the IOProcessor(world) function to print using ONLY the root process */
/*------------------------------------------------------------------------------------*/
return;
}
......@@ -290,8 +309,8 @@ void tumor_cell_phenotype_with_signaling(Cell *pCell, Phenotype &phenotype, doub
tnf_bm_interface_main(pCell, phenotype, dt);
}
//different in physiboss_cell_lines
// cell coloring function for ploting the svg files
std::vector<std::string> my_coloring_function(Cell *pCell)
{
// start with live coloring
......@@ -318,8 +337,8 @@ std::vector<std::string> my_coloring_function(Cell *pCell)
return output;
}
// Funtion to read init files created with PhysiBoSSv2
//unneaded
// Function to read init files created with PhysiBoSSv2
std::vector<init_record> read_init_file(std::string filename, char delimiter, bool header)
{
// File pointer
......@@ -372,16 +391,20 @@ std::vector<init_record> read_init_file(std::string filename, char delimiter, bo
return result;
}
// aded from example
void set_input_nodes(Cell* pCell) {}
void from_nodes_to_cell(Cell* pCell, Phenotype& phenotype, double dt) {}
void color_node(Cell* pCell){
void color_node(Cell* pCell)
{
std::string node_name = parameters.strings("node_to_visualize");
pCell->custom_data[node_name] = pCell->phenotype.intracellular->get_boolean_variable_value(node_name);
}
void inject_density_sphere(int density_index, double concentration, double membrane_lenght)
{
// Inject given concentration on the extremities only
#pragma omp parallel for
for (int n = 0; n < microenvironment.number_of_voxels(); n++)
......@@ -401,39 +424,6 @@ void remove_density(int density_index)
std::cout << "Removal done" << std::endl;
}
// std::string cells_message_builder(std::vector<Cell*> all_cells, double timepoint){
// // implement count different situation of cells
// int alive_no,apoptotic_no,necrotic_no=0;
// // int p_id = all_cells[i]->ID
// static int TUMOR_TYPE=0;
// std::string message;
// for(int i=0;i<all_cells.size();i++)
// {
// if( all_cells[i]->phenotype.cycle.pCycle_Model )
// {
// int code= all_cells[i]->phenotype.cycle.current_phase().code;
// if (code ==PhysiCell_constants::Ki67_positive_premitotic || code==PhysiCell_constants::Ki67_positive_postmitotic || code==PhysiCell_constants::Ki67_positive || code==PhysiCell_constants::Ki67_negative || code==PhysiCell_constants::live)
// // _nameCore="LIVE";
// alive_no +=1;
// else if (code==PhysiCell_constants::apoptotic)
// // _nameCore="APOP";
// apoptotic_no +=1;
// else if (code==PhysiCell_constants::necrotic_swelling || code==PhysiCell_constants::necrotic_lysed || code==PhysiCell_constants::necrotic)
// // _nameCore="NEC";
// necrotic_no += 1;
// // else if (code==PhysiCell_constants::debris)
// // _nameCore="DEBR";
// // else
// // _nameCore="MISC";
// }
// else if(all_cells[i]->type==TUMOR_TYPE)
// // _nameCore="LIVE";
// alive_no+=1;
// }
// // message = std::to_string(p_id) + ';' + std::to_string(timepoint) + ';' + std::to_string(alive_no) + ';' + std::to_string(apoptotic_no) + ';' + std::to_string(necrotic_no) + ';';
// message = std::to_string(timepoint) + ';' + std::to_string(alive_no) + ';' + std::to_string(apoptotic_no) + ';' + std::to_string(necrotic_no) + ';';
// return message;
// }
double total_live_cell_count()
{
......
......@@ -65,12 +65,24 @@
###############################################################################
*/
/*================================================================================
+ If you use PhysiCell-X in your project, we would really appreciate if you can +
+ +
+ [1] Cite the PhysiCell-X repository by giving its URL +
+ +
+ [2] Cite BioFVM-X: +
+ Saxena, Gaurav, Miguel Ponce-de-Leon, Arnau Montagud, David Vicente Dorca, +
+ and Alfonso Valencia. "BioFVM-X: An MPI+ OpenMP 3-D Simulator for Biological +
+ Systems." In International Conference on Computational Methods in Systems +
+ Biology, pp. 266-279. Springer, Cham, 2021. +
=================================================================================*/
#include "../core/PhysiCell.h"
#include "../modules/PhysiCell_standard_modules.h"
#include "./tnf_receptor_dynamics.h"
#include "./tnf_boolean_model_interface.h"
//change
#include "../addons/PhysiBoSS/src/maboss_network.h"
using namespace BioFVM;
......@@ -127,7 +139,7 @@ void inject_density_sphere(int density_index, double concentration, double membr
// helper function to remove a density
void remove_density( int density_index );
//function to create a message
// function to create a message
// std::string cells_message_builder(std::vector<Cell*> all_cells, double timepoint);
double total_live_cell_count();
......
......@@ -85,16 +85,16 @@ void update_cell_from_boolean_model(Cell* pCell, Phenotype& phenotype, double dt
if ( nonACD )
{
pCell->start_death(necrosis_model_index); //In PhysiCell_cell.cpp - no parallelization needed I think
pCell->start_death(necrosis_model_index);
return;
}
if ( survival && pCell->phenotype.cycle.current_phase_index() == PhysiCell_constants::Ki67_negative )
{
pCell->phenotype.cycle.advance_cycle(pCell, phenotype, dt); //In PhysiCell_phenotype.cpp - don't need parallelization I think
pCell->phenotype.cycle.advance_cycle(pCell, phenotype, dt);
}
// // If NFkB node is active produce some TNF
// If NFkB node is active produce some TNF
if ( NFkB )
{
phenotype.secretion.net_export_rates[nTNF_external] = pCell->custom_data[nTNF_export_rate];
......@@ -129,19 +129,19 @@ void tnf_bm_interface_main(Cell* pCell, Phenotype& phenotype, double dt)
return;
}
if ( pCell->phenotype.intracellular->need_update() ) //No parallelization needed I think
if ( pCell->phenotype.intracellular->need_update() )
{
// First we update the Boolean Model inputs
update_boolean_model_inputs(pCell, phenotype, dt ); //No parallelization needed I think
update_boolean_model_inputs(pCell, phenotype, dt );
// Run maboss to update the boolean state of the cell
pCell->phenotype.intracellular->update(); //No parallelization needed I think
pCell->phenotype.intracellular->update();
// update the cell fate based on the boolean outputs
update_cell_from_boolean_model(pCell, phenotype, dt); //No parallelization needed I think
update_cell_from_boolean_model(pCell, phenotype, dt);
// Get track of some boolean node values for debugging
update_monitor_variables(pCell); //No parallelization needed I think
update_monitor_variables(pCell);
}
return;
......
......@@ -27,21 +27,21 @@ Submodel_Information tnf_receptor_info;
void tnf_receptor_model_setup()
{
tnf_receptor_info.name = "TNF trasnporter model";
tnf_receptor_info.name = "TNF trasnporter model";
tnf_receptor_info.version = "0.1.0";
tnf_receptor_info.main_function = tnf_receptor_model;
tnf_receptor_info.main_function = tnf_receptor_model;
// what custom data do I need?
tnf_receptor_info.cell_variables.push_back( "TNFR_activation_threshold" );
tnf_receptor_info.cell_variables.push_back( "unbound_external_TNFR" );
tnf_receptor_info.cell_variables.push_back( "unbound_external_TNFR" );
tnf_receptor_info.cell_variables.push_back( "bound_external_TNFR" );
tnf_receptor_info.cell_variables.push_back( "bound_internal_TNFR" );
tnf_receptor_info.cell_variables.push_back( "TNFR_binding_rate" );
tnf_receptor_info.cell_variables.push_back( "TNFR_binding_rate" );
tnf_receptor_info.cell_variables.push_back( "TNFR_endocytosis_rate" );
tnf_receptor_info.cell_variables.push_back( "TNFR_recycling_rate" );
tnf_receptor_info.cell_variables.push_back( "TNFR_recycling_rate" );
tnf_receptor_info.cell_variables.push_back( "TFN_net_production_rate" );
tnf_receptor_info.register_model();
......@@ -53,11 +53,11 @@ void tnf_receptor_model( Cell* pCell, Phenotype& phenotype, double dt )
{
static int nTNF_external = microenvironment.find_density_index( "tnf" );
static int nR_EU = pCell->custom_data.find_variable_index( "unbound_external_TNFR" );
static int nR_EU = pCell->custom_data.find_variable_index( "unbound_external_TNFR" );
static int nR_EB = pCell->custom_data.find_variable_index( "bound_external_TNFR" );
static int nR_IB = pCell->custom_data.find_variable_index( "bound_internal_TNFR" );
static int nR_bind = pCell->custom_data.find_variable_index( "TNFR_binding_rate" );
static int nR_bind = pCell->custom_data.find_variable_index( "TNFR_binding_rate" );
static double R_binding_rate = pCell->custom_data[nR_bind];
static int nR_endo = pCell->custom_data.find_variable_index( "TNFR_endocytosis_rate" );
......@@ -70,7 +70,7 @@ void tnf_receptor_model( Cell* pCell, Phenotype& phenotype, double dt )
if( phenotype.death.dead == true )
{ return; }
// internalized TNF tells us how many have recently bound to receptors
// internalized TNF tells us how many have recently bound to receptors
// TNF is internalized at:
// phenotype.secretion.uptake_rates[nTNF_external] =
// pCell->custom_data[nR_bind] * pCell->custom_data[nR_EU];
......@@ -105,8 +105,8 @@ void tnf_receptor_model( Cell* pCell, Phenotype& phenotype, double dt )
// Recylcing
// The internalized bounded TNFR release the TNF
// The TNF is instantaneously degraded by the cell
// The TNF receptor is recycled as an unbounded external receptor
// The TNF is instantaneously degraded by the cell
// The TNF receptor is recycled as an unbounded external receptor
double dR_EU = dt * R_recyc_rate * pCell->custom_data[nR_IB];
if( dR_EU > pCell->custom_data[nR_IB] )
{ dR_EU = pCell->custom_data[nR_IB]; }
......
......@@ -65,6 +65,18 @@
###############################################################################
*/
/*================================================================================
+ If you use PhysiCell-X in your project, we would really appreciate if you can +
+ +
+ [1] Cite the PhysiCell-X repository by giving its URL +
+ +
+ [2] Cite BioFVM-X: +
+ Saxena, Gaurav, Miguel Ponce-de-Leon, Arnau Montagud, David Vicente Dorca, +
+ and Alfonso Valencia. "BioFVM-X: An MPI+ OpenMP 3-D Simulator for Biological +
+ Systems." In International Conference on Computational Methods in Systems +
+ Biology, pp. 266-279. Springer, Cham, 2021. +
=================================================================================*/
/*=======================================*/
/* Include mpi.h in the parallel version */
/*=======================================*/
......@@ -81,14 +93,11 @@
#include "./core/PhysiCell.h"
#include "./modules/PhysiCell_standard_modules.h"
#include "./addons/PhysiBoSS/src/maboss_intracellular.h"
// put custom code modules here!
#include "./custom_modules/custom.h"
// #include <cppkafka/utils/buffered_producer.h>
// #include <cppkafka/configuration.h>
using namespace BioFVM;
using namespace PhysiCell;
// using namespace cppkafka;
/*====================================================*/
/* Use DistPhy::mpi namespace in the parallel version */
......@@ -97,14 +106,13 @@ using namespace DistPhy::mpi;
int main( int argc, char* argv[] )
{
// load and parse settings file(s)
/*=======================================================================================*/
/* Create mpi_Environment object, initialize it, then create Cartesian Topology */
/*=======================================================================================*/
mpi_Environment world; //object contains size of communicator, rank of process
world.Initialize(); //Initialize using MPI_THREAD_MULTIPLE, find comm. size and comm. rank
mpi_Environment world; //Object contains size of communicator, rank of process
world.Initialize(); //Initialize using MPI_THREAD_FUNNELED, find comm. size and comm. rank
mpi_Cartesian cart_topo; //Contains dims[3], ccoords[3] array and MPI_Comm mpi_cart_comm
cart_topo.Build_Cartesian_Topology(world); //Create 1-D X decomposition by setting dims[1]=size.
cart_topo.Find_Cartesian_Coordinates(world); //Find Cartesian Topology coordinates of each process
......@@ -130,197 +138,149 @@ int main( int argc, char* argv[] )
exit(-1);
}
// OpenMP setup, commented out for debugging
//omp_set_num_threads(PhysiCell_settings.omp_num_threads);
/*=======================================================================================*/
/* Setting the number of threads here will take precedence over the threads specified in */
/* PhysiCell_settings.xml and the environment variable OMP_NUM_THREADS */
/*=======================================================================================*/
// PNRG setup
SeedRandom( parameters.ints("random_seed") ); // Or a seed can be specified here
//omp_set_num_threads(PhysiCell_settings.omp_num_threads); <--- We use OMP_NUM_THREADS
// time setup
SeedRandom( parameters.ints("random_seed") ); // Or a seed can be specified here
std::string time_units = "min";
/* Microenvironment setup */
/*================================================================================================*/
/* Objects of both DistPhy_Environment and DistPhy_Cartesian are needed in setup microenvironment */
/* These objects are passed to initialize_microenvironment() which calls resize_uniform() */
/*================================================================================================*/
setup_microenvironment(world, cart_topo); // modify this in the custom code
// User parameters
/*========================*/
/* Microenvironment setup */
/*========================*/
setup_microenvironment(world, cart_topo);
//User parameters
double time_add_tnf = parameters.ints("time_add_tnf");
double time_put_tnf = 0;
double duration_add_tnf = parameters.ints("duration_add_tnf");
double time_tnf_next = 0;
double time_remove_tnf = parameters.ints("time_remove_tnf");
// change concentration units too match voxel volume
// double concentration_tnf = parameters.doubles("concentration_tnf") * microenvironment.voxels(0).volume * 0.000001;
double concentration_tnf = parameters.doubles("concentration_tnf") * 0.1;
// radius around which the tnf pulse is injected
//Radius around which the tnf pulse is injected
double membrane_lenght = parameters.ints("membrane_length");
// tnf density index
static int tnf_idx = microenvironment.find_density_index("tnf");
// this is to emulate PhysiBoSSv1 TNF experiment
/*==============================================================================*/
/* seed_tnf = false in the next statement so the if block (compound statement) */
/* is NOT executed. If seed_tnf becomes true THEN we need to check if the */
/* the function called within the if conditions i.e. inject_density_sphere */
/* is to be parallelized or not. I don't think it needs parallelization so only */
/* changing the simulate_diffusion_decay(...) call to parallel version */
/*==============================================================================*/
//TNF density index
static int tnf_idx = microenvironment.find_density_index("tnf");
bool seed_tnf = false;
// do small diffusion steps alone to initialize densities
//Do small diffusion steps to initialize densities in case seed_tnf is true
if ( seed_tnf )
{
inject_density_sphere(tnf_idx, concentration_tnf, membrane_lenght);
for ( int i = 0; i < 25; i ++ )
microenvironment.simulate_diffusion_decay( diffusion_dt, world, cart_topo );
}
/* PhysiCell setup */
/*=============================*/
/* PhysiCell/PhysiCell-X setup */
/*=============================*/
// set mechanics voxel size, and match the data structure to BioFVM
//GS changed it to 25 otherwise it was 30 because 500 - (-500) / 25 = 40
// Mechanical voxel size must be >= Diffusion voxel size (check with PhysiCell_settings.xml)
double mechanics_voxel_size = 20;
/*---------------------------------------------------------*/
/* Calling the parallel version of Cell Container creation */
/*---------------------------------------------------------*/
// Calling the parallel version of Cell Container creation
Cell_Container* cell_container = create_cell_container_for_microenvironment( microenvironment, mechanics_voxel_size, world, cart_topo );
/* Users typically start modifying here. START USERMODS */
//----> Users typically start modifying here. START USERMODS
create_cell_types();
/* Calling the parallel version of setup_tissue(...) */
//Calling the parallel version of setup_tissue(...)
setup_tissue(microenvironment, world, cart_topo);
/* Users typically stop modifying here. END USERMODS */
//----> Users typically stop modifying here. END USERMODS
// set MultiCellDS save options
//Set MultiCellDS save options <--- These MUST all be 'true' here
set_save_biofvm_mesh_as_matlab( true );
set_save_biofvm_data_as_matlab( true );
set_save_biofvm_cell_data( true );
set_save_biofvm_cell_data_as_custom_matlab( true );
// save a simulation snapshot
//Save a simulation snapshot
char filename[1024];
sprintf( filename , "%s/initial" , PhysiCell_settings.folder.c_str() );
/* Use the parallel version of the function now */
//Use the parallel version of the function for XML file writing
save_PhysiCell_to_MultiCellDS_xml_pugi( filename , microenvironment , PhysiCell_globals.current_time, world, cart_topo );
// save a quick SVG cross section through z = 0, after setting its
// length bar to 200 microns
//Save SVG cross section through z = 0, after setting its length bar to 200 microns
PhysiCell_SVG_options.length_bar = 200;
// for simplicity, set a pathology coloring function
//For simplicity, set a pathology coloring function
std::vector<std::string> (*cell_coloring_function)(Cell*) = my_coloring_function;
sprintf( filename , "%s/initial.svg" , PhysiCell_settings.folder.c_str() );
SVG_plot( filename , microenvironment, 0.0 , PhysiCell_globals.current_time, cell_coloring_function, world, cart_topo );
display_citations();
sprintf( filename , "%s/initial.svg" , PhysiCell_settings.folder.c_str() );
// set the performance timers
//Use the parallel version of the function for SVG plot file
SVG_plot( filename , microenvironment, 0.0 , PhysiCell_globals.current_time, cell_coloring_function, world, cart_topo );
//Set the performance timers
BioFVM::RUNTIME_TIC();
BioFVM::TIC();
std::ofstream report_file;
if( PhysiCell_settings.enable_legacy_saves == true )
{
sprintf( filename , "%s/simulation_report.txt" , PhysiCell_settings.folder.c_str() );
report_file.open(filename); // create the data log file
report_file<<"simulated time\tnum cells\tnum division\tnum death\twall time"<<std::endl;
}
// main loop
//Main loop of the program
try
{
while( PhysiCell_globals.current_time < PhysiCell_settings.max_time + 0.1*diffusion_dt )
{
// save data if it's time.
//Save data if it's time.
if( fabs( PhysiCell_globals.current_time - PhysiCell_globals.next_full_save_time ) < 0.01 * diffusion_dt )
{
std::cout << "Time to save" << std::endl;
display_simulation_status( std::cout );
if( PhysiCell_settings.enable_legacy_saves == true )
{
//GS commented out the next statement
//log_output( PhysiCell_globals.current_time , PhysiCell_globals.full_output_index, microenvironment, report_file);
//Count Necrotic Apoptotic Alive cells
// Producer
if(IOProcessor(world))
std::cout << "Time to save" << std::endl;
display_simulation_status( std::cout );
if( PhysiCell_settings.enable_legacy_saves == true )
{
//Count Necrotic, Apoptotic and Alive cells
std::string message;
std::string topic_name = "cells";
double timepoint = PhysiCell_globals.current_time;
int alive_no,necrotic_no,apoptotic_no;
/* Call the parallel versions of the function now which use MPI_Reduce at rank 0 */
//Call the parallel versions of the function now which use MPI_Reduce at rank 0
alive_no = total_live_cell_count(world, cart_topo);
necrotic_no = total_necrosis_cell_count(world, cart_topo);
apoptotic_no = total_dead_cell_count(world, cart_topo);
pid_t pid_var = getpid();
// MessageBuilder builder(topic_name);
if(IOProcessor(world)) //This is rank 0 which will gather no. of cell types from all processes
//Create number of cell types message on root process (use if desired)
if(IOProcessor(world))
message = std::to_string(pid_var) + ';' + std::to_string(timepoint) + ';' + std::to_string(alive_no) + ';' + std::to_string(apoptotic_no) + ';' + std::to_string(necrotic_no) + ';';
// message = '{' + 'process_id' + std::to_string(pid_var) + ',' + 'timepoint' + std::to_string(timepoint) + ','
// + 'alive' + std::to_string(alive_no) + ',' + 'apoptotic' + std::to_string(apoptotic_no) + ','
// + 'necrotic' + std::to_string(necrotic_no) + '}';
// Define the configuration structure
// Configuration config = { { "metadata.broker.list", "localhost:9092" } };