diff --git a/.gitignore b/.gitignore index d5b1d28..8ff13e2 100644 --- a/.gitignore +++ b/.gitignore @@ -117,9 +117,7 @@ doc/latex *.csv *.xyz *.dat - -# Draft code for comparing models -draft_code_my_own_analysis/ +auxiliar_scripts/out/*.png # Python environment .venv/ @@ -129,3 +127,6 @@ cout.txt #sqlite database files *.db + +#auxiliar_scripts +auxiliar_scripts/ \ No newline at end of file diff --git a/README.md b/README.md index 1675b4f..f019077 100644 --- a/README.md +++ b/README.md @@ -33,9 +33,10 @@ 8. [Input Parameters](#input-parameters) 9. [Running the Simulation](#running-the-simulation) 10. [Visualizing Results](#visualizing-results) -11. [Acknowledgments](#acknowledgments) -12. [License](#license) -13. [Author Contact Information](#author-contact-information) +11. [Model Calibration (Bayesian Optimization)](#model-calibration-bayesian-optimization) +12. [Acknowledgments](#acknowledgments) +13. [License](#license) +14. [Author Contact Information](#author-contact-information) --- @@ -70,7 +71,40 @@ Together, these features make **CARTopiaX** a powerful computational tool for in The diagram below illustrates the **agents and main mechanics** interacting within CARTopiaX: -![Model Diagram](./assets/images/model_outline.png) + + +```mermaid +--- +config: + flowchart: + curve: basis +--- +flowchart LR + CART("CAR-T Cell") + IMMUNO(("Immunostimulatory
factor")) + OXY(("Oxygen")) + GLU(("Glucose")) + TUMOR("Tumor Cell") + ONCO(("Oncoprotein
Level")) + NOTE["Represents
the antigen
expression of the
tumor cell.
Cells with higher
oncoprotein levels are
more proliferative but
also more likely to be
killed by CAR-T cells."] + + CART -->|"moved by
chemotaxis
towards"| IMMUNO + IMMUNO <-->|"secretes"| TUMOR + CART -->|"Induces apoptosis to"| TUMOR + OXY -->|"lack of it
induces
necrosis"| TUMOR + OXY <-->|"consumes"| TUMOR + CART -->|"consumes"| OXY + GLU <-->|"consumes"| TUMOR + GLU -->|"lack of it
induces
apoptosis"| TUMOR + CART -->|"consumes"| GLU + TUMOR -->|"expresses"| ONCO + ONCO -.- NOTE + + linkStyle 1,4,6 marker-end:none + classDef cell font-size:22px + class CART,TUMOR cell +``` + The project is organized into the following components: @@ -95,7 +129,7 @@ The project is organized into the following components: ### Configuration Files -- **[`params.json`](params.json)**: JSON-formatted parameter file for configuring simulation runs without recompilation. +- **[`params.json`](params.json)**: JSON-formatted parameter file for configuring simulation runs without recompilation. In [`params4SphericalModel.json`](params4SphericalModel.json) there is a configuration example for an spherical tumor and in [`params4CylindricalModel.json`](params4CylindricalModel.json) the configuration to simulate a cylindrical one. - **[`bdm.toml`](bdm.toml)**: BioDynaMo-specific configuration for visualization settings. @@ -105,6 +139,12 @@ The project is organized into the following components: - **[`CARTopiaX_Simulation_Analysis.ipynb`](CARTopiaX_Simulation_Analysis.ipynb)**: Jupyter notebook for post-processing simulation results, generating plots, and statistical analysis. +### Model Calibration (`abm_calibration/`) + +- **[`optimize.py`](abm_calibration/optimize.py)**: Bayesian optimization workflow to calibrate model parameters against target data. See [Model Calibration](#model-calibration-bayesian-optimization). + +- **[`requirements.txt`](abm_calibration/requirements.txt)**: Python dependencies for the calibration workflow. + --- @@ -342,17 +382,33 @@ These are basic parameters that are commonly changed when designing a treatment | Parameter | Default Value | Units | Description | |-----------|---------------|-------|-------------| -| `seed` | 1 | - | Seed for random number generation to ensure reproducibility | +| `seed` | 42 | - | Seed for random number generation to ensure reproducibility | | `output_performance_statistics` | false | - | Enable/disable performance statistics output | | `total_minutes_to_simulate` | 43200 | minutes | Total simulation time (default: 30 days) | -| `initial_tumor_radius` | 150 | μm | Initial radius of the spherical tumor | +| `tumor_shape` | `"sphere"` | - | Shape of the initial tumor: `"sphere"` or `"cylinder"` | +| `initial_spherical_tumor_radius` | 150 | μm | Initial radius of the spherical tumor (only used if `tumor_shape` is `"sphere"`) | +| `cylindrical_tumor_radius` | 3000 | μm | Initial radius of the cylindrical tumor (only used if `tumor_shape` is `"cylinder"`) | +| `cylindrical_tumor_height` | 100 | μm | Initial height of the cylindrical tumor (only used if `tumor_shape` is `"cylinder"`) | +| `initial_number_of_cylindrical_tumor_cells` | 2800 | cells | Initial number of tumor cells in the cylindrical tumor (only used if `tumor_shape` is `"cylinder"`) | | `bounded_space_length` | 1000 | μm | Length of the cubic simulation domain | +| `bound_space_toplogy` | `"torus"` | - | Topology of the domain boundaries for the agents: `"torus"`, `"open"` or `"closed"` (see BioDynaMo documentation) | | `treatment` | `{0: 3957, 8: 3957}` | day:cells | Map of treatment days to CAR-T cell counts. Key = day, Value = number of cells | ### Specific Parameters: These are all the other hyperparameters that are not adviced to be modified unless having a deeper understanding of the model and code. +#### Domain Boundary Parameters + +These parameters restrict the region where cells can be located, e.g. to model physical barriers in cylindrical tumors. Coordinates are expressed in the simulation coordinate system, where the domain extends from `-bounded_space_length/2` to `+bounded_space_length/2`. By default they do not impose any restriction. + +| Parameter | Default Value | Units | Description | +|-----------|---------------|-------|-------------| +| `bounded_space_min_allowed_z` | `-bounded_space_length/2` (-500) | μm | Minimum allowed z coordinate for cells | +| `bounded_space_max_allowed_z` | `bounded_space_length/2` (500) | μm | Maximum allowed z coordinate for cells | +| `bounded_space_max_allowed_radius` | `bounded_space_length` (1000) | μm | Maximum allowed distance of cells from the center of the tumor (from the center point for spherical tumors and from the axis for cylindrical ones) | + + #### Time Step Parameters | Parameter | Default Value | Units | Description | @@ -361,7 +417,15 @@ These are all the other hyperparameters that are not adviced to be modified unle | `dt_mechanics` | 0.1 | minutes | Time step for mechanical forces between cells | | `dt_cycle` | 6 | minutes | Time step for cell cycle progression | | `dt_step` | 0.1 | minutes | General simulation time step (same as `dt_mechanics`) | + +#### Output Parameters + +| Parameter | Default Value | Units | Description | +|-----------|---------------|-------|-------------| | `output_csv_interval` | 7200 | steps | Interval for writing simulation data to CSV (default: 12 hours) | +| `output_information_dependent_on_radius` | false | - | Outputs an additional CSV (`output/data_dependent_on_radius_tumor.csv`) with information aggregated by distance from the center of the tumor (from the center point for spherical tumors and from the axis for cylindrical ones) | +| `max_radius_analysis_csv_dependent_on_radius` | 3000 | μm | Maximum radius considered in the radius-dependent CSV | +| `num_radius_intervals` | 10 | - | Number of equal radius intervals, from 0 to `max_radius_analysis_csv_dependent_on_radius`, in which the radius-dependent information is aggregated | #### Apoptosis Parameters @@ -378,13 +442,42 @@ These are all the other hyperparameters that are not adviced to be modified unle | Parameter | Default Value | Units | Description | |-----------|---------------|-------|-------------| | `resolution_grid_substances` | 50 | voxels/axis | Number of voxels per axis for diffusion grids | +| `min_initial_z_substances` | `-bounded_space_length/2` (-500) | μm | Minimum height at which the substance grids are initialized to a value different from 0 | +| `max_initial_z_substances` | `bounded_space_length/2` (500) | μm | Maximum height at which the substance grids are initialized to a value different from 0 | + +Oxygen Parameters + +| Parameter | Default Value | Units | Description | +|-----------|---------------|-------|-------------| | `diffusion_coefficient_oxygen` | 100000 | μm²/min | Diffusion coefficient for oxygen | | `decay_constant_oxygen` | 0.1 | min⁻¹ | Decay constant (λ) for oxygen | -| `diffusion_coefficient_immunostimulatory_factor` | 1000 | μm²/min | Diffusion coefficient for immunostimulatory factor | -| `decay_constant_immunostimulatory_factor` | 0.016 | min⁻¹ | Decay constant (λ) for immunostimulatory factor | | `oxygen_reference_level` | 38 | mmHg | Boundary condition value for oxygen concentration | | `initial_oxygen_level` | 38 | mmHg | Initial oxygen concentration in all voxels | -| `oxygen_saturation` | 30 | mmHg | Oxygen saturation level in the microenvironment | +| `lateral_oxygen_production_min_z` | `-bounded_space_length/2` (-500) | μm | Minimum height at which the Dirichlet oxygen boundary condition is applied. If it equals the domain's minimum z, the floor also produces oxygen | +| `lateral_oxygen_production_max_z` | `bounded_space_length/2` (500) | μm | Maximum height at which the Dirichlet oxygen boundary condition is applied. If it equals the domain's maximum z, the roof also produces oxygen | +| `diffuse_oxygen_on_z_axis` | true | - | Whether oxygen also diffuses along the z axis. If false, it only diffuses in the x and y axes (ideal for a 2D model) | + +Immunostimulatory Factor Parameters + +| Parameter | Default Value | Units | Description | +|-----------|---------------|-------|-------------| +| `add_immunostimulatory_factor` | true | - | Whether to add the immunostimulatory factor diffusion grid | +| `diffusion_coefficient_immunostimulatory_factor` | 1000 | μm²/min | Diffusion coefficient for immunostimulatory factor | +| `decay_constant_immunostimulatory_factor` | 0.016 | min⁻¹ | Decay constant (λ) for immunostimulatory factor | +| `diffuse_immunostimulatory_factor_on_z_axis` | true | - | Whether the immunostimulatory factor also diffuses along the z axis. If false, it only diffuses in the x and y axes (ideal for a 2D model) | + +Glucose Parameters + +Glucose is disabled by default. When `add_glucose` is set to `true`, a glucose diffusion grid is added and glucose levels affect tumor cell growth and death (see the glucose parameters in [Tumor Cell Parameters](#tumor-cell-parameters)). + +| Parameter | Default Value | Units | Description | +|-----------|---------------|-------|-------------| +| `add_glucose` | false | - | Whether to add the glucose diffusion grid | +| `diffusion_coefficient_glucose` | 7800 | μm²/min | Diffusion coefficient for glucose | +| `decay_constant_glucose` | 0.01 | min⁻¹ | Decay constant (λ) for glucose | +| `initial_glucose_level` | 24.98 | mmol/L | Initial glucose concentration in each voxel | +| `max_radius_glucose_initialization` | `bounded_space_length` (1000) | μm | Radius of the sphere (spherical tumor) or cylinder (cylindrical tumor) inside which glucose is initialized. Outside of it glucose is initialized to 0 | +| `diffuse_glucose_on_z_axis` | true for `"sphere"`, false for `"cylinder"` | - | Whether glucose also diffuses along the z axis. If false, it only diffuses in the x and y axes (ideal for a 2D model) | #### Mechanical Forces Parameters @@ -402,32 +495,59 @@ These are all the other hyperparameters that are not adviced to be modified unle | `length_box_mechanics` | 22 | μm | Box length for spatial partitioning in force calculations | | `dnew` | 0.15 | - | Adams-Bashforth coefficient for current velocity (dt × 1.5) | | `dold` | -0.05 | - | Adams-Bashforth coefficient for previous velocity (dt × -0.5) | -| `max_squared_distance_cart_moving_towards_tumor_cell` | 317.746 | μm² | Maximum squared distance for CAR-T adhesion-only movement | #### Tumor Cell Parameters | Parameter | Default Value | Units | Description | |-----------|---------------|-------|-------------| -| `rate_secretion_immunostimulatory_factor` | 10 | 1/min | Secretion rate of immunostimulatory factor by tumor cells | -| `saturation_density_immunostimulatory_factor` | 1 | - | Saturation density for immunostimulatory factor secretion | | `oncoprotein_mean` | 1 | - | Mean oncoprotein expression level in tumor cells | | `oncoprotein_standard_deviation` | 0.25 | - | Standard deviation of oncoprotein expression | | `oncoprotein_limit` | 0.5 | - | Minimum oncoprotein level for CAR-T recognition | | `oncoprotein_saturation` | 2.0 | - | Maximum oncoprotein level | -| `oxygen_saturation_for_proliferation` | 38 | mmHg | Oxygen level for maximum proliferation rate | -| `oxygen_limit_for_proliferation` | 10 | mmHg | Minimum oxygen level for cell proliferation | -| `oxygen_limit_for_necrosis` | 5 | mmHg | Oxygen level below which necrosis begins | -| `oxygen_limit_for_necrosis_maximum` | 2.5 | mmHg | Oxygen level for maximum necrosis probability | | `time_lysis` | 86400 | minutes | Time until lysed necrotic cell removal | -| `maximum_necrosis_rate` | 0.00277778 | min⁻¹ | Maximum necrosis rate at 0 oxygen (1/360 min⁻¹) | -| `default_oxygen_consumption_tumor_cell` | 10 | 1/min | Baseline oxygen consumption rate of tumor cells | -| `default_volume_new_tumor_cell` | 2494 | μm³ | Total volume of newly created tumor cell | -| `default_volume_nucleus_tumor_cell` | 540 | μm³ | Nuclear volume of newly created tumor cell | +| `basal_death_probability_cancer_cells` | 0 | min⁻¹ | Basal death probability of tumor cells due to natural random causes | +| `default_volume_new_tumor_cell` | 2494 | μm³ | Mean total volume of newly created tumor cell | +| `std_volume_new_tumor_cell` | 0 | μm³ | Standard deviation of the total volume of newly created tumor cells | +| `min_volume_new_tumor_cell` | 2494 | μm³ | Minimum total volume of newly created tumor cells (sampled volumes are clipped) | +| `max_volume_new_tumor_cell` | 2494 | μm³ | Maximum total volume of newly created tumor cells (sampled volumes are clipped) | +| `default_fraction_of_volume_for_nucleus_tumor_cell` | 0.21652 | - | Nuclear volume of newly created tumor cell | | `default_fraction_fluid_tumor_cell` | 0.75 | - | Fraction of cytoplasmic volume that is fluid | | `average_time_transformation_random_rate` | 38.6 | hours | Mean cell cycle duration | | `standard_deviation_transformation_random_rate` | 3.7 | hours | Standard deviation of cell cycle duration | +| `minimum_tumor_cell_target_volume_fraction_for_division` | 0 | - | Minimum fraction of its target volume a tumor cell must reach to be able to divide | | `adhesion_time` | 60 | minutes | Average time tumor cell remains attached to CAR-T before escaping | +Tumor Cell Oxygen Parameters + +| Parameter | Default Value | Units | Description | +|-----------|---------------|-------|-------------| +| `default_oxygen_consumption_tumor_cell` | 10 | 1/min | Baseline oxygen consumption rate of tumor cells | +| `oxygen_saturation_for_proliferation` | 38 | mmHg | Oxygen level for maximum proliferation rate | +| `oxygen_limit_for_proliferation` | 10 | mmHg | Minimum oxygen level for cell proliferation | +| `oxygen_limit_for_necrosis` | 5 | mmHg | Oxygen level below which necrosis begins | +| `oxygen_limit_for_necrosis_maximum` | 2.5 | mmHg | Oxygen level for maximum necrosis probability | +| `maximum_necrosis_lack_of_oxygen_rate` | 0.00277778 | min⁻¹ | Maximum necrosis rate at 0 oxygen (1/360 min⁻¹) | + +Tumor Cell Immunostimulatory Factor Parameters + +| Parameter | Default Value | Units | Description | +|-----------|---------------|-------|-------------| +| `rate_secretion_immunostimulatory_factor` | 10 | 1/min | Secretion rate of immunostimulatory factor by tumor cells | +| `saturation_density_immunostimulatory_factor` | 1 | - | Saturation density for immunostimulatory factor secretion | + +Tumor Cell Glucose Parameters + +These parameters only have an effect when `add_glucose` is `true`. + +| Parameter | Default Value | Units | Description | +|-----------|---------------|-------|-------------| +| `default_glucose_consumption_tumor_cell` | 0.0007 | 1/min | Baseline glucose consumption rate of tumor cells | +| `glucose_saturation_for_tumor_cell_growth` | 0 | mmol/L | Glucose level above which tumor cells grow at full speed. Below it, growth slows down linearly, which also delays proliferation since cells must grow before dividing | +| `glucose_limit_for_tumor_cell_growth` | 0 | mmol/L | Glucose level below which tumor cells stop growing | +| `glucose_limit_for_death` | 0 | mmol/L | Glucose level below which tumor cells start dying from lack of glucose | +| `glucose_limit_for_death_maximum` | 0 | mmol/L | Glucose level at which the death probability due to lack of glucose is maximum | +| `maximum_death_lack_of_glucose_rate` | 0 | min⁻¹ | Maximum death rate due to lack of glucose | + Tumor Cell Volume Relaxation Rates | Parameter | Default Value | Units | Description | @@ -455,9 +575,15 @@ These are all the other hyperparameters that are not adviced to be modified unle | Parameter | Default Value | Units | Description | |-----------|---------------|-------|-------------| -| `average_maximum_time_untill_apoptosis_cart` | 12342.86 | minutes | Average CAR-T cell lifespan | +| `average_maximum_time_until_apoptosis_cart` | 12342.86 | minutes | Average CAR-T cell lifespan | | `default_oxygen_consumption_cart` | 1 | 1/min | Baseline oxygen consumption rate of CAR-T cells | -| `default_volume_new_cart_cell` | 2494 | μm³ | Total volume of newly created CAR-T cell | +| `default_glucose_consumption_cart` | 0.0007 | 1/min | Baseline glucose consumption rate of CAR-T cells (only used if `add_glucose` is `true`) | +| `default_volume_new_cart_cell` | 2494 | μm³ | Mean total volume of newly created CAR-T cell | +| `std_volume_new_cart_cell` | 0 | μm³ | Standard deviation of the total volume of newly created CAR-T cells | +| `min_volume_new_cart_cell` | 2494 | μm³ | Minimum total volume of newly created CAR-T cells (sampled volumes are clipped) | +| `max_volume_new_cart_cell` | 2494 | μm³ | Maximum total volume of newly created CAR-T cells (sampled volumes are clipped) | +| `default_fraction_of_volume_for_nucleus_cart_cell` | 0.21652 | - | Nuclear volume fraction of newly created CAR-T cell | +| `default_fraction_fluid_cart_cell` | 0.75 | - | Fraction of the CAR-T cell volume that is fluid | | `kill_rate_cart` | 0.06667 | min⁻¹ | Rate at which CAR-T attempts to kill attached tumor cell | | `adhesion_rate_cart` | 0.013 | min⁻¹ | Rate at which CAR-T attempts to attach to tumor cells | | `max_adhesion_distance_cart` | 18 | μm | Maximum distance for CAR-T to tumor cell attachment | @@ -469,7 +595,8 @@ These are all the other hyperparameters that are not adviced to be modified unle | Parameter | Default Value | Units | Description | |-----------|---------------|-------|-------------| | `persistence_time_cart` | 10 | minutes | Average time before CAR-T changes direction | -| `migration_bias_cart` | 0.5 | - | Chemotaxis bias toward immunostimulatory factor (0=random, 1=fully directed) | +| `avg_migration_bias_cart` | 0.5 | - | Mean chemotaxis bias toward immunostimulatory factor (0=random, 1=fully directed) | +| `std_migration_bias_cart` | 0 | - | Standard deviation of the migration bias among CAR-T cells (sampled values are clipped to ≤ 1) | | `migration_speed_cart` | 5 | μm/min | CAR-T cell migration speed | | `elastic_constant_cart` | 0.01 | - | Elastic constant for CAR-T cell motility | @@ -482,7 +609,7 @@ To modify parameters, edit the [`params.json`](params.json) file. Only include p "seed": 1, "output_performance_statistics": true, "total_minutes_to_simulate": 43200, - "initial_tumor_radius": 150.0, + "initial_spherical_tumor_radius": 150.0, "treatment": { "0": 3957, "8": 3957 @@ -523,6 +650,58 @@ paraview ./output/CARTopiaX/CARTopiaX.pvsm --- +## Model Calibration (Bayesian Optimization) + +Many parameters of the model cannot be measured directly in the laboratory. [`optimize.py`](abm_calibration/optimize.py) finds the values of selected parameters that make the simulation reproduce some target data (e.g. experimental measurements) as closely as possible. + +### How it works + +Each simulation is expensive, so trying every combination of parameters (grid search) is not feasible. Instead, the script uses **Bayesian optimization** through [Optuna](https://optuna.org/) and its **Tree-structured Parzen Estimator (TPE)** sampler: + +1. Each *trial* proposes a value for every parameter to be calibrated within its range. +2. The simulation is run with those values `NUMBER_MONTE_CARLO` times with different random seeds, and the error against the target data is averaged. This reduces the effect of the stochastic nature of the model. +3. TPE builds a probabilistic model from all previous trials, separating the parameter values that gave good results from the ones that gave bad results. The next trial is sampled where good values are more likely, balancing *exploration* of new regions and *exploitation* of promising ones. +4. After `NUMBER_OF_TRIALS` trials, the parameters with the lowest error are reported. + +This way good parameters are usually found with far fewer simulations than with a random or grid search. + +### Error modes + +The error is selected with the `MODE` variable: + +| Mode | Target file | Description | +|------|-------------|-------------| +| `total` | `target_data/final_data.csv` | Mean squared error of a metric of the whole tumor over time (e.g. `average_oxygen_all_cells`, `tumor_radius`) | +| `radius` | `target_data/data_dependent_on_radius_tumor.csv` | Mean squared error of the radial profile of a metric at a fixed minute. Requires `output_information_dependent_on_radius: true` | +| `custom` | - | User-defined error implemented in `compute_custom_error()` | + +Target CSVs must have the same format as the files written by the simulation in `output/`, and are compared on the common values of `total_minutes`. + +### Configuration + +The sections of the script marked with `Change this` are meant to be adapted to each experiment: + +- **Settings:** `EXPERIMENT_ID`, `MODE`, `SEED`, `NUMBER_OF_TRIALS`, `NUMBER_MONTE_CARLO` and `BIODYNAMO_DIR` (path to `thisbdm.sh`). +- **`run_ABM()`:** fixed simulation configuration written to `params.json` for each run, together with the parameters being calibrated. +- **`objective()`:** parameters to be calibrated and their search ranges, e.g. `trial.suggest_float("initial_oxygen_level", 30, 40)`. +- **Error functions:** metric to compare and, for `radius` mode, the minute of the radial profile. + +### Running the calibration + +Install the dependencies and run the script from the repository root: +```bash +pip install -r abm_calibration/requirements.txt +python3 abm_calibration/optimize.py +``` + +Results are stored in `abm_calibration/experiment_/`: +- `abm_optuna.db`: SQLite database with all trials. Running again with the same `EXPERIMENT_ID` resumes the study instead of starting from scratch. Setting `NUMBER_OF_TRIALS = 0` just prints the best result stored. +- `optuna_results.csv`: table with the parameters and error of every trial. + +During the calibration ParaView export is disabled, and the original `params.json` and `bdm.toml` are always restored at the end, even if a run fails. + +--- + ## Acknowledgments This project builds upon the BioDynaMo simulation framework. diff --git a/abm_calibration/optimize.py b/abm_calibration/optimize.py index 1d13237..d930d9a 100644 --- a/abm_calibration/optimize.py +++ b/abm_calibration/optimize.py @@ -29,14 +29,19 @@ import pandas as pd # Change this: File Parameters for the desired experiment -EXPERIMENT_ID = 1 +EXPERIMENT_ID = 0 +# "total": compares against target_data/final_data.csv +# "radius": compares against target_data/data_dependent_on_radius_tumor.csv +# "custom": implement your own error function in compute_custom_error() +MODE = "total" + SEED = 42 # You can set the number of trials to 0 to skip the optimization and just load the best result from the database NUMBER_OF_TRIALS = 8 # Number of Monte Carlo simulations to run for each trial. Use one for an aproximation of the error with a single montecarlo run -NUMBER_MONTE_CARLO = 1 +NUMBER_MONTE_CARLO = 3 # BioDynaMo directory to execute the comand source thisbdm.sh, you can change it to your own path -BIODYNAMO_DIR = " /home/usuario/Desktop/biodynamo/build/bin/thisbdm.sh" +BIODYNAMO_DIR = "/home/usuario/Desktop/biodynamo/build/bin/thisbdm.sh" # Other Hyperparameters for the optimization process # LOGGING @@ -47,10 +52,18 @@ # Directory for the parameters file params.json PARAMS_PATH = Path(__file__).resolve().parent.parent / "params.json" +# Directories with the simulation output and the target data +OUTPUT_DIR = Path(__file__).resolve().parent.parent / "output" +TARGET_DIR = Path(__file__).resolve().parent / "target_data" + # Directory for this experiment EXPERIMENT_DIR = Path("abm_calibration") / f"experiment_{EXPERIMENT_ID}" EXPERIMENT_DIR.mkdir(parents=True, exist_ok=True) +####################################################### +# Change this: Simulations to run the ABM with the given parameters +####################################################### + # Function to run the ABM with the given parameters def run_ABM(params, seed): @@ -61,6 +74,7 @@ def run_ABM(params, seed): ] # Change this configuration for the ABM run + # (for MODE = "radius" it must set output_information_dependent_on_radius: True) config = { "seed": seed, "output_performance_statistics": False, @@ -79,45 +93,59 @@ def run_ABM(params, seed): subprocess.run(["bash", "-c", f"source {BIODYNAMO_DIR} && bdm run"], check=True) -# Change This: Function to compute the error between the ABM simulation results and the experimental data -def compute_error(): - target = Path(__file__).resolve().parent / "target_data" / "final_data.csv" - sim_dir = Path(__file__).resolve().parent.parent / "output" / "final_data.csv" +####################################################### +# Change this: Error functions (examples) +####################################################### + +# Error using the data of the whole tumor as a function of time (final_data.csv) +def compute_error_total(): + metric = "average_oxygen_all_cells" # Change this: choose the metric to compare, e.g., "average_oxygen_all_cells", "tumor_radius", etc. + return mse_against_target("final_data.csv", [metric]) + + +# Error using the data as a function of the radius at a fixed time (data_dependent_on_radius_tumor.csv) +def compute_error_radius(): + csv_name = "data_dependent_on_radius_tumor.csv" + metric = "tumor_cells_type5_dead" # Change this: prefix of the metric columns, e.g., "tumor_cells_type5_dead" + minute = 4320 # Change this: fixed time (total_minutes) at which the radial profile is compared + + target = TARGET_DIR / csv_name if not target.exists(): - logging.error("Missing target CSV: %s", target) + logging.error("Missing CSV: %s", target) return float("inf") - if not sim_dir.exists(): - logging.error("Missing simulation CSV: %s", sim_dir) + header = pd.read_csv(target, nrows=0).columns + columns = [c for c in header if c.startswith(f"{metric}_radius_")] + if not columns: + logging.error("No columns for metric %s in %s", metric, target) + return float("inf") + return mse_against_target(csv_name, columns, minute=minute) + + +# Change this: write your own error. Example: dead cells at the center of a tumor +# of 3000 micrometers of radius split into 20 measurement rings at minute 4320 (72 hours) +def compute_custom_error(): + sim = OUTPUT_DIR / "data_dependent_on_radius_tumor.csv" + if not sim.exists(): + logging.error("Missing simulation CSV: %s", sim) return float("inf") - df_t = pd.read_csv(target, usecols=["total_minutes", "average_oxygen_cancer_cells"]) df_s = pd.read_csv( - sim_dir, usecols=["total_minutes", "average_oxygen_cancer_cells"] + sim, usecols=["total_minutes", "tumor_cells_type5_dead_radius_0_to_150"] ) - - # Merge the target and simulation dataframes on the "total_minutes" column to align the data for error computation - merged = pd.merge( - df_t, df_s, on="total_minutes", how="inner", suffixes=("_t", "_s") + row = df_s[df_s["total_minutes"] == 4320].iloc[0] + target_value_center = 21 + return abs( + float(row["tumor_cells_type5_dead_radius_0_to_150"]) - target_value_center ) - if merged.empty: - logging.error("No common minutes between target and simulation") - return float("inf") - # Compute the mean squared error (MSE) between the average oxygen levels in the target and simulation data - mse = ( - ( - merged["average_oxygen_cancer_cells_s"] - - merged["average_oxygen_cancer_cells_t"] - ) - ** 2 - ).mean() - return float(mse) +####################################################### +# Change this: Objective function for the Optuna optimization process +####################################################### -# Objective function for the Optuna optimization process def objective(trial): # Change this: Define the parameters to be optimized and their ranges params = { @@ -144,17 +172,20 @@ def objective(trial): return error -# Fix the random seed for reproducibility -np.random.seed(SEED) +####################################################### +# Auxiliary functions +####################################################### -# Create an Optuna study to optimize the parameters of the ABM -study = optuna.create_study( - study_name="abm_calibration", - storage=f"sqlite:///{EXPERIMENT_DIR / 'abm_optuna.db'}", - load_if_exists=True, - direction="minimize", - sampler=optuna.samplers.TPESampler(seed=SEED), -) + +# Function to compute the error between the ABM simulation results and the target data +def compute_error(): + if MODE == "total": + return compute_error_total() + if MODE == "radius": + return compute_error_radius() + if MODE == "custom": + return compute_custom_error() + raise ValueError(f"Unknown MODE: {MODE}") # Function to set the export value in the bdm.toml file @@ -171,7 +202,51 @@ def set_export_value(path: Path, value: bool): path.write_text("\n".join(new_lines)) +# MSE between simulation and target for the given csv file and columns, merged on total_minutes +# (if minute is given, only that time is compared) +def mse_against_target(csv_name, columns, minute=None): + target = TARGET_DIR / csv_name + sim = OUTPUT_DIR / csv_name + + for path in (target, sim): + if not path.exists(): + logging.error("Missing CSV: %s", path) + return float("inf") + + df_t = pd.read_csv(target, usecols=["total_minutes", *columns]) + df_s = pd.read_csv(sim, usecols=["total_minutes", *columns]) + + if minute is not None: + df_t = df_t[df_t["total_minutes"] == minute] + df_s = df_s[df_s["total_minutes"] == minute] + + merged = pd.merge( + df_t, df_s, on="total_minutes", how="inner", suffixes=("_t", "_s") + ) + if merged.empty: + logging.error("No common minutes between target and simulation") + return float("inf") + + squared = [(merged[f"{c}_s"] - merged[f"{c}_t"]) ** 2 for c in columns] + return float(pd.concat(squared, axis=1).stack().mean()) + + +####################################################### +# Main +####################################################### + if __name__ == "__main__": + # Fix the random seed for reproducibility + np.random.seed(SEED) + + # Create an Optuna study to optimize the parameters of the ABM + study = optuna.create_study( + study_name="abm_calibration", + storage=f"sqlite:///{EXPERIMENT_DIR / 'abm_optuna.db'}", + load_if_exists=True, + direction="minimize", + sampler=optuna.samplers.TPESampler(seed=SEED), + ) # Save the original content of the bdm.toml and params.json files to restore them later original_bdm = BDM_TOML_PATH.read_text() diff --git a/assets/images/model_outline.png b/assets/images/model_outline.png index ae7a43e..458c6d0 100644 Binary files a/assets/images/model_outline.png and b/assets/images/model_outline.png differ diff --git a/bdm.toml b/bdm.toml index d8a720b..a124c49 100644 --- a/bdm.toml +++ b/bdm.toml @@ -1,13 +1,13 @@ [visualization] export = true -interval = 7200 +interval = 1200 [[visualize_agent]] name = "TumorCell" additional_data_members = [ - {name = "diameter_", function = "diameter_"}, - {name = "volume_", function = "volume_"}, - {name = "cell_type", function = "GetTypeAsInt"} + "diameter_", + "volume_", + "cell_type_" ] [[visualize_agent]] @@ -15,7 +15,4 @@ name = "CarTCell" additional_data_members = ["diameter_","volume_"] [[visualize_diffusion]] -name = "oxygen" - -[[visualize_diffusion]] -name = "immunostimulatory_factor" \ No newline at end of file +name = "oxygen" \ No newline at end of file diff --git a/delete_me.py b/delete_me.py deleted file mode 100644 index 86c71fd..0000000 --- a/delete_me.py +++ /dev/null @@ -1,12 +0,0 @@ -import numpy -import pandas - -try: - import optuna - optuna_version = optuna.__version__ -except ImportError: - optuna_version = "NO INSTALADO" - -print("numpy:", numpy.__version__) -print("pandas:", pandas.__version__) -print("optuna:", optuna_version) \ No newline at end of file diff --git a/params.json b/params.json index 89c8178..562de49 100644 --- a/params.json +++ b/params.json @@ -7,4 +7,4 @@ "0": 3957, "8": 3957 } -} +} \ No newline at end of file diff --git a/params4CylindricalModel.json b/params4CylindricalModel.json new file mode 100644 index 0000000..33a4529 --- /dev/null +++ b/params4CylindricalModel.json @@ -0,0 +1,74 @@ +{ + "seed": 42, + "bound_space_toplogy": "closed", + "num_radius_intervals": 100, + "lateral_oxygen_production_min_z": -300.0, + "lateral_oxygen_production_max_z": 300.0, + "min_initial_z_substances": -300.0, + "max_initial_z_substances": 300.0, + "diffuse_oxygen_on_z_axis": false, + "diffuse_glucose_on_z_axis": false, + "output_performance_statistics": false, + "total_minutes_to_simulate": 4320, + "output_csv_interval": 600, + "bounded_space_length": 6500.0, + "tumor_shape": "cylinder", + "output_information_dependent_on_radius": true, + "cylindrical_tumor_radius": 3000.0, + "max_radius_analysis_csv_dependent_on_radius": 3000.0, + "cylindrical_tumor_height": 100.0, + "initial_number_of_cylindrical_tumor_cells": 28000, + "default_volume_new_tumor_cell": 1468.0, + "std_volume_new_tumor_cell": 166, + "min_volume_new_tumor_cell": 1136.0, + "max_volume_new_tumor_cell": 1800.0, + "default_volume_new_cart_cell": 269.0, + "std_volume_new_cart_cell": 33.625, + "min_volume_new_cart_cell": 201.75, + "max_volume_new_cart_cell": 336.25, + "oncoprotein_mean": 1.0, + "oncoprotein_standard_deviation": 0.0, + "initial_oxygen_level": 0.0, + "oxygen_reference_level": 165.0, + "default_oxygen_consumption_tumor_cell": 77.47, + "default_glucose_consumption_tumor_cell": 0.045, + "diffusion_coefficient_oxygen": 180000.0, + "decay_constant_oxygen": 0.006666666666666667, + "time_apoptosis": 6000.0, + "time_lysis": 6000.0, + "average_time_transformation_random_rate": 72, + "standard_deviation_transformation_random_rate": 15.0, + "decay_constant_glucose": 0.0003333333333333333, + "oxygen_saturation_for_proliferation": 20, + "oxygen_limit_for_proliferation": 10, + "oxygen_limit_for_necrosis_maximum": 0.0, + "oxygen_limit_for_necrosis": 7.1, + "maximum_necrosis_lack_of_oxygen_rate": 7.2214e-05, + "reduction_consumption_dead_cells": 0.0, + "basal_death_probability_cancer_cells": 5e-06, + "bounded_space_min_allowed_z": -50.0, + "bounded_space_max_allowed_z": 50.0, + "bounded_space_max_allowed_radius": 3250.0, + "minimum_distance_from_tumor_to_spawn_cart": 0.0, + "add_immunostimulatory_factor": true, + "add_glucose": true, + "diffusion_coefficient_glucose": 7800, + "initial_glucose_level": 24.98, + "treatment": { + "0": 560000 + }, + "default_oxygen_consumption_cart": 16, + "default_glucose_consumption_cart": 0.007, + "adhesion_rate_cart": 8.5e-06, + "glucose_saturation_for_tumor_cell_growth": 24.98, + "glucose_limit_for_tumor_cell_growth": 0, + "glucose_limit_for_death": 5, + "glucose_limit_for_death_maximum": 0, + "maximum_death_lack_of_glucose_rate": 0.000175, + "minimum_tumor_cell_target_volume_fraction_for_division": 0.9, + "avg_migration_bias_cart": 0, + "std_migration_bias_cart": 0.135, + "persistence_time_cart": 0, + "diffusion_coefficient_immunostimulatory_factor": 100000, + "decay_constant_immunostimulatory_factor": 1e-05 +} \ No newline at end of file diff --git a/params4SphericalModel.json b/params4SphericalModel.json new file mode 100644 index 0000000..562de49 --- /dev/null +++ b/params4SphericalModel.json @@ -0,0 +1,10 @@ +{ + "seed": 42, + "output_performance_statistics": true, + "total_minutes_to_simulate": 43200, + "initial_tumor_radius": 150.0, + "treatment": { + "0": 3957, + "8": 3957 + } +} \ No newline at end of file diff --git a/src/agents/cart_cell.cc b/src/agents/cart_cell.cc index 7d7b16b..cfec9f6 100644 --- a/src/agents/cart_cell.cc +++ b/src/agents/cart_cell.cc @@ -48,18 +48,47 @@ CarTCell::CarTCell(const Real3& position) { SetPosition(position); Simulation* sim = Simulation::GetActive(); const auto* sparams = sim->GetParam()->Get(); - SetVolume(sparams->default_volume_new_cart_cell); + // Set default volume + real_t total_volume = SamplePositiveGaussian( + sparams->default_volume_new_cart_cell, sparams->std_volume_new_cart_cell); + // Clip the value between the minimum and maximum allowed values + if (total_volume < sparams->min_volume_new_cart_cell) { + total_volume = sparams->min_volume_new_cart_cell; + } else if (total_volume > sparams->max_volume_new_cart_cell) { + total_volume = sparams->max_volume_new_cart_cell; + } + SetVolume(total_volume); + // Set default fluid fraction + SetFluidFraction(sparams->default_fraction_fluid_cart_cell); + // Set default nuclear volume + SetNuclearVolume(sparams->default_fraction_of_volume_for_nucleus_cart_cell * + total_volume); + + // Set Random Migration Bias + migration_bias_ = SamplePositiveGaussian(sparams->avg_migration_bias_cart, + sparams->std_migration_bias_cart); + // Clip the value <= 1 + if (migration_bias_ > 1.0) { + migration_bias_ = 1.0; + } + const ResourceManager& rm = *sim->GetResourceManager(); oxygen_dgrid_ = rm.GetDiffusionGrid("oxygen"); immunostimulatory_factor_dgrid_ = - rm.GetDiffusionGrid("immunostimulatory_factor"); + sparams->add_immunostimulatory_factor + ? rm.GetDiffusionGrid("immunostimulatory_factor") + : nullptr; + glucose_dgrid_ = + sparams->add_glucose ? rm.GetDiffusionGrid("glucose") : nullptr; SetCurrentLiveTime((sparams->dt_cycle + 1) * - sparams->average_maximum_time_untill_apoptosis_cart); + sparams->average_maximum_time_until_apoptosis_cart); // Add Consumption and Secretion // Set default oxygen consumption rate SetOxygenConsumptionRate(sparams->default_oxygen_consumption_cart); - // Compute constants for all ConsumptionSecretion of Oxygen + // Set default glucose consumption rate + SetGlucoseConsumptionRate(sparams->default_glucose_consumption_cart); + // Compute constants for all ConsumptionSecretion of Oxygen and glucose ComputeConstantsConsumptionSecretion(); } @@ -141,7 +170,7 @@ void CarTCell::ChangeVolumeExponentialRelaxationEquation( // if the volume has changed if (new_volume != current_total_volume) { SetVolume(new_volume); - // Update constants for all ConsumptionSecretion of Oxygen and + // Update constants for all ConsumptionSecretion of Oxygen, glucose and // Immunostimulatory Factors ComputeConstantsConsumptionSecretion(); } @@ -176,15 +205,24 @@ Real3 CarTCell::CalculateDisplacement(const InteractionForce* force, if (rng->Uniform(0.0, 1.0) < sparams->motility_probability_cart) { // random direction as unitary vector const Real3 random_direction = GenerateRandomDirection(); - Real3 direction_to_immunostimulatory_factor; - // returns normalized gradient towards the immunostimulatory factor source - immunostimulatory_factor_dgrid_->GetGradient( - current_position, &direction_to_immunostimulatory_factor, true); - // motility = bias * direction_to_immunostimulatory_factor + - // (1-bias)*random_direction - motility = - sparams->migration_bias_cart * direction_to_immunostimulatory_factor + - sparams->migration_one_minus_bias_cart * random_direction; + // initialize motility with random direction + motility = random_direction; + + if (IsImmunostimulatoryFactorDefined()) { + // if there is an immunostimulatory factor gradient set, the CAR-T cell + // will moves partially towards it, otherwise it will move randomly + Real3 direction_to_immunostimulatory_factor; + // returns normalized gradient towards the immunostimulatory factor + // source + immunostimulatory_factor_dgrid_->GetGradient( + current_position, &direction_to_immunostimulatory_factor, true); + // motility = bias * direction_to_immunostimulatory_factor + + // (1-bias)*random_direction + const real_t bias = GetMigrationBias(); + motility = bias * direction_to_immunostimulatory_factor + + (1 - bias) * random_direction; + } + const real_t motility_norm_squared = motility[0] * motility[0] + motility[1] * motility[1] + motility[2] * motility[2]; @@ -197,7 +235,6 @@ Real3 CarTCell::CalculateDisplacement(const InteractionForce* force, motility * sparams->migration_speed_cart; } } - //-------------------------------------------- // If cell is not apoptotic if (state_ == CarTCellState::kAlive) { @@ -247,8 +284,7 @@ Real3 CarTCell::CalculateDisplacement(const InteractionForce* force, // is the adhesion force to avoid CAR-T non-stop pushing tumor // cells. In case of being closer than // sparams->max_squared_distance_cart_moving_towards_tumor_cell - // there is a probability kProbabilityPushing for the CAR-T to keep - // pushing the tumor cell + // CAR-T to keep stop pushing the tumor cell if (sq_norm_displac > sparams->max_squared_distance_cart_moving_towards_tumor_cell) { translation_velocity_on_point_mass[0] += @@ -280,6 +316,60 @@ Real3 CarTCell::CalculateDisplacement(const InteractionForce* force, older_velocity_ = translation_velocity_on_point_mass; + // Clamp the movement if it surpasses the more restricted specified + // boundaries. + const double min_z = sparams->bounded_space_min_allowed_z; + const double max_z = sparams->bounded_space_max_allowed_z; + const double max_r_sq = sparams->bounded_space_max_allowed_radius_squared; + + const Real3 next_position = current_position + movement_at_next_step; + real_t radi_sq = 0.0; + switch (sparams->tumor_shape) { + case TumorShape::kCylinder: { + // Check the z coordinate + if (next_position[2] < min_z) { + movement_at_next_step[2] = min_z - current_position[2]; + } else if (next_position[2] > max_z) { + movement_at_next_step[2] = max_z - current_position[2]; + } + // Only consider x and y for cylindrical rumor, distance to the axis of + // the cylinder + radi_sq = next_position[0] * next_position[0] + + next_position[1] * next_position[1]; + if (radi_sq > max_r_sq) { + // Project the target position radially onto the allowed boundary + const double scale_factor = std::sqrt(max_r_sq / radi_sq); + movement_at_next_step[0] = + next_position[0] * scale_factor - current_position[0]; + movement_at_next_step[1] = + next_position[1] * scale_factor - current_position[1]; + } + break; + } + + case TumorShape::kSphere: { + // Consider all three dimensions for spherical radius + radi_sq = next_position[0] * next_position[0] + + next_position[1] * next_position[1] + + next_position[2] * next_position[2]; + if (radi_sq > max_r_sq) { + // Project the target position radially onto the allowed boundary + const double scale_factor = std::sqrt(max_r_sq / radi_sq); + movement_at_next_step[0] = + next_position[0] * scale_factor - current_position[0]; + movement_at_next_step[1] = + next_position[1] * scale_factor - current_position[1]; + movement_at_next_step[2] = + next_position[2] * scale_factor - current_position[2]; + } + break; + } + + default: + Log::Error("TumorCell::CalculateDisplacement", + "Unknown tumor shape, please use 'sphere' or 'cylinder'."); + break; + } // Displacement return movement_at_next_step; } @@ -373,7 +463,7 @@ bool CarTCell::TryToInduceApoptosis(bdm::AgentPointer attached_cell, // The CAR-T has succeeded to induce apoptosis on the Cancer Cell if (succeeded) { - attached_cell->StartApoptosis(); + attached_cell->StartApoptosis(TumorCellDeathCause::kCartKill); } return succeeded; @@ -387,10 +477,17 @@ real_t CarTCell::ConsumeSecreteSubstance(int substance_id, if (substance_id == oxygen_dgrid_->GetContinuumId()) { // consuming oxygen res = (old_concentration + constant1_oxygen_) / constant2_oxygen_; - } else if (substance_id == - immunostimulatory_factor_dgrid_->GetContinuumId()) { - // CAR-T do not change immunostimulatory factor levels + } else if (IsImmunostimulatoryFactorDefined() && + substance_id == + immunostimulatory_factor_dgrid_->GetContinuumId()) { + // There is a definde immunostimulatory factor grid and this is the one + // being updated CAR-T do not change immunostimulatory factor levels res = old_concentration; + } else if (IsGlucoseDefined() && + substance_id == glucose_dgrid_->GetContinuumId()) { + // There is a definde glucose grid and this is the one being updated + // CAR-T do not change glucose levels + res = (old_concentration + constant1_glucose_) / constant2_glucose_; } else { throw std::invalid_argument("Unknown substance id: " + std::to_string(substance_id)); @@ -413,7 +510,7 @@ void CarTCell::ComputeConstantsConsumptionSecretion() { const real_t volume = GetVolume(); const auto* sparams = Simulation::GetActive()->GetParam()->Get(); // compute the constants for the differential equation explicit solution: for - // oxygen and immunostimulatory factor + // oxygen and glucose // dt*(cell_volume/voxel_volume)*quantity_secretion*substance_saturation = dt // · (V_k / V_voxel) · S_k · ρ*_k) constant1_oxygen_ = 0.; @@ -423,6 +520,17 @@ void CarTCell::ComputeConstantsConsumptionSecretion() { constant2_oxygen_ = 1 + sparams->dt_substances * (volume / sparams->voxel_volume) * (oxygen_consumption_rate_); + // Compute constants for glucose consumption if glucose is defined + if (IsGlucoseDefined()) { + // dt*(cell_volume/voxel_volume)*quantity_secretion*substance_saturation = + // dt · (V_k / V_voxel) · S_k · ρ*_k) + constant1_glucose_ = 0.; + // 1 + dt*(cell_volume/voxel_volume)*(quantity_secretion + + // quantity_consumption ) = [1 + dt · (V_k / V_voxel) · (S_k + U_k)] + constant2_glucose_ = 1 + sparams->dt_substances * + (volume / sparams->voxel_volume) * + (glucose_consumption_rate_); + } } /// Main behavior executed at each simulation step @@ -450,15 +558,28 @@ void StateControlCart::Run(Agent* agent) { cell->SetState(CarTCellState::kApoptotic); // Reset timer_state, it should be 0 anyway cell->SetTimerState(0); - // Set target volume to 0 (the cell will shrink) - cell->SetTargetCytoplasmSolid(0.0); - cell->SetTargetNucleusSolid(0.0); + // Set target fuid volume to 0 (the cell shrinks loosing all its + // liquids) + const real_t current_total_volume = cell->GetVolume(); + const real_t fluid_fraction = cell->GetFluidFraction(); + const real_t nuclear_volume = cell->GetNuclearVolume(); + const real_t current_cytoplasm_solid = + (current_total_volume - nuclear_volume) * (1 - fluid_fraction); + const real_t current_nuclear_solid = + nuclear_volume * (1 - fluid_fraction); + cell->SetTargetCytoplasmSolid(current_cytoplasm_solid); + cell->SetTargetNucleusSolid(current_nuclear_solid); cell->SetTargetFractionFluid(0.0); - cell->SetTargetRelationCytoplasmNucleus(0.0); + cell->SetTargetRelationCytoplasmNucleus(current_cytoplasm_solid / + current_nuclear_solid); // Reduce oxygen consumption cell->SetOxygenConsumptionRate( cell->GetOxygenConsumptionRate() * sparams->reduction_consumption_dead_cells); + // Reduce glucose consumption + cell->SetGlucoseConsumptionRate( + cell->GetGlucoseConsumptionRate() * + sparams->reduction_consumption_dead_cells); // Update constants for all Consumption of oxygen cell->ComputeConstantsConsumptionSecretion(); // Detach from tumor cell if it was attached diff --git a/src/agents/cart_cell.h b/src/agents/cart_cell.h index f0accc8..efd36f1 100644 --- a/src/agents/cart_cell.h +++ b/src/agents/cart_cell.h @@ -116,6 +116,9 @@ class CarTCell : public Cell, public ISubstanceInteractor { } bool IsAttachedToTumorCell() const { return attached_to_tumor_cell_; } + real_t GetMigrationBias() const { return migration_bias_; } + void SetMigrationBias(real_t bias) { migration_bias_ = bias; } + Real3 GetOlderVelocity() const { return older_velocity_; } void SetOlderVelocity(const Real3& velocity) { older_velocity_ = velocity; } @@ -124,6 +127,11 @@ class CarTCell : public Cell, public ISubstanceInteractor { oxygen_consumption_rate_ = rate; } + real_t GetGlucoseConsumptionRate() const { return glucose_consumption_rate_; } + void SetGlucoseConsumptionRate(real_t rate) { + glucose_consumption_rate_ = rate; + } + real_t GetCurrentLiveTime() const { return current_live_time_; } void SetCurrentLiveTime(real_t time) { current_live_time_ = time; } @@ -145,6 +153,14 @@ class CarTCell : public Cell, public ISubstanceInteractor { DiffusionGrid* GetImmunostimulatoryFactorDiffusionGrid() const { return immunostimulatory_factor_dgrid_; } + /// Returns wether the Immunostimulatory factor is defined + bool IsImmunostimulatoryFactorDefined() const { + return immunostimulatory_factor_dgrid_ != nullptr; + } + /// Returns the diffusion grid for glucose + DiffusionGrid* GetGlucoseDiffusionGrid() const { return glucose_dgrid_; } + /// Returns wether the Immunostimulatory factor is defined + bool IsGlucoseDefined() const { return glucose_dgrid_ != nullptr; } /// Change volume using exponential relaxation equation /// @@ -227,6 +243,9 @@ class CarTCell : public Cell, public ISubstanceInteractor { /// Pointer to the immunostimulatory factor diffusion grid DiffusionGrid* immunostimulatory_factor_dgrid_ = nullptr; + /// Pointer to the glucose diffusion grid + DiffusionGrid* glucose_dgrid_ = nullptr; + /// Flag indicating if the cell is attached to a tumor cell bool attached_to_tumor_cell_ = false; @@ -252,14 +271,17 @@ class CarTCell : public Cell, public ISubstanceInteractor { /// Target relation between cytoplasm and nucleus volumes real_t target_relation_cytoplasm_nucleus_ = 0.0; + /// Migration bias towards the immunostimulatory factor source (range [0,1]) + real_t migration_bias_ = 0.0; + /// Velocity of the cell in the previous time step Real3 older_velocity_ = {0, 0, 0}; /// Rate of oxygen consumption by the cell real_t oxygen_consumption_rate_ = 0.0; - /// Rate of immunostimulatory factor secretion by the cell - real_t immunostimulatory_factor_secretion_rate_ = 0.0; + /// Rate of glucose consumption by the cell + real_t glucose_consumption_rate_ = 0.0; /// Constant 1 for oxygen consumption/secretion differential equation solution real_t constant1_oxygen_ = 0.0; @@ -267,6 +289,14 @@ class CarTCell : public Cell, public ISubstanceInteractor { /// Constant 2 for oxygen consumption/secretion differential equation solution real_t constant2_oxygen_ = 0.0; + /// Constant 1 for glucose consumption/secretion differential equation + /// solution + real_t constant1_glucose_ = 0.0; + + /// Constant 2 for glucose consumption/secretion differential equation + /// solution + real_t constant2_glucose_ = 0.0; + /// Pointer to the attached tumor cell bdm::AgentPointer attached_cell_ptr_; }; diff --git a/src/agents/tumor_cell.cc b/src/agents/tumor_cell.cc index 7651754..aaacdeb 100644 --- a/src/agents/tumor_cell.cc +++ b/src/agents/tumor_cell.cc @@ -36,6 +36,7 @@ #include "core/util/log.h" #include "core/util/random.h" #include +#include #include #include #include @@ -47,25 +48,38 @@ TumorCell::TumorCell(const Real3& position) { const auto* sparams = Simulation::GetActive()->GetParam()->Get(); // volumes // Set default volume - SetVolume(sparams->default_volume_new_tumor_cell); + real_t total_volume = + SamplePositiveGaussian(sparams->default_volume_new_tumor_cell, + sparams->std_volume_new_tumor_cell); + // Clip the value between the minimum and maximum allowed values + if (total_volume < sparams->min_volume_new_tumor_cell) { + total_volume = sparams->min_volume_new_tumor_cell; + } else if (total_volume > sparams->max_volume_new_tumor_cell) { + total_volume = sparams->max_volume_new_tumor_cell; + } + SetVolume(total_volume); // Set default fluid fraction SetFluidFraction(sparams->default_fraction_fluid_tumor_cell); // Set default nuclear volume - SetNuclearVolume(sparams->default_volume_nucleus_tumor_cell); + SetNuclearVolume(sparams->default_fraction_of_volume_for_nucleus_tumor_cell * + total_volume); // target volumes // Set target fraction of fluid SetTargetFractionFluid(sparams->default_fraction_fluid_tumor_cell); // Set target relation between cytoplasm and nucleus + const real_t default_volume_nucleus_tumor_cell = + sparams->default_fraction_of_volume_for_nucleus_tumor_cell * + sparams->default_volume_new_tumor_cell; SetTargetRelationCytoplasmNucleus( (sparams->default_volume_new_tumor_cell - - sparams->default_volume_nucleus_tumor_cell) / - (kEpsilon + sparams->default_volume_nucleus_tumor_cell)); + default_volume_nucleus_tumor_cell) / + (kEpsilon + default_volume_nucleus_tumor_cell)); // Set target nucleus solid volume - SetTargetNucleusSolid(sparams->default_volume_nucleus_tumor_cell * + SetTargetNucleusSolid(default_volume_nucleus_tumor_cell * (1 - sparams->default_fraction_fluid_tumor_cell)); // Set target cytoplasm solid volume SetTargetCytoplasmSolid((sparams->default_volume_new_tumor_cell - - sparams->default_volume_nucleus_tumor_cell) * + default_volume_nucleus_tumor_cell) * (1 - sparams->default_fraction_fluid_tumor_cell)); // Set initial oncoprotein level with a truncated normal distribution @@ -74,11 +88,18 @@ TumorCell::TumorCell(const Real3& position) { ResourceManager* rm = Simulation::GetActive()->GetResourceManager(); // Pointer to oxygen diffusion grid oxygen_dgrid_ = rm->GetDiffusionGrid("oxygen"); - // Pointer to immunostimulatory_factor diffusion grid + // Pointer to immunostimulatory_factor diffusion grid if it is added in the + // simulation immunostimulatory_factor_dgrid_ = - rm->GetDiffusionGrid("immunostimulatory_factor"); + sparams->add_immunostimulatory_factor + ? rm->GetDiffusionGrid("immunostimulatory_factor") + : nullptr; + glucose_dgrid_ = + sparams->add_glucose ? rm->GetDiffusionGrid("glucose") : nullptr; // Set state transition random rate SetTransformationRandomRate(); + // Set basal death probability + SetBasalDeathProbability(sparams->basal_death_probability_cancer_cells); // Add Consumption and Secretion // Set default oxygen consumption rate @@ -86,7 +107,8 @@ TumorCell::TumorCell(const Real3& position) { // Set default immunostimulatory factor secretion rate SetImmunostimulatoryFactorSecretionRate( sparams->rate_secretion_immunostimulatory_factor); - // Compute constants for all ConsumptionSecretion of Oxygen and + SetGlucoseConsumptionRate(sparams->default_glucose_consumption_tumor_cell); + // Compute constants for all ConsumptionSecretion of Oxygen, glucose and // Immunostimulatory Factor ComputeConstantsConsumptionSecretion(); } @@ -106,7 +128,13 @@ void TumorCell::Initialize(const NewAgentEvent& event) { // Pointer to the oxygen diffusion grid oxygen_dgrid_ = mother->oxygen_dgrid_; // Pointer to the immunostimulatory_factor diffusion grid - immunostimulatory_factor_dgrid_ = mother->immunostimulatory_factor_dgrid_; + immunostimulatory_factor_dgrid_ = + mother->IsImmunostimulatoryFactorDefined() + ? mother->immunostimulatory_factor_dgrid_ + : nullptr; + // Pointer to the glucose diffusion grid + glucose_dgrid_ = + mother->IsGlucoseDefined() ? mother->glucose_dgrid_ : nullptr; // inherit oncoprotein level from mother cell this->SetOncoproteinLevel(mother->oncoprotein_level_); // inherit oxygen consumption from mother cell @@ -114,6 +142,8 @@ void TumorCell::Initialize(const NewAgentEvent& event) { // inherit immunostimulatory factor secretion rate from mother cell this->SetImmunostimulatoryFactorSecretionRate( mother->GetImmunostimulatoryFactorSecretionRate()); + // inherit glucose consumption from mother cell + this->SetGlucoseConsumptionRate(mother->GetGlucoseConsumptionRate()); // Update the constants for all ConsumptionSecretion mother->ComputeConstantsConsumptionSecretion(); @@ -140,6 +170,8 @@ void TumorCell::Initialize(const NewAgentEvent& event) { // Set state transition random rate this->SetTransformationRandomRate(); + // Set basal death probability + this->SetBasalDeathProbability(mother->GetBasalDeathProbability()); // Initially not attached to a cart this->attached_to_cart_ = false; } @@ -152,19 +184,19 @@ void TumorCell::SetOncoproteinLevel(real_t level) { // cell type if (level >= sparams->threshold_cancer_cell_type1) { // between 1.5 and 2.0 - type_ = TumorCellType::kType1; + SetType(TumorCellType::kType1); } else if (level >= sparams->threshold_cancer_cell_type2 && level < sparams->threshold_cancer_cell_type1) { - type_ = TumorCellType::kType2; + SetType(TumorCellType::kType2); } else if (level >= sparams->threshold_cancer_cell_type3 && level < sparams->threshold_cancer_cell_type2) { - type_ = TumorCellType::kType3; + SetType(TumorCellType::kType3); } else if (level >= sparams->threshold_cancer_cell_type4 && level < sparams->threshold_cancer_cell_type3) { - type_ = TumorCellType::kType4; + SetType(TumorCellType::kType4); } else { // undefined type - type_ = TumorCellType::kType0; + SetType(TumorCellType::kType0); } } @@ -252,7 +284,7 @@ void TumorCell::ChangeVolumeExponentialRelaxationEquation( // if the volume has changed if (new_volume != current_total_volume) { SetVolume(new_volume); - // Update constants for all ConsumptionSecretion of Oxygen and + // Update constants for all ConsumptionSecretion of Oxygen, Glucose and // Immunostimulatory Factors ComputeConstantsConsumptionSecretion(); } @@ -303,6 +335,62 @@ Real3 TumorCell::CalculateDisplacement(const InteractionForce* force, older_velocity_ = translation_velocity_on_point_mass; + // Clamp the movement if it surpasses the more restricted specified + // boundaries. + Real3 current_position = GetPosition(); + const double min_z = sparams->bounded_space_min_allowed_z; + const double max_z = sparams->bounded_space_max_allowed_z; + const double max_r_sq = sparams->bounded_space_max_allowed_radius_squared; + + const Real3 next_position = current_position + movement_at_next_step; + real_t radi_sq = 0.0; + switch (sparams->tumor_shape) { + case TumorShape::kCylinder: { + // Check the z coordinate + if (next_position[2] < min_z) { + movement_at_next_step[2] = min_z - current_position[2]; + } else if (next_position[2] > max_z) { + movement_at_next_step[2] = max_z - current_position[2]; + } + // Only consider x and y for cylindrical rumor, distance to the axis of + // the cylinder + radi_sq = next_position[0] * next_position[0] + + next_position[1] * next_position[1]; + if (radi_sq > max_r_sq) { + // Project the target position radially onto the allowed boundary + const double scale_factor = std::sqrt(max_r_sq / radi_sq); + movement_at_next_step[0] = + next_position[0] * scale_factor - current_position[0]; + movement_at_next_step[1] = + next_position[1] * scale_factor - current_position[1]; + } + break; + } + + case TumorShape::kSphere: { + // Consider all three dimensions for spherical radius + radi_sq = next_position[0] * next_position[0] + + next_position[1] * next_position[1] + + next_position[2] * next_position[2]; + if (radi_sq > max_r_sq) { + // Project the target position radially onto the allowed boundary + const double scale_factor = std::sqrt(max_r_sq / radi_sq); + movement_at_next_step[0] = + next_position[0] * scale_factor - current_position[0]; + movement_at_next_step[1] = + next_position[1] * scale_factor - current_position[1]; + movement_at_next_step[2] = + next_position[2] * scale_factor - current_position[2]; + } + break; + } + + default: + Log::Error("TumorCell::CalculateDisplacement", + "Unknown tumor shape, please use 'sphere' or 'cylinder'."); + break; + } + // Displacement return movement_at_next_step; } @@ -315,11 +403,18 @@ real_t TumorCell::ConsumeSecreteSubstance(int substance_id, 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 + } else if (IsImmunostimulatoryFactorDefined() && + substance_id == + immunostimulatory_factor_dgrid_->GetContinuumId()) { + // There is a definde immunostimulatory factor grid and this is the one + // being updated secreting immunostimulatory factor res = (old_concentration + constant1_immunostimulatory_factor_) / constant2_immunostimulatory_factor_; + } else if (IsGlucoseDefined() && + substance_id == glucose_dgrid_->GetContinuumId()) { + // There is a defined glucose grid and this is the one being updated + // CAR-T do not change glucose levels + res = (old_concentration + constant1_glucose_) / constant2_glucose_; } else { throw std::invalid_argument("Unknown substance id: " + std::to_string(substance_id)); @@ -340,28 +435,45 @@ void TumorCell::ComputeConstantsConsumptionSecretion() { const auto* sparams = Simulation::GetActive()->GetParam()->Get(); const real_t new_volume = GetVolume(); // compute the constants for the differential equation explicit solution: for - // oxygen and immunostimulatory factor + // oxygen glucose 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_ * - sparams->saturation_density_immunostimulatory_factor * - sparams->dt_substances * (new_volume / sparams->voxel_volume); // 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 + sparams->dt_substances * (new_volume / sparams->voxel_volume) * (oxygen_consumption_rate_); - // Scale by the volume of the cell in the Voxel and time step - constant2_immunostimulatory_factor_ = - 1 + sparams->dt_substances * (new_volume / sparams->voxel_volume) * - (immunostimulatory_factor_secretion_rate_); + if (IsImmunostimulatoryFactorDefined()) { + // dt*(cell_volume/voxel_volume)*quantity_secretion*substance_saturation = + // dt · (V_k / V_voxel) · S_k · ρ*_k) Scale by the volume of the cell in the + // Voxel and time step + constant1_immunostimulatory_factor_ = + immunostimulatory_factor_secretion_rate_ * + sparams->saturation_density_immunostimulatory_factor * + sparams->dt_substances * (new_volume / sparams->voxel_volume); + // 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_immunostimulatory_factor_ = + 1 + sparams->dt_substances * (new_volume / sparams->voxel_volume) * + (immunostimulatory_factor_secretion_rate_); + } + + if (IsGlucoseDefined()) { + // dt*(cell_volume/voxel_volume)*quantity_secretion*substance_saturation = + // dt · (V_k / V_voxel) · S_k · ρ*_k) + constant1_glucose_ = 0.; + // 1 + dt*(cell_volume/voxel_volume)*(quantity_secretion + + // quantity_consumption ) = [1 + dt · (V_k / V_voxel) · (S_k + U_k)] + constant2_glucose_ = 1 + sparams->dt_substances * + (new_volume / sparams->voxel_volume) * + (glucose_consumption_rate_); + } } -void TumorCell::StartApoptosis() { +void TumorCell::StartApoptosis(TumorCellDeathCause cause) { // If the cell is already dead, do nothing if (IsDead()) { return; @@ -370,16 +482,25 @@ void TumorCell::StartApoptosis() { // The cell Dies SetState(TumorCellState::kApoptotic); + // Register the cause of death + SetDeathCause(cause); // Reset timer_state SetTimerState(0); // Set type to 5 to indicate dead cell SetType(TumorCellType::kType5); - // Set target volume to 0 (the cell shrinks) - SetTargetCytoplasmSolid(0.0); - SetTargetNucleusSolid(0.0); + // Set target fluid volume to 0 (the cell shrinks loosing all its liquids) + const real_t current_total_volume = GetVolume(); + const real_t fluid_fraction = GetFluidFraction(); + const real_t nuclear_volume = GetNuclearVolume(); + const real_t current_cytoplasm_solid = + (current_total_volume - nuclear_volume) * (1 - fluid_fraction); + const real_t current_nuclear_solid = nuclear_volume * (1 - fluid_fraction); + SetTargetCytoplasmSolid(current_cytoplasm_solid); + SetTargetNucleusSolid(current_nuclear_solid); SetTargetFractionFluid(0.0); - SetTargetRelationCytoplasmNucleus(0.0); + SetTargetRelationCytoplasmNucleus(current_cytoplasm_solid / + current_nuclear_solid); // Reduce oxygen consumption SetOxygenConsumptionRate(GetOxygenConsumptionRate() * sparams->reduction_consumption_dead_cells); @@ -406,12 +527,6 @@ void StateControlGrowProliferate::Run(Agent* agent) { // If the cell is attached to a cart, skip the state control and growth return; } - // Oxygen levels - const Real3 current_position = cell->GetPosition(); - // Pointer to the oxygen diffusion grid - DiffusionGrid* oxygen_dgrid = cell->GetOxygenDiffusionGrid(); - const real_t oxygen_level = oxygen_dgrid->GetValue(current_position); - switch (cell->GetState()) { case TumorCellState::kAlive: { // the cell is growing to real_t its size before mitosis @@ -419,12 +534,26 @@ void StateControlGrowProliferate::Run(Agent* agent) { // minutes per step) cell->SetTimerState(cell->GetTimerState() + sparams->dt_cycle); - // Enter necrosis if oxygen level is too low - if (ShouldEnterNecrosis(oxygen_level, cell)) { + // Oxygen levels + const Real3 current_position = cell->GetPosition(); + // Pointer to the oxygen diffusion grid + DiffusionGrid* oxygen_dgrid = cell->GetOxygenDiffusionGrid(); + const real_t oxygen_level = oxygen_dgrid->GetValue(current_position); + + // Glucose levels if a glucose diffusion grid is defined + real_t glucose_level = 0.0; + if (cell->IsGlucoseDefined()) { + DiffusionGrid* glucose_dgrid = cell->GetGlucoseDiffusionGrid(); + glucose_level = glucose_dgrid->GetValue(current_position); + } + + // Die if oxygen or glucose levels are too low or because of random + // natural causes + if (ShouldDie(oxygen_level, glucose_level, cell)) { // Exit the function to prevent further processing return; } - ManageLivingCell(cell, oxygen_level); + ManageLivingCell(cell, oxygen_level, glucose_level); break; } case TumorCellState::kNecroticSwelling: { @@ -447,11 +576,24 @@ void StateControlGrowProliferate::Run(Agent* agent) { cell->SetState(TumorCellState::kNecroticLysed); // Reset timer_state cell->SetTimerState(0); - // Set target volume to 0 (the cell will shrink) - cell->SetTargetCytoplasmSolid(0.0); - cell->SetTargetNucleusSolid(0.0); + // Set target fuid volume to 0 (the cell shrinks loosing all its + // liquids) + const real_t current_total_volume = cell->GetVolume(); + const real_t fluid_fraction = cell->GetFluidFraction(); + const real_t nuclear_volume = cell->GetNuclearVolume(); + const real_t current_cytoplasm_solid = + (current_total_volume - nuclear_volume) * (1 - fluid_fraction); + const real_t current_nuclear_solid = + nuclear_volume * (1 - fluid_fraction); + cell->SetTargetCytoplasmSolid(current_cytoplasm_solid); + cell->SetTargetNucleusSolid(current_nuclear_solid); cell->SetTargetFractionFluid(0.0); - cell->SetTargetRelationCytoplasmNucleus(0.0); + cell->SetTargetRelationCytoplasmNucleus(current_cytoplasm_solid / + current_nuclear_solid); + // Reduce oxygen consumption + cell->SetOxygenConsumptionRate( + cell->GetOxygenConsumptionRate() * + sparams->reduction_consumption_dead_cells); // Stop secretion and consumption rate // Stop consumption cell->SetOxygenConsumptionRate(0.0); @@ -516,30 +658,72 @@ void StateControlGrowProliferate::Run(Agent* agent) { // ManageLivingCell function to handle living cell behavior void StateControlGrowProliferate::ManageLivingCell(TumorCell* cell, - real_t oxygen_level) { - // Initialize multiplier - real_t multiplier = 1.0; + real_t oxygen_level, + real_t glucose_level) { const auto* sparams = Simulation::GetActive()->GetParam()->Get(); // volume change - cell->ChangeVolumeExponentialRelaxationEquation( - sparams->volume_relaxation_rate_alive_tumor_cell_cytoplasm, - sparams->volume_relaxation_rate_alive_tumor_cell_nucleus, - sparams->volume_relaxation_rate_alive_tumor_cell_fluid); + // If this tumoral cell is already bigger than the target volume, the volume + // does not change since it is already big enough to divide and it should not + // shrink + if (cell->GetVolume() < cell->GetTargetTotalVolume()) { + // If there is a glucose diffusion grid defined it affects the speed of + // growth of the cell Initialize multiplier + real_t multiplier_glucose = 1.0; + if (cell->IsGlucoseDefined()) { + // if there is a definde glucose gradient in the simulation + // glucose threshold for considering an effect on the growth speed + if (glucose_level < sparams->glucose_saturation_for_tumor_cell_growth) { + multiplier_glucose = + (glucose_level - sparams->glucose_limit_for_tumor_cell_growth) / + (sparams->glucose_saturation_for_tumor_cell_growth - + sparams->glucose_limit_for_tumor_cell_growth); + } + // If glucose is below the limit, set multiplier to 0 + if (glucose_level < sparams->glucose_limit_for_tumor_cell_growth) { + multiplier_glucose = 0.0; + } + } + // change volume change speed depending on the vailable glucose. This also + // affects ploriferation because if the cells do not grow big enough they + // cannot divide + cell->ChangeVolumeExponentialRelaxationEquation( + sparams->volume_relaxation_rate_alive_tumor_cell_cytoplasm * + multiplier_glucose, + sparams->volume_relaxation_rate_alive_tumor_cell_nucleus * + multiplier_glucose, + sparams->volume_relaxation_rate_alive_tumor_cell_fluid * + multiplier_glucose); + } + // cell state control + // The division rate depends on the available resources + // Oxygen + // Initialize multiplier + real_t multiplier_oxygen = 1.0; // oxygen threshold for considering an effect on the proliferation cycle if (oxygen_level < sparams->oxygen_saturation_for_proliferation) { - multiplier = (oxygen_level - sparams->oxygen_limit_for_proliferation) / - (sparams->oxygen_saturation_for_proliferation - - sparams->oxygen_limit_for_proliferation); + multiplier_oxygen = + (oxygen_level - sparams->oxygen_limit_for_proliferation) / + (sparams->oxygen_saturation_for_proliferation - + sparams->oxygen_limit_for_proliferation); } // If oxygen is below the limit, set multiplier to 0 if (oxygen_level < sparams->oxygen_limit_for_proliferation) { - multiplier = 0.0; + multiplier_oxygen = 0.0; } // Calculate the rate of state change based on oxygen level and oncoprotein // (min^-1) - const real_t final_rate_transition = cell->GetTransformationRandomRate() * - multiplier * cell->GetOncoproteinLevel(); + real_t final_rate_transition = cell->GetTransformationRandomRate() * + multiplier_oxygen * + cell->GetOncoproteinLevel(); + + // Check the dimensions of the cell, if it is too small, it cannot divide: it + // needs to grow more before dividing + if (cell->GetVolume() < + cell->GetTargetTotalVolume() * + sparams->minimum_tumor_cell_target_volume_fraction_for_division) { + final_rate_transition = 0; + } // Calculate the time to wait (in minutes) real_t time_to_wait = kTimeTooLarge; @@ -558,13 +742,16 @@ void StateControlGrowProliferate::ManageLivingCell(TumorCell* cell, } } -// computes the probability of the cell entering necrosis -bool StateControlGrowProliferate::ShouldEnterNecrosis(real_t oxygen_level, - TumorCell* cell) { +// computes the probability of the cell of dying due to low oxygen or glucose +// levels or due to random natural causes +bool StateControlGrowProliferate::ShouldDie(real_t oxygen_level, + real_t glucose_level, + TumorCell* cell) { Simulation* sim = Simulation::GetActive(); const auto* sparams = sim->GetParam()->Get(); - // necrosis probability - // Default multiplier for necrosis probability + Random* random = sim->GetRandom(); + // necrosis probability because of lack of Oxygen + // Default multiplier for necrosis probability because of low oxygen level real_t multiplier = 0.0; // oxygen threshold for considering necrosis if (oxygen_level < sparams->oxygen_limit_for_necrosis) { @@ -577,13 +764,15 @@ bool StateControlGrowProliferate::ShouldEnterNecrosis(real_t oxygen_level, multiplier = 1.0; } // Calculate the probability of necrosis based on oxygen level - // and multiply by sparams->dt_cycle since each timestep is sparams->dt_cycle - // minutes - const real_t probability_necrosis = - sparams->dt_cycle * sparams->maximum_necrosis_rate * multiplier; + const real_t maximum_necrosis_rate_oxygen_multiplier = + sparams->maximum_necrosis_lack_of_oxygen_rate * multiplier; - Random* random = sim->GetRandom(); + // Final probability: multiply by sparams->dt_cycle since each timestep is + // sparams->dt_cycle minutes + const real_t probability_necrosis = + sparams->dt_cycle * maximum_necrosis_rate_oxygen_multiplier; const bool enter_necrosis = random->Uniform(0, 1) < probability_necrosis; + // If the random number is less than the probability, enter necrosis if (enter_necrosis) { // If oxygen is too low, enter necrosis @@ -597,20 +786,79 @@ bool StateControlGrowProliferate::ShouldEnterNecrosis(real_t oxygen_level, // Reduce consumption cell->SetOxygenConsumptionRate(cell->GetOxygenConsumptionRate() * sparams->reduction_consumption_dead_cells); + cell->SetGlucoseConsumptionRate(cell->GetGlucoseConsumptionRate() * + sparams->reduction_consumption_dead_cells); // 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); - // Set target fraction of fluid to 1.0 + const real_t current_total_volume = cell->GetVolume(); + const real_t fluid_fraction = cell->GetFluidFraction(); + const real_t nuclear_volume = cell->GetNuclearVolume(); + const real_t current_cytoplasm_solid = + (current_total_volume - nuclear_volume) * (1 - fluid_fraction); + const real_t current_nuclear_solid = nuclear_volume * (1 - fluid_fraction); + cell->SetTargetCytoplasmSolid(current_cytoplasm_solid); + cell->SetTargetNucleusSolid(current_nuclear_solid); cell->SetTargetFractionFluid(1.0); - cell->SetTargetRelationCytoplasmNucleus(0.0); + cell->SetTargetRelationCytoplasmNucleus(current_cytoplasm_solid / + current_nuclear_solid); // Set type to 5 to indicate dead cell cell->SetType(TumorCellType::kType5); + // Register the cause of death + cell->SetDeathCause(TumorCellDeathCause::kLackOfOxygen); + return true; + } + + // If the cell does not eneter necrosis (due to low oxygen), check if it + // should die due to low glucose or random natural causes (apoptosis) Glucose + // (only if glucose is defined in the simulation) Default multiplier for + // apoptosis probability because of low glucose level + multiplier = 0.0; + if (cell->IsGlucoseDefined()) { + // there is a glucose diffusion grid defined in the simulation + // glucose threshold for considering necrosis + if (glucose_level < sparams->glucose_limit_for_death) { + multiplier = (sparams->glucose_limit_for_death - glucose_level) / + (sparams->glucose_limit_for_death - + sparams->glucose_limit_for_death_maximum); + } + // threshold for maximum apoptosis probability + if (glucose_level < sparams->glucose_limit_for_death_maximum) { + multiplier = 1.0; + } + } + // Calculate the probability of apoptosis based on glucose level. If no + // glucose gradient is defined it will be zero + const real_t maximum_death_rate_glucose_multiplier = + sparams->maximum_death_lack_of_glucose_rate * multiplier; + + // Lack of glucose: first independent check. Multiply by sparams->dt_cycle + // since each timestep is sparams->dt_cycle minutes. If there is no glucose + // effect the check is skipped + const real_t probability_glucose_death = + sparams->dt_cycle * maximum_death_rate_glucose_multiplier; + if (probability_glucose_death > 0 && + random->Uniform(0, 1) < probability_glucose_death) { + // If the random number is less than the probability, enter apoptosis + cell->StartApoptosis(TumorCellDeathCause::kLackOfGlucose); + return true; + } + + // Random natural causes: second independent check, only if the cell did not + // die because of lack of glucose. The total probability of apoptosis is + // dt*g + dt*b - dt^2*g*b + const real_t current_basal_death_probability = + cell->GetBasalDeathProbability(); + const real_t probability_natural_death = + sparams->dt_cycle * current_basal_death_probability; + if (random->Uniform(0, 1) < probability_natural_death) { + // If the random number is less than the probability, enter apoptosis + cell->StartApoptosis(TumorCellDeathCause::kRandomNaturalCauses); + return true; } - return enter_necrosis; // Return whether the cell entered necrosis + return false; // Return whether the cell died } } // namespace bdm diff --git a/src/agents/tumor_cell.h b/src/agents/tumor_cell.h index 982f9d2..b903ee5 100644 --- a/src/agents/tumor_cell.h +++ b/src/agents/tumor_cell.h @@ -70,6 +70,17 @@ enum class TumorCellType : int { kType5 = 5 ///< Dead tumor cells }; +/// Enumeration representing the cause of death of a tumor cell +/// +/// This enum class defines why a dead tumor cell (type 5) died. +enum class TumorCellDeathCause : int { + kNone = 0, ///< The cell is alive + kLackOfOxygen = 1, ///< Necrosis due to low oxygen levels + kLackOfGlucose = 2, ///< Apoptosis due to low glucose levels + kRandomNaturalCauses = 3, ///< Apoptosis due to random natural causes + kCartKill = 4 ///< Apoptosis induced by a CAR-T cell +}; + /// Tumor cell class implementation /// /// This class represents a cancer cell that forms a heterogeneous tumor in the @@ -149,11 +160,20 @@ class TumorCell : public Cell, public ISubstanceInteractor { return transformation_random_rate_; } + void SetBasalDeathProbability(real_t basal_death_probability) { + basal_death_probability_ = basal_death_probability; + } + real_t GetBasalDeathProbability() const { return basal_death_probability_; } + void SetAttachedToCart(bool attached) { attached_to_cart_ = attached; } bool IsAttachedToCart() const { return attached_to_cart_; } - void SetType(TumorCellType type) { type_ = type; } + void SetType(TumorCellType type) { + type_ = type; + cell_type_ = static_cast(type_); + } TumorCellType GetType() const { return type_; } + int GetTypeAsInt() const { return cell_type_; } Real3 GetOlderVelocity() const { return older_velocity_; } void SetOlderVelocity(const Real3& velocity) { older_velocity_ = velocity; } @@ -163,6 +183,11 @@ class TumorCell : public Cell, public ISubstanceInteractor { oxygen_consumption_rate_ = rate; } + real_t GetGlucoseConsumptionRate() const { return glucose_consumption_rate_; } + void SetGlucoseConsumptionRate(real_t rate) { + glucose_consumption_rate_ = rate; + } + real_t GetImmunostimulatoryFactorSecretionRate() const { return immunostimulatory_factor_secretion_rate_; } @@ -174,7 +199,8 @@ class TumorCell : public Cell, public ISubstanceInteractor { bool IsDead() const { return type_ == TumorCellType::kType5; } - int GetTypeAsInt() const { return static_cast(type_); } + void SetDeathCause(TumorCellDeathCause cause) { death_cause_ = cause; } + TumorCellDeathCause GetDeathCause() const { return death_cause_; } /// Returns the diffusion grid for oxygen DiffusionGrid* GetOxygenDiffusionGrid() const { return oxygen_dgrid_; } @@ -182,6 +208,14 @@ class TumorCell : public Cell, public ISubstanceInteractor { DiffusionGrid* GetImmunostimulatoryFactorDiffusionGrid() const { return immunostimulatory_factor_dgrid_; } + /// Returns wether the Immunostimulatory factor is defined + bool IsImmunostimulatoryFactorDefined() const { + return immunostimulatory_factor_dgrid_ != nullptr; + } + /// Returns the diffusion grid for glucose + DiffusionGrid* GetGlucoseDiffusionGrid() const { return glucose_dgrid_; } + /// Returns wether the Immunostimulatory factor is defined + bool IsGlucoseDefined() const { return glucose_dgrid_ != nullptr; } /// Change volume using exponential relaxation equation /// @@ -232,8 +266,10 @@ class TumorCell : public Cell, public ISubstanceInteractor { /// Start apoptosis /// /// This function is called when the tumor cell is induced to apoptosis by a - /// CAR-T cell. - void StartApoptosis(); + /// CAR-T cell, by lack of glucose or by random natural causes. + /// + /// @param cause The cause of death of the tumor cell + void StartApoptosis(TumorCellDeathCause cause); private: /// Current state of the tumor cell @@ -248,6 +284,9 @@ class TumorCell : public Cell, public ISubstanceInteractor { /// Pointer to the immunostimulatory factor diffusion grid DiffusionGrid* immunostimulatory_factor_dgrid_ = nullptr; + /// Pointer to the glucose diffusion grid + DiffusionGrid* glucose_dgrid_ = nullptr; + /// Level of oncoprotein expression real_t oncoprotein_level_ = 0.0; @@ -256,6 +295,9 @@ class TumorCell : public Cell, public ISubstanceInteractor { /// cell. This rate is kept constant during the cell's lifetime. real_t transformation_random_rate_ = 0.0; + /// Current probability of death due to natural causes + real_t basal_death_probability_ = 0.0; + /// Flag indicating if the cell is attached to a CAR-T cell bool attached_to_cart_ = false; @@ -285,6 +327,12 @@ class TumorCell : public Cell, public ISubstanceInteractor { /// aggressive. Type 5 means the cell is dead. TumorCellType type_ = TumorCellType::kType0; + /// Cell type as int for exporting to paraview + int cell_type_ = 0; + + /// Cause of death of the cell (kNone while the cell is alive) + TumorCellDeathCause death_cause_ = TumorCellDeathCause::kNone; + /// Velocity of the cell in the previous time step Real3 older_velocity_ = {0, 0, 0}; @@ -294,12 +342,23 @@ class TumorCell : public Cell, public ISubstanceInteractor { /// Rate of immunostimulatory factor secretion by the cell real_t immunostimulatory_factor_secretion_rate_ = 0.0; + /// Rate of glucose consumption by the cell + real_t glucose_consumption_rate_ = 0.0; + /// Constant 1 for oxygen consumption/secretion differential equation solution real_t constant1_oxygen_ = 0.0; /// Constant 2 for oxygen consumption/secretion differential equation solution real_t constant2_oxygen_ = 0.0; + /// Constant 1 for glucose consumption/secretion differential equation + /// solution + real_t constant1_glucose_ = 0.0; + + /// Constant 2 for glucose consumption/secretion differential equation + /// solution + real_t constant2_glucose_ = 0.0; + /// Constant 1 for immunostimulatory factor consumption/secretion differential /// equation solution real_t constant1_immunostimulatory_factor_ = 0.0; @@ -314,7 +373,7 @@ class TumorCell : public Cell, public ISubstanceInteractor { /// 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. +/// death based on oxygen levels and other environmental factors. struct StateControlGrowProliferate : public Behavior { // NOLINTNEXTLINE(llvm-else-after-return, moderinize-type-traits, @@ -337,20 +396,26 @@ struct StateControlGrowProliferate : public Behavior { void Run(Agent* agent) override; private: - /// Compute the probability of the cell entering necrosis + /// Compute the probability of the cell dying due to low oxygen or glucose + /// levels or due to random natural causes /// - /// Determines whether a cell should enter necrosis based on oxygen levels + /// Determines whether a cell should enter necrosis based on oxygen levels or + /// apoptosis because of random causes and, if deffined, glucose levels /// /// @param oxygen_level Current oxygen concentration at the cell's location + /// @param glucose_level Current glucose concentration at the cell's location /// @param cell Pointer to the tumor cell being evaluated - /// @return True if the cell should enter necrosis, false otherwise - static bool ShouldEnterNecrosis(real_t oxygen_level, TumorCell* cell); + /// @return True if the cell should die, false otherwise + static bool ShouldDie(real_t oxygen_level, real_t glucose_level, + TumorCell* cell); /// Manage the behavior of a living tumor cell /// /// @param cell Pointer to the tumor cell being managed /// @param oxygen_level Current oxygen concentration at the cell's location - static void ManageLivingCell(TumorCell* cell, real_t oxygen_level); + /// @param glucose_level Current glucose concentration at the cell's location + static void ManageLivingCell(TumorCell* cell, real_t oxygen_level, + real_t glucose_level); }; } // namespace bdm diff --git a/src/cart_tumor.cc b/src/cart_tumor.cc index 724ed08..f8048ae 100644 --- a/src/cart_tumor.cc +++ b/src/cart_tumor.cc @@ -21,6 +21,7 @@ #include "cart_tumor.h" #include "agents/tumor_cell.h" +#include "diffusion/cylinder_wall_boundary_condition.h" // NOLINT(llvm-include-order) #include "diffusion/diffusion_thomas_algorithm.h" #include "forces/forces_tumor_cart.h" #include "params/hyperparams.h" @@ -36,6 +37,7 @@ #include "core/resource_manager.h" #include "core/scheduler.h" #include "core/simulation.h" +#include "core/util/log.h" #include #include #include @@ -58,7 +60,7 @@ int Simulate(int argc, const char** argv) { auto set_param = [sparam_ref](Param* param) { // Set simulation bounds using the parameters param->random_seed = sparam_ref->seed; - param->bound_space = Param::BoundSpaceMode::kTorus; + param->bound_space = sparam_ref->bound_space_toplogy; param->min_bound = -sparam_ref->bounded_space_length / kHalf; param->max_bound = sparam_ref->bounded_space_length / kHalf; param->simulation_time_step = sparam_ref->dt_step; @@ -89,6 +91,9 @@ int Simulate(int argc, const char** argv) { // Define Substances ResourceManager* rm = Simulation::GetActive()->GetResourceManager(); + // Obtain the tumor shape from the parameters + const TumorShape tumor_shape = sparam->tumor_shape; + // Oxygen // substance_id, name, diffusion_coefficient, decay_constant, resolution, // time_step @@ -97,46 +102,140 @@ int Simulate(int argc, const char** argv) { kOxygen, "oxygen", sparam->diffusion_coefficient_oxygen, sparam->decay_constant_oxygen, sparam->resolution_grid_substances, sparam->dt_substances, - /*dirichlet_border=*/true); + /*dirichlet_border=*/true, sparam->diffuse_oxygen_on_z_axis); rm->AddContinuum(oxygen_grid.release()); - // Immunostimulatory Factor - // substance_id, name, diffusion_coefficient, decay_constant, resolution - std::unique_ptr immunostimulatory_factor_grid = - std::make_unique( - kImmunostimulatoryFactor, "immunostimulatory_factor", - sparam->diffusion_coefficient_immunostimulatory_factor, - sparam->decay_constant_immunostimulatory_factor, - sparam->resolution_grid_substances, sparam->dt_substances, - /*dirichlet_border=*/false); - rm->AddContinuum(immunostimulatory_factor_grid.release()); - // Boundary Conditions Dirichlet: simulating absorption or total loss at the // boundaries of the space. - // Oxygen comming from the borders (capillary vessels) + // Oxygen comming from the borders (capillary vessels). + // The boundary condition secretes oxygen wherever z falls within + // [lateral_oxygen_production_min_z, lateral_oxygen_production_max_z], + // regardless of which face of the domain is being evaluated (x, y, or z + // faces). Depending on how these two parameters are set, this results in + // either: + // - Oxygen produced only from a horizontal strip of the side walls, if + // lateral_oxygen_production_min_z, lateral_oxygen_production_max_z + // are set strictly inside the domain (floor and roof fall outside the + // range and produce no oxygen), or + // - Oxygen also produced from the floor and/or roof, if + // lateral_oxygen_production_min_z equals the domain's minimum z + // (floor) and/or lateral_oxygen_production_max_z equals the + // domain's maximum z (roof). ModelInitializer::AddBoundaryConditions( kOxygen, BoundaryConditionType::kDirichlet, // oxygen_reference_level mmHg is the physiological level of oxygen in // tissues, o2 saturation is 100% at this level - std::make_unique( - sparam->oxygen_reference_level)); - - // This is useless now but should be added this way in a future version of - // BioDynaMo - ModelInitializer::AddBoundaryConditions( - kImmunostimulatoryFactor, BoundaryConditionType::kNeumann, nullptr); + std::make_unique( + sparam->oxygen_reference_level, + sparam->lateral_oxygen_production_min_z, + sparam->lateral_oxygen_production_max_z)); + const real_t min_initial_z_substances = sparam->min_initial_z_substances; + const real_t max_initial_z_substances = sparam->max_initial_z_substances; + const real_t initial_oxygen_level = sparam->initial_oxygen_level; // Initialize oxygen voxels ModelInitializer::InitializeSubstance( - kOxygen, [sparam](real_t /*x*/, real_t /*y*/, real_t /*z*/) { - // Set all voxels to initial_oxygen_level mmHg - return sparam->initial_oxygen_level; + kOxygen, + [initial_oxygen_level, min_initial_z_substances, + max_initial_z_substances](real_t /*x*/, real_t /*y*/, real_t z) { + if (z <= max_initial_z_substances && z >= min_initial_z_substances) { + return initial_oxygen_level; + } + return 0.0; }); - // One spherical tumor of radius initial_tumor_radius in the center of the - // simulation space - const std::vector positions = - CreateSphereOfTumorCells(sparam->initial_tumor_radius); + // Immunostimulatory Factor + if (sparam->add_immunostimulatory_factor) { + // substance_id, name, diffusion_coefficient, decay_constant, resolution + std::unique_ptr immunostimulatory_factor_grid = + std::make_unique( + kImmunostimulatoryFactor, "immunostimulatory_factor", + sparam->diffusion_coefficient_immunostimulatory_factor, + sparam->decay_constant_immunostimulatory_factor, + sparam->resolution_grid_substances, sparam->dt_substances, + /*dirichlet_border=*/false, + sparam->diffuse_immunostimulatory_factor_on_z_axis); + rm->AddContinuum(immunostimulatory_factor_grid.release()); + + // This is useless now but should be added this way in a future version of + // BioDynaMo. Neumann boundary condition: the is not consumed nor produced + // at the boundaries + ModelInitializer::AddBoundaryConditions( + kImmunostimulatoryFactor, BoundaryConditionType::kNeumann, nullptr); + } + // Glucose + // substance_id, name, diffusion_coefficient, decay_constant, resolution, + // time_step + if (sparam->add_glucose) { + std::unique_ptr glucose_grid = + std::make_unique( + kGlucose, "glucose", sparam->diffusion_coefficient_glucose, + sparam->decay_constant_glucose, sparam->resolution_grid_substances, + sparam->dt_substances, + /*dirichlet_border=*/false, sparam->diffuse_glucose_on_z_axis); + rm->AddContinuum(glucose_grid.release()); + + // This is useless now but should be added this way in a future version of + // BioDynaMo. Neumann boundary condition: the is not consumed nor produced + // at the boundaries + ModelInitializer::AddBoundaryConditions( + kGlucose, BoundaryConditionType::kNeumann, nullptr); + const real_t initial_glucose_level = sparam->initial_glucose_level; + const real_t squared_radius_glucose_initialization = + sparam->max_radius_glucose_initialization * + sparam->max_radius_glucose_initialization; + // Initialize oxygen voxels + ModelInitializer::InitializeSubstance( + kGlucose, + [initial_glucose_level, min_initial_z_substances, + max_initial_z_substances, tumor_shape, + squared_radius_glucose_initialization](real_t x, real_t y, real_t z) { + if (z <= max_initial_z_substances && z >= min_initial_z_substances) { + if (tumor_shape == TumorShape::kSphere && + (x * x + y * y + z * z <= + squared_radius_glucose_initialization)) { + // Its spherical and the voxel is within the maximun spherical + // radius + return initial_glucose_level; + } + if (tumor_shape == TumorShape::kCylinder && + (x * x + y * y <= squared_radius_glucose_initialization)) { + // Its cylindrical and the voxel is within the maximun cylindrical + // radius + return initial_glucose_level; + } + } + return 0.0; + }); + } + + // Tumor cells initialization + std::vector positions; + switch (tumor_shape) { + case TumorShape::kSphere: { + // One spherical tumor of radius initial_spherical_tumor_radius in the + // center of the simulation space + positions = + CreateSphereOfTumorCells(sparam->initial_spherical_tumor_radius); + break; + } + + case TumorShape::kCylinder: { + // One cylindrical tumor of radius initial_spherical_tumor_radius and + // height cylindrical_tumor_height in the center of the simulation space + positions = CreateCylinderOfTumorCells( + sparam->cylindrical_tumor_radius, sparam->cylindrical_tumor_height, + sparam->initial_number_of_cylindrical_tumor_cells); + break; + } + + default: + Log::Error("Simulate", + "Unknown tumor shape, please use 'sphere' or 'cylinder'."); + // Exit with an error code + return 1; + } + for (const auto& pos : positions) { std::unique_ptr tumor_cell = std::make_unique(pos); std::unique_ptr state_control = diff --git a/src/cart_tumor.h b/src/cart_tumor.h index a8646fa..6c033a5 100644 --- a/src/cart_tumor.h +++ b/src/cart_tumor.h @@ -24,7 +24,7 @@ namespace bdm { /// List the diffused substances -enum Substances { kImmunostimulatoryFactor, kOxygen }; +enum Substances { kImmunostimulatoryFactor, kOxygen, kGlucose }; /// Function declaration for the main simulation int Simulate(int argc, const char** argv); diff --git a/src/diffusion/cylinder_wall_boundary_condition.cc b/src/diffusion/cylinder_wall_boundary_condition.cc new file mode 100644 index 0000000..54b2683 --- /dev/null +++ b/src/diffusion/cylinder_wall_boundary_condition.cc @@ -0,0 +1,40 @@ +/* + * Copyright 2026 compiler-research.org, Salvador de la Torre Gonzalez + * + * 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/cylinder_wall_boundary_condition.h" +#include "core/real_t.h" + +namespace bdm { + +CylinderWallBoundaryCondition::CylinderWallBoundaryCondition(real_t value, + real_t min_z, + real_t max_z) + : value_(value), min_z_(min_z), max_z_(max_z) {} + +real_t CylinderWallBoundaryCondition::Evaluate(real_t /*x*/, real_t /*y*/, + real_t z, + real_t /*time*/) const { + if (z < min_z_ || z > max_z_) { + return 0.0; + } + return value_; +} + +} // namespace bdm \ No newline at end of file diff --git a/src/diffusion/cylinder_wall_boundary_condition.h b/src/diffusion/cylinder_wall_boundary_condition.h new file mode 100644 index 0000000..fd8a95d --- /dev/null +++ b/src/diffusion/cylinder_wall_boundary_condition.h @@ -0,0 +1,76 @@ +/* + * Copyright 2026 compiler-research.org, Salvador de la Torre Gonzalez + * + * 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 CYLINDER_WALL_BOUNDARY_CONDITION_H_ // NOLINT(llvm-header-guard) +#define CYLINDER_WALL_BOUNDARY_CONDITION_H_ + +#include "core/diffusion/diffusion_grid.h" +#include "core/real_t.h" +#include "core/util/root.h" + +namespace bdm { + +/// @brief Dirichlet boundary condition that secretes a substance wherever +/// z falls within [min_z_, max_z_], regardless of which face of the domain +/// is being evaluated (x, y, or z faces). Outside that range, the boundary +/// condition evaluates to 0. +/// +/// In the typical use case — modeling secretion through the side walls of +/// a cylindrical capillary region — min_z_ and max_z_ are set strictly +/// inside the domain, so the top and bottom faces (z = domain min / domain +/// max) naturally fall outside [min_z_, max_z_] and evaluate to 0. +/// +/// However, this is a consequence of the chosen z-range, not a hardcoded +/// exclusion: if min_z_ is set exactly to the domain's minimum z (the +/// floor), the floor will also secrete value_ at that layer. Likewise, if +/// max_z_ is set exactly to the domain's maximum z (the roof), the roof +/// will also secrete value_ at that layer. Callers who need to guarantee +/// that floor/roof never secrete, regardless of the configured range, +/// should choose min_z_/max_z_ strictly inside the domain bounds. +class CylinderWallBoundaryCondition : public BoundaryCondition { + using BoundaryCondition::BoundaryCondition; + + public: + /// @param value Constant value secreted when min_z <= z <= max_z + /// @param min_z Lower bound (inclusive) of the z-range where value_ is + /// secreted. If set to the domain's minimum z, the floor also secretes. + /// @param max_z Upper bound (inclusive) of the z-range where value_ is + /// secreted. If set to the domain's maximum z, the roof also secretes. + CylinderWallBoundaryCondition(real_t value, real_t min_z, real_t max_z); + + /// @brief see BoundaryCondition::Evaluate() + [[nodiscard]] real_t Evaluate(real_t x, real_t y, real_t z, + real_t time) const final; + + private: + /// Constant value of the boundary condition within [min_z_, max_z_] + real_t value_ = 0.0; + /// Lower z-bound (inclusive) where the boundary condition is active + real_t min_z_ = 0.0; + /// Upper z-bound (inclusive) where the boundary condition is active + real_t max_z_ = 0.0; + + // NOLINTNEXTLINE(modernize-type-traits,llvm-else-after-return,readability-else-after-return,cppcoreguidelines-avoid-non-const-global-variables) + BDM_CLASS_DEF_OVERRIDE(CylinderWallBoundaryCondition, 1); +}; + +} // namespace bdm + +#endif // CYLINDER_WALL_BOUNDARY_CONDITION_H_ diff --git a/src/diffusion/diffusion_thomas_algorithm.cc b/src/diffusion/diffusion_thomas_algorithm.cc index 41b3b26..bb6e495 100644 --- a/src/diffusion/diffusion_thomas_algorithm.cc +++ b/src/diffusion/diffusion_thomas_algorithm.cc @@ -38,11 +38,9 @@ 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) +DiffusionThomasAlgorithm::DiffusionThomasAlgorithm( + int substance_id, std::string substance_name, real_t dc, real_t mu, + int resolution, real_t dt, bool dirichlet_border, bool diffuse_on_z_axis) : DiffusionGrid(substance_id, std::move(substance_name), dc, mu, resolution), resolution_(static_cast(GetResolution())), @@ -52,12 +50,14 @@ DiffusionThomasAlgorithm::DiffusionThomasAlgorithm(int substance_id, ->bounded_space_length) / static_cast(resolution_)), dirichlet_border_(dirichlet_border), + diffuse_on_z_axis_(diffuse_on_z_axis), jump_i_(1), jump_j_(resolution_), jump_(resolution_ * resolution_), spatial_diffusion_coeff_(dc * dt / (d_space_ * d_space_)), neg_diffusion_factor_(-spatial_diffusion_coeff_), - temporal_decay_coeff_(mu * dt / 3.0), + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + temporal_decay_coeff_(mu * dt / (diffuse_on_z_axis ? 3.0 : 2.0)), central_coeff_(1.0 + 2 * spatial_diffusion_coeff_ + temporal_decay_coeff_), edge_coeff_(1.0 + spatial_diffusion_coeff_ + temporal_decay_coeff_), @@ -208,9 +208,11 @@ void DiffusionThomasAlgorithm::DiffuseChemical() { SolveDirectionThomas(1); ApplyBoundaryConditionsIfNeeded(); - // Solve for Z-direction (direction = 2) - SolveDirectionThomas(2); - ApplyBoundaryConditionsIfNeeded(); + // Solve for Z-direction (direction = 2) if diffuse_on_z_axis is true + if (diffuse_on_z_axis_) { + SolveDirectionThomas(2); + ApplyBoundaryConditionsIfNeeded(); + } // Change of concentration levels because of agents ComputeConsumptionsSecretions(); diff --git a/src/diffusion/diffusion_thomas_algorithm.h b/src/diffusion/diffusion_thomas_algorithm.h index 261c304..c9c31fa 100644 --- a/src/diffusion/diffusion_thomas_algorithm.h +++ b/src/diffusion/diffusion_thomas_algorithm.h @@ -42,6 +42,7 @@ class DiffusionThomasAlgorithm : public DiffusionGrid { : resolution_(0), d_space_(0.0), dirichlet_border_(false), + diffuse_on_z_axis_(true), jump_i_(0), jump_j_(0), jump_(0), @@ -53,7 +54,7 @@ class DiffusionThomasAlgorithm : public DiffusionGrid { DiffusionThomasAlgorithm(int substance_id, std::string substance_name, real_t dc, real_t mu, int resolution, real_t dt, - bool dirichlet_border); + bool dirichlet_border, bool diffuse_on_z_axis); /// Concentration setters void SetConcentration(int x, int y, int z, real_t amount) { @@ -140,6 +141,9 @@ class DiffusionThomasAlgorithm : public DiffusionGrid { /// Flag indicating Dirichlet boundary conditions bool dirichlet_border_; + /// Flag indicating whether to diffuse along the z-axis + bool diffuse_on_z_axis_; + /// Index jump for i-direction (x-axis) int jump_i_; diff --git a/src/params/hyperparams.cc b/src/params/hyperparams.cc index 03568f1..252ac9d 100644 --- a/src/params/hyperparams.cc +++ b/src/params/hyperparams.cc @@ -20,8 +20,10 @@ */ #include "params/hyperparams.h" +#include "core/param/param.h" #include "core/param/param_group.h" #include "core/real_t.h" +#include "core/util/log.h" #include "core/util/math.h" #include #include @@ -69,12 +71,85 @@ void SimParam::LoadParams(const std::string& filename) { } }; + auto load_string = [&](const std::string& key, std::string& var) { + if (jfile.contains(key)) { + var = jfile[key].get(); + } + }; + // Load parameters from JSON file load_int("seed", seed); load_bool("output_performance_statistics", output_performance_statistics); load_int("total_minutes_to_simulate", total_minutes_to_simulate); - load_double("initial_tumor_radius", initial_tumor_radius); + + std::string tumor_shape_string; + load_string("tumor_shape", tumor_shape_string); + if (jfile.contains("tumor_shape")) { + tumor_shape_string = jfile["tumor_shape"].get(); + if (tumor_shape_string == "sphere") { + tumor_shape = TumorShape::kSphere; + } else if (tumor_shape_string == "cylinder") { + tumor_shape = TumorShape::kCylinder; + } else { + tumor_shape = TumorShape::kSphere; + Log::Error("SimParam::LoadParams", + "Unknown tumor shape '" + tumor_shape_string + + "'. It should be either 'sphere' or 'cylinder'. Loading " + "default value ('sphere') instead"); + } + } + + load_int("initial_number_of_cylindrical_tumor_cells", + initial_number_of_cylindrical_tumor_cells); load_int("bounded_space_length", bounded_space_length); + std::string bound_space_toplogy_string; + if (jfile.contains("bound_space_toplogy")) { + bound_space_toplogy_string = + jfile["bound_space_toplogy"].get(); + if (bound_space_toplogy_string == "torus") { + bound_space_toplogy = Param::BoundSpaceMode::kTorus; + } else if (bound_space_toplogy_string == "open") { + bound_space_toplogy = Param::BoundSpaceMode::kOpen; + } else if (bound_space_toplogy_string == "closed") { + bound_space_toplogy = Param::BoundSpaceMode::kClosed; + } else { + bound_space_toplogy = Param::BoundSpaceMode::kTorus; + Log::Error("SimParam::LoadParams", + "Unknown bound space topology '" + bound_space_toplogy_string + + "'. It should be either 'torus', 'open' or 'closed'. " + "Loading default value ('torus') instead"); + } + } + + load_double("initial_spherical_tumor_radius", initial_spherical_tumor_radius); + load_double("cylindrical_tumor_radius", cylindrical_tumor_radius); + load_double("cylindrical_tumor_height", cylindrical_tumor_height); + + // Half of the bounded space length (integer division, as in the original + // defaults) + const int half_bounded_space_length = bounded_space_length / 2; + + if (jfile.contains("bounded_space_min_allowed_z")) { + bounded_space_min_allowed_z = + jfile["bounded_space_min_allowed_z"].get(); + } else { + bounded_space_min_allowed_z = -half_bounded_space_length; + } + if (jfile.contains("bounded_space_max_allowed_z")) { + bounded_space_max_allowed_z = + jfile["bounded_space_max_allowed_z"].get(); + } else { + bounded_space_max_allowed_z = half_bounded_space_length; + } + + if (jfile.contains("bounded_space_max_allowed_radius")) { + bounded_space_max_allowed_radius = + jfile["bounded_space_max_allowed_radius"].get(); + } else { + // if it is not defined we assign the whole length of the bounded space to + // avoid any restrictions + bounded_space_max_allowed_radius = bounded_space_length; + } if (jfile.contains("treatment")) { treatment.clear(); @@ -103,6 +178,14 @@ void SimParam::LoadParams(const std::string& filename) { output_csv_interval = static_cast(12 * 60 / dt_step); } + load_bool("output_information_dependent_on_radius", + output_information_dependent_on_radius); + + load_double("max_radius_analysis_csv_dependent_on_radius", + max_radius_analysis_csv_dependent_on_radius); + + load_int("num_radius_intervals", num_radius_intervals); + load_double("volume_relaxation_rate_cytoplasm_apoptotic_cells", volume_relaxation_rate_cytoplasm_apoptotic_cells); load_double("volume_relaxation_rate_nucleus_apoptotic_cells", @@ -115,15 +198,64 @@ void SimParam::LoadParams(const std::string& filename) { reduction_consumption_dead_cells); load_int("resolution_grid_substances", resolution_grid_substances); + if (jfile.contains("min_initial_z_substances")) { + min_initial_z_substances = jfile["min_initial_z_substances"].get(); + } else { + min_initial_z_substances = -half_bounded_space_length; + } + + if (jfile.contains("max_initial_z_substances")) { + max_initial_z_substances = jfile["max_initial_z_substances"].get(); + } else { + max_initial_z_substances = half_bounded_space_length; + } + + if (jfile.contains("lateral_oxygen_production_min_z")) { + lateral_oxygen_production_min_z = + jfile["lateral_oxygen_production_min_z"].get(); + } else { + lateral_oxygen_production_min_z = -half_bounded_space_length; + } + + if (jfile.contains("lateral_oxygen_production_max_z")) { + lateral_oxygen_production_max_z = + jfile["lateral_oxygen_production_max_z"].get(); + } else { + lateral_oxygen_production_max_z = half_bounded_space_length; + } + load_bool("diffuse_oxygen_on_z_axis", diffuse_oxygen_on_z_axis); + load_bool("diffuse_immunostimulatory_factor_on_z_axis", + diffuse_immunostimulatory_factor_on_z_axis); + load_bool("diffuse_glucose_on_z_axis", diffuse_glucose_on_z_axis); + load_double("diffusion_coefficient_oxygen", diffusion_coefficient_oxygen); load_double("decay_constant_oxygen", decay_constant_oxygen); + load_double("oxygen_reference_level", oxygen_reference_level); + load_double("initial_oxygen_level", initial_oxygen_level); + load_bool("add_immunostimulatory_factor", add_immunostimulatory_factor); load_double("diffusion_coefficient_immunostimulatory_factor", diffusion_coefficient_immunostimulatory_factor); load_double("decay_constant_immunostimulatory_factor", decay_constant_immunostimulatory_factor); - load_double("oxygen_reference_level", oxygen_reference_level); - load_double("initial_oxygen_level", initial_oxygen_level); - load_double("oxygen_saturation", oxygen_saturation); + load_bool("add_glucose", add_glucose); + load_double("diffusion_coefficient_glucose", diffusion_coefficient_glucose); + load_double("decay_constant_glucose", decay_constant_glucose); + load_double("initial_glucose_level", initial_glucose_level); + + if (jfile.contains("max_radius_glucose_initialization")) { + max_radius_glucose_initialization = + jfile["max_radius_glucose_initialization"].get(); + } else { + max_radius_glucose_initialization = bounded_space_length; + } + + if (jfile.contains("diffuse_glucose_on_z_axis")) { + diffuse_glucose_on_z_axis = jfile["diffuse_glucose_on_z_axis"].get(); + } else { + // if the tumor shape is cylindrical it should be set to false, otherwise it + // should be set to true + diffuse_glucose_on_z_axis = (tumor_shape != TumorShape::kCylinder); + } load_double("cell_repulsion_between_tumor_tumor", cell_repulsion_between_tumor_tumor); @@ -168,23 +300,43 @@ void SimParam::LoadParams(const std::string& filename) { load_double("oxygen_saturation_for_proliferation", oxygen_saturation_for_proliferation); load_double("oxygen_limit_for_proliferation", oxygen_limit_for_proliferation); + load_double("glucose_saturation_for_tumor_cell_growth", + glucose_saturation_for_tumor_cell_growth); + load_double("glucose_limit_for_tumor_cell_growth", + glucose_limit_for_tumor_cell_growth); load_double("oxygen_limit_for_necrosis", oxygen_limit_for_necrosis); load_double("oxygen_limit_for_necrosis_maximum", oxygen_limit_for_necrosis_maximum); + load_double("maximum_necrosis_lack_of_oxygen_rate", + maximum_necrosis_lack_of_oxygen_rate); + load_double("glucose_limit_for_death", glucose_limit_for_death); + load_double("glucose_limit_for_death_maximum", + glucose_limit_for_death_maximum); + load_double("maximum_death_lack_of_glucose_rate", + maximum_death_lack_of_glucose_rate); + load_double("basal_death_probability_cancer_cells", + basal_death_probability_cancer_cells); load_double("time_lysis", time_lysis); - load_double("maximum_necrosis_rate", maximum_necrosis_rate); load_double("default_oxygen_consumption_tumor_cell", default_oxygen_consumption_tumor_cell); + load_double("default_glucose_consumption_tumor_cell", + default_glucose_consumption_tumor_cell); load_double("default_volume_new_tumor_cell", default_volume_new_tumor_cell); - load_double("default_volume_nucleus_tumor_cell", - default_volume_nucleus_tumor_cell); + load_double("std_volume_new_tumor_cell", std_volume_new_tumor_cell); + load_double("max_volume_new_tumor_cell", max_volume_new_tumor_cell); + load_double("min_volume_new_tumor_cell", min_volume_new_tumor_cell); + load_double("default_fraction_of_volume_for_nucleus_tumor_cell", + default_fraction_of_volume_for_nucleus_tumor_cell); load_double("default_fraction_fluid_tumor_cell", default_fraction_fluid_tumor_cell); load_double("average_time_transformation_random_rate", average_time_transformation_random_rate); load_double("standard_deviation_transformation_random_rate", standard_deviation_transformation_random_rate); + load_double("minimum_tumor_cell_target_volume_fraction_for_division", + minimum_tumor_cell_target_volume_fraction_for_division); + load_double("adhesion_time", adhesion_time); load_double("oncoprotein_limit", oncoprotein_limit); load_double("oncoprotein_saturation", oncoprotein_saturation); @@ -218,18 +370,27 @@ void SimParam::LoadParams(const std::string& filename) { load_double("threshold_cancer_cell_type3", threshold_cancer_cell_type3); load_double("threshold_cancer_cell_type4", threshold_cancer_cell_type4); - if (jfile.contains("average_maximum_time_untill_apoptosis_cart")) { - average_maximum_time_untill_apoptosis_cart = - jfile["average_maximum_time_untill_apoptosis_cart"].get(); + if (jfile.contains("average_maximum_time_until_apoptosis_cart")) { + average_maximum_time_until_apoptosis_cart = + jfile["average_maximum_time_until_apoptosis_cart"].get(); } else { - average_maximum_time_untill_apoptosis_cart = + average_maximum_time_until_apoptosis_cart = // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) dt_cycle * 10.0 * 24.0 * 60.0 / (dt_cycle + 1.0); } load_double("default_oxygen_consumption_cart", default_oxygen_consumption_cart); + load_double("default_glucose_consumption_cart", + default_glucose_consumption_cart); load_double("default_volume_new_cart_cell", default_volume_new_cart_cell); + load_double("std_volume_new_cart_cell", std_volume_new_cart_cell); + load_double("max_volume_new_cart_cell", max_volume_new_cart_cell); + load_double("min_volume_new_cart_cell", min_volume_new_cart_cell); + load_double("default_fraction_fluid_cart_cell", + default_fraction_fluid_cart_cell); + load_double("default_fraction_of_volume_for_nucleus_cart_cell", + default_fraction_of_volume_for_nucleus_cart_cell); load_double("kill_rate_cart", kill_rate_cart); load_double("adhesion_rate_cart", adhesion_rate_cart); load_double("max_adhesion_distance_cart", max_adhesion_distance_cart); @@ -237,7 +398,8 @@ void SimParam::LoadParams(const std::string& filename) { load_double("minimum_distance_from_tumor_to_spawn_cart", minimum_distance_from_tumor_to_spawn_cart); load_double("persistence_time_cart", persistence_time_cart); - load_double("migration_bias_cart", migration_bias_cart); + load_double("avg_migration_bias_cart", avg_migration_bias_cart); + load_double("std_migration_bias_cart", std_migration_bias_cart); load_double("migration_speed_cart", migration_speed_cart); load_double("elastic_constant_cart", elastic_constant_cart); @@ -249,16 +411,20 @@ void SimParam::LoadParams(const std::string& filename) { // Calculate steps per day. This is always calculated here // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) steps_in_one_day = static_cast(24 * 60 / dt_step); - // Calculate the volume of a single mechanical voxel in μm³ + // Calculate the volume of a single substance voxel in μm³ voxel_volume = (static_cast(bounded_space_length) / resolution_grid_substances) * (static_cast(bounded_space_length) / resolution_grid_substances) * (static_cast(bounded_space_length) / resolution_grid_substances); - // 1-migration_bias_cart - migration_one_minus_bias_cart = 1.0 - migration_bias_cart; + + // Squared max allowed radius + bounded_space_max_allowed_radius_squared = + bounded_space_max_allowed_radius * bounded_space_max_allowed_radius; + // Probability of a CAR-T cell to migrate in a given // mechanical time step - motility_probability_cart = dt_mechanics / persistence_time_cart; + motility_probability_cart = + (persistence_time_cart == 0) ? 1.0 : dt_mechanics / persistence_time_cart; // Probability of a Tumor cell to escape in a given // mechanical time step probability_escape_from_cart = @@ -312,10 +478,41 @@ void SimParam::PrintParams() const { << (output_performance_statistics ? "true" : "false") << "\n"; std::cout << "Total simulation time in minutes (30 days): " << total_minutes_to_simulate << "\n"; - std::cout << "Initial radius of the spherical tumor (micrometers): " - << initial_tumor_radius << "\n"; + switch (tumor_shape) { + case TumorShape::kSphere: + std::cout << "Tumor shape: 'sphere' \n"; + std::cout << "Initial radius of the spherical tumor (micrometers): " + << initial_spherical_tumor_radius << "\n"; + break; + case TumorShape::kCylinder: + std::cout << "Tumor shape: 'cylinder' \n"; + std::cout << "Initial radius of the cylindrical tumor (micrometers): " + << cylindrical_tumor_radius << "\n"; + std::cout << "Initial height of the cylindrical tumor (micrometers): " + << cylindrical_tumor_height << "\n"; + std::cout << "Min allowed Z coordinate for cells (micrometers): " + << bounded_space_min_allowed_z << "\n"; + std::cout << "Max allowed Z coordinate for cells (micrometers): " + << bounded_space_max_allowed_z << "\n"; + std::cout << "Max allowed radius for cells (micrometers): " + << bounded_space_max_allowed_radius << "\n"; + std::cout << "Initial number of cylindrical tumor cells: " + << initial_number_of_cylindrical_tumor_cells << "\n"; + break; + } std::cout << "Length of the bounded space (micrometers): " << bounded_space_length << "\n\n"; + switch (bound_space_toplogy) { + case Param::BoundSpaceMode::kTorus: + std::cout << "Bound space topology: 'torus' \n"; + break; + case Param::BoundSpaceMode::kOpen: + std::cout << "Bounded space topology: 'open' \n"; + break; + case Param::BoundSpaceMode::kClosed: + std::cout << "Bounded space topology: 'closed' \n"; + break; + } /// Treatment Dosages std::cout << "/// Treatment Dosages\n"; @@ -340,6 +537,17 @@ void SimParam::PrintParams() const { std::cout << "Time step for the cell cycle (minutes): " << dt_cycle << "\n"; std::cout << "General time step for the simulation: " << dt_step << "\n"; std::cout << "Output CSV interval: " << output_csv_interval << "\n\n"; + std::cout << "Output information dependent on radius: " + << (output_information_dependent_on_radius ? "true" : "false") + << "\n\n"; + if (output_information_dependent_on_radius) { + std::cout << "Number of radius intervals to output the agregated " + "information CSV: " + << num_radius_intervals << "\n\n"; + std::cout << "Maximum radius for the analysis of information dependent on " + "radius: " + << max_radius_analysis_csv_dependent_on_radius << "\n\n"; + } /// Apoptotic cells volume change std::cout << "/// Apoptotic cells volume change\n"; @@ -352,8 +560,7 @@ void SimParam::PrintParams() const { << volume_relaxation_rate_fluid_apoptotic_cells << "\n"; std::cout << "Time in minutes until an apoptotic cell is removed: " << time_apoptosis << "\n"; - std::cout << "Reduction of consumption rate of dead cells when they enter " - "necrosis: " + std::cout << "Reduction of consumption rate of dead cells when they die: " << reduction_consumption_dead_cells << "\n\n"; /// Chemicals @@ -362,20 +569,49 @@ void SimParam::PrintParams() const { std::cout << "///\n\n"; std::cout << "Number of voxels per axis for the substances grid: " << resolution_grid_substances << "\n"; + std::cout << "Minimum initial z-coordinate for substances values different " + "than 0 (micrometers): " + << min_initial_z_substances << "\n"; + std::cout << "Maximum initial z-coordinate for substances values different " + "than 0 (micrometers): " + << max_initial_z_substances << "\n"; + std::cout + << "Minimum z-coordinate for lateral oxygen production (micrometers): " + << lateral_oxygen_production_min_z << "\n"; + std::cout + << "Maximum z-coordinate for lateral oxygen production (micrometers): " + << lateral_oxygen_production_max_z << "\n"; + std::cout << "Whether to diffuse oxygen on the z-axis: " + << (diffuse_oxygen_on_z_axis ? "true" : "false") << "\n"; + std::cout << "Whether to diffuse immunostimulatory factor on the z-axis: " + << (diffuse_immunostimulatory_factor_on_z_axis ? "true" : "false") + << "\n"; + std::cout << "Whether to diffuse glucose on the z-axis: " + << (diffuse_glucose_on_z_axis ? "true" : "false") << "\n"; std::cout << "Diffusion coefficient of oxygen (μm²/min): " << diffusion_coefficient_oxygen << "\n"; std::cout << "Decay constant of oxygen (min⁻¹): " << decay_constant_oxygen << "\n"; - std::cout << "Diffusion coefficient of immunostimulatory factor (μm²/min): " - << diffusion_coefficient_immunostimulatory_factor << "\n"; - std::cout << "Decay constant of immunostimulatory factor (min⁻¹): " - << decay_constant_immunostimulatory_factor << "\n"; std::cout << "Reference level of oxygen at the boundaries (mmHg): " << oxygen_reference_level << "\n"; std::cout << "Initial oxygen concentration in each voxel (mmHg): " << initial_oxygen_level << "\n"; - std::cout << "Oxygen saturation in the microenvironment (mmHg): " - << oxygen_saturation << "\n\n"; + std::cout << "Whether to add immunostimulatory factor: " + << (add_immunostimulatory_factor ? "true" : "false") << "\n"; + std::cout << "Diffusion coefficient of immunostimulatory factor (μm²/min): " + << diffusion_coefficient_immunostimulatory_factor << "\n"; + std::cout << "Decay constant of immunostimulatory factor (min⁻¹): " + << decay_constant_immunostimulatory_factor << "\n"; + std::cout << "Whether to add glucose: " << (add_glucose ? "true" : "false") + << "\n"; + std::cout << "Diffusion coefficient of glucose (μm²/min): " + << diffusion_coefficient_glucose << "\n"; + std::cout << "Decay constant of glucose (min⁻¹): " << decay_constant_glucose + << "\n"; + std::cout << "Initial glucose concentration in each voxel (mM): " + << initial_glucose_level << "\n\n"; + std::cout << "Maximum radius for glucose initialization (micrometers): " + << max_radius_glucose_initialization << "\n\n"; /// Forces /// @@ -431,22 +667,44 @@ void SimParam::PrintParams() const { << oxygen_saturation_for_proliferation << "\n"; std::cout << "Limit of oxygen level for tumor cell proliferation: " << oxygen_limit_for_proliferation << "\n"; + std::cout << "Glucose saturation level in tumor cells for growth: " + << glucose_saturation_for_tumor_cell_growth << "\n"; + std::cout << "Limit of glucose level for tumor cell growth: " + << glucose_limit_for_tumor_cell_growth << "\n"; std::cout << "Limit of oxygen to start causing necrosis: " << oxygen_limit_for_necrosis << "\n"; std::cout << "Limit of oxygen to maximum necrosis probability: " << oxygen_limit_for_necrosis_maximum << "\n"; + std::cout << "Basal death probability for tumor cells: " + << basal_death_probability_cancer_cells << "\n"; std::cout << "Time in minutes until a lysed necrotic cell is removed: " << time_lysis << "\n"; std::cout << "Maximum rate per minute of necrosis for tumor cells: " - << maximum_necrosis_rate << "\n"; + << maximum_necrosis_lack_of_oxygen_rate << "\n"; + std::cout << "Limit of glucose to start causing death: " + << glucose_limit_for_death << "\n"; + std::cout << "Limit of glucose to maximum death probability: " + << glucose_limit_for_death_maximum << "\n"; + std::cout << "Maximum rate per minute of death for tumor cells: " + << maximum_death_lack_of_glucose_rate << "\n"; + std::cout << "Default oxygen consumption rate of tumor cell: " << default_oxygen_consumption_tumor_cell << "\n"; + std::cout << "Default glucose consumption rate of tumor cell: " + << default_glucose_consumption_tumor_cell << "\n"; std::cout << "\nVolume parameters:\n"; std::cout << "Default total volume of a new tumor cell (μm³): " << default_volume_new_tumor_cell << "\n"; - std::cout << "Default volume of the nucleus of a new tumor cell (μm³): " - << default_volume_nucleus_tumor_cell << "\n"; + std::cout + << "Standard deviation of the total volume of a new tumor cell (μm³): " + << std_volume_new_tumor_cell << "\n"; + std::cout << "Maximum total volume of a new tumor cell (μm³): " + << max_volume_new_tumor_cell << "\n"; + std::cout << "Minimum total volume of a new tumor cell (μm³): " + << min_volume_new_tumor_cell << "\n"; + std::cout << "Default fraction of volume for nucleus tumor cell: " + << default_fraction_of_volume_for_nucleus_tumor_cell << "\n"; std::cout << "Default fraction of fluid volume in a new tumor cell: " << default_fraction_fluid_tumor_cell << "\n"; @@ -455,6 +713,8 @@ void SimParam::PrintParams() const { << average_time_transformation_random_rate << "\n"; std::cout << "Standard Deviation for transformation Random Rate (hours): " << standard_deviation_transformation_random_rate << "\n"; + std::cout << "Minimum tumor cell target volume fraction for division: " + << minimum_tumor_cell_target_volume_fraction_for_division << "\n"; std::cout << "Average adhesion time for Tumor Cell under CAR-T attack (minutes): " << adhesion_time << "\n"; @@ -511,13 +771,26 @@ void SimParam::PrintParams() const { std::cout << "///\n\n"; std::cout << "Average time in minutes until a CAR-T cell dies: " - << average_maximum_time_untill_apoptosis_cart << "\n"; + << average_maximum_time_until_apoptosis_cart << "\n"; std::cout << "Default oxygen consumption rate of CAR-T cell: " << default_oxygen_consumption_cart << "\n"; + std::cout << "Default glucose consumption rate of CAR-T cell: " + << default_glucose_consumption_cart << "\n"; std::cout << "\nVolume parameters:\n"; std::cout << "Default total volume of a new CAR-T cell (μm³): " << default_volume_new_cart_cell << "\n"; + std::cout + << "Standard deviation of the total volume of a new CAR-T cell (μm³): " + << std_volume_new_cart_cell << "\n"; + std::cout << "Maximum total volume of a new CAR-T cell (μm³): " + << max_volume_new_cart_cell << "\n"; + std::cout << "Minimum total volume of a new CAR-T cell (μm³): " + << min_volume_new_cart_cell << "\n"; + std::cout << "Default fraction of fluid volume in a new CAR-T cell: " + << default_fraction_fluid_cart_cell << "\n"; + std::cout << "Default fraction of volume for nucleus CAR-T cell: " + << default_fraction_of_volume_for_nucleus_cart_cell << "\n"; std::cout << "\nKilling and adhesion rates:\n"; std::cout << "How often a CAR-T cell tries to kill an attached cancer cell " @@ -541,8 +814,11 @@ void SimParam::PrintParams() const { std::cout << "\nMotility parameters:\n"; std::cout << "Average persistence time before CAR-T cell moves: " << persistence_time_cart << "\n"; - std::cout << "Migration bias (higher values = more directed movement): " - << migration_bias_cart << "\n"; + std::cout + << "Average Migration bias (higher values = more directed movement): " + << avg_migration_bias_cart << "\n"; + std::cout << "Standard deviation migration bias: " << std_migration_bias_cart + << "\n"; std::cout << "Migration speed: " << migration_speed_cart << "\n"; std::cout << "Elastic constant: " << elastic_constant_cart << "\n\n"; diff --git a/src/params/hyperparams.h b/src/params/hyperparams.h index f8711b5..b8d4009 100644 --- a/src/params/hyperparams.h +++ b/src/params/hyperparams.h @@ -22,6 +22,7 @@ #ifndef TUMOR_HYPERPARAMS_H_ #define TUMOR_HYPERPARAMS_H_ +#include "core/param/param.h" #include "core/param/param_group.h" #include "core/real_t.h" #include "core/util/math.h" @@ -52,6 +53,9 @@ constexpr real_t kMinutesInAnHour = 60.0; /// Hours in a day constexpr real_t kHoursInADay = 24.0; +/// Enum for tumor shapes +enum class TumorShape { kSphere, kCylinder }; + /// Contains the default values of the hyperparameters used in the simulation. struct SimParam : public ParamGroup { // NOLINTNEXTLINE(modernize-type-traits,llvm-else-after-return,readability-else-after-return,cppcoreguidelines-owning-memory) @@ -69,21 +73,55 @@ struct SimParam : public ParamGroup { /// Seed for random number generation // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - int seed = 1; + int seed = 42; /// Output Performance Statistics // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) bool output_performance_statistics = false; /// Total simulation time in minutes (30 days by default) // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) int total_minutes_to_simulate = 43200; - /// Initial radius of the spherical tumor (group of cancer cells) in - /// micrometers + /// Bounded space topology for the boundaries it can be "open", "torus" or + /// "closed" (see BioDynaMo documentation for more information) // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t initial_tumor_radius = 150; + Param::BoundSpaceMode bound_space_toplogy = Param::BoundSpaceMode::kTorus; /// Length of the bounded space in micrometers // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) int bounded_space_length = 1000; + /// Tumor topology + /// Tumor shape read as "cylinder" or "sphere" from the JSON file + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + TumorShape tumor_shape = TumorShape::kSphere; + /// Initial radius of the spherical tumor (group of cancer cells) in + /// micrometers + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t initial_spherical_tumor_radius = 150; + /// Initial radius of the cylindrical tumor in micrometers + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t cylindrical_tumor_radius = 3000; + /// Initial height of the cylindrical tumor in micrometers + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t cylindrical_tumor_height = 100; + /// Min allowed Z coordinate for cells in cylindrical tumors with phisical + /// barriers. By default it is automatically set to -bounded_space_length/2 to + /// avoid any restrictions. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t bounded_space_min_allowed_z = -500; + /// Max allowed Z coordinate for cells in cylindrical tumors with phisical + /// barriers. By default it is automatically set to bounded_space_length/2 to + /// avoid any restrictions. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t bounded_space_max_allowed_z = 500; + /// Max allowed radius for cells meassured from the center of the tumor and + /// considering the z-coordinate depending wether it has a cylindrical or + /// spherical shape. By default it is automatically set to + /// bounded_space_length to avoid any restrictions. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t bounded_space_max_allowed_radius = 1000; + /// Initial number of cylindrical tumor cells + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + int initial_number_of_cylindrical_tumor_cells = 2800; + /// Treatment Dosages /// /// Specifies the CAR-T cell infusion schedule as a map where: @@ -114,6 +152,22 @@ struct SimParam : public ParamGroup { // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) int output_csv_interval = 7200; + /// Outputs an additional CSV with information dependent on the radius of the + /// tumor if the tumor is spherical its from its center point and if it is + /// cylindrical its from the axis of the cylinder + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + bool output_information_dependent_on_radius = false; + /// Maximum radius for the analysis of information dependent on radius + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t max_radius_analysis_csv_dependent_on_radius = 3000; + /// Number of radius intervals to output the agregated information CSV. The + /// concatenation of the intervals goes from 0 (the center of the tumor) to + /// max_radius_analysis_csv_dependent_on_radius. The number of intervals is + /// used to divide the radius into equal parts and output the agregated + /// information for each interval. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + int num_radius_intervals = 10; + /// Apoptotic cells volume change // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t volume_relaxation_rate_cytoplasm_apoptotic_cells = 0.0166667; @@ -124,7 +178,7 @@ struct SimParam : public ParamGroup { /// Time in minutes until an apoptotic cell is removed from the simulation // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t time_apoptosis = 516; - /// Reduction of consumption rate of dead cells when they enter necrosis + /// Reduction of consumption rate of dead cells when they die // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t reduction_consumption_dead_cells = 0.1; @@ -132,27 +186,101 @@ struct SimParam : public ParamGroup { /// Number of voxels per axis for the substances grid // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) int resolution_grid_substances = 50; + + /// Minimum height (μm) at which the substances grid is initialized to a value + /// different from 0. Heights are expressed in the simulation coordinate + /// system, where the domain extends from -bounded_space_length / 2 to + /// +bounded_space_length / 2 By default it is automatically set to + /// -bounded_space_length/2 to avoid any restrictions. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t min_initial_z_substances = -500; + /// Maximum height (μm) at which the substances grid is initialized to a value + /// different from 0. Heights are expressed in the simulation coordinate + /// system, where the domain extends from -bounded_space_length / 2 to + /// +bounded_space_length / 2 By default it is automatically set to + /// +bounded_space_length/2 to avoid any restrictions. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t max_initial_z_substances = 500; + + /// Minimum height (μm) at which lateral Dirichlet oxygen boundary conditions + /// are applied. Heights are expressed in the simulation coordinate system, + /// where the domain extends from -bounded_space_length / 2 to + /// +bounded_space_length / 2 along the z-axis. Setting this value equal to + /// the minimum domain height + /// (-bounded_space_length / 2) causes the floor boundary (z = + /// -bounded_space_length / 2) to also act as an oxygen-producing Dirichlet + /// boundary. By default it is automatically set to -bounded_space_length/2 + /// to avoid any restrictions. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t lateral_oxygen_production_min_z = -500.0; + /// Maximum height (μm) at which lateral Dirichlet oxygen boundary conditions + /// are applied. Heights are expressed in the simulation coordinate system, + /// where the domain extends from -bounded_space_length / 2 to + /// +bounded_space_length / 2 along the z-axis. Setting this value equal to + /// the maximum domain height + /// (+bounded_space_length / 2) causes the roof boundary (z = + /// +bounded_space_length / 2) to also act as an oxygen-producing Dirichlet + /// boundary. By default it is automatically set to bounded_space_length/2 to + /// avoid any restrictions. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t lateral_oxygen_production_max_z = 500.0; + /// Wether the Thomas algortithm should also diffuse the oxygen Chemical in + /// the z-axis or not. If false, the diffusion will only be done in the x and + /// y axis, ideal for a 2D model + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + bool diffuse_oxygen_on_z_axis = true; + /// Wether the Thomas algortithm should also diffuse the + /// immunostimulatory_factor Chemical in the z-axis or not. If false, the + /// diffusion will only be done in the x and y axis, ideal for a 2D model + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + bool diffuse_immunostimulatory_factor_on_z_axis = true; + /// Wether the Thomas algortithm should also diffuse the glucose Chemical in + /// the z-axis or not. By default it will be true if the tumor shape is + /// spherical and false if the tumor shape is cylindrical. If false, the + /// diffusion will only be done in the x and y axis, ideal for a 2D model + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + bool diffuse_glucose_on_z_axis = true; /// Diffusion coefficient of oxygen in μm²/min // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t diffusion_coefficient_oxygen = 100000; /// Decay constant of oxygen in min⁻¹ // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t decay_constant_oxygen = 0.1; + /// Reference level of oxygen at the boundaries in mmHg + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t oxygen_reference_level = 38; + /// Initial oxygen concentration in each voxel in mmHg + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t initial_oxygen_level = 38; + /// Wether to add immunostimulatory factor grid or not + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + bool add_immunostimulatory_factor = true; /// Diffusion coefficient of immunostimulatory factor in μm²/min // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t diffusion_coefficient_immunostimulatory_factor = 1000; /// Decay constant of immunostimulatory factor in min⁻¹ // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t decay_constant_immunostimulatory_factor = 0.016; - /// Reference level of oxygen at the boundaries in mmHg + /// Wether to add glucose gradient grid or not // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t oxygen_reference_level = 38; - /// Initial oxygen concentration in each voxel in mmHg + bool add_glucose = false; + /// Diffusion coefficient of immunostimulatory factor in μm²/min // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t initial_oxygen_level = 38; - /// Oxygen saturation in the microenvironment in mmHg + real_t diffusion_coefficient_glucose = 7800; + /// Decay constant of immunostimulatory factor in min⁻¹ + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t decay_constant_glucose = 0.01; + /// Initial glucose concentration in each voxel in mmol/L // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t oxygen_saturation = 30; + real_t initial_glucose_level = 24.98; + /// Maximum radius for glucose initialization in micrometers. If the tumor is + /// spherical, the glucose will be initialized in a sphere of this radius. If + /// the tumor is cylindrical, the glucose will be initialized in a cylinder of + /// this radius. Outside of this radius, the glucose will be initialized to 0. + /// By default it is set to the bounded_space_length to avoid any + /// restrictions. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t max_radius_glucose_initialization = 1000; /// Forces /// Repulsion coeficient between tumor cells @@ -223,6 +351,14 @@ struct SimParam : public ParamGroup { /// Oxygen saturation level in tumor cells for proliferation // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t oxygen_saturation_for_proliferation = 38; + /// Glucose saturation level in tumor cells for cell growth: it affects + /// plorifereation since if cells are not big enough they will not divide + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t glucose_saturation_for_tumor_cell_growth = 0; + /// Limit of glucose level for tumor cell growth: it affects plorifereation + /// since if cells are not big enough they will not divide + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t glucose_limit_for_tumor_cell_growth = 0; /// Limit of oxygen level for tumor cell proliferation // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t oxygen_limit_for_proliferation = 10; @@ -232,23 +368,50 @@ struct SimParam : public ParamGroup { /// Limit of oxygen to maximum necrosis probability // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t oxygen_limit_for_necrosis_maximum = 2.5; + /// Maximum rate per minute of necrosis for tumor cells in case of hypoxia + /// with oxygen_limit_for_necrosis_maximum oxygen + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t maximum_necrosis_lack_of_oxygen_rate = 0.00277778; + /// Limit of glucose to start causing cell death + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t glucose_limit_for_death = 0; + /// Limit of glucose to maximum cell death probability + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t glucose_limit_for_death_maximum = 0; + /// Maximum rate per minute of cell death for tumor cells in case of lack of + /// nutrients with glucose_limit_for_death_maximum glucose + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t maximum_death_lack_of_glucose_rate = 0.0; + /// Basal cell death probability for cancer cells per minute: death due to + /// natural random causes + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t basal_death_probability_cancer_cells = 0.0; /// Time in minutes until a lysed necrotic cell is removed from the simulation // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t time_lysis = 86400; - /// Maximum rate per minute of necrosis for tumor cells in case of hypoxia - /// with 0 oxygen - // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t maximum_necrosis_rate = 0.00277778; /// Default oxygen consumption rate of tumor cell // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t default_oxygen_consumption_tumor_cell = 10; + /// Default glucose consumption rate of tumor cell + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t default_glucose_consumption_tumor_cell = 0.0007; /// Volume parameters /// Default total volume of a new tumor cell in μm³ // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t default_volume_new_tumor_cell = 2494; - /// Default volume of the nucleus of a new tumor cell in μm³ + /// Standard deviation of the total volume of a new tumor cell in μm³ + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t std_volume_new_tumor_cell = 0; + /// Maximum volume of a tumor cell in μm³ // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t default_volume_nucleus_tumor_cell = 540; + real_t max_volume_new_tumor_cell = 2494; + /// Minimum volume of a tumor cell in μm³ + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t min_volume_new_tumor_cell = 2494; + /// Default fraction of volume of the nucleus of a tumor cell. It is the + /// target fraction and the fraction a new born tumor cell will have + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t default_fraction_of_volume_for_nucleus_tumor_cell = 0.21652; /// Default fraction of fluid volume in a new tumor cell // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t default_fraction_fluid_tumor_cell = 0.75; @@ -258,6 +421,11 @@ struct SimParam : public ParamGroup { /// Standard Deviation for transformation Random Rate in hours // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t standard_deviation_transformation_random_rate = 3.7; + /// Minimum tumor cell volume for division. If the tumor cell volume is below + /// this value times the target volume, it will not divide. This is used to + /// avoid division of tumor cells that are too small due to lack of nutrients. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t minimum_tumor_cell_target_volume_fraction_for_division = 0.0; /// Average adhesion time in minutes for Tumor Cell under CAR-T attack before /// escaping // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) @@ -315,15 +483,32 @@ struct SimParam : public ParamGroup { /// Average time in minutes until a CAR-T cell dies // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t average_maximum_time_untill_apoptosis_cart = 12342.86; + real_t average_maximum_time_until_apoptosis_cart = 12342.86; /// Default oxygen consumption rate of CAR-T cell // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t default_oxygen_consumption_cart = 1; + /// Default glucose consumption rate of CAR-T cell + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t default_glucose_consumption_cart = 0.0007; /// Volume parameters /// Default total volume of a new CAR-T cell in μm³ // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t default_volume_new_cart_cell = 2494; - + /// Standard deviation of the total volume of a new CAR-T cell in μm³ + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t std_volume_new_cart_cell = 0; + /// Maximum volume of a new CAR-T cell in μm³ + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t max_volume_new_cart_cell = 2494; + /// Minimum volume of a new CAR-T cell in μm³ + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t min_volume_new_cart_cell = 2494; + /// Default fraction of volume not solid of a CAR-T cell. + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t default_fraction_fluid_cart_cell = 0.75; + /// Default fraction of volume of the nucleus of a CAR-T cell + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t default_fraction_of_volume_for_nucleus_cart_cell = 0.21652; /// How often a CAR-T cell tries to kill an attached cancer cell in 1/min // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t kill_rate_cart = 0.06667; @@ -345,11 +530,16 @@ struct SimParam : public ParamGroup { /// Average persistence time before CAR-T cell moves // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t persistence_time_cart = 10; // 10 minutes - /// Higher bias (\in [0,1]) makes CAR-T movement more directed toward - /// immunostimulatory factor source; while a bias of 0 makes the movement - /// random + /// Average migration bias of CAR-T cells. Higher bias (\in [0,1]) makes CAR-T + /// movement more directed toward immunostimulatory factor source; while a + /// bias of 0 makes the movement random // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t migration_bias_cart = 0.5; + real_t avg_migration_bias_cart = 0.5; + /// Standard deviation of the migration bias for CAR-T cells. Higher values + /// make the migration bias more variable among CAR-T cells + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + real_t std_migration_bias_cart = 0.0; + /// Migration speed // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t migration_speed_cart = 5; @@ -368,12 +558,15 @@ struct SimParam : public ParamGroup { /// Steps in one day // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) size_t steps_in_one_day = 14400; + /// Volume of a single voxel in μm³ // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) real_t voxel_volume = 8000; - /// 1-migration_bias_cart + + /// Squared max allowed radius for bounded space in μm² // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) - real_t migration_one_minus_bias_cart = 0.5; + real_t bounded_space_max_allowed_radius_squared = 1000000; + /// Probability of a CAR-T cell to migrate in a given /// mechanical time step // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) diff --git a/src/utils/utils_aux.cc b/src/utils/utils_aux.cc index 2f4e08b..bbfd3fe 100644 --- a/src/utils/utils_aux.cc +++ b/src/utils/utils_aux.cc @@ -31,6 +31,7 @@ #include "core/resource_manager.h" #include "core/scheduler.h" #include "core/simulation.h" +#include "core/util/log.h" #include "core/util/math.h" #include "core/util/random.h" #include @@ -97,15 +98,60 @@ std::vector CreateSphereOfTumorCells(real_t sphere_radius) { return positions; } +/// Create a cylindrical arrangement of tumor cells randomly distributed within +/// its volume +std::vector CreateCylinderOfTumorCells(real_t cylinder_radius, + real_t cylinder_height, + size_t number_of_cells) { + std::vector positions; + positions.reserve(number_of_cells); + + Random* random = Simulation::GetActive()->GetRandom(); + + for (size_t i = 0; i < number_of_cells; ++i) { + // sqrt ensures uniform distribution over the disk area (area element = r dr + // dθ) + const real_t angle = random->Uniform(0.0, kTwicePi); + const real_t radius = + cylinder_radius * std::sqrt(random->Uniform(0.0, 1.0)); + // NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers) + const real_t height = + random->Uniform(-cylinder_height / 2.0, cylinder_height / 2.0); + + positions.push_back( + {radius * std::cos(angle), radius * std::sin(angle), height}); + } + + return positions; +} + // Function to compute the number of tumor cells of each type, the radius of the -// tumor, the average oncoprotein level and the average oxygen level of the -// cancer cells. -std::tuple -AnalyzeTumor() { +// tumor, the number of alive and dead cart cells, the average oncoprotein +// level, the average oxygen level of the cancer cells, the average oxygen level +// of all cells, the average glucose level of the cancer cells, the average +// glucose level of all cells and the average radius distance from the center of +// the tumor of the living cart cells. Only considering cells within the +// specified inner and outer radius. +std::tuple +// NOLINTNEXTLINE(readability-function-cognitive-complexity) +AnalyzeTumor(real_t inner_radius_considered, real_t outer_radius_considered) { + real_t inner_radius_squared = + inner_radius_considered * inner_radius_considered; + real_t outer_radius_squared = + outer_radius_considered * outer_radius_considered; + Simulation* sim = Simulation::GetActive(); ResourceManager* rm = sim->GetResourceManager(); DiffusionGrid* oxygen_dgrid = rm->GetDiffusionGrid("oxygen"); + const auto* sparams = sim->GetParam()->Get(); + + // Pointer to glucose diffusion grid only if it is added in the simulation + DiffusionGrid* glucose_dgrid = nullptr; + if (sparams->add_glucose) { + glucose_dgrid = rm->GetDiffusionGrid("glucose"); + } int total_num_tumor_cells = 0; int num_tumor_cells_type1 = 0; @@ -113,23 +159,67 @@ AnalyzeTumor() { int num_tumor_cells_type3 = 0; int num_tumor_cells_type4 = 0; int num_tumor_cells_type5_dead = 0; + int num_dead_lack_of_oxygen = 0; + int num_dead_lack_of_glucose = 0; + int num_dead_random_natural_causes = 0; + int num_dead_cartcell_kill = 0; int num_alive_cart = 0; + int num_dead_cart = 0; real_t max_dist_sq = 0.0; real_t acumulator_oncoprotein = 0.0; real_t acumulator_oxygen_cancer_cells = 0.0; + real_t acumulator_oxygen_all_cells = 0.0; + real_t acumulator_glucose_cancer_cells = 0.0; + real_t acumulator_glucose_all_cells = 0.0; + real_t acumulator_squared_radius_distance_living_cart_cells = 0.0; + + TumorShape tumor_shape = sparams->tumor_shape; rm->ForEachAgent([&](const Agent* agent) { + // Compute the distance to the center depending on the tumor shape + const Real3& pos = agent->GetPosition(); + real_t dist_sq = 0.0; + switch (tumor_shape) { + case TumorShape::kCylinder: { + // Only consider x and y for cylindrical rumor, distance to the axis of + // the cylinder + dist_sq = pos[0] * pos[0] + pos[1] * pos[1]; + break; + } + case TumorShape::kSphere: { + // Consider all three dimensions for spherical radius + dist_sq = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]; + break; + } + default: + Log::Error("AnalyzeTumor", + "Unknown tumor shape, please use 'sphere' or 'cylinder'."); + break; + } + + // If the cell is outside the considered radius range, skip it + if (dist_sq < inner_radius_squared || dist_sq >= outer_radius_squared) { + return; + } + + // The agent is within the considered radius range, analyze it if (const auto* tumor_cell = dynamic_cast(agent)) { + // Accumulate oxygen level for average calculation + acumulator_oxygen_all_cells += oxygen_dgrid->GetValue(pos); + total_num_tumor_cells++; - const Real3& pos = agent->GetPosition(); // Accumulate oxygen level for average calculation acumulator_oxygen_cancer_cells += oxygen_dgrid->GetValue(pos); + if (glucose_dgrid != nullptr) { + // Accumulate glucose level for average calculation + acumulator_glucose_all_cells += glucose_dgrid->GetValue(pos); + acumulator_glucose_cancer_cells += glucose_dgrid->GetValue(pos); + } + // computing tumor radius - const real_t dist_sq = - pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]; if (dist_sq > max_dist_sq) { max_dist_sq = dist_sq; } @@ -155,13 +245,40 @@ AnalyzeTumor() { break; case TumorCellType::kType5: num_tumor_cells_type5_dead++; + // Count dead tumor cells by cause of death + switch (tumor_cell->GetDeathCause()) { + case TumorCellDeathCause::kLackOfOxygen: + num_dead_lack_of_oxygen++; + break; + case TumorCellDeathCause::kLackOfGlucose: + num_dead_lack_of_glucose++; + break; + case TumorCellDeathCause::kRandomNaturalCauses: + num_dead_random_natural_causes++; + break; + case TumorCellDeathCause::kCartKill: + num_dead_cartcell_kill++; + break; + default: + break; + } break; default: break; } } else if (const auto* cart_cell = dynamic_cast(agent)) { + // Accumulate oxygen level for average calculation + acumulator_oxygen_all_cells += oxygen_dgrid->GetValue(pos); + if (glucose_dgrid != nullptr) { + // Accumulate glucose level for average calculation + acumulator_glucose_all_cells += glucose_dgrid->GetValue(pos); + } if (cart_cell->GetState() == CarTCellState::kAlive) { num_alive_cart++; + acumulator_squared_radius_distance_living_cart_cells += dist_sq; + } else if (cart_cell->GetState() == CarTCellState::kApoptotic) { + // Dead CART cell + num_dead_cart++; } } }); @@ -174,11 +291,48 @@ AnalyzeTumor() { (total_num_tumor_cells > 0) ? (acumulator_oxygen_cancer_cells / total_num_tumor_cells) : 0.0; - return {total_num_tumor_cells, num_tumor_cells_type1, - num_tumor_cells_type2, num_tumor_cells_type3, - num_tumor_cells_type4, num_tumor_cells_type5_dead, - num_alive_cart, std::sqrt(max_dist_sq), - average_oncoprotein, average_oxygen_cancer_cells}; + + const real_t average_oxygen_all_cells = + (total_num_tumor_cells + num_alive_cart + num_dead_cart > 0) + ? (acumulator_oxygen_all_cells / + (total_num_tumor_cells + num_alive_cart + num_dead_cart)) + : 0.0; + const real_t average_glucose_cancer_cells = + (total_num_tumor_cells > 0 && glucose_dgrid != nullptr) + ? (acumulator_glucose_cancer_cells / total_num_tumor_cells) + : 0.0; + const real_t average_glucose_all_cells = + (total_num_tumor_cells + num_alive_cart + num_dead_cart > 0 && + glucose_dgrid != nullptr) + ? (acumulator_glucose_all_cells / + (total_num_tumor_cells + num_alive_cart + num_dead_cart)) + : 0.0; + + const real_t average_radius_distance_living_cart_cells = + (num_alive_cart > 0) + ? std::sqrt(acumulator_squared_radius_distance_living_cart_cells / + num_alive_cart) + : 0.0; + + return {total_num_tumor_cells, + num_tumor_cells_type1, + num_tumor_cells_type2, + num_tumor_cells_type3, + num_tumor_cells_type4, + num_tumor_cells_type5_dead, + num_dead_lack_of_oxygen, + num_dead_lack_of_glucose, + num_dead_random_natural_causes, + num_dead_cartcell_kill, + num_alive_cart, + num_dead_cart, + std::sqrt(max_dist_sq), + average_oncoprotein, + average_oxygen_cancer_cells, + average_oxygen_all_cells, + average_glucose_cancer_cells, + average_glucose_all_cells, + average_radius_distance_living_cart_cells}; } // Function to output summary CSV @@ -189,27 +343,32 @@ void OutputSummary::operator()() { const uint64_t current_step = scheduler->GetSimulatedSteps(); if (current_step % frequency_ == 0) { + // Calculate time in days, hours, minutes + const double total_minutes = + static_cast(current_step) * sparams->dt_step; + const double total_hours = total_minutes / kMinutesInAnHour; + const double total_days = total_hours / kHoursInADay; + // Delete csv content current_step == 0 to, otherwise append mode std::ofstream file("output/final_data.csv", current_step == 0 ? std::ios::trunc : 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,num_alive_" - "cart,average_oncoprotein,average_oxygen_cancer_cells\n"; // Header - // for - // CSV - // file + 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," + "dead_tumor_cells_lack_of_oxygen,dead_tumor_cells_lack_of_" + "glucose,dead_tumor_cells_random_natural_causes,dead_tumor_" + "cells_cartcell_kill,num_alive_cart,num_dead_cart,average_" + "oncoprotein,average_oxygen_cancer_" + "cells,average_oxygen_all_cells,average_glucose_cancer_cells," + "average_glucose_all_cells,average_radius_distance_living_cart_" + "cells\n"; // Header + // for + // CSV + // file } - // Calculate time in days, hours, minutes - const double total_minutes = - static_cast(current_step) * sparams->dt_step; - const double total_hours = total_minutes / kMinutesInAnHour; - const double total_days = total_hours / kHoursInADay; - // Count total cells, tumor cells of each type and tumor radius int total_num_tumor_cells = 0; int num_tumor_cells_type1 = 0; @@ -217,15 +376,30 @@ void OutputSummary::operator()() { int num_tumor_cells_type3 = 0; int num_tumor_cells_type4 = 0; int num_tumor_cells_type5_dead = 0; + int num_dead_lack_of_oxygen = 0; + int num_dead_lack_of_glucose = 0; + int num_dead_random_natural_causes = 0; + int num_dead_cartcell_kill = 0; int num_alive_cart = 0; + int num_dead_cart = 0; real_t tumor_radius = 0.0; real_t average_oncoprotein = 0.0; real_t average_oxygen_cancer_cells = 0.0; + real_t average_oxygen_all_cells = 0.0; + real_t average_glucose_cancer_cells = 0.0; + real_t average_glucose_all_cells = 0.0; + real_t average_radius_distance_living_cart_cells = 0.0; + // Analyze tumor with no limits on inner and outer 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, - num_alive_cart, tumor_radius, average_oncoprotein, - average_oxygen_cancer_cells) = AnalyzeTumor(); + num_dead_lack_of_oxygen, num_dead_lack_of_glucose, + num_dead_random_natural_causes, num_dead_cartcell_kill, + num_alive_cart, num_dead_cart, tumor_radius, average_oncoprotein, + average_oxygen_cancer_cells, average_oxygen_all_cells, + average_glucose_cancer_cells, average_glucose_all_cells, + average_radius_distance_living_cart_cells) = + AnalyzeTumor(0, sparams->bounded_space_length); size_t total_num_cells = simulation->GetResourceManager()->GetNumAgents(); // If a dosage is administred this exact time the numbers are not seen in @@ -248,12 +422,154 @@ void OutputSummary::operator()() { << total_num_tumor_cells << "," << num_tumor_cells_type1 << "," << num_tumor_cells_type2 << "," << num_tumor_cells_type3 << "," << num_tumor_cells_type4 << "," << num_tumor_cells_type5_dead << "," - << num_alive_cart << "," << average_oncoprotein << "," - << average_oxygen_cancer_cells << "\n"; + << num_dead_lack_of_oxygen << "," << num_dead_lack_of_glucose << "," + << num_dead_random_natural_causes << "," << num_dead_cartcell_kill + << "," << num_alive_cart << "," << num_dead_cart << "," + << average_oncoprotein << "," << average_oxygen_cancer_cells << "," + << average_oxygen_all_cells << "," << average_glucose_cancer_cells + << "," << average_glucose_all_cells << "," + << average_radius_distance_living_cart_cells << "\n"; + } + + if (sparams->output_information_dependent_on_radius) { + OutputInformationBasedOnRadiusCSV(current_step, total_minutes, + total_hours, total_days); } } } +// Function to output information based on radius to CSV +void OutputInformationBasedOnRadiusCSV(const uint64_t current_step, + const real_t total_minutes, + const real_t total_hours, + const real_t total_days) { + const Simulation* simulation = Simulation::GetActive(); + const auto* sparams = simulation->GetParam()->Get(); + const real_t interval_size = + sparams->max_radius_analysis_csv_dependent_on_radius / + sparams->num_radius_intervals; + + // Delete csv content current_step == 0 to, otherwise append mode + std::ofstream file("output/data_dependent_on_radius_tumor.csv", + current_step == 0 ? std::ios::trunc : std::ios::app); + if (file.is_open()) { + if (current_step == 0) { + // Header for CSV file + file << "total_days,total_hours,total_minutes"; + for (int i = 0; i < sparams->num_radius_intervals; ++i) { + file << ",num_alive_cells_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",num_dead_cells_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",total_num_cells_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",num_alive_tumor_cells_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",total_num_tumor_cells_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",num_alive_cart_cells_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",num_dead_cart_cells_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",tumor_cells_type1_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",tumor_cells_type2_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",tumor_cells_type3_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",tumor_cells_type4_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",tumor_cells_type5_dead_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",dead_tumor_cells_lack_of_oxygen_radius_" << i * interval_size + << "_to_" << (i + 1) * interval_size; + file << ",dead_tumor_cells_lack_of_glucose_radius_" << i * interval_size + << "_to_" << (i + 1) * interval_size; + file << ",dead_tumor_cells_random_natural_causes_radius_" + << i * interval_size << "_to_" << (i + 1) * interval_size; + file << ",dead_tumor_cells_cartcell_kill_radius_" << i * interval_size + << "_to_" << (i + 1) * interval_size; + file << ",average_oncoprotein_radius_" << i * interval_size << "_to_" + << (i + 1) * interval_size; + file << ",average_oxygen_cancer_cells_radius_" << i * interval_size + << "_to_" << (i + 1) * interval_size; + file << ",average_oxygen_all_cells_radius_" << i * interval_size + << "_to_" << (i + 1) * interval_size; + file << ",average_glucose_cancer_cells_radius_" << i * interval_size + << "_to_" << (i + 1) * interval_size; + file << ",average_glucose_all_cells_radius_" << i * interval_size + << "_to_" << (i + 1) * interval_size; + file << ",average_radius_distance_living_cart_cells_radius_" + << i * interval_size << "_to_" << (i + 1) * interval_size; + } + // End of header line + file << "\n"; + } + // Write the time information for the current step + file << total_days << "," << total_hours << "," << total_minutes; + + for (int i = 0; i < sparams->num_radius_intervals; ++i) { + const real_t inner_radius_considered = i * interval_size; + const real_t outer_radius_considered = (i + 1) * interval_size; + + // Count total cells, tumor cells of each type and tumor radius + int total_num_tumor_cells = 0; + int num_tumor_cells_type1 = 0; + int num_tumor_cells_type2 = 0; + int num_tumor_cells_type3 = 0; + int num_tumor_cells_type4 = 0; + int num_tumor_cells_type5_dead = 0; + int num_dead_lack_of_oxygen = 0; + int num_dead_lack_of_glucose = 0; + int num_dead_random_natural_causes = 0; + int num_dead_cartcell_kill = 0; + int num_alive_cart = 0; + int num_dead_cart = 0; + real_t tumor_radius = 0.0; + real_t average_oncoprotein = 0.0; + real_t average_oxygen_cancer_cells = 0.0; + real_t average_oxygen_all_cells = 0.0; + real_t average_glucose_cancer_cells = 0.0; + real_t average_glucose_all_cells = 0.0; + real_t average_radius_distance_living_cart_cells = 0.0; + // Analyze tumor with no limits on inner and outer 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, + num_dead_lack_of_oxygen, num_dead_lack_of_glucose, + num_dead_random_natural_causes, num_dead_cartcell_kill, + num_alive_cart, num_dead_cart, tumor_radius, average_oncoprotein, + average_oxygen_cancer_cells, average_oxygen_all_cells, + average_glucose_cancer_cells, average_glucose_all_cells, + average_radius_distance_living_cart_cells) = + AnalyzeTumor(inner_radius_considered, outer_radius_considered); + + const int num_alive_tumor_cells = + num_tumor_cells_type1 + num_tumor_cells_type2 + + num_tumor_cells_type3 + num_tumor_cells_type4; + const int num_alive_cells = num_alive_cart + num_alive_tumor_cells; + const int num_dead_cells = num_dead_cart + num_tumor_cells_type5_dead; + const int total_num_cells = num_alive_cells + num_dead_cells; + + // Write data to CSV file + file << "," << num_alive_cells << "," << num_dead_cells << "," + << total_num_cells << "," << num_alive_tumor_cells << "," + << total_num_tumor_cells << "," << num_alive_cart << "," + << num_dead_cart << "," << num_tumor_cells_type1 << "," + << num_tumor_cells_type2 << "," << num_tumor_cells_type3 << "," + << num_tumor_cells_type4 << "," << num_tumor_cells_type5_dead << "," + << num_dead_lack_of_oxygen << "," << num_dead_lack_of_glucose << "," + << num_dead_random_natural_causes << "," << num_dead_cartcell_kill + << "," << average_oncoprotein << "," << average_oxygen_cancer_cells + << "," << average_oxygen_all_cells << "," + << average_glucose_cancer_cells << "," << average_glucose_all_cells + << "," << average_radius_distance_living_cart_cells; + } + // End of line for the current step + file << "\n"; + } +}; + // Function to spawn CAR-T cell dosages void SpawnCart::operator()() { Simulation* simulation = Simulation::GetActive(); @@ -280,12 +596,31 @@ void SpawnCart::operator()() { // compute tumor radius ResourceManager* rm = simulation->GetResourceManager(); real_t max_dist_sq = 0.0; + TumorShape tumor_shape = sparams->tumor_shape; rm->ForEachAgent([&](const Agent* agent) { if (const auto* cancer_cell = dynamic_cast(agent)) { + // Compute the distance to the center depending on the tumor shape + real_t dist_sq = 0.0; const Real3& pos = cancer_cell->GetPosition(); - const real_t dist_sq = - pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]; + switch (tumor_shape) { + case TumorShape::kCylinder: { + // Only consider x and y for cylindrical rumor, distance to the axis + // of the cylinder + dist_sq = pos[0] * pos[0] + pos[1] * pos[1]; + break; + } + case TumorShape::kSphere: { + // Consider all three dimensions for spherical radius + dist_sq = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]; + break; + } + default: + Log::Error( + "SpawnCart", + "Unknown tumor shape, please use 'sphere' or 'cylinder'."); + break; + } if (dist_sq > max_dist_sq) { max_dist_sq = dist_sq; } @@ -303,8 +638,12 @@ void SpawnCart::operator()() { // for generating car-t positions Random* rng = simulation->GetRandom(); const Param* param = simulation->GetParam(); + const auto* sparams = param->Get(); const real_t min_b = param->min_bound; const real_t max_b = param->max_bound; + const real_t min_bz = sparams->bounded_space_min_allowed_z; + const real_t max_bz = sparams->bounded_space_max_allowed_z; + const real_t max_r_sq = sparams->bounded_space_max_allowed_radius_squared; real_t px = 0.0; real_t py = 0.0; real_t pz = 0.0; @@ -317,9 +656,26 @@ void SpawnCart::operator()() { while (true) { px = rng->Uniform(min_b, max_b); py = rng->Uniform(min_b, max_b); - pz = rng->Uniform(min_b, max_b); - radi_sq = px * px + py * py + pz * pz; - if (radi_sq >= minimum_squared_radius) { + pz = rng->Uniform(min_bz, max_bz); + switch (tumor_shape) { + case TumorShape::kCylinder: { + // Only consider x and y for cylindrical rumor, distance to the axis + // of the cylinder + radi_sq = px * px + py * py; + break; + } + case TumorShape::kSphere: { + // Consider all three dimensions for spherical radius + radi_sq = px * px + py * py + pz * pz; + break; + } + default: + Log::Error( + "SpawnCart", + "Unknown tumor shape, please use 'sphere' or 'cylinder'."); + break; + } + if (radi_sq >= minimum_squared_radius && radi_sq <= max_r_sq) { break; } } diff --git a/src/utils/utils_aux.h b/src/utils/utils_aux.h index 3cc7a53..028ef0f 100644 --- a/src/utils/utils_aux.h +++ b/src/utils/utils_aux.h @@ -1,5 +1,5 @@ /* - * Copyright 2025 compiler-research.org, Salvador de la Torre Gonzalez, Luciana + * Copyright 2026 compiler-research.org, Salvador de la Torre Gonzalez, Luciana * Melina Luque * * Licensed under the Apache License, Version 2.0 (the "License"); @@ -56,11 +56,27 @@ real_t SamplePositiveGaussian(real_t mean, real_t sigma); /// @return Vector of 3D positions where tumor cells should be placed std::vector CreateSphereOfTumorCells(real_t sphere_radius); +/// Create a cylindrical arrangement of tumor cells +/// +/// Generates a vector of 3D positions for tumor cells arranged in a cylindrical +/// pattern with the specified radius and height. The cells are positioned +/// randomly uniformly distributed within the cylinder's volume. +/// +/// @param cylinder_radius Radius of the cylindrical tumor in micrometers +/// @param cylinder_height Height of the cylindrical tumor in micrometers +/// @param number_of_cells Number of tumor cells to place +/// @return Vector of 3D positions where tumor cells should be placed +std::vector CreateCylinderOfTumorCells(real_t cylinder_radius, + real_t cylinder_height, + size_t number_of_cells); + /// Compute tumor statistics and characteristics /// -/// Analyzes the current tumor population to compute the number of tumor cells +/// Analyzes the current tumor population to compute the number of cells /// of each type and the overall radius of the tumor mass. In addition, it -/// computes the average oncoprotein level and oxygen across all tumor cells. +/// computes the average oncoprotein level and oxygen across all tumor cells and +/// all cells. Only cells within the specified inner and outer radius are +/// considered for the analysis. /// /// @return Tuple containing: /// - Total number of tumor cells @@ -69,13 +85,23 @@ std::vector CreateSphereOfTumorCells(real_t sphere_radius); /// - Number of type 3 tumor cells /// - Number of type 4 tumor cells (least aggressive) /// - Number of type 5 tumor cells (dead) +/// - Number of dead tumor cells that died because of lack of oxygen +/// - Number of dead tumor cells that died because of lack of glucose +/// - Number of dead tumor cells that died because of random natural causes +/// - Number of dead tumor cells that were killed by CAR-T cells /// - Number of living CAR-T cells (not apoptotic) +/// - Number of dead CAR-T cells (apoptotic) /// - Current tumor radius in micrometers /// - Average oncoprotein level across all tumor cells /// - Average oxygen level across all tumor cells -std::tuple -AnalyzeTumor(); +/// - Average oxygen level across all cells +/// - Average glucose level across all tumor cells +/// - Average glucose level across all cells +/// - Average radial position of the living cart cells from the tumor center +std::tuple +AnalyzeTumor(real_t inner_radius_considered, real_t outer_radius_considered); /// Generates a random direction unitary vector /// @@ -133,6 +159,16 @@ struct OutputSummary : public StandaloneOperationImpl { // NOLINTNEXTLINE(cppcoreguidelines-owning-memory) inline BDM_REGISTER_OP(OutputSummary, "OutputSummary", kCpu); +/// Function to output information based on radius to CSV, agregate values for +/// each radius from the center of the cylindrical tumor +/// @param current_step Current simulation step +/// @param total_minutes Total simulation time in minutes +/// @param total_hours Total simulation time in hours +/// @param total_days Total simulation time in days +void OutputInformationBasedOnRadiusCSV(uint64_t current_step, + real_t total_minutes, real_t total_hours, + real_t total_days); + } // namespace bdm #endif // CORE_UTIL_UTILS_AUX_H_ diff --git a/test/test-diffusion.cc b/test/test-diffusion.cc index dc5ddf5..7557841 100644 --- a/test/test-diffusion.cc +++ b/test/test-diffusion.cc @@ -26,7 +26,7 @@ class DiffusionThomasAlgorithmTest : public ::testing::Test { // case: valid index computes correctly TEST_F(DiffusionThomasAlgorithmTest, GetBoxIndexValidCoordinates) { - DiffusionThomasAlgorithm grid(0, "test", 1.0, 0.1, 10, 0.01, false); + DiffusionThomasAlgorithm grid(0, "test", 1.0, 0.1, 10, 0.01, false, true); EXPECT_EQ(grid.GetBoxIndex(4, 3, 2), static_cast(234)); // z * res * res + y * res + x = 2*100 + 3*10 + 4 = 234 }