A small network inside a mechanistic solver
The paper does not replace the transport model with a black box. It keeps the population balance equation (PBE), a partial differential equation that tracks how many crystals of each size and shape exist over time, and keeps the solute mass balance that describes the liquid phase. Only the convection term, the rate at which crystals grow or dissolve, is swapped between a published empirical law and a neural-network representation. Because the network sits inside the governing equations, the conservation structure of the model is preserved.
The solver is a two-dimensional finite-volume discretization in the size coordinates with forward Euler time stepping, built in JAX and run on GPUs. Its key property is end-to-end differentiability: automatic differentiation propagates gradients from a trajectory mean-squared-error loss back through the solver to whichever inputs are declared free. The free variables change by task. For model building they are the rate-law parameters plus the initial seed distribution and initial concentration; for process optimization they are the temperature profile, optionally together with the initial conditions.
This matters because the cost of a gradient becomes effectively independent of the number of parameters. The authors report that 500 optimization epochs over roughly 4.7×10^5 parameters, dominated by 11 particle size and shape distributions of 360×120 each, take 5 to 6 minutes. A benchmark MATLAB model using numerical gradients needed upwards of 15 minutes to fit a two-parameter model, more than twice as long while optimizing five orders of magnitude fewer parameters.
Exact initial conditions buy accuracy but break the law
The experiments are batch temperature-driven dissolution of L-glutamic acid in water, monitored online by an instrument that characterizes crystals in two characteristic lengths, length L1 and width L2. The dataset contains 13 real trajectories with varying temperature profiles and initial conditions. In a dynamical system of this kind, trajectories are sensitive to where they start, and the measured seed distribution and concentration carry instrument noise and device limitations.
When the authors treated those measured initial conditions as ground truth and let the network fit only the rate law, predictive accuracy improved but the learned kinetics were physically inconsistent. The learned supersaturation dependence failed to vanish at equilibrium, where supersaturation S equals 1, and an artificial temperature dependence appeared. The paper's explanation is that a flexible network coupled to a solver sensitive to its initial state cannot distinguish error in the data from the underlying dynamics, so it fits both.
That failure is only visible because the learned component is a small, inspectable function embedded in a mechanistic model. A black-box model of the whole system would have absorbed the same inconsistency, and an evaluation based only on predictive performance would not have flagged it.
Assimilating initial conditions recovers a consistent law
The fix was to make the initial conditions free variables and optimize them jointly with the rate law. The recovered kinetics then behaved as the physics requires: the S dependence approaches zero at S = 1, and temperature dependence becomes negligible, consistent with the published empirical baseline. The learned function is of the same order of magnitude as that baseline while showing a non-linear response in S that the linear empirical model cannot describe.
The calibrated initial conditions moved only slightly. In E7, the validation experiment with the largest shift, mean length and width changed by 24.1 µm (4.6%) and 1.4 µm (2.0%), with a smaller shift in concentration. Assimilation adjusts only the shape of the seed distribution, not its mass, because the distribution is renormalized to the measured seed mass during simulation.
These small shifts carried most of the predictive gain. The trajectory loss on the original validation experiments fell to 0.176, a sevenfold improvement over conventional modelling. Applying the published empirical law alone with assimilated initial conditions reached 0.195, so most of the improvement came from correcting the initial conditions rather than from the neural rate law. The authors also note that the non-linearity appears where supersaturation is densely sampled, which suggests it is present in the data rather than imposed by the network; the ensemble spread across 20 random initializations is narrow in that region and wider where sampling is sparse.
Gradients make 300-variable profiles tractable
The optimization tasks parametrize the operating condition as 300 temperatures, far more degrees of freedom than finite-difference or non-differentiable workflows can practically handle. A forward-difference gradient over 300 temperatures would require 301 forward solves, roughly 4 minutes per iteration at the reference solver's simulation time and more than a day for a full optimization. With automatic differentiation, the gradient with respect to the entire temperature profile takes about 0.75 seconds.
The authors first validated trajectory recovery on synthetic targets generated from a stepped ramp and a sigmoid, starting from the E5 initial condition; the original operating condition was reliably recovered across a range of initial guesses. They then targeted the held-out experiment E5. With the initial condition fixed at its measured value, the optimized profile departed from the experimental one (mean absolute error 0.11 °C) and the length and width trajectories were not accurately captured, the same failure mode seen in training. Freeing the initial condition recovered a profile closely matching the experiment (mean absolute error 0.06 °C) and the target trajectories, in under 5 minutes. The profile does not match below the saturation temperature because the hybrid model used here has no growth term.
A simulated design cuts aspect ratio by 25%
The second optimization targets a product specification rather than a known answer. Using a combined size-dependent growth and dissolution model, the authors minimize the final aspect ratio, defined as average length L1 divided by average width L2, subject to the constraint that product volume does not fall below seed volume, so no material is lost. They start from the E5 initial condition and a linear ramp that begins in dissolution and ends in growth.
The optimized profile reduces the aspect ratio from 6.50 to 4.88, a 25% reduction from seed to product; the initial guess gave 6.05. The profile keeps a single cycle: it first dissolves at constant supersaturation, with product volume falling to about 35%, then enters a growth region where volume returns to 100%, staying within supersaturation limits and satisfying the volume constraint. The mechanism the authors give is that length growth scales with L1 to the power 1.2 while width growth is size-independent, so dissolving the population first reduces subsequent length growth relative to width growth and redistributes mass toward width as volume is recovered.
This result is a purely model-generated hypothesis. The authors state that experimental validation is beyond the scope of the work, so the 25% figure should be read as a simulated design target, not a measured outcome.
One material system and thirteen experiments
The demonstration covers a single material system and 13 experiments, and the learned rate law is size-independent. The aspect-ratio optimization relies on published size-dependent growth and dissolution kinetics rather than on the learned hybrid law, because manipulating aspect ratio requires both growth and dissolution while the trained hybrid model is purely size-independent dissolution. The trajectory-recovery task also cannot match the experimental profile below the saturation temperature for the same reason.
Within those bounds, the authors argue there is no additional data requirement compared with empirical modelling, the networks remain small and easy to train, and all tasks cost a fraction of conventional approaches. They present the framework as augmenting rather than replacing mechanistic modelling, and as extendable in principle to size-dependent kinetics, other PBE mechanisms such as nucleation, breakage and agglomeration, and other transport equations.
The priority is separating data error from model error
The paper's most consequential claim is not that a neural network can fit crystallization trajectories. It is that improved predictive performance is not sufficient evidence of a transferable physical law. In this system, the version that fit better while assuming exact initial conditions learned kinetics that violated the equilibrium condition and invented a temperature dependence; the version that fit slightly worse on the rate law alone, but assimilated the initial conditions, recovered a law consistent with the physics and cut trajectory error sevenfold.
For readers evaluating similar hybrid or machine-learning transport models, the practical split is therefore between workflows that treat measured initial conditions as fixed and those that let them be inferred. The paper shows the first can hide a data problem inside a plausible-looking rate law, and that the second is computationally feasible once gradients flow through the solver. The 300-variable optimization results follow from the same capability, but the aspect-ratio design remains a simulation awaiting experiment, whereas the initial-condition diagnosis is supported by the 13 real trajectories.
Sources used in this answer
Differentiable Hybrid Modelling for Learning and Optimising Chemical Transport Processes from Experimental Data
Embedding a neural network in a differentiable GPU-accelerated PBE solver and assimilating measured initial conditions recovered a physically consistent non-linear dissolution law and cut trajectory error sevenfold; the same framework optimized a 300-variable temperature profile that reduced simulated crystal aspect ratio by 25% at constant mass.