SOCR ≫ BPAD Website ≫ BPAD GitHub ≫

1 Introduction

This BPAD section in Chapter 3 is focused on Tumor Growth Simulation and Logistic Growth Modeling. It demonstrates cancer data simulation and modeling of tumor growth under different treatment scenarios (radiation therapy and chemotherapy). Readers are encouraged to review the SOCR - DSPA sections on Dynamical Systems and Stochastic Differential Equation Modeling (in R-Julia) and Cancer Treatment Stochastic Differential Equation (SDE) Modeling, as well as the BPAD section on Exponential Growth and Decay.

The data-driven simulation and modeling accounts for:

  1. Local Control (LC): A measure of tumor response to treatment,
  2. Toxicity (RP2): Radiation pneumonitis Grade 2+ toxicity, a side effect of radiation therapy, and
  3. Additional Biological Variables: Biomarkers, imaging features, and patient-specific factors.

The model allows for forward prediction of outcomes given user-specified time-periods.

2 Experimental settings

2.1 Parameter Configuration

Below we show the core parameter settings, which can be modified to explore different conditions and treatment regimens.

# Set random seed for reproducibility
set.seed(12345)

# Simulation parameters
n_patients <- 50      # Number of patients to simulate
days <- 90            # Days of observation
prediction_days <- 30 # Days to predict forward

# Treatment parameters
radiation_doses <- c("Low" = 60, "Medium" = 70, "High" = 80) # Gy
chemo_regimens <- c("None" = 0, "Standard" = 1, "Intensive" = 2)

# Logistic growth parameters - base values
K_base <- 100        # Carrying capacity (maximum tumor size in mm)
r_base <- 0.08       # Intrinsic growth rate
N0_base <- 10        # Initial tumor size (mm)

# Treatment effect parameters
# How much each treatment reduces the growth rate
radiation_effect <- c("Low" = 0.3, "Medium" = 0.5, "High" = 0.7)
chemo_effect <- c("None" = 0, "Standard" = 0.4, "Intensive" = 0.6)

# Toxicity parameters
toxicity_radiation <- c("Low" = 0.1, "Medium" = 0.25, "High" = 0.4)
toxicity_chemo <- c("None" = 0, "Standard" = 0.15, "Intensive" = 0.35)

# Patient-specific biological factors (influence on growth and toxicity)
biomarker_names <- c("Immune_Score", "Hypoxia", "Proliferation_Index", 
                     "Vascular_Density", "PD_L1_Expression")
imaging_names <- c("SUV_max", "Tumor_Heterogeneity", "Metabolic_Volume")
clinical_names <- c("Age", "Performance_Status", "Smoking_Pack_Years", "Comorbidity_Score")

# Factor effects on carrying capacity, growth rate, and toxicity
bio_K_effect <- c(0.2, -0.15, 0.25, 0.1, -0.2)      # Effects on carrying capacity
bio_r_effect <- c(-0.1, 0.2, 0.15, 0.1, -0.15)      # Effects on growth rate
bio_tox_effect <- c(-0.2, 0.15, 0.1, 0.2, 0.1)      # Effects on toxicity

img_K_effect <- c(0.15, 0.2, 0.25)                  # Effects on carrying capacity
img_r_effect <- c(0.1, 0.15, 0.2)                   # Effects on growth rate
img_tox_effect <- c(0.15, 0.1, 0.2)                 # Effects on toxicity

clin_K_effect <- c(0.1, -0.2, 0.15, 0.2)            # Effects on carrying capacity
clin_r_effect <- c(0.05, -0.15, 0.1, 0.15)          # Effects on growth rate
clin_tox_effect <- c(0.15, 0.2, 0.25, 0.3)          # Effects on toxicity

2.2 Simulation Functions

Here we define the core finctions for simulating realsitic cancer datasets.

