diff --git a/.clang-tidy b/.clang-tidy index f25eb79..68dc3bd 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -11,7 +11,9 @@ Checks: > openmp-* performance-*, portability-*, - readability-* + readability-*, + -modernize-use-trailing-return-type, + -readability-identifier-length CheckOptions: - key: google-readability-braces-around-statements.ShortStatementLines diff --git a/.gitignore b/.gitignore index d4fb281..7d30e11 100644 --- a/.gitignore +++ b/.gitignore @@ -39,3 +39,86 @@ # debug information files *.dwo + +build/ +output/ + +**/.DS_Store + +#----------------- intellij ----------------- +.idea +cx3d-cpp.iml + +#----------------- eclipse ------------------ +.project +.cproject +.settings/ + +#----------------- VS Code ------------------ +.vscode/* + +#----------------- GNU global ------------------ +GPATH +GRTAGS +GTAGS + +#------------------- C++ -------------------- +# Compiled Object files +*.slo +*.lo +*.o +*.obj + +# Precompiled Headers +*.gch +*.pch + +# Compiled Dynamic libraries +*.so +*.dylib +*.dll +*.jnilib + +# Fortran module files +*.mod + +# Compiled Static libraries +*.lai +*.la +*.a +*.lib + +# Executables +*.exe +*.out +*.app + +#------------------- CMake -------------------- +CMakeCache.txt +CMakeFiles +CMakeScripts +Makefile +cmake_install.cmake +install_manifest.txt + +#------------------- Doxygen files -------------------- +Doxyfile +doc/html +doc/latex + +#------------------- Paraview -------------------- +*.vtu +*.vti +*.pvtu +*.pvti + + +#Debugging files +*.csv + +# Result files +*.xyz +*.dat + +# Draft code for comparing models +draft_code_my_own_analysis/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..c56ee3b --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,46 @@ +# ----------------------------------------------------------------------------- +# +# Copyright (C) 2021 CERN & University of Surrey for the benefit of the +# BioDynaMo collaboration. All Rights Reserved. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# +# See the LICENSE file distributed with this work for details. +# See the NOTICE file distributed with this work for additional information +# regarding copyright ownership. +# +# ----------------------------------------------------------------------------- +cmake_minimum_required(VERSION 3.19.3) +project(cart_tumor) + + +# BioDynaMo curretly uses the C++17 standard. +set(CMAKE_CXX_STANDARD 17) + +# Use BioDynaMo in this project. +find_package(BioDynaMo REQUIRED) + +# See UseBioDynaMo.cmake in your BioDynaMo build folder for details. +# Note that BioDynaMo provides gtest header/libraries in its include/lib dir. +include(${BDM_USE_FILE}) + +# Consider all files in src/ for BioDynaMo simulation. +include_directories("src") +file(GLOB_RECURSE PROJECT_HEADERS src/*.h) +file(GLOB_RECURSE PROJECT_SOURCES src/*.cc) + +bdm_add_executable(${CMAKE_PROJECT_NAME} + HEADERS ${PROJECT_HEADERS} + SOURCES ${PROJECT_SOURCES} + LIBRARIES ${BDM_REQUIRED_LIBRARIES}) + +# Consider all files in test/ for GoogleTests. +include_directories("test") +file(GLOB_RECURSE TEST_SOURCES test/*.cc) +file(GLOB_RECURSE TEST_HEADERS test/*.h) + +bdm_add_test(${CMAKE_PROJECT_NAME}-test + SOURCES ${TEST_SOURCES} + HEADERS ${TEST_HEADERS} + LIBRARIES ${BDM_REQUIRED_LIBRARIES} ${CMAKE_PROJECT_NAME}) diff --git a/README.md b/README.md index 97f937b..42c47e8 100644 --- a/README.md +++ b/README.md @@ -1 +1,130 @@ -# CARTopiaX \ No newline at end of file +# CARTopiaX + +This repository provides an **agent-based simulation** of tumour-derived organoids and their interaction with **CAR T-cell therapy**. +Developed as part of Google Summer of Code 2025, the project is released under the **Apache License 2.0**. + +The simulation integrates computational modeling and biological insights to explore tumour–immune dynamics and assess treatment outcomes under various scenarios. + +--- + +## Table of Contents + +1. [Project Overview](#project-overview) +2. [Dependencies](#dependencies) +3. [Installation](#installation) +4. [Building the Simulation](#building-the-simulation) +5. [Running the Simulation](#running-the-simulation) +6. [Visualizing Results](#results) +7. [Acknowledgments](#acknowledgments) +8. [License](#license) + + +--- + +## Project Overview + +This project implements an **agent-based model** to simulate the behavior of *tumour-derived organoids* (lab-grown tumour models) and their response to **CAR T-cell therapy**. + +With this simulation, researchers can: +- Recreate in vitro conditions for tumour growth. +- Introduce CAR T-cells and analyze their effectiveness in solid tumor microenvironments. +- Explore different treatment strategies and parameter variations. +- Evaluate outcomes such as tumour reduction, elimination, or relapse risk. + +By adjusting biological and therapeutic parameters, the model enables **in silico experimentation** to complement laboratory research. + +--- + +## Dependencies + +- [BioDynaMo](https://biodynamo.org/) (tested with version 1.05.132) +- CMake ≥ 3.13 +- GCC or Clang with C++17 support +- GoogleTest (for unit testing) + +**Note:** Ensure BioDynaMo is installed and sourced before running the simulation. + +--- + +## Installation + +Clone the repository: +```bash +git clone https://github.com/compiler-research/CARTopiaX.git +cd CARTopiaX + +``` + +--- + +## Building the Simulation + +**Option 1:** +Use BioDynaMo’s build system: +```bash +biodynamo build +``` + +**Option 2:** +Manual build: +```bash +mkdir build && cd build +cmake .. +make -j +``` + +--- + +## Running the Simulation + +After building, run the simulation using one of the following methods: + +**Option 1:** +With BioDynaMo: +```bash +biodynamo run +``` + +**Option 2:** +Directly from the build directory: +```bash + +./build/CARTopiaX +``` + +--- + +## Results + +Data about tumor growth and diffrent types of cell populations is output in ./output/final_data.csv + +To visualize the results in paraview use: +```bash +paraview ./output/CARTopiaX/CARTopiaX.pvsm + +``` + +--- + +## Acknowledgments + +This project builds upon the BioDynaMo simulation framework. + +> Lukas Breitwieser, Ahmad Hesam, Jean de Montigny, Vasileios Vavourakis, Alexandros Iosif, Jack Jennings, Marcus Kaiser, Marco Manca, Alberto Di Meglio, Zaid Al-Ars, Fons Rademakers, Onur Mutlu, Roman Bauer. +> *BioDynaMo: a modular platform for high-performance agent-based simulation*. +> Bioinformatics, Volume 38, Issue 2, January 2022, Pages 453–460. +> [https://doi.org/10.1093/bioinformatics/btab649](https://doi.org/10.1093/bioinformatics/btab649) + +Some of the mathematical models and solver implementations are based on the research of +Luciana Melina Luque and collaborators, as described in: + +> Luque, L.M., Carlevaro, C.M., Rodriguez-Lomba, E. et al. +> *In silico study of heterogeneous tumour-derived organoid response to CAR T-cell therapy*. +> Scientific Reports 14, 12307 (2024). +> [https://doi.org/10.1038/s41598-024-63125-5](https://doi.org/10.1038/s41598-024-63125-5) + +--- + +## License + +This project is licensed under the Apache License 2.0. See the [LICENSE](LICENSE) file for details. \ No newline at end of file diff --git a/bdm.toml b/bdm.toml new file mode 100644 index 0000000..5557514 --- /dev/null +++ b/bdm.toml @@ -0,0 +1,14 @@ +[visualization] +export = true +interval = 7200 + +[[visualize_agent]] +name = "TumorCell" +additional_data_members = [ + {name = "diameter_", function = "diameter_"}, + {name = "volume_", function = "volume_"}, + {name = "cell_type", function = "GetTypeAsInt"} +] + +[[visualize_diffusion]] +name = "oxygen" diff --git a/src/cart_cell.cc b/src/cart_cell.cc new file mode 100644 index 0000000..b18c66d --- /dev/null +++ b/src/cart_cell.cc @@ -0,0 +1,310 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#include +#include +#include +#include +#include + +#include "core/agent/agent.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/util/log.h" + +#include "cart_cell.h" +#include "hyperparams.h" +#include "tumor_cell.h" +#include "utils_aux.h" + +namespace bdm { + +CartCell::CartCell(const Real3& position) { + SetPosition(position); + SetVolume(kDefaultVolumeNewCartCell); + const ResourceManager& rm = *Simulation::GetActive()->GetResourceManager(); + oxygen_dgrid_ = rm.GetDiffusionGrid("oxygen"); + immunostimulatory_factor_dgrid_ = + rm.GetDiffusionGrid("immunostimulatory_factor"); + ComputeConstantsConsumptionSecretion(); +} + +// Cart cells can move if they are alive and not attached to a tumor cell +bool CartCell::DoesCellMove() { + return (state_ == CartCellState::kAlive && !attached_to_tumor_cell_); +} + +real_t CartCell::GetTargetTotalVolume() const { + return GetTargetNucleusSolid() * (1 + GetTargetRelationCytoplasmNucleus()) / + (1 - GetTargetFractionFluid()); +} + +// This method explicitly solves the system of exponential relaxation +// differential equation using a discrete update step. It is used to shrink the +// volume (and proportions) smoothly toward a desired target volume over time +// whe the cell is apoptotic. The relaxations rate controls the speed of +// convergence +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +void CartCell::ChangeVolumeExponentialRelaxationEquation( + real_t relaxation_rate_cytoplasm, real_t relaxation_rate_nucleus, // NOLINT + real_t relaxation_rate_fluid) { + // Exponential relaxation towards the target volume + const real_t current_total_volume = GetVolume(); + const real_t fluid_fraction = GetFluidFraction(); + const real_t nuclear_volume = GetNuclearVolume(); + + const real_t current_nuclear_solid = nuclear_volume * (1 - fluid_fraction); + const real_t current_cytoplasm_solid = + (current_total_volume - nuclear_volume) * (1 - fluid_fraction); + + const real_t current_fluid = fluid_fraction * current_total_volume; + + // Update fluid volume + real_t new_fluid = + current_fluid + + kDtCycle * relaxation_rate_fluid * + (GetTargetFractionFluid() * current_total_volume - current_fluid); + // Clamp to zero to prevent negative volumes + if (new_fluid < 0.0) { + new_fluid = 0.0; + } + + const real_t nuclear_fluid = + new_fluid * (nuclear_volume / current_total_volume); + // real_t cytoplasm_fluid = new_fluid - nuclear_fluid; + + real_t nuclear_solid = current_nuclear_solid + + kDtCycle * relaxation_rate_nucleus * + (GetTargetNucleusSolid() - current_nuclear_solid); + // Clamp to zero to prevent negative volumes + if (nuclear_solid < 0.0) { + nuclear_solid = 0.0; + } + + const real_t target_cytoplasm_solid = + GetTargetRelationCytoplasmNucleus() * GetTargetNucleusSolid(); + real_t cytoplasm_solid = + current_cytoplasm_solid + + kDtCycle * relaxation_rate_cytoplasm * + (target_cytoplasm_solid - current_cytoplasm_solid); + // Clamp to zero to prevent negative volumes + if (cytoplasm_solid < 0.0) { + cytoplasm_solid = 0.0; + } + + const real_t new_total_solid = nuclear_solid + cytoplasm_solid; + + const real_t total_nuclear = nuclear_solid + nuclear_fluid; + + // real_t total_cytoplasm= cytoplasm_solid + cytoplasm_fluid; + + const real_t new_volume = new_total_solid + new_fluid; + + // Avoid division by zero + const real_t new_fraction_fluid = new_fluid / (kEpsilon + new_volume); + + // Update the cell's properties + // if the volume has changed + if (new_volume != current_total_volume) { + SetVolume(new_volume); + // Update constants for all ConsumptionSecretion of Oxygen and + // Immunostimulatory Factors + ComputeConstantsConsumptionSecretion(); + } + SetFluidFraction(new_fraction_fluid); + SetNuclearVolume(total_nuclear); +} + +// compute Displacement +Real3 CartCell::CalculateDisplacement(const InteractionForce* force, + real_t squared_radius, real_t /*dt*/) { + // 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 + // cleaner code instead of hardcoding it here + squared_radius = gKSquaredMaxDistanceNeighborsForce; + + // the physics force to move the point mass + + 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 + + uint64_t non_zero_neighbor_forces = 0; // NOLINT + if (!IsStatic()) { + auto* ctxt = Simulation::GetActive()->GetExecutionContext(); + auto calculate_neighbor_forces = + 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]; + } + }); + ctxt->ForEachNeighbor(calculate_neighbor_forces, *this, squared_radius); + + if (non_zero_neighbor_forces > 1) { + SetStaticnessNextTimestep(false); + } + } + + // 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; + + older_velocity_ = translation_velocity_on_point_mass; + + // Displacement + return movement_at_next_step; +} + +// Compute new oxygen or immunostimulatory factor concentration after +// consumption/ secretion +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +real_t CartCell::ConsumeSecreteSubstance(int substance_id, + real_t old_concentration) { + real_t res = NAN; // NOLINT + if (substance_id == oxygen_dgrid_->GetContinuumId()) { + // consuming oxygen + res = (old_concentration + constant1_oxygen_) / constant2_oxygen_; + } else if (substance_id == + immunostimulatory_factor_dgrid_->GetContinuumId()) { + // This point should never be reached + res = old_concentration; + } else { + throw std::invalid_argument("Unknown substance id: " + + std::to_string(substance_id)); + } + return res; +} + +// Recompute Consumption constants whenever oxygen_consumption_rate_ or the +// volume changes +void CartCell::ComputeConstantsConsumptionSecretion() { + // constant1_= dt · (V_k / V_voxel) · S_k · ρ*_k) + // constant2_ = [1 + dt · (V_k / V_voxel) · (S_k + U_k)] + // where: + // S_k = secretion rate of cell k + // U_k = uptake (consumption) rate of cell k + // ρ*_k = saturation (target) density for secretion + // V_k = volume of the cell k + // V_voxel = volume of the voxel containing the cell + // dt = simulation time step + const real_t volume = GetVolume(); + // compute the constants for the differential equation explicit solution: for + // oxygen and immunostimulatory factor + // dt*(cell_volume/voxel_volume)*quantity_secretion*substance_saturation = dt + // · (V_k / V_voxel) · S_k · ρ*_k) + constant1_oxygen_ = 0.; + // 1 + dt*(cell_volume/voxel_volume)*(quantity_secretion + + // quantity_consumption ) = [1 + dt · (V_k / V_voxel) · (S_k + U_k)] + // Scale by the volume of the cell in the Voxel and time step + constant2_oxygen_ = + 1 + kDtSubstances * (volume / kVoxelVolume) * (oxygen_consumption_rate_); +} + +/// Main behavior executed at each simulation step +void StateControlCart::Run(Agent* agent) { + auto* sim = Simulation::GetActive(); + if (sim->GetScheduler()->GetSimulatedSteps() % kStepsPerCycle != 0) { + return; + } // Run only every kDtCycle minutes, fmod does not work with the type + // returned by GetSimulatedTime() + + if (auto* cell = dynamic_cast(agent)) { + switch (cell->GetState()) { + case CartCellState::kAlive: { // the cell is growing to real_t its size + // before mitosis + + if (sim->GetRandom()->Uniform(1.0) < + kDtCycle / + std::max( + cell->GetCurrentLiveTime(), + kEpsilon)) { // Probability of death= 1/CurrentLiveTime, + // avoiding division by 0 + // the cell Dies + cell->SetState(CartCellState::kApoptotic); + cell->SetTimerState(0); // Reset timer_state, it should be 0 anyway + // Set target volume to 0 (the cell will shrink) + cell->SetTargetCytoplasmSolid(0.0); + cell->SetTargetNucleusSolid(0.0); + cell->SetTargetFractionFluid(0.0); + cell->SetTargetRelationCytoplasmNucleus(0.0); + // Reduce oxygen consumption + cell->SetOxygenConsumptionRate(cell->GetOxygenConsumptionRate() * + kReductionConsumptionDeadCells); + cell->ComputeConstantsConsumptionSecretion(); // Update constants for + // all Consumption of + // oxygen + if (cell->IsAttachedToTumorCell()) { // Detach from tumor cell if it + // was attached + cell->GetAttachedCell()->SetAttachedToCart(false); + cell->SetAttachedCell(nullptr); + cell->SetAttachedToTumorCell(false); + } + } else { + // decrease current life time + cell->SetCurrentLiveTime( + (cell->GetCurrentLiveTime() - (kDtCycle * kDtCycle))); + } + break; + } + case CartCellState::kApoptotic: { + cell->SetTimerState(static_cast( + static_cast(cell->GetTimerState()) + kDtCycle)); + + cell->ChangeVolumeExponentialRelaxationEquation( + kVolumeRelaxarionRateCytoplasmApoptotic, + kVolumeRelaxarionRateNucleusApoptotic, + kVolumeRelaxarionRateFluidApoptotic); + + if (kTimeApoptosis < + cell->GetTimerState()) { // If the timer_state exceeds the time to + // transition (this is a fixed duration + // transition) + // remove the cell from the simulation + auto* ctxt = sim->GetExecutionContext(); + ctxt->RemoveAgent(agent->GetUid()); + } + break; + } + default: { + Log::Error("StateControlCart::Run", "Unknown CartCellState"); + break; + } + } + } else { + Log::Error("StateControlCart::Run", "SimObject is not a CartCell"); + } +} + +} // namespace bdm \ No newline at end of file diff --git a/src/cart_cell.h b/src/cart_cell.h new file mode 100644 index 0000000..d370054 --- /dev/null +++ b/src/cart_cell.h @@ -0,0 +1,287 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#ifndef CART_CELL_H_ +#define CART_CELL_H_ + +#include "core/agent/agent.h" +#include "core/agent/cell.h" +#include "core/agent/new_agent_event.h" +#include "core/behavior/behavior.h" +#include "core/container/math_array.h" +#include "core/diffusion/diffusion_grid.h" +#include "core/interaction_force.h" +#include "core/real_t.h" + +#include "tumor_cell.h" + +namespace bdm { + +/// Enumeration defining the possible states of a CAR-T cell +/// +/// This enum class represents the different states that a CAR-T cell can be in +/// during its lifecycle in the simulation. +enum class CartCellState : int { + kAlive = + 0, ///< Living cell state - the cell is alive and functioning normally + + kApoptotic = + 1 ///< Apoptotic phase - the cell is undergoing programmed cell death + ///< characterized by cell shrinkage and controlled death +}; + +/// CAR-T cell class implementation +/// +/// This class represents a CAR-T (Chimeric Antigen Receptor T-cell) in the +/// simulation. It inherits from the base Cell class and includes specific +/// behaviors and properties related to CAR-T cell biology, including states, +/// volume dynamics, and interactions with tumor cells. +class CartCell : public Cell { + BDM_AGENT_HEADER(CartCell, Cell, 1); + + public: + CartCell() = default; + explicit CartCell(const Real3& position); + + // Copy and move constructors/destructors/assignment operators + CartCell(const CartCell&) = default; + CartCell(CartCell&&) = default; + ~CartCell() override = default; + + /// Getters and Setters + void SetState(CartCellState state) { state_ = state; } + CartCellState GetState() const { return state_; } + + void SetTimerState(int timer_state) { timer_state_ = timer_state; } + int GetTimerState() const { return timer_state_; } + + void SetFluidFraction(real_t fluid_fraction) { + fluid_fraction_ = fluid_fraction; + } + real_t GetFluidFraction() const { return fluid_fraction_; } + + void SetNuclearVolume(real_t nuclear_volume) { + nuclear_volume_ = nuclear_volume; + } + real_t GetNuclearVolume() const { return nuclear_volume_; } + + void SetTargetCytoplasmSolid(real_t target_cytoplasm_solid) { + target_cytoplasm_solid_ = target_cytoplasm_solid; + } + real_t GetTargetCytoplasmSolid() const { return target_cytoplasm_solid_; } + + void SetTargetNucleusSolid(real_t target_nucleus_solid) { + target_nucleus_solid_ = target_nucleus_solid; + } + real_t GetTargetNucleusSolid() const { return target_nucleus_solid_; } + + void SetTargetFractionFluid(real_t target_fraction_fluid) { + target_fraction_fluid_ = target_fraction_fluid; + } + real_t GetTargetFractionFluid() const { return target_fraction_fluid_; } + + void SetTargetRelationCytoplasmNucleus( + real_t target_relation_cytoplasm_nucleus) { + target_relation_cytoplasm_nucleus_ = target_relation_cytoplasm_nucleus; + } + real_t GetTargetRelationCytoplasmNucleus() const { + return target_relation_cytoplasm_nucleus_; + } + + void SetAttachedToTumorCell(bool attached) { + attached_to_tumor_cell_ = attached; + } + bool IsAttachedToTumorCell() const { return attached_to_tumor_cell_; } + + Real3 GetOlderVelocity() const { return older_velocity_; } + void SetOlderVelocity(const Real3& velocity) { older_velocity_ = velocity; } + + real_t GetOxygenConsumptionRate() const { return oxygen_consumption_rate_; } + void SetOxygenConsumptionRate(real_t rate) { + oxygen_consumption_rate_ = rate; + } + + 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; } + + /// Returns whether the cell moves by its own + bool DoesCellMove(); + + real_t GetTargetTotalVolume() const; + + /// Returns the diffusion grid for oxygen + DiffusionGrid* GetOxygenDiffusionGrid() const { return oxygen_dgrid_; } + /// Returns the diffusion grid for immunostimulatory factors + DiffusionGrid* GetImmunostimulatoryFactorDiffusionGrid() const { + return immunostimulatory_factor_dgrid_; + } + + /// Change volume using exponential relaxation equation + /// + /// This method explicitly solves the system of exponential relaxation + /// differential equations using a discrete update step. It is used to grow or + /// shrink the volume (and proportions) smoothly toward a desired target + /// volume over time. The relaxation rate controls the speed of convergence + /// and dt=1 (the time_step). + /// + /// @param relaxation_rate_cytoplasm Relaxation rate for cytoplasm volume + /// changes + /// @param relaxation_rate_nucleus Relaxation rate for nucleus volume changes + /// @param relaxation_rate_fluid Relaxation rate for fluid volume changes + void ChangeVolumeExponentialRelaxationEquation( + real_t relaxation_rate_cytoplasm, real_t relaxation_rate_nucleus, + real_t + relaxation_rate_fluid); // NOLINT(bugprone-easily-swappable-parameters) + + /// Calculate displacement of the cell + /// + /// Computes the displacement of the cell based on interaction forces. + /// + /// @param force Pointer to the interaction force object + /// @param squared_radius The squared radius of the cell + /// @param dt The time step for the simulation + /// @return The calculated displacement vector + Real3 CalculateDisplacement(const InteractionForce* force, + real_t squared_radius, real_t dt) override; + + /// Consume or secrete substances + /// + /// Computes new oxygen or immunostimulatory factor concentration after + /// consumption or secretion by the cell. + /// + /// @param substance_id The ID of the substance (oxygen or immunostimulatory + /// factor) + /// @param old_concentration The previous concentration of the substance + /// @return The new concentration after consumption/secretion + real_t ConsumeSecreteSubstance(int substance_id, real_t old_concentration); + + /// Compute constants for consumption and secretion + /// + /// Updates constants after the cell's change of volume or quantities. + /// These constants are used in the consumption/secretion differential + /// equations. + void ComputeConstantsConsumptionSecretion(); + + private: + /// Current state of the CAR-T cell + // NOLINTNEXTLINE(readability-identifier-naming) + CartCellState state_ = CartCellState::kAlive; + + /// Timer to track time in the current state (in minutes) + /// Used for apoptotic state timing + // NOLINTNEXTLINE(readability-identifier-naming) + int timer_state_ = 0; + + /// Pointer to the oxygen diffusion grid + // NOLINTNEXTLINE(readability-identifier-naming) + DiffusionGrid* oxygen_dgrid_ = nullptr; + + // NOLINTNEXTLINE(readability-identifier-naming) + /// Pointer to the immunostimulatory factor diffusion grid + // NOLINTNEXTLINE(readability-identifier-naming) + DiffusionGrid* immunostimulatory_factor_dgrid_ = nullptr; + + /// Flag indicating if the cell is attached to a tumor cell + // NOLINTNEXTLINE(readability-identifier-naming) + bool attached_to_tumor_cell_ = false; + + /// Current time until apoptosis + // NOLINTNEXTLINE(readability-identifier-naming) + real_t current_live_time_ = 0.0; + + /// Fluid fraction of the cell volume + // NOLINTNEXTLINE(readability-identifier-naming) + real_t fluid_fraction_ = 0.0; + + /// Volume of the nucleus + // NOLINTNEXTLINE(readability-identifier-naming) + real_t nuclear_volume_ = 0.0; + + /// Target cytoplasm solid volume for exponential relaxation + /// Used during volume changes following exponential relaxation equation + // NOLINTNEXTLINE(readability-identifier-naming) + real_t target_cytoplasm_solid_ = 0.0; + + /// Target nucleus solid volume for exponential relaxation + // NOLINTNEXTLINE(readability-identifier-naming) + real_t target_nucleus_solid_ = 0.0; + + /// Target fluid fraction for exponential relaxation + // NOLINTNEXTLINE(readability-identifier-naming) + real_t target_fraction_fluid_ = 0.0; + + /// Target relation between cytoplasm and nucleus volumes + // NOLINTNEXTLINE(readability-identifier-naming) + real_t target_relation_cytoplasm_nucleus_ = 0.0; + + /// Velocity of the cell in the previous time step + // NOLINTNEXTLINE(readability-identifier-naming) + Real3 older_velocity_ = {0.0, 0.0, 0.0}; + + /// Rate of oxygen consumption by the cell + // NOLINTNEXTLINE(readability-identifier-naming) + real_t oxygen_consumption_rate_ = 0.0; + + /// Rate of immunostimulatory factor secretion by the cell + // NOLINTNEXTLINE(readability-identifier-naming) + real_t immunostimulatory_factor_secretion_rate_ = 0.0; + + /// Constant 1 for oxygen consumption/secretion differential equation solution + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant1_oxygen_ = 0.0; + + /// Constant 2 for oxygen consumption/secretion differential equation solution + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant2_oxygen_ = 0.0; + + /// Pointer to the attached tumor cell + // NOLINTNEXTLINE(readability-identifier-naming) + TumorCell* attached_cell_ = nullptr; +}; + +/// Behavior class for controlling CAR-T cell state transitions +/// +/// This behavior handles the state control logic for CAR-T cells, managing +/// transitions between different cell states: alive and apoptotic phases. +/// It inherits from the base Behavior class and implements the Run method to +/// execute the state control logic during simulation steps. +struct StateControlCart : public Behavior { + BDM_BEHAVIOR_HEADER(StateControlCart, Behavior, 1); + + StateControlCart() { AlwaysCopyToNew(); } + + // Special member functions + StateControlCart(const StateControlCart&) = default; + StateControlCart(StateControlCart&&) = default; + StateControlCart& operator=(const StateControlCart&) = default; + StateControlCart& operator=(StateControlCart&&) = default; + ~StateControlCart() override = default; + + /// Execute the state control behavior + void Run(Agent* agent) override; +}; + +} // namespace bdm + +#endif // CART_CELL_H_ diff --git a/src/cart_tumor.cc b/src/cart_tumor.cc new file mode 100644 index 0000000..10ac19e --- /dev/null +++ b/src/cart_tumor.cc @@ -0,0 +1,155 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#include +#include +#include + +#include "core/container/math_array.h" +#include "core/diffusion/diffusion_grid.h" +#include "core/environment/uniform_grid_environment.h" +#include "core/model_initializer.h" +#include "core/operation/mechanical_forces_op.h" +#include "core/operation/operation.h" +#include "core/param/param.h" +#include "core/real_t.h" +#include "core/simulation.h" + +#include "cart_tumor.h" +#include "diffusion_thomas_algorithm.h" +#include "forces_tumor_cart.h" +#include "hyperparams.h" +#include "tumor_cell.h" +#include "utils_aux.h" + +namespace bdm { + +int Simulate(int argc, const char** argv) { + // Set simulation bounds + auto set_param = [](Param* param) { + param->random_seed = kSeed; // Set a fixed random seed for reproducibility + param->bound_space = Param::BoundSpaceMode::kTorus; // Periodic boundary + param->min_bound = + -kBoundedSpaceLength / + 2.0; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + param->max_bound = + kBoundedSpaceLength / + 2.0; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + // Cube of 1000x1000x1000 centered at origin + param->simulation_time_step = kDt; + }; + + Simulation simulation(argc, argv, set_param); + auto* ctxt = simulation.GetExecutionContext(); + + // Change Forces + auto* scheduler = simulation.GetScheduler(); + + auto* op = scheduler->GetOps("mechanical forces")[0]; + std::unique_ptr interaction_velocity = + std::make_unique(); + op->GetImplementation()->SetInteractionForce( + interaction_velocity.release()); + + auto* env = dynamic_cast( + Simulation::GetActive()->GetEnvironment()); + // Fix the box length for the uniform grid environment + env->SetBoxLength(gKLengthBoxMechanics); + + // Define Substances + auto* rm = Simulation::GetActive()->GetResourceManager(); + + // Oxygen + // substance_id, name, diffusion_coefficient, decay_constant, resolution, + // time_step + std::unique_ptr oxygen_grid = + std::make_unique( + kOxygen, "oxygen", + kDiffusionCoefficientOxygen, // 100000 micrometers^2/minute + kDecayConstantOxygen, // 0.1 minutes^-1 + kResolutionGridSubstances, kDtSubstances, + true); // true indicates Dirichlet border conditions + rm->AddContinuum(oxygen_grid.release()); + + // Immunostimulatory Factor + // substance_id, name, diffusion_coefficient, decay_constant, resolution + std::unique_ptr immunostimulatory_factor_grid = + std::make_unique( + kImmunostimulatoryFactor, "immunostimulatory_factor", + kDiffusionCoefficientImmunostimulatoryFactor, // 1000 + // micrometers^2/minute + kDecayConstantImmunostimulatoryFactor, // 0.016 minutes^-1 + kResolutionGridSubstances, kDtSubstances, + false); // false indicates Neumann border conditions + rm->AddContinuum(immunostimulatory_factor_grid.release()); + + // Boundary Conditions Dirichlet: simulating absorption or total loss at the + // boundaries of the space. + // Oxygen comming from the borders (capillary vessels) + ModelInitializer::AddBoundaryConditions( + kOxygen, BoundaryConditionType::kDirichlet, + std::make_unique( + kOxygenReferenceLevel)); // kOxygenReferenceLevel mmHg is the + // physiological level of oxygen in tissues, + // o2 saturation is 100% at this level + + // This is useless now but should be added this way in a future version of + // BioDynaMo + ModelInitializer::AddBoundaryConditions( + kImmunostimulatoryFactor, BoundaryConditionType::kNeumann, nullptr); + + // Initialize oxygen voxels + // NOLINTNEXTLINE(bugprone-easily-swappable-parameters) + ModelInitializer::InitializeSubstance(kOxygen, [](real_t /*x*/, real_t /*y*/, + real_t /*z*/) { + return kInitialOxygenLevel; // Set all voxels to kInitialOxygenLevel mmHg + }); + + // One spherical tumor of radius kInitialRadiusTumor in the center of the + // simulation space + const std::vector positions = + CreateSphereOfTumorCells(kInitialRadiusTumor); + for (const auto& pos : positions) { + std::unique_ptr tumor_cell = std::make_unique(pos); + std::unique_ptr state_control = + std::make_unique(); + tumor_cell->AddBehavior(state_control.release()); + ctxt->AddAgent(tumor_cell.release()); + } + + // OutputSummary operation + std::unique_ptr summary_op = + std::make_unique("OutputSummary", kOutputCsvInterval); + std::unique_ptr output_summary = + std::make_unique(); + summary_op->AddOperationImpl(bdm::kCpu, output_summary.release()); + scheduler->ScheduleOp(summary_op.release()); + + // Run simulation + // simulate kTotalMinutesToSimulate minutes including the last minute + scheduler->Simulate(1 + kTotalMinutesToSimulate / kDt); + std::cout << "Simulation completed successfully!" << std::endl; + return 0; +} + +} // namespace bdm + +int main(int argc, const char** argv) { return bdm::Simulate(argc, argv); } \ No newline at end of file diff --git a/src/cart_tumor.h b/src/cart_tumor.h new file mode 100644 index 0000000..a8646fa --- /dev/null +++ b/src/cart_tumor.h @@ -0,0 +1,34 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ +#ifndef CART_TUMOR_H_ +#define CART_TUMOR_H_ + +namespace bdm { + +/// List the diffused substances +enum Substances { kImmunostimulatoryFactor, kOxygen }; + +/// Function declaration for the main simulation +int Simulate(int argc, const char** argv); + +} // namespace bdm + +#endif // CART_TUMOR_H_ \ No newline at end of file diff --git a/src/diffusion_thomas_algorithm.cc b/src/diffusion_thomas_algorithm.cc new file mode 100644 index 0000000..2b16aa7 --- /dev/null +++ b/src/diffusion_thomas_algorithm.cc @@ -0,0 +1,354 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#include +#include +#include +#include + +#include "core/agent/agent.h" +#include "core/diffusion/diffusion_grid.h" +#include "core/real_t.h" +#include "core/resource_manager.h" + +#include "cart_cell.h" +#include "diffusion_thomas_algorithm.h" +#include "hyperparams.h" +#include "tumor_cell.h" + +namespace bdm { + +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +DiffusionThomasAlgorithm::DiffusionThomasAlgorithm( + int substance_id, // NOLINT + std::string substance_name, // NOLINT + real_t dc, real_t mu, // NOLINT + real_t resolution, real_t dt, // NOLINT + bool dirichlet_border) // NOLINT + : DiffusionGrid(substance_id, std::move(substance_name), dc, mu, + static_cast( + resolution)), // Added cast for consistency with parent + resolution_(GetResolution()), + d_space_(static_cast(kBoundedSpaceLength) / + static_cast(resolution_)), + dirichlet_border_(dirichlet_border), + jump_i_(1), + jump_j_(static_cast(resolution_)), + jump_k_(static_cast(resolution_ * resolution_)), + constant1_(dc * dt / (d_space_ * d_space_)), + constant1a_(-constant1_), + constant2_( + mu * dt / + 3.0), // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + constant3_(1.0 + 2 * constant1_ + constant2_), + constant3a_(1.0 + constant1_ + constant2_), + thomas_c_x_(resolution_, constant1a_), + thomas_denom_x_(resolution_, constant3_), + thomas_c_y_(resolution_, constant1a_), + thomas_denom_y_(resolution_, constant3_), + thomas_c_z_(resolution_, constant1a_), + thomas_denom_z_(resolution_, constant3_) { + SetTimeStep(dt); + + // Initialize the denominators and coefficients for the Thomas algorithm + InitializeThomasAlgorithmVectors(thomas_denom_x_, thomas_c_x_); + InitializeThomasAlgorithmVectors(thomas_denom_y_, thomas_c_y_); + InitializeThomasAlgorithmVectors(thomas_denom_z_, thomas_c_z_); +} + +void DiffusionThomasAlgorithm::InitializeThomasAlgorithmVectors( + std::vector& thomas_denom, std::vector& thomas_c) const { + thomas_denom[0] = constant3a_; + thomas_denom[resolution_ - 1] = constant3a_; + if (resolution_ == 1) { + thomas_denom[0] = 1.0 + constant2_; + } + thomas_c[0] /= thomas_denom[0]; + for (size_t i = 1; i < resolution_; ++i) { + thomas_denom[i] += constant1_ * thomas_c[i - 1]; + thomas_c[i] /= thomas_denom[i]; + } +} + +// Apply Dirichlet boundary conditions to the grid +void DiffusionThomasAlgorithm::ApplyDirichletBoundaryConditions() { + const auto* dimensions_ptr = GetDimensionsPtr(); + const real_t origin = dimensions_ptr + [0]; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic) + const real_t simulated_time = GetSimulatedTime(); +#pragma omp parallel + { +// We apply the Dirichlet boundary conditions to the first and last layers in +// each direction For z=0 and z=resolution_-1 +#pragma omp for collapse(2) + for (size_t y = 0; y < resolution_; y++) { + for (size_t x = 0; x < resolution_; x++) { + const real_t real_x = origin + static_cast(x) * d_space_; + const real_t real_y = origin + static_cast(y) * d_space_; + // For z=0 + size_t z = 0; + real_t real_z = origin + static_cast(z) * d_space_; + SetConcentration(static_cast(x), static_cast(y), + static_cast(z), + GetBoundaryCondition()->Evaluate( + real_x, real_y, real_z, simulated_time)); + // For z=resolution_-1 + z = resolution_ - 1; + real_z = origin + static_cast(z) * d_space_; + SetConcentration(static_cast(x), static_cast(y), + static_cast(z), + GetBoundaryCondition()->Evaluate( + real_x, real_y, real_z, simulated_time)); + } + } +// For y=0 and y=resolution_-1 +#pragma omp for collapse(2) + for (size_t z = 0; z < resolution_; z++) { + for (size_t x = 0; x < resolution_; x++) { + const real_t real_x = origin + static_cast(x) * d_space_; + const real_t real_z = origin + static_cast(z) * d_space_; + // For y=0 + size_t y = 0; + real_t real_y = origin + static_cast(y) * d_space_; + SetConcentration(static_cast(x), static_cast(y), + static_cast(z), + GetBoundaryCondition()->Evaluate( + real_x, real_y, real_z, simulated_time)); + // For y=resolution_-1 + y = resolution_ - 1; + real_y = origin + static_cast(y) * d_space_; + SetConcentration(static_cast(x), static_cast(y), + static_cast(z), + GetBoundaryCondition()->Evaluate( + real_x, real_y, real_z, simulated_time)); + } + } +// For x=0 and x=resolution_-1 +#pragma omp for collapse(2) + for (size_t z = 0; z < resolution_; z++) { + for (size_t y = 0; y < resolution_; y++) { + const real_t real_y = origin + static_cast(y) * d_space_; + const real_t real_z = origin + static_cast(z) * d_space_; + // For x=0 + size_t x = 0; + real_t real_x = origin + static_cast(x) * d_space_; + SetConcentration(static_cast(x), static_cast(y), + static_cast(z), + GetBoundaryCondition()->Evaluate( + real_x, real_y, real_z, simulated_time)); + // For x=resolution_-1 + x = resolution_ - 1; + real_x = origin + static_cast(x) * d_space_; + SetConcentration(static_cast(x), static_cast(y), + static_cast(z), + GetBoundaryCondition()->Evaluate( + real_x, real_y, real_z, simulated_time)); + } + } + } +} + +// Sets the concentration at a specific voxel +void DiffusionThomasAlgorithm::SetConcentration(size_t idx, real_t amount) { + const auto* all_concentrations = GetAllConcentrations(); + const real_t current_concentration = all_concentrations[idx]; + ChangeConcentrationBy(idx, amount - current_concentration, + InteractionMode::kAdditive, false); +} + +// Flattens the 3D coordinates (x, y, z) into a 1D index +size_t DiffusionThomasAlgorithm::GetBoxIndex(size_t x, size_t y, + size_t z) const { + return z * resolution_ * resolution_ + y * resolution_ + x; +} + +void DiffusionThomasAlgorithm::Step(real_t /*dt*/) { + // check if diffusion coefficient and decay constant are 0 + // i.e. if we don't need to calculate diffusion update + if (IsFixedSubstance()) { + return; + } + DiffuseChemical(); + + // This should be done considering different border cases instead of using the + // dirichlet_border_ flag. However, there is a bug in BioDynaMo that makes + // bc_type be "Neumann" no matter what. In future versions of BioDynaMo this + // should be fixed +} + +// This method solves the Diffusion Diferential equation using the Alternating +// Direction Implicit approach +void DiffusionThomasAlgorithm::DiffuseChemical() { + ApplyBoundaryConditionsIfNeeded(); + + // Solve for X-direction (direction = 0) + SolveDirectionThomas(0); + ApplyBoundaryConditionsIfNeeded(); + + // Solve for Y-direction (direction = 1) + SolveDirectionThomas(1); + ApplyBoundaryConditionsIfNeeded(); + + // Solve for Z-direction (direction = 2) + SolveDirectionThomas(2); + ApplyBoundaryConditionsIfNeeded(); + + // Change of concentration levels because of agents + ComputeConsumptionsSecretions(); +} + +void DiffusionThomasAlgorithm::ApplyBoundaryConditionsIfNeeded() { + if (dirichlet_border_) { + ApplyDirichletBoundaryConditions(); + } +} + +void DiffusionThomasAlgorithm::SolveDirectionThomas(unsigned int direction) { + const auto& thomas_denom = [this, direction]() -> const std::vector& { + if (direction == 0) { + return thomas_denom_x_; + } + if (direction == 1) { + return thomas_denom_y_; + } + // direction == 2 + return thomas_denom_z_; + }(); + + const auto& thomas_c = [this, direction]() -> const std::vector& { + if (direction == 0) { + return thomas_c_x_; + } + if (direction == 1) { + return thomas_c_y_; + } + // direction == 2 + return thomas_c_z_; + }(); + + const unsigned int jump = [this, direction]() -> unsigned int { + if (direction == 0) { + return static_cast(jump_i_); + } + if (direction == 1) { + return static_cast(jump_j_); + } + // direction == 2 + return static_cast(jump_k_); + }(); + +#pragma omp parallel for collapse(2) + for (unsigned int outer = 0; outer < resolution_; outer++) { + for (unsigned int middle = 0; middle < resolution_; middle++) { + // Forward elimination step + ForwardElimination(direction, outer, middle, thomas_denom, jump); + + // Back substitution step + BackSubstitution(direction, outer, middle, thomas_c, jump); + } + } +} + +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +void DiffusionThomasAlgorithm::ForwardElimination( + unsigned int direction, unsigned int outer, unsigned int middle, + const std::vector& thomas_denom, unsigned int jump) { + // Get initial index based on direction + size_t ind = GetLoopIndex(direction, outer, middle, 0); + const auto* all_concentrations = GetAllConcentrations(); + const real_t initial_concentration = all_concentrations[ind]; + SetConcentration(ind, initial_concentration / thomas_denom[0]); + + // Forward elimination loop + for (unsigned int inner = 1; inner < resolution_; inner++) { + ind = GetLoopIndex(direction, outer, middle, inner); + const real_t current_concentration = all_concentrations[ind]; + const real_t prev_concentration = all_concentrations[ind - jump]; + SetConcentration(ind, + (current_concentration + constant1_ * prev_concentration) / + thomas_denom[inner]); + } +} + +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +void DiffusionThomasAlgorithm::BackSubstitution( + unsigned int direction, unsigned int outer, unsigned int middle, + const std::vector& thomas_c, unsigned int jump) { + const auto* all_concentrations = GetAllConcentrations(); + + // Back substitution loop + for (int inner = static_cast(resolution_) - 2; inner >= 0; inner--) { + const size_t ind = GetLoopIndex(direction, outer, middle, + static_cast(inner)); + const real_t current_concentration = all_concentrations[ind]; + const real_t next_concentration = all_concentrations[ind + jump]; + SetConcentration( + ind, current_concentration - + thomas_c[static_cast(inner)] * next_concentration); + } +} + +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +size_t DiffusionThomasAlgorithm::GetLoopIndex(unsigned int direction, + unsigned int outer, + unsigned int middle, + unsigned int inner) const { + switch (direction) { + case 0: // X-direction: outer=k, middle=j, inner=i + return GetBoxIndex(inner, middle, outer); + case 1: // Y-direction: outer=k, middle=i, inner=j + return GetBoxIndex(middle, inner, outer); + case 2: // Z-direction: outer=j, middle=i, inner=k + return GetBoxIndex(middle, outer, inner); + default: + return 0; + } +} + +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(); + // in a future version of BioDynaMo this should be parallelized getting the + // agents inside each chemical voxel and treating each voxel independently. + rm->ForEachAgent([this](bdm::Agent* agent) { + if (auto* cell = dynamic_cast(agent)) { + // Handle TumorCell agents + const auto& pos = cell->GetPosition(); + const real_t conc = this->GetValue(pos); + const real_t new_conc = + cell->ConsumeSecreteSubstance(GetContinuumId(), conc); + this->ChangeConcentrationBy(pos, new_conc - conc, + InteractionMode::kAdditive, false); + } else if (auto* cell = dynamic_cast(agent)) { + // Handle CartCell agents + const auto& pos = cell->GetPosition(); + const real_t conc = GetValue(pos); + const real_t new_conc = + cell->ConsumeSecreteSubstance(GetContinuumId(), conc); + ChangeConcentrationBy(pos, new_conc - conc, InteractionMode::kAdditive, + false); + } + }); +} + +} // namespace bdm \ No newline at end of file diff --git a/src/diffusion_thomas_algorithm.h b/src/diffusion_thomas_algorithm.h new file mode 100644 index 0000000..94b2119 --- /dev/null +++ b/src/diffusion_thomas_algorithm.h @@ -0,0 +1,266 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#ifndef DIFFUSION_THOMAS_ALGORITHM_H_ +#define DIFFUSION_THOMAS_ALGORITHM_H_ + +#include +#include +#include + +#include "core/diffusion/diffusion_grid.h" +#include "core/real_t.h" +#include "core/util/root.h" + +namespace bdm { + +/// Continuum model for the 3D heat equation with exponential decay +/// +/// Implements the diffusion equation, solved implicitly: ∂t u = ∇D∇u - μu +/// Uses the Thomas algorithm for solving tridiagonal systems efficiently. +class DiffusionThomasAlgorithm : public DiffusionGrid { + public: + DiffusionThomasAlgorithm() + : resolution_(0), + d_space_(0.0), + dirichlet_border_(false), + jump_i_(0), + jump_j_(0), + jump_k_(0), + constant1_(0.0), + constant1a_(0.0), + constant2_(0.0), + constant3_(0.0), + constant3a_(0.0) {} + + DiffusionThomasAlgorithm(int substance_id, std::string substance_name, + real_t dc, real_t mu, real_t resolution, real_t dt, + bool dirichlet_border); + + /// Concentration setters + void SetConcentration(real_t x, real_t y, real_t z, real_t amount) { + SetConcentration(GetBoxIndex(static_cast(x), static_cast(y), + static_cast(z)), + amount); + } + + void SetConcentration(size_t idx, real_t amount); + + /// These methods are overridden but empty because they are not used. + /// This should be fixed in future versions of BioDynaMo. + void DiffuseWithClosedEdge(real_t /*dt*/) override {} + void DiffuseWithOpenEdge(real_t /*dt*/) override {} + void DiffuseWithNeumann(real_t /*dt*/) override {} + void DiffuseWithPeriodic(real_t /*dt*/) override {} + void DiffuseWithDirichlet(real_t /*dt*/) override {} + + /// Perform chemical diffusion using Thomas algorithm + /// + /// Computes the diffusion of the substance using the Thomas algorithm + /// for solving tridiagonal systems efficiently. + void DiffuseChemical(); + + /// Execute one simulation step + /// + /// Main stepping function that performs one time step of the simulation, + /// including diffusion and cellular consumption/secretion. + /// + /// @param dt Time step for the simulation + void Step(real_t dt) override; + + /// Compute cellular consumption and secretion effects + /// + /// Handles secretion or consumption of substances following the differential + /// equation: + /// + /// ∂ρ/∂t = ∇·(D ∇ρ) − λ · ρ + sum_{cells in voxel}((V_k / V_voxel) · [ S_k · + /// ( ρ*_k − ρ ) − (S_k + U_k) · ρ ]) + /// + /// Where: + /// - ρ = concentration of the substance in the microenvironment + /// - S_k = secretion rate of cell k + /// - U_k = uptake (consumption) rate of cell k + /// - ρ*_k = saturation (target) density for secretion + /// - V_k = volume of cell k (approximated to default volume of new tumor + /// cell) + /// - V_voxel = volume of the voxel containing the cell + /// - dt = simulation time step + /// + /// In this class, we only model the secretion and consumption of the + /// substance, not its diffusion, which follows: (ρ − σ)/dt = sum_{cells in + /// voxel}((V_k / V_voxel) · [ S_k · ( ρ*_k − ρ ) − (S_k + U_k) · ρ ]) + /// + /// Where σ is the concentration at the previous time step (may include + /// diffusion term). The solution is: ρⁿ⁺¹ = (ρⁿ + dt · (V_k / V_voxel) · S_k + /// · ρ*_k) / [1 + dt · (V_k / V_voxel) · (S_k + U_k)] + /// + /// Where: + /// - ρⁿ = current concentration + /// - ρⁿ⁺¹ = updated concentration after secretion/uptake + /// + /// This assumes secretion is toward a saturation level, and uptake is + /// proportional to ρ. + /// + /// In a future version, consider using a Behavior associated to each agent + /// but controlling the time in which it is applied so that it is executed + /// always after the diffusion module + void ComputeConsumptionsSecretions(); + + private: + /// Number of voxels in each spatial direction + // NOLINTNEXTLINE(readability-identifier-naming) + size_t resolution_; + + /// Voxel side length in micrometers + // NOLINTNEXTLINE(readability-identifier-naming) + real_t d_space_; + + /// Flag indicating Dirichlet boundary conditions + // NOLINTNEXTLINE(readability-identifier-naming) + bool dirichlet_border_; + + /// Index jump for i-direction (x-axis) + // NOLINTNEXTLINE(readability-identifier-naming) + int jump_i_; + + /// Index jump for j-direction (y-axis) + // NOLINTNEXTLINE(readability-identifier-naming) + int jump_j_; + + /// Index jump for k-direction (z-axis) + // NOLINTNEXTLINE(readability-identifier-naming) + int jump_k_; + + /// First diffusion constant + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant1_; + + /// Alternative first diffusion constant + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant1a_; + + /// Second diffusion constant + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant2_; + + /// Third diffusion constant + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant3_; + + /// Alternative third diffusion constant + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant3a_; + + /// Coefficients for x-direction Thomas algorithm + // NOLINTNEXTLINE(readability-identifier-naming) + std::vector thomas_c_x_; + + /// Denominators for x-direction Thomas algorithm + // NOLINTNEXTLINE(readability-identifier-naming) + std::vector thomas_denom_x_; + + /// Coefficients for y-direction Thomas algorithm + // NOLINTNEXTLINE(readability-identifier-naming) + std::vector thomas_c_y_; + + /// Denominators for y-direction Thomas algorithm + // NOLINTNEXTLINE(readability-identifier-naming) + std::vector thomas_denom_y_; + + /// Coefficients for z-direction Thomas algorithm + // NOLINTNEXTLINE(readability-identifier-naming) + std::vector thomas_c_z_; + + /// Denominators for z-direction Thomas algorithm + // NOLINTNEXTLINE(readability-identifier-naming) + std::vector thomas_denom_z_; + + /// Initialize Thomas algorithm coefficient vectors + /// + /// Sets up the precomputed coefficients for efficient Thomas algorithm + /// execution in the specified direction. + /// + /// @param thomas_denom Reference to denominator vector to initialize + /// @param thomas_c Reference to coefficient vector to initialize + void InitializeThomasAlgorithmVectors(std::vector& thomas_denom, + std::vector& thomas_c) const; + + /// Apply Dirichlet boundary conditions to the diffusion grid + /// + /// Sets the boundary values according to Dirichlet boundary conditions, + /// maintaining constant values at the grid boundaries. + void ApplyDirichletBoundaryConditions(); + + /// Convert 3D coordinates to linear index + /// + /// @param x X-coordinate in voxel space + /// @param y Y-coordinate in voxel space + /// @param z Z-coordinate in voxel space + /// @return Linear index in the flattened 3D array + size_t GetBoxIndex(size_t x, size_t y, size_t z) const; + + /// Apply boundary conditions if Dirichlet boundaries are enabled + void ApplyBoundaryConditionsIfNeeded(); + + /// Solve Thomas algorithm for a specific direction + /// + /// @param direction Direction to solve (0=X, 1=Y, 2=Z) + void SolveDirectionThomas(unsigned int direction); + + /// Perform forward elimination step of Thomas algorithm + /// + /// @param direction Direction being solved (0=X, 1=Y, 2=Z) + /// @param outer Outer loop index + /// @param middle Middle loop index + /// @param thomas_denom Precomputed denominators for this direction + /// @param jump Index jump value for this direction + void ForwardElimination(unsigned int direction, unsigned int outer, + unsigned int middle, + const std::vector& thomas_denom, + unsigned int jump); + + /// Perform back substitution step of Thomas algorithm + /// + /// @param direction Direction being solved (0=X, 1=Y, 2=Z) + /// @param outer Outer loop index + /// @param middle Middle loop index + /// @param thomas_c Precomputed coefficients for this direction + /// @param jump Index jump value for this direction + void BackSubstitution(unsigned int direction, unsigned int outer, + unsigned int middle, + const std::vector& thomas_c, unsigned int jump); + + /// Get the linear index for given direction and loop indices + /// + /// @param direction Direction (0=X, 1=Y, 2=Z) + /// @param outer Outer loop index + /// @param middle Middle loop index + /// @param inner Inner loop index + /// @return Linear index in the flattened array + size_t GetLoopIndex(unsigned int direction, unsigned int outer, + unsigned int middle, unsigned int inner) const; + + BDM_CLASS_DEF_OVERRIDE(DiffusionThomasAlgorithm, 1); +}; + +} // namespace bdm + +#endif // DIFFUSION_THOMAS_ALGORITHM_H_ \ No newline at end of file diff --git a/src/forces_tumor_cart.cc b/src/forces_tumor_cart.cc new file mode 100644 index 0000000..9dc9359 --- /dev/null +++ b/src/forces_tumor_cart.cc @@ -0,0 +1,176 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#include +#include +#include + +#include "core/agent/agent.h" +#include "core/agent/cell.h" +#include "core/container/math_array.h" +#include "core/interaction_force.h" +#include "core/real_t.h" + +#include "forces_tumor_cart.h" +#include "hyperparams.h" +#include "tumor_cell.h" +#include "utils_aux.h" + +namespace bdm { + +Real4 InteractionVelocity::Calculate(const Agent* lhs, const Agent* rhs) const { + const auto* a = dynamic_cast(lhs); + const auto* b = dynamic_cast(rhs); + + // Ignore self-interaction + if (a->GetUid() == b->GetUid()) { + return {0.0, 0.0, 0.0, 0.0}; + } + + Real3 displacement = a->GetPosition() - b->GetPosition(); + + // For periodic boundary conditions, we need to adjust the displacement + displacement[0] = + displacement[0] - + (kBoundedSpaceLength)*round(displacement[0] / (kBoundedSpaceLength)); + displacement[1] = + displacement[1] - + (kBoundedSpaceLength)*round(displacement[1] / (kBoundedSpaceLength)); + displacement[2] = + displacement[2] - + (kBoundedSpaceLength)*round(displacement[2] / (kBoundedSpaceLength)); + + const real_t dist_sq = displacement[0] * displacement[0] + + displacement[1] * displacement[1] + + displacement[2] * displacement[2]; + const real_t distance = std::max(std::sqrt(dist_sq), kEpsilonDistance); + + constexpr real_t kHalf = 2.0; + const real_t radius_a = a->GetDiameter() / kHalf; + const real_t radius_b = b->GetDiameter() / kHalf; + const real_t combined_radius = radius_a + radius_b; + // combined_radius=16.8254;//Debug + // std::cout << "Debug: combined_radius = " << combined_radius << ", distance + // = " << distance << std::endl;// Debug output + real_t temp_r = 0.0; + + const auto* a_tumor = dynamic_cast(a); + const auto* b_tumor = dynamic_cast(b); + + if (distance < combined_radius) { + // 1 - d/combined_radius + temp_r = 1.0 - distance / combined_radius; + // (1 - d/combined_radius)^2 + temp_r *= temp_r; + + real_t repulsion = NAN; + // std::cout << "temp_r = " << temp_r<< std::endl;// Debug output + + if ((a_tumor != nullptr) && (b_tumor != nullptr)) { // two tumor cells + repulsion = kRepulsionTumorTumor; // std::sqrt(kRepulsionTumorTumor * + // kRepulsionTumorTumor); + } else if ((a_tumor == nullptr) && + (b_tumor == nullptr)) { // two CAR-T cells + repulsion = + kRepulsionCartCart; // std::sqrt(kRepulsionCartCart*kRepulsionCartCart); + } else { // one tumor cell and one CAR-T + repulsion = std::sqrt(kRepulsionCartTumor * kRepulsionTumorCart); + } + + temp_r *= repulsion; + } + + // std::cout << "temp_r after repulsion = " << temp_r<< std::endl;// Debug + // output + + // Adhesion + const real_t max_interaction_distance = + kMaxRelativeAdhesionDistance * combined_radius; + // max_interaction_distance=21.0318;//Debug + // std::cout << "max_interaction_distance = " << max_interaction_distance << + // std::endl;// Debug output + + if (distance < max_interaction_distance) { + // 1 - d/S + real_t temp_a = 1.0 - distance / max_interaction_distance; + // (1-d/S)^2 + temp_a *= temp_a; + + // std::cout << "temp_a = " << temp_a << std::endl;// Debug output + + real_t adhesion = NAN; // Initialize to NAN + if ((a_tumor != nullptr) && (b_tumor != nullptr)) { // two tumor cells + adhesion = kAdhesionTumorTumor; + } else if ((a_tumor == nullptr) && + (b_tumor == nullptr)) { // two CAR-T cells + adhesion = kAdhesionCartCart; + } else { // one tumor cell and one CAR-T + adhesion = std::sqrt(kAdhesionCartTumor * kAdhesionTumorCart); + } + + // std::cout << "adhesion = " << adhesion << std::endl;// Debug output + + temp_a *= adhesion; + temp_r -= temp_a; + + // std::cout << "temp_a after adhesion= " << temp_a << std::endl;// Debug + // output + } + + if (std::abs(temp_r) < kEpsilon) { + return {0.0, 0.0, 0.0, 0.0}; + } + const real_t force_magnitude = temp_r / distance; + + // Debug Output volcities + // std::ofstream file("output/intercation_velocities.csv", std::ios::app); + // if (file.is_open()) { + + // real_t total_minutes = + // Simulation::GetActive()->GetScheduler()->GetSimulatedTime(); Real3 + // position = a->GetPosition(); + // // Write data to CSV file + // file << total_minutes << ",position" + // << position[0] << "," + // << position[1] << "," + // << position[2] << ",displacement" + // << displacement[0] << "," + // << displacement[1] << "," + // << displacement[2] << ",distance" + // << distance << ",force_magnitude" + // << force_magnitude << ",temp_r" + // << temp_r << "\n"; + // } + // End Debug Output + + // return{0.,0.,0.,0.};//debug + + return {2 * force_magnitude * displacement[0], + 2 * force_magnitude * displacement[1], + 2 * force_magnitude * displacement[2], + 0.0}; // 4th component is unused +} + +InteractionForce* InteractionVelocity::NewCopy() const { + return std::make_unique().release(); +} + +} // namespace bdm \ No newline at end of file diff --git a/src/forces_tumor_cart.h b/src/forces_tumor_cart.h new file mode 100644 index 0000000..3da9d28 --- /dev/null +++ b/src/forces_tumor_cart.h @@ -0,0 +1,64 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#ifndef FORCES_TUMOR_CART_H_ +#define FORCES_TUMOR_CART_H_ + +#include "core/container/math_array.h" +#include "core/interaction_force.h" + +namespace bdm { + +/// Custom interaction force implementation for velocity-based cell interactions +/// +/// This class implements a specialized interaction force that takes into +/// account the velocity of cells when calculating forces between agents (tumor +/// cells and CAR-T cells). It extends the base InteractionForce class to +/// provide custom force calculations specific to the tumor-CAR-T cell +/// interaction simulation. +class InteractionVelocity : public InteractionForce { + public: + InteractionVelocity() = default; + InteractionVelocity(const InteractionVelocity&) = default; + InteractionVelocity& operator=(const InteractionVelocity&) = default; + InteractionVelocity(InteractionVelocity&&) = default; + InteractionVelocity& operator=(InteractionVelocity&&) = default; + + ~InteractionVelocity() override = default; + + /// Calculate interaction force between two agents + /// + /// Computes the force vector between two agents (cells) based on their + /// positions, properties, and velocities. This method is called by the + /// mechanical forces operation during each simulation step. + /// + /// @param lhs Pointer to the first agent (left-hand side) + /// @param rhs Pointer to the second agent (right-hand side) + /// @return Real4 vector containing the force components (fx, fy, fz, + /// magnitude) + Real4 Calculate(const Agent* lhs, const Agent* rhs) const override; + + [[nodiscard]] InteractionForce* NewCopy() const override; +}; + +} // namespace bdm + +#endif // FORCES_TUMOR_CART_H_ \ No newline at end of file diff --git a/src/hyperparams.h b/src/hyperparams.h new file mode 100644 index 0000000..c7b609e --- /dev/null +++ b/src/hyperparams.h @@ -0,0 +1,245 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#ifndef TUMOR_HYPERPARAMS_H_ +#define TUMOR_HYPERPARAMS_H_ + +#include +#include + +#include "core/real_t.h" +#include "core/util/math.h" + +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 +/// + +/// Rate of secretion of immunostimulatory factor of tumor cells per minute +constexpr real_t kRateSecretionImmunostimulatoryFactor = 10.0; +/// Saturation density of immunostimulatory factor for tumor cells +constexpr real_t kSaturationDensityImmunostimulatoryFactor = 1.0; +/// Mean level of oncoprotein expression in tumor cells +constexpr real_t kOncoproteinMean = 1.0; +/// Standard deviation of oncoprotein expression in tumor cells +constexpr real_t kOncoproteinStandardDeviation = 0.25; +/// Oxygen saturation level in tumor cells for proliferation +constexpr real_t kOxygenSaturationInProliferation = 38.0; +/// Limit of oxygen level for tumor cell proliferation +constexpr real_t kOxygenLimitForProliferation = 10.0; +/// Limit of oxygen to start causing necrosis +constexpr real_t kOxygenLimitForNecrosis = 5.0; +/// Limit of oxygen to maximum necrosis probability +constexpr real_t kOxygenLimitForNecrosisMaximum = 2.5; +/// Time in minutes until a lysed necrotic cell is removed from the simulation +constexpr real_t kTimeLysis = 60 * 24 * 60.; +/// Rate of cell division in min**-1 +constexpr real_t kDivisionRate = 0.02717 / 60.0; +/// Maximum rate per minute of necrosis for tumor cells in case of hypoxia with +/// 0 oxygen +constexpr real_t kMaximumNecrosisRate = 1.0 / (6.0 * 60.0); +/// Default oxygen consumption rate of tumor cell +constexpr real_t kDefaultOxygenConsumption = 10.0; +/// Volume parameters +/// Default total volume of a new tumor cell in μm³ +constexpr real_t kDefaultVolumeNewTumorCell = 2494.0; +/// Default volume of the nucleus of a new tumor cell in μm³ +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; + +/// volume relaxation rate (min^-1) for each state +constexpr real_t kVolumeRelaxarionRateAliveCytoplasm = + 0.13 / 60.; // 0.27/ 60.0; +constexpr real_t kVolumeRelaxarionRateAliveNucleus = 0.22 / 60.; // 0.33/60. +constexpr real_t kVolumeRelaxarionRateAliveFluid = 1.3 / 60.; // 3.0/60. + +constexpr real_t kVolumeRelaxarionRateCytoplasmNecroticSwelling = 0.0032 / 60.0; +constexpr real_t kVolumeRelaxarionRateNucleusNecroticSwelling = 0.013 / 60.; +constexpr real_t kVolumeRelaxarionRateFluidNecroticSwelling = 0.050 / 60.0; + +constexpr real_t kVolumeRelaxarionRateCytoplasmNecroticLysed = 0.0032 / 60.00; +constexpr real_t kVolumeRelaxarionRateNucleusNecroticLysed = 0.013 / 60.; +constexpr real_t kVolumeRelaxarionRateFluidNecroticLysed = 0.050 / 60.0; + +/// Thresholds in oncoprotein levels for differentiating 4 cancer cell types +constexpr real_t kThresholdCancerCellType1 = 1.5; +constexpr real_t kThresholdCancerCellType2 = 1.0; +constexpr real_t kThresholdCancerCellType3 = 0.5; +constexpr real_t kThresholdCancerCellType4 = 0.0; + +/// +/// General Hyperparameters +/// +/// Seed for random number generation +constexpr int kSeed = 3; + +/// 0.01 minutes time step for substances secretion/consumption +constexpr real_t kDtSubstances = 0.01; +/// 0.1 minutes time step for the cell mechanics +constexpr real_t kDtMechanics = 0.1; +/// 6 minutes time step for the cell cycle +constexpr real_t kDtCycle = 6.0; + +/// General time step for the simulation: it is the same as kDtMechanics, do not +/// modify this line +constexpr real_t kDt = kDtMechanics; +/// Number of steps per cycle step, do not modify this line. Needs to be +/// computed to avoid errors with fmod +constexpr int kStepsPerCycle = kDtCycle / kDt; + +/// 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; + +constexpr real_t kVolumeRelaxarionRateCytoplasmApoptotic = 1.0 / 60.0; +constexpr real_t kVolumeRelaxarionRateNucleusApoptotic = 0.35 / 60.0; +constexpr real_t kVolumeRelaxarionRateFluidApoptotic = 0.0; +/// Time in minutes until an apoptotic cell is removed from the simulation +constexpr real_t kTimeApoptosis = 8.6 * 60; +/// Reduction of consumption rate of dead cells when they enter necrosis +constexpr real_t kReductionConsumptionDeadCells = 0.1; + +/// Chemicals +/// Number of voxels per axis for the substances grid +constexpr int kResolutionGridSubstances = 50; // 50 // voxels per axis +/// Volume of a single voxel in μm³ (do not modify this line) +constexpr real_t kVoxelVolume = + (static_cast(kBoundedSpaceLength) / kResolutionGridSubstances) * + (static_cast(kBoundedSpaceLength) / kResolutionGridSubstances) * + (static_cast(kBoundedSpaceLength) / + kResolutionGridSubstances); // Do not modify this line +/// Diffusion coefficient of oxygen in μm²/min +constexpr real_t kDiffusionCoefficientOxygen = + 100000; // 100000 micrometers^2/minute +/// Decay constant of oxygen in min⁻¹ +constexpr real_t kDecayConstantOxygen = 0.1; // 0.1 minutes^-1 +/// Diffusion coefficient of immunostimulatory factor in μm²/min +constexpr real_t kDiffusionCoefficientImmunostimulatoryFactor = + 1000; // 1000 micrometers^2/minute +/// Decay constant of immunostimulatory factor in min⁻¹ +constexpr real_t kDecayConstantImmunostimulatoryFactor = + 0.016; // 0.016 minutes^-1 +/// Time step for oxygen diffusion in minutes +constexpr real_t kTimeStepOxygen = 0.0005; // 0.001 minutes CHANGE +/// Time step for immunostimulatory factor diffusion in minutes +constexpr real_t kTimeStepImmunostimulatoryFactor = 0.01; // 0.01 minutes +/// Reference level of oxygen at the boundaries in mmHg +constexpr real_t kOxygenReferenceLevel = + 38.; // Reference level of oxygen at the boundaries of the simulation space + // in mmHg +/// Initial oxygen concentration in each voxel in mmHg +constexpr real_t kInitialOxygenLevel = + 38.0; // Initial voxel concentration of oxygen in mmHg +/// Oxygen saturation in the microenvironment in mmHg +constexpr real_t kOxygenSaturation = + 30.0; // 30.0 // Oxygen saturation in mmHg in microenvironment +/// Forces +/// Repulsion coeficient between tumor cells +constexpr real_t kRepulsionTumorTumor = 10.0; +/// Repulsion coeficient between CAR-T cells +constexpr real_t kRepulsionCartCart = 50.0; +/// Repulsion coeficient between CAR-T cells and tumor cells +constexpr real_t kRepulsionCartTumor = 50.0; +/// Repulsion coeficient between tumor cells and CAR-T cells +constexpr real_t kRepulsionTumorCart = 10.0; +/// Maximum relative adhesion distance for cell interactions +constexpr real_t kMaxRelativeAdhesionDistance = 1.25; +/// Adhesion coeficient between tumor cells +constexpr real_t kAdhesionTumorTumor = 0.4; +/// Adhesion coeficient between CAR-T cells +constexpr real_t kAdhesionCartCart = 0.0; +/// Adhesion coeficient between CAR-T cells and tumor cells +constexpr real_t kAdhesionCartTumor = 0.0; +/// Adhesion coeficient between tumor cells and CAR-T cells +constexpr real_t kAdhesionTumorCart = 0.0; + +/// Do not change +// coefficientes for the 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) ] +/// Coefficient for the current time step in the Adams-Bashforth method (dt +/// * 1.5) +constexpr real_t kDnew = 1.5 * kDtMechanics; +/// Coefficient for the previous time step in the Adams-Bashforth method (dt * +/// -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; // Length of the box for mechanics in micrometers + +/// Max Distance for considering two cells as neighbours for force calculations +/// in μm Do not change this line +const real_t gKSquaredMaxDistanceNeighborsForce = std::pow( + 0.1 + std::cbrt(kDefaultVolumeNewTumorCell * 6 / Math::kPi) * + kMaxRelativeAdhesionDistance, + 2); // (twice biggest cell radius (in case to cells tha maximum size + // encounter each other) times kMaxRelativeAdhesionDistance + 0.1 to + // avoid mismatch because of numerical errors)**2 + +/// +/// CAR-T Cell Hyperparameters +/// +constexpr real_t kAverageMaximumTimeUntillApoptosisCart = + kDtCycle * 10.0 * 24.0 * 60.0; +/// Volume parameters +/// Default total volume of a new CAR-T cell in μm³ +constexpr real_t kDefaultVolumeNewCartCell = 2494.0; +/// Default volume of the nucleus of a new CAR-T cell in μm³ +constexpr real_t kDefaultVolumeNucleusCartCell = 540.0; +/// Default fraction of fluid volume in a new CAR-T cell +constexpr real_t kDefaultFractionFluidCartCell = 0.75; + +} // namespace bdm + +#endif // TUMOR_HYPERPARAMS_H_ diff --git a/src/tumor_cell.cc b/src/tumor_cell.cc new file mode 100644 index 0000000..4c4f718 --- /dev/null +++ b/src/tumor_cell.cc @@ -0,0 +1,656 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#include +#include +#include +#include + +#include "core/agent/agent.h" +#include "core/agent/cell_division_event.h" +#include "core/agent/new_agent_event.h" +#include "core/container/math_array.h" +#include "core/functor.h" +#include "core/interaction_force.h" +#include "core/real_t.h" +#include "core/resource_manager.h" +#include "core/util/log.h" + +#include "hyperparams.h" +#include "tumor_cell.h" +#include "utils_aux.h" + +namespace bdm { + +TumorCell::TumorCell(const Real3& position) { + SetPosition(position); + + // volumes + SetVolume(kDefaultVolumeNewTumorCell); // Set default volume + SetFluidFraction( + kDefaultFractionFluidTumorCell); // Set default fluid fraction + SetNuclearVolume( + kDefaultVolumeNucleusTumorCell); // Set default nuclear volume + // target volumes + SetTargetFractionFluid( + kDefaultFractionFluidTumorCell); // Set target fraction of fluid + SetTargetRelationCytoplasmNucleus( + (kDefaultVolumeNewTumorCell - kDefaultVolumeNucleusTumorCell) / + (kEpsilon + + kDefaultVolumeNucleusTumorCell)); // Set target relation between + // cytoplasm and nucleus + SetTargetNucleusSolid( + kDefaultVolumeNucleusTumorCell * + (1 - kDefaultFractionFluidTumorCell)); // Set target nucleus solid volume + // to real_t + SetTargetCytoplasmSolid( + (kDefaultVolumeNewTumorCell - kDefaultVolumeNucleusTumorCell) * + (1 - kDefaultFractionFluidTumorCell)); // Set target cytoplasm solid + // volume to real_t + + SetOncoproteineLevel(SamplePositiveGaussian( + kOncoproteinMean, + kOncoproteinStandardDeviation)); // Set initial oncoproteine level with a + // truncated normal distribution + // SetOncoproteineLevel(1.); //Debug + oxygen_dgrid_ = + Simulation::GetActive()->GetResourceManager()->GetDiffusionGrid( + "oxygen"); // Pointer to oxygen diffusion grid + immunostimulatory_factor_dgrid_ = + Simulation::GetActive()->GetResourceManager()->GetDiffusionGrid( + "immunostimulatory_factor"); // Pointer to immunostimulatory_factor + // diffusion grid + SetTransformationRandomRate(); // Set state transition random rate + + // Add Consumption and Secretion + SetOxygenConsumptionRate( + kDefaultOxygenConsumption); // Set default oxygen consumption rate + SetImmunostimulatoryFactorSecretionRate( + kRateSecretionImmunostimulatoryFactor); // Set default immunostimulatory + // factor secretion rate + ComputeConstantsConsumptionSecretion(); // Compute constants for all + // ConsumptionSecretion of Oxygen and + // Immunostimulatory Factors +} + +/// Called when a new agent is created (e.g., after cell division) +void TumorCell::Initialize(const NewAgentEvent& event) { + Base::Initialize(event); + + if (auto* mother = dynamic_cast( + event.existing_agent)) { // if the cell is created from division + if (event.GetUid() == CellDivisionEvent::kUid) { + // Initialize daughter cell from mother cell + state_ = TumorCellState::kAlive; // state after division + timer_state_ = 0; + // diffusion grids + oxygen_dgrid_ = + mother->oxygen_dgrid_; // Pointer to the oxygen diffusion grid + immunostimulatory_factor_dgrid_ = + mother->immunostimulatory_factor_dgrid_; // Pointer to the + // immunostimulatory_factor + // diffusion grid + this->SetOncoproteineLevel( + mother->oncoproteine_level_); // inherit oncoproteine level from + // mother cell + this->SetOxygenConsumptionRate( + mother->GetOxygenConsumptionRate()); // inherit oxygen consumption + // rate from mother cell + this->SetImmunostimulatoryFactorSecretionRate( + mother + ->GetImmunostimulatoryFactorSecretionRate()); // inherit + // immunostimulatory + // factor secretion + // rate from mother + // cell + + // Update the constants for all ConsumptionSecretion + mother->ComputeConstantsConsumptionSecretion(); + this->ComputeConstantsConsumptionSecretion(); + + // divde mother's nuclear volume by 2 + const real_t new_nuclear_volume = + mother->GetNuclearVolume() / + 2.0; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + mother->SetNuclearVolume( + new_nuclear_volume); // Set mother's nuclear volume to the new volume + this->SetNuclearVolume(new_nuclear_volume); + + // Inherit mother's fluid fraction and velocity + this->SetFluidFraction( + mother->GetFluidFraction()); // Set fluid fraction to mother's fluid + // fraction + this->SetOlderVelocity( + mother->GetOlderVelocity()); // Copy velocity from mother cell + + // inherit target volumes of the daughter cell + this->SetTargetFractionFluid(mother->GetTargetFractionFluid()); + this->SetTargetRelationCytoplasmNucleus( + mother->GetTargetRelationCytoplasmNucleus()); + this->SetTargetNucleusSolid(mother->GetTargetNucleusSolid()); + this->SetTargetCytoplasmSolid(mother->GetTargetCytoplasmSolid()); + + this->SetTransformationRandomRate(); // Set state transition random rate + this->attached_to_cart_ = false; // Initially not attached to a cart + } + } +} + +void TumorCell::SetOncoproteineLevel(real_t level) { + oncoproteine_level_ = level; // oncoproteine_level_ + // cell type + if (level >= kThresholdCancerCellType1) { // between 1.5 and 2.0 + type_ = TumorCellType::kType1; + } else if (level >= kThresholdCancerCellType2 && + level < kThresholdCancerCellType1) { + type_ = TumorCellType::kType2; + } else if (level >= kThresholdCancerCellType3 && + level < kThresholdCancerCellType2) { + type_ = TumorCellType::kType3; + } else if (level >= kThresholdCancerCellType4 && + level < kThresholdCancerCellType3) { + type_ = TumorCellType::kType4; + } else { // undefined type + type_ = TumorCellType::kType0; + } +} + +void TumorCell::SetTransformationRandomRate() { + // avoid division by zero + transformation_random_rate_ = + 1 / (std::max(SamplePositiveGaussian( + kAverageTimeTransformationRandomRate, + kStandardDeviationTransformationRandomRate) * + kMinutesInAnHour, + kEpsilon)); +} + +real_t TumorCell::GetTargetTotalVolume() const { + return GetTargetNucleusSolid() * (1 + GetTargetRelationCytoplasmNucleus()) / + (1 - GetTargetFractionFluid()); +} + +// This method explicitly solves the system of exponential relaxation +// differential equation using a discrete update step. It is used to grow or +// shrink the volume (and proportions) smoothly toward a desired target volume +// over time. The relaxations rate controls the speed of convergence +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +void TumorCell::ChangeVolumeExponentialRelaxationEquation( + real_t relaxation_rate_cytoplasm, real_t relaxation_rate_nucleus, // NOLINT + real_t relaxation_rate_fluid) { + // Exponential relaxation towards the target volume + const real_t current_total_volume = GetVolume(); + const real_t fluid_fraction = GetFluidFraction(); + const real_t nuclear_volume = GetNuclearVolume(); + + const real_t current_nuclear_solid = nuclear_volume * (1 - fluid_fraction); + const real_t current_cytoplasm_solid = + (current_total_volume - nuclear_volume) * (1 - fluid_fraction); + + // std::cout << "time=" << + // Simulation::GetActive()->GetScheduler()->GetSimulatedTime() << + // ", current_total_volume=" << current_total_volume << + // ", current_nuclear_volume=" << nuclear_volume << + // ", current_cytoplasm_solid=" << current_cytoplasm_solid << + // std::endl; + + const real_t current_fluid = fluid_fraction * current_total_volume; + + // Update fluid volume + real_t new_fluid = + current_fluid + + kDtCycle * relaxation_rate_fluid * + (GetTargetFractionFluid() * current_total_volume - current_fluid); + // Clamp to zero to prevent negative volumes + if (new_fluid < 0.0) { + new_fluid = 0.0; + } + + const real_t nuclear_fluid = + new_fluid * (nuclear_volume / current_total_volume); + // real_t cytoplasm_fluid = new_fluid - nuclear_fluid; + + real_t nuclear_solid = current_nuclear_solid + + kDtCycle * relaxation_rate_nucleus * + (GetTargetNucleusSolid() - current_nuclear_solid); + // Clamp to zero to prevent negative volumes + if (nuclear_solid < 0.0) { + nuclear_solid = 0.0; + } + + const real_t target_cytoplasm_solid = + GetTargetRelationCytoplasmNucleus() * GetTargetNucleusSolid(); + real_t cytoplasm_solid = + current_cytoplasm_solid + + kDtCycle * relaxation_rate_cytoplasm * + (target_cytoplasm_solid - current_cytoplasm_solid); + // Clamp to zero to prevent negative volumes + if (cytoplasm_solid < 0.0) { + cytoplasm_solid = 0.0; + } + + const real_t new_total_solid = nuclear_solid + cytoplasm_solid; + + const real_t total_nuclear = nuclear_solid + nuclear_fluid; + + // real_t total_cytoplasm= cytoplasm_solid + cytoplasm_fluid; + + const real_t new_volume = new_total_solid + new_fluid; + + // Avoid division by zero + const real_t new_fraction_fluid = new_fluid / (kEpsilon + new_volume); + + // Debug Debug Output params + // std::ofstream file("output/volumes.csv", std::ios::app); + // if (file.is_open()) { + + // // Write data to CSV file + // file << Simulation::GetActive()->GetScheduler()->GetSimulatedTime() << + // ",cytoplasm," + // << new_volume-total_nuclear << ",nuclear," + // << total_nuclear <<",fraction fluid," + // << new_fraction_fluid<< "cytoplasm solid: " + // <GetExecutionContext(); + auto calculate_neighbor_forces = + 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]; + } + }); + ctxt->ForEachNeighbor(calculate_neighbor_forces, *this, squared_radius); + + if (non_zero_neighbor_forces > 1) { + SetStaticnessNextTimestep(false); + } + } + + // 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; + + older_velocity_ = translation_velocity_on_point_mass; + + // Displacement + return movement_at_next_step; +} + +// Compute new oxygen or immunostimulatory factor concentration after +// consumption/ secretion +// NOLINTNEXTLINE(bugprone-easily-swappable-parameters) +real_t TumorCell::ConsumeSecreteSubstance(int substance_id, + real_t old_concentration) { + // constant1_oxygen_ = 0; // Debug + // constant2_oxygen_ = 1.3; // Debug + real_t res = 0.0; + if (substance_id == oxygen_dgrid_->GetContinuumId()) { + // consuming oxygen + res = (old_concentration + constant1_oxygen_) / constant2_oxygen_; + } else if (substance_id == + immunostimulatory_factor_dgrid_->GetContinuumId()) { + // secreting immunostimulatory factor + res = (old_concentration + constant1_immunostimulatory_factor_) / + constant2_immunostimulatory_factor_; + } else { + throw std::invalid_argument("Unknown substance id: " + + std::to_string(substance_id)); + } + return res; +} + +void TumorCell::ComputeConstantsConsumptionSecretion() { + // constant1_= dt · (V_k / V_voxel) · S_k · ρ*_k) + // constant2_ = [1 + dt · (V_k / V_voxel) · (S_k + U_k)] + // where: + // S_k = secretion rate of cell k + // U_k = uptake (consumption) rate of cell k + // ρ*_k = saturation (target) density for secretion + // V_k = volume of the cell k + // V_voxel = volume of the voxel containing the cell + // dt = simulation time step + + const real_t new_volume = GetVolume(); + // compute the constants for the differential equation explicit solution: for + // oxygen and immunostimulatory factor + // dt*(cell_volume/voxel_volume)*quantity_secretion*substance_saturation = dt + // · (V_k / V_voxel) · S_k · ρ*_k) + constant1_oxygen_ = 0.; + // Scale by the volume of the cell in the Voxel and time step + constant1_immunostimulatory_factor_ = + immunostimulatory_factor_secretion_rate_ * + kSaturationDensityImmunostimulatoryFactor * kDtSubstances * + (new_volume / kVoxelVolume); + // 1 + dt*(cell_volume/voxel_volume)*(quantity_secretion + + // quantity_consumption ) = [1 + dt · (V_k / V_voxel) · (S_k + U_k)] + // Scale by the volume of the cell in the Voxel and time step + constant2_oxygen_ = 1 + kDtSubstances * (new_volume / kVoxelVolume) * + (oxygen_consumption_rate_); + // Scale by the volume of the cell in the Voxel and time step + constant2_immunostimulatory_factor_ = + 1 + kDtSubstances * (new_volume / kVoxelVolume) * + (immunostimulatory_factor_secretion_rate_); +} + +/// Main behavior executed at each simulation step +void StateControlGrowProliferate::Run(Agent* agent) { + auto* sim = Simulation::GetActive(); + if (sim->GetScheduler()->GetSimulatedSteps() % kStepsPerCycle != 0) { + return; + } // Run only every kDtCycle minutes, fmod does not work with the type + // returned by GetSimulatedTime() + // Debug + // // Print simulation minute and number of TumorCell agents + // int num_steps = sim->GetScheduler()->GetSimulatedSteps(); + // int current_minute = 6 * num_steps; + // size_t num_cells = sim->GetResourceManager()->GetNumAgents(); + // int current_hour = current_minute / 60; + // int current_day = current_hour / 24; + // std::cout << "Dia: " << current_day << " Hora: " << (current_hour % 24) + // << " Minuto: " << (current_minute % 60) + // << " Numero de celulas: " << num_cells << std::endl; + // // ---------------------------------------- + // // End Debug + + if (auto* cell = dynamic_cast(agent)) { + if (cell->IsAttachedToCart()) { + // If the cell is attached to a cart, skip the state control and growth + return; + } + // Oxygen levels + const Real3 current_position = cell->GetPosition(); + auto* oxygen_dgrid = + cell->GetOxygenDiffusionGrid(); // Pointer to the oxygen diffusion grid + const real_t oxygen_level = oxygen_dgrid->GetValue(current_position); + // oxygen_level = 30.; // Debug + + // Debug + // std::cout << oxygen_level << std::endl; + + switch (cell->GetState()) { + case TumorCellState::kAlive: { // the cell is growing to real_t its size + // before mitosis + cell->SetTimerState( + cell->GetTimerState() + + kDtCycle); // Increase timer_state to track time in this state + // (kDtCycle minutes per step) + + if (ShouldEnterNecrosis( + oxygen_level, + cell)) { // Enter necrosis if oxygen level is too low + return; // Exit the function to prevent further processing + } + ManageLivingCell(cell, oxygen_level); + break; + } + case TumorCellState::kNecroticSwelling: { // the cell is swelling before + // lysing + cell->SetTimerState( + cell->GetTimerState() + + kDtCycle); // Increase timer_state to track time in this state + // (kDtCycle minutes per step) + // volume change + // The cell swells + cell->ChangeVolumeExponentialRelaxationEquation( + kVolumeRelaxarionRateCytoplasmNecroticSwelling, + kVolumeRelaxarionRateNucleusNecroticSwelling, + kVolumeRelaxarionRateFluidNecroticSwelling); + if (cell->GetVolume() >= + 2 * kDefaultVolumeNewTumorCell) { // If the cell has swollen to 2 + // times its original volume, it + // lyses + cell->SetState(TumorCellState::kNecroticLysed); // Change state to + // necrotic lysed + cell->SetTimerState(0); // Reset timer_state + // Set target volume to 0 (the cell will shrink) + cell->SetTargetCytoplasmSolid(0.0); + cell->SetTargetNucleusSolid(0.0); + cell->SetTargetFractionFluid(0.0); + cell->SetTargetRelationCytoplasmNucleus(0.0); + // Stop secretion and consumption rate + // Stop consumption + cell->SetOxygenConsumptionRate(0.0); + // Stop secretion + cell->SetImmunostimulatoryFactorSecretionRate(0.0); + // Update constants for all ConsumptionSecretion of Oxygen and + // Immunostimulatory Factors + cell->ComputeConstantsConsumptionSecretion(); + } + break; + } + case TumorCellState::kNecroticLysed: { // the cell is shirinking and will + // be removed after a certain time + cell->SetTimerState( + cell->GetTimerState() + + kDtCycle); // Increase timer_state to track time in this state + // (kDtCycle minutes per step) + // volume change + // The cell shrinks + cell->ChangeVolumeExponentialRelaxationEquation( + kVolumeRelaxarionRateCytoplasmNecroticLysed, + kVolumeRelaxarionRateNucleusNecroticLysed, + kVolumeRelaxarionRateFluidNecroticLysed); + if (kTimeLysis < + cell->GetTimerState()) { // If the timer_state exceeds the time to + // transition (this is a fixed duration + // transition) + // remove the cell from the simulation + auto* ctxt = sim->GetExecutionContext(); + ctxt->RemoveAgent(agent->GetUid()); + } + break; + } + case TumorCellState::kApoptotic: { + // CHANGE write this in the function that causes apoptosis + // //Stop Secretion and reduce consumption + // for (auto* beh : cell->GetAllBehaviors()) { + // if (auto* c = dynamic_cast(beh)) + // {c->SetQuantity(c->GetQuantity()*kReductionConsumptionDeadCells);}// + // Reduce consumption rate else if (auto* s = + // dynamic_cast(beh)) + // {s->SetQuantity(0);} // Stop secretion of immunostimulatory + // factor + // } + // cell->SetType(5); // Set type to 5 to indicate dead cell + + cell->SetTimerState( + cell->GetTimerState() + + kDtCycle); // Increase timer_state to track time in this state + // (kDtCycle minutes per step) + // volume change CHANGe check if it should indeed be reduced to 0 + // The cell shrinks + cell->ChangeVolumeExponentialRelaxationEquation( + kVolumeRelaxarionRateCytoplasmApoptotic, + kVolumeRelaxarionRateNucleusApoptotic, + kVolumeRelaxarionRateFluidApoptotic); + if (kTimeApoptosis < + cell->GetTimerState()) { // If the timer_state exceeds the time to + // transition (this is a fixed duration + // transition) + // remove the cell from the simulation + auto* ctxt = sim->GetExecutionContext(); + ctxt->RemoveAgent(agent->GetUid()); + } + break; + } + default: { + Log::Error("StateControlGrowProliferate::Run", + "Unknown TumorCellState"); + break; + } + } + } else { + Log::Error("StateControlGrowProliferate::Run", + "SimObject is not a TumorCell"); + } +} + +// ManageLivingCell function to handle living cell behavior +void StateControlGrowProliferate::ManageLivingCell(TumorCell* cell, + real_t oxygen_level) { + real_t multiplier = 1.0; // Initialize multiplier + // volume change + cell->ChangeVolumeExponentialRelaxationEquation( + kVolumeRelaxarionRateAliveCytoplasm, kVolumeRelaxarionRateAliveNucleus, + kVolumeRelaxarionRateAliveFluid); // The cell grows to real_t its + // size + // cell state control + if (oxygen_level < + kOxygenSaturationInProliferation) { // oxygen threshold for + // considering an effect on the + // proliferation cycle + multiplier = + (oxygen_level - kOxygenLimitForProliferation) / + (kOxygenSaturationInProliferation - kOxygenLimitForProliferation); + } + if (oxygen_level < kOxygenLimitForProliferation) { + multiplier = 0.0; // If oxygen is below the limit, set multiplier to 0 + } + // double multiplier1 = multiplier; //Debug + + const real_t final_rate_transition = + cell->GetTransformationRandomRate() * multiplier * + cell->GetOncoproteineLevel(); // Calculate the rate of state change + // based on oxygen level and + // oncoproteine (min^-1) + + // Debug + // int current_time = sim->GetScheduler()->GetSimulatedSteps()* kDt; // + // Get the current time step in minutes std::ofstream + // file("output/simulation_data_mine" + + // std::to_string(current_time/(12*60)) + ".csv", std::ios::app); if + // (file.is_open()) { file << oxygen_level << "," + // << cell->GetOncoproteineLevel() << "," + // <GetTransformationRandomRate()<< "," + // << final_rate_transition << "\n"; + // } + // End Debug + + real_t time_to_wait = + kTimeTooLarge; // Set a very large time to avoid division by zero + if (final_rate_transition > 0) { + time_to_wait = 1. / final_rate_transition; // Calculate the time to + // transition (in minutes ) + } + if (time_to_wait < + cell->GetTimerState()) { // If the timer_state exceeds the time to + // transition, change state (this is a + // fixed duration transition) + // mitosis: cell divides + cell->SetState(TumorCellState::kAlive); + cell->Divide(); + cell->SetTimerState(0); // Reset timer_state + } +} + +// computes the probability of the cell entering necrosis +bool StateControlGrowProliferate::ShouldEnterNecrosis(real_t oxygen_level, + TumorCell* cell) { + // necrosis probability + real_t multiplier = 0.0; // Default multiplier for necrosis probability + + if (oxygen_level < + kOxygenLimitForNecrosis) { // oxygen threshold for considering necrosis + multiplier = (kOxygenLimitForNecrosis - oxygen_level) / + (kOxygenLimitForNecrosis - kOxygenLimitForNecrosisMaximum); + } + if (oxygen_level < kOxygenLimitForNecrosisMaximum) { // threshold for maximum + // necrosis probability + multiplier = 1.0; + } + + const real_t probability_necrosis = + kDtCycle // multiply by kDtCycle since each timestep is kDtCycle minutes + * kMaximumNecrosisRate * multiplier; // Calculate the probability of + // necrosis based on oxygen level + + auto* sim = Simulation::GetActive(); + auto* random = sim->GetRandom(); + const bool enter_necrosis = random->Uniform(0, 1) < probability_necrosis; + if (enter_necrosis) { // If the random number is less than the probability, + // enter necrosis + cell->SetState(TumorCellState::kNecroticSwelling); // If oxygen is too low, + // enter necrosis + cell->SetTimerState(0); // Reset timer_state + + // Stop Secretion and reduce consumption + // Stop secretion + cell->SetImmunostimulatoryFactorSecretionRate(0.0); + // Reduce consumption + cell->SetOxygenConsumptionRate(cell->GetOxygenConsumptionRate() * + kReductionConsumptionDeadCells); + // Update constants for all ConsumptionSecretion of Oxygen and + // Immunostimulatory Factors + cell->ComputeConstantsConsumptionSecretion(); + + // The cell will swell getting filled with fluid + cell->SetTargetCytoplasmSolid(0); + cell->SetTargetNucleusSolid(0); + cell->SetTargetFractionFluid(1.0); // Set target fraction of fluid to 1.0 + cell->SetTargetRelationCytoplasmNucleus(0.0); + // Set type to 5 to indicate dead cell + cell->SetType(TumorCellType::kType5); + } + return enter_necrosis; // Return whether the cell entered necrosis +} + +} // namespace bdm \ No newline at end of file diff --git a/src/tumor_cell.h b/src/tumor_cell.h new file mode 100644 index 0000000..82edde3 --- /dev/null +++ b/src/tumor_cell.h @@ -0,0 +1,368 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#ifndef TUMOR_CELL_H_ +#define TUMOR_CELL_H_ + +#include "core/agent/agent.h" +#include "core/agent/cell.h" +#include "core/agent/new_agent_event.h" +#include "core/behavior/behavior.h" +#include "core/container/math_array.h" +#include "core/diffusion/diffusion_grid.h" +#include "core/interaction_force.h" +#include "core/real_t.h" +#include "core/scheduler.h" + +namespace bdm { + +/// Enumeration representing the different states of a tumor cell +/// +/// This enum class defines the various states a tumor cell can be in during its +/// lifecycle, and various death pathways (necrosis and apoptosis). +enum class TumorCellState : int { + kAlive = + 0, ///< Living cell state - cell is alive and can potentially proliferate + + kNecroticSwelling = + 1, ///< Necrotic swelling phase: The cell loses membrane integrity and + ///< starts absorbing fluid, swelling abnormally, in volume before + ///< rupture. This is the first phase of necrotic cell death. + kNecroticLysed = 2, ///< Necrotic lysed phase: The cell membrane breaks + ///< apart, releasing its contents. The cell will be + ///< removed from the simulation after a defined time. + + kApoptotic = 3 ///< Apoptotic phase: The cell is undergoing programmed cell + ///< death characterized by cell shrinkage. This is a + ///< controlled form of cell death. +}; + +/// Enumeration representing the different types of tumor cells +/// +/// This enum class defines the various types of tumor cells depending on their +/// agressiveness from 1 (most agressive) until 4 (least agressive) based on the +/// expressed oncoprotein level. Type 5 cells are dead cells +enum class TumorCellType : int { + kType0 = 0, ///< Unclassified tumor cells + kType1 = 1, ///< Most aggressive tumor cells + kType2 = 2, ///< Moderately aggressive tumor cells + kType3 = 3, ///< Less aggressive tumor cells + kType4 = 4, ///< Least aggressive tumor cells + kType5 = 5 ///< Dead tumor cells +}; + +/// Tumor cell class implementation +/// +/// This class represents a cancer cell that forms a heterogeneous tumor in the +/// simulation. The class includes capabilities for: +/// - Different cellular states (alive, necrotic, apoptotic) +/// - Volume dynamics with exponential relaxation +/// - Cell division for tumor proliferation +/// - Oxygen consumption and immunostimulatory factor secretion +/// - Displacement computation applying pushing/adhesive forces between cells +/// - Oncoprotein expression levels +/// - Interactions with CAR-T cells + +class TumorCell : public Cell { + BDM_AGENT_HEADER(TumorCell, Cell, 1); + + public: + TumorCell() = default; + + explicit TumorCell(const Real3& position); + + // Special member functions + TumorCell(const TumorCell&) = default; + TumorCell(TumorCell&&) = default; + ~TumorCell() override = default; + + // NOLINTNEXTLINE(cppcoreguidelines-special-member-functions) + // Assignment operators are implicitly deleted due to base class (Cell/Agent) + + /// Called when a new agent is created (after cell division) + /// @param event The new agent event containing initialization data + void Initialize(const NewAgentEvent& event) override; + + /// Getters and Setters + void SetState(TumorCellState state) { state_ = state; } + TumorCellState GetState() const { return state_; } + + 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 SetFluidFraction(real_t fluid_fraction) { + fluid_fraction_ = fluid_fraction; + } + real_t GetFluidFraction() const { return fluid_fraction_; } + + void SetNuclearVolume(real_t nuclear_volume) { + nuclear_volume_ = nuclear_volume; + } + real_t GetNuclearVolume() const { return nuclear_volume_; } + + void SetTargetCytoplasmSolid(real_t target_cytoplasm_solid) { + target_cytoplasm_solid_ = target_cytoplasm_solid; + } + real_t GetTargetCytoplasmSolid() const { return target_cytoplasm_solid_; } + + void SetTargetNucleusSolid(real_t target_nucleus_solid) { + target_nucleus_solid_ = target_nucleus_solid; + } + real_t GetTargetNucleusSolid() const { return target_nucleus_solid_; } + + void SetTargetFractionFluid(real_t target_fraction_fluid) { + target_fraction_fluid_ = target_fraction_fluid; + } + real_t GetTargetFractionFluid() const { return target_fraction_fluid_; } + + void SetTargetRelationCytoplasmNucleus( + real_t target_relation_cytoplasm_nucleus) { + target_relation_cytoplasm_nucleus_ = target_relation_cytoplasm_nucleus; + } + real_t GetTargetRelationCytoplasmNucleus() const { + return target_relation_cytoplasm_nucleus_; + } + + void SetTransformationRandomRate(); + real_t GetTransformationRandomRate() const { + return transformation_random_rate_; + } + + void SetAttachedToCart(bool attached) { attached_to_cart_ = attached; } + bool IsAttachedToCart() const { return attached_to_cart_; } + + void SetType(TumorCellType type) { type_ = type; } + TumorCellType GetType() const { return type_; } + + Real3 GetOlderVelocity() const { return older_velocity_; } + void SetOlderVelocity(const Real3& velocity) { older_velocity_ = velocity; } + + real_t GetOxygenConsumptionRate() const { return oxygen_consumption_rate_; } + void SetOxygenConsumptionRate(real_t rate) { + oxygen_consumption_rate_ = rate; + } + + real_t GetImmunostimulatoryFactorSecretionRate() const { + return immunostimulatory_factor_secretion_rate_; + } + void SetImmunostimulatoryFactorSecretionRate(real_t rate) { + immunostimulatory_factor_secretion_rate_ = rate; + } + + real_t GetTargetTotalVolume() const; + + int GetTypeAsInt() const { return static_cast(type_); } + + /// Returns the diffusion grid for oxygen + DiffusionGrid* GetOxygenDiffusionGrid() const { return oxygen_dgrid_; } + /// Returns the diffusion grid for immunostimulatory factors + DiffusionGrid* GetImmunostimulatoryFactorDiffusionGrid() const { + return immunostimulatory_factor_dgrid_; + } + + /// Change volume using exponential relaxation equation + /// + /// This method explicitly solves the system of exponential relaxation + /// differential equations using a discrete update step. It is used to grow or + /// shrink the volume (and proportions) smoothly toward a desired target + /// volume over time. The relaxation rate controls the speed of convergence + /// and dt=1 (the time_step). + /// + /// @param relaxation_rate_cytoplasm Relaxation rate for cytoplasm volume + /// changes + /// @param relaxation_rate_nucleus Relaxation rate for nucleus volume changes + /// @param relaxation_rate_fluid Relaxation rate for fluid volume changes + void ChangeVolumeExponentialRelaxationEquation( + real_t relaxation_rate_cytoplasm, real_t relaxation_rate_nucleus, + real_t relaxation_rate_fluid); + + /// Calculate displacement of the cell + /// + /// Computes the displacement of the cell based on interaction forces. + /// + /// @param force Pointer to the interaction force object + /// @param squared_radius The squared radius of the cell + /// @param dt The time step for the simulation + /// @return The calculated displacement vector + Real3 CalculateDisplacement(const InteractionForce* force, + real_t squared_radius, real_t dt) override; + + /// Consume or secrete substances + /// + /// Computes new oxygen or immunostimulatory factor concentration after + /// consumption or secretion by the cell. + /// + /// @param substance_id The ID of the substance (oxygen or immunostimulatory + /// factor) + /// @param old_concentration The previous concentration of the substance + /// @return The new concentration after consumption/secretion + real_t ConsumeSecreteSubstance(int substance_id, real_t old_concentration); + + /// Compute constants for consumption and secretion + /// + /// Updates constants after the cell's change of volume or quantities. + /// These constants are used in the consumption/secretion differential + /// equations. + void ComputeConstantsConsumptionSecretion(); + + private: + /// Current state of the tumor cell + // NOLINTNEXTLINE(readability-identifier-naming) + TumorCellState state_ = TumorCellState::kAlive; + + /// Timer to track time in the current state (in minutes) + // NOLINTNEXTLINE(readability-identifier-naming) + real_t timer_state_ = 0; + + /// Pointer to the oxygen diffusion grid + // NOLINTNEXTLINE(readability-identifier-naming) + DiffusionGrid* oxygen_dgrid_ = nullptr; + + /// Pointer to the immunostimulatory factor diffusion grid + // NOLINTNEXTLINE(readability-identifier-naming) + DiffusionGrid* immunostimulatory_factor_dgrid_ = nullptr; + + /// Level of oncoprotein expression + // NOLINTNEXTLINE(readability-identifier-naming) + real_t oncoproteine_level_ = 0.0; + + /// Transition random rate between states: + /// Affects the probability of transitioning and depends on the individual + /// cell. This rate is kept constant during the cell's lifetime. + // NOLINTNEXTLINE(readability-identifier-naming) + real_t transformation_random_rate_ = 0.0; + + /// Flag indicating if the cell is attached to a CAR-T cell + // NOLINTNEXTLINE(readability-identifier-naming) + bool attached_to_cart_ = false; + + /// Fluid fraction of the cell volume + // NOLINTNEXTLINE(readability-identifier-naming) + real_t fluid_fraction_ = 0.0; + + /// Volume of the nucleus + // NOLINTNEXTLINE(readability-identifier-naming) + real_t nuclear_volume_ = 0.0; + + /// Target cytoplasm solid volume for exponential relaxation + /// + /// Used for growing or shrinking tumor cells. The volume change follows + /// an exponential relaxation equation toward this target volume. + // NOLINTNEXTLINE(readability-identifier-naming) + real_t target_cytoplasm_solid_ = 0.0; + + /// Target nucleus solid volume for exponential relaxation + // NOLINTNEXTLINE(readability-identifier-naming) + real_t target_nucleus_solid_ = 0.0; + + /// Target fluid fraction for exponential relaxation + // NOLINTNEXTLINE(readability-identifier-naming) + real_t target_fraction_fluid_ = 0.0; + + /// Target relation between cytoplasm and nucleus volumes + // NOLINTNEXTLINE(readability-identifier-naming) + real_t target_relation_cytoplasm_nucleus_ = 0.0; + + /// Cell type according to oncoprotein level: + /// Types 1-4: 1 is the most mutated and proliferative type, 4 is the least + /// aggressive. Type 5 means the cell is dead. + // NOLINTNEXTLINE(readability-identifier-naming) + TumorCellType type_ = TumorCellType::kType0; + + /// Velocity of the cell in the previous time step + // NOLINTNEXTLINE(readability-identifier-naming) + Real3 older_velocity_; + + /// Rate of oxygen consumption by the cell + // NOLINTNEXTLINE(readability-identifier-naming) + real_t oxygen_consumption_rate_ = 0.0; + + /// Rate of immunostimulatory factor secretion by the cell + // NOLINTNEXTLINE(readability-identifier-naming) + real_t immunostimulatory_factor_secretion_rate_ = 0.0; + + /// Constant 1 for oxygen consumption/secretion differential equation solution + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant1_oxygen_ = 0.0; + + /// Constant 2 for oxygen consumption/secretion differential equation solution + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant2_oxygen_ = 0.0; + + /// Constant 1 for immunostimulatory factor consumption/secretion differential + /// equation solution + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant1_immunostimulatory_factor_ = 0.0; + + /// Constant 2 for immunostimulatory factor consumption/secretion differential + /// equation solution + // NOLINTNEXTLINE(readability-identifier-naming) + real_t constant2_immunostimulatory_factor_ = 0.0; +}; + +/// Behavior class for controlling tumor cell state transitions and growth +/// +/// This behavior handles the state control logic for tumor cells, managing +/// transitions between different cell states, growth, proliferation, and death +/// processes. It includes logic for determining when cells should enter +/// necrosis based on oxygen levels and other environmental factors. + +struct StateControlGrowProliferate : public Behavior { + BDM_BEHAVIOR_HEADER(StateControlGrowProliferate, Behavior, 1); + + StateControlGrowProliferate() { AlwaysCopyToNew(); } + + // Special member functions + StateControlGrowProliferate(const StateControlGrowProliferate&) = default; + StateControlGrowProliferate& operator=(const StateControlGrowProliferate&) = + default; + StateControlGrowProliferate(StateControlGrowProliferate&&) = default; + StateControlGrowProliferate& operator=(StateControlGrowProliferate&&) = + default; + + ~StateControlGrowProliferate() override = default; + + /// Execute the state control and growth behavior + void Run(Agent* agent) override; + + private: + /// Compute the probability of the cell entering necrosis + /// + /// Determines whether a cell should enter necrosis based on oxygen levels + /// + /// @param oxygen_level Current oxygen concentration at the cell's location + /// @param cell Pointer to the tumor cell being evaluated + /// @return True if the cell should enter necrosis, false otherwise + static bool ShouldEnterNecrosis(real_t oxygen_level, TumorCell* cell); + + /// Manage the behavior of a living tumor cell + /// + /// @param cell Pointer to the tumor cell being managed + /// @param oxygen_level Current oxygen concentration at the cell's location + static void ManageLivingCell(TumorCell* cell, real_t oxygen_level); +}; + +} // namespace bdm + +#endif // TUMOR_CELL_H_ \ No newline at end of file diff --git a/src/utils_aux.cc b/src/utils_aux.cc new file mode 100644 index 0000000..a24138b --- /dev/null +++ b/src/utils_aux.cc @@ -0,0 +1,211 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#include +#include +#include +#include +#include +#include + +#include "core/agent/agent.h" +#include "core/container/math_array.h" +#include "core/real_t.h" +#include "core/resource_manager.h" +#include "core/util/math.h" + +#include "hyperparams.h" +#include "tumor_cell.h" +#include "utils_aux.h" + +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(); + real_t value = random->Gaus(mean, sigma); + if (value < 0.) { + value = 0.; + } + return value; +} + +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; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + + std::vector positions; + + // Hexagonal close-packing spacing + const real_t spacing_x = + cell_radius * + std::sqrt( + 3.0); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + const real_t spacing_y = + cell_radius * + 2.0; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + const real_t spacing_z = + cell_radius * + std::sqrt( + 3.0); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + + // Use integer counters instead of floating-point loop variables + const int z_steps = static_cast((2 * sphere_radius) / spacing_z) + 1; + const int x_steps = static_cast((2 * sphere_radius) / spacing_x) + 1; + const int y_steps = static_cast((2 * sphere_radius) / spacing_y) + 1; + + for (int zi = 0; zi < z_steps; ++zi) { + const real_t z = -sphere_radius + zi * spacing_z; + for (int xi = 0; xi < x_steps; ++xi) { + const real_t x = -sphere_radius + xi * spacing_x; + for (int yi = 0; yi < y_steps; ++yi) { + const real_t y = -sphere_radius + yi * spacing_y; + + // Compute cell center with HCP offset + const real_t px = + x + + (zi % 2) * 0.5 * + cell_radius; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + const real_t py = + y + + (xi % 2) * + cell_radius; // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + const real_t pz = z; + + const real_t dist = std::sqrt(px * px + py * py + pz * pz); + + if (dist <= sphere_radius) { + positions.push_back({px, py, pz}); + } + } + } + } + + 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; + + real_t max_dist_sq = 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 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 + switch (tumor_cell->GetType()) { + case TumorCellType::kType1: + num_tumor_cells_type1++; + break; + case TumorCellType::kType2: + num_tumor_cells_type2++; + break; + case TumorCellType::kType3: + num_tumor_cells_type3++; + break; + case TumorCellType::kType4: + num_tumor_cells_type4++; + break; + case TumorCellType::kType5: + num_tumor_cells_type5_dead++; + break; + default: + break; + } + } + }); + 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)}; +} + +// Function to output summary CSV +void OutputSummary::operator()() { + auto* scheduler = Simulation::GetActive()->GetScheduler(); + auto current_step = scheduler->GetSimulatedSteps(); + + if (current_step % frequency_ == 0) { + std::ofstream file("output/final_data.csv", std::ios::app); + if (file.is_open()) { + 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"; // Header + // for + // CSV + // file + } + + // Calculate time in days, hours, minutes + const double total_minutes = + Simulation::GetActive()->GetScheduler()->GetSimulatedTime(); + 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; + real_t tumor_radius = 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(); + + // Write data to CSV file + file << total_days << "," << total_hours << "," << total_minutes << "," + << tumor_radius << "," + << Simulation::GetActive()->GetResourceManager()->GetNumAgents() + << "," // number of 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 + << "\n"; + } + } +} + +} // namespace bdm \ No newline at end of file diff --git a/src/utils_aux.h b/src/utils_aux.h new file mode 100644 index 0000000..f94486e --- /dev/null +++ b/src/utils_aux.h @@ -0,0 +1,107 @@ +/* + * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * 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 obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * SPDX-License-Identifier: Apache-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * This file contains a model developed under Google Summer of Code (GSoC) + * for the compiler-research.org organization. + */ + +#ifndef CORE_UTIL_UTILS_AUX_H_ +#define CORE_UTIL_UTILS_AUX_H_ + +#include +#include +#include +#include + +#include "core/container/math_array.h" +#include "core/operation/operation.h" +#include "core/operation/operation_registry.h" +#include "core/real_t.h" + +namespace bdm { +/// Forward declaration of TumorCell class +class TumorCell; + +/// Sample a positive Gaussian value +/// +/// Samples a Gaussian value with given mean and standard deviation. +/// All negative values are mapped to zero to ensure positive results. +/// +/// @param mean Mean value of the Gaussian distribution +/// @param sigma Standard deviation of the Gaussian distribution +/// @return Sampled positive value (negative values mapped to zero) +real_t SamplePositiveGaussian(float mean, float sigma); + +/// Create a spherical arrangement of tumor cells +/// +/// Generates a vector of 3D positions for tumor cells arranged in a spherical +/// pattern with the specified radius. The cells are positioned to form a +/// forces-stable structure. The achieved tumor mimics a heterogeneous organoid +/// thanks to the different cancer cell types that can be randomly generated +/// +/// @param sphere_radius Radius of the spherical tumor in micrometers +/// @return Vector of 3D positions where tumor cells should be placed +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. +/// +/// @return Tuple containing: +/// - 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 +/// - Current tumor radius in micrometers +std::tuple +ComputeNumberTumorCellsAndRadius(); + +/// Operation for outputting simulation summary data to CSV files +/// +/// 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. +struct OutputSummary : public StandaloneOperationImpl { + BDM_OP_HEADER(OutputSummary); + + public: + void SetFrequency(uint64_t frequency) { frequency_ = frequency; } + [[nodiscard]] uint64_t GetFrequency() const { return frequency_; } + + private: + /// Frequency of output (every N simulation steps) + // NOLINTNEXTLINE(readability-identifier-naming) + 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; +}; + +/// Register OutputSummary operation with CPU as compute target +inline BDM_REGISTER_OP(OutputSummary, "OutputSummary", kCpu); + +} // namespace bdm + +#endif // CORE_UTIL_UTILS_AUX_H_ \ No newline at end of file diff --git a/test/test-main.cc b/test/test-main.cc new file mode 100644 index 0000000..c47c7cd --- /dev/null +++ b/test/test-main.cc @@ -0,0 +1,21 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- + +#include + +int main(int argc, char** argv) { + ::testing::FLAGS_gtest_death_test_style = "threadsafe"; + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/test/test-suite.cc b/test/test-suite.cc new file mode 100644 index 0000000..b671c8a --- /dev/null +++ b/test/test-suite.cc @@ -0,0 +1,60 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2021 CERN & University of Surrey for the benefit of the +// BioDynaMo collaboration. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// +// See the LICENSE file distributed with this work for details. +// See the NOTICE file distributed with this work for additional information +// regarding copyright ownership. +// +// ----------------------------------------------------------------------------- + +#include +#include "biodynamo.h" + +// Googletest in combination with the provided CMakeLists.txt allows you to +// define tests in arbitrary .cc files in the `test/` folder. This file should +// serve as an inspiration for testing user-defined, custom behaviors, basic as +// well as compicated functions, or similar things. If you wish to add tests for +// specific aspects, you can either add them to the existing test-suite.cc file +// or create a new *.cc file in the `test/` folder. CMake will handle it +// automatically. For more information regarding testing with Googletest, please +// consider the following sources: +// * https://google.github.io/googletest/primer.html +// * https://github.com/google/googletest + +#define TEST_NAME typeid(*this).name() + +namespace bdm { + +// A function to test +int Compute42() { return 6 * 7; }; + +// Show how to compare two numbers +TEST(UtilTest, NumberTest) { + // Expect equality + EXPECT_EQ(Compute42(), 42); +} + +// Test if we can add agents to the simulation +TEST(AgentTest, AddAgentsToSimulation) { + // Create simulation + Simulation simulation(TEST_NAME); + + // Add some cells to the simulation + auto* rm = simulation.GetResourceManager(); + uint8_t expected_no_cells{20}; + for (int i = 0; i < expected_no_cells; i++) { + auto* cell = new Cell(30); + rm->AddAgent(cell); + } + + // Test if all 20 cells are in the simulation + auto no_cells = rm->GetNumAgents(); + EXPECT_EQ(expected_no_cells, no_cells); +} + +} // namespace bdm