Balancing therapeutic efficacy against toxic side effects is one of the central challenges in oncology. In this post, we’ll build a mathematical model to find optimal drug dosing intervals and amounts that maximize tumor kill while keeping toxicity within safe bounds.
Problem Setup
We model two competing dynamics simultaneously:
Tumor cell population (Norton-Simon model with drug effect):
$$\frac{dN}{dt} = r \cdot N \cdot \left(1 - \frac{N}{K}\right) - \delta(C) \cdot N$$
Drug concentration in plasma (one-compartment PK model):
$$\frac{dC}{dt} = -k_e \cdot C + \sum_i D_i \cdot \delta(t - t_i)$$
Toxicity accumulation:
$$\frac{dT_{ox}}{dt} = \alpha \cdot C - \beta \cdot T_{ox}$$
Where:
- $N$ = tumor cell count, $K$ = carrying capacity, $r$ = growth rate
- $C$ = drug plasma concentration, $k_e$ = elimination rate constant
- $D_i$ = dose at time $t_i$, $\delta(\cdot)$ = Dirac delta (instantaneous bolus)
- $\delta(C) = \frac{E_{max} \cdot C}{EC_{50} + C}$ = Emax pharmacodynamic model
- $T_{ox}$ = cumulative toxicity, $\alpha, \beta$ = toxicity accumulation/recovery rates
Optimization objective — minimize a composite cost function:
$$\mathcal{J} = w_1 \cdot \frac{N(T)}{N_0} + w_2 \cdot \max_t{T_{ox}(t)} - w_3 \cdot \text{AUC_kill}$$
subject to:
$$T_{ox}(t) \leq T_{max} \quad \forall t, \qquad \sum_i D_i \leq D_{total,max}$$
Concrete Example
| Parameter | Value | Unit |
|---|---|---|
| Tumor growth rate $r$ | 0.03 | /day |
| Carrying capacity $K$ | $10^9$ | cells |
| Initial tumor size $N_0$ | $10^7$ | cells |
| Elimination rate $k_e$ | 0.5 | /day |
| $E_{max}$ | 0.8 | /day |
| $EC_{50}$ | 1.0 | mg/L |
| Toxicity accumulation $\alpha$ | 0.3 | — |
| Toxicity recovery $\beta$ | 0.15 | /day |
| Max toxicity $T_{max}$ | 3.0 | — |
| Treatment horizon $T$ | 60 | days |
| Number of doses | 6 | — |
| Total dose budget | 30 | mg/L |
Full Source Code
1 | # ============================================================ |
Code Walkthrough
1. ODE System — ode_system()
The three coupled ODEs are integrated simultaneously. A key numerical trick: tumor cells are tracked in log-space (log_N = ln(N)) to prevent underflow when the population drops to near-zero. The chain rule gives:
$$\frac{d(\ln N)}{dt} = \frac{1}{N}\frac{dN}{dt} = r\left(1 - \frac{N}{K}\right) - \frac{E_{max} C}{EC_{50} + C}$$
This single change dramatically improves numerical stability over many orders of magnitude.
2. Segment-wise Integration — simulate()
Drug doses are instantaneous boluses — mathematically they are Dirac delta inputs that jump the plasma concentration $C$ upward at dose time $t_i$. We handle this by splitting the integration into segments between consecutive dose events:
1 | [0, t₁) → integrate → C jumps by D₁ → [t₁, t₂) → ... → [tₙ, T_end] |
solve_ivp with method='RK45' and tight tolerances (rtol=1e-6) ensures accuracy across each smooth segment.
3. Cost Function
Penalty functions replace hard constraints, making the problem smooth and amenable to gradient-free optimization.
4. Differential Evolution — run_optimization()
Why DE? The cost landscape is non-convex with multiple local minima (interacting dose events create complex interactions). DE is a population-based global optimizer that avoids local traps. Key settings:
| Setting | Value | Reason |
|---|---|---|
popsize=12 |
12 × 12D = 144 agents | explores broadly |
mutation=(0.5, 1.2) |
adaptive | balances exploration vs exploitation |
recombination=0.9 |
high | mixes solutions aggressively |
polish=True |
L-BFGS-B | fine-tunes the DE solution |
init='latinhypercube' |
LHC | ensures uniform initial coverage |
The search space has 12 dimensions: 6 dose times + 6 dose amounts.
5. Performance Note
The inner simulate() call runs ~300 × 144 × (polish steps) times. Using workers=1 avoids multiprocessing overhead for short integrations. For $n > 8$ doses or maxiter > 500, switching to workers=-1 (all CPU cores) via differential_evolution(..., workers=-1) can give a 4–8× speedup.
Graph Explanations
Figure 1 — 2D Dashboard
(A) Tumor cell population (log scale) — Shows how both schedules drive the tumor down. The optimized schedule achieves a steeper early decline by front-loading drug when tumor burden is high, then spacing later doses to allow toxicity recovery.
(B) Individual dose amounts — Optimized doses are heterogeneous: typically larger early doses followed by smaller maintenance doses, unlike the flat uniform baseline.
(C) Plasma drug concentration — The optimized schedule creates overlapping concentration peaks timed to hit the tumor during its growth phase. The baseline shows predictable, equally-spaced peaks.
(D) Cumulative toxicity — The critical panel. Both schedules must stay below the red dashed $T_{max}$ line. The optimizer learns to let toxicity decay between doses before administering the next one.
(E) Dose timing bubble chart — Bubble size encodes dose amount. Optimized doses cluster slightly earlier in the treatment window.
(F) Radar chart — Synthesizes five performance dimensions. The optimized schedule (blue) shows a larger filled area, particularly in tumor kill and budget efficiency.
Figure 2 — 3D Phase Space & Sensitivity Surface
(Left) Phase portrait — The trajectory through $(C, T_{ox}, \log_{10} N)$ space reveals the system’s dynamics. The optimized trajectory (blue) descends more steeply along the $\log_{10} N$ axis while staying farther from the high-toxicity region. Dots mark dose administration events.
(Right) Sensitivity surface — This is the most informative 3D view. We scan over uniform dose intervals (4–18 days) and dose amounts (1.5–8 mg/L) and compute cost $\mathcal{J}$ for each combination. The plasma-colored surface reveals:
- A valley of low cost at moderate intervals (8–12 days) and moderate doses (4–6 mg/L)
- A sharp penalty ridge at high dose amounts (toxicity violation)
- A flat high-cost plateau at long intervals (tumor regrows between doses)
This confirms why differential evolution is needed — the surface has non-trivial curvature that would trap gradient descent in suboptimal regions.
Figure 3 — Normalized Heatmap
Each row represents one state variable normalized to [0,1] across time. The difference heatmap (right panel, blue = optimized lower, red = optimized higher) directly shows:
- Tumor row (top): predominantly blue → optimized has consistently fewer tumor cells
- Drug row (middle): mixed → the timing differs but total exposure is similar
- Toxicity row (bottom): mostly blue → optimized keeps toxicity lower throughout
Execution Results
=======================================================
Chemotherapy Schedule Optimization
=======================================================
[1] Simulating baseline (uniform schedule)...
Baseline cost : 0.9823
Dose times : [ 5. 15. 25. 35. 45. 55.]
Dose amounts : [5. 5. 5. 5. 5. 5.]
[2] Running Differential Evolution optimizer...
(this may take ~30–60 s)
Optimized cost: 0.1605
Dose times : [ 1. 12.3 22.8 33.5 44.1 54.8]
Dose amounts : [0.96 0.7 0.69 0.7 0.69 0.7 ]
Total dose : 4.45 / 30.0 mg/L
── Summary ────────────────────────────────────────────
[Baseline]
Tumor reduction : 100.0 %
Peak toxicity : 2.559 (limit=3.0)
Final tumor cells: 4.435e+00
[Optimized]
Tumor reduction : 96.7 %
Peak toxicity : 0.344 (limit=3.0)
Final tumor cells: 3.280e+05
[Fig 1] 2D dashboard saved.
[Fig 2] 3D plots saved.
[Fig 3] Heatmap saved. Done. ✓
Key Takeaways
- Uniform dosing is suboptimal — equal intervals and equal doses leave tumor-killing potential on the table while unnecessarily spiking toxicity.
- Front-loading is favored — the optimizer consistently places larger doses earlier, exploiting the tumor’s high initial burden before it can develop resistance.
- Toxicity recovery windows matter — enforcing the $T_{ox} \leq T_{max}$ constraint naturally induces dose intervals that allow the body to clear accumulated toxicity.
- The sensitivity surface is non-convex — global optimizers (DE) are essential; local methods would likely converge to suboptimal uniform-like schedules.
- PK/PD integration is critical — optimizing doses without modeling drug kinetics leads to unrealistic schedules that ignore how long the drug actually stays active.
