diff --git a/.gitignore b/.gitignore index 7d30e11..2881546 100644 --- a/.gitignore +++ b/.gitignore @@ -121,4 +121,7 @@ doc/latex *.dat # Draft code for comparing models -draft_code_my_own_analysis/ \ No newline at end of file +draft_code_my_own_analysis/ + +# Temporary files +cout.txt \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 4d2fdcd..d8b452d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.19.3) -project(cart_tumor) +project(CARTopiaX) # BioDynaMo curretly uses the C++17 standard. diff --git a/bdm.toml b/bdm.toml index 5557514..d681e58 100644 --- a/bdm.toml +++ b/bdm.toml @@ -10,5 +10,12 @@ additional_data_members = [ {name = "cell_type", function = "GetTypeAsInt"} ] +[[visualize_agent]] +name = "CarTCell" +additional_data_members = ["diameter_","volume_"] + [[visualize_diffusion]] name = "oxygen" + +[[visualize_diffusion]] +name = "immunostimulatory_factor" diff --git a/src/cart_cell.cc b/src/cart_cell.cc index bab5c30..745d994 100644 --- a/src/cart_cell.cc +++ b/src/cart_cell.cc @@ -24,6 +24,7 @@ #include "tumor_cell.h" #include "utils_aux.h" #include "core/agent/agent.h" +#include "core/agent/agent_pointer.h" #include "core/agent/new_agent_event.h" #include "core/container/math_array.h" #include "core/diffusion/diffusion_grid.h" @@ -31,7 +32,9 @@ #include "core/interaction_force.h" #include "core/real_t.h" #include "core/resource_manager.h" +#include "core/simulation.h" #include "core/util/log.h" +#include "core/util/random.h" #include #include #include @@ -47,6 +50,12 @@ CarTCell::CarTCell(const Real3& position) { oxygen_dgrid_ = rm.GetDiffusionGrid("oxygen"); immunostimulatory_factor_dgrid_ = rm.GetDiffusionGrid("immunostimulatory_factor"); + + SetCurrentLiveTime(kAverageMaximumTimeUntillApoptosisCart); + // Add Consumption and Secretion + // Set default oxygen consumption rate + SetOxygenConsumptionRate(kDefaultOxygenConsumptionCarT); + // Compute constants for all ConsumptionSecretion of Oxygen ComputeConstantsConsumptionSecretion(); } @@ -138,6 +147,8 @@ void CarTCell::ChangeVolumeExponentialRelaxationEquation( // compute Displacement Real3 CarTCell::CalculateDisplacement(const InteractionForce* force, real_t squared_radius, real_t /*dt*/) { + Simulation* sim = Simulation::GetActive(); + // real_t h = dt; Real3 movement_at_next_step{0, 0, 0}; // this should be chaged in a future version of BioDynaMo in order to have a @@ -148,33 +159,114 @@ Real3 CarTCell::CalculateDisplacement(const InteractionForce* force, Real3 translation_velocity_on_point_mass{0, 0, 0}; - // We check for every neighbor object if they touch us, i.e. push us - // away and agreagate the velocities + //-------------------------------------------- + // CAR-T self motility (in case of migration) + //-------------------------------------------- + Real3 current_position = GetPosition(); + ExecutionContext* ctxt = sim->GetExecutionContext(); + Random* rng = sim->GetRandom(); + Real3 motility; + if (DoesCellMove()) { + // compute motility + if (rng->Uniform(0.0, 1.0) < kMotilityProbability) { + // random direction as unitary vector + const Real3 random_direction = GenerateRandomDirection(); + Real3 direction_to_immunostimulatory_factor; + // returns normalized gradient towards the immunostimulatory factor source + immunostimulatory_factor_dgrid_->GetGradient( + current_position, &direction_to_immunostimulatory_factor, true); + // motility = bias * direction_to_immunostimulatory_factor + + // (1-bias)*random_direction + motility = kMigrationBiasCart * direction_to_immunostimulatory_factor + + kMigrationOneMinusBiasCart * random_direction; + const real_t motility_norm_squared = motility[0] * motility[0] + + motility[1] * motility[1] + + motility[2] * motility[2]; + // Convert to unit direction + if (motility_norm_squared > 0) { + motility.Normalize(); + } + // Scale by migration speed and add to the velocity + translation_velocity_on_point_mass += motility * kMigrationSpeedCart; + } + } + + //-------------------------------------------- + // CAR-T killing or victim cell escaping + //-------------------------------------------- + // If cell is not apoptotic + if (state_ == CarTCellState::kAlive) { + // if it is attached to tumor cell + if (attached_to_tumor_cell_) { + // try to kill the cancer cell and in case of failure see if it manages to + // scape the order needs to be this one because it should try to kill + // before seeing if it scapes + if (TryToInduceApoptosis(attached_cell_ptr_, rng) || + rng->Uniform(0.0, 1.0) < kProbabilityEscape) { + // the cancer cell is detached + attached_cell_ptr_->SetAttachedToCart(false); + // empty ID + attached_cell_ptr_ = nullptr; + attached_to_tumor_cell_ = false; + } + } - uint64_t non_zero_neighbor_forces = 0; - if (!IsStatic()) { - auto* ctxt = Simulation::GetActive()->GetExecutionContext(); - auto calculate_neighbor_forces = + //-------------------------------------------- + // CAR-T adhesion to victim cell + //-------------------------------------------- + // Compute forces between the cells and check for a new attachment + auto calculate_forces_and_elastic_displacement = L2F([&](Agent* neighbor, real_t /*squared_distance*/) { - auto neighbor_force = force->Calculate(this, neighbor); - if (neighbor_force[0] != 0 || neighbor_force[1] != 0 || - neighbor_force[2] != 0) { - non_zero_neighbor_forces++; - translation_velocity_on_point_mass[0] += neighbor_force[0]; - translation_velocity_on_point_mass[1] += neighbor_force[1]; - translation_velocity_on_point_mass[2] += neighbor_force[2]; + // Adhesion repulsion forces between cells + // We check for every neighbor object if they touch us, i.e. push us + // away and aggregate the velocities + Real4 neighbor_force = force->Calculate(this, neighbor); + translation_velocity_on_point_mass[0] += neighbor_force[0]; + translation_velocity_on_point_mass[1] += neighbor_force[1]; + translation_velocity_on_point_mass[2] += neighbor_force[2]; + + // CAR-T adhesion to new victim cell + Real3 displac = neighbor->GetPosition() - current_position; + + if (auto* cancer_cell = dynamic_cast(neighbor)) { + // movement towards the tumor cells + const real_t sq_norm_displac = displac[0] * displac[0] + + displac[1] * displac[1] + + displac[2] * displac[2]; + + // The cart moves towards the tumor cell only if they are not + // touching already If they are too close the only force affecting + // is the adhesion force to avoid CAR-T non-stop pushing tumor + // cells. In case of being closer than + // gKMaxSquaredDistanceCartMovingTowardsTumorCell there is a + // probability kProbabilityPushing for the CAR-T to keep pushing the + // tumor cell + if (sq_norm_displac > + gKMaxSquaredDistanceCartMovingTowardsTumorCell) { + translation_velocity_on_point_mass[0] += + displac[0] * kElasticConstantCart; + translation_velocity_on_point_mass[1] += + displac[1] * kElasticConstantCart; + translation_velocity_on_point_mass[2] += + displac[2] * kElasticConstantCart; + } + + // If the CAR-T has not succeeded in attaching to a tumor cell yet, + // it tries again + if (!attached_to_tumor_cell_) { + TryToGetAttachedTo(cancer_cell, sq_norm_displac, rng); + } } }); - ctxt->ForEachNeighbor(calculate_neighbor_forces, *this, squared_radius); - - if (non_zero_neighbor_forces > 1) { - SetStaticnessNextTimestep(false); - } + ctxt->ForEachNeighbor(calculate_forces_and_elastic_displacement, *this, + squared_radius); } + //-------------------------------------------- // Two step Adams-Bashforth approximation of the time derivative for position // position(t + dt) ≈ position(t) + dt * [ 1.5 * velocity(t) - 0.5 * // velocity(t - dt) ] + //-------------------------------------------- movement_at_next_step += translation_velocity_on_point_mass * kDnew + older_velocity_ * kDold; @@ -184,6 +276,97 @@ Real3 CarTCell::CalculateDisplacement(const InteractionForce* force, return movement_at_next_step; } +// Try to get attached to a tumor cell +void CarTCell::TryToGetAttachedTo(TumorCell* victim, real_t squared_distance, + Random* rng) { + // If the tumor cell is not already attached to a CAR-T cell, is not dead and + // is not too far away. + if (!victim->IsAttachedToCart() && !victim->IsDead() && + squared_distance < kSquaredMaxAdhesionDistanceCart) { + // factor of how high is the oncoprotein level of the cancer cell + real_t oncoprotein_scale_factor = + (victim->GetOncoproteinLevel() - kOncoproteinLimit) / + kOncoproteinDifference; + // Clamp scale_factor to be in [0,1] + if (oncoprotein_scale_factor > 1.0) { + oncoprotein_scale_factor = 1.0; + } + // If oncoprotein level is lower than the limit the cancer cell does not get + // detected + if (oncoprotein_scale_factor <= 0.0) { + // oncoprotein_scale_factor = 0.0; the probability is going to be 0 so + // return the function is the most efficient + return; + } + + // factor of how far the cancer cell is + real_t distance_scale_factor = + (kMaxAdhesionDistanceCart - std::sqrt(squared_distance)) / + kDifferenceCartAdhesionDistances; + // Clamp scale_factor to be in [0,1]. We already checked that it is > 0 + // because squared_distance < kSquaredMaxAdhesionDistanceCart + if (distance_scale_factor > 1.0) { + distance_scale_factor = 1.0; + } + + // It tries to attach the CAR-T cell to the tumor cell with probability + // kAdhesionRateCart * oncoprotein_scale_factor * distance_scale_factor * + // kDtMechanics + if (rng->Uniform(0.0, 1.0) < kAdhesionRateCart * oncoprotein_scale_factor * + distance_scale_factor * kDtMechanics) { +// avoid race condition. Only one cell can be attached to the tumor cell. +#pragma omp critical + { + // We need to check again if the victim is not attached to a CAR-T cell + // yet. This could be made more efficiently with a semaphore for each + // cancer cell + if (!victim->IsAttachedToCart()) { + attached_to_tumor_cell_ = true; + attached_cell_ptr_ = victim->GetAgentPtr(); + victim->SetAttachedToCart(true); + } + } + } + } +} + +// Try to induce apoptosis +bool CarTCell::TryToInduceApoptosis(bdm::AgentPointer attached_cell, + Random* rng) const { + // If there is no attached cell (this should not happen) + + if (!attached_to_tumor_cell_) { + return false; + } + + // factor of how high is the oncoprotein levelof the cancer cell + real_t scale_factor = + (attached_cell->GetOncoproteinLevel() - kOncoproteinLimit) / + kOncoproteinDifference; + // Clamp scale_factor to be in [0,1] + if (scale_factor > 1.0) { + scale_factor = 1.0; + } + // If oncoprotein level is lower than the limit the cancer cell does not + // become apoptotic + if (scale_factor < 0.0) { + // scale_factor = 0.0; the probability is going to be 0 so return the + // function is the most efficient + return false; + } + // CAR-T success of killing probability: aggressive cancer cells (high + // oncoprotein level) are more likely to be killed + const bool succeeded = + rng->Uniform(0.0, 1.0) < kKillRateCart * scale_factor * kDtMechanics; + + // The CAR-T has succeeded to induce apoptosis on the Cancer Cell + if (succeeded) { + attached_cell->StartApoptosis(); + } + + return succeeded; +} + // Compute new oxygen or immunostimulatory factor concentration after // consumption/ secretion real_t CarTCell::ConsumeSecreteSubstance(int substance_id, @@ -194,7 +377,7 @@ real_t CarTCell::ConsumeSecreteSubstance(int substance_id, res = (old_concentration + constant1_oxygen_) / constant2_oxygen_; } else if (substance_id == immunostimulatory_factor_dgrid_->GetContinuumId()) { - // This point should never be reached + // CAR-T do not change immunostimulatory factor levels res = old_concentration; } else { throw std::invalid_argument("Unknown substance id: " + @@ -230,7 +413,7 @@ void CarTCell::ComputeConstantsConsumptionSecretion() { /// Main behavior executed at each simulation step void StateControlCart::Run(Agent* agent) { - auto* sim = Simulation::GetActive(); + Simulation* sim = Simulation::GetActive(); // Run only every kDtCycle minutes, fmod does not work with the type // returned by GetSimulatedTime() if (sim->GetScheduler()->GetSimulatedSteps() % kStepsPerCycle != 0) { @@ -241,7 +424,7 @@ void StateControlCart::Run(Agent* agent) { switch (cell->GetState()) { case CarTCellState::kAlive: { // the cell is growing to real_t its size before mitosis - // Probability of death= 1/CurrentLiveTime, division by 0 + // Probability of death= 1/CurrentLiveTime,avoiding division by 0 if (sim->GetRandom()->Uniform(1.0) < kDtCycle / std::max(cell->GetCurrentLiveTime(), kEpsilon)) { // the cell Dies @@ -260,8 +443,8 @@ void StateControlCart::Run(Agent* agent) { cell->ComputeConstantsConsumptionSecretion(); // Detach from tumor cell if it was attached if (cell->IsAttachedToTumorCell()) { - cell->GetAttachedCell()->SetAttachedToCart(false); - cell->SetAttachedCell(nullptr); + cell->GetAttachedCellPointer()->SetAttachedToCart(false); + cell->SetAttachedCellPointer(nullptr); cell->SetAttachedToTumorCell(false); } } else { @@ -284,7 +467,7 @@ void StateControlCart::Run(Agent* agent) { // duration transition) if (kTimeApoptosis < cell->GetTimerState()) { // remove the cell from the simulation - auto* ctxt = sim->GetExecutionContext(); + ExecutionContext* ctxt = sim->GetExecutionContext(); ctxt->RemoveAgent(agent->GetUid()); } break; diff --git a/src/cart_cell.h b/src/cart_cell.h index 7a92ed6..54a6c16 100644 --- a/src/cart_cell.h +++ b/src/cart_cell.h @@ -24,6 +24,7 @@ #include "tumor_cell.h" #include "core/agent/agent.h" +#include "core/agent/agent_pointer.h" #include "core/agent/cell.h" #include "core/agent/new_agent_event.h" #include "core/behavior/behavior.h" @@ -31,6 +32,7 @@ #include "core/diffusion/diffusion_grid.h" #include "core/interaction_force.h" #include "core/real_t.h" +#include "core/util/random.h" namespace bdm { @@ -124,8 +126,12 @@ class CarTCell : public Cell { real_t GetCurrentLiveTime() const { return current_live_time_; } void SetCurrentLiveTime(real_t time) { current_live_time_ = time; } - TumorCell* GetAttachedCell() const { return attached_cell_; } - void SetAttachedCell(TumorCell* cell) { attached_cell_ = cell; } + bdm::AgentPointer GetAttachedCellPointer() const { + return attached_cell_ptr_; + } + void SetAttachedCellPointer(bdm::AgentPointer cell_ptr) { + attached_cell_ptr_ = cell_ptr; + } /// Returns whether the cell moves by its own bool DoesCellMove(); @@ -184,6 +190,27 @@ class CarTCell : public Cell { /// equations. void ComputeConstantsConsumptionSecretion(); + /// Try to get attached to a tumor cell + /// + /// @param victim The tumor cell to which the CAR-T cell tries to attach + /// @param squared_distance The squared distance between the CAR-T cell and + /// the tumor cell + /// @param rng Pointer to the random number generator + /// Attempts to attach the CAR-T cell to a selected tumor cell. + void TryToGetAttachedTo(TumorCell* victim, real_t squared_distance, + Random* rng); + + /// Try to induce apoptosis + /// + /// Tries stochastically to induce apoptosis in the attached tumor cell and in + /// case of success induces the apoptosis + /// + /// @param attached_cell The tumor cell to which the CAR-T cell is attached + /// @param rng Pointer to the random number generator + /// @return true if apoptosis was induced, false otherwise + bool TryToInduceApoptosis(bdm::AgentPointer attached_cell, + Random* rng) const; + private: /// Current state of the CAR-T cell CarTCellState state_ = CarTCellState::kAlive; @@ -224,7 +251,7 @@ class CarTCell : public Cell { real_t target_relation_cytoplasm_nucleus_ = 0.0; /// Velocity of the cell in the previous time step - Real3 older_velocity_ = {0.0, 0.0, 0.0}; + Real3 older_velocity_ = {0, 0, 0}; /// Rate of oxygen consumption by the cell real_t oxygen_consumption_rate_ = 0.0; @@ -239,7 +266,7 @@ class CarTCell : public Cell { real_t constant2_oxygen_ = 0.0; /// Pointer to the attached tumor cell - TumorCell* attached_cell_ = nullptr; + bdm::AgentPointer attached_cell_ptr_; }; /// Behavior class for controlling CAR-T cell state transitions diff --git a/src/cart_tumor.cc b/src/cart_tumor.cc index 1097509..6dade31 100644 --- a/src/cart_tumor.cc +++ b/src/cart_tumor.cc @@ -33,6 +33,8 @@ #include "core/operation/operation.h" #include "core/param/param.h" #include "core/real_t.h" +#include "core/resource_manager.h" +#include "core/scheduler.h" #include "core/simulation.h" #include #include @@ -48,18 +50,20 @@ int Simulate(int argc, const char** argv) { // Periodic boundary conditions param->bound_space = Param::BoundSpaceMode::kTorus; // Cube of kBoundedSpaceLength³ centered at origin - param->min_bound = -kBoundedSpaceLength / 2.0; - param->max_bound = kBoundedSpaceLength / 2.0; + param->min_bound = -kBoundedSpaceLength / kHalf; + param->max_bound = kBoundedSpaceLength / kHalf; param->simulation_time_step = kDt; + // for outputing performance statistics + param->statistics = kOutputPerformanceStatistics; }; Simulation simulation(argc, argv, set_param); - auto* ctxt = simulation.GetExecutionContext(); + ExecutionContext* ctxt = simulation.GetExecutionContext(); // Change Forces - auto* scheduler = simulation.GetScheduler(); + Scheduler* scheduler = simulation.GetScheduler(); - auto* op = scheduler->GetOps("mechanical forces")[0]; + Operation* op = scheduler->GetOps("mechanical forces")[0]; std::unique_ptr interaction_velocity = std::make_unique(); op->GetImplementation()->SetInteractionForce( @@ -71,7 +75,7 @@ int Simulate(int argc, const char** argv) { env->SetBoxLength(gKLengthBoxMechanics); // Define Substances - auto* rm = Simulation::GetActive()->GetResourceManager(); + ResourceManager* rm = Simulation::GetActive()->GetResourceManager(); // Oxygen // substance_id, name, diffusion_coefficient, decay_constant, resolution, @@ -80,7 +84,7 @@ int Simulate(int argc, const char** argv) { std::make_unique( kOxygen, "oxygen", kDiffusionCoefficientOxygen, kDecayConstantOxygen, kResolutionGridSubstances, kDtSubstances, - true); // true indicates Dirichlet border conditions + /*dirichlet_border=*/true); rm->AddContinuum(oxygen_grid.release()); // Immunostimulatory Factor @@ -91,8 +95,7 @@ int Simulate(int argc, const char** argv) { kDiffusionCoefficientImmunostimulatoryFactor, kDecayConstantImmunostimulatoryFactor, kResolutionGridSubstances, kDtSubstances, - // false indicates Neumann border conditions - false); + /*dirichlet_border=*/false); rm->AddContinuum(immunostimulatory_factor_grid.release()); // Boundary Conditions Dirichlet: simulating absorption or total loss at the @@ -128,6 +131,14 @@ int Simulate(int argc, const char** argv) { ctxt->AddAgent(tumor_cell.release()); } + // Treatment administration operation + std::unique_ptr treatment_op = + std::make_unique("SpawnCart", kStepsOneDay); + std::unique_ptr spawn_cart = + std::make_unique(); + treatment_op->AddOperationImpl(bdm::kCpu, spawn_cart.release()); + scheduler->ScheduleOp(treatment_op.release()); + // OutputSummary operation std::unique_ptr summary_op = std::make_unique("OutputSummary", kOutputCsvInterval); diff --git a/src/diffusion_thomas_algorithm.cc b/src/diffusion_thomas_algorithm.cc index efd03b7..3afc242 100644 --- a/src/diffusion_thomas_algorithm.cc +++ b/src/diffusion_thomas_algorithm.cc @@ -24,10 +24,12 @@ #include "hyperparams.h" #include "tumor_cell.h" #include "core/agent/agent.h" +#include "core/container/math_array.h" #include "core/diffusion/diffusion_grid.h" #include "core/real_t.h" #include "core/resource_manager.h" #include +#include #include #include #include @@ -85,7 +87,7 @@ void DiffusionThomasAlgorithm::InitializeThomasAlgorithmVectors( // Apply Dirichlet boundary conditions to the grid void DiffusionThomasAlgorithm::ApplyDirichletBoundaryConditions() { // FIXME: Fix BioDynaMo by returning a view or c++20 std::span. - const auto* dimensions_ptr = GetDimensionsPtr(); + const int32_t* dimensions_ptr = GetDimensionsPtr(); // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) const real_t origin = dimensions_ptr[0]; const real_t simulated_time = GetSimulatedTime(); @@ -158,7 +160,7 @@ void DiffusionThomasAlgorithm::ApplyDirichletBoundaryConditions() { // Sets the concentration at a specific voxel void DiffusionThomasAlgorithm::SetConcentration(size_t idx, real_t amount) { // NOLINTBEGIN(cppcoreguidelines-pro-bounds-pointer-arithmetic) - const auto* all_concentrations = GetAllConcentrations(); + const real_t* all_concentrations = GetAllConcentrations(); const real_t current_concentration = all_concentrations[idx]; ChangeConcentrationBy(idx, amount - current_concentration, InteractionMode::kAdditive, false); @@ -264,7 +266,7 @@ void DiffusionThomasAlgorithm::ForwardElimination( // Get initial index based on direction size_t ind = GetLoopIndex(direction, outer, middle, 0); // NOLINTBEGIN(cppcoreguidelines-pro-bounds-pointer-arithmetic) - const auto* all_concentrations = GetAllConcentrations(); + const real_t* all_concentrations = GetAllConcentrations(); const real_t initial_concentration = all_concentrations[ind]; SetConcentration(ind, initial_concentration / thomas_denom[0]); @@ -284,7 +286,7 @@ void DiffusionThomasAlgorithm::BackSubstitution( int direction, int outer, int middle, const std::vector& thomas_c, int jump) { // NOLINTBEGIN(cppcoreguidelines-pro-bounds-pointer-arithmetic) - const auto* all_concentrations = GetAllConcentrations(); + const real_t* all_concentrations = GetAllConcentrations(); // Back substitution loop for (int inner = resolution_ - 2; inner >= 0; inner--) { @@ -315,14 +317,14 @@ void DiffusionThomasAlgorithm::ComputeConsumptionsSecretions() { // This method is called to compute the consumptions and secretions of // substances by the tumor cells. It iterates over all agents and applies the // consumption and secretion behaviors defined in the TumorCell class. - auto* rm = bdm::Simulation::GetActive()->GetResourceManager(); + ResourceManager* rm = bdm::Simulation::GetActive()->GetResourceManager(); // in a future version of BioDynaMo this should be parallelized getting the // agents inside each chemical voxel and treating each voxel independently. // Fixme:this dynamic casting could be optimized in a future version rm->ForEachAgent([this](bdm::Agent* agent) { if (auto* cell = dynamic_cast(agent)) { // Handle TumorCell agents - const auto& pos = cell->GetPosition(); + const Real3& pos = cell->GetPosition(); const real_t conc = this->GetValue(pos); const real_t new_conc = cell->ConsumeSecreteSubstance(GetContinuumId(), conc); @@ -330,7 +332,7 @@ void DiffusionThomasAlgorithm::ComputeConsumptionsSecretions() { InteractionMode::kAdditive, false); } else if (auto* cell = dynamic_cast(agent)) { // Handle CarTCell agents - const auto& pos = cell->GetPosition(); + const Real3& pos = cell->GetPosition(); const real_t conc = GetValue(pos); const real_t new_conc = cell->ConsumeSecreteSubstance(GetContinuumId(), conc); diff --git a/src/forces_tumor_cart.cc b/src/forces_tumor_cart.cc index 157a56f..850093c 100644 --- a/src/forces_tumor_cart.cc +++ b/src/forces_tumor_cart.cc @@ -124,9 +124,8 @@ Real4 InteractionVelocity::Calculate(const Agent* lhs, const Agent* rhs) const { } const real_t force_magnitude = temp_r / distance; - return {2 * force_magnitude * displacement[0], - 2 * force_magnitude * displacement[1], - 2 * force_magnitude * displacement[2], + return {force_magnitude * displacement[0], force_magnitude * displacement[1], + force_magnitude * displacement[2], // 4th component is unused 0.0}; } diff --git a/src/hyperparams.h b/src/hyperparams.h index 57ae8b1..8d2aa3f 100644 --- a/src/hyperparams.h +++ b/src/hyperparams.h @@ -26,13 +26,13 @@ #include "core/util/math.h" #include #include +#include namespace bdm { /// This file contains hyperparameters used in the simulation. Change: In a /// future version it needs to be changed into a params file with no need to be /// recompiled - /// /// TumorCell Hyperparameters /// @@ -69,12 +69,20 @@ constexpr real_t kDefaultVolumeNewTumorCell = 2494.0; constexpr real_t kDefaultVolumeNucleusTumorCell = 540.0; /// Default fraction of fluid volume in a new tumor cell constexpr real_t kDefaultFractionFluidTumorCell = 0.75; - /// Average time for transformation Random Rate in hours constexpr real_t kAverageTimeTransformationRandomRate = 38.6; - /// Standard Deviation for transformation Random Rate in hours constexpr real_t kStandardDeviationTransformationRandomRate = 3.7; +/// Average adhesion time in minutes for Tumor Cell under CAR-T attack before +/// escaping +constexpr real_t kAdhesionTime = 60.0; +/// Min oncoprotein level to be killed by a CAR-T cell +constexpr real_t kOncoproteinLimit = 0.5; +/// Max oncoprotein level +constexpr real_t kOncoproteinSaturation = 2.0; +/// Do not modify this line: difference between saturation and limit +constexpr real_t kOncoproteinDifference = + kOncoproteinSaturation - kOncoproteinLimit; /// volume relaxation rate (min^-1) for each state constexpr real_t kVolumeRelaxarionRateAliveCytoplasm = 0.13 / 60.; @@ -98,8 +106,12 @@ constexpr real_t kThresholdCancerCellType4 = 0.0; /// /// General Hyperparameters /// + /// Seed for random number generation -constexpr int kSeed = 3; +constexpr int kSeed = 1; + +/// Output Performance Statistics +constexpr bool kOutputPerformanceStatistics = true; /// 0.01 minutes time step for substances secretion/consumption constexpr real_t kDtSubstances = 0.01; @@ -115,22 +127,26 @@ constexpr real_t kDt = kDtMechanics; /// computed to avoid errors with fmod constexpr int kStepsPerCycle = kDtCycle / kDt; +/// Output little summary each half a day +constexpr int kOutputCsvInterval = 12 * 60 / kDt; + +/// Total simulation time in minutes (30 days) +constexpr int kTotalMinutesToSimulate = 30 * 24 * 60; + /// Epsilon for avoiding division by 0 constexpr real_t kEpsilon = 1e-10; /// kEpsilon distance constexpr real_t kEpsilonDistance = 1e-5; -/// Output little summary each half a day -constexpr int kOutputCsvInterval = 12 * 60 / kDt; - -/// Total simulation time in minutes (30 days) -constexpr int kTotalMinutesToSimulate = 30 * 24 * 60; /// Length of the bounded space in micrometers constexpr int kBoundedSpaceLength = 1000; /// Initial radius of the spherical tumor (group of cancer cells) in micrometers constexpr real_t kInitialRadiusTumor = 150; +/// Do not modify this line: Twice Pi +constexpr real_t kTwicePi = 2. * Math::kPi; + constexpr real_t kVolumeRelaxarionRateCytoplasmApoptotic = 1.0 / 60.0; constexpr real_t kVolumeRelaxarionRateNucleusApoptotic = 0.35 / 60.0; constexpr real_t kVolumeRelaxarionRateFluidApoptotic = 0.0; @@ -195,16 +211,8 @@ constexpr real_t kDnew = 1.5 * kDtMechanics; /// -0.5) constexpr real_t kDold = -0.5 * kDtMechanics; -/// Large time to avoid division by 0 -constexpr real_t kTimeTooLarge = 1e100; - -/// Minutes in an hour -constexpr real_t kMinutesInAnHour = 60.0; -/// Hours in a day -constexpr real_t kHoursInADay = 24.0; - /// Do not change this line -const size_t gKLengthBoxMechanics = 22; +const real_t gKLengthBoxMechanics = 22; /// Max Distance for considering two cells as neighbours for force calculations /// in μm Do not change this line @@ -216,11 +224,21 @@ const real_t gKSquaredMaxDistanceNeighborsForce = kMaxRelativeAdhesionDistance, 2); +/// Large time to avoid division by 0 +constexpr real_t kTimeTooLarge = 1e100; + +/// Minutes in an hour +constexpr real_t kMinutesInAnHour = 60.0; +/// Hours in a day +constexpr real_t kHoursInADay = 24.0; + /// /// CAR-T Cell Hyperparameters /// constexpr real_t kAverageMaximumTimeUntillApoptosisCart = kDtCycle * 10.0 * 24.0 * 60.0; +/// Default oxygen consumption rate of CAR-T cell +constexpr real_t kDefaultOxygenConsumptionCarT = 1.0; /// Volume parameters /// Default total volume of a new CAR-T cell in μm³ constexpr real_t kDefaultVolumeNewCarTCell = 2494.0; @@ -229,6 +247,75 @@ constexpr real_t kDefaultVolumeNucleusCarTCell = 540.0; /// Default fraction of fluid volume in a new CAR-T cell constexpr real_t kDefaultFractionFluidCarTCell = 0.75; +/// How often a CAR-T cell tries to kill an attached cancer cell +constexpr real_t kKillRateCart = 0.06667; // 1/min +/// How often a CAR-T cell tries to attach to a cancer cell +constexpr real_t kAdhesionRateCart = 0.013; // 1/min +/// Maximum adhesion distance between CAR-T and tumor cells +constexpr real_t kMaxAdhesionDistanceCart = 18.0; // micrometers +/// Minimum adhesion distance between CAR-T and tumor cells +constexpr real_t kMinAdhesionDistanceCart = 14.0; // micrometers + +/// Minimum distance in micrometers from the tumor for spawning CAR-T cells +constexpr real_t kMinimumDistanceCarTFromTumor = 50.0; + +/// Motility parameters +/// Average persistence time before CAR-T cell moves +constexpr real_t kPersistenceTimeCart = 10; // 10 minutes +/// Higher bias (\in [0,1]) makes CAR-T movement more directed toward +/// immunostimulatory factor source; while a bias of 0 makes the movement random +constexpr real_t kMigrationBiasCart = 0.5; +/// Migration speed +constexpr real_t kMigrationSpeedCart = 5.0; +/// Elastic constant +constexpr real_t kElasticConstantCart = 0.01; + +/// Treatment Dosages +/// +/// Specifies the CAR-T cell infusion schedule as a map where: +/// - The key represents the day of treatment (starting from day 0). +/// - The value represents the number of CAR-T cells administered on that day. +/// Example: On day 0 and day 8, 3957 CAR-T cells are introduced (matching the +/// initial tumor cell count). +// FIXME: add this parameters in an external no recompiled file which will also +// avoid cppcoreguidelines-avoid-magic-numbers warning +inline const std::map gKTreatment = { + {0, 3957}, + {8, + 3957}}; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + +/// Do not modify this line: 1-kMigrationBiasCart +constexpr real_t kMigrationOneMinusBiasCart = 1.0 - kMigrationBiasCart; +/// Do not modify this line: probability of a CAR-T cell to migrate in a given +/// mechanical time step +constexpr real_t kMotilityProbability = kDtMechanics / kPersistenceTimeCart; +/// Do not modify this line: probability of a Tumor cell to escape in a given +/// mechanical time step +constexpr real_t kProbabilityEscape = kDtMechanics / (kAdhesionTime + 1e-15); +/// Do not modify this line: Steps in one day +constexpr size_t kStepsOneDay = 24 * 60 / kDt; +/// Do not modify this line: maximum adhesion distance squared +constexpr real_t kSquaredMaxAdhesionDistanceCart = + kMaxAdhesionDistanceCart * kMaxAdhesionDistanceCart; +/// Do not modify this line: difference between min and max adhesion distance +constexpr real_t kDifferenceCartAdhesionDistances = + kMaxAdhesionDistanceCart - kMinAdhesionDistanceCart; + +/// Do not change this line: radius tumor cell +const real_t gKRadiusTumorCell = + std::cbrt(kDefaultVolumeNewTumorCell * 3. / (4. * Math::kPi)); + +/// Do not change this line: radius cart cell +const real_t gKRadiusCarTCell = + std::cbrt(kDefaultVolumeNewCarTCell * 3. / (4. * Math::kPi)); + +// Do not change this line: maximum squared distance to avoid CAR-T pushing +// tumor cells If a CAR-T and a Tumor Cell are closer than this distance, the +// CAR-T cell will only move to the tumor cell with the adhesion forces +//(radiusCART + radiusTumorCell + 0.1 to avoid numerical errors)**2 +const real_t gKMaxSquaredDistanceCartMovingTowardsTumorCell = + std::pow(gKRadiusCarTCell + gKRadiusTumorCell + 1, 2); + } // namespace bdm #endif // TUMOR_HYPERPARAMS_H_ diff --git a/src/tumor_cell.cc b/src/tumor_cell.cc index d0a5baa..b401876 100644 --- a/src/tumor_cell.cc +++ b/src/tumor_cell.cc @@ -26,11 +26,14 @@ #include "core/agent/cell_division_event.h" #include "core/agent/new_agent_event.h" #include "core/container/math_array.h" +#include "core/diffusion/diffusion_grid.h" #include "core/functor.h" #include "core/interaction_force.h" #include "core/real_t.h" #include "core/resource_manager.h" +#include "core/simulation.h" #include "core/util/log.h" +#include "core/util/random.h" #include #include #include @@ -63,16 +66,15 @@ TumorCell::TumorCell(const Real3& position) { (kDefaultVolumeNewTumorCell - kDefaultVolumeNucleusTumorCell) * (1 - kDefaultFractionFluidTumorCell)); - // Set initial oncoproteine level with a truncated normal distribution - SetOncoproteineLevel( + // Set initial oncoprotein level with a truncated normal distribution + SetOncoproteinLevel( SamplePositiveGaussian(kOncoproteinMean, kOncoproteinStandardDeviation)); + ResourceManager* rm = Simulation::GetActive()->GetResourceManager(); // Pointer to oxygen diffusion grid - oxygen_dgrid_ = - Simulation::GetActive()->GetResourceManager()->GetDiffusionGrid("oxygen"); + oxygen_dgrid_ = rm->GetDiffusionGrid("oxygen"); // Pointer to immunostimulatory_factor diffusion grid immunostimulatory_factor_dgrid_ = - Simulation::GetActive()->GetResourceManager()->GetDiffusionGrid( - "immunostimulatory_factor"); + rm->GetDiffusionGrid("immunostimulatory_factor"); // Set state transition random rate SetTransformationRandomRate(); @@ -103,8 +105,8 @@ void TumorCell::Initialize(const NewAgentEvent& event) { oxygen_dgrid_ = mother->oxygen_dgrid_; // Pointer to the immunostimulatory_factor diffusion grid immunostimulatory_factor_dgrid_ = mother->immunostimulatory_factor_dgrid_; - // inherit oncoproteine level from mother cell - this->SetOncoproteineLevel(mother->oncoproteine_level_); + // inherit oncoprotein level from mother cell + this->SetOncoproteinLevel(mother->oncoprotein_level_); // inherit oxygen consumption from mother cell this->SetOxygenConsumptionRate(mother->GetOxygenConsumptionRate()); // inherit immunostimulatory factor secretion rate from mother cell @@ -142,8 +144,8 @@ void TumorCell::Initialize(const NewAgentEvent& event) { } } -void TumorCell::SetOncoproteineLevel(real_t level) { - oncoproteine_level_ = level; +void TumorCell::SetOncoproteinLevel(real_t level) { + oncoprotein_level_ = level; // cell type if (level >= kThresholdCancerCellType1) { // between 1.5 and 2.0 @@ -164,7 +166,7 @@ void TumorCell::SetOncoproteineLevel(real_t level) { } void TumorCell::SetTransformationRandomRate() { - // avoid division by zero + // avoid division by zero with kEpsilon transformation_random_rate_ = 1 / (std::max(SamplePositiveGaussian( kAverageTimeTransformationRandomRate, @@ -267,10 +269,10 @@ Real3 TumorCell::CalculateDisplacement(const InteractionForce* force, uint64_t non_zero_neighbor_forces = 0; if (!IsStatic()) { - auto* ctxt = Simulation::GetActive()->GetExecutionContext(); + ExecutionContext* ctxt = Simulation::GetActive()->GetExecutionContext(); auto calculate_neighbor_forces = L2F([&](Agent* neighbor, real_t /*squared_distance*/) { - auto neighbor_force = force->Calculate(this, neighbor); + Real4 neighbor_force = force->Calculate(this, neighbor); if (neighbor_force[0] != 0 || neighbor_force[1] != 0 || neighbor_force[2] != 0) { non_zero_neighbor_forces++; @@ -351,9 +353,36 @@ void TumorCell::ComputeConstantsConsumptionSecretion() { (immunostimulatory_factor_secretion_rate_); } +void TumorCell::StartApoptosis() { + // If the cell is already dead, do nothing + if (IsDead()) { + return; + } + + // The cell Dies + SetState(TumorCellState::kApoptotic); + + // Reset timer_state + SetTimerState(0); + // Set type to 5 to indicate dead cell + SetType(TumorCellType::kType5); + // Set target volume to 0 (the cell shrinks) + SetTargetCytoplasmSolid(0.0); + SetTargetNucleusSolid(0.0); + SetTargetFractionFluid(0.0); + SetTargetRelationCytoplasmNucleus(0.0); + // Reduce oxygen consumption + SetOxygenConsumptionRate(GetOxygenConsumptionRate() * + kReductionConsumptionDeadCells); + // Stop Immunostimulatory Factor Secretion + SetImmunostimulatoryFactorSecretionRate(0.0); + // Update constants for consumption/secretion differential equation solving + ComputeConstantsConsumptionSecretion(); +} + /// Main behavior executed at each simulation step void StateControlGrowProliferate::Run(Agent* agent) { - auto* sim = Simulation::GetActive(); + Simulation* sim = Simulation::GetActive(); // Run only every kDtCycle minutes, fmod does not work with the type // returned by GetSimulatedTime() if (sim->GetScheduler()->GetSimulatedSteps() % kStepsPerCycle != 0) { @@ -368,7 +397,7 @@ void StateControlGrowProliferate::Run(Agent* agent) { // Oxygen levels const Real3 current_position = cell->GetPosition(); // Pointer to the oxygen diffusion grid - auto* oxygen_dgrid = cell->GetOxygenDiffusionGrid(); + DiffusionGrid* oxygen_dgrid = cell->GetOxygenDiffusionGrid(); const real_t oxygen_level = oxygen_dgrid->GetValue(current_position); switch (cell->GetState()) { @@ -435,7 +464,7 @@ void StateControlGrowProliferate::Run(Agent* agent) { // duration transition) if (kTimeLysis < cell->GetTimerState()) { // remove the cell from the simulation - auto* ctxt = sim->GetExecutionContext(); + ExecutionContext* ctxt = sim->GetExecutionContext(); ctxt->RemoveAgent(agent->GetUid()); } break; @@ -454,7 +483,7 @@ void StateControlGrowProliferate::Run(Agent* agent) { // duration transition) if (kTimeApoptosis < cell->GetTimerState()) { // remove the cell from the simulation - auto* ctxt = sim->GetExecutionContext(); + ExecutionContext* ctxt = sim->GetExecutionContext(); ctxt->RemoveAgent(agent->GetUid()); } break; @@ -491,11 +520,10 @@ void StateControlGrowProliferate::ManageLivingCell(TumorCell* cell, if (oxygen_level < kOxygenLimitForProliferation) { multiplier = 0.0; } - // Calculate the rate of state change based on oxygen level and oncoproteine + // Calculate the rate of state change based on oxygen level and oncoprotein // (min^-1) const real_t final_rate_transition = cell->GetTransformationRandomRate() * - multiplier * - cell->GetOncoproteineLevel(); + multiplier * cell->GetOncoproteinLevel(); // Calculate the time to wait (in minutes) real_t time_to_wait = kTimeTooLarge; @@ -534,8 +562,8 @@ bool StateControlGrowProliferate::ShouldEnterNecrosis(real_t oxygen_level, const real_t probability_necrosis = kDtCycle * kMaximumNecrosisRate * multiplier; - auto* sim = Simulation::GetActive(); - auto* random = sim->GetRandom(); + Simulation* sim = Simulation::GetActive(); + Random* random = sim->GetRandom(); const bool enter_necrosis = random->Uniform(0, 1) < probability_necrosis; // If the random number is less than the probability, enter necrosis if (enter_necrosis) { diff --git a/src/tumor_cell.h b/src/tumor_cell.h index d4b6dea..55b7fb3 100644 --- a/src/tumor_cell.h +++ b/src/tumor_cell.h @@ -107,8 +107,8 @@ class TumorCell : public Cell { void SetTimerState(real_t timer_state) { timer_state_ = timer_state; } real_t GetTimerState() const { return timer_state_; } - void SetOncoproteineLevel(real_t level); - real_t GetOncoproteineLevel() const { return oncoproteine_level_; } + void SetOncoproteinLevel(real_t level); + real_t GetOncoproteinLevel() const { return oncoprotein_level_; } void SetFluidFraction(real_t fluid_fraction) { fluid_fraction_ = fluid_fraction; @@ -171,6 +171,8 @@ class TumorCell : public Cell { real_t GetTargetTotalVolume() const; + bool IsDead() const { return type_ == TumorCellType::kType5; } + int GetTypeAsInt() const { return static_cast(type_); } /// Returns the diffusion grid for oxygen @@ -225,6 +227,12 @@ class TumorCell : public Cell { /// equations. void ComputeConstantsConsumptionSecretion(); + /// Start apoptosis + /// + /// This function is called when the tumor cell is induced to apoptosis by a + /// CAR-T cell. + void StartApoptosis(); + private: /// Current state of the tumor cell TumorCellState state_ = TumorCellState::kAlive; @@ -239,7 +247,7 @@ class TumorCell : public Cell { DiffusionGrid* immunostimulatory_factor_dgrid_ = nullptr; /// Level of oncoprotein expression - real_t oncoproteine_level_ = 0.0; + real_t oncoprotein_level_ = 0.0; /// Transition random rate between states: /// Affects the probability of transitioning and depends on the individual @@ -276,7 +284,7 @@ class TumorCell : public Cell { TumorCellType type_ = TumorCellType::kType0; /// Velocity of the cell in the previous time step - Real3 older_velocity_; + Real3 older_velocity_ = {0, 0, 0}; /// Rate of oxygen consumption by the cell real_t oxygen_consumption_rate_ = 0.0; diff --git a/src/utils_aux.cc b/src/utils_aux.cc index 1aed07e..36906ef 100644 --- a/src/utils_aux.cc +++ b/src/utils_aux.cc @@ -3,7 +3,13 @@ * Melina Luque * * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. + * you may not use this file e // Count total cells, tumor cells of each + type and tumor radius size_t total_num_tumor_cells = 0; size_t + num_tumor_cells_type1 = 0; size_t num_tumor_cells_type2 = 0; size_t + num_tumor_cells_type3 = 0; size_t num_tumor_cells_type4 = 0; size_t + num_tumor_cells_type5_dead = 0; size_t num_alive_cart = 0; real_t tumor_radius + = 0.0; real_t average_oncoprotein = 0.0; real_t average_oxygen_cancer_cells = + 0.0;in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 @@ -20,17 +26,25 @@ */ #include "utils_aux.h" +#include "cart_cell.h" #include "hyperparams.h" #include "tumor_cell.h" #include "core/agent/agent.h" #include "core/container/math_array.h" +#include "core/diffusion/diffusion_grid.h" +#include "core/param/param.h" #include "core/real_t.h" #include "core/resource_manager.h" +#include "core/scheduler.h" +#include "core/simulation.h" #include "core/util/math.h" +#include "core/util/random.h" #include #include +#include #include #include +#include #include #include @@ -39,7 +53,7 @@ namespace bdm { // Samples a Gaussian value with given mean and standard deviation but all // negative values are mapped to zero real_t SamplePositiveGaussian(float mean, float sigma) { - auto* random = Simulation::GetActive()->GetRandom(); + Random* random = Simulation::GetActive()->GetRandom(); real_t value = random->Gaus(mean, sigma); if (value < 0.) { value = 0.; @@ -50,7 +64,7 @@ real_t SamplePositiveGaussian(float mean, float sigma) { std::vector CreateSphereOfTumorCells(real_t sphere_radius) { // V = (4/3)*pi*r^3 = (pi/6)*diameter^3 const real_t cell_radius = - std::cbrt(kDefaultVolumeNewTumorCell * 6 / Math::kPi) / 2; + std::cbrt(kDefaultVolumeNewTumorCell * 6 / Math::kPi) / kHalf; std::vector positions; @@ -88,42 +102,59 @@ std::vector CreateSphereOfTumorCells(real_t sphere_radius) { return positions; } -// Function to compute the number of tumor cells of each type and the radius of -// the tumor -std::tuple -ComputeNumberTumorCellsAndRadius() { - auto* rm = Simulation::GetActive()->GetResourceManager(); - size_t total_num_tumor_cells = 0; - size_t num_tumor_cells_type1 = 0; - size_t num_tumor_cells_type2 = 0; - size_t num_tumor_cells_type3 = 0; - size_t num_tumor_cells_type4 = 0; - size_t num_tumor_cells_type5_dead = 0; +// Function to compute the number of tumor cells of each type, the radius of the +// tumor, the average oncoprotein level and the average oxygen level of the +// cancer cells. +std::tuple +AnalyzeTumor() { + ResourceManager* rm = Simulation::GetActive()->GetResourceManager(); + DiffusionGrid* oxygen_dgrid = rm->GetDiffusionGrid("oxygen"); + + int total_num_tumor_cells = 0; + int num_tumor_cells_type1 = 0; + int num_tumor_cells_type2 = 0; + int num_tumor_cells_type3 = 0; + int num_tumor_cells_type4 = 0; + int num_tumor_cells_type5_dead = 0; + int num_alive_cart = 0; real_t max_dist_sq = 0.0; + real_t acumulator_oncoprotein = 0.0; + real_t acumulator_oxygen_cancer_cells = 0.0; rm->ForEachAgent([&](const Agent* agent) { if (const auto* tumor_cell = dynamic_cast(agent)) { total_num_tumor_cells++; - const auto& pos = agent->GetPosition(); + const Real3& pos = agent->GetPosition(); + + // Accumulate oxygen level for average calculation + acumulator_oxygen_cancer_cells += oxygen_dgrid->GetValue(pos); + + // computing tumor radius const real_t dist_sq = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]; if (dist_sq > max_dist_sq) { max_dist_sq = dist_sq; } - // Count tumor cells by type + // Count tumor cells by type and accumulate oncoprotein levels if the cell + // is alive switch (tumor_cell->GetType()) { case TumorCellType::kType1: + acumulator_oncoprotein += tumor_cell->GetOncoproteinLevel(); num_tumor_cells_type1++; break; case TumorCellType::kType2: + acumulator_oncoprotein += tumor_cell->GetOncoproteinLevel(); num_tumor_cells_type2++; break; case TumorCellType::kType3: + acumulator_oncoprotein += tumor_cell->GetOncoproteinLevel(); num_tumor_cells_type3++; break; case TumorCellType::kType4: + acumulator_oncoprotein += tumor_cell->GetOncoproteinLevel(); num_tumor_cells_type4++; break; case TumorCellType::kType5: @@ -132,59 +163,182 @@ ComputeNumberTumorCellsAndRadius() { default: break; } + } else if (const auto* cart_cell = dynamic_cast(agent)) { + if (cart_cell->GetState() == CarTCellState::kAlive) { + num_alive_cart++; + } } }); + const int total_alive_cells = num_tumor_cells_type1 + num_tumor_cells_type2 + + num_tumor_cells_type3 + num_tumor_cells_type4; + const real_t average_oncoprotein = + (total_alive_cells > 0) ? (acumulator_oncoprotein / total_alive_cells) + : 0.0; + const real_t average_oxygen_cancer_cells = + (total_num_tumor_cells > 0) + ? (acumulator_oxygen_cancer_cells / total_num_tumor_cells) + : 0.0; return {total_num_tumor_cells, num_tumor_cells_type1, num_tumor_cells_type2, num_tumor_cells_type3, num_tumor_cells_type4, num_tumor_cells_type5_dead, - std::sqrt(max_dist_sq)}; + num_alive_cart, std::sqrt(max_dist_sq), + average_oncoprotein, average_oxygen_cancer_cells}; } // Function to output summary CSV void OutputSummary::operator()() { - auto* scheduler = Simulation::GetActive()->GetScheduler(); - auto current_step = scheduler->GetSimulatedSteps(); + Simulation* simulation = Simulation::GetActive(); + Scheduler* scheduler = simulation->GetScheduler(); + const uint64_t current_step = scheduler->GetSimulatedSteps(); if (current_step % frequency_ == 0) { std::ofstream file("output/final_data.csv", std::ios::app); if (file.is_open()) { - // Header for the CSV file if (current_step == 0) { - file << "total_days,total_hours,total_minutes,tumor_radius,num_cells," - "num_tumor_cells,tumor_cells_type1,tumor_cells_type2,tumor_" - "cells_type3,tumor_cells_type4,tumor_cells_type5_dead\n"; + file + << "total_days,total_hours,total_minutes,tumor_radius,num_cells," + "num_tumor_cells,tumor_cells_type1,tumor_cells_type2,tumor_" + "cells_type3,tumor_cells_type4,tumor_cells_type5_dead,num_alive_" + "cart,average_oncoprotein,average_oxygen_cancer_cells\n"; // Header + // for + // CSV + // file } // Calculate time in days, hours, minutes - const double total_minutes = - Simulation::GetActive()->GetScheduler()->GetSimulatedTime(); + const double total_minutes = static_cast(current_step) * kDt; const double total_hours = total_minutes / kMinutesInAnHour; const double total_days = total_hours / kHoursInADay; // Count total cells, tumor cells of each type and tumor radius - size_t total_num_tumor_cells = 0; - size_t num_tumor_cells_type1 = 0; - size_t num_tumor_cells_type2 = 0; - size_t num_tumor_cells_type3 = 0; - size_t num_tumor_cells_type4 = 0; - size_t num_tumor_cells_type5_dead = 0; + int total_num_tumor_cells = 0; + int num_tumor_cells_type1 = 0; + int num_tumor_cells_type2 = 0; + int num_tumor_cells_type3 = 0; + int num_tumor_cells_type4 = 0; + int num_tumor_cells_type5_dead = 0; + int num_alive_cart = 0; real_t tumor_radius = 0.0; - + real_t average_oncoprotein = 0.0; + real_t average_oxygen_cancer_cells = 0.0; std::tie(total_num_tumor_cells, num_tumor_cells_type1, num_tumor_cells_type2, num_tumor_cells_type3, num_tumor_cells_type4, num_tumor_cells_type5_dead, - tumor_radius) = ComputeNumberTumorCellsAndRadius(); + num_alive_cart, tumor_radius, average_oncoprotein, + average_oxygen_cancer_cells) = AnalyzeTumor(); + size_t total_num_cells = simulation->GetResourceManager()->GetNumAgents(); + + // If a dosage is administred this exact time the numbers are not seen in + // the resource manager yet because of how BioDynaMo is built + // therefore we need to add the just added new cells to the statistics + // here. + const auto current_day_int = static_cast(total_days); + if (current_step % kStepsOneDay == 0 && + gKTreatment.find(current_day_int) != gKTreatment.end()) { + const size_t just_spawned_cells = gKTreatment.at(current_day_int); + total_num_cells += just_spawned_cells; + num_alive_cart += static_cast(just_spawned_cells); + } // Write data to CSV file file << total_days << "," << total_hours << "," << total_minutes << "," - << tumor_radius << "," - << Simulation::GetActive()->GetResourceManager()->GetNumAgents() - << "," << total_num_tumor_cells << "," << num_tumor_cells_type1 - << "," << num_tumor_cells_type2 << "," << num_tumor_cells_type3 - << "," << num_tumor_cells_type4 << "," << num_tumor_cells_type5_dead - << "\n"; + << tumor_radius << "," << total_num_cells << "," + << total_num_tumor_cells << "," << num_tumor_cells_type1 << "," + << num_tumor_cells_type2 << "," << num_tumor_cells_type3 << "," + << num_tumor_cells_type4 << "," << num_tumor_cells_type5_dead << "," + << num_alive_cart << "," << average_oncoprotein << "," + << average_oxygen_cancer_cells << "\n"; } } } +// Function to spawn CAR-T cell dosages +void SpawnCart::operator()() { + Simulation* simulation = Simulation::GetActive(); + Scheduler* scheduler = simulation->GetScheduler(); + const uint64_t current_step = scheduler->GetSimulatedSteps(); + // This function only executes each day + if (current_step % frequency_ != 0) { + return; + } + // See if there is any dosage to apply in this day + const auto current_day = + static_cast(static_cast(current_step) * kDt / + (kMinutesInAnHour * kHoursInADay)); + size_t cells_to_spawn = 0; + + if (gKTreatment.find(current_day) != gKTreatment.end()) { + cells_to_spawn = gKTreatment.at(current_day); + } + + // if there are cells to spawn in the treatment + if (cells_to_spawn > 0) { + // compute tumor radius + ResourceManager* rm = simulation->GetResourceManager(); + real_t max_dist_sq = 0.0; + + rm->ForEachAgent([&](const Agent* agent) { + if (const auto* cancer_cell = dynamic_cast(agent)) { + const Real3& pos = cancer_cell->GetPosition(); + const real_t dist_sq = + pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]; + if (dist_sq > max_dist_sq) { + max_dist_sq = dist_sq; + } + } + }); + + // the car-t spawns at least kMinimumDistanceCarTFromTumor micrometers away + // from the tumor + real_t minimum_squared_radius = + std::sqrt(max_dist_sq) + kMinimumDistanceCarTFromTumor; + minimum_squared_radius *= minimum_squared_radius; + + // for generating car-t positions + Random* rng = simulation->GetRandom(); + const Param* param = simulation->GetParam(); + const real_t min_b = param->min_bound; + const real_t max_b = param->max_bound; + real_t px = 0.0; + real_t py = 0.0; + real_t pz = 0.0; + real_t radi_sq = 0.0; + + ExecutionContext* ctxt = simulation->GetExecutionContext(); + + for (unsigned int i = 0; i < cells_to_spawn; i++) { + // look for a valid position + while (true) { + px = rng->Uniform(min_b, max_b); + py = rng->Uniform(min_b, max_b); + pz = rng->Uniform(min_b, max_b); + radi_sq = px * px + py * py + pz * pz; + if (radi_sq >= minimum_squared_radius) { + break; + } + } + // spawn CAR-T + std::unique_ptr cart = + std::make_unique(Real3{px, py, pz}); + std::unique_ptr state_control = + std::make_unique(); + cart->AddBehavior(state_control.release()); + ctxt->AddAgent(cart.release()); + } + } +} + +// Function to generate a random direction unitary vector +Real3 GenerateRandomDirection() { + Random* rnd = Simulation::GetActive()->GetRandom(); + // Generate a random direction vector + const double theta = kTwicePi * rnd->Uniform(0.0, 1.0); + const double phi = Math::kPi * rnd->Uniform(0.0, 1.0); + + const double sin_phi = std::sin(phi); + const double cos_phi = std::cos(phi); + + return {sin_phi * std::cos(theta), sin_phi * std::sin(theta), cos_phi}; +} + } // namespace bdm diff --git a/src/utils_aux.h b/src/utils_aux.h index 2a08fa9..85e4b8c 100644 --- a/src/utils_aux.h +++ b/src/utils_aux.h @@ -59,24 +59,61 @@ std::vector CreateSphereOfTumorCells(real_t sphere_radius); /// Compute tumor statistics and characteristics /// /// Analyzes the current tumor population to compute the number of tumor cells -/// of each type and the overall radius of the tumor mass. +/// of each type and the overall radius of the tumor mass. In addition, it +/// computes the average oncoprotein level and oxygen across all tumor cells. /// /// @return Tuple containing: +/// - Total number of tumor cells /// - Number of type 1 tumor cells (most aggressive) /// - Number of type 2 tumor cells /// - Number of type 3 tumor cells /// - Number of type 4 tumor cells (least aggressive) /// - Number of type 5 tumor cells (dead) -/// - Total number of tumor cells +/// - Number of living CAR-T cells (not apoptotic) /// - Current tumor radius in micrometers -std::tuple -ComputeNumberTumorCellsAndRadius(); +/// - Average oncoprotein level across all tumor cells +/// - Average oxygen level across all tumor cells +std::tuple +AnalyzeTumor(); + +/// Generates a random direction unitary vector +/// +/// @return A 3D vector representing a random direction +Real3 GenerateRandomDirection(); + +/// Spawns a Dosage of CAR-T cells around the tumor +/// +/// Called automatically by the simulation scheduler at the specified frequency. +/// Spawns a dosage of CAR-T cells around the tumor following the map +/// gKTreatment where +/// - The key represents the day of treatment (starting from day 0). +/// - The value represents the number of CAR-T cells administered on that day. +struct SpawnCart : public StandaloneOperationImpl { + // NOLINTNEXTLINE(cppcoreguidelines-owning-memory) + BDM_OP_HEADER(SpawnCart); + + public: + void SetFrequency(uint64_t frequency) { frequency_ = frequency; } + [[nodiscard]] uint64_t GetFrequency() const { return frequency_; } + + private: + /// Frequency of output (every N simulation steps) + uint64_t frequency_ = 1; + + void operator()() override; +}; + +/// Register SpawnCart operation with CPU as compute target +// NOLINTNEXTLINE(cppcoreguidelines-owning-memory) +inline BDM_REGISTER_OP(SpawnCart, "SpawnCart", kCpu); -/// Operation for outputting simulation summary data to CSV files +/// Operation for outputting simulation summary data to a CSV file /// /// This operation collects and outputs summary statistics about the simulation -/// state to CSV files for post-processing and analysis. It includes information -/// about cell populations, tumor characteristics, and other relevant metrics. +/// state to output/final_data.csv for post-processing and analysis. It includes +/// information about cell populations, tumor characteristics, and other +/// relevant metrics. struct OutputSummary : public StandaloneOperationImpl { // NOLINTNEXTLINE(cppcoreguidelines-owning-memory) BDM_OP_HEADER(OutputSummary); @@ -89,12 +126,6 @@ struct OutputSummary : public StandaloneOperationImpl { /// Frequency of output (every N simulation steps) uint64_t frequency_ = 1; - /// Collects current simulation data and writes it to CSV files - /// - /// Called automatically by the simulation scheduler at the specified - /// frequency. Gathers statistics about cell populations, tumor radius, and - /// other metrics, then outputs them to appropriately named CSV files for - /// analysis. void operator()() override; }; diff --git a/test/test-suite.cc b/test/test-suite.cc index e33f20d..60d3e8b 100644 --- a/test/test-suite.cc +++ b/test/test-suite.cc @@ -20,6 +20,7 @@ */ #include "biodynamo.h" +#include "core/resource_manager.h" #include // Googletest in combination with the provided CMakeLists.txt allows you to @@ -52,7 +53,7 @@ TEST(AgentTest, AddAgentsToSimulation) { Simulation simulation(TEST_NAME); // Add some cells to the simulation - auto* rm = simulation.GetResourceManager(); + ResourceManager* rm = simulation.GetResourceManager(); uint8_t expected_no_cells{20}; for (int i = 0; i < expected_no_cells; i++) { auto* cell = new Cell(30);