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..2db1355 --- /dev/null +++ b/bdm.toml @@ -0,0 +1,10 @@ +[visualization] +export = true +interval = 7200 + +[[visualize_agent]] +name = "TumorCell" +additional_data_members = ["diameter_","volume_", "type_"] + +[[visualize_diffusion]] +name = "oxygen" diff --git a/src/cart_cell.cc b/src/cart_cell.cc new file mode 100644 index 0000000..c94961a --- /dev/null +++ b/src/cart_cell.cc @@ -0,0 +1,275 @@ +/* + * 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 "cart_cell.h" + +namespace bdm { + +CartCell::CartCell(const Real3& position) { + SetPosition(position); + // Default state for new cells + state_ = CartCellState::kAlive; + // Initial timer_state for apoptotic state + timer_state_ = 0; + + //volumes + // Set default volume + SetVolume(kDefaultVolumeNewCartCell); + // Set default fluid fraction + SetFluidFraction(kDefaultFractionFluidCartCell); + // Set default nuclear volume + SetNuclearVolume(kDefaultVolumeNucleusCartCell); + + + + ResourceManager &rm = *Simulation::GetActive()->GetResourceManager(); + // Pointer to oxygen diffusion grid + oxygen_dgrid_ = rm.GetDiffusionGrid("oxygen"); + // Pointer to immunostimulatory_factor diffusion grid + immunostimulatory_factor_dgrid_ = rm.GetDiffusionGrid("immunostimulatory_factor"); + // Initially not attached to a tumor cell + attached_to_tumor_cell_ = false; + // Initialize attached cell pointer to null + attached_cell_ = nullptr; + + // Initialize the velocity of the cell in the previous step to zero + older_velocity_ = {0, 0, 0}; + + + SetCurrentLiveTime(kAverageMaximumTimeUntillApoptosisCart); + + //Add Consumption and Secretion + // Set default oxygen consumption rate + SetOxygenConsumptionRate(kDefaultOxygenConsumption); + // Compute constants for all ConsumptionSecretion of Oxygen + 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() { + 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 +void CartCell::ChangeVolumeExponentialRelaxationEquation(real_t relaxation_rate_cytoplasm, real_t relaxation_rate_nucleus, real_t relaxation_rate_fluid) { + // Exponential relaxation towards the target volume + real_t current_total_volume = GetVolume(); + real_t fluid_fraction= GetFluidFraction(); + real_t nuclear_volume = GetNuclearVolume(); + + real_t current_nuclear_solid = nuclear_volume * (1 - fluid_fraction); + real_t current_cytoplasm_solid = (current_total_volume - nuclear_volume) * (1-fluid_fraction); + + 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; } + + 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; } + + 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; } + + real_t new_total_solid= nuclear_solid + cytoplasm_solid; + + real_t total_nuclear= nuclear_solid + nuclear_fluid; + + // real_t total_cytoplasm= cytoplasm_solid + cytoplasm_fluid; + + real_t new_volume = new_total_solid + new_fluid; + + // Avoid division by zero + real_t new_fraction_fluid = new_fluid / (1e-16 + 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=kSquaredMaxDistanceNeighborsForce; + + // 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; + 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 +real_t CartCell::ConsumeSecreteSubstance(int substance_id, real_t old_concentration) { + real_t res; + 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 + 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(), 1e-10)) { // 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(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 diff --git a/src/cart_cell.h b/src/cart_cell.h new file mode 100644 index 0000000..06e00b0 --- /dev/null +++ b/src/cart_cell.h @@ -0,0 +1,224 @@ +/* + * 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 "biodynamo.h" +#include "core/util/log.h" +#include "core/util/root.h" +#include "utils_aux.h" +#include "hyperparams.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() {} + explicit CartCell(const Real3& position); + virtual ~CartCell() {} + + ///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(); + + /// 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 CAR-T cell + CartCellState state_; + + /// Timer to track time in the current state (in minutes) + /// Used for apoptotic state timing + int timer_state_; + + /// Pointer to the oxygen diffusion grid + DiffusionGrid* oxygen_dgrid_; + + /// Pointer to the immunostimulatory factor diffusion grid + DiffusionGrid* immunostimulatory_factor_dgrid_; + + /// Flag indicating if the cell is attached to a tumor cell + bool attached_to_tumor_cell_; + + /// Current time until apoptosis + real_t current_live_time_; + + /// Fluid fraction of the cell volume + real_t fluid_fraction_; + + /// Volume of the nucleus + real_t nuclear_volume_; + + /// Target cytoplasm solid volume for exponential relaxation + /// Used during volume changes following exponential relaxation equation + real_t target_cytoplasm_solid_; + + /// Target nucleus solid volume for exponential relaxation + real_t target_nucleus_solid_; + + /// Target fluid fraction for exponential relaxation + real_t target_fraction_fluid_; + + /// Target relation between cytoplasm and nucleus volumes + real_t target_relation_cytoplasm_nucleus_; + + /// Velocity of the cell in the previous time step + Real3 older_velocity_; + + /// Rate of oxygen consumption by the cell + real_t oxygen_consumption_rate_; + + /// Rate of immunostimulatory factor secretion by the cell + real_t immunostimulatory_factor_secretion_rate_; + + /// Constant 1 for oxygen consumption/secretion differential equation solution + real_t constant1_oxygen_; + + /// Constant 2 for oxygen consumption/secretion differential equation solution + real_t constant2_oxygen_; + + /// Pointer to the attached tumor cell + TumorCell* attached_cell_; + +}; + +/// 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(); } + + virtual ~StateControlCart() {} + + /// 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..7e30490 --- /dev/null +++ b/src/cart_tumor.cc @@ -0,0 +1,117 @@ +/* + * 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 "cart_tumor.h" +#include "tumor_cell.h" +#include "cart_cell.h" +#include "diffusion_thomas_algorithm.h" +#include "forces_tumor_cart.h" +#include "core/environment/uniform_grid_environment.h" +#include "core/operation/mechanical_forces_op.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; + param->max_bound = kBoundedSpaceLength/2; // 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]; + op->GetImplementation()->SetInteractionForce(new InteractionVelocity()); + + auto* env = dynamic_cast(Simulation::GetActive()->GetEnvironment()); + // Fix the box length for the uniform grid environment + env->SetBoxLength(kLengthBoxMechanics); + + // Define Substances + auto* rm = Simulation::GetActive()->GetResourceManager(); + + // Oxygen + // substance_id, name, diffusion_coefficient, decay_constant, resolution, time_step + auto* oxygen_grid = new DiffusionThomasAlgorithm( + kOxygen, "oxygen", + kDiffusionCoefficientOxygen,// 100000 micrometers^2/minute + kDecayConstantOxygen, // 0.1 minutes^-1 + kResolutionGridSubstances, + kDtSubstances, + true); // true indicates Dirichlet border conditions + rm->AddContinuum(oxygen_grid); + + // Immunostimulatory Factor + // substance_id, name, diffusion_coefficient, decay_constant, resolution + auto* immunostimulatory_factor_grid = new DiffusionThomasAlgorithm( + 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); + + // 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 + 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 + std::vector positions=CreateSphereOfTumorCells(kInitialRadiusTumor);//positions of the cells + for (const auto& pos : positions) { + TumorCell* tumor_cell = new TumorCell(pos); + tumor_cell->AddBehavior(new StateControlGrowProliferate()); + ctxt->AddAgent(tumor_cell); + } + + //OutputSummary operation + auto* summary_op = new bdm::Operation("OutputSummary"); + // Set the interval for outputting CSV files + summary_op->frequency_ = kOutputCsvInterval; + summary_op->AddOperationImpl(bdm::kCpu, new bdm::OutputSummary()); + scheduler->ScheduleOp(summary_op); + + // 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); } diff --git a/src/cart_tumor.h b/src/cart_tumor.h new file mode 100644 index 0000000..a214602 --- /dev/null +++ b/src/cart_tumor.h @@ -0,0 +1,35 @@ +/* + * 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_ + +#include "biodynamo.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..6c6962c --- /dev/null +++ b/src/diffusion_thomas_algorithm.cc @@ -0,0 +1,282 @@ +/* + * 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 "diffusion_thomas_algorithm.h" +#include "core/resource_manager.h" +#include "core/simulation.h" +#include "hyperparams.h" +#include "tumor_cell.h" +#include "cart_cell.h" + +namespace bdm { + + +DiffusionThomasAlgorithm::DiffusionThomasAlgorithm(int substance_id, std::string substance_name, real_t dc, real_t mu,int resolution, real_t dt, bool dirichlet_border)//time step + : DiffusionGrid(substance_id, std::move(substance_name), dc, mu, resolution) { + + SetTimeStep(dt); + //num of voxels in each direction + resolution_ = GetResolution(); + // Voxel side length in micrometers + d_space_ = kBoundedSpaceLength / resolution_; + + dirichlet_border_ = dirichlet_border; + + jump_i_ = 1; + jump_j_ = resolution_; + jump_k_ = resolution_ * resolution_; + + //all diffusion coefficients are the same for all directions (isotropic) + constant1_ = dc; + constant1_ *=dt/(d_space_ * d_space_); + constant1a_ = -constant1_; + //decay constant + constant2_ = mu; + // Divide by 3 for the three directions + constant2_ *= dt / 3.0; + + constant3_ = 1.0 + 2 * constant1_ + constant2_; + constant3a_ = 1.0 + constant1_ + constant2_; + + // Initialize the denominators and coefficients for the Thomas algorithm + + thomas_c_x_ = std::vector(resolution_, constant1a_); + thomas_denom_x_ = std::vector(resolution_, constant3_); + InitializeThomasAlgorithmVectors(thomas_denom_x_, thomas_c_x_); + + thomas_c_y_ = std::vector(resolution_, constant1a_); + thomas_denom_y_ = std::vector(resolution_, constant3_); + InitializeThomasAlgorithmVectors(thomas_denom_y_, thomas_c_y_); + + thomas_c_z_ = std::vector(resolution_, constant1a_); + thomas_denom_z_ = std::vector(resolution_, constant3_); + InitializeThomasAlgorithmVectors(thomas_denom_z_, thomas_c_z_); + +} + +void DiffusionThomasAlgorithm::InitializeThomasAlgorithmVectors(std::vector& thomas_denom, std::vector& thomas_c) { + 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 (unsigned int 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() { + real_t origin= GetDimensionsPtr()[0]; + 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++) { + real_t real_x = origin + x * d_space_; + real_t real_y = origin + y * d_space_; + //For z=0 + size_t z=0; + real_t real_z = origin + z * d_space_; + SetConcentration(x, y, z, GetBoundaryCondition()->Evaluate(real_x, real_y, real_z, simulated_time)); + //For z=resolution_-1 + z = resolution_ - 1; + real_z = origin + z * d_space_; + SetConcentration(x, y, 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++) { + real_t real_x = origin + x * d_space_; + real_t real_z = origin + z * d_space_; + //For y=0 + size_t y=0; + real_t real_y = origin + y * d_space_; + SetConcentration(x, y, z, GetBoundaryCondition()->Evaluate(real_x, real_y, real_z, simulated_time)); + //For y=resolution_-1 + y = resolution_ - 1; + real_y = origin + y * d_space_; + SetConcentration(x, y, 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++) { + real_t real_y = origin + y * d_space_; + real_t real_z = origin + z * d_space_; + //For x=0 + size_t x=0; + real_t real_x = origin + x * d_space_; + SetConcentration(x, y, z, GetBoundaryCondition()->Evaluate(real_x, real_y, real_z, simulated_time)); + //For x=resolution_-1 + x = resolution_ - 1; + real_x = origin + x * d_space_; + SetConcentration(x, y, 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){ + ChangeConcentrationBy(idx, amount - GetAllConcentrations()[idx], 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) {//instead of overwriting Step, in future versions of BioDynaMo, we should overwrite CheckParameters + // check if diffusion coefficient and decay constant are 0 + // i.e. if we don't need to calculate diffusion update + if (IsFixedSubstance()) { + return; + } + DiffuseChemical(dt); + + //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(real_t dt) { + + + // Change for the future: to add double buffer for paralelization + + if (dirichlet_border_) { ApplyDirichletBoundaryConditions();} + + //X-direction + #pragma omp parallel for collapse(2) + for( unsigned int k=0; k < resolution_; k++) { + for( unsigned int j=0; j < resolution_; j++) { + int ind = GetBoxIndex(0, j, k); + + SetConcentration(ind, GetAllConcentrations()[ind]/thomas_denom_x_[0]); + // Forward elimination step for x direction + for (unsigned int i = 1; i < resolution_ ; i++) { + ind = GetBoxIndex(i, j, k); + auto* all_concentrations = GetAllConcentrations(); + SetConcentration(ind, (all_concentrations[ind] + constant1_ * all_concentrations[ind-jump_i_]) / thomas_denom_x_[i]); + } + // Back substitution step for x direction + for (int i = resolution_ - 2; i >= 0; i--) { + ind = GetBoxIndex(i, j, k); + auto* all_concentrations = GetAllConcentrations(); + SetConcentration(ind, all_concentrations[ind] - thomas_c_x_[i] * all_concentrations[ind + jump_i_]); + } + } + } + + if (dirichlet_border_) { ApplyDirichletBoundaryConditions();} + + //Y-direction + #pragma omp parallel for collapse(2) + for( unsigned int k=0; k < resolution_; k++) { + for( unsigned int i=0; i < resolution_; i++) { + int ind = GetBoxIndex(i, 0, k); + + SetConcentration(ind, GetAllConcentrations()[ind]/thomas_denom_y_[0]); + // Forward elimination step for y direction + for (unsigned int j = 1; j < resolution_ ; j++) { + ind = GetBoxIndex(i, j, k); + auto* all_concentrations = GetAllConcentrations(); + SetConcentration(ind, (all_concentrations[ind] + constant1_ * all_concentrations[ind-jump_j_]) / thomas_denom_y_[j]); + } + // Back substitution step for y direction + for (int j = resolution_ - 2; j >= 0; j--) { + ind = GetBoxIndex(i, j, k); + auto* all_concentrations = GetAllConcentrations(); + SetConcentration(ind, all_concentrations[ind] - thomas_c_y_[j] * all_concentrations[ind + jump_j_]); + } + } + } + + if (dirichlet_border_) { ApplyDirichletBoundaryConditions();} + + //Z-direction + #pragma omp parallel for collapse(2) + for( unsigned int j=0; j < resolution_; j++) { + for( unsigned int i=0; i < resolution_; i++) { + int ind = GetBoxIndex(i, j, 0); + SetConcentration(ind, GetAllConcentrations()[ind]/thomas_denom_z_[0]); + // Forward elimination step for z direction + for (unsigned int k = 1; k < resolution_ ; k++) { + ind = GetBoxIndex(i, j, k); + auto* all_concentrations = GetAllConcentrations(); + SetConcentration(ind, (all_concentrations[ind] + constant1_ * all_concentrations[ind-jump_k_]) / thomas_denom_z_[k]); + } + // Back substitution step for z direction + for (int k = resolution_ - 2; k >= 0; k--) { + ind = GetBoxIndex(i, j, k); + auto* all_concentrations = GetAllConcentrations(); + SetConcentration(ind, all_concentrations[ind] - thomas_c_z_[k] * all_concentrations[ind + jump_k_]); + } + } + } + if (dirichlet_border_) { ApplyDirichletBoundaryConditions(); } + + //Change of concentration levels because of agents + ComputeConsumptionsSecretions(); + + return; + + +} + + +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(); + real_t current_time = GetSimulatedTime(); + //in a future version of BioDynaMo this should be parallelized getting the agents inside each chemical voxel and trating each voxel independently. + rm->ForEachAgent([this, current_time](bdm::Agent* agent) { + if (auto* cell = dynamic_cast(agent)) { + // Handle TumorCell agents + const auto& pos = cell->GetPosition(); + real_t conc = this->GetValue(pos); + 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(); + real_t conc = GetValue(pos); + real_t new_conc = cell->ConsumeSecreteSubstance(GetContinuumId(),conc); + ChangeConcentrationBy(pos, new_conc - conc, InteractionMode::kAdditive, false); + } + + }); + + return; +} + +} // namespace bdm diff --git a/src/diffusion_thomas_algorithm.h b/src/diffusion_thomas_algorithm.h new file mode 100644 index 0000000..32c92d4 --- /dev/null +++ b/src/diffusion_thomas_algorithm.h @@ -0,0 +1,193 @@ +/* + * 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 "core/diffusion/diffusion_grid.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() = default; + + DiffusionThomasAlgorithm(int substance_id, + std::string substance_name, + real_t dc, + real_t mu, + int 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(x, y, 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. + /// + /// @param dt Time step for the diffusion computation + void DiffuseChemical(real_t dt); + + /// 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 + size_t resolution_; + + /// Voxel side length in micrometers + real_t d_space_; + + /// Denominators for x-direction Thomas algorithm + std::vector thomas_denom_x_; + + /// Coefficients for x-direction Thomas algorithm + std::vector thomas_c_x_; + + /// Denominators for y-direction Thomas algorithm + std::vector thomas_denom_y_; + + /// Coefficients for y-direction Thomas algorithm + std::vector thomas_c_y_; + + /// Denominators for z-direction Thomas algorithm + std::vector thomas_denom_z_; + + /// Coefficients for z-direction Thomas algorithm + std::vector thomas_c_z_; + + /// Index jump for i-direction (x-axis) + int jump_i_; + + /// Index jump for j-direction (y-axis) + int jump_j_; + + /// Index jump for k-direction (z-axis) + int jump_k_; + + /// First diffusion constant + real_t constant1_; + + /// Alternative first diffusion constant + real_t constant1a_; + + /// Second diffusion constant + real_t constant2_; + + /// Third diffusion constant + real_t constant3_; + + /// Alternative third diffusion constant + real_t constant3a_; + + /// Flag indicating Dirichlet boundary conditions + bool dirichlet_border_; + + /// 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); + + /// 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; + + + BDM_CLASS_DEF_OVERRIDE(DiffusionThomasAlgorithm, 1); +}; + +} // namespace bdm + +#endif // DIFFUSION_THOMAS_ALGORITHM_H_ diff --git a/src/forces_tumor_cart.cc b/src/forces_tumor_cart.cc new file mode 100644 index 0000000..463124e --- /dev/null +++ b/src/forces_tumor_cart.cc @@ -0,0 +1,160 @@ +/* + * 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 "forces_tumor_cart.h" + + +namespace bdm { + +Real4 InteractionVelocity::Calculate(const Agent* lhs, const Agent* rhs) const { + + auto* a = dynamic_cast(lhs); + 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)); + + double dist_sq = displacement[0] * displacement[0] + + displacement[1] * displacement[1] + + displacement[2] * displacement[2]; + double distance = std::max(std::sqrt(dist_sq), 1e-5); + + double radius_a = a->GetDiameter() / 2.0; + double radius_b = b->GetDiameter() / 2.0; + double R = radius_a + radius_b; + // R=16.8254;//Debug + // std::cout << "Debug: R = " << R << ", distance = " << distance << std::endl;// Debug output + double temp_r = 0.0; + + const TumorCell* a_tumor = dynamic_cast(a); + const TumorCell* b_tumor = dynamic_cast(b); + + if (distance < R) { + + // 1 - d/R + temp_r = 1.0 - distance / R; + // (1 - d/R)^2 + temp_r *= temp_r; + + double repulsion; + // std::cout << "temp_r = " << temp_r<< std::endl;// Debug output + + + if (a_tumor && b_tumor) {// two tumor cells + repulsion = kRepulsionTumorTumor;//std::sqrt(kRepulsionTumorTumor * kRepulsionTumorTumor); + } else if (!a_tumor && !b_tumor) {// two CAR-T cells + repulsion = kRepulsionCartCart;//std::sqrt(kRepulsionCartCart*kRepulsionCartCart); + } else {// one tumor cell and one CAR-T + repulsion = std::sqrt(kRepulsionCartTumor * + kRepulsionTumorCart); + } + // std::cout << "repulsion = " << repulsion<< std::endl;// Debug output + + temp_r *= repulsion; + } + + // std::cout << "temp_r after repulsion = " << temp_r<< std::endl;// Debug output + + + // Adhesion + double max_interaction_distance = kMaxRelativeAdhesionDistance * R; + // 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 + double 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 + + + double adhesion; + if (a_tumor && b_tumor) {// two tumor cells + adhesion = kAdhesionTumorTumor; + } else if (!a_tumor && !b_tumor) {// 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) < 1e-16) { + return {0.0, 0.0, 0.0, 0.0}; + } + double force_magnitude = temp_r / distance; + + + + //Debug Output volcities + // std::ofstream file("output/intercation_velocities.csv", std::ios::app); + // if (file.is_open()) { + + // double 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 new InteractionVelocity(); +} + +} // 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..b167b11 --- /dev/null +++ b/src/forces_tumor_cart.h @@ -0,0 +1,62 @@ +/* + * 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/interaction_force.h" +#include "core/operation/mechanical_forces_op.h" +#include "biodynamo.h" +#include "core/util/log.h" +#include "core/util/root.h" +#include "hyperparams.h" +#include "utils_aux.h" +#include "tumor_cell.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() 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; + + 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..0c47ecc --- /dev/null +++ b/src/hyperparams.h @@ -0,0 +1,195 @@ +/* + * 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 + +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; + + +///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; + +/// +/// 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; + +/// 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 = (kBoundedSpaceLength / kResolutionGridSubstances)*(kBoundedSpaceLength / kResolutionGridSubstances)*(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; + +///Do not change this line +const real_t kLengthBoxMechanics =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 kSquaredMaxDistanceNeighborsForce = 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..3f76e9f --- /dev/null +++ b/src/tumor_cell.cc @@ -0,0 +1,509 @@ +/* + * 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 "tumor_cell.h" + +namespace bdm { + +TumorCell::TumorCell(const Real3& position) { + SetPosition(position); + state_ = TumorCellState::kAlive; // Default state for new cells + timer_state_ = 0; // Initial timer_state + + //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) / ( 1e-16 + 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 + attached_to_cart_ = false; // Initially not attached to a cart + + older_velocity_ = {0, 0, 0}; // Initialize the velocity of the cell in the previous step to zero + + //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 + real_t new_nuclear_volume = mother->GetNuclearVolume() / 2.0; // Divide mother's nuclear volume by 2 + 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 >= 1.5) {//between 1.5 and 2.0 + type_ = 1; + } else if (level >= 1.0 && level < 1.5) { + type_ = 2; + } else if (level >= 0.5 && level < 1.0) { + type_ = 3; + } else {//between 0.0 and 0.5 + type_ = 4; + } +} + +void TumorCell::SetTransformationRandomRate() { + // avoid division by zero + transformation_random_rate_ = 1/(std::max(SamplePositiveGaussian(38.6,3.7)*60., 1e-16)); +} + +real_t TumorCell::GetTargetTotalVolume() { + 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 +void TumorCell::ChangeVolumeExponentialRelaxationEquation(real_t relaxation_rate_cytoplasm, real_t relaxation_rate_nucleus, real_t relaxation_rate_fluid) { + // Exponential relaxation towards the target volume + real_t current_total_volume = GetVolume(); + real_t fluid_fraction= GetFluidFraction(); + real_t nuclear_volume = GetNuclearVolume(); + + real_t current_nuclear_solid = nuclear_volume * (1 - fluid_fraction); + 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; + + 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; } + + 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; } + + 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; } + + real_t new_total_solid= nuclear_solid + cytoplasm_solid; + + real_t total_nuclear= nuclear_solid + nuclear_fluid; + + // real_t total_cytoplasm= cytoplasm_solid + cytoplasm_fluid; + + real_t new_volume = new_total_solid + new_fluid; + + // Avoid division by zero + real_t new_fraction_fluid = new_fluid / (1e-16 + 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 +real_t TumorCell::ConsumeSecreteSubstance(int substance_id, real_t old_concentration) { + // constant1_oxygen_ = 0; // Debug + // constant2_oxygen_ = 1.3; // Debug + real_t res; + 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 + + 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 + Real3 current_position = cell->GetPosition(); + auto* oxygen_dgrid = cell->GetOxygenDiffusionGrid(); // Pointer to the oxygen diffusion grid + real_t oxygen_level = oxygen_dgrid->GetValue(current_position); + // oxygen_level = 30.; // Debug + + // Debug + // std::cout << oxygen_level << std::endl; + + real_t multiplier; + real_t final_rate_transition; + + 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 + } + + //volume change + cell->ChangeVolumeExponentialRelaxationEquation(kVolumeRelaxarionRateAliveCytoplasm, + kVolumeRelaxarionRateAliveNucleus, + kVolumeRelaxarionRateAliveFluid); // The cell grows to real_t its size + //cell state control + multiplier = 1.0; // Default multiplier for transition cycle + 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 + + + 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 + + //Debug Debug Output params + // std::ofstream file2("output/params_o2_oncoproteine.csv", std::ios::app); + // if (file2.is_open()) { + + // // Write data to CSV file + // file2 << currennt_time << ",multiplier1," + // << multiplier1 << ",multiplier2," + // << multiplier2 << ",transition_rate," + // << final_rate_transition + // <<"\n"; + // } + // End Debug Output + //End Debug + + // //volume change + // cell->ChangeVolumeExponentialRelaxationEquation(kVolumeRelaxarionRateAliveCytoplasm, + // kVolumeRelaxarionRateAliveNucleus, + // kVolumeRelaxarionRateAliveFluid); // The cell grows to real_t its size + // //cell state control + // multiplier = 1.0; // Default multiplier for transition cycle + // 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 + // } + + + // final_rate_transition= cell->GetTransformationRandomRate() * multiplier * cell->GetOncoproteineLevel(); // Calculate the rate of state change based on oxygen level and oncoproteine (min^-1) + + real_t time_to_wait=1e100; // 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 + } + 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"); + } +} + +// computes the probability of the cell entering necrosis +bool StateControlGrowProliferate::ShouldEnterNecrosis(real_t oxygen_level,TumorCell* cell) const { + //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; + } + + 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(); + 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); + cell->SetType(5); // Set type to 5 to indicate dead cell + } + 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..156f11e --- /dev/null +++ b/src/tumor_cell.h @@ -0,0 +1,264 @@ +/* + * 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 "biodynamo.h" +#include "core/util/log.h" +#include "core/util/root.h" +#include "hyperparams.h" +#include "utils_aux.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. +}; + +/// 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() {} + + explicit TumorCell(const Real3& position); + + virtual ~TumorCell() {} + + /// 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(int timer_state) { timer_state_ = timer_state; } + int 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(int type) { type_ = type; } + int 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(); + + /// 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 + TumorCellState state_; + + /// Timer to track time in the current state (in minutes) + int timer_state_; + + /// Pointer to the oxygen diffusion grid + DiffusionGrid* oxygen_dgrid_; + + /// Pointer to the immunostimulatory factor diffusion grid + DiffusionGrid* immunostimulatory_factor_dgrid_; + + /// Level of oncoprotein expression + real_t oncoproteine_level_; + + /// 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. + real_t transformation_random_rate_; + + /// Flag indicating if the cell is attached to a CAR-T cell + bool attached_to_cart_; + + /// Fluid fraction of the cell volume + real_t fluid_fraction_; + + /// Volume of the nucleus + real_t nuclear_volume_; + + /// 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. + real_t target_cytoplasm_solid_; + + /// Target nucleus solid volume for exponential relaxation + real_t target_nucleus_solid_; + + /// Target fluid fraction for exponential relaxation + real_t target_fraction_fluid_; + + /// Target relation between cytoplasm and nucleus volumes + real_t target_relation_cytoplasm_nucleus_; + + /// 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. + int type_; + + /// Velocity of the cell in the previous time step + Real3 older_velocity_; + + /// Rate of oxygen consumption by the cell + real_t oxygen_consumption_rate_; + + /// Rate of immunostimulatory factor secretion by the cell + real_t immunostimulatory_factor_secretion_rate_; + + /// Constant 1 for oxygen consumption/secretion differential equation solution + real_t constant1_oxygen_; + + /// Constant 2 for oxygen consumption/secretion differential equation solution + real_t constant2_oxygen_; + + /// Constant 1 for immunostimulatory factor consumption/secretion differential equation solution + real_t constant1_immunostimulatory_factor_; + + /// Constant 2 for immunostimulatory factor consumption/secretion differential equation solution + real_t constant2_immunostimulatory_factor_; + +}; + +/// 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(); } + + virtual ~StateControlGrowProliferate() {} + + /// 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 + bool ShouldEnterNecrosis(real_t oxygen_level,TumorCell* cell) const; +}; + +} // 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..15594b4 --- /dev/null +++ b/src/utils_aux.cc @@ -0,0 +1,141 @@ +/* + * 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 "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 + real_t cell_radius = std::cbrt(kDefaultVolumeNewTumorCell * 6 / Math::kPi)/2; + + std::vector positions; + + // Hexagonal close-packing spacing + real_t spacing_x = cell_radius * std::sqrt(3.0); + real_t spacing_y = cell_radius * 2.0; + real_t spacing_z = cell_radius * std::sqrt(3.0); + + int zc = 0; + for (real_t z = -sphere_radius; z < sphere_radius; z += spacing_z, ++zc) { + int xc = 0; + for (real_t x = -sphere_radius; x < sphere_radius; x += spacing_x, ++xc) { + int yc = 0; + for (real_t y = -sphere_radius; y < sphere_radius; y += spacing_y, ++yc) { + + // Compute cell center with HCP offset + real_t px = x + (zc % 2) * 0.5 * cell_radius; + real_t py = y + (xc % 2) * cell_radius; + real_t pz = z; + + 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 (auto* tumor_cell = dynamic_cast(agent)) { + total_num_tumor_cells++; + const auto& pos = agent->GetPosition(); + 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 1: num_tumor_cells_type1++; break; + case 2: num_tumor_cells_type2++; break; + case 3: num_tumor_cells_type3++; break; + case 4: num_tumor_cells_type4++; break; + case 5: 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 + double total_minutes = Simulation::GetActive()->GetScheduler()->GetSimulatedTime(); + double total_hours = total_minutes / 60.0; + double total_days = total_hours / 24.0; + + // Count total cells, tumor cells of each type and tumor radius + size_t total_num_tumor_cells; + size_t num_tumor_cells_type1, num_tumor_cells_type2, num_tumor_cells_type3, num_tumor_cells_type4, num_tumor_cells_type5_dead; + real_t tumor_radius; + 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..8bb5f3e --- /dev/null +++ b/src/utils_aux.h @@ -0,0 +1,93 @@ +/* + * 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 "biodynamo.h" +#include "core/util/random.h" +#include "hyperparams.h" +#include "tumor_cell.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); + + /// Frequency of output (every N simulation steps) + uint64_t frequency_ = 1; + + /// Collects current simulation data and writes it to CSV files + /// + /// Called automatically by the simulation scheduler at the specified frequency. + /// Gathers statistics about cell populations, tumor radius, and other metrics, + /// then outputs them to appropriately named CSV files for analysis. + void operator()() override; + }; + + /// 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