This DSPA Appendix shows examples of dynamic stochastic differential equation (SDE) models of cancer-treatment. The main objective is to develop an advanced generative stochastic differential equation (SDE) model proposing a new approach of cancer-treatment prediction. Ideally, the model should balance the relations between the primary pair of clinical outcome measures - intensity of the cancer-treatment intervention (by surgery, chemotherapy and radiation, manifesting as tumor local-control, \(LC\)) and side-effects (\(RP2\)), radiation pneumonitis grade 2 or higher, a measure of significant side effect in radiation therapy, see DOI: 10.1259/bjr.20220239.

Review the earlier Volterra-Lotka preditor-pray model and Hodgkin-Huxley (HH) model of neuronal neuron action potential propagation, which provide an oversimplified motivation for building a more advanced SDE cancer-model, which includes many controllable parameters accounting for the univariate distributions of dozens of observable clinical variables (biomarkers) and multivariate relations. The model (hyper)parameters can be tuned by Bayesian optimization using very sparsely sampled longitudinal clinical cancer-treatment data that includes the temporal dynamics of \(LC\) and \(RP2\) clinical outcomes from the point of initial diagnosis, through cancer-treatment intervention using alternative treatment strategies, to the final outcome stages. We’ll include a simulation of the proposed model and visualize SDE model instantiations that shows the time dynamics of \(LC\) and \(RP2\).

This advanced stochastic differential equation (SDE) model to predict and simulate cancer-treatment outcomes will track

Local Control (LC) of the tumor (a proxy for treatment intensity/efficacy),
Toxicity or side effects (RP2),
Additional clinical/biological variables (e.g., biomarkers, imaging signals, performance status), and
Patient survival or other endpoints.

Such a generative SDE framework can be viewed either as a Lotka-Volterra extension, where \(LC\) and \(RP2\) interact like predator-prey, or as an Hodgkin-Huxley-like system (with gating-type variables controlling “firing” or “flare-up” episodes of toxicity vs. tumor response). The key is to incorporate both deterministic dynamics (e.g., standard dose-response models) and stochastic components, e.g., random fluctuations reflecting inter-patient heterogeneity, measurement noise, or temporal variability in a patient’s biology.

The following R chunk includes knitr-specific parallel configuration specification.

library(knitr)
knitr::opts_chunk$set(
  cache = FALSE,        # Disable caching for parallel code
  autodep = FALSE,      # Disable automatic dependency detection
  message = FALSE,      # Suppress messages
  warning = FALSE       # Suppress warnings
)

1 First-Generation SDE Model

We will start with an oversimplified mathematical model of tumor dynamics without utilizing advanced biophysical principles. A more sophisticated biophysics SDE model will be presented later.

1.1 SDE Model Design

Let’s denote the state variables by \[\mathbf{X}(t) = \begin{pmatrix} LC(t) \\ RP2(t) \\ \mathbf{B}(t) \end{pmatrix},\] where \(LC(t)\) measures local tumor control (e.g., size reduction, imaging readouts), \(RP2(t)\) measures radiation- or drug-related side effects or toxicity, and \(\mathbf{B}(t)\) is a vector of relevant biomarkers or other continuous clinical variables (e.g., lab values, immunological markers), potentially of dimension 10-30 depending on data availability.

The treatment inputs can be denoted by \(\mathbf{u}(t)\) the treatment inputs (e.g., dose schedules for surgery/chemotherapy/radiation, or immunotherapy cycles). These external controls can drive the tumor and/or side-effect trajectories.

Next, we need to formulate the time evolution of the SDE model

\[d \mathbf{X}(t) = \mathbf{f}\bigl(\mathbf{X}(t), \mathbf{u}(t), \boldsymbol{\theta}\bigr)\,dt \;+\; \mathbf{g}\bigl(\mathbf{X}(t), \mathbf{u}(t), \boldsymbol{\theta}\bigr)\,d\mathbf{W}(t),\] where \(\mathbf{f}(\cdot)\) is the drift (deterministic) part, modeling how tumor control and toxicity progress on average, \(\mathbf{g}(\cdot)\) encodes the diffusion (noise) part, allowing for random fluctuations in responses and side effects, \(\boldsymbol{\theta}\) is a parameter set controlling tumor growth rates, clearance rates, synergy between multiple therapies, and so on, and \(\mathbf{W}(t)\) is a vector Wiener process (Brownian motion).

Using intuition from the Lotka-Volterra model, we may treat \(LC\) as a “prey” variable (grows unless suppressed) and \(RP2\) as a “predator” (side effects escalate and feed off continuing treatment). Let’s consider a possible drift in the form

