New model GSoC2026

.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/

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()

bdm.toml

@@ -1,21 +1,18 @@
 [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]]
 name = "CarTCell"
 additional_data_members = ["diameter_","volume_"]
 
 [[visualize_diffusion]]
-name = "oxygen"
-
-[[visualize_diffusion]]
-name = "immunostimulatory_factor"
+name = "oxygen"

params.json

@@ -7,4 +7,4 @@
     "0": 3957,
     "8": 3957
   }
-}
+}

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
+}

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
+  }
+}