Automated Design for Polypills

• The paper introduces an automated framework that co-optimizes pill geometry and material layout to achieve precise drug-release profiles.
• It employs a supershape parameterization combined with a neural network-based mapping for excipient distribution to navigate high-dimensional design spaces.
• Integrating differentiable simulation with gradient-based optimization, the framework achieves low mean-squared error in meeting prescribed release kinetics.

Automated design frameworks for polypills employ computational methods to systematically co-optimize the geometry and spatial material distribution of single-dose oral formulations containing multiple active pharmaceutical ingredients and excipients. These frameworks transcend traditional ad hoc parameter-sweep approaches by enabling the systematic exploration and gradient-based optimization of high-dimensional design spaces defined by pill shape, excipient placement, and multi-phase drug release kinetics. Recent advances in additive manufacturing permit the physical realization of multi-material, patient-specific polypills, provided that their compositional and geometrical complexity can be robustly connected to desired therapeutic release profiles via physics-informed computational pipelines (Padhy et al., 9 Dec 2025).

1. Problem Motivation and Conventional Limitations

Modern polypills target coordinated, fixed-dose combination therapies and patient-customized pharmacokinetic profiles by integrating multiple APIs and excipients within a single oral dosage form. Additive manufacturing (AM) now allows for spatially heterogeneous, multi-material fabrication, greatly expanding the feasible design space for both geometry and composition. However, prevailing design practices rely on low-dimensional parameter sweeps, iteratively adjusting geometric or compositional variables to approximate target dissolution and release behaviors.

The primary limitations of these traditional approaches include:

• Scalability: Ad hoc parameter sweeps are not tractable in the intrinsically high-dimensional space of possible shapes and material layouts.
• Lack of formal guarantees: These methods provide no assurance that a prescribed, often tightly constrained release curve can be matched with available materials and shapes.
• Limited co-optimization: Current workflows rarely support the simultaneous optimization of geometry and material distribution at fine spatial resolution.

This context motivates the development of automated design frameworks, exemplified by PILLTOP, which couple physics-informed, differentiable simulation with topology optimization (TO) to robustly meet arbitrarily specified drug-release kinetics (Padhy et al., 9 Dec 2025).

2. Supershape Parametrization of Pill Geometry

The geometric freedom required for advanced polypill optimization is achieved with a generalized Gielis “supershape” function. The outer 2D topology is defined in local polar coordinates by the boundary radius

$$R(\varphi) = \left[ \left| \frac{1}{a} \cos\left(\frac{m\varphi}{4}\right) \right|^n + \left| \frac{1}{b} \sin\left(\frac{m\varphi}{4}\right) \right|^n \right]^{-1/n}$$

with the parameters $a, b > 0$ (axes scales), $m \in \mathbb{R}$ (rotational symmetry), and $n > 0$ (curvature exponent). The supershape is then affine-transformed (translation, rotation) with centroid $c=(c_x, c_y)$ and angle $\theta$.

A smooth phase-field “distance” function is constructed as

$$\Phi(x) = \|x_{\text{loc}}\| - R(\arg(x_{\text{loc}})), \quad x_{\text{loc}} = R(\theta) \cdot (x - c)$$

which is mapped to a sharp indicator via $\phi(x) = \frac{1}{2}[1 - \tanh(\Phi(x)/\mu)]$ ($\mu \ll 1$), resulting in $\phi \approx 1$ inside the pill and $\phi \approx 0$ outside.

The set of seven supershape design variables

$\zeta = [c_x, c_y, \theta, a, b, n, m]^T$

are bounded and optimized via a latent-space to sigmoid transformation. This parameterization provides strong geometric expressivity while remaining low-dimensional and differentiable.

3. Neural Network-Based Excipient Distribution

Material distribution within the supershape-defined pill domain is encoded using a coordinate-based multi-layer perceptron (MLP). For $S$ excipient types, the architecture takes $(x, y)$ coordinates (within the pill, $\phi \approx 1$), applies a Fourier feature mapping to the input, and employs two hidden ReLU layers with 40 neurons each, culminating in a Softmax output of dimension $S$. The output function

$$\gamma_s(x) \in [0, 1], \quad \sum_{s=1}^S \gamma_s(x) = 1$$

enforces non-negativity and the partition of unity, while the neural weights $w$ serve as design variables. Initialization is uniform: $\gamma_s \equiv 1/S$.

This continuous, mesh-independent representation supports analytic gradient computation $\partial \gamma / \partial w$, allowing end-to-end optimization. Discrete material assignment is encouraged using grayness constraints (see Section 5), ensuring near-discrete allocation of excipient types.

4. Coupled Dissolution–Diffusion Simulation

Release kinetics are governed by a coupled system of transient partial differential equations (PDEs) defined on the pill and solvent domain $\Omega$:

a) Phase-Field Dissolution (modified Allen–Cahn):