\[\begin{aligned} dLC &= \Bigl[\alpha\,LC - \beta\,LC \cdot RP2 \;-\;\phi_{\text{decay}}\,LC \Bigr] dt \;+\; \sigma_{LC}(LC,RP2)\, dW_1,\\[6pt] dRP2 &= \Bigl[\gamma\,LC \cdot RP2 - \delta\,RP2 \;-\;\psi_{\text{treatment}}\, \mathbf{u}(t)\,RP2 \Bigr] dt \;+\; \sigma_{RP2}(LC,RP2)\, dW_2,\\[6pt] d\mathbf{B} &= \mathbf{F}_{\mathbf{B}}(\mathbf{B}, LC, RP2, \mathbf{u})\,dt \;+\; \boldsymbol{\Sigma}_{\mathbf{B}}(\mathbf{B},LC,RP2)\,d\mathbf{W}_{\mathbf{B}}, \end{aligned}\] where \(\alpha, \beta, \gamma, \delta, \dots\) are parameters, and \(\sigma_{LC}, \sigma_{RP2}, \boldsymbol{\Sigma}_{\mathbf{B}}\) are noise intensities (possibly state-dependent).

Similarly to the Hodgkin-Huxley model, \(LC\) is analogous to a “membrane voltage,” while \(RP2\) (and possibly some gating variables) evolves more slowly, controlling when toxicity “fires.” The drift terms would mimic the classic HH model formulation:

\[\begin{aligned} d LC &= \text{(nonlinear function of LC, gating vars, u(t))}\,dt \;+\; \ldots,\\ d RP2 &= \text{(slower gating dynamics, driven by LC transitions)}\,dt \;+\; \ldots \end{aligned}.\]

Such a SDE model captures bursts of toxicity or short-lived flare-ups.

The (hyper)parameter vector estimation can be achieved by using Bayesian optimization. Because clinical data are often limited and noisy, we can use Bayesian methods (e.g., approximate Bayesian computation in Julia, Stan, or specialized packages) or utilize Bayesian optimization to tune the parameter vector \(\boldsymbol{\theta}\).

In clinical setting, naturally sparse sampling (few time points per patient) is handled by fitting population-level priors, then refining posteriors for subgroups or individuals. We can define a likelihood linking the SDE outputs to observed data (LC, toxicities, biomarkers). Markov chain Monte Carlo (MCMC) or gradient-based Bayesian approaches can then infer the posterior distribution for \(\boldsymbol{\theta}\).

1.2 `R` Implementation

1.2.1 Simplified 2D SDE Modeling

We’ll start with a simpler first-generation SDE 2D model, using the snssde2d() function for the \(LC\)-\(RP2\) subsystem (as higher-dimensional SDEs require specialized handling in R). We’ll focus just on \(LC\) and \(RP2\) dynamics, with treatment as an external parameter.

The state variables are \[\mathbf{X}(t) = \begin{pmatrix} LC(t) \\ RP2(t) \end{pmatrix}. \]

The 2D SDE model system is \[\begin{aligned} dLC &= \left[\alpha LC \left(1 - \frac{LC}{K}\right) - \beta LC \cdot RP2 - \phi_{\text{treat}} u(t) \right] dt + \sigma_{LC} LC \, dW_1, \\ dRP2 &= \left[\gamma LC \cdot RP2 - \delta RP2 + \psi_{\text{treat}} u(t) \right] dt + \sigma_{RP2} \sqrt{RP2} \, dW_2. \end{aligned}\]

library(Sim.DiffProc)
library(plotly)

# Parameters
alphaP <- 0.3    # Tumor growth rate
KP <- 1.5        # Tumor carrying capacity
betaP <- 0.02    # Toxicity suppression rate
gammaP <- 0.015  # Toxicity activation rate
deltaP <- 0.4    # Toxicity clearance rate
phi_treat <- 0.1 # Treatment efficacy on LC
psi_treat <- 0.05 # Treatment-induced toxicity
sigma_LC <- 0.05 # Noise for LC
sigma_RP2 <- 0.1 # Noise for RP2

# Treatment schedule (example: weekly cycles)
u_t <- function(t) {
  ifelse(t %% 7 < 2, 1.0, 0.0)  # Treatment active 2 days/week
}

# Define drift and diffusion for snssde2d()
fx <- expression(
  alphaP * x * (1 - x/KP) - betaP * x * y - phi_treat * 1.0, # u(t) simplified to 1 when active
  gammaP * x * y - deltaP * y + psi_treat * 1.0
)

gx <- expression(
  sigma_LC * x,
  sigma_RP2 * sqrt(y)
)

