From 7d53c634136bab33b3e12c38794f00a4de461bc5 Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Sun, 24 Nov 2024 15:34:39 -0500 Subject: [PATCH 1/5] refactor update_all_cells - fix bug in neighbor list - phagocytosis and fusion could result in removal of a cell - but the neighbor list was not updated - resolve interactions after phenotype updates - this allows the update velocity step to update the neighbors list one time - does move interactions before several steps: - gradient computation - evaluate interactions (so fewer interactions to evaluate) - custom rule - update velocity - springs - do not allow spontaneous spring detachments caused by attack - refactor how timing is done for simplicity --- core/PhysiCell_cell_container.cpp | 89 ++++++++++++++---------------- core/PhysiCell_cell_container.h | 1 - core/PhysiCell_standard_models.cpp | 50 +++++++++-------- 3 files changed, 66 insertions(+), 74 deletions(-) diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index 17851a5ea..496c683a0 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -124,6 +124,14 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me { // secretions and uptakes. Syncing with BioFVM is automated. + static double time_since_last_cycle = phenotype_dt_; + static double time_since_last_mechanics = mechanics_dt_; + static double phenotype_threshold = phenotype_dt_ - 0.5 * diffusion_dt_; + static double mechanics_threshold = mechanics_dt_ - 0.5 * diffusion_dt_; + + bool time_for_phenotype = time_since_last_cycle > phenotype_threshold; + bool time_for_mechanics = time_since_last_mechanics > mechanics_threshold; + #pragma omp parallel for for( int i=0; i < (*all_cells).size(); i++ ) { @@ -133,18 +141,12 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me } } - //if it is the time for running cell cycle, do it! - double time_since_last_cycle= t- last_cell_cycle_time; - - static double phenotype_dt_tolerance = 0.001 * phenotype_dt_; - static double mechanics_dt_tolerance = 0.001 * mechanics_dt_; - // intracellular update. called for every diffusion_dt, but actually depends on the intracellular_dt of each cell (as it can be noisy) #pragma omp parallel for for( int i=0; i < (*all_cells).size(); i++ ) { - if( (*all_cells)[i]->is_out_of_domain == false && initialzed ) { + if( (*all_cells)[i]->is_out_of_domain == false ) { if( (*all_cells)[i]->phenotype.intracellular != NULL && (*all_cells)[i]->phenotype.intracellular->need_update()) { @@ -159,17 +161,12 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me } } - if( time_since_last_cycle > phenotype_dt_ - 0.5 * diffusion_dt_ || !initialzed ) + if( time_for_phenotype ) { // Reset the max_radius in each voxel. It will be filled in set_total_volume // It might be better if we calculate it before mechanics each time // std::fill(max_cell_interactive_distance_in_voxel.begin(), max_cell_interactive_distance_in_voxel.end(), 0.0); - if(!initialzed) - { - time_since_last_cycle = phenotype_dt_; - } - // new as of 1.2.1 -- bundles cell phenotype parameter update, volume update, geometry update, // checking for death, and advancing the cell cycle. Not motility, though. (that's in mechanics) #pragma omp parallel for @@ -180,33 +177,46 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me (*all_cells)[i]->advance_bundled_phenotype_functions( time_since_last_cycle ); } } - // process divides / removes for( int i=0; i < cells_ready_to_divide.size(); i++ ) { cells_ready_to_divide[i]->divide(); } + // remove cells here (and not only below after interactions) in case phenotype and mechanics time steps don't line up for( int i=0; i < cells_ready_to_die.size(); i++ ) { cells_ready_to_die[i]->die(); } - num_divisions_in_current_step+= cells_ready_to_divide.size(); num_deaths_in_current_step+= cells_ready_to_die.size(); + num_divisions_in_current_step+= cells_ready_to_divide.size(); + cells_ready_to_die.clear(); cells_ready_to_divide.clear(); - last_cell_cycle_time= t; - } - double time_since_last_mechanics= t- last_mechanics_time; - - // if( time_since_last_mechanics>= mechanics_dt || !initialzed) - if( time_since_last_mechanics > mechanics_dt_ - 0.5 * diffusion_dt_ || !initialzed ) + time_since_last_cycle = 0.0; + } + else + { time_since_last_cycle += diffusion_dt_; } + + if( time_for_mechanics ) { - if(!initialzed) + // new March 2022: + // run standard interactions (phagocytosis, attack, fusion) here + #pragma omp parallel for + for( int i=0; i < (*all_cells).size(); i++ ) { - time_since_last_mechanics = mechanics_dt_; + Cell* pC = (*all_cells)[i]; + standard_cell_cell_interactions(pC,pC->phenotype,time_since_last_mechanics); } + + for( int i=0; i < cells_ready_to_die.size(); i++ ) + { + cells_ready_to_die[i]->die(); + } + num_deaths_in_current_step+= cells_ready_to_die.size(); + + cells_ready_to_die.clear(); // new February 2018 // if we need gradients, compute them @@ -271,30 +281,6 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me } } - // new March 2022: - // run standard interactions (phagocytosis, attack, fusion) here - #pragma omp parallel for - for( int i=0; i < (*all_cells).size(); i++ ) - { - Cell* pC = (*all_cells)[i]; - standard_cell_cell_interactions(pC,pC->phenotype,time_since_last_mechanics); - } - // super-critical to performance! clear the "dummy" cells from phagocytosis / fusion - // otherwise, comptuational cost increases at polynomial rate VERY fast, as O(10,000) - // dummy cells of size zero are left ot interact mechanically, etc. - if( cells_ready_to_die.size() > 0 ) - { - /* - std::cout << "\tClearing dummy cells from phagocytosis and fusion events ... " << std::endl; - std::cout << "\t\tClearing " << cells_ready_to_die.size() << " cells ... " << std::endl; - // there might be a lot of "dummy" cells ready for removal. Let's do it. - */ - for( int i=0; i < cells_ready_to_die.size(); i++ ) - { cells_ready_to_die[i]->die(); } - cells_ready_to_die.clear(); - } - - // update positions #pragma omp parallel for @@ -309,12 +295,17 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me // Update cell indices in the container for( int i=0; i < (*all_cells).size(); i++ ) + { if(!(*all_cells)[i]->is_out_of_domain && (*all_cells)[i]->is_movable) + { (*all_cells)[i]->update_voxel_in_container(); - last_mechanics_time=t; + } + } + time_since_last_mechanics = 0.0; } + else + { time_since_last_mechanics += diffusion_dt_; } - initialzed=true; return; } diff --git a/core/PhysiCell_cell_container.h b/core/PhysiCell_cell_container.h index c1c7243c1..04b7d8b98 100644 --- a/core/PhysiCell_cell_container.h +++ b/core/PhysiCell_cell_container.h @@ -84,7 +84,6 @@ class Cell_Container : public BioFVM::Agent_Container std::vector cells_ready_to_divide; // the index of agents ready to divide std::vector cells_ready_to_die; int boundary_condition_for_pushed_out_agents; // what to do with pushed out cells - bool initialzed = false; public: BioFVM::Cartesian_Mesh underlying_mesh; diff --git a/core/PhysiCell_standard_models.cpp b/core/PhysiCell_standard_models.cpp index 2a38d3bb9..2e34aa37d 100644 --- a/core/PhysiCell_standard_models.cpp +++ b/core/PhysiCell_standard_models.cpp @@ -1317,36 +1317,36 @@ void standard_cell_cell_interactions( Cell* pCell, Phenotype& phenotype, double // move effector attack here. - if( pCell->phenotype.cell_interactions.pAttackTarget != NULL ) - { - Cell* pTarget = pCell->phenotype.cell_interactions.pAttackTarget; + if( pCell->phenotype.cell_interactions.pAttackTarget != NULL ) + { + Cell* pTarget = pCell->phenotype.cell_interactions.pAttackTarget; - pCell->attack_cell(pTarget,dt); - attacked = true; // attacked at least one cell in this time step + pCell->attack_cell(pTarget,dt); + attacked = true; // attacked at least one cell in this time step - // attack_cell + // attack_cell - // probability of ending attack - // end attack if target is dead - probability = dt / (1e-15 + pCell->phenotype.cell_interactions.attack_duration); + // probability of ending attack + // end attack if target is dead + probability = dt / (1e-15 + pCell->phenotype.cell_interactions.attack_duration); - if( UniformRandom() < probability || pTarget->phenotype.death.dead ) - { - /* - std::cout << "********* ********* ******** attack done **** " << PhysiCell_globals.current_time << " " - << probability << " " - << "attack time: " << pTarget->state.total_attack_time << " " - << "damage: " << pTarget->phenotype.cell_integrity.damage << " " - << "dead? " << (int) pTarget->phenotype.death.dead << " " - << "damage delivered: " << pCell->phenotype.cell_interactions.total_damage_delivered << std::endl; - */ - - detach_cells_as_spring(pCell,pTarget); - - pCell->phenotype.cell_interactions.pAttackTarget = NULL; - } + if( UniformRandom() < probability || pTarget->phenotype.death.dead ) + { + /* + std::cout << "********* ********* ******** attack done **** " << PhysiCell_globals.current_time << " " + << probability << " " + << "attack time: " << pTarget->state.total_attack_time << " " + << "damage: " << pTarget->phenotype.cell_integrity.damage << " " + << "dead? " << (int) pTarget->phenotype.death.dead << " " + << "damage delivered: " << pCell->phenotype.cell_interactions.total_damage_delivered << std::endl; + */ + + detach_cells_as_spring(pCell,pTarget); + + pCell->phenotype.cell_interactions.pAttackTarget = NULL; } + } // move decision if to end attack here. @@ -1468,6 +1468,8 @@ void dynamic_spring_attachments( Cell* pCell , Phenotype& phenotype, double dt ) for( int j=0; j < pCell->state.spring_attachments.size(); j++ ) { Cell* pTest = pCell->state.spring_attachments[j]; + if (phenotype.cell_interactions.pAttackTarget==pTest || pTest->phenotype.cell_interactions.pAttackTarget==pCell) // do not let attackers detach randomly + { continue; } if( UniformRandom() <= detachment_probability ) { detach_cells_as_spring( pCell , pTest ); } } From 2e8fb1f8a4a12ff14e0aa472b0e3ee6016e51155 Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Tue, 26 Nov 2024 09:06:07 -0500 Subject: [PATCH 2/5] fix logic in timing of phenotype/mechanics updates --- core/PhysiCell_cell_container.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index 496c683a0..7900e18b6 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -194,7 +194,7 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me cells_ready_to_die.clear(); cells_ready_to_divide.clear(); - time_since_last_cycle = 0.0; + time_since_last_cycle = diffusion_dt_; // setting it for next cycle } else { time_since_last_cycle += diffusion_dt_; } @@ -301,7 +301,7 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me (*all_cells)[i]->update_voxel_in_container(); } } - time_since_last_mechanics = 0.0; + time_since_last_mechanics = diffusion_dt_; // setting it for next cycle } else { time_since_last_mechanics += diffusion_dt_; } From 7ac3e3967074cbf2d83a96b210a15c2923b376c5 Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Tue, 26 Nov 2024 20:32:09 -0500 Subject: [PATCH 3/5] remove change to detach --- core/PhysiCell_cell_container.cpp | 20 ++++---------------- core/PhysiCell_standard_models.cpp | 2 -- 2 files changed, 4 insertions(+), 18 deletions(-) diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index 7900e18b6..94a1e26bd 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -173,24 +173,17 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me for( int i=0; i < (*all_cells).size(); i++ ) { if( (*all_cells)[i]->is_out_of_domain == false ) - { - (*all_cells)[i]->advance_bundled_phenotype_functions( time_since_last_cycle ); - } + { (*all_cells)[i]->advance_bundled_phenotype_functions( time_since_last_cycle ); } } // process divides / removes for( int i=0; i < cells_ready_to_divide.size(); i++ ) - { - cells_ready_to_divide[i]->divide(); - } + { cells_ready_to_divide[i]->divide(); } // remove cells here (and not only below after interactions) in case phenotype and mechanics time steps don't line up for( int i=0; i < cells_ready_to_die.size(); i++ ) - { - cells_ready_to_die[i]->die(); - } + { cells_ready_to_die[i]->die(); } num_deaths_in_current_step+= cells_ready_to_die.size(); num_divisions_in_current_step+= cells_ready_to_divide.size(); - cells_ready_to_die.clear(); cells_ready_to_divide.clear(); @@ -209,13 +202,8 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me Cell* pC = (*all_cells)[i]; standard_cell_cell_interactions(pC,pC->phenotype,time_since_last_mechanics); } - for( int i=0; i < cells_ready_to_die.size(); i++ ) - { - cells_ready_to_die[i]->die(); - } - num_deaths_in_current_step+= cells_ready_to_die.size(); - + { cells_ready_to_die[i]->die(); } cells_ready_to_die.clear(); // new February 2018 diff --git a/core/PhysiCell_standard_models.cpp b/core/PhysiCell_standard_models.cpp index 2e34aa37d..50884880f 100644 --- a/core/PhysiCell_standard_models.cpp +++ b/core/PhysiCell_standard_models.cpp @@ -1468,8 +1468,6 @@ void dynamic_spring_attachments( Cell* pCell , Phenotype& phenotype, double dt ) for( int j=0; j < pCell->state.spring_attachments.size(); j++ ) { Cell* pTest = pCell->state.spring_attachments[j]; - if (phenotype.cell_interactions.pAttackTarget==pTest || pTest->phenotype.cell_interactions.pAttackTarget==pCell) // do not let attackers detach randomly - { continue; } if( UniformRandom() <= detachment_probability ) { detach_cells_as_spring( pCell , pTest ); } } From 9abcf80ca55447421222df773417889dddd7d329 Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Tue, 26 Nov 2024 20:34:05 -0500 Subject: [PATCH 4/5] edit whitespaces --- core/PhysiCell_cell_container.cpp | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index 94a1e26bd..ffdfabb88 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -173,16 +173,22 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me for( int i=0; i < (*all_cells).size(); i++ ) { if( (*all_cells)[i]->is_out_of_domain == false ) - { (*all_cells)[i]->advance_bundled_phenotype_functions( time_since_last_cycle ); } + { + (*all_cells)[i]->advance_bundled_phenotype_functions(time_since_last_cycle); + } } // process divides / removes for( int i=0; i < cells_ready_to_divide.size(); i++ ) - { cells_ready_to_divide[i]->divide(); } + { + cells_ready_to_divide[i]->divide(); + } // remove cells here (and not only below after interactions) in case phenotype and mechanics time steps don't line up for( int i=0; i < cells_ready_to_die.size(); i++ ) - { cells_ready_to_die[i]->die(); } - num_deaths_in_current_step+= cells_ready_to_die.size(); + { + cells_ready_to_die[i]->die(); + } num_divisions_in_current_step+= cells_ready_to_divide.size(); + num_deaths_in_current_step+= cells_ready_to_die.size(); cells_ready_to_die.clear(); cells_ready_to_divide.clear(); @@ -203,7 +209,9 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me standard_cell_cell_interactions(pC,pC->phenotype,time_since_last_mechanics); } for( int i=0; i < cells_ready_to_die.size(); i++ ) - { cells_ready_to_die[i]->die(); } + { + cells_ready_to_die[i]->die(); + } cells_ready_to_die.clear(); // new February 2018 @@ -279,8 +287,6 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me { pC->update_position(time_since_last_mechanics); } } - // When somebody reviews this code, let's add proper braces for clarity!!! - // Update cell indices in the container for( int i=0; i < (*all_cells).size(); i++ ) { From 686f16be9b1f418479bd2b16f388cfc491c03a73 Mon Sep 17 00:00:00 2001 From: Daniel Bergman Date: Mon, 2 Dec 2024 11:06:44 -0500 Subject: [PATCH 5/5] rename var to time_since_last_phenotype --- core/PhysiCell_cell_container.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/core/PhysiCell_cell_container.cpp b/core/PhysiCell_cell_container.cpp index ffdfabb88..cec097d9d 100644 --- a/core/PhysiCell_cell_container.cpp +++ b/core/PhysiCell_cell_container.cpp @@ -124,12 +124,12 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me { // secretions and uptakes. Syncing with BioFVM is automated. - static double time_since_last_cycle = phenotype_dt_; + static double time_since_last_phenotype = phenotype_dt_; static double time_since_last_mechanics = mechanics_dt_; static double phenotype_threshold = phenotype_dt_ - 0.5 * diffusion_dt_; static double mechanics_threshold = mechanics_dt_ - 0.5 * diffusion_dt_; - bool time_for_phenotype = time_since_last_cycle > phenotype_threshold; + bool time_for_phenotype = time_since_last_phenotype > phenotype_threshold; bool time_for_mechanics = time_since_last_mechanics > mechanics_threshold; #pragma omp parallel for @@ -174,7 +174,7 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me { if( (*all_cells)[i]->is_out_of_domain == false ) { - (*all_cells)[i]->advance_bundled_phenotype_functions(time_since_last_cycle); + (*all_cells)[i]->advance_bundled_phenotype_functions(time_since_last_phenotype); } } // process divides / removes @@ -193,10 +193,10 @@ void Cell_Container::update_all_cells(double t, double phenotype_dt_ , double me cells_ready_to_die.clear(); cells_ready_to_divide.clear(); - time_since_last_cycle = diffusion_dt_; // setting it for next cycle + time_since_last_phenotype = diffusion_dt_; // setting it for next cycle } else - { time_since_last_cycle += diffusion_dt_; } + { time_since_last_phenotype += diffusion_dt_; } if( time_for_mechanics ) {