Machine learning of the anomalous diffusion of branched polymers in crosslinked networks

Abstract

Diffusion processes in complex environments, such as the extracellular matrix, are crucial in drug delivery. Common analytical theories developed for diffusion in such environments assume spherical, rigid particles. However, polymeric nanoparticles can often be aspherical and highly deformable, which introduces complexity beyond the spherical and rigid body assumptions. Moreover, it is challenging to measure classical diffusion coefficients under strong confinement or pronounced sub-diffusive conditions. We theoretically investigate the diffusion of branched polymers (bottlebrushes and stars) in polymeric mesh networks using coarse-grained molecular dynamics simulations. We introduce the Debye–Waller factor, a metric of confined mobility that we prove predicts long-time diffusion. We show that in relevant confinement regimes, elongated bottlebrushes have higher mobility than spherical stars. We can reliably predict the Debye–Waller factor from particle and network descriptors using Gaussian process regression. These results characterise the diffusion of arbitrary branched polymer nanoparticles and provide new, easily obtained metrics and protocols to design more efficient drug delivery carriers based on simple physical principles.

---

Introduction

Understanding the diffusion mechanism in confined environments is essential not only for fundamental transport theory^1–3 but also for designing effective drug delivery systems (DDS)^4–9 that operate within dense biological tissues. Nanomedicine offers powerful new ways to deliver therapeutic molecules, but its success ultimately depends on how well drug-loaded nanoparticles can navigate the body's crowded and complex interior.⁴ After entering the bloodstream, these tiny carriers must evade rapid clearance, reach the target tissue, and then thread their way through the narrow pores of the surrounding extracellular matrix (ECM).^10–12 The composition and pore size of ECM vary considerably among tissues, with approximately 85% of pores within the brain interstitial matrix being smaller than 125 nm.¹³ Therefore, how nanoparticles diffuse through ECM-like networks is a central question for modern drug-delivery research. Although the diffusion of rigid nanoparticles with well-defined shapes (spheres or rods) is now relatively well understood,^14–16 the transport mechanisms governing anisotropic and deformable polymeric nanoparticles remain largely unexplored. Understanding these mechanisms is especially important because particles used in practical DDS are rarely perfectly rigid and spherical.¹⁷ Beyond shape and deformability, the degree of branching directly tunes hydrodynamic size, segmental mobility, and steric interactions, making architecture a key design parameter for controlling nanoparticle diffusion. Most notably, experiments by Rabanel et al. showed that soft, anisotropic bottlebrush polymers cross the blood–brain barrier and penetrate the ECM-like gels far more efficiently than equal-size nanospheres, highlighting the practical importance of deformability and anisotropy in real DDS.¹⁸ Recent studies have also shown that soft nanoparticles can shrink in confined environments, boosting their diffusion compared to hard, non-deformable ones.¹⁹

The theoretical framework for understanding rigid spherical nanoparticle diffusion in polymer networks was substantially advanced by Cai et al., who developed a hopping diffusion model for large nanoparticles confined within network meshes.¹⁵ According to this model, particles larger than the mesh size of unentangled polymer networks become trapped in confinement cells.^15,20 At longer time scales, these particles can diffuse by overcoming the free energy barrier between neighbouring cells through a hopping process, so that the terminal diffusion coefficient depends on the particle-to-mesh size ratio. Subsequent simulation studies have refined and broadened this picture. Sorichetti et al. showed that spherical nanoparticles in polydisperse networks evolve from free diffusion to activated hopping and finally to trapping as confinement increases.²¹ Xu et al. introduced semiflexible network strands and found that matching the persistence length to the mesh size lowers the hopping barrier and intensifies subdiffusive behaviour.²² Recent work has further demonstrated that nanoparticle-network interactions can induce non-monotonic diffusivity through intermittent hopping,²³ and that increasing network stiffness raises the corresponding hopping barriers.²⁴ Shape effects add an additional layer of complexity: rigid rod-like systems exhibit direction-dependent diffusion scaling¹⁶ and can even outpace Stokes–Einstein predictions.²⁵ Lin et al. recently mapped the sphere-to-rod crossover, demonstrating that subtle changes in aspect ratio strongly tune confinement coupling,²⁶ thereby underscoring the profound role of morphology in confined transport. Yet, systematic work on non-spherical, deformable particles is still missing in this body of literature. Additionally, lower-resolution models of spherical or cylindrical nanoparticles lack information about the structure–property relation of polymeric micelles (such as the link between polymer/micelle architecture and diffusion) that is crucial to quantitatively inform experimental synthesis of new DDS.

In this work, we combine extensive coarse-grained molecular dynamics (CGMD) simulations with a Gaussian-process regression (GPR) framework to analyse the diffusion of deformable, anisotropic polymeric nanoparticles, bottlebrushes and stars, through flexible polymer networks of varying mesh size. We first introduce the Debye–Waller (DW) factor, which is extensively used in glass physics to quantify cage-scale vibrations^27–31 and is later shown to predict long-time diffusion³² in polymer melts. It is a practical predictor when direct diffusion coefficient measurements are prohibitively expensive, and we show here that it correlates well with the diffusion coefficients of our nanoparticles. Its value spans two orders of magnitude across our parameter space, showing the diverse dynamics within different systems. A GPR surrogate model, trained using the simulation input parameters, accurately predicts the DW factor with tight confidence intervals, transforming this observable into a powerful design metric. In the diffusive regime, bottlebrush carriers outperform their star-shaped counterparts with the same molecular weight across all mesh sizes. They remain diffusive even at the highest molecular weights examined. Stars, by contrast, become subdiffusive except under the weakest confinement and at the lowest molecular weight. This confirms the usefulness of anisotropic micelles as drug delivery carriers. In every system, the diffusion coefficient decreases monotonically with the confinement ratio, confirming geometric confinement as the primary control parameter,¹⁵ but with important secondary effects coming from polymer architecture. Collectively, our results go beyond existing hopping theories and incorporate fine-resolution molecular details into the diffusion of polymer nanoparticles, highlighting a new, simple metric of mobility of great use for the design of new drug nanocarriers and other transport applications.