# Simulate
set.seed(123)
mod <- snssde2d(
  drift = fx,
  diffusion = gx,
  x0 = c(x = 0.8, y = 0.1),
  t0 = 0,
  T = 60,
  N = 500,
  M = 50
)

# Visualize
fig <- plot_ly() 
for (i in 1:5) {  # Plot first 5 paths
  fig <- fig %>%
    add_trace(x = mod$time, y = mod$X[, i], name = paste("LC", i),
              type = "scatter", mode = "lines", line = list(color = "blue")) %>%
    add_trace(x = mod$time, y = mod$Y[, i], name = paste("RP2", i),
              line = list(color = "red"), yaxis = "y2")
}
fig <- fig %>% layout(
  title = "LC (Tumor Control) vs. RP2 (Toxicity)",
  yaxis2 = list(overlaying = "y", side = "right", title = "RP2"),
  xaxis = list(title = "Time (Days)")
)
fig

To match the existing R package capabilities, in this simplified 2D model the state space only covers \(LC\) and \(RP2\), but excludes the vector of biomarkers \(B\). We’ll revisit this later. The treatment schedule uses treatment cycles directly in the drift terms; phi_treat * 1.0 assumes full dose when active. For time-dependent u(t), we should precompute its values over the time grid and reference them in the drift expressions. An extension accounting for observable vectors of biomarkers, which are omitted in the 2D SDE model, may require simulating each biomarker as a 1D SDE or use alternative packages like pomp.

Below we show the protocol for R-based parameter estimation using rBayesianOptimization.

library(rBayesianOptimization)

# Simulate observed data (ground truth)
observed_data <- data.frame(
  time = seq(from=mod$t0, to=mod$T, by=mod$Dt),
  LC = rowMeans(mod$X),
  RP2 = rowMeans(mod$Y)
)

# Objective function to minimize MSE
objective_func <- function(alphaP, betaP, gammaP) {
  sim_mod <- snssde2d(
    drift = expression(
      alphaP * x * (1 - x/1.5) - betaP * x * y - 0.1,
      gammaP * x * y - 0.4 * y + 0.05
    ),
    diffusion = expression(0.05*x, 0.1*sqrt(y)),
    x0 = c(0.8, 0.1),
    t0 = 0,
    T = 60,
    N = 500
  )
  mse <- mean((sim_mod$X - observed_data$LC)^2 + (sim_mod$Y - observed_data$RP2)^2)
  return(list(Score = -mse))
}

# Bayesian optimization
opt_result <- BayesianOptimization(
  objective_func,
  bounds = list(
    alphaP = c(0.1, 0.5),
    betaP = c(0.01, 0.05),
    gammaP = c(0.01, 0.03)
  ),
  init_points = 5,
  n_iter = 20
)