$$\frac{\partial \phi}{\partial t} + M_\phi [\psi'(\phi) - W\epsilon_t^2 \nabla^2 \phi] - f_{\text{diss}}(\phi, C) = 0$$

with $\psi(\phi) = W\phi^2(1 - \phi)^2$ and $$f_{\text{diss}} = -k(x)/[\rho_s(C_{\text{sat}} - C)|\nabla \phi|]$$. Here $k(x) = \sum_s \gamma_s(x) k^{(s)}$, $\rho_s$ is the solid density, and $C_{\text{sat}}$ saturation concentration. Initial $\phi(x,0)$ is the projected pill shape.

b) API Diffusion in Solvent (Fick's law):

$$\frac{\partial C}{\partial t} - \nabla \cdot [D(\phi)\nabla C] - S_{\text{source}}(\phi, C) = 0$$

with boundary condition $C=0$ (infinite sink) and initial zero concentration. The effective diffusivity is $$D(\phi) = D_{\text{solvent}} + (D_{\text{solid}} - D_{\text{solvent}}) h(\phi)$$, $h(\phi) = \phi^3 (10 - 15\phi + 6\phi^2)$. Source term is $$S_{\text{source}} = k(x)(C_{\text{sat}} - C)|\nabla \phi|$$.

Spatial discretization uses bilinear quadrilateral finite elements, temporal discretization is backward Euler. The resulting nonlinear system is solved by Newton-Raphson iterations at each time step (Padhy et al., 9 Dec 2025).

5. Optimization Problem Formulation

The design optimization aims to match a prescribed instantaneous mass-release rate $\dot{m}^*(t)$. The joint design variable set comprises supershape latents $\tilde{\zeta}$ and NN weights $w$. At discrete times $n$:

• Mass release rate:

$$\dot{m}_n = \frac{\rho_s}{\Delta t} \sum_e [\phi_n(x_e) - \phi_{n-1}(x_e)] v_e$$

• Objective: Minimize the mean-squared error

$$J = \frac{1}{N_t} \sum_{n=1}^{N_t} [\dot{m}_n - \dot{m}^*_n]^2$$

Subject to:

• Grayness suppression: encourages discrete $\gamma_s$ assignments,

$$g_r = \frac{1}{S N_e} \sum_{e,s} \gamma_s (1 - \gamma_s) - \xi \leq 0$$

• Minimum excipient volume fractions:

$$\lambda_s = \frac{\sum_e \gamma_s \phi v_e}{\lambda^*_s V_{\text{solid}}} - 1 \geq 0$$

constraints are aggregated via smooth-min $g_v \leq 0$.

These constraints are incorporated via a log-barrier:

$\mathcal{L} = J/J^0 + \psi(g_r) + \psi(g_v)$

with $\psi(g) = -\frac{1}{\tau}\log(-g)$ and continuation schemes for $\xi$, $\tau$.

Optimization is performed by unconstrained gradient descent (Adam), co-optimizing both the geometry and material variables.

6. Computational Implementation with Automatic Differentiation

The entire PILLTOP pipeline operates within JAX, supporting high-performance computation and end-to-end differentiability. Pipeline components include:

• Supershape projection ($\tilde{\zeta} \to \phi(x)$)
• Neural material mapping ($x \to \gamma(x)$)
• FE mesh mapping and assembly
• Time-stepped, Newton-solved PDEs
• Mass release and loss computation

Reverse-mode automatic differentiation (AD) is employed throughout, with the Implicit Function Theorem (IFT) used to differentiate through the final converged nonlinear system (not through all Newton iterations). Transient adjoint gradient computation exploits state checkpointing for efficient memory use. Hyperparameters typical of the framework include Adam with learning rate $\sim$8×10⁻³ and gradient clipping, with grayness slack and barrier weight scheduled over the iterations. Optimization of 100 steps typically completes in 30 minutes on a standard laptop, averaging about 10 Newton solves per time step over 100 time steps (Padhy et al., 9 Dec 2025).

7. Empirical Validation and Case Studies

The framework demonstrates its capabilities across a series of case studies:

• Single-excipient polypill: For monotonic release objectives, optimization of geometry alone (no NN) achieved $<1\%$ error to the target release, yielding a smooth, capsule-like supershape.
• Multi-material, non-monotonic release: Single-material TO failed on non-monotonic profiles (e.g., slow–spike–decay), whereas supershape + NN topology optimization (TO) achieved $<2\%$ mean-squared error (MSE) by spatially arranging fast and slow excipients.
• Sensitivity to initialization: Diverse initial supershape guesses ("spike," "circle," "sunflower") converged to distinct, locally optimal topologies but produced nearly identical release curves (MSE $\lesssim 10^{-3}$), allowing post hoc selection for secondary design criteria (e.g., comfort).
• Robustness to material degradation: With three grades of excipients (fresh, moderate, heavily aged, with decreasing $k$), and minimum-use constraints (e.g., $\geq 10\%$ aged), the optimizer adjusted allocations to approximately match the prescribed release, demonstrating resilience to on-demand manufacturing constraints.

In all scenarios, the framework consistently produced solutions aligning with the prescribed drug-release profiles within practical computational budgets (Padhy et al., 9 Dec 2025).