Mathematical oncology · Preprint

Some tumor cells survive chemotherapy
because the drug never arrives.

A patient-specific reaction-diffusion model of triple-negative breast cancer, proving that spatial drug-insufficient zones exist before treatment even begins, and locating them directly from a baseline MRI scan.

Koushik Chowdhury Universität des Saarlandes ORCID 0009-0007-6319-3769
∂Ω Ω_h

Drug enters at the tumor boundary ∂Ω and decays over a characteristic length L_d = √(D_d/γ) ≈ 0.29 mm. Every tumor in this cohort was larger than that, so a core region Ω_h never saw a therapeutic dose.

46.6%
of mean tumor volume received sub-therapeutic drug under standard dosing
0.29mm
computed drug penetration depth, matching published doxorubicin measurements
15
TNBC patients, I-SPY2 trial, reconstructed from raw DICOM to 3D tumor geometry
03
theorems proving failure zones must exist, and when they can be dosed away

Abstract

The physics of the tumor, not just its biology, decides who responds.

Fewer than three in ten TNBC patients see a complete response to standard chemotherapy. The usual explanation is molecular resistance. This paper asks a simpler, geometric question first: does the drug even reach the cell?

The model

Three coupled fields, tumor cell density, oxygen, and drug concentration, evolve by reaction-diffusion over the patient's actual tumor shape, reconstructed from MRI. A failure zone is defined as any point where oxygen and drug both stay below their effective thresholds at once.

The claim

Given a tumor's size and shape alone, before a single dose is given, it is possible to prove whether a drug-insufficient region must appear, how large it will get, and the dosing frequency that would make it disappear.

Why geometry matters

A pharmacological sanctuary, made precise.

Experimentally, cytotoxic drugs are known to penetrate only a few cell layers past the nearest vessel before their concentration drops too low to matter. This paper turns that observation into a length scale that falls directly out of the drug's own diffusion equation.

At steady state, drug supplied at the tumor boundary decays over a characteristic distance L_d = √(D_d / γ), where D_d is the drug's diffusivity through tissue and γ is its tissue clearance rate. Any point further than L_d · log(d_peak / d_thr) from the boundary will sit below the therapeutic threshold no matter how much drug is given at the surface.

With baseline doxorubicin parameters, that works out to roughly 0.29 mm. The smallest tumor in this cohort had a radius of about 8 mm. Every patient's geometry, by construction, contained a region the drug could not reach.

∂Ω tumor core d(x) d_thr L_d ≈ 0.29 mm Ω_h — below threshold

Drug concentration decays exponentially from the tumor boundary. Past L_d, concentration is already below d_thr — the shaded region receives no effective dose regardless of the dose given at the surface.

Main results

Three theorems, one on existence, one on why, one on what to do about it.

Each proof sketch reduces to a physical argument: a maximum principle, a comparison of sub- and supersolutions, and a contraction on the dosing period.

01

Well-posedness on patient-specific geometry

The coupled system has a unique, bounded, non-negative weak solution on any bounded Lipschitz domain, which is exactly the regularity that an MRI-reconstructed tumor surface actually has.

→ the model is mathematically sound on real, irregular tumor shapes, not just idealized spheres.
02

Existence and persistence of failure zones

If a tumor contains a point far enough from its boundary, and its local cell density is high enough to out-consume oxygen diffusion, then a small ball around that point becomes and stays a failure zone.

|Ω_h(t)| ≥ ω₃ρ³ for all t ≥ t*

Compact, convex, low-density tumors are unlikely to satisfy this. Large, irregular tumors with necrotic cores are much more likely to.

03

A dosing frequency that erases the failure zone

Because drug supply is periodic, its long-run concentration converges to a fixed periodic orbit. There is a critical period P* below which that orbit stays above threshold everywhere, and above which a failure zone reappears and grows with P − P*.

P* = (1/γ)·log(d_peak/d_thr) − (L_d/γ)·log(d_peak/d_thr)·dist(x*,∂Ω)/L_d
→ P* is computable from one baseline MRI and the drug's own pharmacokinetics, before treatment starts.

From scanner to simulation

A full pipeline, from raw DICOM to a failure zone map.

STEP 1

Cohort selection

360 TNBC patients filtered from I-SPY2 by hormone receptor and HER2 status. 16 had complete baseline DCE-MRI available in TCIA; one produced no usable mesh, leaving 15.

STEP 2

Mask extraction

Tumor label identified from the multi-label segmentation DICOM by smallest non-background voxel count. Computed volumes matched the trial's radiomic volumes within 15%.

STEP 3

Mesh reconstruction

Marching cubes at isosurface 0.5, then Laplacian smoothing, turns each binary mask into the 3D tumor surface the PDE solver runs on.

STEP 4

PDE simulation

Finite-difference solve of all three fields over an 84-day, four-cycle AC-T schedule, on a grid matched to the scan's own voxel spacing.

DCE-MRI baseline scan for patient ISPY2-164468 with tumor mask overlay

Fig. 1 — Baseline DCE-MRI for patient ISPY2-164468. Left: raw scan. Centre: extracted binary tumor mask. Right: mask overlaid, confirming it lines up with the contrast-enhancing tissue.

Why the shape itself is the point

A sphere is the best possible shape for drug delivery: every interior point is as close to the surface as geometry allows. Real TNBC tumors are not spheres. They branch, they have thin connected clusters, they have deep interior folds where the distance to the nearest supplying boundary is far larger than the tumor's nominal radius would suggest.