## elapsed = 0.02   Round = 1   alphaP = 0.340964   betaP = 0.04482155  gammaP = 0.02975993 Value = -0.01000155 
## elapsed = 0.00   Round = 2   alphaP = 0.1091407  betaP = 0.01837854  gammaP = 0.02581208 Value = -0.01493705 
## elapsed = 0.02   Round = 3   alphaP = 0.428225   betaP = 0.03315608  gammaP = 0.01017334 Value = -0.005088571 
## elapsed = 0.00   Round = 4   alphaP = 0.1146278  betaP = 0.02775184  gammaP = 0.02165458 Value = -0.005338672 
## elapsed = 0.00   Round = 5   alphaP = 0.1940198  betaP = 0.0246884   gammaP = 0.01649646 Value = -0.007238786 
## elapsed = 0.02   Round = 6   alphaP = 0.4813248  betaP = 0.03072383  gammaP = 0.02149576 Value = -0.1544111 
## elapsed = 0.01   Round = 7   alphaP = 0.1550263  betaP = 0.03071254  gammaP = 0.02586016 Value = -0.009795527 
## elapsed = 0.00   Round = 8   alphaP = 0.1000 betaP = 0.0500  gammaP = 0.0100 Value = -0.04291667 
## elapsed = 0.00   Round = 9   alphaP = 0.3880017  betaP = 0.02188004  gammaP = 0.01164873 Value = -0.01688399 
## elapsed = 0.01   Round = 10  alphaP = 0.1963528  betaP = 0.01123368  gammaP = 0.0300 Value = -0.007565522 
## elapsed = 0.02   Round = 11  alphaP = 0.282257   betaP = 0.04799122  gammaP = 0.01532491 Value = -0.009646454 
## elapsed = 0.00   Round = 12  alphaP = 0.245288   betaP = 0.04974854  gammaP = 0.02067502 Value = -0.02056674 
## elapsed = 0.00   Round = 13  alphaP = 0.2989052  betaP = 0.0500  gammaP = 0.0100 Value = -0.01580561 
## elapsed = 0.00   Round = 14  alphaP = 0.1000 betaP = 0.0500  gammaP = 0.0300 Value = -0.02294954 
## elapsed = 0.00   Round = 15  alphaP = 0.1000 betaP = 0.0500  gammaP = 0.01689369 Value = -0.02679663 
## elapsed = 0.00   Round = 16  alphaP = 0.5000 betaP = 0.0100  gammaP = 0.0300 Value = -0.006871026 
## elapsed = 0.00   Round = 17  alphaP = 0.5000 betaP = 0.03723532  gammaP = 0.0100 Value = -0.02225222 
## elapsed = 0.00   Round = 18  alphaP = 0.1854054  betaP = 0.02100875  gammaP = 0.01888987 Value = -0.009016901 
## elapsed = 0.00   Round = 19  alphaP = 0.2725012  betaP = 0.02431722  gammaP = 0.0274043  Value = -0.01257695 
## elapsed = 0.00   Round = 20  alphaP = 0.2398215  betaP = 0.03820745  gammaP = 0.01359014 Value = -0.01843583 
## elapsed = 0.00   Round = 21  alphaP = 0.5000 betaP = 0.03124543  gammaP = 0.02835901 Value = -0.01552359 
## elapsed = 0.00   Round = 22  alphaP = 0.1661688  betaP = 0.03703741  gammaP = 0.02313721 Value = -0.01509109 
## elapsed = 0.00   Round = 23  alphaP = 0.4314345  betaP = 0.01653852  gammaP = 0.02955281 Value = -0.009896026 
## elapsed = 0.00   Round = 24  alphaP = 0.1721756  betaP = 0.02145591  gammaP = 0.02817582 Value = -0.01390439 
## elapsed = 0.00   Round = 25  alphaP = 0.1000 betaP = 0.02164099  gammaP = 0.01373757 Value = -0.01346684 
## 
##  Best Parameters Found: 
## Round = 3    alphaP = 0.428225   betaP = 0.03315608  gammaP = 0.01017334 Value = -0.005088571

print(opt_result$Best_Par)

##     alphaP      betaP     gammaP 
## 0.42822499 0.03315608 0.01017334

# Set the 3 parameters to the BO estimates and plot 50 SDE Model solutions
alphaP <- opt_result$Best_Par[[1]]    # Tumor growth rate
betaP  <- opt_result$Best_Par[[2]]    # Toxicity suppression rate
gammaP <- opt_result$Best_Par[[3]]  # Toxicity activation rate

# Define drift and diffusion for snssde2d()
fx <- expression(
  alphaP * x * (1 - x/KP) - betaP * x * y - phi_treat * 1.0, # u(t) simplified to 1 when active
  gammaP * x * y - deltaP * y + psi_treat * 1.0
)

gx <- expression(
  sigma_LC * x,
  sigma_RP2 * sqrt(y)
)

# Simulate
set.seed(123)
mod <- snssde2d(
  drift = fx,
  diffusion = gx,
  x0 = c(x = 0.8, y = 0.1),
  t0 = 0,
  T = 60,
  N = 500,
  M = 50
)

# Visualize
fig <- plot_ly() 
for (i in 6:10){    # 1:mod$M) {  # Plot just 5 paths or all of them
  fig <- fig %>%
    add_trace(x = mod$time, y = mod$X[, i], name = paste("LC", i),
              type = "scatter", mode = "lines", line = list(color = "blue")) %>%
    add_trace(x = mod$time, y = mod$Y[, i], name = paste("RP2", i),
              line = list(color = "red"), yaxis = "y2")
}
fig <- fig %>% layout(
  title = "LC (Tumor Control) vs. RP2 (Toxicity)",
  yaxis2 = list(overlaying = "y", side = "right", title = "RP2"),
  xaxis = list(title = "Time (Days)")
)
fig

DSPA2: Data Science and Predictive Analytics (UMich HS650)

Appendix 19: Cancer Treatment SDE-Modeling

SOCR (Ivo Dinov and SOCR)

September 16, 2025

1 First-Generation SDE Model

1.1 SDE Model Design

1.2 `R` Implementation

1.2.1 Simplified 2D SDE Modeling

DSPA2: Data Science and Predictive Analytics (UMich HS650)

Appendix 19: Cancer Treatment SDE-Modeling

SOCR (Ivo Dinov and SOCR)

September 16, 2025

1 First-Generation SDE Model

1.1 SDE Model Design

1.2 R Implementation

1.2.1 Simplified 2D SDE Modeling

1.2 `R` Implementation