Surgery, Chemotherapy, and Radiation as a Combinatorial Optimization Problem
When a patient is diagnosed with a tumor that requires multiple treatment modalities, doctors don’t just choose which treatments to use — they also have to choose the order. Should surgery come first, followed by radiation to clean up the margins? Or should chemotherapy shrink the tumor first, making surgery safer and less invasive? Each ordering carries a different balance of risk, healing time, and toxicity.
This is, at its core, a sequencing optimization problem — mathematically very close to the famous Traveling Salesman Problem (TSP), except instead of minimizing travel distance between cities, we minimize a “treatment cost” that depends on which therapy follows which.
In this article, I’ll build a concrete, solvable example with three treatments — Surgery (S), Chemotherapy (C), and Radiation (R) — model the cost of every possible ordering, find the optimal sequence with Python, visualize it (including in 3D), and then scale the problem up to nine treatment modalities to show why a smarter algorithm (dynamic programming) becomes essential as the number of choices grows.
Note: the numbers used below are illustrative, synthetic values chosen to make the optimization clear — not real clinical data. The point of this article is the algorithm, not medical advice.
Formulating the Problem
Let treatments be indexed $1, 2, \dots, n$. We define two cost components:
- $B(i)$: the baseline cost of starting the entire treatment plan with treatment $i$ (e.g., the inherent risk of operating on an untreated tumor).
- $C(i, j)$: the transition cost of performing treatment $j$ immediately after treatment $i$ (e.g., the risk of operating on tissue that was just irradiated).
For a treatment order $\pi = (\pi_1, \pi_2, \dots, \pi_n)$, the total cost is:
$$
\text{TotalCost}(\pi) = B(\pi_1) + \sum_{k=1}^{n-1} C(\pi_k, \pi_{k+1})
$$
We want to find the permutation $\pi^*$ that minimizes this:
$$
\pi^* = \arg\min_{\pi \in S_n} \text{TotalCost}(\pi)
$$
This is exactly the structure of the shortest Hamiltonian path problem. For small $n$, we can simply check every permutation (there are $n!$ of them). For larger $n$, we’ll use the Held-Karp dynamic programming algorithm, which solves it in:
$$
O(n^2 \cdot 2^n) \quad \text{instead of} \quad O(n!)
$$
The DP recurrence builds up the optimal cost of visiting a subset of treatments $S$ and ending at treatment $j$:
$$
dp[S][j] = \min_{i \in S \setminus {j}} \Big( dp[S \setminus {j}][i] + C(i,j) \Big), \qquad dp[{j}][j] = B(j)
$$
The Python Implementation
1 | import numpy as np |
Execution Results
All possible treatment sequences and their total cost: Surgery -> Chemotherapy -> Radiation cost = 9.50 Surgery -> Radiation -> Chemotherapy cost = 14.00 Chemotherapy -> Surgery -> Radiation cost = 10.50 Chemotherapy -> Radiation -> Surgery cost = 9.50 Radiation -> Surgery -> Chemotherapy cost = 10.00 Radiation -> Chemotherapy -> Surgery cost = 8.50 <== OPTIMAL Optimal treatment order: Radiation -> Chemotherapy -> Surgery Minimum total cost: 8.5
Brute force : best cost = 19.502, time = 0.823 s Held-Karp DP : best cost = 19.502, time = 0.0095 s Results match : True Speed-up : 86.4x faster with dynamic programming Optimal 9-step order (DP): Chemotherapy -> Photodynamic Therapy -> Radiation -> Surgery -> Brachytherapy -> Immunotherapy -> Hormone Therapy -> Cryotherapy -> Targeted Therapy
Walking Through the Code
Section 1 — Defining the problem. B3 holds the baseline cost of starting with each treatment, and C3 is a 3×3 matrix of transition costs. Row $i$, column $j$ of C3 represents “the extra cost of doing treatment $j$ right after treatment $i$.” For instance, C3[0][2] = 6.0 (Surgery → Radiation) is deliberately set high, modeling the realistic concern that irradiating a fresh surgical wound increases complication risk. C3[1][0] = 1.5 (Chemo → Surgery) is low, reflecting that a tumor downstaged by chemotherapy is often easier and safer to remove.
total_cost() simply walks through a permutation and sums the baseline cost plus every transition cost along the path — a direct implementation of the formula above.
evaluate_all_permutations() uses Python’s itertools.permutations to brute-force all $3! = 6$ orderings, which is trivial at this scale, and tracks the minimum.
Section 2 — Bar chart. This plots the total cost of all six possible orderings side by side, with the optimal one highlighted in red, making it visually obvious which choice wins and by how much.
Section 3 — 3D transition matrix. Using ax.bar3d, we render the transition cost matrix as a literal 3D bar chart: the x-axis is the “from” treatment, the y-axis is the “to” treatment, and the height (z-axis) is the cost of that transition. This makes asymmetries immediately visible — for example, you can see that the Surgery → Radiation bar towers over Chemotherapy → Surgery.
Section 4 — 3D cumulative cost trajectories. Each of the six possible sequences is drawn as a 3D line: the x-axis is the treatment step (1st, 2nd, 3rd), the y-axis separates the six sequences into “lanes,” and the z-axis shows the cumulative cost as it builds up step by step. The optimal sequence is drawn thicker and in red, so you can see exactly where it pulls ahead of the alternatives.
Section 5 — Scaling to 9 treatments. A real hospital might consider far more than three modalities — immunotherapy, hormone therapy, targeted therapy, and so on. With 9 treatments, brute force must check $9! = 362{,}880$ orderings. brute_force_best_only() is a leaner version that discards results as soon as the running cost exceeds the current best (a simple branch-and-bound pruning trick), avoiding the memory overhead of storing all permutations.
held_karp() is the optimized algorithm. It uses a bitmask mask to represent “the set of treatments already scheduled,” and dp[mask][j] stores the minimum cost of a sequence that uses exactly the treatments in mask and ends at treatment j. By building this table up from single-treatment subsets to the full set, it explores only $O(n^2 2^n)$ states instead of $O(n!)$ permutations — for $n=9$ that’s roughly 41,000 operations versus 362,880 permutations (each requiring multiple additions), a dramatic reduction. The parent table lets us reconstruct the actual optimal path afterward, not just its cost.
Section 6 — Runtime comparison. Finally, we time both approaches directly and plot the result, along with a sanity check (math.isclose) confirming both methods agree on the optimal cost — proof that the dynamic programming version isn’t an approximation, just a faster exact algorithm.
Interpreting the Results
In the 3-treatment example, the cost structure favors Radiation → Chemotherapy → Surgery as the lowest-cost path. This pattern — delivering non-surgical therapy before the operation — mirrors a real and well-established clinical strategy: neoadjuvant treatment, used for example in locally advanced rectal cancer, where shrinking the tumor first can make surgery less extensive and more effective. Our toy model arrives at a structurally similar conclusion purely from the cost numbers, which is a nice demonstration of how combinatorial optimization can echo real-world clinical reasoning — though actual treatment decisions are always made by multidisciplinary medical teams based on clinical guidelines and individual patient factors, not a cost matrix.
The 9-treatment scaling experiment is the real takeaway for anyone building decision-support tools: as soon as the number of options grows past a handful, brute-force enumeration becomes computationally infeasible ($15! \approx 1.3$ trillion), while the Held-Karp dynamic programming approach remains practical up to roughly 20–25 items, since its complexity grows exponentially in $n$ but only polynomially in the cost-table size — a far gentler curve than factorial growth.
