A Two-Compartment Pharmacokinetic Approach
Why particle size and infusion rate matter
Nanoparticle-based cancer therapeutics rely on the Enhanced Permeability and Retention (EPR) effect: tumor blood vessels are leaky and poorly drained, so nanoparticles of the “right” size accumulate preferentially inside tumor tissue. But size cuts both ways. Particles that are too small are filtered out by the kidneys before they can accumulate anywhere. Particles that are too large are swept up by the reticuloendothelial system (liver, spleen) and cleared from circulation quickly. There is a sweet spot, typically somewhere in the 80–150 nm range, where a particle circulates long enough and permeates the tumor efficiently enough to deliver a meaningful dose.
At the same time, the infusion rate is not a free parameter — push too much drug into the blood and you hit systemic toxicity limits before the tumor ever gets a useful concentration.
This post builds a small, fully quantitative model that captures both effects, then numerically finds the particle size and infusion rate that maximize tumor drug exposure without violating a plasma toxicity constraint.
The model
We treat the body as two compartments: plasma and tumor tissue. Drug crosses between them at a rate set by the permeability–surface area product $PS$, which depends on particle diameter $d$ through the EPR effect:
$$PS(d) = PS_{max}\exp!\left(-\frac{(d-d_{opt})^2}{\sigma_d^2}\right)$$
Clearance from plasma also depends on particle size — small particles are cleared fast by renal filtration, large ones by the RES, and clearance is minimized near some optimal circulation diameter $d_c$:
$$CL(d) = CL_0\left(1+\left(\frac{d-d_c}{w_c}\right)^2\right)$$
With these two size-dependent parameters, the plasma concentration $C_p$ and tumor concentration $C_t$ evolve as:
$$\frac{dC_p}{dt} = \frac{R_{in}}{V_p} - \frac{CL(d)+PS(d)}{V_p}C_p + \frac{PS(d)}{V_p}C_t$$
$$\frac{dC_t}{dt} = \frac{PS(d)}{V_t}\left(C_p-C_t\right) - k_{deg}C_t$$
where $R_{in}$ is the (zero-order) infusion rate, $V_p$ and $V_t$ are the plasma and tumor volumes, and $k_{deg}$ is the tumor drug elimination rate.
This is a linear system $\dot{\mathbf{x}} = A\mathbf{x} + \mathbf{b}$, so it has a closed-form solution — no numerical ODE integration needed at all. We use this to evaluate the whole model for a full grid of particle sizes at once, using nothing but NumPy array broadcasting.
The metric we optimize is the tumor drug exposure, i.e. the area under the tumor concentration curve:
$$AUC_{tumor} = \int_0^{T_{end}} C_t(t),dt$$
subject to a plasma toxicity ceiling:
$$\max_{d,,R_{in}}\ AUC_{tumor}(d,R_{in}) \quad \text{s.t.} \quad \max_t C_p(t) \le C_{tox}$$
Why we can avoid a slow parameter sweep
A naive implementation would loop over every (particle size, infusion rate) pair — say 300 × 200 = 60,000 combinations — and call scipy.integrate.solve_ivp for each one. That’s 60,000 separate numerical integrations, which is slow and wasteful.
We avoid this in two ways:
- Analytical matrix exponential. For a 2×2 system, $e^{At}$ has a closed form in terms of the eigenvalues of $A$ (Sylvester’s formula), so we never call a generic matrix-exponential routine — we compute it directly with array operations, broadcast simultaneously over every particle size and every time point.
- Linearity in the infusion rate. Because the system is linear, the entire trajectory scales linearly with $R_{in}$. We only need to solve the model once per particle size at a unit infusion rate; every other infusion rate is just a multiplication. This collapses a genuinely 2D sweep into a 1D computation plus simple algebra.
The result: the full parameter sweep over 300 particle sizes and 200 infusion rates runs in a fraction of a second, with zero explicit loops over the parameter grid.
Source code
1 | import numpy as np |
Code walkthrough
Parameters (section 1): All values are illustrative, chosen to produce realistic-looking pharmacokinetic behavior rather than sourced from a specific real drug. Vt is deliberately tiny (10 mL) compared to Vp (5 L), which is what makes drug exchange with the tumor compartment fast relative to systemic elimination — a common feature of small, well-perfused tumor lesions.
Size-dependent physiology (section 2): PS_of_d encodes the EPR effect as a Gaussian peaking at 100 nm — permeability collapses quickly for particles that are too small or too large. CL_of_d encodes the trade-off on the clearance side: a parabola with a minimum at 130 nm, reflecting that very small particles are lost through renal filtration and very large ones are captured by the liver/spleen. These two peaks don’t coincide, which is exactly what creates a genuine, non-trivial optimum — not just an assumption stated up front.
Closed-form linear solver (section 3): Instead of stepping the ODE forward in time with an integrator, eig_2x2 and expm_At compute the matrix exponential $e^{At}$ analytically via Sylvester’s formula for a 2×2 matrix. All of the array shapes carry a leading “particle size” axis, so expm_At returns the exponential for all 300 particle sizes and every time point in one shot — this is the core trick that keeps everything fast. simulate_unit_response then assembles the full concentration trajectories for both the infusion phase and the post-infusion decay phase, again for every particle size simultaneously.
The sweep (section 4): We only ever solve the model once, at a unit infusion rate. Because the whole system is linear, the response to any infusion rate is just that unit response scaled by R0. This lets us derive the toxicity-constrained optimal infusion rate for every particle size analytically: R0_star_d = C_tox / Cp_max_unit, with no additional simulation required. Maximizing AUC_star_d over the particle-size axis then gives the overall optimum.
Figures 3 and 4 still sweep a full 200-point grid of infusion rates for the visualization — this is just an outer product (np.outer), essentially free computationally, and it lets the reader see the entire feasible operating region, not just the single optimal point.
Reading the results
Running the script gives:
1 | Optimal particle diameter : 127.4 nm |
The optimal size sits between the pure permeability peak (100 nm) and the pure circulation-time peak (130 nm) — the model is genuinely balancing two competing effects rather than just reproducing one of the input assumptions. About 82% of the sampled (size, dose) grid turns out to be within the toxicity constraint, but only a narrow ridge along that constraint boundary actually maximizes tumor exposure — which is exactly what Figures 3 and 4 are meant to make visible.
Figure 1 shows why 127 nm is the answer: permeability is already past its peak and declining, but clearance is still improving (still falling), and the net effect nets out slightly past the permeability optimum.
Figure 2 shows the actual concentration-time trace at the optimum — plasma concentration rises during the 2-hour infusion and decays afterward, while tumor concentration lags slightly behind and clears more slowly, since it’s driven by the smaller kdeg rate rather than plasma clearance.
Figure 3 is the main result: a 3D surface of tumor AUC over the full (particle size, infusion rate) plane, with everything above the toxicity ceiling masked out (NaN) so only physically realizable treatment plans are shown. The red marker sits on the ridge of the feasible region where AUC is maximized.
Figure 4 is the same surface from directly above, as a contour plot, with the toxicity boundary drawn explicitly in red — this is the easiest view for reading off exact coordinates of the optimum.
Figure 5 collapses the whole problem down to a single curve: for each particle size, what is the best tumor exposure achievable without crossing the toxicity limit? This makes the interior optimum at ~127 nm unambiguous.
Takeaway
Framing drug delivery as a small linear compartmental model, rather than reaching for a general-purpose ODE solver, buys two things at once: a rigorous, closed-form answer to “what particle size and dose are optimal under a safety constraint,” and a computation cheap enough to sweep and visualize the entire design space rather than just reporting a single number.
