diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index ca40b72f7..136cd7941 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -64,4 +64,14 @@ jobs: - name: Check ${{ matrix.projects.name }} project simulation results if: matrix.projects.output_folder != '' run: | - python beta/test_diff_svg.py ${{ matrix.projects.output_folder }} tests/cases/output_${{ matrix.projects.project }} \ No newline at end of file + python beta/test_diff_svg.py ${{ matrix.projects.output_folder }} tests/cases/output_${{ matrix.projects.project }} + + start_and_stop_tests: + runs-on: ubuntu-latest + steps: + - name: Checkout repository + uses: actions/checkout@v4 + + - name: Run start and stop test + run: | + python beta/start_stop_tests.py \ No newline at end of file diff --git a/.gitignore b/.gitignore index c6e78d6e8..25744e2cf 100644 --- a/.gitignore +++ b/.gitignore @@ -19,3 +19,4 @@ project studio_debug.log user_projects/* !user_projects/empty.txt +output/* diff --git a/BioFVM/BioFVM_microenvironment.cpp b/BioFVM/BioFVM_microenvironment.cpp index 47e3d5621..c3b9942f4 100644 --- a/BioFVM/BioFVM_microenvironment.cpp +++ b/BioFVM/BioFVM_microenvironment.cpp @@ -50,6 +50,11 @@ #include "BioFVM_solvers.h" #include "BioFVM_vector.h" #include +#include +#include +#include +#include +#include "../core/PhysiCell_phenotype.h" #include "BioFVM_basic_agent.h" @@ -1613,4 +1618,861 @@ void get_row_from_substrate_initial_condition_csv(std::vector &voxel_set, c } voxel_set.push_back(voxel_ind); } + +std::ostream& operator<<(std::ostream& os, Microenvironment& S) +{ + os << std::endl << "Microenvironment summary: " << std::endl; + os << S.name << ": " << std::endl; + + os << "Microenvironment_private_data:" << std::endl; + //temporary_density_vectors1 + os << "temporary_density_vectors1:" << std::endl; + for (const auto& temporary_density_vectors1_v : S.temporary_density_vectors1) { + for( const auto& temporary_density_vectors1 : temporary_density_vectors1_v){ + os << temporary_density_vectors1 << " "; + }os << std::endl;} + os << std::endl; + + //temporary_density_vectors + os << "temporary_density_vectors2:" << std::endl; + for (const auto& temporary_density_vectors2_v : S.temporary_density_vectors2) { + for( const auto& temporary_density_vectors2 : temporary_density_vectors2_v){ + os << temporary_density_vectors2 << " "; + }os << std::endl;} + os << std::endl; + + //bulk_source_sink_solver_temp1 + os << "bulk_source_sink_solver_temp1:" << std::endl; + for (const auto& bulk_source_sink_solver_temp1_v : S.bulk_source_sink_solver_temp1) { + for( const auto& bulk_source_sink_solver_temp1 : bulk_source_sink_solver_temp1_v){ + os << bulk_source_sink_solver_temp1 << " "; + }os << std::endl;} + os << std::endl; + + //bulk_source_sink_solver_temp2 + os << "bulk_source_sink_solver_temp2:" << std::endl; + for (const auto& bulk_source_sink_solver_temp2_v : S.bulk_source_sink_solver_temp2) { + for( const auto& bulk_source_sink_solver_temp2 : bulk_source_sink_solver_temp2_v){ + os << bulk_source_sink_solver_temp2 << " "; + }os << std::endl;} + os << std::endl; + + //bulk_source_sink_solver_tem3 + os << "bulk_source_sink_solver_temp3:" << std::endl; + for (const auto& bulk_source_sink_solver_temp3_v : S.bulk_source_sink_solver_temp3) { + for( const auto& bulk_source_sink_solver_temp3 : bulk_source_sink_solver_temp3_v){ + os << bulk_source_sink_solver_temp3 << " "; + }os << std::endl;} + os << std::endl; + + //bulk_source_sink_solver_setup_done + os << "bulk_source_sink_solver_setup_done: " << (S.bulk_source_sink_solver_setup_done ? "true" : "false") << std::endl; + + //gradient_vectors + os << "gradient_vectors:" << std::endl; + for (const auto& gradient_vectors_v : S.gradient_vectors) { + for( const auto& gradient_vectors_n : gradient_vectors_v){ + for( const auto& gradient_vectors : gradient_vectors_n){ + os << gradient_vectors << " ";} + }os << std::endl;} + os << std::endl; + + //gradient_vector_computed + os << "gradient_vector_computed: "; + for (const auto& gradient_vector_computed : S.gradient_vector_computed) { + os << (gradient_vector_computed ? "true" : "false") << " "; + } + os << std::endl; + + //one + os << "one: "; + for (const auto& one : S.one) { + os << one << " "; + } + os << std::endl; + + //zero + os << "zero: "; + for (const auto& zero : S.zero) { + os << zero << " "; + } + os << std::endl; + + //one_half + os << "one_half: "; + for (const auto& one_half : S.one_half) { + os << one_half << " "; + } + os << std::endl; + + //one_third + os << "one_third: "; + for (const auto& one_third : S.one_third) { + os << one_third << " "; + } + os << std::endl; + + //thomas_temp1 + os << "thomas_temp1:" << std::endl; + for (const auto& thomas_temp1_v : S.thomas_temp1) { + for( const auto& thomas_temp1 : thomas_temp1_v){ + os << thomas_temp1 << " "; + }os << std::endl;} + os << std::endl; + + //thomas_temp2 + os << "thomas_temp2:" << std::endl; + for (const auto& thomas_temp2_v : S.thomas_temp2) { + for( const auto& thomas_temp2 : thomas_temp2_v){ + os << thomas_temp2 << " "; + }os << std::endl;} + os << std::endl; + + //thomas_constant1x + os << "thomas_constant1x: "; + for (const auto& thomas_constant1x : S.thomas_constant1x) { + os << thomas_constant1x << " "; + } + os << std::endl; + + //thomas_constant1y + os << "thomas_constant1y: "; + for (const auto& thomas_constant1y : S.thomas_constant1y) { + os << thomas_constant1y << " "; + } + os << std::endl; + + //thomas_constant1z + os << "thomas_constant1z: "; + for (const auto& thomas_constant1z : S.thomas_constant1z) { + os << thomas_constant1z << " "; + } + os << std::endl; + + //thomas_neg_constant1x + os << "thomas_neg_constant1x: "; + for (const auto& thomas_neg_constant1x : S.thomas_neg_constant1x) { + os << thomas_neg_constant1x << " "; + } + os << std::endl; + + //thomas_neg_constant1y + os << "thomas_neg_constant1y: "; + for (const auto& thomas_neg_constant1y : S.thomas_neg_constant1y) { + os << thomas_neg_constant1y << " "; + } + os << std::endl; + + //thomas_neg_constant1z + os << "thomas_neg_constant1z: "; + for (const auto& thomas_neg_constant1z : S.thomas_neg_constant1z) { + os << thomas_neg_constant1z << " "; + } + os << std::endl; + + //thomas_setup_done + os << "thomas_setup_done: " << (S.thomas_setup_done ? "true" : "false") << std::endl; + + //thomas_i_jump + os << "thomas_i_jump: " << S.thomas_i_jump << std::endl; + + //thomas_i_jump + os << "thomas_j_jump: " << S.thomas_j_jump << std::endl; + + //thomas_i_jump + os << "thomas_k_jump: " << S.thomas_k_jump << std::endl; + + //thomas_constant1 + os << "thomas_constant1: "; + for (const auto& thomas_constant1 : S.thomas_constant1) { + os << thomas_constant1 << " "; + } + os << std::endl; + + //thomas_constant1a + os << "thomas_constant1a: "; + for (const auto& thomas_constant1a : S.thomas_constant1a) { + os << thomas_constant1a << " "; + } + os << std::endl; + + //thomas_constant2 + os << "thomas_constant2: "; + for (const auto& thomas_constant2 : S.thomas_constant2) { + os << thomas_constant2 << " "; + } + os << std::endl; + + //thomas_constant3 + os << "thomas_constant3: "; + for (const auto& thomas_constant3 : S.thomas_constant3) { + os << thomas_constant3 << " "; + } + os << std::endl; + + //thomas_constant3a + os << "thomas_constant3a: "; + for (const auto& thomas_constant3a : S.thomas_constant3a) { + os << thomas_constant3a << " "; + } + os << std::endl; + + //name + os << "name: " << S.name << std::endl; + + //supply_rates + os << "supply_rates:" << std::endl; + for (const auto& supply_rates_v : S.supply_rates) { + for( const auto& supply_rates : supply_rates_v){ + os << supply_rates << " "; + }os << std::endl;} + os << std::endl; + + //uptake_rates + os << "uptake_rates:" << std::endl; + for (const auto& uptake_rates_v : S.uptake_rates) { + for( const auto& uptake_rates : uptake_rates_v){ + os << uptake_rates << " "; + }os << std::endl;} + os << std::endl; + + os << "Densities: " << S.number_of_densities() << std::endl; + for (unsigned int i = 0; i < S.density_names.size(); i++) + { + os << S.density_names[i] << std::endl + << "units: " << S.density_units[i] << std::endl + << "diffusion_coefficient: " << S.diffusion_coefficients[i] + << " " << S.spatial_units << " " << S.time_units << std::endl + << "decay_rate: " << S.decay_rates[i] + << " " << S.time_units << std::endl; + } + os << "\n"; + // Display all the densities for each voxel. + os << "Values for each voxel: " << std::endl; + for (unsigned int n = 0; n < S.number_of_voxels(); n++) + { + os << S.density_vector(n) << std::endl; + } + os << std::endl; + + os << "centers of voxels: " << std::endl; + for (unsigned int n = 0; n < S.number_of_voxels(); n++) + { + os << S.voxels(n).center << std::endl; + } + os << std::endl; + + // dirichlet_value_vectors + os << "dirichlet_value_vectors:" << std::endl; + for (const auto& vec : S.dirichlet_value_vectors) { + for (const auto& val : vec) { + os << val << " "; + } + os << std::endl; + } + os << std::endl; + + // dirichlet_activation_vector (singolo) + os << "dirichlet_activation_vector: "; + for (bool val : S.dirichlet_activation_vector) { + os << (val ? "true" : "false") << " "; + } + os << std::endl << std::endl; + + // dirichlet_activation_vectors (matrice) + os << "dirichlet_activation_vectors:" << std::endl; + for (const auto& vec : S.dirichlet_activation_vectors) { + for (bool val : vec) { + os << (val ? "true" : "false") << " "; + } + os << std::endl; + } + os << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Microenvironment& microenv){ + + std::string dummy; + double value; + int j = 0; + std::string value_2; + + // read headers + std::getline(is, dummy); + std::getline(is, dummy); + std::getline(is, dummy); + std::getline(is, dummy); + + //temporary_density_vectors1 + std::getline(is, dummy); + + unsigned int n_lines = 0; + + for( unsigned int n=0; n < microenv.temporary_density_vectors1.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.temporary_density_vectors1[n][j] = value; + j = j+1; + } + } + + + + assert(n_lines == microenv.temporary_density_vectors1.size()); + //read empty line + std::getline(is, dummy); + + //temporary_density_vectors2 + std::getline(is, dummy); + + n_lines = 0; + for( unsigned int n=0; n < microenv.temporary_density_vectors2.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.temporary_density_vectors2[n][j] = value; + j = j+1; + } + } + + + assert(n_lines == microenv.temporary_density_vectors2.size()); + //read empty line + std::getline(is, dummy); + + //bulk_source_sink_solver_temp1 + std::getline(is, dummy); + n_lines = 0; + for( unsigned int n=0; n < microenv.bulk_source_sink_solver_temp1.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.bulk_source_sink_solver_temp1[n][j] = value; + j = j+1; + } + } + + assert(n_lines == microenv.bulk_source_sink_solver_temp1.size()); + //read empty line + std::getline(is, dummy); + + //bulk_source_sink_solver_temp2 + std::getline(is, dummy); + n_lines = 0; + for( unsigned int n=0; n < microenv.bulk_source_sink_solver_temp2.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.bulk_source_sink_solver_temp2[n][j] = value; + j = j+1; + } + } + assert(n_lines == microenv.bulk_source_sink_solver_temp2.size()); + //read empty line + std::getline(is, dummy); + + //bulk_source_sink_solver_temp3 + std::getline(is, dummy); + n_lines = 0; + for( unsigned int n=0; n < microenv.bulk_source_sink_solver_temp3.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.bulk_source_sink_solver_temp3[n][j] = value; + j = j+1; + } + } + assert(n_lines == microenv.bulk_source_sink_solver_temp3.size()); + + //read empty line + std::getline(is, dummy); + + //bulk_source_sink_solver_setup_done + std::getline(is, dummy); + microenv.bulk_source_sink_solver_setup_done = read_number_in_line_bool(dummy); + + + //gradient_vectors + std::getline(is, dummy); + + n_lines = 0; + for (unsigned int n = 0; n < microenv.gradient_vectors.size(); n++) + { std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + for (unsigned int l = 0; l < microenv.gradient_vectors[n].size(); l++) + { + + double value; + for (int j = 0; j < microenv.gradient_vectors[n][l].size(); j++) + { + stream4 >> value; + microenv.gradient_vectors[n][l][j] = value; + } + } + } + assert(n_lines == microenv.gradient_vectors.size()); + + + //skip empty line + std::getline(is, dummy); + + //gradient_vector_computed + std::getline(is, dummy); + bool gradient_vector_computed; + std::istringstream stream_gradient_vector_computed(dummy); + std::string key_gradient_vector_computed; + stream_gradient_vector_computed >> key_gradient_vector_computed; + for (int i = 0; i < microenv.gradient_vector_computed.size(); i++){ + stream_gradient_vector_computed >> value; + if(value_2 == "true"){ + gradient_vector_computed = true; + } + else if(value_2 == "false"){ + gradient_vector_computed = false; + } + microenv.gradient_vector_computed[i] = gradient_vector_computed; + } + + // one + std::getline(is, dummy); + std::istringstream stream_one(dummy); + std::string key_one; + stream_one >> key_one; + std::vector one; + while(stream_one >> value) { + one.push_back(value); + } + microenv.one.clear(); + microenv.one = one; + + // zero + std::getline(is, dummy); + std::istringstream stream_zero(dummy); + std::string key_zero; + stream_zero >> key_zero; + std::vector zero; + while(stream_zero >> value) { + zero.push_back(value); + } + microenv.zero.clear(); + microenv.zero = zero; + + // one_half + std::getline(is, dummy); + std::istringstream stream_one_half(dummy); + std::string key_one_half; + stream_one_half >> key_one_half; + std::vector one_half; + while(stream_one_half >> value) { + one_half.push_back(value); + } + microenv.one_half.clear(); + microenv.one_half = one_half; + + // one_third + std::getline(is, dummy); + std::istringstream stream_one_third(dummy); + std::string key_one_third; + stream_one_third >> key_one_third; + std::vector one_third; + while(stream_one_third >> value) { + one_third.push_back(value); + } + microenv.one_third.clear(); + microenv.one_third = one_third; + + //thomas_temp1 + std::getline(is, dummy); + n_lines = 0; + for( unsigned int n=0; n < microenv.thomas_temp1.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.thomas_temp1[n][j] = value; + j = j+1; + } + } + + assert(n_lines == microenv.thomas_temp1.size()); + //read empty line + std::getline(is, dummy); + + //thomas_temp2 + std::getline(is, dummy); + n_lines = 0; + for( unsigned int n=0; n < microenv.thomas_temp2.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.thomas_temp2[n][j] = value; + j = j+1; + } + } + assert(n_lines == microenv.thomas_temp2.size()); + //read empty line + std::getline(is, dummy); + + // thomas_constant1x + std::getline(is, dummy); + std::istringstream stream_thomas_constant1x(dummy); + std::string key_thomas_constant1x; + stream_thomas_constant1x >> key_thomas_constant1x; + std::vector thomas_constant1x; + while(stream_thomas_constant1x >> value) { + thomas_constant1x.push_back(value); + } + microenv.thomas_constant1x.clear(); + microenv.thomas_constant1x = thomas_constant1x; + + // thomas_constant1y + std::getline(is, dummy); + std::istringstream stream_thomas_constant1y(dummy); + std::string key_thomas_constant1y; + stream_thomas_constant1y >> key_thomas_constant1y; + std::vector thomas_constant1y; + while(stream_thomas_constant1y >> value) { + thomas_constant1y.push_back(value); + } + microenv.thomas_constant1y.clear(); + microenv.thomas_constant1y = thomas_constant1y; + + // thomas_constant1z + std::getline(is, dummy); + std::istringstream stream_thomas_constant1z(dummy); + std::string key_thomas_constant1z; + stream_thomas_constant1z >> key_thomas_constant1z; + std::vector thomas_constant1z; + while(stream_thomas_constant1z >> value) { + thomas_constant1z.push_back(value); + } + microenv.thomas_constant1z.clear(); + microenv.thomas_constant1z = thomas_constant1z; + + // thomas_neg_constant1x + std::getline(is, dummy); + std::istringstream stream_thomas_neg_constant1x(dummy); + std::string key_thomas_neg_constant1x; + stream_thomas_neg_constant1x >> key_thomas_neg_constant1x; + std::vector thomas_neg_constant1x; + while(stream_thomas_neg_constant1x >> value) { + thomas_neg_constant1x.push_back(value); + } + microenv.thomas_neg_constant1x.clear(); + microenv.thomas_neg_constant1x = thomas_neg_constant1x; + + // thomas_neg_constant1y + std::getline(is, dummy); + std::istringstream stream_thomas_neg_constant1y(dummy); + std::string key_thomas_neg_constant1y; + stream_thomas_neg_constant1y >> key_thomas_neg_constant1y; + std::vector thomas_neg_constant1y; + while(stream_thomas_neg_constant1y >> value) { + thomas_neg_constant1y.push_back(value); + } + microenv.thomas_neg_constant1y.clear(); + microenv.thomas_neg_constant1y = thomas_neg_constant1y; + + // thomas_neg_constant1y + std::getline(is, dummy); + std::istringstream stream_thomas_neg_constant1z(dummy); + std::string key_thomas_neg_constant1z; + stream_thomas_neg_constant1z >> key_thomas_neg_constant1z; + std::vector thomas_neg_constant1z; + while(stream_thomas_neg_constant1z >> value) { + thomas_neg_constant1z.push_back(value); + } + microenv.thomas_neg_constant1z.clear(); + microenv.thomas_neg_constant1z = thomas_neg_constant1z; + + //thomas_setup_done + std::getline(is, dummy); + microenv.thomas_setup_done = read_number_in_line_bool(dummy); + + //thomas_i_jump + std::getline(is, dummy); + microenv.thomas_i_jump = read_number_in_line_int(dummy); + + //thomas_j_jump + std::getline(is, dummy); + microenv.thomas_j_jump = read_number_in_line_int(dummy); + + //thomas_k_jump + std::getline(is, dummy); + microenv.thomas_k_jump = read_number_in_line_int(dummy); + + // thomas_constant1 + std::getline(is, dummy); + std::istringstream stream_thomas_constant1(dummy); + std::string key_thomas_constant1; + stream_thomas_constant1 >> key_thomas_constant1; + std::vector thomas_constant1; + while(stream_thomas_constant1 >> value) { + thomas_constant1.push_back(value); + } + microenv.thomas_constant1.clear(); + microenv.thomas_constant1 = thomas_constant1; + + // thomas_constant1a + std::getline(is, dummy); + std::istringstream stream_thomas_constant1a(dummy); + std::string key_thomas_constant1a; + stream_thomas_constant1a >> key_thomas_constant1a; + std::vector thomas_constant1a; + while(stream_thomas_constant1a >> value) { + thomas_constant1a.push_back(value); + } + microenv.thomas_constant1a.clear(); + microenv.thomas_constant1a = thomas_constant1a; + + // thomas_constant2 + std::getline(is, dummy); + std::istringstream stream_thomas_constant2(dummy); + std::string key_thomas_constant2; + stream_thomas_constant2 >> key_thomas_constant2; + std::vector thomas_constant2; + while(stream_thomas_constant2 >> value) { + thomas_constant2.push_back(value); + } + microenv.thomas_constant2.clear(); + microenv.thomas_constant2 = thomas_constant2; + + // thomas_constant3 + std::getline(is, dummy); + std::istringstream stream_thomas_constant3(dummy); + std::string key_thomas_constant3; + stream_thomas_constant3 >> key_thomas_constant3; + std::vector thomas_constant3; + while(stream_thomas_constant3 >> value) { + thomas_constant3.push_back(value); + } + microenv.thomas_constant3.clear(); + microenv.thomas_constant3 = thomas_constant3; + + // thomas_constant3a + std::getline(is, dummy); + std::istringstream stream_thomas_constant3a(dummy); + std::string key_thomas_constant3a; + stream_thomas_constant3a >> key_thomas_constant3a; + std::vector thomas_constant3a; + while(stream_thomas_constant3a >> value) { + thomas_constant3a.push_back(value); + } + microenv.thomas_constant3a.clear(); + microenv.thomas_constant3a = thomas_constant3a; + + //name + std::getline(is, dummy); + std::istringstream stream_name(dummy); + std::string key_name; + stream_name >> key_name >> value_2; + microenv.name = value_2; + + //supply_rates + std::getline(is, dummy); + n_lines = 0; + for( unsigned int n=0; n < microenv.supply_rates.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.supply_rates[n][j] = value; + j = j+1; + } + } + assert(n_lines == microenv.supply_rates.size()); + //read empty line + std::getline(is, dummy); + + //uptake_rates + std::getline(is, dummy); + n_lines = 0; + for( unsigned int n=0; n < microenv.uptake_rates.size() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + j = 0; + while(stream4 >> value) { + microenv.uptake_rates[n][j] = value; + j = j+1; + } + } + assert(n_lines == microenv.uptake_rates.size()); + //read empty line + std::getline(is, dummy); + + //number_of_densities + std::getline(is, dummy); + + double number_of_densities; + std::istringstream stream(dummy); + std::string key; + stream >> key; + double number; + while(stream >> number) { + number_of_densities = number; + } + + for (size_t i = 0; i < number_of_densities; ++i) { + //get name + std::string name; + std::getline(is, name); + microenv.density_names[i] = name; + + + // get unit + std::getline(is, dummy); + std::istringstream stream(dummy); + std::string key; + std::string unit; + stream >> key; + stream >> unit; + microenv.density_units[i] = unit; + + // get diffusion_coefficient + std::string spatial_unit; + std::string time_unit; + std::getline(is, dummy); + std::istringstream stream2(dummy); + double diffusion_coefficient; + stream2 >> key; + stream2 >> diffusion_coefficient; + stream2 >> spatial_unit; + stream2 >> time_unit; + microenv.diffusion_coefficients[i] = diffusion_coefficient; + microenv.spatial_units = spatial_unit; + microenv.time_units = time_unit; + + + //get decay rate + std::getline(is, dummy); + std::istringstream stream3(dummy); + double decay_rate; + stream3 >> key; + stream3 >> decay_rate; + microenv.decay_rates[i] = decay_rate; + } + + // skip empty line + std::getline(is, dummy); + + // skip header + std::getline(is, dummy); + + n_lines = 0; + //reassign all values for all the voxels. + for( unsigned int n=0; n < microenv.number_of_voxels() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + double value; + int j = 0; + while(stream4 >> value) { + microenv.density_vector(n)[j] = value; + j = j+1; + } + } + assert(n_lines == microenv.number_of_voxels()); + + std::getline(is, dummy); + std::getline(is, dummy); + + //reassign all values for all the voxels. + n_lines = 0; + for( unsigned int n=0; n < microenv.number_of_voxels() ; n++ ) + { + std::getline(is, dummy); n_lines++; + std::istringstream stream4(dummy); + double value; + int j = 0; + while(stream4 >> value) { + microenv.voxels(n).center[j] = value; + j = j+1; + } + } + assert(n_lines == microenv.number_of_voxels()); + + // skip empty line + std::getline(is, dummy); + std::getline(is, dummy); + + // dirichlet_value_vectors + n_lines = 0; + for (unsigned int n = 0; n < microenv.dirichlet_value_vectors.size(); n++) { + std::getline(is, dummy); + n_lines++; + std::istringstream stream(dummy); + size_t j = 0; + while (stream >> value) { + microenv.dirichlet_value_vectors[n][j] = value; + j++; + } + } + assert(n_lines == microenv.dirichlet_value_vectors.size()); + + // read empty line + std::getline(is, dummy); + // dirichlet_activation_vector + n_lines = 0; + std::getline(is, dummy); + bool dirichlet_activation_vector; + std::istringstream stream_dirichlet_activation_vector(dummy); + std::string key_dirichlet_activation_vector; + stream_dirichlet_activation_vector >> key_dirichlet_activation_vector; + for (int i = 0; i < microenv.dirichlet_activation_vector.size(); i++){ + stream_dirichlet_activation_vector >> value_2; + if(value_2 == "true"){ + dirichlet_activation_vector = true; + } + else if(value_2 == "false"){ + dirichlet_activation_vector = false; + } + microenv.dirichlet_activation_vector[i] = dirichlet_activation_vector; + } + + // read empty line + std::getline(is, dummy); + std::getline(is, dummy); + + // dirichlet_activation_vectors + n_lines = 0; + for (unsigned int n = 0; n < microenv.dirichlet_activation_vectors.size(); n++) { + std::getline(is, dummy); + n_lines++; + std::istringstream stream(dummy); + size_t j = 0; + while (stream >> value_2) { + if(value_2 == "true"){ + microenv.dirichlet_activation_vectors[n][j] = true; + } + else if(value_2 == "false"){ + microenv.dirichlet_activation_vectors[n][j] = false; + } + j++; + } + } + + // verifica che il numero di righe lette corrisponda + assert(n_lines == microenv.dirichlet_activation_vectors.size()); + + + return is; +} }; diff --git a/BioFVM/BioFVM_microenvironment.h b/BioFVM/BioFVM_microenvironment.h index 8b894e880..c8593051d 100644 --- a/BioFVM/BioFVM_microenvironment.h +++ b/BioFVM/BioFVM_microenvironment.h @@ -237,6 +237,8 @@ class Microenvironment void simulate_cell_sources_and_sinks( double dt ); void display_information( std::ostream& os ); + friend std::ostream& operator<<(std::ostream& os, Microenvironment& S); + friend std::istream& operator>>(std::istream& is, Microenvironment& microenv); void add_dirichlet_node( int voxel_index, std::vector& value ); void update_dirichlet_node( int voxel_index , std::vector& new_value ); diff --git a/Makefile b/Makefile index 76ad9b813..fc6c6816a 100644 --- a/Makefile +++ b/Makefile @@ -77,7 +77,7 @@ list-projects: @echo " asymmetric-division-sample" @echo "" @echo "Sample intracellular projects: template_BM ode-energy-sample physiboss-cell-lines-sample" - @echo " cancer-metabolism-sample physiboss-tutorial physiboss-tutorial-invasion" + @echo " physiboss-tnf-model cancer-metabolism-sample physiboss-tutorial physiboss-tutorial-invasion" @echo "" template: @@ -235,6 +235,16 @@ physiboss-cell-lines-sample: cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml cp ./sample_projects_intracellular/boolean/physiboss_cell_lines/config/* ./config/ +physiboss-tnf-model: + cp ./sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/* ./custom_modules/ + touch main.cpp && cp main.cpp main-backup.cpp + cp ./sample_projects_intracellular/boolean/spheroid_tnf_model/main-spheroid_TNF.cpp ./main.cpp + cp Makefile Makefile-backup + cp ./sample_projects_intracellular/boolean/spheroid_tnf_model/Makefile . + cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml + cp ./sample_projects_intracellular/boolean/spheroid_tnf_model/config/* ./config/ + cp -r ./sample_projects_intracellular/boolean/spheroid_tnf_model/scripts ./ + physiboss-tutorial: cp ./sample_projects_intracellular/boolean/tutorial/custom_modules/* ./custom_modules/ touch main.cpp && cp main.cpp main-backup.cpp diff --git a/addons/PhysiBoSS/src/maboss_intracellular.h b/addons/PhysiBoSS/src/maboss_intracellular.h index 06242a8c7..52d60ea35 100644 --- a/addons/PhysiBoSS/src/maboss_intracellular.h +++ b/addons/PhysiBoSS/src/maboss_intracellular.h @@ -117,6 +117,231 @@ class MaBoSSIntracellular : public PhysiCell::Intracellular { void print_current_nodes(){ this->maboss.print_nodes(); } + void save_current_nodes(std::ostream& out_stream){ + this->maboss.save_nodes(out_stream); + } + + void save_current_parameters_maboss(std::ostream& out_stream){ + this->maboss.save_current_parameters(out_stream); + } + void read_current_parameter_maboss(std::ifstream& in_stream){ + this->maboss.read_current_parameter(in_stream); + } + int get_number_of_nodes(){ + return this->maboss.get_nodes_number().size(); + } + + void reinit_maboss(std::string networkFile, std::string configFile){ + this->maboss.init_maboss(networkFile, configFile); + } + + void save_current_parameters(std::ostream& out_stream){ + + out_stream << "maboss_intracellular_parameters:" << std::endl; + out_stream << "time_step: " << this->time_step << std::endl; + out_stream << "discrete_time: " << (this->discrete_time ? "true" : "false") << std::endl; + out_stream << "time_tick: " << this->time_tick << std::endl; + out_stream << "scaling: " << this->scaling << std::endl; + out_stream << "time_stochasticity: " << this->time_stochasticity << std::endl; + out_stream << "inherit_state: " << (this->inherit_state ? "true" : "false") << std::endl; + out_stream << "start_time: " << this->start_time << std::endl; + out_stream << "next_physiboss_run: " << this->next_physiboss_run << std::endl; + + // Stampa inherit_nodes + out_stream << "inherit_nodes_size: " << inherit_nodes.size() << std::endl; + for (const auto& pair : inherit_nodes) { + out_stream << "Key: " << pair.first << " Value: " << (pair.second ? "true" : "false") << std::endl; + } + + // Stampa initial_values + out_stream << "initial_values_size: " << initial_values.size() << std::endl; + for (const auto& pair : initial_values) { + out_stream << "Key: " << pair.first << " Value: " << pair.second << std::endl; + } + + // Stampa mutations + out_stream << "mutations_size: " << mutations.size() << std::endl; + for (const auto& pair : mutations) { + out_stream << "Key: " << pair.first << " Value: " << pair.second << std::endl; + } + + // Stampa parameters + out_stream << "parameters_size: " << parameters.size() << std::endl; + for (const auto& pair : parameters) { + out_stream << "Key: " << pair.first << " Value: " << pair.second << std::endl; + } + + // Stampa indicesOfInputs + out_stream << "indicesOfInputs_size: " << indicesOfInputs.size() << std::endl; + for (size_t i = 0; i < indicesOfInputs.size(); ++i) { + out_stream << "Index_" << i << ": " << indicesOfInputs[i] << std::endl; + } + + // Stampa indicesOfOutputs + out_stream << "indicesOfOutputs_size: " << indicesOfOutputs.size() << std::endl; + for (size_t i = 0; i < indicesOfOutputs.size(); ++i) { + out_stream << "Index_" << i << ": " << indicesOfOutputs[i] << std::endl; + } + out_stream << "listOfInputs_size: " << listOfInputs.size() << std::endl; + for (const auto& pair : listOfInputs) { + out_stream << pair.first << std::endl; + pair.second.save_maboss_input(out_stream); + } + + out_stream << "listOfoUTPUTS_size: " << listOfOutputs.size() << std::endl; + for (const auto& pair : listOfOutputs) { + out_stream << pair.first << std::endl; + pair.second.save_maboss_output(out_stream); + } + + out_stream << std::endl; + + } + + void read_current_parameter(std::ifstream& in_stream){ + std::string dummy; + + //skip header + std::getline(in_stream, dummy); + + //time_step + std::getline(in_stream, dummy); + this->time_step = read_number_in_line(dummy); + + //discrete_time + std::getline(in_stream, dummy); + this->discrete_time = read_number_in_line_bool(dummy); + + //time_tick + std::getline(in_stream, dummy); + this->time_tick = read_number_in_line(dummy); + + //scaling + std::getline(in_stream, dummy); + this->scaling = read_number_in_line(dummy); + + //time_stochasticity + std::getline(in_stream, dummy); + this->time_stochasticity = read_number_in_line(dummy); + + //inherit_state + std::getline(in_stream, dummy); + this->inherit_state = read_number_in_line_bool(dummy); + + //start_time + std::getline(in_stream, dummy); + this->start_time = read_number_in_line(dummy); + + //next_physiboss_run + std::getline(in_stream, dummy); + this->next_physiboss_run = read_number_in_line(dummy); + + //inherit_nodes + int inherit_nodes_size; + std::getline(in_stream, dummy); + inherit_nodes_size = read_number_in_line_int(dummy); + for (int i = 0; i < inherit_nodes_size; ++i) { + std::getline(in_stream, dummy); + std::istringstream stream(dummy); + std::string key, value; + stream >> key >> value; + + if (value == "true") { + inherit_nodes[key] = true; + } else if (value == "false") { + inherit_nodes[key] = false; + } else { + throw std::runtime_error("Invalid boolean value in line: " + dummy); + } + } + // initial_values + int initial_values_size; + std::getline(in_stream, dummy); + initial_values_size = read_number_in_line_int(dummy); + for (int i = 0; i < initial_values_size; ++i) { + std::getline(in_stream, dummy); + std::istringstream stream(dummy); + std::string key; + double value; + stream >> key >> value; + initial_values[key] = value; + } + + // mutations + int mutations_size; + std::getline(in_stream, dummy); + mutations_size = read_number_in_line_int(dummy); + for (int i = 0; i < mutations_size; ++i) { + std::getline(in_stream, dummy); + std::istringstream stream(dummy); + std::string key; + double value; + stream >> key >> value; + mutations[key] = value; + } + + // parameters + int parameters_size; + std::getline(in_stream, dummy); + parameters_size = read_number_in_line_int(dummy); + for (int i = 0; i < parameters_size; ++i) { + std::getline(in_stream, dummy); + std::istringstream stream(dummy); + std::string key; + double value; + stream >> key >> value; + parameters[key] = value; + } + + // indicesOfInputs + int indicesOfInputs_size; + std::getline(in_stream, dummy); + indicesOfInputs_size = read_number_in_line_int(dummy); + indicesOfInputs.clear(); + for (int i = 0; i < indicesOfInputs_size; ++i) { + std::getline(in_stream, dummy); + + int new_value = read_number_in_line_int(dummy); + + indicesOfInputs.push_back(new_value); + } + + // indicesOfOutputs + int indicesOfOutputs_size; + std::getline(in_stream, dummy); + indicesOfOutputs_size = read_number_in_line_int(dummy); + indicesOfOutputs.clear(); // Assicuriamoci che sia vuoto prima di riempirlo + for (int i = 0; i < indicesOfOutputs_size; ++i) { + std::getline(in_stream, dummy); + indicesOfOutputs.push_back(read_number_in_line_int(dummy)); + } + + int listOfInputs_size; + std::getline(in_stream, dummy); + listOfInputs_size = read_number_in_line_int(dummy); + listOfInputs.clear(); + for (int i = 0; i < listOfInputs_size; ++i) { + std::getline(in_stream, dummy); + std::string key = dummy; + MaBoSSInput input; + input.read_maboss_input(in_stream); + listOfInputs[key] = input; + } + + int listOfOutputs_size; + std::getline(in_stream, dummy); + listOfOutputs_size = read_number_in_line_int(dummy); + listOfOutputs.clear(); + for (int i = 0; i < listOfOutputs_size; ++i) { + std::getline(in_stream, dummy); + std::string key = dummy; + MaBoSSOutput output; + output.read_maboss_output(in_stream); + listOfOutputs[key] = output; + } + // skip empty line + std::getline(in_stream, dummy); + } void display(std::ostream& os); diff --git a/addons/PhysiBoSS/src/maboss_network.cpp b/addons/PhysiBoSS/src/maboss_network.cpp index ebcee8be3..b4561381e 100644 --- a/addons/PhysiBoSS/src/maboss_network.cpp +++ b/addons/PhysiBoSS/src/maboss_network.cpp @@ -164,4 +164,21 @@ void MaBoSSNetwork::print_nodes() i++; } std::cout << std::endl; +} + +void MaBoSSNetwork::save_nodes(std::ostream& out_stream) +{ + int i = 0; + //out_stream << this->get_state(); + + std::vector nodes = this->network->getNodes(); + for ( auto node: nodes ) + { + out_stream << node->getLabel() << "=" << state.getNodeState(node) << " "; + i++; + } +} + +std::vector MaBoSSNetwork::get_nodes_number(){ + return this->network->getNodes(); } \ No newline at end of file diff --git a/addons/PhysiBoSS/src/maboss_network.h b/addons/PhysiBoSS/src/maboss_network.h index 5aace16d3..a2d0aed6b 100644 --- a/addons/PhysiBoSS/src/maboss_network.h +++ b/addons/PhysiBoSS/src/maboss_network.h @@ -5,6 +5,7 @@ #include "BooleanNetwork.h" #include "RunConfig.h" #include "../../../core/PhysiCell_utilities.h" +#include "../../../core/PhysiCell_phenotype.h" /** * \class MaBoSSNetwork @@ -111,6 +112,7 @@ class MaBoSSNetwork bool has_node( std::string name ); void set_node_value(std::string name, bool value); bool get_node_value(std::string name); + std::vector get_nodes_number(); std::string get_state(); /** \brief Return update time value */ @@ -141,6 +143,89 @@ class MaBoSSNetwork */ void print_nodes(); + /** + * \brief Save current state of all the nodes of the network + * \param node_values Boolean vector mapping a boolean network + * \param out_stream Output file used to save data + */ + void save_nodes(std::ostream& out_stream); + +void save_current_parameters(std::ostream& out_stream){ + + out_stream << "maboss_network_parameters:" << std::endl; + out_stream << "update_time_step: " << this->update_time_step << std::endl; + out_stream << "time_to_update: " << this->time_to_update << std::endl; + out_stream << "scaling: " << this->scaling << std::endl; + out_stream << "time_stochasticity: " << this->time_stochasticity << std::endl; + + out_stream << "initial_values_size: " << initial_values.size() << std::endl; + for (const auto& pair : initial_values) { + out_stream << "Key: " << pair.first << " Value: " << pair.second << std::endl; + } + + out_stream << "mutations_size: " << mutations.size() << std::endl; + for (const auto& pair : mutations) { + out_stream << "Key: " << pair.first << " Value: " << pair.second << std::endl; + } + + + //out_stream << "seed: " << this->seed << std::endl; + out_stream << std::endl; + } + + void read_current_parameter(std::ifstream& in_stream){ + std::string dummy; + + //skip header + std::getline(in_stream, dummy); + + // update_time_step + std::getline(in_stream, dummy); + this->update_time_step = read_number_in_line(dummy); + + //time_to_update + std::getline(in_stream, dummy); + this->time_to_update = read_number_in_line(dummy); + + //scaling + std::getline(in_stream, dummy); + this->scaling = read_number_in_line(dummy); + + //time_stochasticity + std::getline(in_stream, dummy); + this->time_stochasticity = read_number_in_line(dummy); + + // initial_values + int initial_values_size; + std::getline(in_stream, dummy); + initial_values_size = read_number_in_line_int(dummy); + for (int i = 0; i < initial_values_size; ++i) { + std::getline(in_stream, dummy); + std::istringstream stream(dummy); + std::string key; + double value; + stream >> key >> value; + initial_values[key] = value; + } + + // mutations + int mutations_size; + std::getline(in_stream, dummy); + mutations_size = read_number_in_line_int(dummy); + for (int i = 0; i < mutations_size; ++i) { + std::getline(in_stream, dummy); + std::istringstream stream(dummy); + std::string key; + double value; + stream >> key >> value; + mutations[key] = value; + } + + // skip empty line + std::getline(in_stream, dummy); + + } + void set_state(NetworkState _state) { state = NetworkState(_state.getState()); } NetworkState get_maboss_state() { return state;} void inherit_state(NetworkState mother, bool inherit_state, std::map& inherit_nodes) { diff --git a/addons/PhysiBoSS/src/utils.h b/addons/PhysiBoSS/src/utils.h index 332eacb84..9ba4e6a45 100644 --- a/addons/PhysiBoSS/src/utils.h +++ b/addons/PhysiBoSS/src/utils.h @@ -18,6 +18,75 @@ class MaBoSSInput int smoothing; double smoothed_value; bool use_for_dead; + + MaBoSSInput() + : physicell_name(""), + type(NODE), + intracellular_name(""), + intracellular_parameter(""), + action(""), + threshold(0.0), + inact_threshold(0.0), + scaling(1.0), + smoothing(1), + smoothed_value(0.0), + use_for_dead(false) {} + + void save_maboss_input(std::ostream& out_stream) const { + out_stream << "MaBoSSInput: " << std::endl; + out_stream << this->physicell_name << std::endl; + out_stream << "type: " << this->type << std::endl; + out_stream << this->intracellular_name << std::endl; + out_stream << this->intracellular_parameter << std::endl; + out_stream << this->action << std::endl; + out_stream << "threshold: " << this->threshold << std::endl; + out_stream << "inact_threshold: " << this->inact_threshold << std::endl; + out_stream << "scaling: " << this->scaling << std::endl; + out_stream << "smoothing: " << this->smoothing << std::endl; + out_stream << "smoothed_value: " << this->smoothed_value << std::endl; + out_stream << "use_for_dead: " << (this->use_for_dead ? "true" : "false") << std::endl; + } + + void read_maboss_input(std::ifstream& in_stream) { + std::string dummy; + std::getline(in_stream, dummy); + + std::getline(in_stream, dummy); + physicell_name = dummy; + + std::getline(in_stream, dummy); + type = read_number_in_line_int(dummy); + + std::getline(in_stream, dummy); + intracellular_name = dummy; + + std::getline(in_stream, dummy); + intracellular_parameter = dummy; + + std::getline(in_stream, dummy); + action = dummy; + + std::getline(in_stream, dummy); + threshold = read_number_in_line(dummy); + + std::getline(in_stream, dummy); + inact_threshold = read_number_in_line(dummy); + + std::getline(in_stream, dummy); + //std::cout << "dummy: " << dummy << std::endl; + //std::cout << "scaling pre: " << scaling << std::endl; + scaling = read_number_in_line(dummy); + //std::cout << "scaling post: " << scaling << std::endl; + + std::getline(in_stream, dummy); + smoothing = read_number_in_line_int(dummy); + + std::getline(in_stream, dummy); + smoothed_value = read_number_in_line(dummy); + + std::getline(in_stream, dummy); + use_for_dead = read_number_in_line_bool(dummy); + } MaBoSSInput(std::string physicell_name, std::string intracellular_name, std::string action, double threshold, double inact_threshold, int smoothing, bool use_for_dead) : physicell_name(physicell_name), intracellular_name(intracellular_name), action(action), threshold(threshold), inact_threshold(inact_threshold), smoothing(smoothing), use_for_dead(use_for_dead){ type = NODE; @@ -88,6 +157,67 @@ class MaBoSSOutput int steepness; bool use_for_dead; + MaBoSSOutput() + : physicell_name(""), + intracellular_name(""), + action(""), + value(0.0), + base_value(0.0), + smoothing(1), + probability(0.5), + initialized(false), + steepness(1), + use_for_dead(false) {} + + void save_maboss_output(std::ostream& out_stream) const { + out_stream << "MaBoSSOutput: " << std::endl; + out_stream << this->physicell_name << std::endl; + out_stream << this->intracellular_name << std::endl; + out_stream << this->action << std::endl; + out_stream << "value: " << this->value << std::endl; + out_stream << "base_value: " << this->base_value << std::endl; + out_stream << "smoothing: " << this->smoothing << std::endl; + out_stream << "probability: " << this->probability << std::endl; + out_stream << "initialized: " << (this->initialized ? "true" : "false") << std::endl; + out_stream << "steepness: " << this->steepness << std::endl; + out_stream << "use_for_dead: " << (this->use_for_dead ? "true" : "false") << std::endl; + } + + void read_maboss_output(std::ifstream& in_stream) { + std::string dummy; + std::getline(in_stream, dummy); // Legge l'intestazione (da ignorare) + + std::getline(in_stream, dummy); + physicell_name = dummy; + + std::getline(in_stream, dummy); + intracellular_name = dummy; + + std::getline(in_stream, dummy); + action = dummy; + + std::getline(in_stream, dummy); + value = read_number_in_line(dummy); + + std::getline(in_stream, dummy); + base_value = read_number_in_line(dummy); + + std::getline(in_stream, dummy); + smoothing = read_number_in_line_int(dummy); + + std::getline(in_stream, dummy); + probability = read_number_in_line(dummy); + + std::getline(in_stream, dummy); + initialized = read_number_in_line_bool(dummy); + + std::getline(in_stream, dummy); + steepness = read_number_in_line_int(dummy); + + std::getline(in_stream, dummy); + use_for_dead = read_number_in_line_bool(dummy); + } + MaBoSSOutput(std::string physicell_name, std::string intracellular_name, std::string action, double value, double base_value, int smoothing, int steepness, bool use_for_dead) diff --git a/addons/start_and_stop/start_and_stop.cpp b/addons/start_and_stop/start_and_stop.cpp new file mode 100644 index 000000000..ce3e660cc --- /dev/null +++ b/addons/start_and_stop/start_and_stop.cpp @@ -0,0 +1,501 @@ +#include "./start_and_stop.h" +#include +#include + +using namespace std; + +// reset microenv +void reset_microenv(){ + + string filename = "start_and_stop_saving_files/microenv_data.txt"; + + // Open the file for input + ifstream input_file(filename); + if (!input_file.is_open()) { + cerr << "Error: Unable to open file " << filename << endl; + return; + } + + // Call the operator>> function to read data from the file into the Microenvironment instance + input_file >> microenvironment; + + setup_signal_behavior_dictionaries(); + + // Close the file + input_file.close(); + return; +} + +// Function to save cell data +void save_cell_microenv_data(Cell_Container* cell_container) +{ + + // decide if restart from a previous point or from the beginning + string filename_cells = "start_and_stop_saving_files/cell_data.txt"; + string filename_bool = "start_and_stop_saving_files/bool_data.txt"; + string filename_microenv = "start_and_stop_saving_files/microenv_data.txt"; + string cell_pos_filename = "start_and_stop_saving_files/cells.csv"; + string global_param_filename = "start_and_stop_saving_files/global_param.txt"; + string random_counters_filename = "start_and_stop_saving_files/random_counters.txt"; + + // Open the files in write mode + ofstream file_cells(filename_cells); + ofstream file_bool(filename_bool); + ofstream file_microenv(filename_microenv); + ofstream cell_pos_file(cell_pos_filename); + ofstream global_param_file(global_param_filename); + ofstream random_counters_file(random_counters_filename); + + // Check if the files are successfully opened + if (!file_cells.is_open()) + { + // Handle the case where the file cannot be opened + cerr << "Error opening file " << filename_cells << endl; + return; + } + if (!file_bool.is_open()) + { + // Handle the case where the file cannot be opened + cerr << "Error opening file " << filename_bool << endl; + return; + } + if (!file_microenv.is_open()) + { + // Handle the case where the file cannot be opened + cerr << "Error opening file " << filename_microenv << endl; + return; + } + if (!cell_pos_file.is_open()) + { + // Handle the case where the file cannot be opened + cerr << "Error opening file " << cell_pos_filename << endl; + return; + } + if (!global_param_file.is_open()) + { + // Handle the case where the file cannot be opened + cerr << "Error opening file " << global_param_filename << endl; + return; + } + if (!random_counters_file.is_open()) + { + // Handle the case where the file cannot be opened + cerr << "Error opening file " << random_counters_filename << endl; + return; + } + + // display microenvironment + + file_microenv << microenvironment; + + // add header to pos file + cell_pos_file << "x,y,z,type" << endl; + + // Iterate through all cells + for (int i = 0; i < (*all_cells).size(); i++) + { + // Access the current cell + Cell *pCell = (*all_cells)[i]; + + // Write the cell data to the file using the insertion operator + file_cells << *pCell << endl; + + if (pCell->phenotype.intracellular){ + //save parameters maboss_intracellular + pCell->phenotype.intracellular->save_current_parameters(file_cells); + + //save parameters maboss network + pCell->phenotype.intracellular->save_current_parameters_maboss(file_cells); + } + + // current_voxel_index + file_cells << "current_voxel_index: " << pCell->get_current_voxel_index() << std::endl; + + // nearest_density_vector + file_cells << "nearest_density_vector: " << pCell->nearest_density_vector() << std::endl; + + + // empty line + file_cells << std::endl; + + // Write the positions in the cell file + cell_pos_file << pCell->position[0] << ","; + cell_pos_file << pCell->position[1] << ","; + cell_pos_file << pCell->position[2] << ","; + cell_pos_file << pCell->type_name; + + + if (pCell->phenotype.intracellular){ + // Write the output for Boolean network + // iterate among the nodes and write the name of all the active nodes for each cell, a row = a cell. + pCell->phenotype.intracellular->save_current_nodes(file_bool); + + + + // Write newline except for the last cell + if (i != (*all_cells).size() - 1) { + cell_pos_file << "\n"; + file_bool << "\n"; + } + } + + + } + // Save global parameters and cell container + global_param_file << "current_time: " << PhysiCell_globals.current_time << endl; + global_param_file << "full_output_index: " <phenotype.intracellular){ + global_param_file << "intracellular_dt: " << intracellular_dt << std::endl; + } + cell_container->save_data(global_param_file); + + //save random counters + save_counters(random_counters_file); + + //close all the files + file_cells.close(); + file_bool.close(); + file_microenv.close(); + global_param_file.close(); + random_counters_file.close(); + + return; +} +using namespace std; + + +void reset_cell(double last_cell_cycle_time) +{ + + // Specify the path to the input files + std::string input_file_cells_path = "start_and_stop_saving_files/cell_data.txt"; + std::string input_file_bool_path = "start_and_stop_saving_files/bool_data.txt"; + + // Open the files in input mode using std::ifstream + std::ifstream input_file_cells(input_file_cells_path); + std::ifstream input_file_bool(input_file_bool_path); + + + + // Check if the files was opened successfully + if (!input_file_cells.is_open()) + { + std::cerr << "Error opening the file: " << input_file_cells_path << std::endl; + return; + } + if (!input_file_bool.is_open()) + { + std::cerr << "Error opening the file: " << input_file_bool_path << std::endl; + return; + } + + // Iterate through all cells + for (int i = 0; i < (*all_cells).size(); i++) + { + // Use the extraction operator >> to read data from the file line into the cell + + // Access the current cell + Cell *pCell = (*all_cells)[i]; + //read header + string header; + getline(input_file_cells, header); + // ID: + getline(input_file_cells, header); + pCell->ID = read_number_in_line_int(header); + // index: + getline(input_file_cells, header); + pCell->index = read_number_in_line_int(header); + // type: + getline(input_file_cells, header); + pCell->type = read_number_in_line_int(header); + + // is_active: + std::getline(input_file_cells, header); + pCell->is_active = read_number_in_line_bool(header); + // type_name + getline(input_file_cells >> ws, pCell->type_name); + // is_out_of_domain + std::getline(input_file_cells, header); + pCell->is_out_of_domain = read_number_in_line_bool(header); + + // is_movable + std::getline(input_file_cells, header); + pCell->is_movable = read_number_in_line_bool(header); + + // state + input_file_cells >> pCell->state; + + + //parameters + input_file_cells >> pCell->parameters; + + + + string line; + + // flagged_for_division: + std::getline(input_file_cells, line); + pCell->phenotype.flagged_for_division = read_number_in_line_bool(line); + + + // flagged_for_removal: + std::getline(input_file_cells, line); + pCell->phenotype.flagged_for_removal = read_number_in_line_bool(line); + + + //death + input_file_cells >> pCell->phenotype.death; + + if(pCell->phenotype.death.dead == true){ + pCell->phenotype.cycle.sync_to_cycle_model(pCell->phenotype.death.current_model()); + } + //cycle + input_file_cells >> pCell->phenotype.cycle; + if( pCell->phenotype.cycle.current_phase().entry_function ) + { + pCell->phenotype.cycle.current_phase().entry_function( pCell, pCell->phenotype, PhysiCell_globals.current_time - last_cell_cycle_time); + } + + //volume + input_file_cells >> pCell->phenotype.volume; + pCell->set_total_volume(pCell->phenotype.volume.total); + //geometry + input_file_cells >> pCell->phenotype.geometry; + //mechanics + input_file_cells >> pCell->phenotype.mechanics; + //motility + input_file_cells >> pCell->phenotype.motility; + + //secretion + input_file_cells >> pCell->phenotype.secretion; + + //molecular + input_file_cells >> pCell->phenotype.molecular; + + //Cell_Interactions + input_file_cells >> pCell->phenotype.cell_interactions; + + //Cell_Transformations + input_file_cells >> pCell->phenotype.cell_transformations; + + //custom_data + input_file_cells >> pCell->custom_data; + + getline(input_file_cells, header); + getline(input_file_cells, header); + + //intracellular + if (pCell->phenotype.intracellular){ + pCell->phenotype.intracellular->read_current_parameter(input_file_cells); + + pCell->phenotype.intracellular->read_current_parameter_maboss(input_file_cells); + + std::string current_bnd_filename; + std::string current_cfg_filename; + + // load and parse settings file(s) + pugi::xml_document physicell_xml; + pugi::xml_node node; + pugi::xml_parse_result result = physicell_xml.load_file( "./config/PhysiCell_settings.xml" ); + + if( result.status != pugi::xml_parse_status::status_ok ) + { + std::cout << "Error loading " << "./config/PhysiCell_settings.xml" << "!" << std::endl; + } + + // navigate to the desired value + pugi::xml_node physicell_settings = physicell_xml.child("PhysiCell_settings"); + pugi::xml_node cell_defs = physicell_settings.child("cell_definitions"); + + for (pugi::xml_node cell_def = cell_defs.child("cell_definition"); cell_def; cell_def = cell_def.next_sibling("cell_definition")) + { + if (std::string(cell_def.attribute("name").value()) == pCell->type_name) + { + pugi::xml_node phenotype = cell_def.child("phenotype"); + + pugi::xml_node intracellular = phenotype.child("intracellular"); + + if(intracellular && std::string(intracellular.attribute("type").value()) == "maboss") + { + pugi::xml_node bnd_filename = intracellular.child("bnd_filename"); + pugi::xml_node cfg_filename = intracellular.child("cfg_filename"); + + current_bnd_filename = bnd_filename.text().get(); + current_cfg_filename = cfg_filename.text().get(); + + } + } + } + // reinit maboss + //pCell->phenotype.intracellular->reinit_maboss(current_bnd_filename, current_cfg_filename); + } + + // current_voxel_index + std::getline(input_file_cells, header); + + // nearest_density_vector + std::getline(input_file_cells, header); + + // empty line + std::getline(input_file_cells, header); + + if (pCell->phenotype.intracellular){ + // Boolean network reinitialization + string dummy; + getline(input_file_bool, dummy); + + // Use a stringstream to extract tokens from the line + istringstream iss(dummy); + + // Loop to extract node names and values from the line + string token; + unsigned int n_line = 0; + while (getline(iss, token, ' ')) { // Read the string until the ';' character + // Find the position of the '=' character + size_t pos = token.find('='); + if (pos != string::npos) { // If the '=' character is found + // Extract the node name and its value from the substring + string node_name = token.substr(0, pos); + string node_value_str = token.substr(pos + 1); + bool node_value = (node_value_str == "1"); + + pCell->phenotype.intracellular->set_boolean_variable_value(node_name, node_value); + n_line++; + + } + + + } + // Check if the number of nodes is the same + + assert(n_line == pCell->phenotype.intracellular->get_number_of_nodes()); + //pCell->phenotype.intracellular->set_state_int(); + //pCell->phenotype.intracellular->update(pCell, pCell->phenotype, phenotype_dt); + } + + } + input_file_cells.close(); + input_file_bool.close(); + return; + + } + + void reset_global_parameters(Cell_Container* cell_container) + { + string global_param_filename = "start_and_stop_saving_files/global_param.txt"; + + ifstream global_param_file(global_param_filename); + + if (!global_param_file.is_open()) + { + // Handle the case where the file cannot be opened + cerr << "Error opening file " << global_param_filename << endl; + return; + } + + // current_time: + string line_current_time; + getline(global_param_file, line_current_time); + PhysiCell_globals.current_time = read_number_in_line(line_current_time); + + // full_output_index: + string line_full_output_index; + getline(global_param_file, line_full_output_index); + PhysiCell_globals.full_output_index = read_number_in_line_int(line_full_output_index); + + // intracellular_output_index: + string line_intracellular_output_index; + getline(global_param_file, line_intracellular_output_index); + PhysiCell_globals.intracellular_output_index = read_number_in_line_int(line_intracellular_output_index); + + // next_full_save_time: + string line_next_full_save_time; + getline(global_param_file, line_next_full_save_time); + PhysiCell_globals.next_full_save_time = read_number_in_line(line_next_full_save_time); + + // next_intracellular_save_time: + string line_next_intracellular_save_time; + getline(global_param_file, line_next_intracellular_save_time); + PhysiCell_globals.next_intracellular_save_time = read_number_in_line(line_next_intracellular_save_time); + + // next_SVG_save_time: + string line_next_SVG_save_time; + getline(global_param_file, line_next_SVG_save_time); + PhysiCell_globals.next_SVG_save_time = read_number_in_line(line_next_SVG_save_time); + + // SVG_output_index: + string line_SVG_output_index; + getline(global_param_file, line_SVG_output_index); + PhysiCell_globals.SVG_output_index = read_number_in_line_int(line_SVG_output_index); + + // diffusion_dt: + string line_diffusion_dt; + getline(global_param_file, line_diffusion_dt); + diffusion_dt = read_number_in_line(line_diffusion_dt); + + // diffusion_dt: + string line_mechanics_dt; + getline(global_param_file, line_mechanics_dt); + mechanics_dt = read_number_in_line(line_mechanics_dt); + + // phenotype_dt: + string line_phenotype_dt; + getline(global_param_file, line_phenotype_dt); + phenotype_dt = read_number_in_line(line_phenotype_dt); + + // intracellular_dt: + if (((*all_cells)[0])->phenotype.intracellular){ + string line_intracellular_dt; + getline(global_param_file, line_intracellular_dt); + intracellular_dt = read_number_in_line(line_intracellular_dt); + } + + //reset cell_container + cell_container->reset_data(global_param_file); + + // Close the file after finishing reading + global_param_file.close(); + + return; + } + + void reset_randomness() + { + + string random_counters_filename = "start_and_stop_saving_files/random_counters.txt"; + + ifstream random_counters_file(random_counters_filename); + + if (!random_counters_file.is_open()) + { + // Handle the case where the file cannot be opened + cerr << "Error opening file " << random_counters_filename << endl; + return; + } + + // counter_double_random: + string line_counter_double_random; + getline(random_counters_file, line_counter_double_random); + set_counter_double_random(read_number_in_line(line_counter_double_random)); + + // counter_int_random: + string line_counter_int_random; + getline(random_counters_file, line_counter_int_random); + set_counter_int_random(read_number_in_line(line_counter_int_random)); + + // counter_normal_random: + string line_counter_normal_random; + getline(random_counters_file, line_counter_normal_random); + set_counter_normal_random(read_number_in_line(line_counter_normal_random)); + + random_counters_file.close(); + + } diff --git a/addons/start_and_stop/start_and_stop.h b/addons/start_and_stop/start_and_stop.h new file mode 100644 index 000000000..c9bee9c07 --- /dev/null +++ b/addons/start_and_stop/start_and_stop.h @@ -0,0 +1,17 @@ + +#include "../../core/PhysiCell.h" +#include "../../modules/PhysiCell_standard_modules.h" + +using namespace BioFVM; +using namespace PhysiCell; + +void reset_microenv(); + +// Function to save cell data +void save_cell_microenv_data(Cell_Container* cell_container); + +void reset_cell(double last_cell_cycle_time); + +void reset_global_parameters(Cell_Container* cell_container); + +void reset_randomness(); \ No newline at end of file diff --git a/beta/start_stop_tests.py b/beta/start_stop_tests.py new file mode 100644 index 000000000..954576196 --- /dev/null +++ b/beta/start_stop_tests.py @@ -0,0 +1,61 @@ +import subprocess +import xml.etree.ElementTree as ET +import os +import time + +user_proj = ["biorobots-sample", "physiboss-tnf-model"] + +model_execs = ["./biorobots", "./spheroid_TNF_model"] + + +# Using dummy max_time values of 99 for many projects; 61 for the more time-consuming cancer_immune_3D. +# Users can change them as they wish. +max_times = [20, 1440] + +for (uproj, myexec, max_time) in zip(user_proj, model_execs, max_times): + times = [max_time/2, max_time] + for current_time in times: + print("\n\n------------ ",uproj,myexec, " ----------------------------------") + if current_time == max_time: + # execute everything but do not reload the project + tree = ET.parse('config/PhysiCell_settings.xml') + root = tree.getroot() + root.find(".//max_time").text = str(current_time) + root.find(".//omp_num_threads").text = "1" + new_output_dir = "output_" + uproj + root.find(".//save//folder").text = new_output_dir + root.find(".//user_parameters//start_stop").text = "true" + + # set the start and stop to true for the second execution + tree.write('./config/PhysiCell_settings.xml') + try: + os.makedirs(new_output_dir) + except: + pass + time.sleep(1) + + subprocess.run([myexec]) + else: + subprocess.run(["make","reset"]) + subprocess.run(["make","clean"]) + # make load PROJ=template + subprocess.run(["make",uproj]) + subprocess.run(["make"]) + # update max_time and omp_num_threads (=1) + tree = ET.parse('config/PhysiCell_settings.xml') + root = tree.getroot() + root.find(".//max_time").text = str(current_time) + root.find(".//omp_num_threads").text = "1" + new_output_dir = "output_" + uproj + root.find(".//save//folder").text = new_output_dir + root.find(".//user_parameters//start_stop").text = "false" + + # set the start and stop to true for the second execution + tree.write('./config/PhysiCell_settings.xml') + try: + os.makedirs(new_output_dir) + except: + pass + time.sleep(1) + + subprocess.run([myexec]) \ No newline at end of file diff --git a/beta/test_build_samples.sh b/beta/test_build_samples.sh old mode 100644 new mode 100755 diff --git a/config/PhysiCell_settings.xml b/config/PhysiCell_settings.xml index e3cd03b4a..a32e7c30d 100644 --- a/config/PhysiCell_settings.xml +++ b/config/PhysiCell_settings.xml @@ -1,271 +1,225 @@ + + -500 + 500 + -500 + 500 + -10 + 10 + 20 + 20 + 20 + true + + + + 180 + min + micron + 0.02 + 0.1 + 2 + + + + 8 + + + + output + + 30 + true + + + 30 + true + + + false + + + + + + + 100000.0 + .1 + + 38.0 + 38.0 + + + + + 1200.0 + .0275 + + 0.0 + 0.0 + - - -500 - 500 - -500 - 500 - -10 - 10 - 20 - 20 - 20 - true - + + false + true + + + ./config/initial.mat + + + + ./config/dirichlet.mat + + + + + + + + + + + + 0.0011 + + + + + + 5.31667e-05 + + 516 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + - - 14400 - min - micron - 0.01 - 0.1 - 6 - + + 0.0 + + 0 + 86400 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + + + 2494 + 0.75 + 540 + + 0.05 + 0.0045 + 0.0055 + + 0 + 0 + + 2.0 + - - 6 - + + + 0 + 38 + 10 + 0 + + + 0 + 0 + 0 + 0 + + + - - output - - 60 - true - - - 60 - true - - - false - - + + ./config/TNF_nodes.bnd + ./config/TNF_conf.cfg + + 10.0 + 0.01 + + + + + + 0.243 + 0.128 + 0.293 + + 0.50 + 1.0 + 0 - - false - true - + + 0 + 0 + 0 + 0 + + + + 0.0035 + 0.0035 + 1.0 + + + - - - - 100000.0 - 0.1 - - 38 - 0 - - 38 - 10 - 10 - 38 - 0 - 0 - - - - - 10 - .1 - - 0 - 0 - - 0 - 0 - 0 - 0 - 0 - 0 - - - - - 10 - 0.1 - - 0 - 0 - - 0 - 0 - 0 - 0 - 0 - 0 - - - - true - true - - ./config/initial.mat - - - ./config/dirichlet.mat - - - + + 0 + false + ./start_and_stop_saving_files/initial.tsv + 100 + - - - - - - 0.000 - - - - - 5.31667e-05 - - 516 - - - 0.05 - 0 - 1.66667e-02 - 5.83333e-03 - 0 - 2.0 - - - - 2.80E-03 - - 0 - 86400 - - - 1.11667e-2 - 8.33333e-4 - 5.33333e-5 - 2.16667e-3 - 0 - 2.0 - - - - - 2494 - 0.75 - 540 - 0.05 - 0.0045 - 0.0055 - 0 - 0 - 2.0 - - - 0.4 - 10.0 - 1.25 - - 1 - - - 1.8 - 15.12 - - 4.0 - 10.0 - 0.01 - 0.0 - 0.0 - - - 1 - 1 - .5 - - false - true - - false - oxygen - 1 - - - false - false - - 0.0 - 0.0 - 0.0 - - - - - - - 0 - 1 - 10 - 0 - - - 0.0 - 1 - 0.0 - 0.0 - - - 0.0 - 1 - 0.0 - 0.0 - - - - 0 - - 0 - 0 - - - 0 - 0 - - 1 - - 0 - 0 - - - - - 0 - 0 - - - - - 1.0 - - - + 120 + + false + false + false - - - ./config - cells.csv - - + true + 150 + 10 + 0.005 + 999999 + + + + 0.0 + 0.05 + 0.95 - - - - ./config - cell_rules_v2.csv - - - - + + 0.0 + 0.05 + 0.95 + + 0.0 + 0.05 + 0.95 + - - 0 - 0 - + + + + \ No newline at end of file diff --git a/config/initial.tsv b/config/initial.tsv new file mode 100644 index 000000000..02e46f221 --- /dev/null +++ b/config/initial.tsv @@ -0,0 +1,126 @@ +x y z +-93.7351 -16.7883 0.0 +-95.7338 2.90696 0.0 +-95.9399 18.5505 0.0 +-79.8482 -55.7662 0.0 +-82.7174 -41.7014 0.0 +-80.7886 -22946 0.0 +-79.0268 -6.05337 0.0 +-82.0067 10.4321 0.0 +-80.6049 25.4132 0.0 +-82.8366 44.0516 0.0 +-78.9526 61124 0.0 +-64102 -65.5156 0.0 +-64.2225 -46.9066 0.0 +-68468 -33.3416 0.0 +-68.4907 -14.3586 0.0 +-67762 -0.0687358 0.0 +-64.3596 17.0034 0.0 +-65.1953 36802 0.0 +-65.7418 51.8249 0.0 +-64.8725 70.5213 0.0 +-49.7993 -75.1665 0.0 +-50.1662 -59.6431 0.0 +-51.6159 -42.4253 0.0 +-52.2097 -24.4339 0.0 +-49.6998 -5.49337 0.0 +-52.8067 10.2768 0.0 +-53.0525 25.5959 0.0 +-50.1145 42.8981 0.0 +-50.8895 63.7084 0.0 +-51.5362 77.8187 0.0 +-35.1206 -84.1124 0.0 +-39.0413 -67.4176 0.0 +-36.7631 -50.2966 0.0 +-35.5977 -32.8759 0.0 +-38.2835 -13.9819 0.0 +-36793 3.2073 0.0 +-38.2603 20.0752 0.0 +-38.5707 36.7987 0.0 +-35.5583 50.5694 0.0 +-35.4349 71.8035 0.0 +-34977 87.9757 0.0 +-20.4709 -93.9292 0.0 +-20396 -76.3795 0.0 +-23.7124 -57.6671 0.0 +-23908 -40.2429 0.0 +-21.0551 -22.9972 0.0 +-20.3947 -5.09517 0.0 +-20212 9.72254 0.0 +-24.1934 27.5216 0.0 +-19.9637 47.0325 0.0 +-22.1828 63.9466 0.0 +-22.5229 76.7938 0.0 +-23.5922 96.0568 0.0 +-7.4651 -98934 0.0 +-5.63692 -80.5666 0.0 +-9.19887 -66.9041 0.0 +-7.83299 -49.2237 0.0 +-7.58047 -29.8289 0.0 +-9.17806 -14.3046 0.0 +-8.47321 1.79866 0.0 +-6.22558 17232 0.0 +-6.1979 34592 0.0 +-6.06438 51.8374 0.0 +-7.99579 71.8464 0.0 +-9.34279 88.3144 0.0 +8.09689 -92.1459 0.0 +8.56241 -76.1541 0.0 +5.75799 -58.7743 0.0 +6.80561 -42.8107 0.0 +7.29498 -24.8781 0.0 +9.49315 -8.67633 0.0 +5.48934 10.7607 0.0 +7.48964 25.2213 0.0 +8.37825 44.0642 0.0 +6.63344 60.3534 0.0 +8.82467 76.4799 0.0 +5.50615 97.7667 0.0 +20.0973 -85.3025 0.0 +23.6045 -66.3319 0.0 +21.6682 -31.0106 0.0 +20.8476 -14.3643 0.0 +23354 4.3938 0.0 +23.2571 21.0108 0.0 +21.9115 37.8318 0.0 +23.0517 50.9854 0.0 +20875 68.0303 0.0 +37.9323 -76.0558 0.0 +38.8808 -59.6001 0.0 +34.7439 -39.8988 0.0 +34.7224 -26.0126 0.0 +38066 -8.93449 0.0 +36.8513 9.09404 0.0 +36.4099 27.3056 0.0 +35.4412 44.6587 0.0 +36.1567 62.8726 0.0 +37.3943 79938 0.0 +50.5785 -83.4395 0.0 +51.7827 -64.4546 0.0 +50186 -48.5523 0.0 +53.3226 -33.7905 0.0 +51.1424 -15.7653 0.0 +52.0446 1.17038 0.0 +52.8293 16.8413 0.0 +49.7305 34.3809 0.0 +51.0108 55.0179 0.0 +49.4942 67.6309 0.0 +68.0616 -56.1252 0.0 +67.5741 -38.8549 0.0 +64.4535 -24.3364 0.0 +67.9933 -6.46241 0.0 +65.0907 12.6911 0.0 +65.1839 25.4711 0.0 +66.8215 44.4639 0.0 +67.4613 60.1319 0.0 +82.3886 -46.9536 0.0 +82.0147 -33.4872 0.0 +81.3058 -15.8358 0.0 +79.4432 -0.165153 0.0 +81181 18.2513 0.0 +82.8324 35.9978 0.0 +82.6803 55.3748 0.0 +95.6087 -23.2664 0.0 +96.4069 -7.70411 0.0 +94.1872 12.5324 0.0 +93.7707 26.4937 0.0 diff --git a/core/PhysiCell_cell.cpp b/core/PhysiCell_cell.cpp index 3780f0ea1..9253777ee 100644 --- a/core/PhysiCell_cell.cpp +++ b/core/PhysiCell_cell.cpp @@ -3525,6 +3525,171 @@ int find_cell_definition_index( int search_type ) return -1; } +// Definition of the output operator for parameters class for now i haven't included the pointer to phenotype. +std::ostream& operator<<(std::ostream& os, const Cell_Parameters& params) +{ + os << "o2_hypoxic_threshold: " << params.o2_hypoxic_threshold << std::endl; + os << "o2_hypoxic_response: " << params.o2_hypoxic_response << std::endl; + os << "o2_hypoxic_saturation: " << params.o2_hypoxic_saturation << std::endl; + os << "o2_proliferation_saturation: " << params.o2_proliferation_saturation << std::endl; + os << "o2_proliferation_threshold: " << params.o2_proliferation_threshold << std::endl; + os << "o2_reference: " << params.o2_reference << std::endl; + os << "o2_necrosis_threshold: " << params.o2_necrosis_threshold << std::endl; + os << "o2_necrosis_max: " << params.o2_necrosis_max << std::endl; + os << "max_necrosis_rate: " << params.max_necrosis_rate << std::endl; + os << "necrosis_type: " << params.necrosis_type << std::endl; + + // Add other class members if necessary + + return os; +} + +// Definition of the input operator for parameters class +std::istream& operator>>(std::istream& is, Cell_Parameters& params) +{ std::string dummy; + //o2_hypoxic_threshold + std::getline(is, dummy); + params.o2_hypoxic_threshold = read_number_in_line(dummy); + + //o2_hypoxic_response + std::getline(is, dummy); + params.o2_hypoxic_response = read_number_in_line(dummy); + + //o2_hypoxic_saturation + std::getline(is, dummy); + params.o2_hypoxic_saturation = read_number_in_line(dummy); + + //o2_proliferation_saturation + std::getline(is, dummy); + params.o2_proliferation_saturation = read_number_in_line(dummy); + + //o2_proliferation_threshold + std::getline(is, dummy); + params.o2_proliferation_threshold = read_number_in_line(dummy); + + //o2_reference + std::getline(is, dummy); + params.o2_reference = read_number_in_line(dummy); + + //o2_necrosis_threshold + std::getline(is, dummy); + params.o2_necrosis_threshold = read_number_in_line(dummy); + + //o2_necrosis_max + std::getline(is, dummy); + params.o2_necrosis_max = read_number_in_line(dummy); + + //max_necrosis_rate + std::getline(is, dummy); + params.max_necrosis_rate = read_number_in_line(dummy); + + //necrosis_type + std::getline(is, dummy); + params.necrosis_type = read_number_in_line_int(dummy); + + // skip empty rows + std::getline(is, dummy); + return is; +} + + +// Implementation of stream operator for cell_state class but no neighbors since not yet tracked (comment in the .h file) + +std::ostream& operator<<(std::ostream& os, const Cell_State& cellState) +{ + + // Output orientation + os << "Orientation: "; + for (double angle : cellState.orientation) + { + os << angle << " "; + } + os << std::endl; + + // Output simple pressure + os << "Simple_Pressure: " << cellState.simple_pressure << std::endl; + + // number of nuclei + os << "number_of_nuclei: " << cellState.number_of_nuclei<< std::endl; + + // total_attack_time + //added this if since sometimes the values were very low but not zero and they were creating problems + if (cellState.total_attack_time < std::numeric_limits::min()) { + os << "total_attack_time: " << 0.0<< std::endl; + } else { + os << "total_attack_time: " << cellState.total_attack_time<< std::endl; + } + + // contact_with_basement_membrane + os << "contact_with_basement_membrane: " << cellState.contact_with_basement_membrane; + + return os; +} + +std::istream& operator>>(std::istream& is, Cell_State& cellState) +{ + // Input orientation + std::string dummy; + std::getline(is, dummy); + + // Use a stringstream to parse the string + std::istringstream stream(dummy); + + // Ignore the "Orientation:" part + std::string key_1; + stream >> key_1; // Ignore the "Orientation:" + + double number; + + // Clear orientation + cellState.orientation.clear(); + + // Update orientation + while (stream >> number) { + cellState.orientation.push_back(number); + } + + + //Simple_Pressure + std::getline(is, dummy); + cellState.simple_pressure = read_number_in_line(dummy); + + + //number_of_nuclei + std::getline(is, dummy); + cellState.number_of_nuclei = read_number_in_line_int(dummy); + + + //total_attack_time + std::getline(is, dummy); + cellState.total_attack_time = read_number_in_line(dummy); + + //contact_with_basement_membrane + std::getline(is, dummy); + cellState.contact_with_basement_membrane = read_number_in_line(dummy); + + return is; +} + + +// Implementation of the insertion operator for cell class +std::ostream& operator<<(std::ostream& os, const Cell& cell) { + os << "Cell" << std::endl; + // Write the cell's attributes to 'os' + os << "ID: " << cell.ID << std::endl; + os << "index: " << cell.index << std::endl; + os << "type: " << cell.type << std::endl; + os << "is_active: " << (cell.is_active ? "true" : "false") << std::endl; + os << cell.type_name << "\n"; + os << "is_out_of_domain: " << (cell.is_out_of_domain ? "true" : "false") << std::endl; + os << "is_movable: " << (cell.is_movable ? "true" : "false") << std::endl; + os << cell.state << "\n"; + os << cell.parameters << "\n"; + os << cell.phenotype << "\n"; + os << cell.custom_data << "\n"; + + return os; +} }; diff --git a/core/PhysiCell_cell.h b/core/PhysiCell_cell.h index 573d93704..ef2f0cda0 100644 --- a/core/PhysiCell_cell.h +++ b/core/PhysiCell_cell.h @@ -111,6 +111,10 @@ class Cell_Parameters // necrosis parameters (may evenually be moved into a reference necrotic phenotype double max_necrosis_rate; // deprecate int necrosis_type; // deprecate + + // Declaration of input/output operators + friend std::ostream& operator<<(std::ostream& os, const Cell_Parameters& params); + friend std::istream& operator>>(std::istream& is, Cell_Parameters& params); Cell_Parameters(); }; @@ -159,6 +163,9 @@ class Cell_State bool contact_with_basement_membrane; // not implemented yet Cell_State(); + friend std::ostream& operator<<(std::ostream& os, const Cell_State& cellState); + friend std::istream& operator>>(std::istream& is, Cell_State& cellState); + }; class Cell : public Basic_Agent @@ -249,6 +256,9 @@ class Cell : public Basic_Agent std::vector nearby_interacting_cells( void ); // new in 1.8.0 void convert_to_cell_definition( Cell_Definition& cd ); + + // stream operator to save cells + friend std::ostream& operator<<(std::ostream& os, const Cell& cell); }; Cell* create_cell( Cell* (*custom_instantiate)() = NULL ); diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index 17851a5ea..49f12a11e 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -72,6 +72,7 @@ #include #include +#include using namespace BioFVM; @@ -425,4 +426,121 @@ Cell_Container* create_cell_container_for_microenvironment( BioFVM::Microenviron return cell_container; } -}; +void Cell_Container::save_data(std::ofstream& os){ + os << "Cell_container:" << std::endl; + + //Cells_ready_to_divide + os << "Cells_ready_to_divide: " << cells_ready_to_divide.size() << std::endl; + + //cells_ready_to_die + os << "cells_ready_to_die: " << cells_ready_to_die.size() << std::endl; + + //boundary_condition_for_pushed_out_agents + os << "boundary_condition_for_pushed_out_agents: " << boundary_condition_for_pushed_out_agents << std::endl; + + //initialzed + os << "initialzed: " << (initialzed ? "true" : "false") << std::endl; + + //boundary_condition_for_pushed_out_agents + os << "num_divisions_in_current_step: " << num_divisions_in_current_step << std::endl; + + //boundary_condition_for_pushed_out_agents + os << "num_deaths_in_current_step: " << num_deaths_in_current_step << std::endl; + + //boundary_condition_for_pushed_out_agents + os << "last_diffusion_time: " << last_diffusion_time << std::endl; + + //boundary_condition_for_pushed_out_agents + os << "last_cell_cycle_time: " << last_cell_cycle_time << std::endl; + + //boundary_condition_for_pushed_out_agents + os << "last_mechanics_time: " << last_mechanics_time << std::endl; + + //max_cell_interactive_distance_in_voxel + os << "max_cell_interactive_distance_in_voxel: "; + for (const auto& value : max_cell_interactive_distance_in_voxel) { + os << value << " "; + } + os << std::endl; + //cells_ready_to_die + os << "cells_ready_to_die: " << cells_ready_to_die.size() << std::endl; + + os << "cells_ready_to_divide: " << cells_ready_to_divide.size() << std::endl; +/* + //cells_ready_to_die + os << "cells_ready_to_die: "; + for (const auto& cell_ptr : cells_ready_to_die) { + os << cell_ptr << " "; + } + os << std::endl; + + //cells_ready_to_divide + os << "cells_ready_to_divide: "; + for (const auto& cell_ptr : cells_ready_to_divide) { + os << cell_ptr << " "; + } + os << std::endl;*/ +} + +void Cell_Container::reset_data(std::ifstream& is){ + std::string dummy; + + //skip first 3 lines + std::getline(is, dummy); + std::getline(is, dummy); + std::getline(is, dummy); + + //cell_cell_adhesion_strength + std::getline(is, dummy); + boundary_condition_for_pushed_out_agents = read_number_in_line(dummy); + + //initialzed + std::getline(is, dummy); + initialzed = read_number_in_line_bool(dummy); + + //num_divisions_in_current_step + std::getline(is, dummy); + num_divisions_in_current_step = read_number_in_line_int(dummy); + + //num_deaths_in_current_step + std::getline(is, dummy); + num_deaths_in_current_step = read_number_in_line_int(dummy); + + //last_diffusion_time + std::getline(is, dummy); + last_diffusion_time = read_number_in_line(dummy); + + //last_cell_cycle_time + std::getline(is, dummy); + last_cell_cycle_time = read_number_in_line(dummy); + + //last_mechanics_time + std::getline(is, dummy); + last_mechanics_time = read_number_in_line(dummy); + + // max_cell_interactive_distance_in_voxel + static int n_line = 0; + std::getline(is, dummy); + std::istringstream stream_max_cell_interactive_distance_in_voxel(dummy); + std::string key_max_cell_interactive_distance_in_voxel; + stream_max_cell_interactive_distance_in_voxel >> key_max_cell_interactive_distance_in_voxel; + double number; + std::vector max_cell_interactive_distance_in_voxel_new; + while(stream_max_cell_interactive_distance_in_voxel >> number) { + max_cell_interactive_distance_in_voxel_new.push_back(number); + n_line++; + } + assert(n_line == max_cell_interactive_distance_in_voxel.size()); + max_cell_interactive_distance_in_voxel.clear(); + max_cell_interactive_distance_in_voxel = max_cell_interactive_distance_in_voxel_new; + + + // cells_ready_to_die + std::getline(is, dummy); + + // cells_ready_to_divide + std::getline(is, dummy); + + +} +}; \ No newline at end of file diff --git a/core/PhysiCell_cell_container.h b/core/PhysiCell_cell_container.h index c1c7243c1..191ec7159 100644 --- a/core/PhysiCell_cell_container.h +++ b/core/PhysiCell_cell_container.h @@ -115,6 +115,9 @@ class Cell_Container : public BioFVM::Agent_Container void flag_cell_for_division( Cell* pCell ); void flag_cell_for_removal( Cell* pCell ); bool contain_any_cell(int voxel_index); + + void save_data(std::ofstream& os); + void reset_data(std::ifstream& is); }; int find_escaping_face_index(Cell* agent); diff --git a/core/PhysiCell_custom.cpp b/core/PhysiCell_custom.cpp index ffe23ccb6..a1a44e02c 100644 --- a/core/PhysiCell_custom.cpp +++ b/core/PhysiCell_custom.cpp @@ -69,6 +69,7 @@ #include #include #include +#include #include namespace PhysiCell @@ -85,7 +86,13 @@ Variable::Variable() std::ostream& operator<<(std::ostream& os, const Variable& v) { - os << v.name << ": " << v.value << " " << v.units; + if (v.units.empty()){ + os << v.name << ": " << v.value << " no_dim"; + } + else{ + os << v.name << ": " << v.value << " " << v.units; + } + return os; } @@ -238,21 +245,137 @@ double& Custom_Cell_Data::operator[]( std::string name ) return variables[ name_to_index_map[name] ].value; } +//std::ostream& operator<<(std::ostream& os, const Custom_Cell_Data& ccd) +//{ +// os << "Custom data (scalar): " << std::endl; +// for( int i=0 ; i < ccd.variables.size() ; i++ ) +// { +// os << i << ": " << ccd.variables[i] << std::endl; +// } +// +// os << "Custom data (vector): " << std::endl; +// for( int i=0 ; i < ccd.vector_variables.size() ; i++ ) +// { +// os << i << ": " << ccd.vector_variables[i] << std::endl; +// } +// +// return os; +//} + +//Operators for variable class only is because os was already defined +std::istream& operator>>(std::istream& is, PhysiCell::Variable& v) +{ + // Read name, value, and units directly from the line + is >> v.name >> v.value >> v.units; + //std::cout << "ayooo " << v.name << std::endl; + // Check if the input was successful + if (is.fail()) + { + std::cerr << v.name << " Error: Unable to read variable data." << std::endl; + } + + return is; +} + + +//Operators for Vector_variable class only is because os was already defined +std::istream& operator>>(std::istream& is, PhysiCell::Vector_Variable& v) +{ + // Read name + is >> v.name; + + // Read values until a non-numeric character is encountered + double value; + while (is >> value) + { + v.value.push_back(value); + } + + // Clear the failbit if it was set due to encountering a non-numeric character + is.clear(); + + // Read units + is >> v.units; + + return is; +} + + +//Operators for custom data class std::ostream& operator<<(std::ostream& os, const Custom_Cell_Data& ccd) { - os << "Custom data (scalar): " << std::endl; - for( int i=0 ; i < ccd.variables.size() ; i++ ) - { - os << i << ": " << ccd.variables[i] << std::endl; - } + for (int i = 0; i < ccd.variables.size(); i++) + { - os << "Custom data (vector): " << std::endl; - for( int i=0 ; i < ccd.vector_variables.size() ; i++ ) - { - os << i << ": " << ccd.vector_variables[i] << std::endl; - } - - return os; + os << ccd.variables[i] << std::endl; + } + + for (int i = 0; i < ccd.vector_variables.size(); i++) + { + os << ccd.vector_variables[i] << std::endl; + } + + os << "End of custom data" << std::endl; + return os; } +std::istream& operator>>(std::istream& is, PhysiCell::Custom_Cell_Data& ccd) +{ + // Read scalar data until a line starts with a space + + std::string line; + for (;;) + { + // Read the whole line + std::getline(is, line); + + + // Check if the first character is a space or the line is empty + if (line.find("End of custom data") != std::string::npos) + { + break; + } + + // Use a stringstream to extract values from the line + std::istringstream line_stream(line); + PhysiCell::Variable v; + line_stream >> v; + + + // Add the variable to ccd in this version it is saving the : in name. i have to remove it + std::string correct_name = v.name.substr(0, v.name.size() - 1); + + int index = ccd.find_variable_index(correct_name); + ccd[index] = v.value; + } + + // Read vector data until a line starts with a space + //std::cout << "Enter Custom data (vector): " << std::endl; + /* + for (;;) + { + // Read the whole line + std::string line; + std::getline(is, line); + + // Check if the first character is a space or the line is empty + if (line.empty() || std::isspace(line[0])) + { + // Break the loop if the first character is a space + break; + } + + // Use a stringstream to extract values from the line + std::istringstream line_stream(line); + PhysiCell::Vector_Variable vv; + line_stream >> vv; + + // Add the vector variable to ccd + ccd.add_vector_variable(vv); + } +*/ + return is; +} + + }; diff --git a/core/PhysiCell_custom.h b/core/PhysiCell_custom.h index aea0743db..c2619593f 100644 --- a/core/PhysiCell_custom.h +++ b/core/PhysiCell_custom.h @@ -81,6 +81,7 @@ class Variable { private: friend std::ostream& operator<<(std::ostream& os, const Variable& v); // done + friend std::istream& operator>>(std::istream& is, Variable& v); public: std::string name; double value; @@ -94,6 +95,7 @@ class Vector_Variable { private: friend std::ostream& operator<<(std::ostream& os, const Vector_Variable& v); // done + friend std::istream& operator>>(std::istream& is, Vector_Variable& v); public: std::string name; @@ -133,6 +135,10 @@ class Custom_Cell_Data Custom_Cell_Data(); // done Custom_Cell_Data( const Custom_Cell_Data& ccd ); + + //stream operator declaration + friend std::ostream& operator<<(std::ostream& os, const Custom_Cell_Data& ccd); + friend std::istream& operator>>(std::istream& is, Custom_Cell_Data& ccd); }; }; diff --git a/core/PhysiCell_phenotype.cpp b/core/PhysiCell_phenotype.cpp index a98559b7a..c7c0a4c6e 100644 --- a/core/PhysiCell_phenotype.cpp +++ b/core/PhysiCell_phenotype.cpp @@ -15,11 +15,75 @@ # able software used in your PhysiCell application. # # # # Because PhysiCell extensively uses BioFVM, we suggest you also cite BioFVM # -# as below: # -# # -# We implemented and solved the model using PhysiCell (Version x.y.z) [1], # -# with BioFVM [2] to solve the transport equations. # -# # +# as below: #double read_number_in_line(const std::string& line) { + std::istringstream stream(line); + std::string key; + stream >> key; + std::string numberStr; + + if (stream >> numberStr) { + try { + long double value = std::stold(numberStr); + + // Se il valore è fuori dal range di un double, lancia un'eccezione + if (value > std::numeric_limits::max() || value < -std::numeric_limits::max()) { + throw std::out_of_range("Value out of range for double: " + numberStr); + } + + return static_cast(value); + } catch (const std::invalid_argument&) { + throw std::runtime_error("Error: Failed to parse value from line: " + line); + } catch (const std::out_of_range&) { + throw std::runtime_error("Error: Value out of range for double in line: " + line); + } + } + + throw std::runtime_error("Error: Failed to parse value from line: " + line); +} + +int read_number_in_line_int(const std::string& line) { + std::istringstream stream(line); + std::string key; + stream >> key; + std::string numberStr; + + if (stream >> numberStr) { + try { + long long value = std::stoll(numberStr); // Usa `stoll` per numeri grandi + + // Se il valore è fuori dal range di un int, lancia un'eccezione + if (value > std::numeric_limits::max() || value < std::numeric_limits::min()) { + throw std::out_of_range("Value out of range for int: " + numberStr); + } + + return static_cast(value); + } catch (const std::invalid_argument&) { + throw std::runtime_error("Error: Failed to parse value from line: " + line); + } catch (const std::out_of_range&) { + throw std::runtime_error("Error: Value out of range for int in line: " + line); + } + } + + throw std::runtime_error("Error: Failed to parse value from line: " + line); +} + +bool read_number_in_line_bool(const std::string& line){ + std::istringstream stream(line); + std::string key, value; + stream >> key >> value; + + + if(value == "true"){ + return true; + } + else if(value == "false"){ + return false; + } + else{ + throw std::runtime_error("Failed to parse value from line"); + } +}; + # [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # # PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # # lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # @@ -1462,5 +1526,1069 @@ void Cell_Integrity::advance_damage( double dt ) return; } +// streaming operator for phase +std::ostream& operator<<(std::ostream& os, const Phase& phase) { + os << "index: " << phase.index << std::endl; + os << "code: " << phase.code << std::endl; + os << "name: " << phase.name << std::endl; + os << "division_at_phase_exit: " << (phase.division_at_phase_exit ? "true" : "false") << std::endl; + os << "removal_at_phase_exit: " << (phase.removal_at_phase_exit ? "true" : "false") << std::endl; + os << "entry_function: " << phase.entry_function << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Phase& phase){ + + std::string dummy; + + //index + std::getline(is, dummy); + phase.index = read_number_in_line_int(dummy); + + //code + std::getline(is, dummy); + phase.code = read_number_in_line_int(dummy); + + // name + std::getline(is, dummy); + std::istringstream stream(dummy); + std::string key, name; + stream >> key; + std::getline(stream, name); + name.erase(0, name.find_first_not_of(" ")); + phase.name = name; + + // division_at_phase_exit + std::getline(is, dummy); + phase.division_at_phase_exit = read_number_in_line_bool(dummy); + + // removal_at_phase_exit + std::getline(is, dummy); + phase.removal_at_phase_exit = read_number_in_line_bool(dummy); + + // entry_function + std::getline(is, dummy); + + //empty_row + std::getline(is, dummy); + + return is; +} + +// streaming operator for phase link +std::ostream& operator<<(std::ostream& os, const Phase_Link& phase_link) { + os << "start_phase_index: " << phase_link.start_phase_index << std::endl; + os << "end_phase_index: " << phase_link.end_phase_index << std::endl; + os << "fixed_duration: " << (phase_link.fixed_duration ? "true" : "false") << std::endl; + os << "arrest_function: " << (phase_link.arrest_function ? "true" : "false") << std::endl; + //not boolean need fixing + os << "exit_function: " << phase_link.exit_function << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Phase_Link& phase_link){ + + std::string dummy; + + + //start_phase_index + std::getline(is, dummy); + phase_link.start_phase_index = read_number_in_line_int(dummy); + + //end_phase_index + std::getline(is, dummy); + phase_link.end_phase_index = read_number_in_line_int(dummy); + + // fixed_duration + std::getline(is, dummy); + phase_link.fixed_duration = read_number_in_line_bool(dummy); + + // arrest_function + std::getline(is, dummy); + + // exit_function + std::getline(is, dummy); + + //empty_row + std::getline(is, dummy); + + return is; +} + +// Streaming operator for cycle model +void Cycle_Model::save_data(std::ostream& os ) +{ + + //name + os << "name: " << name << std::endl; + + //code + os << "code: " << code << std::endl; + + //phases + os << "phases:" << std::endl; + for (size_t i = 0; i < phases.size(); ++i) { + os << phases[i] << std::endl;; + } + os << std::endl; + + //phase_links + os << "phase_links:" << std::endl; + for (const auto& phase_links_v : phase_links) { + for( const auto& phase_links_current : phase_links_v){ + os << phase_links_current << std::endl; + }} + os << std::endl; + + //default_phase_index + os << "default_phase_index: " << default_phase_index << std::endl; + +} + +void Cycle_Model::reset_data(std::istream& is) +{ + + std::string dummy; + + // Skip "Cycle" and read subsequent data + std::getline(is, dummy); + + + // name + std::getline(is, dummy); + std::istringstream stream(dummy); + std::string key, name_new; + stream >> key >> name_new; + name_new = name; + + //code + std::getline(is, dummy); + code = read_number_in_line_int(dummy); + + + //phases + std::getline(is, dummy); + + for (size_t i = 0; i < phases.size(); ++i) { + is >> phases[i]; + }std::getline(is, dummy); + + //phase_links + std::getline(is, dummy); + + for (auto& phase_links_v : phase_links) { + for( auto& phase_links_current : phase_links_v){ + is >> phase_links_current; + }}std::getline(is, dummy); + + //default_phase_index + std::getline(is, dummy); + default_phase_index = read_number_in_line_int(dummy); + + +} +// Streaming operator for cycle data +std::ostream& operator<<(std::ostream& os, const Cycle_Data& cycleData) { + os << "Time_Units: " << cycleData.time_units << std::endl; + os << "Current_Phase_Index: " << cycleData.current_phase_index << std::endl; + os << "Elapsed_Time_in_Phase: " << cycleData.elapsed_time_in_phase << std::endl; + + os << "Transition_Rates:" << std::endl; + int num_phases = cycleData.transition_rates.size(); + os << "num_Phases: " << num_phases << std::endl; + for (std::size_t i = 0; i < cycleData.transition_rates.size(); ++i) { + for (std::size_t j = 0; j < cycleData.transition_rates[i].size(); ++j) { + os << "i" << i << "j" << j << ": " << cycleData.transition_rates[i][j] << std::endl; + } + } + os << "End of cycle data" << std::endl; + return os; +} + +std::istream& operator>>(std::istream& is, Cycle_Data& cycleData) { + std::string dummy; + + // Read Time Units + std::getline(is, dummy); + std::istringstream stream(dummy); + std::string key, time_units; + stream >> key >> time_units; + cycleData.time_units = time_units; + + // current_phase_index + std::getline(is, dummy); + cycleData.current_phase_index = read_number_in_line_int(dummy); + + // elapsed_time_in_phase + std::getline(is, dummy); + cycleData.elapsed_time_in_phase = read_number_in_line(dummy); + + // Read "Transition Rates:" + std::getline(is, dummy); + + //Read num phases + std::getline(is, dummy); + double num_phases; + num_phases = read_number_in_line(dummy); + + // no need to reupdate since now the dimension should be good. + + // Read the transition rates until "End of cycle data" is encountered + while (std::getline(is, dummy)) { + // Check if the line contains "End of cycle data" + if (dummy.find("End of cycle data") != std::string::npos) { + break; + } + //read the phase + int i,j; + + size_t pos_i = dummy.find('i'); + size_t pos_j = dummy.find('j'); + + std::string value_i = dummy.substr(pos_i + 1, 1); + std::string value_j = dummy.substr(pos_j + 1, 1); + + i = std::stoi(value_i); + j = std::stoi(value_j); + + std::istringstream stream_case1(dummy); + std::string key_case1; + stream_case1 >> key_case1; + double p1; + while (stream_case1 >> p1) { + cycleData.transition_rates[i][j] = p1; + } + } + + // discard empty rows; + + std::getline(is, dummy); + std::getline(is, dummy); + + return is; +} + +// Stream operator for Cycle +std::ostream& operator<<(std::ostream& os, const Cycle& cycle) { + //for now no pointers in the saved model + cycle.pCycle_Model->save_data(os); + os << cycle.data << std::endl; + return os; +} + +std::istream& operator>>(std::istream& is, Cycle& cycle) { + // read string "Cycle Data:" + cycle.pCycle_Model->reset_data(is); + is >> cycle.data; + + return is; +} + +// stream operator for death parameter class +std::ostream& operator<<(std::ostream& os, const Death_Parameters& DParameters) { + os << "Death_Parameters:" << std::endl; + os << "Time_Units: " << DParameters.time_units << std::endl; + os << "Unlysed_Fluid_Change_Rate: " << DParameters.unlysed_fluid_change_rate << std::endl; + os << "Lysed_Fluid_Change_Rate: " << DParameters.lysed_fluid_change_rate << std::endl; + os << "Cytoplasmic_Biomass_Change_Rate: " << DParameters.cytoplasmic_biomass_change_rate << std::endl; + os << "Nuclear_Biomass_Change_Rate: " << DParameters.nuclear_biomass_change_rate << std::endl; + os << "Calcification_Rate: " << DParameters.calcification_rate << std::endl; + os << "Relative_Rupture_Volume: " << DParameters.relative_rupture_volume << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Death_Parameters& DParameters) { + std::string dummy; + + // skip the header death parameters and parameters: + std::getline(is, dummy); + std::getline(is, dummy); + + // time_units + std::getline(is, dummy); + std::istringstream stream(dummy); + std::string key, time_units; + stream >> key >> time_units; + DParameters.time_units = time_units; + + //unlysed_fluid_change_rate + std::getline(is, dummy); + DParameters.unlysed_fluid_change_rate = read_number_in_line(dummy); + + // lysed_fluid_change_rate + std::getline(is, dummy); + DParameters.lysed_fluid_change_rate = read_number_in_line(dummy); + + // cytoplasmic_biomass_change_rate + std::getline(is, dummy); + DParameters.cytoplasmic_biomass_change_rate = read_number_in_line(dummy); + + // nuclear_biomass_change_rate + std::getline(is, dummy); + DParameters.nuclear_biomass_change_rate = read_number_in_line(dummy); + + //calcification_rate + std::getline(is, dummy); + DParameters.calcification_rate = read_number_in_line(dummy); + + // relative_rupture_volume + std::getline(is, dummy); + DParameters.relative_rupture_volume = read_number_in_line(dummy); + + + //skip empty row + std::getline(is, dummy); + + return is; +} + + +//Stream operator for death class +std::ostream& operator<<(std::ostream& os, const Death& death) { + os << "Death_Information:" << std::endl; + os << "Current_Death_Model_Index: " << death.current_death_model_index << std::endl; + os << "Dead: " << (death.dead ? "true" : "false") << std::endl; + + os << "Death_Models: " << death.models.size() << std::endl; + for (size_t i = 0; i < death.models.size(); ++i) { + os << "Model " << i << ":" << std::endl; + os << "Rate: " << death.rates[i] << std::endl; + os << "Parameters:\n" << death.parameters[i] << std::endl; + } + + return os; +} + +std::istream& operator>>(std::istream& is, Death& death) { + std::string dummy; + // Read the header "Death:" and "Death Information:" + std::getline(is, dummy); + std::getline(is, dummy); + + // current_death_model_index + std::getline(is, dummy); + death.current_death_model_index = read_number_in_line_int(dummy); + + // Dead status + std::getline(is, dummy); + death.dead = read_number_in_line_bool(dummy); + + // num_Death_Models + std::getline(is, dummy); + double num_Death_Models; + num_Death_Models = read_number_in_line(dummy); + + for (int i = 0; i < num_Death_Models; ++i) { + // skip the header Model: i + std::getline(is, dummy); + + // rate + std::getline(is, dummy); + death.rates[i] = read_number_in_line(dummy); + + // Parameters + is >> death.parameters[i]; + } + + // skip empty rows + std::getline(is, dummy); + return is; +} + +// Stream operator for volume class +std::ostream& operator<<(std::ostream& os, const Volume& volume) { + os << "Total: " << volume.total << std::endl; + os << "Solid: " << volume.solid << std::endl; + os << "Fluid: " << volume.fluid << std::endl; + os << "Fluid_Fraction: " << volume.fluid_fraction << std::endl; + os << "Nuclear: " << volume.nuclear << std::endl; + os << "Nuclear_Fluid: " << volume.nuclear_fluid << std::endl; + os << "Nuclear_Solid: " << volume.nuclear_solid << std::endl; + os << "Cytoplasmic: " << volume.cytoplasmic << std::endl; + os << "Cytoplasmic_Fluid: " << volume.cytoplasmic_fluid << std::endl; + os << "Cytoplasmic_Solid: " << volume.cytoplasmic_solid << std::endl; + os << "Calcified_Fraction: " << volume.calcified_fraction << std::endl; + os << "Cytoplasmic_to_Nuclear_Ratio: " << volume.cytoplasmic_to_nuclear_ratio << std::endl; + os << "Rupture_Volume: " << volume.rupture_volume << std::endl; + os << "Cytoplasmic_Biomass_Change_Rate: " << volume.cytoplasmic_biomass_change_rate << std::endl; + os << "Nuclear_Biomass_Change_Rate: " << volume.nuclear_biomass_change_rate << std::endl; + os << "Fluid_Change_Rate: " << volume.fluid_change_rate << std::endl; + os << "Calcification_Rate: " << volume.calcification_rate << std::endl; + os << "Target_Solid_Cytoplasmic: " << volume.target_solid_cytoplasmic << std::endl; + os << "Target_Solid_Nuclear: " << volume.target_solid_nuclear << std::endl; + os << "Target_Fluid_Fraction: " << volume.target_fluid_fraction << std::endl; + os << "Target_Cytoplasmic_to_Nuclear_Ratio: " << volume.target_cytoplasmic_to_nuclear_ratio << std::endl; + os << "Relative_Rupture_Volume: " << volume.relative_rupture_volume << std::endl; + + return os; +} +std::istream& operator>>(std::istream& is, Volume& volume) { + std::string dummy; + + //skip first 2 lines line + std::getline(is, dummy); + + //total + std::getline(is, dummy); + volume.total = read_number_in_line(dummy); + + //solid + std::getline(is, dummy); + volume.solid = read_number_in_line(dummy); + + //fluid + std::getline(is, dummy); + volume.fluid = read_number_in_line(dummy); + + //fluid_fraction + std::getline(is, dummy); + volume.fluid_fraction = read_number_in_line(dummy); + + //nuclear + std::getline(is, dummy); + volume.nuclear = read_number_in_line(dummy); + + //nuclear_fluid + std::getline(is, dummy); + volume.nuclear_fluid = read_number_in_line(dummy); + + //nuclear_solid + std::getline(is, dummy); + volume.nuclear_solid = read_number_in_line(dummy); + + //cytoplasmic + std::getline(is, dummy); + volume.cytoplasmic = read_number_in_line(dummy); + + //cytoplasmic_fluid + std::getline(is, dummy); + volume.cytoplasmic_fluid = read_number_in_line(dummy); + + //cytoplasmic_solid + std::getline(is, dummy); + volume.cytoplasmic_solid = read_number_in_line(dummy); + + //calcified_fraction + std::getline(is, dummy); + volume.calcified_fraction = read_number_in_line(dummy); + + //cytoplasmic_to_nuclear_ratio + std::getline(is, dummy); + volume.cytoplasmic_to_nuclear_ratio = read_number_in_line(dummy); + + //rupture_volume + std::getline(is, dummy); + volume.rupture_volume = read_number_in_line(dummy); + + //cytoplasmic_biomass_change_rate + std::getline(is, dummy); + volume.cytoplasmic_biomass_change_rate = read_number_in_line(dummy); + + //nuclear_biomass_change_rate + std::getline(is, dummy); + volume.nuclear_biomass_change_rate = read_number_in_line(dummy); + + //fluid_change_rate + std::getline(is, dummy); + volume.fluid_change_rate = read_number_in_line(dummy); + + //calcification_rate + std::getline(is, dummy); + volume.calcification_rate = read_number_in_line(dummy); + + //target_solid_cytoplasmic + std::getline(is, dummy); + volume.target_solid_cytoplasmic = read_number_in_line(dummy); + + //target_solid_nuclear + std::getline(is, dummy); + volume.target_solid_nuclear = read_number_in_line(dummy); + + //target_fluid_fraction + std::getline(is, dummy); + volume.target_fluid_fraction = read_number_in_line(dummy); + + //target_cytoplasmic_to_nuclear_ratio + std::getline(is, dummy); + volume.target_cytoplasmic_to_nuclear_ratio = read_number_in_line(dummy); + + //relative_rupture_volume + std::getline(is, dummy); + volume.relative_rupture_volume = read_number_in_line(dummy); + + //skip empty lines + std::getline(is, dummy); + return is; +} + +// Stream operator for geometry class +std::ostream& operator<<(std::ostream& os, const Geometry& geometry) { + os << "Radius: " << geometry.radius << std::endl; + os << "Nuclear_radius: " << geometry.nuclear_radius << std::endl; + os << "Surface_area: " << geometry.surface_area << std::endl; + os << "Polarity: " << geometry.polarity << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Geometry& geometry) { + std::string dummy; + + //skip first line + std::getline(is, dummy); + + //radius + std::getline(is, dummy); + geometry.radius = read_number_in_line(dummy); + + //nuclear_radius + std::getline(is, dummy); + geometry.nuclear_radius = read_number_in_line(dummy); + + //surface_area + std::getline(is, dummy); + geometry.surface_area = read_number_in_line(dummy); + + //polarity + std::getline(is, dummy); + geometry.polarity = read_number_in_line(dummy); + + + //skip last rows + std::getline(is, dummy); + return is; +} + +// Stream operator for Mechanics class +std::ostream& operator<<(std::ostream& os, const Mechanics& mechanics) { + os << "Cell_cell_adhesion_strength: " << mechanics.cell_cell_adhesion_strength << std::endl; + os << "Cell_BM_adhesion_strength: " << mechanics.cell_BM_adhesion_strength << std::endl; + os << "Cell_cell_repulsion_strength: " << mechanics.cell_cell_repulsion_strength << std::endl; + os << "Cell_BM_repulsion_strength: " << mechanics.cell_BM_repulsion_strength << std::endl; + os << "cell_adhesion_affinities: "; + for (const auto& affinity : mechanics.cell_adhesion_affinities) { + os << affinity << " "; + } + os << std::endl; + os << "Relative_maximum_adhesion_distance " << mechanics.relative_maximum_adhesion_distance << std::endl; + os << "maximum_number_of_attachments " << mechanics.maximum_number_of_attachments << std::endl; + os << "attachment_elastic_constant " << mechanics.attachment_elastic_constant << std::endl; + os << "attachment_rate " << mechanics.attachment_rate << std::endl; + os << "detachment_rate " << mechanics.detachment_rate << std::endl; + os << "relative_maximum_attachment_distance " << mechanics.relative_maximum_attachment_distance << std::endl; + os << "relative_detachment_distance " << mechanics.relative_detachment_distance << std::endl; + os << "maximum_attachment_rate " << mechanics.maximum_attachment_rate << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Mechanics& mechanics) { + std::string dummy; + + //skip first line + std::getline(is, dummy); + + //cell_cell_adhesion_strength + std::getline(is, dummy); + mechanics.cell_cell_adhesion_strength = read_number_in_line(dummy); + + //cell_BM_adhesion_strength + std::getline(is, dummy); + mechanics.cell_BM_adhesion_strength = read_number_in_line(dummy); + + //cell_cell_repulsion_strength + std::getline(is, dummy); + mechanics.cell_cell_repulsion_strength = read_number_in_line(dummy); + + //cell_BM_repulsion_strength + std::getline(is, dummy); + mechanics.cell_BM_repulsion_strength = read_number_in_line(dummy); + + // cell_adhesion_affinities + std::getline(is, dummy); + std::istringstream stream_cell_adhesion_affinities(dummy); + std::string key_cell_adhesion_affinities; + stream_cell_adhesion_affinities >> key_cell_adhesion_affinities; + double number; + std::vector cell_adhesion_affinities; + while(stream_cell_adhesion_affinities >> number) { + cell_adhesion_affinities.push_back(number); + } + mechanics.cell_adhesion_affinities.clear(); + mechanics.cell_adhesion_affinities = cell_adhesion_affinities; + + //relative_maximum_adhesion_distance + std::getline(is, dummy); + mechanics.relative_maximum_adhesion_distance = read_number_in_line(dummy); + + //maximum_number_of_attachments + std::getline(is, dummy); + mechanics.maximum_number_of_attachments = read_number_in_line_int(dummy); + + //attachment_elastic_constant + std::getline(is, dummy); + mechanics.attachment_elastic_constant = read_number_in_line(dummy); + + //attachment_rate + std::getline(is, dummy); + mechanics.attachment_rate = read_number_in_line(dummy); + + //detachment_rate + std::getline(is, dummy); + mechanics.detachment_rate = read_number_in_line(dummy); + + //relative_maximum_attachment_distance + std::getline(is, dummy); + mechanics.relative_maximum_attachment_distance = read_number_in_line(dummy); + //relative_detachment_distance + std::getline(is, dummy); + mechanics.relative_detachment_distance = read_number_in_line(dummy); + + //maximum_attachment_rate + std::getline(is, dummy); + mechanics.maximum_attachment_rate = read_number_in_line(dummy); + + //skip empty rows + std::getline(is, dummy); + return is; +} + +// stream operators for Motility class +std::ostream& operator<<(std::ostream& os, const Motility& motility) { + os << "Motility_Information:" << std::endl; + os << "Is_Motile: " << (motility.is_motile ? "true" : "false") << std::endl; + os << "Persistence_Time: " << motility.persistence_time << std::endl; + os << "Migration_Speed: " << motility.migration_speed << std::endl; + + os << "Migration_Bias_Direction: "; + for (const auto& direction : motility.migration_bias_direction) { + os << direction << " "; + } + os << std::endl; + + os << "Migration_Bias: " << motility.migration_bias << std::endl; + os << "Restrict_to_2D: " << (motility.restrict_to_2D ? "true" : "false") << std::endl; + + os << "Motility_Vector: "; + for (const auto& component : motility.motility_vector) { + os << component << " "; + } + os << std::endl; + + os << "Chemotaxis_Index: " << motility.chemotaxis_index << std::endl; + os << "Chemotaxis_Direction: " << motility.chemotaxis_direction << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Motility& motilityparams) { + std::string dummy; + + // read headers + std::getline(is, dummy); + std::getline(is, dummy); + + //is motile + std::getline(is, dummy); + motilityparams.is_motile = read_number_in_line_bool(dummy); + + //persistence_time + std::getline(is, dummy); + motilityparams.persistence_time = read_number_in_line(dummy); + + //migration_speed + std::getline(is, dummy); + motilityparams.migration_speed = read_number_in_line(dummy); + + // migration_bias_direction + std::getline(is, dummy); + std::istringstream stream_migration_bias_direction(dummy); + std::string key_migration_bias_direction; + stream_migration_bias_direction >> key_migration_bias_direction; + double number; + std::vector migration_bias_direction; + while(stream_migration_bias_direction >> number) { + migration_bias_direction.push_back(number); + } + motilityparams.migration_bias_direction.clear(); + motilityparams.migration_bias_direction = migration_bias_direction; + + //migration_bias + std::getline(is, dummy); + motilityparams.migration_bias = read_number_in_line(dummy); + + //restrict_to_2D + std::getline(is, dummy); + motilityparams.restrict_to_2D = read_number_in_line_bool(dummy); + + // motility_vector + std::getline(is, dummy); + std::istringstream stream_motility_vector(dummy); + std::string key_motility_vector; + stream_motility_vector >> key_motility_vector; + double number_2; + std::vector motility_vector; + while(stream_motility_vector >> number_2) { + motility_vector.push_back(number_2); + } + motilityparams.motility_vector.clear(); + motilityparams.motility_vector = motility_vector; + + //chemotaxis_index + std::getline(is, dummy); + motilityparams.chemotaxis_index = read_number_in_line_int(dummy); + + //chemotaxis_direction + std::getline(is, dummy); + motilityparams.chemotaxis_direction = read_number_in_line_int(dummy); + + //skip empty rows + std::getline(is, dummy); + + return is; +} + +// Streaming operators for Secretion class +std::ostream& operator<<(std::ostream& os, const Secretion& secretionparams) { + os << "Secretion_Information:" << std::endl; + os << "Number_of_Secretions: " << secretionparams.secretion_rates.size() << std::endl; + + for (size_t i = 0; i < secretionparams.secretion_rates.size(); ++i) { + os << "Secretion " << i << ":" << std::endl; + os << "Secretion_Rate: " << secretionparams.secretion_rates[i] << std::endl; + os << "Uptake_Rate: " << secretionparams.uptake_rates[i] << std::endl; + os << "Saturation_Density: " << secretionparams.saturation_densities[i] << std::endl; + os << "Net_Export_Rate: " << secretionparams.net_export_rates[i] << "\n" << std::endl; + } + + return os; +} + +std::istream& operator>>(std::istream& is, Secretion& secretionparams) { + std::string dummy; + + // read headers + std::getline(is, dummy); + std::getline(is, dummy); + + //number_of_secretions + std::getline(is, dummy); + double number_of_secretions; + number_of_secretions = read_number_in_line(dummy); + + for (size_t i = 0; i < number_of_secretions; ++i) { + + //skip header + std::getline(is, dummy); + + //secretion_rate + std::getline(is, dummy); + secretionparams.secretion_rates[i] = read_number_in_line(dummy); + + //uptake_rate + std::getline(is, dummy); + secretionparams.uptake_rates[i] = read_number_in_line(dummy); + + //saturation_density + std::getline(is, dummy); + secretionparams.saturation_densities[i] = read_number_in_line(dummy); + + //net_export_rate + std::getline(is, dummy); + secretionparams.net_export_rates[i] = read_number_in_line(dummy); + + //skip empty rows + std::getline(is, dummy); + + } + + //skip last rows + std::getline(is, dummy); + + return is; +} + +// Streaming operators for class Molecular +std::ostream& operator<<(std::ostream& os, const Molecular& molecular) { + os << "Molecular_Information:" << std::endl; + os << "number_of_substrates: " << molecular.internalized_total_substrates.size() << std::endl; + + for (size_t i = 0; i < molecular.internalized_total_substrates.size(); ++i){ + + os << "Substrate: " << i << std::endl; + os << "Internalized_Total_Substrates: " << molecular.internalized_total_substrates[i] << std::endl; + os << "Fraction_Released_at_Death: " << molecular.fraction_released_at_death[i] << std::endl; + os << "Fraction_Transferred_When_Ingested: " << molecular.fraction_transferred_when_ingested[i] << "\n" << std::endl; + + } + + return os; +} + +std::istream& operator>>(std::istream& is, Molecular& molecular) { + std::string dummy; + + // read headers + std::getline(is, dummy); + std::getline(is, dummy); + + //number_of_secretions + std::getline(is, dummy); + double number_of_substrates; + number_of_substrates = read_number_in_line(dummy); + + for (size_t i = 0; i < number_of_substrates; ++i) { + + //skip header + std::getline(is, dummy); + + //Internalized_Total_Substrates + std::getline(is, dummy); + molecular.internalized_total_substrates[i] = read_number_in_line(dummy); + + //Fraction_Released_at_Death + std::getline(is, dummy); + molecular.fraction_released_at_death[i] = read_number_in_line(dummy); + + //Fraction_Transferred_When_Ingested + std::getline(is, dummy); + molecular.fraction_released_at_death[i] = read_number_in_line(dummy); + //skip empty rows + std::getline(is, dummy); + + } + + //skip last rows + std::getline(is, dummy); + std::getline(is, dummy); + return is; +} + +std::ostream& operator<<(std::ostream& os, const Cell_Interactions& cell_interactions) { + os << "Cell_Interactions_Information:" << std::endl; + os << "Dead_Phagocytosis_Rate: " << cell_interactions.other_dead_phagocytosis_rate << std::endl; + os << "Live_Phagocytosis_Rates: "; + for (size_t i = 0; i < cell_interactions.live_phagocytosis_rates.size(); ++i) { + os << cell_interactions.live_phagocytosis_rates[i] << " "; + } + os << std::endl; + os << "Attack_Rates: "; + for (size_t i = 0; i < cell_interactions.attack_rates.size(); ++i) { + os << cell_interactions.attack_rates[i] << " "; + } + os << std::endl; + os << "Immunogenicities: "; + for (size_t i = 0; i < cell_interactions.immunogenicities.size(); ++i) { + os << cell_interactions.immunogenicities[i] << " "; + } + os << std::endl; + + os << "Damage_Rate: " << cell_interactions.attack_damage_rate << std::endl; + os << "Fusion_Rates: "; + for (size_t i = 0; i < cell_interactions.fusion_rates.size(); ++i) { + os << cell_interactions.fusion_rates[i] << " "; + } + os << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Cell_Interactions& cell_interactions){ + std::string dummy; + + // read headers + std::getline(is, dummy); + + //other_dead_phagocytosis_rate + std::getline(is, dummy); + cell_interactions.other_dead_phagocytosis_rate = read_number_in_line(dummy); + + // Live_Phagocytosis_Rates + std::getline(is, dummy); + std::istringstream stream_Live_Phagocytosis_Rates(dummy); + std::string key_Live_Phagocytosis_Rates; + stream_Live_Phagocytosis_Rates >> key_Live_Phagocytosis_Rates; + double number1; + std::vector Live_Phagocytosis_Rates; + while(stream_Live_Phagocytosis_Rates >> number1) { + Live_Phagocytosis_Rates.push_back(number1); + } + cell_interactions.live_phagocytosis_rates.clear(); + cell_interactions.live_phagocytosis_rates = Live_Phagocytosis_Rates; + + // Attack_Rates + std::getline(is, dummy); + std::istringstream stream_Attack_Rates(dummy); + std::string key_Attack_Rates; + stream_Attack_Rates >> key_Attack_Rates; + double number2; + std::vector Attack_Rates; + while(stream_Attack_Rates >> number2) { + Attack_Rates.push_back(number2); + } + cell_interactions.attack_rates.clear(); + cell_interactions.attack_rates = Attack_Rates; + + // Immunogenicities + std::getline(is, dummy); + std::istringstream stream_Immunogenicities(dummy); + std::string key_Immunogenicities; + stream_Immunogenicities >> key_Immunogenicities; + double number3; + std::vector Immunogenicities; + while(stream_Immunogenicities >> number3) { + Immunogenicities.push_back(number3); + } + cell_interactions.immunogenicities.clear(); + cell_interactions.immunogenicities = Immunogenicities; + + //Damage_Rate + std::getline(is, dummy); + cell_interactions.attack_damage_rate = read_number_in_line(dummy); + + // Fusion_Rates + std::getline(is, dummy); + std::istringstream stream_Fusion_Rates(dummy); + std::string key_Fusion_Rates; + stream_Fusion_Rates >> key_Fusion_Rates; + double number4; + std::vector Fusion_Rates; + while(stream_Fusion_Rates >> number4) { + Fusion_Rates.push_back(number4); + } + cell_interactions.fusion_rates.clear(); + cell_interactions.fusion_rates = Fusion_Rates; + + std::getline(is, dummy); + + + return is; +} + +std::ostream& operator<<(std::ostream& os, const Cell_Transformations& cell_transformations) { + os << "Cell_Transformations_Information:" << std::endl; + os << "Transformation_Rates: "; + for (size_t i = 0; i < cell_transformations.transformation_rates.size(); ++i) { + os << cell_transformations.transformation_rates[i] << " "; + } + os << std::endl; + + return os; +} + +std::istream& operator>>(std::istream& is, Cell_Transformations& cell_transformations) { + std::string dummy; + + // read headers + std::getline(is, dummy); + std::getline(is, dummy); + + // Transformation_Rates + std::getline(is, dummy); + std::istringstream stream_Transformation_Rates(dummy); + std::string key_Transformation_Rates; + stream_Transformation_Rates >> key_Transformation_Rates; + double number; + std::vector Transformation_Rates; + while(stream_Transformation_Rates >> number) { + Transformation_Rates.push_back(number); + } + cell_transformations.transformation_rates.clear(); + cell_transformations.transformation_rates = Transformation_Rates; + + std::getline(is, dummy); + std::getline(is, dummy); + return is; + +} + +// Stream operator for Phenotyope class +std::ostream& operator<<(std::ostream& os, const Phenotype& phenotype) { + os << "Flagged_for_division: " << (phenotype.flagged_for_division ? "true" : "false") << std::endl; + os << "Flagged_for_removal: " << (phenotype.flagged_for_removal ? "true" : "false") << std::endl; + os << "Death:\n" << phenotype.death << std::endl; + os << "Cycle:\n" << phenotype.cycle << std::endl; + os << "Volume:\n" << phenotype.volume << std::endl; + os << "Geometry:\n" << phenotype.geometry << std::endl; + os << "Mechanics:\n" << phenotype.mechanics << std::endl; + os << "Motility:\n" << phenotype.motility << std::endl; + os << "Secretion:\n" << phenotype.secretion << std::endl; + os << "Molecular:\n" << phenotype.molecular << std::endl; + os << "Cell_Interactions:\n" << phenotype.cell_interactions << std::endl; + os << "Cell_Transformations:\n" << phenotype.cell_transformations << std::endl; + return os; + } + +}; +double read_number_in_line(const std::string& line) { + std::istringstream stream(line); + std::string key; + stream >> key; + std::string numberStr; + + if (stream >> numberStr) { + try { + long double value = std::stold(numberStr); + + // Se il valore è fuori dal range di un double, lancia un'eccezione + if (value > std::numeric_limits::max() || value < -std::numeric_limits::max()) { + throw std::out_of_range("Value out of range for double: " + numberStr); + } + + return static_cast(value); + } catch (const std::invalid_argument&) { + throw std::runtime_error("Error: Failed to parse value from line: " + line); + } catch (const std::out_of_range&) { + throw std::runtime_error("Error: Value out of range for double in line: " + line); + } + } + + throw std::runtime_error("Error: Failed to parse value from line: " + line); +} + +int read_number_in_line_int(const std::string& line) { + std::istringstream stream(line); + std::string key; + stream >> key; + std::string numberStr; + + if (stream >> numberStr) { + try { + long long value = std::stoll(numberStr); // Usa `stoll` per numeri grandi + + // Se il valore è fuori dal range di un int, lancia un'eccezione + if (value > std::numeric_limits::max() || value < std::numeric_limits::min()) { + throw std::out_of_range("Value out of range for int: " + numberStr); + } + + return static_cast(value); + } catch (const std::invalid_argument&) { + throw std::runtime_error("Error: Failed to parse value from line: " + line); + } catch (const std::out_of_range&) { + throw std::runtime_error("Error: Value out of range for int in line: " + line); + } + } + + throw std::runtime_error("Error: Failed to parse value from line: " + line); +} + +bool read_number_in_line_bool(const std::string& line){ + std::istringstream stream(line); + std::string key, value; + stream >> key >> value; + + + if(value == "true"){ + return true; + } + else if(value == "false"){ + return false; + } + else{ + throw std::runtime_error("Failed to parse value from line"); + } }; diff --git a/core/PhysiCell_phenotype.h b/core/PhysiCell_phenotype.h index c69b19400..4eece7b4c 100644 --- a/core/PhysiCell_phenotype.h +++ b/core/PhysiCell_phenotype.h @@ -72,6 +72,7 @@ #include #include #include +#include #include "../BioFVM/BioFVM.h" @@ -115,6 +116,9 @@ class Phase void (*entry_function)( Cell* pCell, Phenotype& phenotype, double dt ); Phase(); // done + + friend std::ostream& operator<<(std::ostream& os, const Phase& phase); + friend std::istream& operator>>(std::istream& is, Phase& phase); }; class Phase_Link @@ -132,6 +136,9 @@ class Phase_Link // function to be excecuted when completing the phase transition Phase_Link(); // done + + friend std::ostream& operator<<(std::ostream& os, const Phase_Link& phase_link); + friend std::istream& operator>>(std::istream& is, Phase_Link& phase_link); }; class Cycle_Data @@ -169,6 +176,8 @@ 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). + friend std::ostream& operator<<(std::ostream& os, const Cycle_Data& cycleData); + friend std::istream& operator>>(std::istream& is, Cycle_Data& cycleData); }; class Cycle_Model @@ -211,6 +220,9 @@ class Cycle_Model Phase_Link& phase_link(int start_index,int end_index ); // done std::ostream& display( std::ostream& os ); // done + + void save_data(std::ostream& os ); + void reset_data(std::istream& is ); }; class Asymmetric_Division @@ -247,6 +259,8 @@ class Cycle int& current_phase_index( void ); // done void sync_to_cycle_model( Cycle_Model& cm ); // done + friend std::ostream& operator<<(std::ostream& os, const Cycle& cycle_params); + friend std::istream& operator>>(std::istream& is, Cycle& cycle_params); Asymmetric_Division asymmetric_division; }; @@ -267,6 +281,9 @@ class Death_Parameters double relative_rupture_volume; Death_Parameters(); // done + + friend std::ostream& operator<<(std::ostream& os, const Death_Parameters& DParameters); + friend std::istream& operator>>(std::istream& is, Death_Parameters& DParameters); }; class Death @@ -297,6 +314,9 @@ class Death // ease of access double& apoptosis_rate(void); double& necrosis_rate(void); + + friend std::ostream& operator<<(std::ostream& os, const Death& deathparams); + friend std::istream& operator>>(std::istream& is, Death& deathparams); }; class Volume @@ -355,6 +375,8 @@ class Volume void divide( void ); // done void multiply_by_ratio(double); // done + friend std::ostream& operator<<(std::ostream& os, const Volume& volumeparams); + friend std::istream& operator>>(std::istream& is, Volume& volumeparams); }; class Geometry @@ -373,6 +395,8 @@ class Geometry void update_surface_area( Cell* pCell, Phenotype& phenotype, double dt ); // done void update( Cell* pCell, Phenotype& phenotype, double dt ); // done + friend std::ostream& operator<<(std::ostream& os, const Geometry& geometryparams); + friend std::istream& operator>>(std::istream& is, Geometry& geometryparams); }; class Mechanics @@ -416,7 +440,8 @@ class Mechanics void set_absolute_equilibrium_distance( Phenotype& phenotype, double new_value ); // done - + friend std::ostream& operator<<(std::ostream& os, const Mechanics& mechanicsparams); + friend std::istream& operator>>(std::istream& is, Mechanics& mechanicsparams); }; class Motility @@ -451,6 +476,9 @@ class Motility Motility(); // done + + friend std::ostream& operator<<(std::ostream& os, const Motility& motilityparams); + friend std::istream& operator>>(std::istream& is, Motility& motilityparams); }; class Secretion @@ -486,7 +514,10 @@ class Secretion double& secretion_rate( std::string name ); double& uptake_rate( std::string name ); double& saturation_density( std::string name ); - double& net_export_rate( std::string name ); + double& net_export_rate( std::string name ); + + friend std::ostream& operator<<(std::ostream& os, const Secretion& secretionparams); + friend std::istream& operator>>(std::istream& is, Secretion& secretionparams); }; class Cell_Functions @@ -528,7 +559,9 @@ class Cell_Functions void (*plot_agent_SVG)(std::ofstream& os, Cell* pCell, double z_slice, std::vector (*cell_coloring_function)(Cell*), double X_lower, double Y_lower); void (*plot_agent_legend)(std::ofstream& os, Cell_Definition* cell_def, double& cursor_x, double& cursor_y, std::vector (*cell_coloring_function)(Cell*), double temp_cell_radius); - + friend std::ostream& operator<<(std::ostream& os, const Cell_Functions& cellFunctions); + friend std::istream& operator>>(std::istream& is, Cell_Functions& cellFunctions); + Cell_Functions(); // done }; @@ -610,7 +643,10 @@ class Molecular void sync_to_cell( Basic_Agent* pCell ); // ease of access - double& internalized_total_substrate( std::string name ); + double& internalized_total_substrate( std::string name ); + + friend std::ostream& operator<<(std::ostream& os, const Molecular& molecular); + friend std::istream& operator>>(std::istream& is, Molecular& molecular); }; @@ -661,9 +697,16 @@ class Intracellular virtual bool has_variable(std::string name) = 0; virtual bool get_boolean_variable_value(std::string name) = 0; virtual void set_boolean_variable_value(std::string name, bool value) = 0; + virtual int get_number_of_nodes() = 0; // virtual bool get_double_variable_value(std::string name) = 0; // virtual void set_double_variable_value(std::string name, bool value) = 0; virtual void print_current_nodes() = 0; + virtual void save_current_nodes(std::ostream& out_stream) = 0; + virtual void save_current_parameters(std::ostream& out_stream) = 0; + virtual void read_current_parameter(std::ifstream& in_stream) = 0; + virtual void save_current_parameters_maboss(std::ostream& out_stream) = 0; + virtual void read_current_parameter_maboss(std::ifstream& in_stream) = 0; + virtual void reinit_maboss(std::string networkFile, std::string configFile) = 0; // ================ specific to "roadrunner" ================ @@ -716,7 +759,10 @@ class Cell_Interactions double& immunogenicity( std::string type_name ); // done // automated cell phagocytosis, attack, and fusion -// void perform_interactions( Cell* pCell, Phenotype& phenotype, double dt ); +// void perform_interactions( Cell* pCell, Phenotype& phenotype, double dt ); + + friend std::ostream& operator<<(std::ostream& os, const Cell_Interactions& cell_interactions); + friend std::istream& operator>>(std::istream& is, Cell_Interactions& cell_interactions); }; class Cell_Transformations @@ -735,6 +781,9 @@ class Cell_Transformations // automated cell transformations // void perform_transformations( Cell* pCell, Phenotype& phenotype, double dt ); + + friend std::ostream& operator<<(std::ostream& os, const Cell_Transformations& cell_transformations); + friend std::istream& operator>>(std::istream& is, Cell_Transformations& cell_transformations); }; // pre-beta functionality in 1.10.3 @@ -766,6 +815,9 @@ class Cell_Integrity Cell_Integrity(); void advance_damage( double dt ); + + friend std::ostream& operator<<(std::ostream& os, const Cell_Integrity& integrity); + friend std::istream& operator>>(std::istream& is, Cell_Integrity& integrity); }; class Phenotype @@ -805,8 +857,13 @@ class Phenotype // make sure cycle, death, etc. are synced to the defaults. void sync_to_default_functions( void ); // done + + friend std::ostream& operator<<(std::ostream& os, const Phenotype& pheno_params); }; }; #endif +double read_number_in_line(const std::string& line); +int read_number_in_line_int(const std::string& line); +bool read_number_in_line_bool(const std::string& line); diff --git a/core/PhysiCell_utilities.cpp b/core/PhysiCell_utilities.cpp index 97875642a..490597b7c 100644 --- a/core/PhysiCell_utilities.cpp +++ b/core/PhysiCell_utilities.cpp @@ -72,6 +72,10 @@ namespace PhysiCell{ +int counter_normal_random = 0; +int counter_int_random = 0; +int counter_double_random = 0; + thread_local std::mt19937_64 physicell_PRNG_generator; thread_local bool local_pnrg_setup_done = false; @@ -359,6 +363,53 @@ int choose_event( std::vector& probabilities ) return probabilities.size(); } +void save_counters(std::ostream& out_stream) +{ + out_stream << "counter_double_random: " << counter_double_random << std::endl; + out_stream << "counter_int_random: " << counter_int_random << std::endl; + out_stream << "counter_normal_random: " << counter_normal_random << std::endl; +} + +void print_counters() +{ + std::cout << "counter_double_random: " << counter_double_random << std::endl; + std::cout << "counter_int_random: " << counter_int_random << std::endl; + std::cout << "counter_normal_random: " << counter_normal_random << std::endl; +} + +void set_counter_double_random(int value) +{ + if (counter_double_random > value){ + std::cout << "UniformRandom counter still greater than value" << std::endl; + return; + } + while (counter_double_random != value){ + UniformRandom(); + } +} + +void set_counter_int_random(int value) +{ + if (counter_int_random > value){ + std::cout << "UniformInt counter still greater than value" << std::endl; + return; + } + while (counter_int_random != value){ + UniformInt(); + } +} + +void set_counter_normal_random(int value) +{ + if (counter_normal_random > value){ + std::cout << "NormalRandom counter still greater than value" << std::endl; + return; + } + while (counter_normal_random != value){ + NormalRandom(0.5, 1); + } +} + void copy_file_to_output(std::string filename) { std::cout << "Copying " << filename << " to output folder." << std::endl; diff --git a/core/PhysiCell_utilities.h b/core/PhysiCell_utilities.h index c17dd9217..77282ce84 100644 --- a/core/PhysiCell_utilities.h +++ b/core/PhysiCell_utilities.h @@ -111,7 +111,18 @@ void add_software_citation( std::string name , std::string version, std::string int choose_event( std::vector& probabilities ); -void copy_file_to_output( std::string filename ); +void save_counters(std::ostream& out_stream); + +void set_counter_double_random(int value); + +void set_counter_int_random(int value); + +void set_counter_normal_random(int value); + +void print_counters(); + +void copy_file_to_output(std::string filename); + }; #endif diff --git a/output/empty.txt b/output/empty.txt deleted file mode 100644 index e69de29bb..000000000 diff --git a/sample_projects/Makefile-default b/sample_projects/Makefile-default index 76ad9b813..fc6c6816a 100644 --- a/sample_projects/Makefile-default +++ b/sample_projects/Makefile-default @@ -77,7 +77,7 @@ list-projects: @echo " asymmetric-division-sample" @echo "" @echo "Sample intracellular projects: template_BM ode-energy-sample physiboss-cell-lines-sample" - @echo " cancer-metabolism-sample physiboss-tutorial physiboss-tutorial-invasion" + @echo " physiboss-tnf-model cancer-metabolism-sample physiboss-tutorial physiboss-tutorial-invasion" @echo "" template: @@ -235,6 +235,16 @@ physiboss-cell-lines-sample: cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml cp ./sample_projects_intracellular/boolean/physiboss_cell_lines/config/* ./config/ +physiboss-tnf-model: + cp ./sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/* ./custom_modules/ + touch main.cpp && cp main.cpp main-backup.cpp + cp ./sample_projects_intracellular/boolean/spheroid_tnf_model/main-spheroid_TNF.cpp ./main.cpp + cp Makefile Makefile-backup + cp ./sample_projects_intracellular/boolean/spheroid_tnf_model/Makefile . + cp ./config/PhysiCell_settings.xml ./config/PhysiCell_settings-backup.xml + cp ./sample_projects_intracellular/boolean/spheroid_tnf_model/config/* ./config/ + cp -r ./sample_projects_intracellular/boolean/spheroid_tnf_model/scripts ./ + physiboss-tutorial: cp ./sample_projects_intracellular/boolean/tutorial/custom_modules/* ./custom_modules/ touch main.cpp && cp main.cpp main-backup.cpp diff --git a/sample_projects/biorobots/Makefile b/sample_projects/biorobots/Makefile index 26314225c..4f149880b 100644 --- a/sample_projects/biorobots/Makefile +++ b/sample_projects/biorobots/Makefile @@ -62,8 +62,10 @@ PhysiCell_custom_module_OBJECTS := custom.o pugixml_OBJECTS := pugixml.o +start_and_stop_OBJECTS := start_and_stop.o + PhysiCell_OBJECTS := $(BioFVM_OBJECTS) $(pugixml_OBJECTS) $(PhysiCell_core_OBJECTS) $(PhysiCell_module_OBJECTS) -ALL_OBJECTS := $(PhysiCell_OBJECTS) $(PhysiCell_custom_module_OBJECTS) +ALL_OBJECTS := $(PhysiCell_OBJECTS) $(PhysiCell_custom_module_OBJECTS) $(start_and_stop_OBJECTS) # compile the project @@ -165,7 +167,11 @@ PhysiCell_basic_signaling.o: ./core/PhysiCell_basic_signaling.cpp PhysiCell_geometry.o: ./modules/PhysiCell_geometry.cpp $(COMPILE_COMMAND) -c ./modules/PhysiCell_geometry.cpp - + +# start_and_stop +start_and_stop.o: ./addons/start_and_stop/start_and_stop.cpp + $(COMPILE_COMMAND) $(INC) -c ./addons/start_and_stop/start_and_stop.cpp + # user-defined PhysiCell modules custom.o: ./custom_modules/custom.cpp diff --git a/sample_projects/biorobots/config/PhysiCell_settings.xml b/sample_projects/biorobots/config/PhysiCell_settings.xml index 2d5f7e46b..844095d1b 100644 --- a/sample_projects/biorobots/config/PhysiCell_settings.xml +++ b/sample_projects/biorobots/config/PhysiCell_settings.xml @@ -551,6 +551,7 @@ + false 0 0 red diff --git a/sample_projects/biorobots/custom_modules/custom.h b/sample_projects/biorobots/custom_modules/custom.h index 13c67b857..294027d16 100644 --- a/sample_projects/biorobots/custom_modules/custom.h +++ b/sample_projects/biorobots/custom_modules/custom.h @@ -67,6 +67,8 @@ #include "../core/PhysiCell.h" #include "../modules/PhysiCell_standard_modules.h" +//include the start and stop add-on +#include "../addons/start_and_stop/start_and_stop.h" using namespace BioFVM; using namespace PhysiCell; diff --git a/sample_projects/biorobots/main.cpp b/sample_projects/biorobots/main.cpp index 2e5a302f0..b6fa3f08b 100644 --- a/sample_projects/biorobots/main.cpp +++ b/sample_projects/biorobots/main.cpp @@ -115,6 +115,8 @@ int main( int argc, char* argv[] ) setup_microenvironment(); // modify this in the custom code + bool start_stop = parameters.bools("start_stop"); + /* PhysiCell setup */ // set mechanics voxel size, and match the data structure to BioFVM @@ -125,7 +127,21 @@ int main( int argc, char* argv[] ) create_cell_types(); - setup_tissue(); + if( start_stop ){ + + // reset cells as they were in the previous simulation + setup_tissue(); + + reset_cell(cell_container->last_cell_cycle_time); + //exit(-1); + + reset_global_parameters(cell_container); + + reset_microenv(); + + } else{ + setup_tissue(); //death model index = 1 == necrotic...= 0 == apoptotic. + } /* Users typically stop modifying here. END USERMODS */ @@ -173,6 +189,9 @@ int main( int argc, char* argv[] ) report_file<<"simulated time\tnum cells\tnum division\tnum death\twall time"< - 2880.0 + 2880 min micron 0.01 @@ -29,11 +29,11 @@ output - 10 + 30 true - 10 + 30 true false @@ -670,6 +670,9 @@ + false + false + 0 30.0 30 diff --git a/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom.cpp b/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom.cpp index d36cd42f8..c8fdd6355 100644 --- a/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom.cpp +++ b/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom.cpp @@ -399,3 +399,19 @@ std::string my_coloring_function_for_stroma( double concentration, double max_co return paint_by_density_percentage( concentration, max_conc, min_conc); } + +bool auto_stop() { + + bool condition = false; + bool stop; + float stop_time = 1440; + if ((fabs( PhysiCell_globals.current_time - PhysiCell_globals.next_full_save_time ) < 0.01 * diffusion_dt) && (PhysiCell_globals.next_full_save_time == stop_time)) { + condition = true; + } + if (condition) { + stop = true; + } else { + stop = false; + } + return stop; +} diff --git a/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom.h b/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom.h index ccc4e80e0..abde164fc 100644 --- a/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom.h +++ b/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom.h @@ -101,4 +101,5 @@ void pre_update_intracellular(Cell* pCell, Phenotype& phenotype, double dt); void post_update_intracellular(Cell* pCell, Phenotype& phenotype, double dt); std::string my_coloring_function_for_stroma( double concentration, double max_conc, double min_conc ); -void color_node(Cell* pCell); \ No newline at end of file +void color_node(Cell* pCell); +bool auto_stop(); diff --git a/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom_main.h b/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom_main.h index 24db740f8..8e94d8b43 100755 --- a/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom_main.h +++ b/sample_projects_intracellular/boolean/cancer_invasion/custom_modules/custom_main.h @@ -2,6 +2,8 @@ #define __custom_main_h__ #include "../core/PhysiCell.h" +#include "../addons/start_and_stop/start_and_stop.h" + inline double current_value( double min, double max, double percent ) { return (min + (max-min) * percent); }; diff --git a/sample_projects_intracellular/boolean/cancer_invasion/main.cpp b/sample_projects_intracellular/boolean/cancer_invasion/main.cpp index eac3e89cf..73fdd42b8 100644 --- a/sample_projects_intracellular/boolean/cancer_invasion/main.cpp +++ b/sample_projects_intracellular/boolean/cancer_invasion/main.cpp @@ -86,7 +86,13 @@ using namespace PhysiCell; int main( int argc, char* argv[] ) { + clock_t T_save_start, T_save_stop, T_reload_start, T_reload_stop, T_total_start, T_total_stop, T_main_start, T_main_stop; + T_total_start = clock(); + T_reload_start = clock(); // load and parse settings file(s) + + std::ofstream file_times("output/interesting_times.txt", std::ios::app); + bool XML_status = false; char copy_command [1024]; @@ -122,6 +128,8 @@ int main( int argc, char* argv[] ) /* Microenvironment setup */ setup_microenvironment(); // modify this in the custom code + + bool start_stop = parameters.bools("start_stop"); /* PhysiCell setup */ @@ -132,8 +140,26 @@ int main( int argc, char* argv[] ) /* Users typically start modifying here. START USERMODS */ create_cell_types(); - - setup_tissue(); + + if( start_stop ){ + + + setup_tissue(); + + reset_cell(cell_container->last_cell_cycle_time); + + //exit(-1); + + + reset_global_parameters(cell_container); + + reset_microenv(); + + + + } else{ + setup_tissue(); //death model index = 1 == necrotic...= 0 == apoptotic. + } /* Users typically stop modifying here. END USERMODS */ @@ -183,15 +209,25 @@ int main( int argc, char* argv[] ) report_file.open(filename); // create the data log file report_file<<"simulated time\tnum cells\tnum division\tnum death\twall time"< PhysiCell_globals.next_full_save_time - 0.5 * diffusion_dt ) + if( fabs( PhysiCell_globals.current_time - PhysiCell_globals.next_full_save_time ) < 0.01 * diffusion_dt ) { display_simulation_status( std::cout ); if( PhysiCell_settings.enable_legacy_saves == true ) @@ -201,17 +237,22 @@ int main( int argc, char* argv[] ) if( PhysiCell_settings.enable_full_saves == true ) { + sprintf( filename , "%s/output%08u" , PhysiCell_settings.folder.c_str(), PhysiCell_globals.full_output_index ); save_PhysiCell_to_MultiCellDS_v2( filename , microenvironment , PhysiCell_globals.current_time ); + if ( parameters.bools("auto_stop") ){ + stop = auto_stop(); + } } PhysiCell_globals.full_output_index++; PhysiCell_globals.next_full_save_time += PhysiCell_settings.full_save_interval; + } // save SVG plot if it's time - if( PhysiCell_globals.current_time > PhysiCell_globals.next_SVG_save_time - 0.5 * diffusion_dt ) + if( fabs( PhysiCell_globals.current_time - PhysiCell_globals.next_SVG_save_time ) < 0.01 * diffusion_dt ) { if( PhysiCell_settings.enable_SVG_saves == true ) { @@ -246,6 +287,7 @@ int main( int argc, char* argv[] ) { // reference to the base of a polymorphic object std::cout << e.what(); // information from length_error printed } + T_main_stop = clock(); // save a final simulation snapshot @@ -254,11 +296,26 @@ int main( int argc, char* argv[] ) sprintf( filename , "%s/final.svg" , PhysiCell_settings.folder.c_str() ); SVG_plot( filename , microenvironment, 0.0 , PhysiCell_globals.current_time, cell_coloring_function, ECM_coloring_function); + + // save the necessary files for the start and stop add-on + T_save_start = clock(); + save_cell_microenv_data(cell_container); + std::cout << "cells data saved succesfully" << std::endl; + T_save_stop = clock(); + // timer std::cout << std::endl << "Total simulation runtime: " << std::endl; BioFVM::display_stopwatch_value( std::cout , BioFVM::runtime_stopwatch_value() ); + T_total_stop = clock(); + double T_save, T_reload, T_total, T_main; + T_save = (double)(T_save_stop - T_save_start)/CLOCKS_PER_SEC; + T_reload = (double)(T_reload_stop - T_reload_start)/CLOCKS_PER_SEC; + T_total = (double)(T_total_stop - T_total_start)/CLOCKS_PER_SEC; + T_main = (double)(T_main_stop - T_main_start)/CLOCKS_PER_SEC; + file_times << T_save << " " << T_reload << " " << T_total << " " << T_main << std::endl; + file_times.close(); return 0; -} +} \ No newline at end of file diff --git a/sample_projects_intracellular/boolean/physiboss_cell_lines/Makefile b/sample_projects_intracellular/boolean/physiboss_cell_lines/Makefile index b6205918e..a260f52f6 100644 --- a/sample_projects_intracellular/boolean/physiboss_cell_lines/Makefile +++ b/sample_projects_intracellular/boolean/physiboss_cell_lines/Makefile @@ -95,8 +95,10 @@ PhysiCell_custom_module_OBJECTS := custom.o pugixml_OBJECTS := pugixml.o +start_and_stop_OBJECTS := start_and_stop.o + PhysiCell_OBJECTS := $(BioFVM_OBJECTS) $(pugixml_OBJECTS) $(PhysiCell_core_OBJECTS) $(PhysiCell_module_OBJECTS) -ALL_OBJECTS := $(PhysiCell_OBJECTS) $(PhysiCell_custom_module_OBJECTS) $(PhysiBoSS_OBJECTS) +ALL_OBJECTS := $(PhysiCell_OBJECTS) $(PhysiCell_custom_module_OBJECTS) $(PhysiBoSS_OBJECTS) $(start_and_stop_OBJECTS) # compile the project @@ -202,6 +204,10 @@ PhysiCell_basic_signaling.o: ./core/PhysiCell_basic_signaling.cpp PhysiCell_geometry.o: ./modules/PhysiCell_geometry.cpp $(COMPILE_COMMAND) -c ./modules/PhysiCell_geometry.cpp +# start_and_stop +start_and_stop.o: ./addons/start_and_stop/start_and_stop.cpp + $(COMPILE_COMMAND) $(INC) -c ./addons/start_and_stop/start_and_stop.cpp + # user-defined PhysiCell modules Compile_MaBoSS: ./addons/PhysiBoSS/MaBoSS/engine/src/BooleanNetwork.h diff --git a/sample_projects_intracellular/boolean/physiboss_cell_lines/config/PhysiCell_settings.xml b/sample_projects_intracellular/boolean/physiboss_cell_lines/config/PhysiCell_settings.xml index b41a3fc53..e15ae442c 100644 --- a/sample_projects_intracellular/boolean/physiboss_cell_lines/config/PhysiCell_settings.xml +++ b/sample_projects_intracellular/boolean/physiboss_cell_lines/config/PhysiCell_settings.xml @@ -1062,6 +1062,7 @@ + false 0 C diff --git a/sample_projects_intracellular/boolean/physiboss_cell_lines/custom_modules/custom.h b/sample_projects_intracellular/boolean/physiboss_cell_lines/custom_modules/custom.h index 44cd833d9..006c3b409 100644 --- a/sample_projects_intracellular/boolean/physiboss_cell_lines/custom_modules/custom.h +++ b/sample_projects_intracellular/boolean/physiboss_cell_lines/custom_modules/custom.h @@ -70,6 +70,8 @@ #include "../core/PhysiCell.h" #include "../modules/PhysiCell_standard_modules.h" +#include "../addons/start_and_stop/start_and_stop.h" + using namespace BioFVM; using namespace PhysiCell; diff --git a/sample_projects_intracellular/boolean/physiboss_cell_lines/main.cpp b/sample_projects_intracellular/boolean/physiboss_cell_lines/main.cpp index 002ed6092..8a851c7fe 100644 --- a/sample_projects_intracellular/boolean/physiboss_cell_lines/main.cpp +++ b/sample_projects_intracellular/boolean/physiboss_cell_lines/main.cpp @@ -118,6 +118,13 @@ int main( int argc, char* argv[] ) setup_microenvironment(); // modify this in the custom code + bool start_stop = parameters.bools("start_stop"); + if( start_stop ){ + + // reset microenvironment and cells as they were in the previous simulation + reset_microenv(); + } + /* PhysiCell setup */ // set mechanics voxel size, and match the data structure to BioFVM @@ -128,7 +135,19 @@ int main( int argc, char* argv[] ) create_cell_types(); - setup_tissue(); + if( start_stop ){ + + // reset cells as they were in the previous simulation + setup_tissue(); + + reset_cell(cell_container->last_cell_cycle_time); + //exit(-1); + + reset_global_parameters(cell_container); + + } else{ + setup_tissue(); //death model index = 1 == necrotic...= 0 == apoptotic. + } /* Users typically stop modifying here. END USERMODS */ @@ -182,6 +201,10 @@ int main( int argc, char* argv[] ) report_file<<"simulated time\tnum cells\tnum division\tnum death\twall time"< 64, change lib path +ifeq ($(shell expr $(MABOSS_MAX_NODES) '>' 64), 1) +LIB := -L$(CUR_DIR)/$(MABOSS_DIR)/lib -lMaBoSS_$(MABOSS_MAX_NODES)n-static $(LDL_FLAG) +endif + +ARCH := native # best auto-tuning +# ARCH := core2 # a reasonably safe default for most CPUs since 2007 +# ARCH := corei7 +# ARCH := corei7-avx # earlier i7 +# ARCH := core-avx-i # i7 ivy bridge or newer +# ARCH := core-avx2 # i7 with Haswell or newer +# ARCH := nehalem +# ARCH := westmere +# ARCH := sandybridge # circa 2011 +# ARCH := ivybridge # circa 2012 +# ARCH := haswell # circa 2013 +# ARCH := broadwell # circa 2014 +# ARCH := skylake # circa 2015 +# ARCH := bonnell +# ARCH := silvermont +# ARCH := skylake-avx512 +# ARCH := nocona #64-bit pentium 4 or later + +# CFLAGS := -march=$(ARCH) -Ofast -s -fomit-frame-pointer -mfpmath=both -fopenmp -m64 -std=c++11 +CFLAGS := -march=$(ARCH) -O3 -fomit-frame-pointer -mfpmath=both -fopenmp -m64 -std=c++11 + +ifeq ($(OS),Windows_NT) +else + UNAME_S := $(shell uname -s) + ifeq ($(UNAME_S),Darwin) + UNAME_P := $(shell uname -p) + var := $(shell which $(CC) | xargs file) + ifeq ($(lastword $(var)),arm64) + CFLAGS := -march=$(ARCH) -O3 -fomit-frame-pointer -fopenmp -m64 -std=c++11 + endif + endif +endif + +COMPILE_COMMAND := $(CC) $(CFLAGS) + +BioFVM_OBJECTS := BioFVM_vector.o BioFVM_mesh.o BioFVM_microenvironment.o BioFVM_solvers.o BioFVM_matlab.o \ +BioFVM_utilities.o BioFVM_basic_agent.o BioFVM_MultiCellDS.o BioFVM_agent_container.o + +PhysiCell_core_OBJECTS := PhysiCell_phenotype.o PhysiCell_cell_container.o PhysiCell_standard_models.o \ +PhysiCell_cell.o PhysiCell_custom.o PhysiCell_utilities.o PhysiCell_constants.o PhysiCell_basic_signaling.o \ +PhysiCell_signal_behavior.o PhysiCell_rules.o + +PhysiCell_module_OBJECTS := PhysiCell_SVG.o PhysiCell_pathology.o PhysiCell_MultiCellDS.o PhysiCell_various_outputs.o \ +PhysiCell_pugixml.o PhysiCell_settings.o PhysiCell_geometry.o + +# put your custom objects here (they should be in the custom_modules directory) +MaBoSS := ./addons/PhysiBoSS/MaBoSS/engine/src/BooleanNetwork.h + +PhysiBoSS_OBJECTS := maboss_network.o maboss_intracellular.o + +PhysiCell_custom_module_OBJECTS := custom.o submodel_data_structures.o tnf_receptor_dynamics.o tnf_boolean_model_interface.o + +pugixml_OBJECTS := pugixml.o + +PhysiCell_OBJECTS := $(BioFVM_OBJECTS) $(pugixml_OBJECTS) $(PhysiCell_core_OBJECTS) $(PhysiCell_module_OBJECTS) + +start_and_stop_OBJECTS := start_and_stop.o +ALL_OBJECTS := $(PhysiCell_OBJECTS) $(PhysiCell_custom_module_OBJECTS) $(PhysiBoSS_OBJECTS) $(PhysiBoSS_module_OBJECTS) $(start_and_stop_OBJECTS) + +# compile the project + +all: main.cpp $(ALL_OBJECTS) $(MaBoSS) + $(COMPILE_COMMAND) $(INC) -o $(PROGRAM_NAME) $(ALL_OBJECTS) main.cpp $(LIB) + @echo "" + @echo "check for $(PROGRAM_NAME)" + +# PhysiCell core components + +PhysiCell_phenotype.o: ./core/PhysiCell_phenotype.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_phenotype.cpp + +PhysiCell_digital_cell_line.o: ./core/PhysiCell_digital_cell_line.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_digital_cell_line.cpp + +PhysiCell_cell.o: ./core/PhysiCell_cell.cpp $(MaBoSS) + $(COMPILE_COMMAND) $(INC) -c ./core/PhysiCell_cell.cpp + +PhysiCell_cell_container.o: ./core/PhysiCell_cell_container.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_cell_container.cpp + +PhysiCell_standard_models.o: ./core/PhysiCell_standard_models.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_standard_models.cpp + +PhysiCell_utilities.o: ./core/PhysiCell_utilities.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_utilities.cpp + +PhysiCell_custom.o: ./core/PhysiCell_custom.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_custom.cpp + +PhysiCell_constants.o: ./core/PhysiCell_constants.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_constants.cpp + +PhysiCell_signal_behavior.o: ./core/PhysiCell_signal_behavior.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_signal_behavior.cpp + +PhysiCell_rules.o: ./core/PhysiCell_rules.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_rules.cpp + +# BioFVM core components (needed by PhysiCell) + +BioFVM_vector.o: ./BioFVM/BioFVM_vector.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_vector.cpp + +BioFVM_agent_container.o: ./BioFVM/BioFVM_agent_container.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_agent_container.cpp + +BioFVM_mesh.o: ./BioFVM/BioFVM_mesh.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_mesh.cpp + +BioFVM_microenvironment.o: ./BioFVM/BioFVM_microenvironment.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_microenvironment.cpp + +BioFVM_solvers.o: ./BioFVM/BioFVM_solvers.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_solvers.cpp + +BioFVM_utilities.o: ./BioFVM/BioFVM_utilities.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_utilities.cpp + +BioFVM_basic_agent.o: ./BioFVM/BioFVM_basic_agent.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_basic_agent.cpp + +BioFVM_matlab.o: ./BioFVM/BioFVM_matlab.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_matlab.cpp + +BioFVM_MultiCellDS.o: ./BioFVM/BioFVM_MultiCellDS.cpp + $(COMPILE_COMMAND) -c ./BioFVM/BioFVM_MultiCellDS.cpp + +pugixml.o: ./BioFVM/pugixml.cpp + $(COMPILE_COMMAND) -c ./BioFVM/pugixml.cpp + +# standard PhysiCell modules + +PhysiCell_SVG.o: ./modules/PhysiCell_SVG.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_SVG.cpp + +PhysiCell_pathology.o: ./modules/PhysiCell_pathology.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_pathology.cpp + +PhysiCell_MultiCellDS.o: ./modules/PhysiCell_MultiCellDS.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_MultiCellDS.cpp + +PhysiCell_various_outputs.o: ./modules/PhysiCell_various_outputs.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_various_outputs.cpp + + +PhysiCell_pugixml.o: ./modules/PhysiCell_pugixml.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_pugixml.cpp + +PhysiCell_settings.o: ./modules/PhysiCell_settings.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_settings.cpp + +PhysiCell_basic_signaling.o: ./core/PhysiCell_basic_signaling.cpp + $(COMPILE_COMMAND) -c ./core/PhysiCell_basic_signaling.cpp + +PhysiCell_geometry.o: ./modules/PhysiCell_geometry.cpp + $(COMPILE_COMMAND) -c ./modules/PhysiCell_geometry.cpp + +# start_and_stop +start_and_stop.o: ./addons/start_and_stop/start_and_stop.cpp + $(COMPILE_COMMAND) $(INC) -c ./addons/start_and_stop/start_and_stop.cpp + +# user-defined PhysiCell modules + +Compile_MaBoSS: ./addons/PhysiBoSS/MaBoSS/engine/src/BooleanNetwork.h + cd ./addons/PhysiBoSS/MaBoSS/engine/src;make CXX=$(CC) install_alib;make clean; cd ../../../../.. + +$(MaBoSS): +ifeq ($(OS), Windows_NT) + python addons/PhysiBoSS/setup_libmaboss.py +else + python3 addons/PhysiBoSS/setup_libmaboss.py +endif + +maboss_network.o: ./addons/PhysiBoSS/src/maboss_network.cpp $(MaBoSS) + $(COMPILE_COMMAND) $(INC) -c ./addons/PhysiBoSS/src/maboss_network.cpp + +maboss_intracellular.o: ./addons/PhysiBoSS/src/maboss_intracellular.cpp $(MaBoSS) + $(COMPILE_COMMAND) $(INC) -c ./addons/PhysiBoSS/src/maboss_intracellular.cpp + +custom.o: ./custom_modules/custom.cpp $(MaBoSS) + $(COMPILE_COMMAND) $(INC) -c ./custom_modules/custom.cpp + +submodel_data_structures.o: ./custom_modules/submodel_data_structures.cpp + $(COMPILE_COMMAND) $(INC) -c ./custom_modules/submodel_data_structures.cpp + +tnf_receptor_dynamics.o: ./custom_modules/tnf_receptor_dynamics.cpp + $(COMPILE_COMMAND) $(INC) -c ./custom_modules/tnf_receptor_dynamics.cpp + +tnf_boolean_model_interface.o: ./custom_modules/tnf_boolean_model_interface.cpp $(MaBoSS) + $(COMPILE_COMMAND) $(INC) -c ./custom_modules/tnf_boolean_model_interface.cpp + @echo "compiling" +# cleanup + +reset: + rm -f *.cpp + cp ./sample_projects/Makefile-default Makefile + rm -f ./custom_modules/* + touch ./custom_modules/empty.txt + rm ALL_CITATIONS.txt + rm -f ./config/PhysiCell_settings.xml + rm -f ./config/PhysiCell_settings_2D.xml + rm -f ./config/PhysiCell_settings_3D.xml + rm -f ./config/init.tsv + rm -f ./config/TNF_conf.cfg + rm -f ./config/TNF_nodes.bnd + rm -fr ./scripts + cp ./config/PhysiCell_settings-backup.xml ./config/PhysiCell_settings.xml + +MaBoSS-clean: + rm -fr addons/PhysiBoSS/MaBoSS + +clean: MaBoSS-clean + rm -f *.o + rm -f $(PROGRAM_NAME)* + +data-cleanup: + rm -f *.mat + rm -f *.xml + rm -f *.svg + rm -f ./output/* + touch ./output/empty.txt + +# archival + +checkpoint: + zip -r $$(date +%b_%d_%Y_%H%M).zip Makefile *.cpp *.h config/*.xml custom_modules/* + +zip: + zip -r latest.zip Makefile* *.cpp *.h BioFVM/* config/* core/* custom_modules/* matlab/* modules/* sample_projects/* + cp latest.zip $$(date +%b_%d_%Y_%H%M).zip + cp latest.zip VERSION_$(VERSION).zip + mv *.zip archives/ + +tar: + tar --ignore-failed-read -czf latest.tar Makefile* *.cpp *.h BioFVM/* config/* core/* custom_modules/* matlab/* modules/* sample_projects/* + cp latest.tar $$(date +%b_%d_%Y_%H%M).tar + cp latest.tar VERSION_$(VERSION).tar + mv *.tar archives/ + +unzip: + cp ./archives/latest.zip . + unzip latest.zip + +untar: + cp ./archives/latest.tar . + tar -xzf latest.tar diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings.xml b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings.xml new file mode 100644 index 000000000..a32e7c30d --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings.xml @@ -0,0 +1,225 @@ + + + -500 + 500 + -500 + 500 + -10 + 10 + 20 + 20 + 20 + true + + + + 180 + min + micron + 0.02 + 0.1 + 2 + + + + 8 + + + + output + + 30 + true + + + 30 + true + + + false + + + + + + + 100000.0 + .1 + + 38.0 + 38.0 + + + + + 1200.0 + .0275 + + 0.0 + 0.0 + + + + false + true + + + ./config/initial.mat + + + + ./config/dirichlet.mat + + + + + + + + + + + + 0.0011 + + + + + + 5.31667e-05 + + 516 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + + 0.0 + + 0 + 86400 + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + + + 2494 + 0.75 + 540 + + 0.05 + 0.0045 + 0.0055 + + 0 + 0 + + 2.0 + + + + + 0 + 38 + 10 + 0 + + + 0 + 0 + 0 + 0 + + + + + + ./config/TNF_nodes.bnd + ./config/TNF_conf.cfg + + 10.0 + 0.01 + + + + + + 0.243 + 0.128 + 0.293 + + 0.50 + 1.0 + 0 + + + 0 + 0 + 0 + 0 + + + + 0.0035 + 0.0035 + 1.0 + + + + + + 0 + false + ./start_and_stop_saving_files/initial.tsv + 100 + + + 120 + + false + false + false + + true + 150 + 10 + 0.005 + 999999 + + + + 0.0 + 0.05 + 0.95 + + + 0.0 + 0.05 + 0.95 + + + 0.0 + 0.05 + 0.95 + + + + + + + \ No newline at end of file diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings_2D.xml b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings_2D.xml new file mode 100644 index 000000000..a6764a6f4 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings_2D.xml @@ -0,0 +1,286 @@ + + + + + + + + -500 + 500 + -500 + 500 + -10 + 10 + 20 + 20 + 20 + true + + + + 1440 + min + micron + 0.02 + 0.1 + 2 + + + + 8 + + + + output + + 30 + true + + + 30 + true + + + false + + + + + + + 100000.0 + .1 + + 38.0 + 38.0 + + + + + 1200.0 + .0275 + + 0.0 + 0.0 + + + + false + true + + + ./config/initial.mat + + + + ./config/dirichlet.mat + + + + + + + + + + + + 0.00075 + + + + + + 5.31667e-05 + + 516 + + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + + 0.0 + + 0 + 86400 + + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + + + 2494 + 0.75 + 540 + + 0.05 + 0.0045 + 0.0055 + + 0 + 0 + + 2.0 + + + + + + 0 + 38 + 10 + 0 + + + 0 + 0 + 0 + 0 + + + + + + ./config/TNF_nodes.bnd + ./config/TNF_conf.cfg + + 10 + 1 + + + + + + 0.243 + 0.128 + 0.293 + + 0.50 + 1.0 + 0 + + + 0 + 0 + 0 + 0 + + + 0.0027 + 0.0055 + + + + + + 0 + false + ./config/init.tsv + 100.0 + + 0.0 + 0.05 + 0.95 + + 0.0 + 0.05 + 0.95 + + 0.0 + 0.05 + 0.95 + + 150.0 + 10.0 + 0.0015 + 999999 + 120 + + + + diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings_3D.xml b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings_3D.xml new file mode 100644 index 000000000..dab033eb4 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/PhysiCell_settings_3D.xml @@ -0,0 +1,287 @@ + + + + + + + + -200 + 200 + -200 + 200 + -200 + 200 + 20 + 20 + 20 + false + + + + 2160 + min + micron + 0.02 + 0.1 + 2 + + + + 8 + + + + output + + 30 + true + + + 30 + true + + + false + + + + + + + 100000.0 + .1 + + 38.0 + 38.0 + + + + + 1200.0 + .0275 + + 0.0 + 0.0 + + + + false + true + + + ./config/initial.mat + + + + ./config/dirichlet.mat + + + + + + + + + + + + 0.00075 + + + + + + 5.31667e-05 + + 516 + + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + + 0.0 + + 0 + 86400 + + + + 0.05 + 0 + 1.66667e-02 + 5.83333e-03 + 0 + 2.0 + + + + + + 2494 + 0.75 + 540 + + 0.05 + 0.0045 + 0.0055 + + 0 + 0 + + 2.0 + + + + + + 0 + 38 + 10 + 0 + + + 0 + 0 + 0 + 0 + + + + + + ./config/TNF_nodes.bnd + ./config/TNF_conf.cfg + + 10 + 1 + + + + + + 0.243 + 0.128 + 0.293 + + 0.50 + 1.0 + 0 + + + 0 + 0 + 0 + 0 + + + 0.0027 + 0.0055 + 1e-7 + + + + + + 0 + ./config/init.tsv + false + 50.0 + + 0.0 + 0.05 + 0.95 + + 0.0 + 0.05 + 0.95 + + 0.0 + 0.05 + 0.95 + + 150.0 + 10.0 + 0.002 + 999999 + 80 + + + + \ No newline at end of file diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/config/TNF_conf.cfg b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/TNF_conf.cfg new file mode 100644 index 000000000..5e9c9892a --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/TNF_conf.cfg @@ -0,0 +1,101 @@ + +// input +FADD.istate = 0; +TNF.istate = 0; +FASL.istate = 0; + +// output +NonACD.istate = 0; +Apoptosis.istate = 0 ; +Survival.istate = 0 ; + +// specific state +ATP.istate = 1; +cIAP.istate = 1; + +TNFR.istate = 0 ; +DISC_TNF.istate = 0 ; +DISC_FAS.istate = 0 ; +RIP1.istate = 0 ; +RIP1ub.istate = 0 ; +RIP1K.istate = 0 ; +IKK.istate = 0 ; +NFkB.istate = 0 ; +CASP8.istate = 0 ; +BAX.istate = 0 ; +BCL2.istate = 0 ; +ROS.istate = 0 ; +mROS.istate = 0 ; +MPT.istate = 0 ; +MOMP.istate = 0 ; +SMAC.istate = 0 ; +mcIAP.istate = 0 ; +Cyt_c.istate = 0 ; +XIAP.istate = 0 ; +mXIAP.istate = 0 ; +apoptosome.istate = 0 ; +CASP3.istate = 0 ; +cFLIP.istate = 0 ; + +//$ProdTNF_NFkB = 1; +$TransRate = 1/24; +//$DivRate = 1/24; +//$Degr_TNF = 1/6; +//$TNF_induc = 0; + +//[TNF].istate = 1 [1] , 0 [0]; + + +//TNF.is_internal = TRUE ; +//ATP.is_internal = TRUE ; +//FADD.is_internal = TRUE ; +//cIAP.is_internal = TRUE ; +// +//FASL.is_internal = TRUE ; +//TNFR.is_internal = TRUE ; +//DISC_TNF.is_internal = TRUE ; +//DISC_FAS.is_internal = TRUE ; +//RIP1.is_internal = TRUE ; +//RIP1ub.is_internal = TRUE ; +//RIP1K.is_internal = TRUE ; +//IKK.is_internal = TRUE ; +//NFkB.is_internal = TRUE ; +//CASP8.is_internal = TRUE ; +//BAX.is_internal = TRUE ; +//BCL2.is_internal = TRUE ; +//ROS.is_internal = TRUE ; +//MPT.is_internal = TRUE ; +//MOMP.is_internal = TRUE ; +//SMAC.is_internal = TRUE ; +//Cyt_c.is_internal = TRUE ; +//XIAP.is_internal = TRUE ; +//apoptosome.is_internal = TRUE ; +//CASP3.is_internal = TRUE ; +//cFLIP.is_internal = TRUE ; + + +sample_count = 1; +max_time = 1; // 10 min: 1 in maboss = 10 min in Physicell +time_tick = .01; +discrete_time = 0; +use_physrandgen = FALSE; +seed_pseudorandom = 37; +display_traj = FALSE; + +thread_count = 1; + +statdist_traj_count = 1; +statdist_cluster_threshold = 0.8; + +$Low_CASP3 = 0.0; +$High_CASP3 = 0.0; +$High_Cytc = 0.0; +$Low_CASP8 = 0.0; +$High_CASP8 = 0.0; +$High_IKK = 0.0; +$High_cFLIP = 0.0; +$Low_cIAP = 0.0; +$High_mROS = 0.0; +$Low_RIP1 = 0.0; +$High_NFkB = 0.0; +$Low_NFkB = 0.0; diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/config/TNF_nodes.bnd b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/TNF_nodes.bnd new file mode 100644 index 000000000..d91d6f8dc --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/TNF_nodes.bnd @@ -0,0 +1,184 @@ +node FASL +{ + rate_up = 0.0; + rate_down = 0.0; +} +node TNF +{ + logic = TNF; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node TNFR +{ + logic = TNF; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node DISC_TNF +{ + logic = FADD & TNFR; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node DISC_FAS +{ + logic = FASL & FADD; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node FADD +{ + rate_up = 0.0; + rate_down = 0.0; +} +node RIP1 +{ + logic = (DISC_FAS | TNFR) & (!CASP8); + rate_up = ($Low_RIP1 ? 0.0 : ( @logic ? 1.0: 0.0 )); + rate_down = ($Low_RIP1 ? 1E+100 : ( @logic ? 0.0 : 1.0 )); +} +node RIP1ub +{ + logic = cIAP & RIP1; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node RIP1K +{ + logic = RIP1; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node IKK +{ + logic = RIP1ub; + rate_up = ($High_IKK ? 1E+100 : (@logic ? 1.0 : 0.0) ); + rate_down = ($High_IKK ? 0.0 : (@logic ? 0.0 : 1.0) ); +} +node NFkB +{ + logic = IKK & (!CASP3); + rate_up = ($Low_NFkB ? 0.0 : ($High_NFkB ? 1E+100 : ( @logic ? 1.0 : 0.0) )); + rate_down = ($Low_NFkB ? 1E+100 : ($High_NFkB ? 0.0 : (@logic ? 0.0 : 1.0) )); +} +node CASP8 +{ + logic = (DISC_TNF | (DISC_FAS | CASP3) ) & (!cFLIP); + rate_up = ($Low_CASP8 ? 0.0 : ( $High_CASP8 ? 1E+100 : ( @logic ? 1.0: 0.0 )) ); + rate_down = ($Low_CASP8 ? 1E+100 : ( $High_CASP8 ? 0.0 : (@logic ? 0.0 : 1.0 )) ); +} +node BAX +{ + logic = CASP8 & (!BCL2); + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node BCL2 +{ + logic = NFkB; + rate_up = (@logic ? $TransRate : 0.0); + rate_down = @logic ? 0.0 : 1.0; +} +node mROS +{ + logic = (!NFkB); + rate_up = ($High_mROS ? 1E+100 : (@logic ? $TransRate : 0.0) ); + rate_down = ($High_mROS ? 0.0: (@logic ? 0.0 : 1.0) ); +} +node ROS +{ + logic = (mROS) & (MPT | RIP1K ); + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node ATP +{ + logic = !MPT; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node MPT +{ + logic = (!BCL2) & ROS; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node MOMP +{ + logic = BAX | MPT; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node SMAC +{ + logic = MOMP; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node mcIAP +{ + logic = (NFkB); + rate_up = @logic ? $TransRate : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node cIAP +{ + logic = (mcIAP & (!SMAC)); + rate_up = ($Low_cIAP ? 0.0 : ( @logic ? $TransRate: 0.0 )); + rate_down = ($Low_cIAP ? 1E+100 : ( (SMAC) ? 1.0 : 0.0 )); +} +node Cyt_c +{ + logic = MOMP; + rate_up = ($High_Cytc ? 1E+100 : (@logic ? 1.0 : 0.0) ); + rate_down = ($High_Cytc ? 0.0 : (@logic ? 0.0 : 1.0) ); +} +node mXIAP +{ + logic = (NFkB); + rate_up = @logic ? $TransRate : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node XIAP +{ + logic = (!SMAC) & mXIAP; + rate_up = (@logic) ? 1.0 : 0.0; + rate_down = (@logic ? 0.0 : 1.0); +} +node apoptosome +{ + logic = Cyt_c & (ATP & (!XIAP)); + rate_up = (@logic) ? 1.0 : 0.0; + rate_down = (@logic ? 0.0 : 1.0); +} +node CASP3 +{ + logic = apoptosome & (!XIAP); + rate_up = ( $Low_CASP3 ? 0.0 : ($High_CASP3 ? 1E+100 : ((@logic ? 1.0 : 0.0)) ) ); + rate_down = ($Low_CASP3 ? 1E+100 : ($High_CASP3 ? 0.0 : ( @logic ? 0.0 : 1.0) ) ); +} +node cFLIP +{ + logic = NFkB ; + rate_up = ($High_cFLIP ? 1E+100 : ((@logic) ? $TransRate : 0.0) ); + rate_down = ($High_cFLIP ? 0.0 : (@logic ? 0.0 : 1.0) ); +} +node NonACD +{ + logic = !ATP; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node Apoptosis +{ + logic = CASP3; + rate_up = @logic ? 1.0 : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} +node Survival +{ + logic = NFkB; + rate_up = @logic ? $TransRate : 0.0; + rate_down = @logic ? 0.0 : 1.0; +} diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/config/init.tsv b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/init.tsv new file mode 100644 index 000000000..02e46f221 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/config/init.tsv @@ -0,0 +1,126 @@ +x y z +-93.7351 -16.7883 0.0 +-95.7338 2.90696 0.0 +-95.9399 18.5505 0.0 +-79.8482 -55.7662 0.0 +-82.7174 -41.7014 0.0 +-80.7886 -22946 0.0 +-79.0268 -6.05337 0.0 +-82.0067 10.4321 0.0 +-80.6049 25.4132 0.0 +-82.8366 44.0516 0.0 +-78.9526 61124 0.0 +-64102 -65.5156 0.0 +-64.2225 -46.9066 0.0 +-68468 -33.3416 0.0 +-68.4907 -14.3586 0.0 +-67762 -0.0687358 0.0 +-64.3596 17.0034 0.0 +-65.1953 36802 0.0 +-65.7418 51.8249 0.0 +-64.8725 70.5213 0.0 +-49.7993 -75.1665 0.0 +-50.1662 -59.6431 0.0 +-51.6159 -42.4253 0.0 +-52.2097 -24.4339 0.0 +-49.6998 -5.49337 0.0 +-52.8067 10.2768 0.0 +-53.0525 25.5959 0.0 +-50.1145 42.8981 0.0 +-50.8895 63.7084 0.0 +-51.5362 77.8187 0.0 +-35.1206 -84.1124 0.0 +-39.0413 -67.4176 0.0 +-36.7631 -50.2966 0.0 +-35.5977 -32.8759 0.0 +-38.2835 -13.9819 0.0 +-36793 3.2073 0.0 +-38.2603 20.0752 0.0 +-38.5707 36.7987 0.0 +-35.5583 50.5694 0.0 +-35.4349 71.8035 0.0 +-34977 87.9757 0.0 +-20.4709 -93.9292 0.0 +-20396 -76.3795 0.0 +-23.7124 -57.6671 0.0 +-23908 -40.2429 0.0 +-21.0551 -22.9972 0.0 +-20.3947 -5.09517 0.0 +-20212 9.72254 0.0 +-24.1934 27.5216 0.0 +-19.9637 47.0325 0.0 +-22.1828 63.9466 0.0 +-22.5229 76.7938 0.0 +-23.5922 96.0568 0.0 +-7.4651 -98934 0.0 +-5.63692 -80.5666 0.0 +-9.19887 -66.9041 0.0 +-7.83299 -49.2237 0.0 +-7.58047 -29.8289 0.0 +-9.17806 -14.3046 0.0 +-8.47321 1.79866 0.0 +-6.22558 17232 0.0 +-6.1979 34592 0.0 +-6.06438 51.8374 0.0 +-7.99579 71.8464 0.0 +-9.34279 88.3144 0.0 +8.09689 -92.1459 0.0 +8.56241 -76.1541 0.0 +5.75799 -58.7743 0.0 +6.80561 -42.8107 0.0 +7.29498 -24.8781 0.0 +9.49315 -8.67633 0.0 +5.48934 10.7607 0.0 +7.48964 25.2213 0.0 +8.37825 44.0642 0.0 +6.63344 60.3534 0.0 +8.82467 76.4799 0.0 +5.50615 97.7667 0.0 +20.0973 -85.3025 0.0 +23.6045 -66.3319 0.0 +21.6682 -31.0106 0.0 +20.8476 -14.3643 0.0 +23354 4.3938 0.0 +23.2571 21.0108 0.0 +21.9115 37.8318 0.0 +23.0517 50.9854 0.0 +20875 68.0303 0.0 +37.9323 -76.0558 0.0 +38.8808 -59.6001 0.0 +34.7439 -39.8988 0.0 +34.7224 -26.0126 0.0 +38066 -8.93449 0.0 +36.8513 9.09404 0.0 +36.4099 27.3056 0.0 +35.4412 44.6587 0.0 +36.1567 62.8726 0.0 +37.3943 79938 0.0 +50.5785 -83.4395 0.0 +51.7827 -64.4546 0.0 +50186 -48.5523 0.0 +53.3226 -33.7905 0.0 +51.1424 -15.7653 0.0 +52.0446 1.17038 0.0 +52.8293 16.8413 0.0 +49.7305 34.3809 0.0 +51.0108 55.0179 0.0 +49.4942 67.6309 0.0 +68.0616 -56.1252 0.0 +67.5741 -38.8549 0.0 +64.4535 -24.3364 0.0 +67.9933 -6.46241 0.0 +65.0907 12.6911 0.0 +65.1839 25.4711 0.0 +66.8215 44.4639 0.0 +67.4613 60.1319 0.0 +82.3886 -46.9536 0.0 +82.0147 -33.4872 0.0 +81.3058 -15.8358 0.0 +79.4432 -0.165153 0.0 +81181 18.2513 0.0 +82.8324 35.9978 0.0 +82.6803 55.3748 0.0 +95.6087 -23.2664 0.0 +96.4069 -7.70411 0.0 +94.1872 12.5324 0.0 +93.7707 26.4937 0.0 diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/custom.cpp b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/custom.cpp new file mode 100644 index 000000000..287d9700e --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/custom.cpp @@ -0,0 +1,545 @@ +/* +############################################################################### +# 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 x.y.z) [1]. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# See VERSION.txt or call get_PhysiCell_version() to get the current version # +# x.y.z. Call display_citations() to get detailed information on all cite-# +# able software used in your PhysiCell application. # +# # +# Because PhysiCell extensively uses BioFVM, we suggest you also cite BioFVM # +# as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1], # +# with BioFVM [2] to solve the transport equations. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# [2] A Ghaffarizadeh, SH Friedman, and P Macklin, BioFVM: an efficient para- # +# llelized diffusive transport solver for 3-D biological simulations, # +# Bioinformatics 32(8): 1256-8, 2016. DOI: 10.1093/bioinformatics/btv730 # +# # +############################################################################### +# # +# BSD 3-Clause License (see https://opensource.org/licenses/BSD-3-Clause) # +# # +# Copyright (c) 2015-2018, Paul Macklin and the PhysiCell Project # +# All rights reserved. # +# # +# Redistribution and use in source and binary forms, with or without # +# modification, are permitted provided that the following conditions are met: # +# # +# 1. Redistributions of source code must retain the above copyright notice, # +# this list of conditions and the following disclaimer. # +# # +# 2. Redistributions in binary form must reproduce the above copyright # +# notice, this list of conditions and the following disclaimer in the # +# documentation and/or other materials provided with the distribution. # +# # +# 3. Neither the name of the copyright holder nor the names of its # +# contributors may be used to endorse or promote products derived from this # +# software without specific prior written permission. # +# # +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE # +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR # +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF # +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS # +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN # +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) # +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # +# POSSIBILITY OF SUCH DAMAGE. # +# # +############################################################################### +*/ + +#include "./custom.h" +#include +#include +#include +#include +#include +#include +#include +#include + + +// declare cell definitions here +void create_cell_types(void) +{ + SeedRandom(); + + initialize_default_cell_definition(); + + /* This parses the cell definitions in the XML config file. */ + initialize_cell_definitions_from_pugixml(); + + // This sets the pre and post intracellular update functions + cell_defaults.functions.pre_update_intracellular = update_boolean_model_inputs; + cell_defaults.functions.post_update_intracellular = update_behaviors; + cell_defaults.functions.update_phenotype = NULL; + + // This initializes the the TNF receptor model + tnf_receptor_model_setup(); + tnf_boolean_model_interface_setup(); + submodel_registry.display(std::cout); + + // 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"]; + cell_defaults.custom_data["bound_external_TNFR"] = 0; + cell_defaults.custom_data["bound_internal_TNFR"] = 0; + + build_cell_definitions_maps(); + + setup_signal_behavior_dictionaries(); + + display_cell_definitions(std::cout); + + return; +} + + +void setup_microenvironment(void) +{ + initialize_microenvironment(); + return; +} + +void setup_tissue(void) +{ + std::vector> positions; + + if ( parameters.bools("read_init") ) + { + std::string csv_fname = parameters.strings("init_cells_filename"); + positions = read_cells_positions(csv_fname, ',', true); + + } + else + { + double cell_radius = cell_defaults.phenotype.geometry.radius; + double tumor_radius = parameters.doubles("tumor_radius"); + if (default_microenvironment_options.simulate_2D == true) + positions = create_cell_disc_positions(cell_radius,tumor_radius); + else + positions = create_cell_sphere_positions(cell_radius,tumor_radius); + + } + + Cell* pCell = NULL; + for (int i = 0; i < positions.size(); i++) + { + pCell = create_cell(get_cell_definition("default")); + pCell->assign_position(positions[i]); + + static int idx_bind_rate = pCell->custom_data.find_variable_index( "TNFR_binding_rate" ); + static float mean_bind_rate = pCell->custom_data[idx_bind_rate]; + static float std_bind_rate = parameters.doubles("TNFR_binding_rate_std"); + static float min_bind_rate = parameters.doubles("TNFR_binding_rate_min"); + static float max_bind_rate = parameters.doubles("TNFR_binding_rate_max"); + + if(std_bind_rate > 0 ) + { + pCell->custom_data[idx_bind_rate] = NormalRandom(mean_bind_rate, std_bind_rate); + if (pCell->custom_data[idx_bind_rate] < min_bind_rate) + { pCell->custom_data[idx_bind_rate] = min_bind_rate; } + if (pCell->custom_data[idx_bind_rate] > max_bind_rate) + { pCell->custom_data[idx_bind_rate] = max_bind_rate; } + } + + + static int idx_endo_rate = pCell->custom_data.find_variable_index( "TNFR_endocytosis_rate" ); + static float mean_endo_rate = pCell->custom_data[idx_endo_rate]; + static float std_endo_rate = parameters.doubles("TNFR_endocytosis_rate_std"); + static float min_endo_rate = parameters.doubles("TNFR_endocytosis_rate_min"); + static float max_endo_rate = parameters.doubles("TNFR_endocytosis_rate_max"); + + if(std_endo_rate > 0) + { + pCell->custom_data[idx_endo_rate] = NormalRandom(mean_endo_rate, std_endo_rate); + if (pCell->custom_data[idx_endo_rate] < min_endo_rate) + { pCell->custom_data[idx_endo_rate] = min_endo_rate; } + if (pCell->custom_data[idx_endo_rate] > max_endo_rate) + { pCell->custom_data[idx_endo_rate] = max_endo_rate; } + } + + static int idx_recycle_rate = pCell->custom_data.find_variable_index( "TNFR_recycling_rate" ); + static float mean_recycle_rate = pCell->custom_data[idx_recycle_rate]; + static float std_recycle_rate = parameters.doubles("TNFR_recycling_rate_std"); + static float min_recycle_rate = parameters.doubles("TNFR_recycling_rate_min"); + static float max_recycle_rate = parameters.doubles("TNFR_recycling_rate_max"); + + if(std_recycle_rate > 0) + { + pCell->custom_data[idx_recycle_rate] = NormalRandom(mean_recycle_rate, std_recycle_rate); + if (pCell->custom_data[idx_recycle_rate] < min_recycle_rate) + { pCell->custom_data[idx_recycle_rate] = min_recycle_rate; } + if (pCell->custom_data[idx_recycle_rate] > max_recycle_rate) + { pCell->custom_data[idx_recycle_rate] = max_recycle_rate; } + } + + update_monitor_variables(pCell); + } + + return; +} + +void update_variables_monitor() +{ + for (int i = 0; i < (*all_cells).size(); i++) + { + // Access the current cell + Cell *pCell = (*all_cells)[i]; + + update_monitor_variables(pCell); + } + + +} + +std::vector> read_cells_positions(std::string filename, char delimiter, bool header) +{ + // File pointer + std::fstream fin; + std::vector> positions; + + // Open an existing file + fin.open(filename, std::ios::in); + + // Read the Data from the file + // as String Vector + std::vector row; + std::string line, word; + + if (header) + { getline(fin, line); } + + do + { + row.clear(); + + // read an entire row and + // store it in a string variable 'line' + getline(fin, line); + + // used for breaking words + std::stringstream s(line); + + while (getline(s, word, delimiter)) + { + row.push_back(word); + } + + std::vector tempPoint(3,0.0); + tempPoint[0]= std::stof(row[0]); + tempPoint[1]= std::stof(row[1]); + tempPoint[2]= std::stof(row[2]); + + positions.push_back(tempPoint); + } while (!fin.eof()); + + return positions; +} + +std::vector> create_cell_sphere_positions(double cell_radius, double sphere_radius) +{ + std::vector> cells; + int xc=0,yc=0,zc=0; + double x_spacing= cell_radius*sqrt(3); + double y_spacing= cell_radius*2; + double z_spacing= cell_radius*sqrt(3); + + std::vector tempPoint(3,0.0); + // std::vector cylinder_center(3,0.0); + + for(double z=-sphere_radius;z> create_cell_disc_positions(double cell_radius, double disc_radius) +{ + double cell_spacing = 0.95 * 2.0 * cell_radius; + + double x = 0.0; + double y = 0.0; + double x_outer = 0.0; + + std::vector> positions; + std::vector tempPoint(3,0.0); + + int n = 0; + while( y < disc_radius ) + { + x = 0.0; + if( n % 2 == 1 ) + { x = 0.5 * cell_spacing; } + x_outer = sqrt( disc_radius*disc_radius - y*y ); + + while( x < x_outer ) + { + tempPoint[0]= x; tempPoint[1]= y; tempPoint[2]= 0.0; + positions.push_back(tempPoint); + if( fabs( y ) > 0.01 ) + { + tempPoint[0]= x; tempPoint[1]= -y; tempPoint[2]= 0.0; + positions.push_back(tempPoint); + } + if( fabs( x ) > 0.01 ) + { + tempPoint[0]= -x; tempPoint[1]= y; tempPoint[2]= 0.0; + positions.push_back(tempPoint); + if( fabs( y ) > 0.01 ) + { + tempPoint[0]= -x; tempPoint[1]= -y; tempPoint[2]= 0.0; + positions.push_back(tempPoint); + } + } + x += cell_spacing; + } + y += cell_spacing * sqrt(3.0)/2.0; + n++; + } + return positions; +} + +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++) + { + auto current_voxel = microenvironment.voxels(n); + std::vector cent = {current_voxel.center[0], current_voxel.center[1], current_voxel.center[2]}; + + if ((membrane_lenght - norm(cent)) <= 0) + microenvironment.density_vector(n)[density_index] = concentration; + } +} + +void remove_density(int density_index) +{ + for (int n = 0; n < microenvironment.number_of_voxels(); n++) + microenvironment.density_vector(n)[density_index] = 0; +} + +std::vector my_coloring_function(Cell *pCell) +{ + // start with live coloring + std::vector output = false_cell_coloring_live_dead(pCell); + + // dead cells + if (pCell->phenotype.death.dead == false) + { + static int nR_EB = pCell->custom_data.find_variable_index("bound external TNFR"); + float activation_threshold = pCell->custom_data.find_variable_index("TNFR activation threshold"); + + int bounded_tnf = (int)round((pCell->custom_data[nR_EB] / activation_threshold) * 255.0); + if (bounded_tnf > 0) + { + char szTempString[128]; + sprintf(szTempString, "rgb(%u,%u,%u)", bounded_tnf, bounded_tnf, 255 - bounded_tnf); + output[0].assign("black"); + output[1].assign(szTempString); + output[2].assign("black"); + output[3].assign(szTempString); + } + } + + return output; +} + + +double total_live_cell_count() +{ + double out = 0.0; + + for( int i=0; i < (*all_cells).size() ; i++ ) + { + if( (*all_cells)[i]->phenotype.death.dead == false && (*all_cells)[i]->type == 0 ) + { out += 1.0; } + } + + return out; +} + +double total_dead_cell_count() +{ + double out = 0.0; + + for( int i=0; i < (*all_cells).size() ; i++ ) + { + if( (*all_cells)[i]->phenotype.death.dead == true && (*all_cells)[i]->phenotype.death.current_death_model_index == 0 ) + { out += 1.0; } + } + + return out; +} + +double total_necrosis_cell_count() +{ + double out = 0.0; + + for( int i=0; i < (*all_cells).size() ; i++ ) + { + if( (*all_cells)[i]->phenotype.death.dead == true && (*all_cells)[i]->phenotype.death.current_death_model_index == 1 ) + { out += 1.0; } + } + + return out; +} + +using namespace std; + +vector vector_alives; + +bool auto_stop_resistance(int alive_cells, int resistant_cells) { + // Define stable states and nodes as vectors instead of unordered_sets + + vector_alives.push_back(alive_cells); + std::cout << "Steps: " << vector_alives.size() << std::endl; + + double threshold = 0.8; + double percentage_of_resistant = static_cast(resistant_cells) / alive_cells; + std::cout << "Number of resistant: " << resistant_cells << std::endl; + + bool stop; + bool condition = false; + if (vector_alives.size() >= 4){ + condition = percentage_of_resistant >= threshold; + } + + if (condition) { + stop = true; + } else { + stop = false; + } + return stop; +} + +bool auto_stop_alive(int alive_cells) { + //concatenate the number of alive cells to the vector + vector_alives.push_back(alive_cells); + std::cout << "Steps: " << vector_alives.size() << std::endl; + + bool condition = false; + // check the number of elements inside the vector to decide if process it and compute the derivative + if (vector_alives.size() >= 8) { + std::vector derivative; + + // compute the derivative only for the last three steps + for (size_t i = vector_alives.size() - 4; i < vector_alives.size(); ++i) { + double slope = vector_alives[i] - vector_alives[i - 1]; + derivative.push_back(slope); + } + + condition = true; + for (double slope : derivative) { + // if the slope is less than or equal to zero, set condition to false + if (slope < 0) { + condition = false; + break; + } + } + + } else { + condition = false; + } + + bool stop; + + if (condition) { + stop = true; + } else { + stop = false; + } + return stop; +} + +int save_resistant_cells(ofstream& file_resistant){ + //count the number of resistant cells for the plots + vector> stable_states = { + {"TNF", "TNFR", "RIP1", "RIP1ub", "RIP1K", "IKK", "NFkB", "BCL2", "ATP", "cIAP", "XIAP", "cFLIP", "Survival"}, + {"FASL", "TNF", "TNFR", "RIP1", "RIP1ub", "RIP1K", "IKK", "NFkB", "BCL2", "ATP", "cIAP", "XIAP", "cFLIP", "Survival"}, + {"TNF", "TNFR", "DISC_TNF", "FADD", "RIP1", "RIP1ub", "RIP1K", "IKK", "NFkB", "BCL2", "ATP", "cIAP", "XIAP", "cFLIP", "Survival"}, + {"FASL", "DISC_FAS", "FADD", "RIP1", "RIP1ub", "RIP1K", "IKK", "NFkB", "BCL2", "ATP", "cIAP", "XIAP", "cFLIP", "Survival"}, + {"FASL", "TNF", "TNFR", "DISC_TNF", "DISC_FAS", "FADD", "RIP1", "RIP1ub", "RIP1K", "IKK", "NFkB", "BCL2", "ATP", "cIAP", "XIAP", "cFLIP", "Survival"} + }; + + // Read the bool_data.txt file + + string line; + int counter_stable = 0; + + for (int i = 0; i < (*all_cells).size(); i++) + { + // Access the current cell + Cell *pCell = (*all_cells)[i]; + string node; + vector nodes; // Define nodes as a vector + bool is_resistant = false; + int num_stable_states = stable_states.size(); + int j = 0; + while (j < num_stable_states && !is_resistant) { + const auto& state = stable_states[j]; + int counter = 0; + for (const auto& node : state) { + if (pCell->phenotype.intracellular->get_boolean_variable_value(node) == 1) { + counter++; + } + } + if (counter == state.size()) { + is_resistant = true; + counter_stable++; + } + j++; + } + } + // save the number of resistant cells on a file + + file_resistant << counter_stable << endl; + + + return counter_stable; +} + +bool auto_stop() { + + bool condition = false; + bool stop; + // implement here your condition to stop the simulation + + if (condition) { + stop = true; + } else { + stop = false; + } + return stop; +} diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/custom.h b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/custom.h new file mode 100644 index 000000000..ee32d2cc8 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/custom.h @@ -0,0 +1,121 @@ +/* +############################################################################### +# 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 x.y.z) [1]. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# See VERSION.txt or call get_PhysiCell_version() to get the current version # +# x.y.z. Call display_citations() to get detailed information on all cite-# +# able software used in your PhysiCell application. # +# # +# Because PhysiCell extensively uses BioFVM, we suggest you also cite BioFVM # +# as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1], # +# with BioFVM [2] to solve the transport equations. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/jourgs.full_save_interval; # +############################################################################### +# # +# BSD 3-Clause License (see https://opensource.org/licenses/BSD-3-Clause) # +# # +# Copyright (c) 2015-2018, Paul Macklin and the PhysiCell Project # +# All rights reserved. # +# # +# Redistribution and use in source and binary forms, with or without # +# modification, are permitted provided that the following conditions are met: # +# # +# 1. Redistributions of source code must retain the above copyright notice, # +# this list of conditions and the following disclaimer. # +# # +# 2. Redistributions in binary form must reproduce the above copyright # +# notice, this list of conditions and the following disclaimer in the # +# documentation and/or other materials provided with the distribution. # +# # +# 3. Neither the name of the copyright holder nor the names of its # +# contributors may be used to endorse or promote products derived from this # +# software without specific prior written permission. # +# # +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE # +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR # +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF # +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS # +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN # +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) # +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # +# POSSIBILITY OF SUCH DAMAGE. # +# # +############################################################################### +*/ + +#include + +#include "../core/PhysiCell.h" +#include "../modules/PhysiCell_standard_modules.h" + +#include "./tnf_receptor_dynamics.h" +#include "./tnf_boolean_model_interface.h" +#include "../addons/start_and_stop/start_and_stop.h" + +using namespace BioFVM; +using namespace PhysiCell; + +// setup functions to help us along +void create_cell_types( void ); +void setup_tissue( void ); +void setup_tissue_input(void); + +// set up the BioFVM microenvironment +void setup_microenvironment( void ); + +// helper function to read init files +std::vector> read_cells_positions(std::string filename, char delimiter, bool header); + +// helper function to create a sphere of cells of a given radius +std::vector> create_cell_sphere_positions(double cell_radius, double sphere_radius); + +// helper function to create a disc of cells of a given radius +std::vector> create_cell_disc_positions(double cell_radius, double disc_radius); + +// helper function that calculates phere volume +inline float sphere_volume_from_radius(float radius) {return 4/3 * PhysiCell_constants::pi * std::pow(radius, 3);} + +// helper function to inject density surrounding a spheroid +void inject_density_sphere(int density_index, double concentration, double membrane_lenght); + +// helper function to remove a density +void remove_density( int density_index ); + +// custom pathology coloring function +std::vector my_coloring_function( Cell* ); + +double total_live_cell_count(); + +// count the number of total dead cells at current time step +double total_dead_cell_count(); + +// count the number of necrotic cells at current time step +double total_necrosis_cell_count(); + +// function to auto stop the simulation +bool auto_stop_resistance(int alive_cells, int resistant_cells); +bool auto_stop_alive(int alive_cells); +int save_resistant_cells(std::ofstream& file_resistant); +void update_variables_monitor(); +bool auto_stop(); + + + + diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/submodel_data_structures.cpp b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/submodel_data_structures.cpp new file mode 100644 index 000000000..22630e3f8 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/submodel_data_structures.cpp @@ -0,0 +1,114 @@ +/* + * submodel_data_structures.cpp + * + * Created on: ? + * Author: Paul Macklin + * + * Description: + * Auxiliary class to manage submodels + * + * Source: https://github.com/pc4covid19/COVID19/blob/master/PhysiCell/custom_modules/submodel_data_structures.cpp +*/ + +#include "./submodel_data_structures.h" + +using namespace PhysiCell; + +Submodel_Registry submodel_registry; + +Submodel_Information::Submodel_Information( void ) +{ + name = "none"; + version = "-1"; + main_function = NULL; + + microenvironment_variables.resize(0); + cell_variables.resize(0); // custom data + + return; +} + +void Submodel_Information::register_model( void ) +{ + // make sure the cell defaults "know" about each custom variable + // Make sure it's there, or add if it's not. + + for( int n = 0 ; n < cell_variables.size() ; n++ ) + { + // let's do this a bit manually (and inefficiently), but safely. + bool found_it = false; + for( unsigned int i=0; i < cell_defaults.custom_data.variables.size() ; i++ ) + { + if( cell_defaults.custom_data.variables[i].name == cell_variables[n] ) + { found_it = true; } + } + if( found_it == false ) + { + cell_defaults.custom_data.add_variable( cell_variables[n] , "none", 0.0 ); + } + + } + + // add the model to the registry of submodels + + submodel_registry.register_model( *this ); +} + +void Submodel_Information::display( std::ostream& os ) +{ + os << "Submodel: " << name << " (Version " << version << ")" << std::endl ; + + os << "\tcell variables: " << std::endl; + for( int n = 0 ; n < cell_variables.size(); n++ ) + { + os << "\t\t" << cell_variables[n] << std::endl; + } + + os << "\tfunction: " ; + if( main_function ) + { os << (long long int) main_function; } + else + { os << "NULL"; } + os << std::endl; + + return; +} + +void Submodel_Registry::register_model( Submodel_Information& model ) +{ + #pragma omp critical + { + // already registered? + bool found = false; + for( int n = 0; n < submodels.size() ; n++ ) + { + if( submodels[n] == &model ) + { found = true; } + } + + if( found == false ) + { + submodels.push_back( &model ); +// add_software_citation(); +// void add_software_citation( std::string name , std::string version, std::string DOI , std::string URL ) + } + } + + return; +} + +void Submodel_Registry::display( std::ostream& os ) +{ + os << "The following submodels are registered: " << std::endl; + os << "=======================================" << std::endl; + for( int n = 0 ; n < submodels.size(); n++ ) + { + submodels[n]->display( os ); + } + os << std::endl; + + return; +} + + + diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/submodel_data_structures.h b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/submodel_data_structures.h new file mode 100644 index 000000000..f69e34102 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/submodel_data_structures.h @@ -0,0 +1,53 @@ +/* + * submodel_data_structures.h + * + * Created on: ? + * Author: Paul Macklin + * + * Description: + * Auxiliary class to manage submodels + * + * Source: https://github.com/pc4covid19/COVID19/blob/master/PhysiCell/custom_modules/submodel_data_structures.cpp +*/ + + +#include "../core/PhysiCell.h" +#include "../modules/PhysiCell_standard_modules.h" + +using namespace BioFVM; +using namespace PhysiCell; + +#ifndef __submodel_data__ +#define __submodel_data__ + +class Submodel_Information +{ + private: + public: + std::string name; + std::string version; + void(*main_function)(Cell*,Phenotype&,double); + + std::vector< std::string > microenvironment_variables; + std::vector< std::string > cell_variables; // custom data + + Submodel_Information( void ); + + void register_model( void ); + + void display( std::ostream& os ); +}; + +class Submodel_Registry +{ + private: + std::vector submodels; + public: + void register_model( Submodel_Information& model ); + void display( std::ostream& os ); + +}; + +extern Submodel_Registry submodel_registry; + +#endif \ No newline at end of file diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_boolean_model_interface.cpp b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_boolean_model_interface.cpp new file mode 100644 index 000000000..2a4f4abaa --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_boolean_model_interface.cpp @@ -0,0 +1,189 @@ +/* + * tnf_boolean_model_interface.cpp + * + * Created on: 15 jun. 2020 + * Author: Miguel Ponce-de-Leon (miguel.ponce@bsc.es) + * Contributor: Gerard Pradas + * Contributor: Arnau Montagud + * Contributor: Thalia Diniaco + * Cite as: arXiv:2103.14132 [q-bio.QM] + * Description: + * Submodel that work as an interface + * between the Boolean Network (BN) and PhysiCell (PC). + * The submodel run the following steps: + * 1- updates BN input nodes based on custom cell variables (see receptor model) + * 2- updates the BN intracellular model by running MaBoSS + * 3- updates cell state/behaviour based on the state of the BN readout nodes + * + * The update_monitor_variables funtion is just used to keep track of some relevand + * BN nodes' values that are stored as custom variables + */ + + +#include "./tnf_boolean_model_interface.h" + +using namespace PhysiCell; + +Submodel_Information tnf_bm_interface_info; + +void tnf_boolean_model_interface_setup() +{ + tnf_bm_interface_info.name = "TNF Boolean model interface"; + tnf_bm_interface_info.version = "0.2.0"; + + // tnf_bm_interface_info.main_function= update_phenotype_with_signaling; + + // These are just auxiliary variables to keep track of some BN nodes + tnf_bm_interface_info.cell_variables.push_back( "tnf_node" ); + tnf_bm_interface_info.cell_variables.push_back( "fadd_node" ); + tnf_bm_interface_info.cell_variables.push_back( "nfkb_node" ); + + tnf_bm_interface_info.register_model(); +} + +void update_boolean_model_inputs( Cell* pCell, Phenotype& phenotype, double dt ) +{ + static int nR_EB = pCell->custom_data.find_variable_index( "bound_external_TNFR" ); + static int nTNF_threshold = pCell->custom_data.find_variable_index( "TNFR_activation_threshold" ); + + // This if the step transfer function used to update the state of boolean model inputs + // using the state of the receptor dynamics model. The continuos value thresholded is + // the total TNF-recptor complex (doi:10.1016/j.cellsig.2010.08.016) + + if ( pCell->custom_data[nR_EB] >= pCell->custom_data[nTNF_threshold] ){ + + pCell->phenotype.intracellular->set_boolean_variable_value("TNF", 1);} + else + pCell->phenotype.intracellular->set_boolean_variable_value("TNF", 0); + + return; +} + + +/** +void update_cell_from_boolean_model(Cell* pCell, Phenotype& phenotype, double dt) +{ + static int nTNF_external = microenvironment.find_density_index( "tnf" ); + static int nTNF_export_rate = pCell->custom_data.find_variable_index( "TNF_net_production_rate" ); + + static int apoptosis_model_index = phenotype.death.find_death_model_index( "Apoptosis" ); + static int necrosis_model_index = phenotype.death.find_death_model_index( "Necrosis" ); + + // Getting the state of the boolean model readouts (Readout can be in the XML) + bool apoptosis = pCell->phenotype.intracellular->get_boolean_variable_value( "Apoptosis" ); + bool nonACD = pCell->phenotype.intracellular->get_boolean_variable_value( "NonACD" ); + bool survival = pCell->phenotype.intracellular->get_boolean_variable_value( "Survival" ); + bool NFkB = pCell->phenotype.intracellular->get_boolean_variable_value( "NFkB" ); + + if ( apoptosis ) { + pCell->start_death(apoptosis_model_index); + return; + } + + if ( nonACD ) { + 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); + } + + // If NFkB node is active produce some TNF + if ( NFkB ) { + double tnf_export_rate = pCell->custom_data[nTNF_export_rate]; + phenotype.secretion.net_export_rates[nTNF_external] = tnf_export_rate; + } else { + phenotype.secretion.net_export_rates[nTNF_external] = 0; + } + + return; +} +*/ + + + +void update_cell_from_boolean_model(Cell* pCell, Phenotype& phenotype, double dt) +{ + + + static int necrosis_index = phenotype.death.find_death_model_index( PhysiCell_constants::necrosis_death_model ); + static int apoptosis_index = phenotype.death.find_death_model_index( PhysiCell_constants::apoptosis_death_model ); + static int cycle_start_index = live.find_phase_index( PhysiCell_constants::live ); + static int cycle_end_index = live.find_phase_index( PhysiCell_constants::live ); + + static int nTNF_external = microenvironment.find_density_index( "tnf" ); + static int nTNF_export_rate = pCell->custom_data.find_variable_index( "TNF_net_production_rate" ); + static int death_decay_idx = pCell->custom_data.find_variable_index( "death_commitment_decay" ); + + + static float necrosis_rate = pCell->custom_data["necrosis_rate"]; + static float apoptosis_rate = pCell->custom_data["apoptosis_rate"]; + static float death_commitment_decay = pCell->custom_data["death_decay_idx"]; + + // Getting the state of the boolean model readouts (Readout can be in the XML) + bool apoptosis = pCell->phenotype.intracellular->get_boolean_variable_value( "Apoptosis" ); + bool nonACD = pCell->phenotype.intracellular->get_boolean_variable_value( "NonACD" ); + bool survival =pCell->phenotype.intracellular->get_boolean_variable_value( "Survival" ); + bool NFkB = pCell->phenotype.intracellular->get_boolean_variable_value( "NFkB" ); + + + if ( apoptosis ) { + // pCell->start_death(apoptosis_index); + phenotype.death.rates[apoptosis_index] = apoptosis_rate; + } else { + phenotype.death.rates[apoptosis_index] -= apoptosis_rate * death_commitment_decay; + if (phenotype.death.rates[apoptosis_index] < 0) + phenotype.death.rates[apoptosis_index] = 0; + } + + if ( nonACD ) { + //myCounter++; + // pCell->start_death(necrosis_index); + phenotype.death.rates[necrosis_index] = necrosis_rate; + } + else { + phenotype.death.rates[necrosis_index] -= necrosis_rate * death_commitment_decay; + if (phenotype.death.rates[necrosis_index] < 0) + phenotype.death.rates[necrosis_index] = 0; + } + + if ( survival && pCell->phenotype.cycle.current_phase_index() == PhysiCell_constants::Ki67_negative ) + { + pCell->phenotype.cycle.advance_cycle(pCell, phenotype, dt); + } + + // If NFkB node is active produce some TNF + if ( NFkB ) + { + phenotype.secretion.net_export_rates[nTNF_external] = pCell->custom_data[nTNF_export_rate]; + } else + { + phenotype.secretion.net_export_rates[nTNF_external] = 0; + } + // update_cell_and_death_parameters_O2_based(pCell, phenotype, dt); + return; +} + + +void update_behaviors(Cell* pCell, Phenotype& phenotype, double dt) +{ + // update the cell fate based on the boolean outputs + update_cell_from_boolean_model(pCell, phenotype, dt); + + // Get track of some boolean node values for debugging + update_monitor_variables(pCell); +} + +void update_monitor_variables(Cell* pCell ) +{ + static int index_tnf_node = pCell->custom_data.find_variable_index("tnf_node"); + static int index_fadd_node = pCell->custom_data.find_variable_index("fadd_node"); + static int index_nfkb_node = pCell->custom_data.find_variable_index("nfkb_node"); + + pCell->custom_data[index_tnf_node] = pCell->phenotype.intracellular->get_boolean_variable_value("TNF"); + pCell->custom_data[index_fadd_node] = pCell->phenotype.intracellular->get_boolean_variable_value("FADD"); + pCell->custom_data[index_nfkb_node] = pCell->phenotype.intracellular->get_boolean_variable_value( "NFkB" ); + + return; +} diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_boolean_model_interface.h b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_boolean_model_interface.h new file mode 100644 index 000000000..9f1432b32 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_boolean_model_interface.h @@ -0,0 +1,42 @@ +/* + * tnf_boolean_model_interface.cpp + * + * Created on: 15 jun. 2020 + * Author: Miguel Ponce-de-Leon (miguel.ponce@bsc.es) + * Contributor: Gerard Pradas + * Contributor: Arnau Montagud + * Contributor: Thalia Diniaco + * Cite as: arXiv:2103.14132 [q-bio.QM] + * Description: + * Submodel that work as an interface + * between the Boolean Network (BN) and PhysiCell (PC). + * The submodel run the following steps: + * 1- updates BN input nodes based on custom cell variables (see receptor model) + * 2- updates the BN intracellular model by running MaBoSS + * 3- updates cell state/behaviour based on the state of the BN readout nodes + * + * The update_monitor_variables funtion is just used to keep track of some relevand + * BN nodes' values that are stored as custom variables + */ + +#include "../core/PhysiCell.h" +#include "../modules/PhysiCell_standard_modules.h" +#include "../addons/PhysiBoSS/src/maboss_intracellular.h" + +using namespace BioFVM; +using namespace PhysiCell; + +#include "./submodel_data_structures.h" + +void tnf_boolean_model_interface_setup(); + +void update_boolean_model_inputs( Cell* pCell, Phenotype& phenotype, double dt ); + +void update_behaviors(Cell* pCell, Phenotype& phenotype, double dt); + +void update_cell_from_boolean_model(Cell* pCell, Phenotype& phenotype, double dt); + +// helper function to keep updated some cell custom variables +void update_monitor_variables( Cell* pCell ); + +extern int myCounter; diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_receptor_dynamics.cpp b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_receptor_dynamics.cpp new file mode 100644 index 000000000..438095bbb --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_receptor_dynamics.cpp @@ -0,0 +1,143 @@ +/* + * tnf_receptor_dynamics.cpp + * + * Created on: 15 jun. 2020 + * Author: Miguel Ponce-de-Leon (miguel.ponce@bsc.es) + * Contributor: Gerard Pradas + * Contributor: Arnau Montagud + * Contributor: Thalia Diniaco + * Cite as: arXiv:2103.14132 [q-bio.QM] + * Description: + * Simplified model of the TNF receptor (TNFR) dynamics. The models has threed states: + * 1- TNFR unbounded (exposed to the extracellular space) + * 2- TNFR with a TNF bounded + * 3- TNFR with a TNF bounded internalized to be recycled by removeing the TNF ligand. + * + * It is an ODE-based model that has the following equations + * d[unbound_external_TNFR]/dt = TNFR_recycling_rate * [bound_internal_TNFR] - TNFR_binding_rate * [unbound_external_TNFR] * [TNF] + * d[bound_external_TNFR]/dt = NFR_binding_rate * [unbound_external_TNFR] * [TNF] - TNFR_endocytosis_rate * [bound_external_TNFR] + * d[bound_internal_TNFR]/dt = TNFR_endocytosis_rate * [bound_external_TNFR] - TNFR_recycling_rate * [bound_internal_TNFR] + */ + +#include "./tnf_receptor_dynamics.h" + +using namespace PhysiCell; + +Submodel_Information tnf_receptor_info; + +void tnf_receptor_model_setup() +{ + tnf_receptor_info.name = "TNF trasnporter model"; + tnf_receptor_info.version = "0.1.0"; + + 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( "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_endocytosis_rate" ); + tnf_receptor_info.cell_variables.push_back( "TNFR_recycling_rate" ); + tnf_receptor_info.cell_variables.push_back( "TNF_net_production_rate" ); + + tnf_receptor_info.register_model(); + + return; +} + +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_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 double R_binding_rate = pCell->custom_data[nR_bind]; + + static int nR_endo = pCell->custom_data.find_variable_index( "TNFR_endocytosis_rate" ); + static double R_endo_rate = pCell->custom_data[nR_endo]; + + static int nR_recycle = pCell->custom_data.find_variable_index( "TNFR_recycling_rate" ); + static double R_recyc_rate = pCell->custom_data[nR_recycle]; + + + if( phenotype.death.dead == true ) + { return; } + + // 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]; + + // The internalization is only used to track the TNF + // The following part of the code takes care of correcly managed + double dR_EB = phenotype.molecular.internalized_total_substrates[nTNF_external]; + //std::cout << "nTNF_external: " << nTNF_external << std::endl; + //std::cout << "dR_EB: " << dR_EB << std::endl; + + // if it tries to bind more TNF than there are receptors, compensate + if( dR_EB > pCell->custom_data[nR_EU] ) + { + double excess_binding = dR_EB - pCell->custom_data[nR_EU]; + dR_EB = pCell->custom_data[nR_EU]; + // dump any excess back into the microenvironment + static double to_density = 1.0 / microenvironment.mesh.dV; + // this needs omp critical because 2 cells writing to 1 voxel is not thread safe + #pragma omp critical + { pCell->nearest_density_vector()[nTNF_external] += excess_binding * to_density; } + } + + // Remove all the internalized TNF from cell + phenotype.molecular.internalized_total_substrates[nTNF_external] = 0.0; + + // Endocytosis + // The bounded receptor is internalized at a rate R_endo_rate + double dR_IB = dt * R_endo_rate * pCell->custom_data[nR_EB]; + if( dR_IB > pCell->custom_data[nR_EB] ) + { dR_IB = pCell->custom_data[nR_EB]; } + + // 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 + 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]; } + + pCell->custom_data[nR_EU] -= dR_EB; // remove newly bound receptor from R_EU + pCell->custom_data[nR_EB] += dR_EB; // add newly bound receptor to R_EB + + pCell->custom_data[nR_EB] -= dR_IB; // move from external bound + pCell->custom_data[nR_IB] += dR_IB; // move to internal bound + + pCell->custom_data[nR_IB] -= dR_EU; // move from internal unbound + pCell->custom_data[nR_EU] += dR_EU; // move to external unbound + + // update the TNF uptake rate + //std::cout << "uptake_pre: " << phenotype.secretion.uptake_rates[nTNF_external] << std::endl; + phenotype.secretion.uptake_rates[nTNF_external] = R_binding_rate * pCell->custom_data[nR_EU]; + //std::cout << "uptake_post: " << phenotype.secretion.uptake_rates[nTNF_external] << std::endl; + + return; +} + +void tnf_receptor_model_main( double dt ) +{ + #pragma omp parallel for + for( int i=0; i < (*all_cells).size() ; i++ ) + { + Cell* pC = (*all_cells)[i]; + if ( pC->is_out_of_domain ) + { continue; } + if( pC->phenotype.death.dead ) + { continue; } + tnf_receptor_model( pC, pC->phenotype , dt ); + } + return; +} \ No newline at end of file diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_receptor_dynamics.h b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_receptor_dynamics.h new file mode 100644 index 000000000..962ff19bc --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/custom_modules/tnf_receptor_dynamics.h @@ -0,0 +1,39 @@ +/* + * tnf_receptor_dynamics.cpp + * + * Created on: 15 jun. 2020 + * Author: Miguel Ponce-de-Leon + * Contributor: Thalia Diniaco + * Created on: 15 jun. 2020 + * Author: Miguel Ponce-de-Leon (miguel.ponce@bsc.es) + * Contributor: Gerard Pradas + * Contributor: Arnau Montagud + * Contributor: Thalia Diniaco + * Cite as: arXiv:2103.14132 [q-bio.QM] + * Description: + * SSimplified model of the TNF receptor (TNFR) dynamics. The models has threed states: + * 1- TNFR unbounded (exposed to the extracellular space) + * 2- TNFR with a TNF bounded + * 3- TNFR with a TNF bounded internalized to be recycled by removeing the TNF ligand. + * + * It is an ODE-based model that has the following equations + * d[unbound_external_TNFR]/dt = TNFR_recycling_rate * [bound_internal_TNFR] - TNFR_binding_rate * [unbound_external_TNFR] * [TNF] + * d[bound_external_TNFR]/dt = NFR_binding_rate * [unbound_external_TNFR] * [TNF] - TNFR_endocytosis_rate * [bound_external_TNFR] + * d[bound_internal_TNFR]/dt = TNFR_endocytosis_rate * [bound_external_TNFR] - TNFR_recycling_rate * [bound_internal_TNFR] + * + * Cite as: arXiv:2103.14132 [q-bio.QM] + */ + +#include "../core/PhysiCell.h" +#include "../modules/PhysiCell_standard_modules.h" + +using namespace BioFVM; +using namespace PhysiCell; + +#include "./submodel_data_structures.h" + +void tnf_receptor_model_setup(); + +void tnf_receptor_model( Cell* pCell, Phenotype& phenotype, double dt ); + +void tnf_receptor_model_main( double dt ); diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/main-spheroid_TNF.cpp b/sample_projects_intracellular/boolean/spheroid_tnf_model/main-spheroid_TNF.cpp new file mode 100644 index 000000000..d3b452346 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/main-spheroid_TNF.cpp @@ -0,0 +1,410 @@ +/* +############################################################################### +# 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 x.y.z) [1]. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# See VERSION.txt or call get_PhysiCell_version() to get the current version # +# x.y.z. Call display_citations() to get detailed information on all cite-# +# able software used in your PhysiCell application. # +# # +# Because PhysiCell extensively uses BioFVM, we suggest you also cite BioFVM # +# as below: # +# # +# We implemented and solved the model using PhysiCell (Version x.y.z) [1], # +# with BioFVM [2] to solve the transport equations. # +# # +# [1] A Ghaffarizadeh, R Heiland, SH Friedman, SM Mumenthaler, and P Macklin, # +# PhysiCell: an Open Source Physics-Based Cell Simulator for Multicellu- # +# lar Systems, PLoS Comput. Biol. 14(2): e1005991, 2018 # +# DOI: 10.1371/journal.pcbi.1005991 # +# # +# [2] A Ghaffarizadeh, SH Friedman, and P Macklin, BioFVM: an efficient para- # +# llelized diffusive transport solver for 3-D biological simulations, # +# Bioinformatics 32(8): 1256-8, 2016. DOI: 10.1093/bioinformatics/btv730 # +# # +############################################################################### +# # +# BSD 3-Clause License (see https://opensource.org/licenses/BSD-3-Clause) # +# # +# Copyright (c) 2015-2018, Paul Macklin and the PhysiCell Project # +# All rights reserved. # +# # +# Redistribution and use in source and binary forms, with or without # +# modification, are permitted provided that the following conditions are met: # +# # +# 1. Redistributions of source code must retain the above copyright notice, # +# this list of conditions and the following disclaimer. # +# # +# 2. Redistributions in binary form must reproduce the above copyright # +# notice, this list of conditions and the following disclaimer in the # +# documentation and/or other materials provided with the distribution. # +# # +# 3. Neither the name of the copyright holder nor the names of its # +# contributors may be used to endorse or promote products derived from this # +# software without specific prior written permission. # +# # +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE # +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR # +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF # +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS # +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN # +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) # +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # +# POSSIBILITY OF SUCH DAMAGE. # +# # +############################################################################### +*/ + +#include +#include +#include +#include +#include +#include +#include + +#include "./core/PhysiCell.h" +#include "./core/PhysiCell_utilities.h" +#include "./modules/PhysiCell_standard_modules.h" +#include "./addons/PhysiBoSS/src/maboss_intracellular.h" +#include "./custom_modules/custom.h" + +using namespace BioFVM; +using namespace PhysiCell; + +int main( int argc, char* argv[] ) +{ + clock_t T_save_start, T_save_stop, T_reload_start, T_reload_stop, T_total_start, T_total_stop, T_main_start, T_main_stop; + T_total_start = clock(); + T_reload_start = clock(); + // load and parse settings file(s) + std::ofstream file_resistant("output/resistant_cells.txt", std::ios::app); + std::ofstream file_times("output/interesting_times.txt", std::ios::app); + + bool XML_status = false; + char copy_command [1024]; + if( argc > 1 ) + { + XML_status = load_PhysiCell_config_file( argv[1] ); + sprintf( copy_command , "cp %s %s" , argv[1] , PhysiCell_settings.folder.c_str() ); + } + else + { + XML_status = load_PhysiCell_config_file( "./config/PhysiCell_settings.xml" ); + sprintf( copy_command , "cp ./config/PhysiCell_settings.xml %s" , PhysiCell_settings.folder.c_str() ); + } + if( !XML_status ) + { exit(-1); } + + // copy config file to output directry + system( copy_command ); + + // OpenMP setup + omp_set_num_threads(PhysiCell_settings.omp_num_threads); + + // PNRG setup + SeedRandom(); // or specify a seed here + + // time setup + std::string time_units = "min"; + + /* Microenvironment setup */ + + setup_microenvironment(); // modify this in the custom code + + bool start_stop = parameters.bools("start_stop"); + + // User parameters + + double tnf_pulse_period = parameters.doubles("tnf_pulse_period"); + double tnf_pulse_duration = parameters.doubles("tnf_pulse_duration"); + double tnf_pulse_concentration = parameters.doubles("tnf_pulse_concentration"); + double time_remove_tnf = parameters.doubles("time_remove_tnf"); + double membrane_lenght = parameters.doubles("membrane_length"); // radious around which the tnf pulse is injected + + + double tnf_pulse_timer = tnf_pulse_period; + double tnf_pulse_injection_timer = tnf_pulse_duration; // tnf_pulse_duration; // -1; + static int tnf_idx = microenvironment.find_density_index("tnf"); + int resistant_cells = 0; + + + /* PhysiCell setup */ + + // set mechanics voxel size, and match the data structure to BioFVM + double mechanics_voxel_size = 30; + Cell_Container* cell_container = create_cell_container_for_microenvironment( microenvironment, mechanics_voxel_size ); + + /* Users typically start modifying here. START USERMODS */ + create_cell_types(); + + + if( start_stop ){ + + + parameters.bools("read_init") = true; + + // reset cells as they were in the previous simulation + setup_tissue(); + + reset_cell(cell_container->last_cell_cycle_time); + + //exit(-1); + + + reset_global_parameters(cell_container); + + update_variables_monitor(); + + reset_microenv(); + + + } else{ + setup_tissue(); //death model index = 1 == necrotic...= 0 == apoptotic. + } + + + // check if we want to start injecting tnf + if(parameters.bools("if_start_inj")){ + + tnf_pulse_timer = PhysiCell_globals.current_time; + + }else{ + + tnf_pulse_timer = PhysiCell_globals.current_time + tnf_pulse_period; + } + + + // set MultiCellDS save options + + 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 + + char filename[1024]; + sprintf( filename , "%s/initial" , PhysiCell_settings.folder.c_str() ); + + save_PhysiCell_to_MultiCellDS_v2( filename , microenvironment , PhysiCell_globals.current_time ); + + sprintf( filename , "%s/states_initial.csv", PhysiCell_settings.folder.c_str()); + MaBoSSIntracellular::save(filename); + + // save a quick 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 + + std::vector (*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 ); + + sprintf( filename , "%s/legend.svg" , PhysiCell_settings.folder.c_str() ); + create_plot_legend( filename , cell_coloring_function ); + + add_software_citation( "PhysiBoSS" , PhysiBoSS_Version , PhysiBoSS_DOI, PhysiBoSS_URL); + + display_citations(); + + // set the performance timers + BioFVM::RUNTIME_TIC(); + BioFVM::TIC(); + T_reload_stop = clock(); + + 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"<= tnf_pulse_timer ) + { + tnf_pulse_injection_timer = PhysiCell_globals.current_time + tnf_pulse_duration; + tnf_pulse_timer += tnf_pulse_period; + } + + if ( PhysiCell_globals.current_time <= tnf_pulse_injection_timer ) + { + + inject_density_sphere(tnf_idx, tnf_pulse_concentration, membrane_lenght); + + } + + if ( PhysiCell_globals.current_time >= time_remove_tnf ) + { + + remove_density(tnf_idx); + time_remove_tnf += PhysiCell_settings.max_time; + } + + + + // update the microenvironment + microenvironment.simulate_diffusion_decay( diffusion_dt ); + + // update te TNF receptor model of each cell + tnf_receptor_model_main( diffusion_dt ); + + // run PhysiCell + ((Cell_Container *)microenvironment.agent_container)->update_all_cells( PhysiCell_globals.current_time ); + + PhysiCell_globals.current_time += diffusion_dt; + } + + if( PhysiCell_settings.enable_legacy_saves == true ) + { + log_output(PhysiCell_globals.current_time, PhysiCell_globals.full_output_index, microenvironment, report_file); + report_file.close(); + } + } + catch( const std::exception& e ) + { // reference to the base of a polymorphic object + std::cout << e.what(); // information from length_error printed + } + + T_main_stop = clock(); + // save a final simulation snapshot + + sprintf( filename , "%s/final" , PhysiCell_settings.folder.c_str() ); + save_PhysiCell_to_MultiCellDS_v2( filename , microenvironment , PhysiCell_globals.current_time ); + + sprintf( filename , "%s/states_final.csv", PhysiCell_settings.folder.c_str()); + MaBoSSIntracellular::save(filename); + + sprintf( filename , "%s/final.svg" , PhysiCell_settings.folder.c_str() ); + SVG_plot( filename , microenvironment, 0.0 , PhysiCell_globals.current_time, cell_coloring_function ); + + // Save all the files needed for Start & Stop at the right point. + T_save_start = clock(); + save_cell_microenv_data(cell_container); + std::cout << "cells data saved successfully" << std::endl; + T_save_stop = clock(); + + // timer + + std::cout << std::endl << "Total simulation runtime: " << std::endl; + BioFVM::display_stopwatch_value( std::cout , BioFVM::runtime_stopwatch_value() ); + + file_resistant.close(); + T_total_stop = clock(); + double T_save, T_reload, T_total, T_main; + T_save = (double)(T_save_stop - T_save_start)/CLOCKS_PER_SEC; + T_reload = (double)(T_reload_stop - T_reload_start)/CLOCKS_PER_SEC; + T_total = (double)(T_total_stop - T_total_start)/CLOCKS_PER_SEC; + T_main = (double)(T_main_stop - T_main_start)/CLOCKS_PER_SEC; + file_times << T_save << " " << T_reload << " " << T_total <<" " << T_main << std::endl; + file_times.close(); + + return 0; +} diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/params/cell_phases_dict.json b/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/params/cell_phases_dict.json new file mode 100644 index 000000000..cf11a01b2 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/params/cell_phases_dict.json @@ -0,0 +1,23 @@ +{ + "0": "Ki67_positive_premitotic", + "1": "Ki67_positive_postmitotic", + "2": "Ki67_positive", + "3": "Ki67_negative", + "4": "G0G1_phase", + "5": "G0_phase", + "6": "G1_phase", + "7": "G1a_phase", + "8": "G1b_phase", + "9": "G1c_phase", + "10": "S_phase", + "11": "G2M_phase", + "12": "G2_phase", + "13": "M_phase", + "14": "live", + "100": "apoptotic", + "101": "necrotic_swelling", + "102": "necrotic_lysed", + "103": "necrotic", + "104": "debris" + } + \ No newline at end of file diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/params/phases_grouping_dict.json b/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/params/phases_grouping_dict.json new file mode 100644 index 000000000..665b84844 --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/params/phases_grouping_dict.json @@ -0,0 +1,20 @@ +{ + "Ki67_positive_premitotic": "alive", + "Ki67_positive_postmitotic": "alive", + "Ki67_positive": "alive", + "Ki67_negative": "alive", + "G0G1_phase": "alive", + "G0_phase": "alive", + "G1_phase": "alive", + "G1a_phase": "alive", + "G1b_phase": "alive", + "G1c_phase": "alive", + "S_phase": "alive", + "G2M_phase": "alive", + "G2_phase": "alive", + "M_phase": "alive", + "live": "alive", + "apoptotic": "apoptotic", + "necrotic_lysed": "necrotic", + "necrotic_swelling": "necrotic" +} diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/plot_time_course.py b/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/plot_time_course.py new file mode 100644 index 000000000..496a9e3ad --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/plot_time_course.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + + +import sys +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import seaborn as sns + + + + +from pctk.multicellds import MultiCellDS + + +def labelsubplot(ax, idx, xpos=-0.1, ypos=1.05, weight="bold", fontsize=12, color='#434343'): + a_ascci_idx = ord('a') + subplot_label = chr(a_ascci_idx+idx) + ')' + box_aspect = ax.get_box_aspect() + if box_aspect is not None and len(box_aspect) == 3: + ax.text2D(xpos, ypos, subplot_label, weight="bold", ha='left', va='center', + fontsize=fontsize, color=color, transform=ax.transAxes) + else: + ax.text(xpos, ypos, subplot_label, weight="bold", ha='left', va='center', + fontsize=fontsize, color=color, transform=ax.transAxes) + + + + +def get_timeserie_mean(mcds, filter_alive=True): + time = [] + values = [] + filter_alive = True + for t, df in mcds.cells_as_frames_iterator(): + time.append(t) + df = df.iloc[:,3:] + if filter_alive: + mask = df['current_phase'] <= 14 + df = df[mask] + values.append(df.mean(axis=0).values) + + cell_columns = df.columns.tolist() + df = pd.DataFrame(values, columns=cell_columns) + df['time'] = time + return df[['time'] + cell_columns] + + +def get_timeserie_density(mcds): + data = [] + for t,m in mcds.microenvironment_as_matrix_iterator(): + data.append((t, m[5,:].sum())) + df = pd.DataFrame(data=data, columns=['time', 'tnf']) + return df + +def plot_molecular_model(df_cell_variables, list_of_variables, ax1): + + for label in list_of_variables: + y = df_cell_variables[label] + time = df_cell_variables["time"] + ax1.plot(time, y, label="% " + label) + + ax1.set_ylabel("% X") + ax1.yaxis.grid(True) + ax1.set_xlim((0, time.values[-1])) + ax1.set_ylim((0, 1.05)) + + +def plot_cells(df_time_course, color_dict, ax): + + # Alive/Apoptotic/Necrotic vs Time + for k in color_dict: + ax.plot(df_time_course.time, df_time_course[k], "-", c=color_dict[k], label=k) + + # setting axes labels + # ax.set_xlabel("time (min)") + ax.set_ylabel("N of cells") + + # Showing legend + ax.legend() + ax.yaxis.grid(True) + + +def plot_time_course(df_time_course, df_time_tnf, df_cell_variables, list_of_variables): + fig, axes = plt.subplots(2, 1, figsize=(10,4), dpi=300, sharex=True) + + custom_palette = sns.color_palette("deep") + color_dict = {"live": custom_palette[2], "apoptotic": custom_palette[3], "necrotic":custom_palette[5]} + + plot_cells(df_time_course, color_dict, axes[0]) + + for i,ax in enumerate(axes): + labelsubplot(ax, i, xpos=-0.065) + + ax2 = axes[0].twinx() + ax2.plot(df_time_tnf.time, df_time_tnf['tnf'], 'lightgrey') + + ax2.set_yticks([]) + ax2.fill_between(df_time_tnf.time, df_time_tnf['tnf'], color='lightgrey', alpha=0.3) + + axes[0].set_zorder(ax2.get_zorder()+1) + axes[0].patch.set_visible(False) + + plot_molecular_model(df_cell_variables, list_of_variables, axes[1]) + threshold = 0.5 + + axes[1].hlines(threshold, 0, df_time_course.time.iloc[-1]) + axes[1].legend(loc="upper left") + axes[1].set_xlabel("Time (min)") + + + fig.tight_layout() + sns.despine(fig) + + return fig + +def load_datasets(instance_folder, labels_dict): + + + mcds = MultiCellDS(instance_folder) + + df_time_course = mcds.get_cells_summary_frame() + df_time_tnf = get_timeserie_density(mcds) + df_cell_variables = get_timeserie_mean(mcds) + df_cell_variables = df_cell_variables.rename(labels_dict, axis=1) + + if (df_time_course["live"] == 0).sum() > 0: + idx = df_time_course.index[df_time_course["live"] == 0][0] + df_time_course = df_time_course.iloc[:idx] + df_cell_variables = df_cell_variables.iloc[:idx] + df_time_tnf = df_time_tnf.iloc[:idx] + + return df_time_course, df_time_tnf, df_cell_variables + + +sns.set_style("white") +sns.set_palette("deep") + + +def main(): + if len(sys.argv) == 1: + output_folder = "output" + else: + output_folder = sys.argv[1] + + + labels_dict = {} + labels_dict['bound_external_TNFR'] = "TNFR-TNF[e]" + labels_dict['unbound_external_TNFR'] = "TNFR[e]" + labels_dict['bound_internal_TNFR'] = "TNFR-TNF[i]" + df_time_course, df_time_tnf, df_cell_variables = load_datasets(output_folder, labels_dict) + + list_of_variables = list(labels_dict.values()) + fig = plot_time_course(df_time_course, df_time_tnf, df_cell_variables, list_of_variables) + + fig.savefig(f"{output_folder}/Time_course.png") + + +main() \ No newline at end of file diff --git a/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/summarize_simulation.py b/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/summarize_simulation.py new file mode 100755 index 000000000..98a0a160f --- /dev/null +++ b/sample_projects_intracellular/boolean/spheroid_tnf_model/scripts/summarize_simulation.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python +# coding: utf-8 + +import re +import os +import sys +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + + +from pctk.multicellds import MultiCellDS + +def get_timeserie_mean(mcds, filter_alive=True): + time = [] + values = [] + filter_alive = True + for t, df in mcds.cells_as_frames_iterator(): + time.append(t) + df = df.iloc[:,3:] + if filter_alive: + mask = df['current_phase'] <= 14 + df = df[mask] + values.append(df.mean(axis=0).values) + + cell_columns = df.columns.tolist() + df = pd.DataFrame(values, columns=cell_columns) + df['time'] = time + return df[['time'] + cell_columns] + + +def get_timeserie_density(mcds): + data = [] + for t,m in mcds.microenvironment_as_matrix_iterator(): + data.append((t, m[5,:].sum())) + df = pd.DataFrame(data=data, columns=['time', 'tnf']) + return df + +def plot_molecular_model(df_cell_variables, list_of_variables, ax1): + + threshold = 0.5 + + for label in list_of_variables: + y = df_cell_variables[label] + time = df_cell_variables["time"] + ax1.plot(time, y, label="% X " + label) + + ax1.set_ylabel("% X") + ax1.yaxis.grid(True) + ax1.set_xlim((0,time.values[-1])) + ax1.set_ylim((0,1)) + # ax1.set_xlabel("time (min)") + +def plot_cells(df_time_course, color_dict, ax): + + # Alive/Apoptotic/Necrotic vs Time + for k in color_dict: + ax.plot(df_time_course.time, df_time_course[k], "-", c=color_dict[k], label=k) + + # setting axes labels + # ax.set_xlabel("time (min)") + ax.set_ylabel("Nº of cells") + + # Showing legend + ax.legend() + ax.yaxis.grid(True) + +def main(): + + color_dict = {"alive": "g", "apoptotic": "r", "necrotic":"k"} + + + output_folder = sys.argv[1] + + mcds = MultiCellDS(output_folder=output_folder) + + df_time_course = mcds.get_cells_summary_frame() + df_cell_variables = get_timeserie_mean(mcds) + df_time_tnf = get_timeserie_density(mcds) + + # df_time_course.to_csv(instance_folder + "time_course.tsv", sep="\t") + # df_cell_variables.to_csv(instance_folder + "cell_variables.tsv", sep="\t") + # df_time_tnf.to_csv(instance_folder + "tnf_time.tsv", sep="\t") + + fig, axes = plt.subplots(3, 1, figsize=(12,12), dpi=150, sharex=True) + plot_cells(df_time_course, color_dict, axes[0]) + + list_of_variables = ['bound_external_TNFR', 'unbound_external_TNFR', 'bound_internal_TNFR'] + plot_molecular_model(df_cell_variables, list_of_variables, axes[1]) + threshold = 0.5 + + axes[1].hlines(threshold, 0, df_time_course.time.iloc[-1], label="Activation threshold") + ax2 = axes[1].twinx() + ax2.plot(df_time_tnf.time, df_time_tnf['tnf'], 'r', label="[TNF]") + ax2.set_ylabel("[TNF]") + # ax2.set_ylim([0, 1000]) + axes[1].legend(loc="upper left") + ax2.legend(loc="upper right") + + list_of_variables = ['tnf_node', 'nfkb_node', 'fadd_node'] + plot_molecular_model(df_cell_variables, list_of_variables, axes[2]) + axes[2].set_xlabel("time (min)") + ax2 = axes[2].twinx() + ax2.plot(df_time_tnf.time, df_time_tnf['tnf'], 'r', label="[TNF]") + ax2.set_ylabel("[TNF]") + ax2.set_ylim([0, 1000]) + axes[2].legend(loc="upper left") + ax2.legend(loc="upper right") + + fig.tight_layout() + fig.savefig('variables_vs_time.png') + +main() diff --git a/sample_projects_intracellular/boolean/tutorial/Makefile b/sample_projects_intracellular/boolean/tutorial/Makefile index ad562dfa2..425ef4db7 100644 --- a/sample_projects_intracellular/boolean/tutorial/Makefile +++ b/sample_projects_intracellular/boolean/tutorial/Makefile @@ -100,7 +100,9 @@ PhysiCell_custom_module_OBJECTS := custom.o pugixml_OBJECTS := pugixml.o PhysiCell_OBJECTS := $(BioFVM_OBJECTS) $(pugixml_OBJECTS) $(PhysiCell_core_OBJECTS) $(PhysiCell_module_OBJECTS) -ALL_OBJECTS := $(PhysiCell_OBJECTS) $(PhysiCell_custom_module_OBJECTS) $(PhysiBoSS_OBJECTS) + +start_and_stop_OBJECTS := start_and_stop.o +ALL_OBJECTS := $(PhysiCell_OBJECTS) $(PhysiCell_custom_module_OBJECTS) $(PhysiBoSS_OBJECTS) $(PhysiBoSS_module_OBJECTS) $(start_and_stop_OBJECTS) # compile the project @@ -206,6 +208,10 @@ PhysiCell_basic_signaling.o: ./core/PhysiCell_basic_signaling.cpp PhysiCell_geometry.o: ./modules/PhysiCell_geometry.cpp $(COMPILE_COMMAND) -c ./modules/PhysiCell_geometry.cpp +# start_and_stop +start_and_stop.o: ./addons/start_and_stop/start_and_stop.cpp + $(COMPILE_COMMAND) $(INC) -c ./addons/start_and_stop/start_and_stop.cpp + # user-defined PhysiCell modules Compile_MaBoSS: ./addons/PhysiBoSS/MaBoSS/engine/src/BooleanNetwork.h cd ./addons/PhysiBoSS/MaBoSS/engine/src;make CXX=$(CC) MAXNODES=$(MABOSS_MAX_NODES) install_alib;make clean; cd ../../../../.. diff --git a/sample_projects_intracellular/boolean/tutorial/config/simple_tnf/2_Short_TNF.xml b/sample_projects_intracellular/boolean/tutorial/config/simple_tnf/2_Short_TNF.xml index d7b2184b7..2fa17a1ed 100644 --- a/sample_projects_intracellular/boolean/tutorial/config/simple_tnf/2_Short_TNF.xml +++ b/sample_projects_intracellular/boolean/tutorial/config/simple_tnf/2_Short_TNF.xml @@ -255,6 +255,11 @@ + + false + false + false + 0 0 True diff --git a/sample_projects_intracellular/boolean/tutorial/custom_modules/custom.cpp b/sample_projects_intracellular/boolean/tutorial/custom_modules/custom.cpp index f06fa23e1..6d32c77a3 100644 --- a/sample_projects_intracellular/boolean/tutorial/custom_modules/custom.cpp +++ b/sample_projects_intracellular/boolean/tutorial/custom_modules/custom.cpp @@ -213,8 +213,7 @@ void contact_function( Cell* pMe, Phenotype& phenoMe , Cell* pOther, Phenotype& void treatment_function () { if (PhysiCell::parameters.bools.find_index("treatment") != -1) - { - int treatment_substrate_index = BioFVM::microenvironment.find_density_index(PhysiCell::parameters.strings("treatment_substrate")); + { int treatment_substrate_index = BioFVM::microenvironment.find_density_index(PhysiCell::parameters.strings("treatment_substrate")); if (PhysiCell::parameters.bools("treatment")){ diff --git a/sample_projects_intracellular/boolean/tutorial/custom_modules/custom.h b/sample_projects_intracellular/boolean/tutorial/custom_modules/custom.h index f541082b6..5adb5b9f5 100644 --- a/sample_projects_intracellular/boolean/tutorial/custom_modules/custom.h +++ b/sample_projects_intracellular/boolean/tutorial/custom_modules/custom.h @@ -67,6 +67,8 @@ #include "../core/PhysiCell.h" #include "../modules/PhysiCell_standard_modules.h" +#include "../addons/start_and_stop/start_and_stop.h" + using namespace BioFVM; using namespace PhysiCell; @@ -89,4 +91,4 @@ void phenotype_function( Cell* pCell, Phenotype& phenotype, double dt ); void custom_function( Cell* pCell, Phenotype& phenotype , double dt ); void contact_function( Cell* pMe, Phenotype& phenoMe , Cell* pOther, Phenotype& phenoOther , double dt ); -void treatment_function (); \ No newline at end of file +void treatment_function (); diff --git a/sample_projects_intracellular/boolean/tutorial/main.cpp b/sample_projects_intracellular/boolean/tutorial/main.cpp index 0cfa59427..cb88165fb 100644 --- a/sample_projects_intracellular/boolean/tutorial/main.cpp +++ b/sample_projects_intracellular/boolean/tutorial/main.cpp @@ -122,6 +122,8 @@ int main( int argc, char* argv[] ) /* Microenvironment setup */ setup_microenvironment(); // modify this in the custom code + + bool start_stop = parameters.bools("start_stop"); /* PhysiCell setup */ @@ -133,7 +135,27 @@ int main( int argc, char* argv[] ) create_cell_types(); - setup_tissue(); + if( start_stop ){ + + + // reset cells as they were in the previous simulation + setup_tissue(); + + reset_cell(cell_container->last_cell_cycle_time); + + //exit(-1); + + + reset_global_parameters(cell_container); + + + reset_microenv(); + + + } else{ + setup_tissue(); //death model index = 1 == necrotic...= 0 == apoptotic. + } + /* Users typically stop modifying here. END USERMODS */ @@ -184,6 +206,9 @@ int main( int argc, char* argv[] ) report_file<<"simulated time\tnum cells\tnum division\tnum death\twall time"<