Methods

Model

The translational diffusion of neutral polymeric nanoparticles within a neutral polymer network is investigated using CGMD simulations in implicit solvent. Both the diffusing nanoparticle and the network are modelled by the bead-spring model, each bead corresponding to a single Kuhn segment.³³ Non-bonded interactions between two beads are described by a truncated, purely repulsive Lennard-Jones (LJ) potential,where ε is the interaction strength, σ the monomer diameter, r the distance between two beads, and r_c = 2^1/6σ ≈ 1.12σ the interaction cutoff. All fundamental units, energy (ε), mass (m), and length (σ) are set to unity, leading to a unity simulation time (τ) defined as . For convenience, the Boltzmann constant is likewise set to unity (k_B = 1). Beads constituting the nanoparticle and the network are assigned identical diameter σ and mass m. Bonded interactions are modelled with a harmonic potential

U_H(r) = k(r − r₀)², (2)

with spring constant k = 1000ε/σ², and equilibrium length r₀ = 0.97σ for all bonds. The polymer network is constructed as a cubic lattice of beads. Each bond crosses the periodic boundaries, causing the lattice to repeat indefinitely and form an infinite three-dimensional mesh.

Simulation protocol

A nanoparticle and the polymer network are generated independently before coupling. The simulation cell coincides with the cubic mesh, and periodic boundary conditions are imposed along each Cartesian axis. The nanoparticle is inserted into a vacant lattice cell, after which stochastic velocities are assigned to all beads, sampled from a Maxwell–Boltzmann distribution at temperature T = 1.0. An energy minimisation stage precedes a soft relaxation executed in the canonical (NVT) ensemble for 5 × 10²τ, ensuring no bead overlaps. Subsequent equilibration employs a Langevin thermostat in the microcanonical (NVE) ensemble for an additional 5 × 10³τ. Finally, the data collection is done during a production run for 2.5 × 10⁶τ. Collective properties are measured at the last frame of this run. To elucidate the influence of nanoparticle architecture, we consider four different molecular weights, defined as the total number of beads in a nanoparticle. For each of them, we change the architecture, maintaining the molecular weight constant. Four different mesh sizes are considered to probe the confinement effects. This parametrisation yields 1024 state points. To improve the quality of dynamical properties, we run ten statistically independent replicas with different initial particle velocities for each state point, culminating in 10 240 simulations. The measured quantities are averaged over these runs, and their means are reported. Simulations are carried out with the LAMMPS (large-scale atomic/molecular massively parallel simulator) software with periodic boundary conditions in all three dimensions (https://www.lammps.org/).³⁴ Snapshots are rendered with OVITO.³⁵

Dynamical analyses

Mean squared displacement. The MSD³⁶ of the nanoparticle's centre of mass iswhere δt is the MD time step and N_t the number of stored frames. The average, therefore, exploits all overlapping time origins and is evaluated by averaging over multiple, overlapping time origins^37–39 which reduces statistical noise by roughly N_t/t compared with a single-origin estimate. We save trajectory frames in a hybrid linear-log scheme: dense, uniform sampling at short times and brief linear bursts around decade boundaries, followed by progressively coarser sampling within each decade. This yields many identical lag times, and we average MSD over overlapping time origins that are randomly sampled. For every state point, the resulting MSD trajectory is ensemble-averaged across ten statistically independent replicas that differ solely in their Maxwell–Boltzmann velocity realisations at T = 1.0. To classify dynamical regimes, we also evaluate the instantaneous logarithmic derivative (sometimes called the local exponent)⁴⁰ of the MSDwhose minima define the caging interval of the nanoparticle.

To determine the diffusive and subdiffusive regimes, we fit each MSD to the power law

〈r²(t)〉 = D_αt^α (5)

for t/τ > 10² where D_α is the generalised diffusion coefficient, and α is the diffusion exponent. Trajectories with α = 1 ± 0.05 are classified as diffusive, for which D_α becomes the standard diffusion coefficient, whereas those with α < 0.95 are classified as subdiffusive.

Debye–Waller factor. The instantaneous logarithmic derivative of the MSD, β(t), is the orthodox indicator of cage formation. Yet, it is notoriously susceptible to statistical volatility as it requires a numerical differentiation of the MSD. This operation magnifies high-frequency noise, thereby obscuring the true location of the global minimum. In several high confinement cases explored in our work, β(t) decays to its first minimum and fluctuates around that value for the remainder of the trajectory. Automated algorithms frequently misidentify a later, noise-corrupted lag time as the global minimum. To avoid this artefact, we locate the minimum by selecting a single global short-time caging window that we apply uniformly to all systems. In every data set, the minimum falls between lag times t/τ = 2 × 10² ∓ 10², so we pick the DW time (t_DW) within this range. We then calculate the corresponding Debye–Waller factor

〈u²〉 = 〈r²(t_DW)〉 (6)

fixing t_DW for all systems for consistency across the full data set. Five candidate windows are selected from the set

t_DW/τ ∈ {100, 150, 200, 250, 300}

each situated within the ballistic-to-caged (or diffusive) crossover. The DW factor reported for a given state point is the arithmetic mean of 〈u²〉 evaluated at these windows. Since 〈u²〉 is obtained directly from the MSD, it offers a more reliable measure of the nanoparticle's vibrational amplitude inside the cage. The DW factor has been previously established as an effective predictor of long-time cage relaxation dynamics in glass-forming systems and complex liquids.^41,42 Thus, it provides an unambiguous operational marker of caging dynamics even when the derivative criterion fails owing to persistent noise or extended plateaux in β(t). Averaging suppresses residual noise yet preserves the underlying plateau value. Since diffusive trajectories show an increasing MSD over the measured DW time window, their five samples spread further, giving a larger standard error of the mean (SEM). The SEM is stored for every state point and passed to the machine-learning model as an observation-noise estimate. This procedure replaces automated minima detection, which proved unreliable for trajectories that remain caged and led to noise-induced late time minima.

Predicting the dynamics by machine learning. GPR has emerged as a powerful tool for creating surrogate models of complex systems, with proven efficacy in MD applications.⁴³ Preliminary experiments showed that the raw DW factors span more than two orders of magnitude and exhibit heteroscedastic noise (the SEM grows as 〈u²〉 increases). Transforming the target to y = log〈u²〉 renders the noise nearly homoscedastic and linearises the dependence on the features, which improves both the fit quality and the calibration of predictive uncertainty. All statistics reported for the GPR (MSE, MAE, LML, and R²) therefore refer to log〈u²〉. It is then straightforward to do a back-transform of the quantities.
Gaussian-process model. We employ a zero-mean Gaussian process with an isotropic squared-exponential covariance functionwhere the length-scale hyper-parameters {[small script l]_i}, the signal variance σ_f² and the noise variance σ_n² are optimised by maximising the log-marginal likelihood (LML) over the training set. The entire GPR framework is implemented through the scikit-learn library.⁴⁴
Feature representation. The GPR model operates on a six-dimensional, normalised input vector,

x = (x₀, L, f, M, N, z)^T, (8)

whose components are defined as follows:

All predictors are linearly scaled to the closed interval [0, 1]. Parameters that are not physically meaningful for a given topology are set to the neutral midpoint 0.5 of the range [0, 1]. Choosing this central value minimises unintended trends and avoids boundary bias, so that the missing features do not create false patterns in the similarity measure. The model compares the architectures in equal terms, with predictions driven by the genuinely informative descriptors.

Training and validation. The dataset comprises the Debye–Waller amplitudes log〈u²〉 corresponding to the 1024 unique state points defined in the simulation settings. The model is trained on a randomly selected 80% subset of the data, which spans the full range of Debye–Waller factors, and evaluated on the remaining 20%. Predictive accuracy is quantified by the coefficient of determination R² and the mean squared error (MSE). Posterior standard deviations supplied by the GPR give an intrinsic estimate of predictive uncertainty for every state point. This protocol enables us to map the high-dimensional design space x onto a continuous estimate of the Debye–Waller factor, thereby providing rapid, uncertainty-aware predictions of nanoparticle caging dynamics without recourse to further molecular dynamics simulations.

Results and discussion

Studied system

We embed a neutral, cubic polymer network of mesh size L into the simulation box and introduce a single soft nanoparticle in one of the cells (Fig. 1a). Here, L denotes the distance between two crosslinking points in bonds and defines the mesh size, so that the number of beads between two crosslinks is L − 1. Six arms are connected to each crosslinking point. We explore four levels of confinement, L/σ = 7, 9, 11, 13, that span tight to loose meshes. We investigate two nanoparticle topologies (Fig. 1b and c):

Fig. 1 Schematic representation of the studied system. (a) A polymeric particle, either star-shaped or bottlebrush-like, diffusing within a polymeric network of mesh size L. We consider four different mesh sizes L/σ = {7, 9, 11, 13} (b) bottlebrush polymer characterised by backbone length N, side chain length M, and grafting interval m, i.e., the spacing between consecutive grafting points. The grafting density is defined as z = 1/m. We explore all topologies satisfying the constraints N > 3, M ≥ 3, and m ≤ 4, at fixed molecular weight M_w = N(1 + M/m). Four distinct molecular weights are examined M^BB_w = {120, 240, 360, 480} (c) star polymer with f arms of length M. The arm length is restricted by the condition 3 ≤ f ≤ 16, with fixed total molecular weight given by M_w = fM + 1, where the constant offset accounts for the central core bead. Four star molecular weights are studied: M^star_w = {121, 241, 361, 481}, this ensures approximate cross-compatibility between star and bottlebrush polymers at fixed size, with a minor offset reflecting their structural distinction.

• Bottlebrush polymers.

We model a linear backbone of N beads and graft flexible side chains of length M at every mth bead, so the grafting density reads z = 1/m. We impose N > 3, M ≥ 3, m ≤ 4, and keep the total molecular weight

at one of four values, M^BB_w = 120, 240, 360, 480.

• Star polymers.

We attach f arms of length M to a central core bead, with 3 ≤ f ≤ 16. The molecular weight becomes

M^star_w = fM + 1

which we choose as M^star_w = 121, 241, 361, 481 to align closely with the bottlebrush values; the extra ‘1’ accounts for the core bead. All beads, including the central core, have the same diameter.

Combining the four mesh sizes with all admissible nanoparticle architectures yields 1024 distinct state points. Since the nanoparticle beads share the same diameter σ and the interaction parameters as the network, any change in the dynamics originates solely from topology, confinement, or their interplay, not from chemical differences.

Architecture-dependent anomalous diffusion

Fig. 2 shows the MSD for all 1024 state-points, each coloured by its DW factor. At very short times, every nanoparticle follows the ballistic regime (〈r²〉 ∝ t²), indicated by the black t² guide-line. Beyond t/τ ≈ 10², the trajectories depart from ballistic motion; depending on the specific architecture and level of confinement, they either pause in a subdiffusive regime (caged), where the MSD forms a plateau even at longer times, or progress towards a diffusive regime (black t guide-line). There are 180 diffusive and 844 subdiffusive state points in our dataset, meaning that 180 systems reach the diffusive regime within the simulated time window. The plateau height is the value of the DW factor, so the colour scale by DW factor value reflects how strongly each nanoparticle is confined by the surrounding mesh. The inset illustrates the determination of the DW factor for an example system: bottlebrush-like polymer with backbone length N = 120, side chain length M = 6 and grafting density z = 0.33 in a network of mesh size L = 13σ. Black symbols show the MSD, and red symbols its logarithmic derivative. Since measurement noise makes it difficult to locate the true minimum of the numerical derivative automatically (see Methods for details), we fix a set of short-time windows t_DW/τ = 100, 150, 200, 250, 300. We compute the DW factor at each of these windows and report their average, propagating the associated standard error so that the upcoming Gaussian-process regression can make use of this uncertainty.

Fig. 2 Ensemble-averaged mean-squared displacement (MSD) across all systems. Each curve is colour-coded by the corresponding carrier's Debye–Waller (DW) factor, (〈u²〉). The black lines denote the ballistic (∝t²) and symptotic diffusive regime (∝t). The red-shaded band (100 ≤ t/τ ≤ 300) highlights the beginning of the intermediate “caging” point, wherein carriers experience transient confinement imposed by the surrounding network. A diverse spectrum of dynamic behaviour is observed, ranging from normal diffusion to pronounced subdiffusion, with diffusion exponents spanning a broad interval (0.14 ≤ α ≤ 1.04). The inset illustrates the MSD (top) and its logarithmic derivative (bottom) for a representative system. The selection of the caging interval is informed by visual inspection across all datasets and fixed for all systems with appropriate error bars, see the Methods section for details.

The whole ensemble displays a broad spectrum of dynamical responses, with some nanoparticles diffusing rapidly while others remain strongly caged. Most trajectories never attain a linear slope within the accessible simulation time window, confirming sustained subdiffusive behaviour. Curves with small plateau values (shaded towards orange) linger longer in the caged region before any power-law growth resumes; those with larger plateaux (shaded towards purple) never enter a caging regime (or escape it earlier) and approach a diffusive (or nearly diffusive) regime.

Fig. 3 reports all DW factor amplitudes as a function of the level of confinement of the polymers, as well as the correlation between DW factor and diffusion coefficient for diffusive systems. We define the confinement ratio where R_h is the hydrodynamic radius given by

where r⃑_i and r⃑_j are all possible positions of chain monomers with i ≠ j.⁴⁵ In hopping theory, is the ratio of the diffusing particle's diameter to the network mesh size.¹⁵ Here, we adopt a similar approach by choosing 2R_h as the nanoparticle's diameter. The R_H range is split into thirty equal-width bins for every state point, and the DW factors are averaged inside each bin. Black symbols are these bin-wise means. The pink envelope shows the standard error of the mean for each bin.

Fig. 3 Dependence of the Debye–Waller factor on the confinement ratio (all systems) and on the diffusion coefficient (diffusive systems only). (a) Binned Debye–Waller factor as a function of confinement ratio. The Debye–Waller (DW) factor is computed for all systems and represented as a function of the confinement ratio 2R_h/L. Each data point corresponds to the mean DW factor within one of 30 equally spaced bins along the horizontal axis. The shaded red envelope denotes the propagated error associated with each measurement. At low confinement ratios, particles exhibit higher mobility. (b) Scatter plot of DW factor versus the diffusion coefficient D_α for systems classified as diffusive (α ≈ 1). Each point denotes a state point; the strong positive correlation (Spearman r = 0.98) indicates that larger 〈u²〉 coincides with larger D_α, substantiating the DW factor as a robust early-time predictor of diffusion.

The uncertainty band in Fig. 3a widens markedly at small confinement ratios since, in that regime, nanoparticles leave the ballistic regime and enter a diffusive one without becoming caged. When we evaluate the DW factor at five, short, fixed time windows, the particle has already travelled a noticeable distance between successive windows; the resulting spread increases the standard error on the mean. As the confinement ratio increases, the network increasingly restrains the particle, the MSD plateau stabilises across all five windows, and the standard error decreases. Hence, the measurement grows more precise for strongly confined systems and less precise for weakly confined, diffusive ones. Using the DW factor in the weak confinement regime is crucial: it supplies a single, universally defined dynamical descriptor that we can evaluate for every trajectory, even those that never reach the long-time diffusive limit. A conventional diffusion coefficient would be undefined in these subdiffusive cases, leaving gaps in the data or requiring much longer simulations. The DW factor gives us a continuous metric across the entire spectrum from freely diffusing to strongly caged nanoparticles, enabling coherent trend analysis and machine-learning prediction for the whole dataset.

To fully isolate the effect of architecture at fixed molecular weight, we show DW factor as a function of a different confinement metric, which we define as the volumetric confinement ratio

where the numerator defines the excluded volume of the nanoparticle. Note that M_w corresponds to the number of beads, and that CR_vol is uniquely defined for each different [M_w,L] set of parameters studied. The denominator is the volume of one network cell. With this definition, we show the distribution of 〈u²〉 as a function of confinement and highlight the variations between different architectures in Fig. 4.

Fig. 4 Statistical distribution of the Debye–Waller factor as a function of volumetric confinement ratio. A violin plot depicts the full distribution of the DW factor 〈u²〉 obtained from all simulated systems at fixed molecular weight and mesh size. The width of the violin plot gives a kernel-density estimation of the frequency of 〈u²〉. The black rectangle corresponds to the inter-quantile range (25th and 75th percentiles). The white horizontal line is the sample median. The black line around the coloured violin is the 95% empirical range, truncated at the extreme observed values. The colour map (legend at right) encodes the logarithmic scale value of 〈u²〉. The snapshots correspond to architectures at the lowest confinement ratio. The top snapshot is a bottlebrush with N = 20, M = 5, m = 1 (the fastest within its set), whereas the bottom snapshot is a star with N = 4, M = 30 (the slowest within its set). In both cases, the molecular weight and the mesh size are fixed.

The DW factor decays steeply with increasing confinement. At weak confinement (leftmost violins), the distributions are broad, suggesting the possibility of achieving faster carriers by changing the architecture at fixed M_w. In the strongest confinement regime (rightmost violins), the distributions collapse to nearly point-like shapes, indicating that particle motion is highly restricted and insensitive to further architectural changes. Therefore, the progressive narrowing of the violins highlights a confinement-dominated transition from heterogeneous mobile systems to uniformly caged ones. For weak confinement, we generally observe that elongated bottlebrushes have higher DW factors than stars with the same M_w. A representative example is shown with the two snapshots in Fig. 4. The same conclusion can be drawn from the architecture-resolved D_α.

Fig. 5 shows the effective diffusion coefficient D_α for every diffusive system as a function of confinement ratio. Symbol colours encode the nanoparticle molecular weight M_w (runs from pale to dark progressively), while the shape distinguishes its topology (circles: bottlebrush, stars: star). To compare the star and bottlebrush topologies, we exclude the bottlebrush-like systems with grafting density z < 1 and only consider those with z = 1 here.

Fig. 5 The effect of architecture on the diffusivity in the diffusive regime. The diffusion coefficient D_α obtained from the power-law fits of the MSD is plotted against the confinement ratio. Each point represents a bottlebrush (circles) or a star (stars) polymer, colour-coded by molecular weight M_w. A clear inverse relationship is observed: the effective diffusivity decreases systematically as the confinement ratio increases. At low confinement, systems with lower molecular weights exhibit higher diffusivities, reflecting reduced size-induced constraints. Notably, a broader spread in diffusivity is observed at low M_w where architectural features more strongly influence the dynamics, beyond the simple predictions of the hopping model.

The primary control is the confinement effect. Regardless of topology or molecular weight, D_α rapidly decays as the confinement ratio increases. Confinement remains the dominant control parameter even when particles differ by shape or mass. At a given confinement ratio, the data spread by up to half an order of magnitude, showing the secondary effect of architecture. Most of this spread appears at the weak confinement range where .⁴⁶ This is particularly significant because in molecular transport, the confinement typically does not exceed unity.⁴⁶ Our results, therefore, highlight the weak confinement regime as the relevant window where architecture-dependent differences emerge. We note the absence of star topology at high confinement. Since star polymers radiate from a central core, their size grows rapidly with molecular weight, causing them to become caged at lower confinement ratios than bottlebrushes. In our dataset, all star topologies with M_w > 120 fall into the subdiffusive class, whereas we can identify at least one bottlebrush architecture that remains in the diffusive regime across the full M_w range examined. This contrast suggests that, once other factors are equalised, bottlebrush topology favours faster long-time transport than its star-shaped counterpart. This trend aligns with a recent experiment showing bottlebrushes as faster diffusing drug carriers than stars.¹⁸

To pinpoint the contribution of nanoparticle architecture, we compute the Spearman correlation coefficient over all simulated data points between 〈u²〉 and some of the nanoparticle structural descriptors in Fig. 6.

Fig. 6 Spearman rank correlations between the Debye–Waller factor and structural descriptors. Bars give the Spearman correlation coefficient between the target observable and each structural variable: confinement ratio , mesh size L, molecular weight (M_w), hydrodynamic radius (R_h), and the component of the gyration tensor (R_g,i). The strongest anti-correlation is observed for , confirming it as the primary mobility limiter. Overall, the ranking highlights confinement geometry as dominant, with nanoparticle size/shape parameters having secondary yet inhibitory effects on diffusion. Correlations are calculated over the full dataset.

Here, and where λ_i² is the gyration tensor eigenvalues of the nanoparticle. The components are sorted in value R_g,3 ≥ R_g,2 ≥ R_g,1. The trace of the gyration tensor gives the radius of gyration .

Particle-level confinement emerges as the dominant determinant of anomalous diffusion because it couples directly to the particle's hydrodynamic size and thus to the free volume available for motion. Nevertheless, architecture-induced effects are evident in the full distribution of 〈u²〉 (Fig. 4) where systems of identical confinement ratio but different topologies display markedly different mobilities. Note also that the correlation between the Debye–Waller factor and most structural features is never very high (lower than 0.5 in absolute value for all molecular descriptors). This disparity signals that the correlations are too coarse to capture the complex interplay between size and shape, motivating the adoption of a more predictive and general framework.

Statistical learning of nanoparticle mobility

GPR is particularly suitable in our case for multiple reasons. First, it gives state-of-the-art predictive accuracy by learning arbitrary nonlinear relations.^47,48 Second, its kernel hyperparameters and automatic-relevance-determination length scales provide quantitative, physically meaningful measures of each feature's influence.⁴⁹ Finally, the posterior variance supplies uncertainty estimates.⁵⁰ To evaluate the performance of our feature selection in GPR, we apply the leave-one-feature-out (LOFO) protocol and plot the inverse of the fitted kernel length scale in Fig. 7.

Fig. 7 Relative importance of input features in the Gaussian-process model for the Debye–Waller factor. (a) Inverse length-scale 1/[small script l] extracted from the fitted kernel; larger values denote stronger local sensitivity of prediction to that feature. The values are normalised within the feature vector so that the sum of all bars adds up to 1. (b) Decrease in the coefficient of determination, ΔR², obtained by the leave-one-feature-out protocol. Each bar shows the loss in predictive power when the indicated variable is omitted from model training.

The LOFO implies how much a given feature improves the model accuracy globally. It shows whether the model can still predict well without the omitted feature. Thanks to the feature encoding we employ in Table 1, we can omit the categorical variable x₀. However, it still increases the log marginal likelihood, so we keep it in our feature vector.

Table 1 Normalised input features for the GPR model

Symbol	Description	Encoding for absent parameter
x 0	Topology identifier (star = 1; bottlebrush = 0)	—
L	Network mesh size	—
f	Number of arms (stars)	f = 0.5 for bottlebrushes
M	Arm length (stars) or side-chain length (bottlebrushes)	—
N	Bottlebrush backbone length	N = 0.5 for stars
z	Bottlebrush grafting density	z = 0.5 for stars

---

Having established the relative importance of the input features, we next examine the overall predictive power of the model. Fig. 8 shows how well GPR performs on our simulation data. We trained the model on 80% of the whole data set, and present the predictions on the remaining test set here. Each circle represents one of those unseen test data points. The horizontal axis is the actual measured DW factor with horizontal error bars showing the measurement error. The vertical axis shows the GPR predictions with the model's 95% confidence interval for each prediction. The red dashed line is the identity line y = x.

Fig. 8 Comparison between predicted and actual values of the Debye–Waller (DW) factor for the test dataset. This figure demonstrates the predictive capability of Gaussian process regression (GPR) in estimating the DW factor based solely on the input parameters of the simulation. The horizontal axis denotes the ground-truth DW factor as computed from the simulation trajectories, while the vertical axis corresponds to the GPR-predicted values. All values have been normalised to the interval [0, 1] for the purposes of model training and evaluation. The red dashed line represents the identity relation y = x, corresponding to perfect predictive performance. Each black data point corresponds to one test instance, with vertical and horizontal grey error bars indicating the 95% confidence interval (CI) of the GPR prediction and the propagated uncertainty of the actual DW factor, respectively. The close alignment of the data with the identity line and the relatively narrow confidence intervals indicate that the GPR model can accurately infer the DW factor from simulation inputs without requiring access to the particle trajectories.

We emphasise the overall agreement between the measured values and the GPR predictions. The test points cluster tightly around the identity line, indicating that the GPR recovers DW factor of previously unseen state points without noticeable bias. Optimisation yields the kernel

with length scales [small script l] = (4.32, 1.72, 2.44, 0.251, 0.232, 3.60) in the order (x₀, L, f, M, N, z). The corresponding signal variance is σ_f² = 8.98² and the white-noise variance is σ_n² = 0.0064. On the 20% test set the model attains a mean-squared error MSE = 5.5 × 10⁻⁴, mean absolute error MAE = 1.72 × 10⁻², coefficient of determination R² = 0.989, and a log-marginal likelihood . These scores refer to the target y = log〈u²〉. These metrics confirm that six normalised input features in eqn (8) encode the essential variables governing the cage amplitude. The 95% confidence intervals widen at small 〈u²〉 (weak confinement, larger measurement noise) and tighten strongly for caged systems with more precise measurements. 99.5% of the test points fall inside their predicted 95% intervals, demonstrating that the model's uncertainty estimates are well calibrated. This accuracy and reliable error estimation allow GPR to replace time-consuming MD simulation. It delivers predictions and confidence bounds in milliseconds, making rapid screening of new particle designs possible, whereas a single MD run can require several days to weeks of computation. We note that, as with all Gaussian process models, predictive accuracy is limited to interpolation within the training range; extrapolation beyond it reverts to the prior with broad confidence intervals.⁴⁹

Limitations of the hopping model

Fig. 9 compares the diffusion coefficient predicted by the hopping model D_hop with the effective coefficient D_α obtained from our MSD analysis for the subset of systems that are diffusive (α ≈ 1). We compute the hopping coefficient according towhere b is the Kuhn length of a polymer network strand, τ_x the Rouse relaxation time of a network strand, and the confinement ratio.¹⁵ Based on the phantom network model,⁵¹b is defined as
where R is the distance between two cross-linking points in the network, and L is the network strand length. Based on the same model, the Rouse relaxation τ_x is proportional to the number of monomers on the strand τ_x ∼ N_x². Since we have a regularly spaced monodisperse network R ≃ L and N_x ∼ L, hence, we approximate D_hop for our systems by

Fig. 9 Comparison between the hopping diffusion coefficient and the effective diffusion coefficient extracted from the power-law fit r² = D_αt^α and reported for the cases α ≈ 1. (a) Correlation between the diffusion coefficient estimated from hopping model (D_hop) and the effective long-time diffusion coefficient (D_α), extracted from the MSD. Each point represents a distinct system, and the red dashed line corresponds to the identity D_hop = D_α, indicating perfect agreement. While D_hop correlates with D_α with a Spearman coefficient 0.91, we observe some systematic non-linearities for systems with higher mobility. A slightly larger Spearman coefficient (0.98, see Fig. 3) between D_α and 〈u²〉 supports that the DW factor is indeed a better metric in predicting the long time diffusion. (b) The residual D_hop − D_α as a function of confinement ratio. The difference converges to zero at high confinement, suggesting that discrete hopping events predominantly govern the dynamics under strong spatial restriction. In contrast, the increased mobility and potential for continuous motion at low confinement lead to larger discrepancies between the two diffusion estimates.

Each point in Fig. 9a shows one state point in the diffusive set. The x-axis gives MSD based D_α and the y-axis shows the hopping estimate D_hop. The data points cluster close to the red identity line, but sit slightly above it, indicating that the hopping model overestimates the diffusion coefficient. To make this clearer, we plot the difference D_hop − D_α as a function of confinement ratio 2R_h/L in Fig. 9b. The residuals converge to zero as CR increases. In the weak confinement regime , the residual becomes nonzero and scattered, showing that the simple hopping dynamics cannot capture the effect. Deviations from the hopping model emerge as the diffusing particle is not trapped within the network cells. By contrast, in the strongly confined regime , the transport mechanism is well captured by the hopping model, since the nanoparticle relies on pores opening with the thermal fluctuations. In this case, most of our diffusion coefficients overlap with the predictions of D_hop, consistent with the model's validity under strong confinement. The hopping model, therefore, provides an accurate description only in the strongly confined regime where the network forces the particle to move by well-defined jumps. A key premise of the hopping model is that the diffusing entity is a rigid, spherical particle whose geometry remains fixed between cage escapes. Our nanoparticles violate this premise: they are anisotropic, possess internal flexibility, and can deform in response to the local network environment. Under strong confinement, the mesh constrains both translation and shape, so motion still resembles discrete hops, and the model performs well. In looser meshes, however, the particle has room to elongate, contract, or reorient while drifting between network strands. These shape fluctuations add a continuous component to the overall displacement that the rigid-sphere hopping picture cannot capture.

Conclusions

Our combined MD and GPR study clarifies how deformability, anisotropy, and branching hierarchy determine polymer nanoparticle motion through polymer-like ECM. We first validate the Debye–Waller factor, a measure of mobility during caging, as an efficient and powerful predictor of long-time diffusion. Its wide variation across state points captures the full spectrum from trapped to freely diffusing nanoparticles, and shows large variations in mobility for systems at fixed size and confinement degree based on their shape. A Gaussian process surrogate then predicts DW factor directly from molecular parameters, thereby predicting both caging strength and long-time diffusion. This turns a noisy observable into a reliable design metric. Most architectures studied here exhibit anomalous diffusion. When normal diffusion emerges, the measured diffusion coefficients follow the hopping model at high confinement ratios, but diverge at weak confinement, highlighting the model's limitations for soft, anisotropic particles. Bottlebrush carriers remain diffusive at higher molecular weight and tighter meshes than stars, confirming that branching can mitigate confinement penalties. Based on similar previous work²² where crosslink connectivity is systematically varied, we do not expect qualitative changes to our results if the network topology is changed by varying the number of connecting strands to a crosslink point. More interesting would be to introduce pore size dispersity, a realistic and important feature of biological networks for which we cannot easily predict the results based on our current data. Overall, these results bridge theoretical models with the practical design of branched nanocarriers with tailored transport properties. Our findings provide practical metrics for engineering polymeric nanoparticles with enhanced transport capabilities in confined environments such as tumour ECM, mucus membranes, and other complex biological gels by elucidating the relationship between molecular architecture and diffusion behaviour.

Conflicts of interest

There are no conflicts to declare.

Data availability

In an effort to promote Open Science practices, the LAMMPS code used to perform our simulations is available on the GitHub page of our group https://github.com/giuntoli-group/confined-polymer-dynamics.

Acknowledgements

We thank the Center for Information Technology of the University of Groningen for their support and for providing access to the Hábrók high performance computing cluster. We thank Anton van Beek for the useful discussions and Alessia Vanni for helping with the design of our TOC graphics.

References

1. P. S. Burada, P. Hänggi, F. Marchesoni, G. Schmid and P. Talkner, Diffusion in confined geometries, ChemPhysChem, 2009, 10, 45–54.
2. M. Bruna and S. J. Chapman, Diffusion of finite-size particles in confined geometries, Bull. Math. Biol., 2014, 76, 947–982.
3. Y. Li, R. Mei, Y. Xu, J. Kurths, J. Duan and R. Metzler, Particle dynamics and transport enhancement in a confined channel with position-dependent diffusivity, New J. Phys., 2020, 22, 053016.
4. E. Blanco, H. Shen and M. Ferrari, Principles of nanoparticle design for overcoming biological barriers to drug delivery, Nat. Biotechnol., 2015, 33, 941–951.
5. S. Mitragotri, P. A. Burke and R. Langer, Overcoming the challenges in administering biopharmaceuticals: formulation and delivery strategies, Nat. Rev. Drug Discovery, 2014, 13, 655–672.
6. M. W. Tibbitt, J. E. Dahlman and R. Langer, Emerging frontiers in drug delivery, J. Am. Chem. Soc., 2016, 138, 704–717.
7. S. Adepu and S. Ramakrishna, Controlled drug delivery systems: current status and future directions, Molecules, 2021, 26, 5905.
8. T. C. Ezike, U. S. Okpala, U. L. Onoja, C. P. Nwike, E. C. Ezeako, O. J. Okpara, C. C. Okoroafor, S. C. Eze, O. L. Kalu and E. C. Odoh, et al., Advances in drug delivery systems, challenges and future directions, Heliyon, 2023, 9, e17488.
9. E. Dong, Q. Huo, J. Zhang, H. Han, T. Cai and D. Liu, Advancements in nanoscale delivery systems: optimizing intermolecular interactions for superior drug encapsulation and precision release, Drug Delivery Transl. Res., 2025, 15, 7–25.
10. F. Alexis, E. Pridgen, L. K. Molnar and O. C. Farokhzad, Factors affecting the clearance and biodistribution of polymeric nanoparticles, Mol. Pharmaceutics, 2008, 5, 505–515.
11. S. Barua and S. Mitragotri, Challenges associated with penetration of nanoparticles across cell and tissue barriers: a review of current status and future prospects, Nano Today, 2014, 9, 223–243.
12. N. Hoshyar, S. Gray, H. Han and G. Bao, The effect of nanoparticle size on in vivo pharmacokinetics and cellular interaction, Nanomedicine, 2016, 11, 673–692.
13. E. A. Nance, G. F. Woodworth, K. A. Sailor, T.-Y. Shih, Q. Xu, G. Swaminathan, D. Xiang, C. Eberhart and J. Hanes, A dense poly (ethylene glycol) coating improves penetration of large polymeric nanoparticles within brain tissue, Sci. Transl. Med., 2012, 4, 149ra119.
14. U. Yamamoto and K. S. Schweizer, Theory of nanoparticle diffusion in unentangled and entangled polymer melts, J. Chem. Phys., 2011, 135, 224902.
15. L.-H. Cai, S. Panyukov and M. Rubinstein, Hopping diffusion of nanoparticles in polymer matrices, Macromolecules, 2015, 48, 847–862.
16. Y. Chen, Z. Xiang, H. Ren, F. Guo, V. Ganesan and J. Liu, Anisotropic and Non-Gaussian Diffusion of Thin Nanorods in Polymer Networks, Macromolecules, 2024, 57, 5105–5118.
17. Q. Yu, M. G. Roberts, L. Houdaihed, Y. Liu, K. Ho, G. Walker, C. Allen, R. M. Reilly, I. Manners and M. A. Winnik, Investigating the influence of block copolymer micelle length on cellular uptake and penetration in a multicellular tumor spheroid model, Nanoscale, 2021, 13, 280–291.
18. J. M. Rabanel, et al., Deep Tissue Penetration of Bottle-Brush Polymers via Cell Capture Evasion and Fast Diffusion, ACS Nano, 2022, 16, 21583–21599,  DOI:10.1021/acsnano.2c10554.
19. P.-L. Latreille, V. Adibnia, A. Nour, J.-M. Rabanel, A. Lalloz, J. Arlt, W. C. Poon, P. Hildgen, V. A. Martinez and X. Banquy, Spontaneous shrinking of soft nanoparticles boosts their diffusion in confined media, Nat. Commun., 2019, 10, 4294.
20. Z. E. Dell and K. S. Schweizer, Theory of localization and activated hopping of nanoparticles in cross-linked networks and entangled polymer melts, Macromolecules, 2014, 47, 405–414.
21. V. Sorichetti, V. Hugouvieux and W. Kob, Dynamics of nanoparticles in polydisperse polymer networks: From free diffusion to hopping, Macromolecules, 2021, 54, 8575–8589.
22. Z. Xu, X. Dai, X. Bu, Y. Yang, X. Zhang, X. Man, X. Zhang, M. Doi and L.-T. Yan, Enhanced heterogeneous diffusion of nanoparticles in semiflexible networks, ACS Nano, 2021, 15, 4608–4616.
23. B.-R. Zhao and B. Li, Molecular simulation of hopping mechanisms of nanoparticles in regular cross-linked polymer networks, J. Chem. Phys., 2022, 157, 104901.
24. Y. Lu, X.-Y. Liu and G.-H. Hu, Double-spring model for nanoparticle diffusion in a polymer network, Macromolecules, 2022, 55, 4548–4556.
25. M. E. Sudakova, V. M. Cherepanov, T. N. Ponamareva, A. G. Savchenko, M. A. Abakumov and A. A. Nikitin, Unconventional Fast Nanorod Diffusion in Entangled Solutions of PEO Revealed by Moössbauer Spectroscopy, Macromolecules, 2024, 57, 11429–11437.
26. T.-W. Lin and C. E. Sing, Effect of penetrant–polymer interactions and shape on the motion of molecular penetrants in dense polymer networks, J. Chem. Phys., 2024, 160, 114905.
27. F. W. Starr, S. Sastry, J. F. Douglas and S. C. Glotzer, What do we learn from the local geometry of glass-forming liquids?, Phys. Rev. Lett., 2002, 89, 125501.
28. D. S. Simmons, M. T. Cicerone, Q. Zhong, M. Tyagi and J. F. Douglas, Generalized localization model of relaxation in glass-forming liquids, Soft Matter, 2012, 8, 11455–11461.
29. B. A. Pazmiño Betancourt, P. Z. Hanakata, F. W. Starr and J. F. Douglas, Quantitative relations between cooperative motion, emergent elasticity, and free volume in model glass-forming polymer materials, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 2966–2971.
30. J. F. Douglas, B. A. P. Betancourt, X. Tong and H. Zhang, Localization model description of diffusion and structural relaxation in glass-forming Cu–Zr alloys, J. Stat. Mech.:Theory Exp., 2016, 2016, 054048.
31. M. Becchi, A. Giuntoli and D. Leporini, Molecular layers in thin supported films exhibit the same scaling as the bulk between slow relaxation and vibrational dynamics, Soft Matter, 2018, 14, 8814–8820.
32. W. Xia, N. K. Hansoge, W.-S. Xu, F. R. Phelan Jr, S. Keten and J. F. Douglas, Energy renormalization for coarse-graining polymers having different segmental structures, Sci. Adv., 2019, 5, eaav4683.
33. K. Kremer and G. S. Grest, Dynamics of entangled linear polymer melts: A moleculardynamics simulation, J. Chem. Phys., 1990, 92, 5057–5086.
34. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al., LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
35. A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool, Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012.
36. D. Frenkel and B. SmitUnderstanding molecular simulation: from algorithms to applications, Elsevier, 2023.
37. H. Qian, M. P. Sheetz and E. L. Elson, Single particle tracking. Analysis of diffusion and flow in two-dimensional systems, Biophys. J., 1991, 60, 910–921.
38. Y. He, S. Burov, R. Metzler and E. Barkai, Random time-scale invariant diffusion and transport coefficients, Phys. Rev. Lett., 2008, 101, 058101.
39. X. Michalet, Mean square displacement analysis of single-particle trajectories with localization error: Brownian motion in an isotropic medium, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2010, 82, 041914.
40. A. Nandi, D. Heinrich and B. Lindner, Distributions of diffusion measures from a local mean-square displacement analysis, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2012, 86, 021926.
41. A. Widmer-Cooper and P. Harrowell, Predicting the Long-Time Dynamic Heterogeneity in a Supercooled Liquid¡? format?¿ on the Basis of Short-Time Heterogeneities, Phys. Rev. Lett., 2006, 96, 185701.
42. L. Larini, A. Ottochian, C. De Michele and D. Leporini, Universal scaling between structural relaxation and vibrational dynamics in glass-forming liquids and polymers, Nat. Phys., 2008, 4, 42–45.
43. F. Noé, A. Tkatchenko, K.-R. Müller and C. Clementi, Machine learning for molecular simulation, Annu. Rev. Phys. Chem., 2020, 71, 361–390.
44. F. Pedregosa, et al., Scikit-learn: Machine Learning in Python, J. Mach. Learn. Res., 2011, 12, 2825–2830.
45. I. Teraoka Polymer Solutions, John Wiley & Sons, Ltd, 2002, ch. 3, pp. 167–275.
46. B. Mei, T.-W. Lin, G. S. Sheridan, C. M. Evans, C. E. Sing and K. S. Schweizer, How segmental dynamics and mesh confinement determine the selective diffusivity of molecules in cross-linked dense polymer networks, ACS Cent. Sci., 2023, 9, 508–518.
47. D. Ye and M. Guo, Gaussian process learning of nonlinear dynamics, Commun. Nonlinear Sci. Numer. Simul., 2024, 138, 108184.
48. J. Wang, An intuitive tutorial to Gaussian process regression, Comput. Sci. Eng., 2023, 25, 4–11.
49. C. K. Williams and C. E. RasmussenGaussian processes for machine learning, MIT Press, Cambridge, MA, 2006, vol. 2.
50. J. Wenger, G. Pleiss, M. Pförtner, P. Hennig and J. P. Cunningham, Posterior and computational uncertainty in Gaussian processes, Adv. Neural Inf. Process. Syst., 2022, 35, 10876–10890.
51. M. Rubinstein, R. H. Colby, et al., Polymer physics, Oxford University Press, New York, 2003, vol. 23.

---