# Function to generate patient-specific biological factors
generate_patient_factors <- function() {
  # Generate standardized biomarkers (mean 0, sd 1)
  biomarkers <- rnorm(length(biomarker_names), mean = 0, sd = 1)
  names(biomarkers) <- biomarker_names
  
  # Generate imaging features (mean 0, sd 1)
  imaging <- rnorm(length(imaging_names), mean = 0, sd = 1)
  names(imaging) <- imaging_names
  
  # Generate clinical factors (appropriate distributions)
  clinical <- c(
    Age = rnorm(1, mean = 65, sd = 10),                          # Age ~ N(65, 10)
    Performance_Status = sample(0:2, 1, prob = c(0.3, 0.5, 0.2)), # 0-2 scale
    Smoking_Pack_Years = max(0, rnorm(1, mean = 30, sd = 15)),     # Pack years (>=0)
    Comorbidity_Score = sample(0:4, 1, prob = c(0.3, 0.3, 0.2, 0.1, 0.1)) # 0-4 scale
  )
  
  # Normalize clinical factors
  clinical_norm <- c(
    Age = (clinical["Age"] - 65) / 10,
    Performance_Status = clinical["Performance_Status"] / 2,
    Smoking_Pack_Years = (clinical["Smoking_Pack_Years"] - 30) / 15,
    Comorbidity_Score = clinical["Comorbidity_Score"] / 4
  )
  
  # Combined factors (original values for reporting, normalized for calculations)
  return(list(
    biomarkers = biomarkers,
    imaging = imaging,
    clinical = clinical,
    clinical_norm = clinical_norm
  ))
}

# Function to calculate patient-specific growth parameters
calculate_growth_parameters <- function(factors, rad_dose_level, chemo_regimen_level) {
  # Base values
  K <- K_base
  r <- r_base
  N0 <- N0_base
  
  # Apply biological factor effects - biomarkers
  for (i in 1:length(biomarker_names)) {
    K <- K * (1 + bio_K_effect[i] * factors$biomarkers[i])
    r <- r * (1 + bio_r_effect[i] * factors$biomarkers[i])
  }
  
  # Apply biological factor effects - imaging
  for (i in 1:length(imaging_names)) {
    K <- K * (1 + img_K_effect[i] * factors$imaging[i])
    r <- r * (1 + img_r_effect[i] * factors$imaging[i])
  }
  
  # Apply biological factor effects - clinical
  for (i in 1:length(clinical_names)) {
    K <- K * (1 + clin_K_effect[i] * factors$clinical_norm[i])
    r <- r * (1 + clin_r_effect[i] * factors$clinical_norm[i])
  }
  
  # Apply treatment effects
  r_reduced <- r * (1 - radiation_effect[rad_dose_level] - chemo_effect[chemo_regimen_level])
  
  # Ensure logical ranges
  K <- max(K, 1.2 * N0)  # Carrying capacity should be larger than initial size
  r_reduced <- max(r_reduced, 0.001)  # Growth rate should be positive but can be very small
  
  return(list(K = K, r = r_reduced, N0 = N0))
}

# Function to calculate toxicity
calculate_toxicity <- function(factors, rad_dose_level, chemo_regimen_level) {
  # Base toxicity from treatments
  tox_base <- toxicity_radiation[rad_dose_level] + toxicity_chemo[chemo_regimen_level]
  
  # Apply biological factor effects - biomarkers
  for (i in 1:length(biomarker_names)) {
    tox_base <- tox_base * (1 + bio_tox_effect[i] * factors$biomarkers[i])
  }
  
  # Apply biological factor effects - imaging
  for (i in 1:length(imaging_names)) {
    tox_base <- tox_base * (1 + img_tox_effect[i] * factors$imaging[i])
  }
  
  # Apply biological factor effects - clinical
  for (i in 1:length(clinical_names)) {
    tox_base <- tox_base * (1 + clin_tox_effect[i] * factors$clinical_norm[i])
  }
  
  # Add random variation and bound between 0 and 1
  toxicity <- pmin(1, pmax(0, tox_base + rnorm(1, 0, 0.05)))
  
  # Convert to binary outcome with probability
  toxicity_bin <- rbinom(1, 1, toxicity)
  
  return(list(
    toxicity_prob = toxicity,
    toxicity_outcome = toxicity_bin
  ))
}

# Logistic growth function
logistic_growth <- function(t, K, r, N0) {
  return(K / (1 + ((K / N0) - 1) * exp(-r * t)))
}

# Add noise to the growth curve to simulate measurement error and biological variation
add_noise <- function(value, cv = 0.05) {
  return(value * (1 + rnorm(1, 0, cv)))
}