That is the entire reason a patient-specific geometry, not a generic spherical approximation, is required to compute where a failure zone will actually sit.

Reconstructed 3D tumor surface mesh for patient ISPY2-164468

Fig. 2 — Reconstructed 3D surface for the same tumor. The irregular, non-convex, branching shape is typical of the cohort, and it is exactly this kind of geometry Theorem 1 is built to handle.

What the solver actually resolves

Degree-two finite elements resolve the curved oxygen and drug profiles near the boundary with far fewer nodes than a linear mesh would need. A semi-implicit Crank-Nicolson scheme keeps the nonlinear kill and consumption terms stable, and an adaptive time step tightens automatically around each of the four dosing days.

Simulation output

Watching the failure zone appear.

For patient ISPY2-164468, the failure zone was empty at day 0. It grew to 1,096 voxels by day 21, to 9,394 by day 42, and to 14,211 by day 84, the persistence predicted by Theorem 2 playing out voxel by voxel.

By day 84 the mean oxygen concentration inside the tumor mask had fallen to 0.68 mmHg, far under the 5 mmHg hypoxic threshold, while tissue just outside the tumor stayed near the 40 mmHg boundary value. That contrast is driven by how much oxygen a dense tumor cell population consumes, not by anything about the drug at all.

Sensitivity testing later in the paper shows this is the dominant effect: hypoxia, not drug transport, decides where the failure zone sits.

PDE simulation output at day 84 showing oxygen field, drug field, and computed failure zone

Fig. 3 — Day 84, middle axial slice. Left: oxygen field (red = deeply hypoxic). Centre: drug field, showing the diffusion gradient from the boundary inward. Right: the computed failure zone in white.

Cohort results

All 15 patients. No exceptions.

Failure zone fraction ranged from 0.378 to 0.527, mean 0.466. Every tumor in this cohort, complete responder or not, carried a substantial region no dose could reach under the standard three-weekly schedule.

Patient IDTumor vol. (cc)FZ vol. (cc)FZ fractionpCR
ISPY2-16446831.9915.460.4830
ISPY2-19602417.247.920.4590
ISPY2-2083036.242.550.4091
ISPY2-25274818.068.970.4971
ISPY2-48950430.3915.610.5141
ISPY2-5357792.751.040.3780
ISPY2-56423413.416.560.4891
ISPY2-57897533.9017.870.5271
ISPY2-62231538.0015.180.4000
ISPY2-6258542.141.040.4871
ISPY2-69709824.6211.050.4491
ISPY2-72780424.7710.960.4421
ISPY2-8294915.452.760.5060
ISPY2-91070640.4619.650.4861
ISPY2-9349065.822.710.4650

pCR = pathological complete response. Failure zone fraction alone does not separate responders from non-responders, which is expected: response also depends on molecular factors this physical model does not capture.

Scatter plot of failure zone volume against total tumor volume for 15 patients, coloured by treatment response

Fig. 4 — Failure zone volume scales almost linearly with total tumor volume, regardless of eventual response.

Reading the scatter

The near-linear relationship is itself informative. It says the failure zone fraction is close to a geometric constant for this drug and this disease, set mostly by the mismatch between tumor size and penetration depth, rather than by anything specific to an individual patient's biology.

Bar chart showing change in failure zone fraction when drug diffusivity, decay rate, and proliferation rate are each varied by 20 percent

Fig. 5 — Varying drug diffusivity or decay rate by 20% moves the failure zone fraction by at most 0.02%. Varying proliferation rate moves it by up to 0.37%.

Sensitivity analysis

It is the oxygen, not the drug.

By day 84, oxygen inside the tumor mask has already collapsed to 0.68 mmHg, far under the 5 mmHg hypoxic threshold. Whether the drug diffuses slightly faster or slower barely changes which voxels are hypoxic, because hypoxia is already the binding constraint almost everywhere.

The practical consequence: a drug engineered to diffuse further or decay slower would not, on its own, close this failure zone. Something that addresses the hypoxia directly, a hypoxia-activated prodrug, or a vascular-normalizing agent given ahead of chemotherapy, is more likely to help.

Does the failure zone predict outcome?

A proof of concept

Leave-one-out cross validation on 15 patients, random forest and logistic regression, three feature groups: radiomics, PDE-derived features, and the two combined.

Radiomics only

0.741 Random forest AUC
0.815 Logistic regression AUC

PDE features only

0.278 Random forest AUC
0.333 Logistic regression AUC

Combined

0.759 Random forest AUC
0.796 Logistic regression AUC

Radiomics dominate on this cohort. Adding PDE features nudges random forest AUC up by 0.018 and logistic regression AUC down by 0.019. Fifteen patients is not enough to call this a real effect either way, but the failure zone fraction still ranks third by feature importance in the ensemble model.

ROC curves for the combined feature set classifier

Fig. 6 — ROC curves, combined feature set, leave-one-out cross validation.

Feature importance ranking from the random forest classifier

Fig. 7 — Random forest feature importance. Failure zone fraction (fz_fraction) sits third, behind early functional tumor volume change and baseline functional tumor volume.

What this could change, and what it can't yet claim

Three questions a clinician could act on.

Limitations: single homogeneous cell population, fixed tumor geometry over the full 84-day course, literature-derived rather than per-patient transport parameters, and a 15-patient cohort too small to validate any clinical claim on its own. Parameters were calibrated to doxorubicin and may need recalibration for paclitaxel or carboplatin.