# Function to generate tumor growth data for one patient
generate_patient_data <- function(patient_id, days, rad_dose_level, chemo_regimen_level) {
  # Generate patient factors
  factors <- generate_patient_factors()
  
  # Calculate growth parameters
  params <- calculate_growth_parameters(factors, rad_dose_level, chemo_regimen_level)
  
  # Calculate toxicity
  toxicity <- calculate_toxicity(factors, rad_dose_level, chemo_regimen_level)
  
  # Generate tumor size at each time point
  data <- data.frame(
    patient_id = patient_id,
    day = 0:days,
    radiation_dose = rad_dose_level,
    radiation_dose_value = radiation_doses[rad_dose_level],
    chemo_regimen = chemo_regimen_level,
    chemo_regimen_value = chemo_regimens[chemo_regimen_level]
  )
  
  # Add predicted tumor size based on logistic growth
  data$tumor_size_pred <- sapply(data$day, function(t) logistic_growth(t, params$K, params$r, params$N0))
  
  # Add noise to create observed tumor size
  data$tumor_size <- sapply(data$tumor_size_pred, add_noise)
  
  # Add growth parameters
  data$K <- params$K
  data$r <- params$r
  data$N0 <- params$N0
  
  # Add toxicity data
  data$toxicity_prob <- toxicity$toxicity_prob
  data$toxicity_outcome <- toxicity$toxicity_outcome
  
  # Add all biological factors
  for (name in biomarker_names) {
    data[[name]] <- factors$biomarkers[name]
  }
  
  for (name in imaging_names) {
    data[[name]] <- factors$imaging[name]
  }
  
  for (name in clinical_names) {
    data[[name]] <- factors$clinical[name]
  }
  
  # Calculate Local Control (LC) - inverse relationship with tumor size relative to baseline
  data$local_control <- with(data, 100 * (1 - (tumor_size / tumor_size[1]) / (K / tumor_size[1])))
  
  # Ensure LC is between 0 and 100
  data$local_control <- pmin(100, pmax(0, data$local_control))
  
  return(data)
}

2.3 Generate Simulated Patient Data

Next we simulate cancer patient datasets.

# Assign treatments to patients
patient_treatments <- data.frame(
  patient_id = 1:n_patients,
  rad_dose_level = sample(names(radiation_doses), n_patients, replace = TRUE),
  chemo_regimen_level = sample(names(chemo_regimens), n_patients, replace = TRUE)
)

# Generate data for all patients
patient_data_list <- map2(
  patient_treatments$patient_id,
  patient_treatments$rad_dose_level,
  function(id, rad) {
    map(
      patient_treatments$chemo_regimen_level[patient_treatments$patient_id == id],
      function(chemo) {
        generate_patient_data(id, days, rad, chemo)
      }
    )
  }
) %>% flatten()

# Combine all patient data
all_patient_data <- bind_rows(patient_data_list)

# Display summary of the data
summary_data <- all_patient_data %>%
  filter(day == days) %>%
  group_by(radiation_dose, chemo_regimen) %>%
  summarize(
    n_patients = n(),
    mean_tumor_size = mean(tumor_size),
    mean_local_control = mean(local_control),
    toxicity_rate = mean(toxicity_outcome),
    .groups = 'drop'
  )

knitr::kable(summary_data, caption = "Summary of Patient Outcomes at End of Observation Period")
Summary of Patient Outcomes at End of Observation Period
radiation_dose chemo_regimen n_patients mean_tumor_size mean_local_control toxicity_rate
High Intensive 7 10.64339 70.77827 0.7142857
High None 4 40.97723 51.81287 0.5000000
High Standard 6 10.96377 89.20703 0.5000000
Low Intensive 4 17.68211 69.86646 0.5000000
Low None 5 59.98504 17.98194 0.0000000
Low Standard 5 57.71604 47.91041 0.6000000
Medium Intensive 10 11.16218 88.53072 0.6000000
Medium None 4 50.70305 17.38160 0.2500000
Medium Standard 5 14.93581 75.98805 0.6000000

2.4 Visualize Tumor Growth Curves

Let’s show the tumor growth/decay curves for different scenarios.

# Plot average tumor growth by treatment combination
avg_growth <- all_patient_data %>%
  group_by(day, radiation_dose, chemo_regimen) %>%
  summarize(
    mean_size = mean(tumor_size),
    lower_ci = quantile(tumor_size, 0.25),
    upper_ci = quantile(tumor_size, 0.75),
    .groups = 'drop'
  )

ggplot(avg_growth, aes(x = day, y = mean_size, color = radiation_dose, linetype = chemo_regimen)) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci, fill = radiation_dose), alpha = 0.2, colour = NA) +
  labs(
    title = "Average Tumor Growth Curves by Treatment",
    x = "Day",
    y = "Tumor Size (mm)",
    color = "Radiation Dose",
    linetype = "Chemotherapy Regimen",
    fill = "Radiation Dose"
  ) +
  theme_minimal() +
  facet_grid(chemo_regimen ~ radiation_dose)