In-silico Quality-by-Design characterization of a cooling-crystallization unit operation using an agentic cloud-science environment

Zoe Huiying Tian¹

¹ Independent researcher, Chicago, IL, USA. Correspondence: zoetian@engineering.upenn.edu.

Preprint / independent study — not peer reviewed. All data in this work are synthetic, generated by a mechanistic simulator written for this study. No proprietary, company, or experimental data were used, referenced, or derived. This manuscript is a personal learning project; it makes a methodological contribution (how to run a full ICH Q8/Q11 characterization entirely in-silico using agentic tooling), not a claim about any real drug substance.

Abstract

Quality-by-Design (QbD) as codified in ICH Q8(R2) and Q11 asks the developer to identify critical process parameters (CPPs), relate them to critical quality attributes (CQAs), and define a design space within which quality is assured [1,4,5]. In practice this requires many experiments, and the underlying structure — where the process changes mechanism, where quoted ranges over- or under-state the truly safe region — is often invisible to the empirical models that summarize the campaign. Here I report a fully in-silico ("dry-run") characterization of a model cooling-crystallization unit operation (designated CR-201), carried out end-to-end inside an agentic cloud-science environment (Claude Science, Anthropic) that couples a large-language-model agent to a live scientific compute kernel [50,51,52]. A first-principles crystallization model — population balance by the method of characteristics, classical nucleation and size-dependent growth, agitation-dependent attrition and suspension, and a PSD-dependent isolation yield [6,7,8,9,17,18,22] — was constructed and subjected to a falsifiable test suite (34 assertions) before any characterization was run. A 54-batch design-of-experiments campaign [28,29,30], a variance-based (Sobol) global sensitivity analysis [23,24,25,26,27], a grid-resolved design space with Proven Acceptable Ranges (PARs) and a Normal Operating Range (NOR), and two empirical surrogate models (PLS regression [31,32,33] and a gradient-boosted design-space classifier [34,35,36]) were then produced by the agent.

Three results are worth stating up front. First, all four studied parameters are CPPs under a total-effect-index threshold of 0.10, but the ranking is CQA-dependent: cooling rate dominates purity (S_T = 0.64), agitation dominates yield (0.59), seed loading dominates span (0.41), seed-point temperature dominates d₅₀ (0.34). Second, the PAR rectangle is only 81% releasable: 19% of the box formed by the four individually-proven ranges lies outside the joint design space, a direct consequence of quoting rectangular ranges over a non-rectangular region [40,41,42]. Third, the empirical surrogates disagree with the generating mechanism most often exactly where it matters — near the metastable limit, the PLS release call is wrong on 23.6% of fresh batches versus 20.7% away from it, and the classifier 17.6% versus 12.2%. Because the mechanism is known (it was authored), these disagreements can be measured rather than assumed — an audit that a real campaign, lacking ground truth, cannot perform. Fourth, a scale-up analysis across 2 L, 200 L, and 2000 L shows that the release outcome is decided by the agitation transfer basis, not the recipe: transferring on constant rpm collapses the product above 2 L (fines rising from 0.5 % to 58 %), whereas transferring on constant power per unit volume reproduces the bench batch exactly at every scale. The study is offered as a worked template for in-silico process understanding and as a demonstration that agentic cloud-science tooling can execute the full QbD workflow — modelling, DoE, sensitivity, design space, surrogate audit, and scale-up — as a single reproducible thread.

Keywords: Quality by Design; cooling crystallization; population balance; design space; global sensitivity analysis; partial least squares; ICH Q8; ICH Q11; agentic AI; in-silico process development.

1. Introduction

The regulatory framework for pharmaceutical process development shifted decisively with ICH Q8(R2), which reframed development around understanding: the systematic identification of the parameters that drive product quality and the definition of a multidimensional design space — "the multidimensional combination and interaction of input variables … that have been demonstrated to provide assurance of quality" [1]. Q9 and Q10 place this in a quality-risk-management and quality-system context [2,3]; Q11 extends the same logic to the manufacture of drug substances, where crystallization is very often the CQA-defining step [4,5,47]. Within that framework a handful of quantities recur: critical process parameters (CPPs), the knobs an operator can set; critical quality attributes (CQAs), the product properties a patient ultimately depends on; Proven Acceptable Ranges (PARs), the range over which a single parameter has been shown acceptable with the others fixed; and the Normal Operating Range (NOR), the tighter box a plant actually runs in [1,5].

Crystallization is the archetypal hard case. The particle-size distribution (PSD), polymorphic form, purity, and yield of an isolated solid are set by a competition between nucleation and growth that is itself governed by the supersaturation trajectory — and that trajectory is a nonlinear function of cooling rate, seeding, temperature, and hydrodynamics [6,7,12,13]. Small changes in how one approaches the metastable zone can flip a process between a seeded, size-controlled regime and an uncontrolled nucleation burst that produces fines and bimodal distributions [9,10,11,48]. The mechanistic backbone for describing this — the population balance equation coupled to nucleation and growth kinetics — is well established [8,17,18,19] and has been used extensively for model-based crystallization control [13,14,15,16,44]. Yet in day-to-day development the design space is usually mapped empirically, by fitting a response surface or a classifier to a DoE campaign [28,29,40]. Those surrogates are indispensable for interpolation, but they cannot tell the developer where they themselves break down, because the truth they are approximating is unknown.

This paper takes the opposite starting point. Rather than fit a surrogate to unknown physics, I author the physics as a mechanistic simulator, treat it as the ground truth of a synthetic unit operation, and then run the entire QbD workflow against it. Because the mechanism is known, every downstream artifact — the sensitivity ranking, the design-space boundary, the PAR/NOR boxes, and the empirical surrogates — can be checked against it. The design space is not estimated; it is resolved. The surrogate error is not assumed; it is measured, and located. This is a "dry run" in the literal sense: an end-to-end rehearsal of process characterization with the answer key in hand.

The second contribution is methodological in a different sense. The whole study — building the model, writing and running its test suite, designing and executing the campaign, the sensitivity and design-space analyses, the surrogate fits, and the figures — was carried out inside an agentic cloud-science environment (Claude Science), in which a large-language-model agent [50] drives a persistent scientific compute kernel with package management, a filesystem, and artifact tracking. Recent work has shown large language models acting as autonomous or semi-autonomous experimental agents in chemistry [51,52] and has argued for AI as a first-class instrument of scientific discovery [53]. The present study is a controlled demonstration of that idea in a domain where the ground truth is fully specified, so the agent's output can be judged against a known answer. I report both the science and, in the Methods and Discussion, an honest account of what the agent found — including two modelling errors it caught and corrected mid-study — because the reproducibility of an agentic workflow is itself part of the result.

The remainder of the paper is organized conventionally. Section 2 describes the agentic environment, the mechanistic model and its verification, the DoE campaign, the sensitivity and design-space methods, and the surrogate models. Section 3 reports results. Section 4 discusses scale-up implications, limitations of a synthetic study, and the role of agentic tooling. Section 5 concludes.

2. Methods

2.1 Agentic cloud-science environment

All computation was performed inside Claude Science (Anthropic), an environment in which a large-language-model agent [50] is coupled to a persistent Python compute kernel with NumPy/SciPy/scikit-learn/Matplotlib [36,37,38,39], on-demand package and environment management, a shared filesystem, and an artifact store that versions every output with its generating code. The agent authored the model and analysis scripts, executed them, inspected intermediate results, produced the figures, and versioned each deliverable; the human author set the goals, the framing, and the acceptance criteria, and reviewed the output. This division — human intent, machine execution, human review — mirrors the way autonomous LLM agents have recently been deployed for chemical synthesis planning and execution [51,52], but here it is applied to process characterization against a known ground truth so that the agent's conclusions are independently checkable. The complete code (src/cr201_model.py, the analysis scripts under scripts/, and the test suite) and all synthetic data are retained as versioned artifacts.

2.2 Mechanistic crystallization model

CR-201 is a seeded batch cooling crystallizer for a small-molecule API. The model integrates one batch by the method of characteristics applied to the population balance equation [8,17,18]: rather than discretize the PSD on a fixed grid [19], each nucleation event spawns a Lagrangian cohort that is then advected in size space by the growth rate, which keeps the numerically-delicate nucleation front sharp. The state variables are the liquid-phase solute concentration and a growing list of (size, count) cohorts.

Supersaturation. The solubility follows a van 't Hoff-type curve; the driving force is the supersaturation ratio S = C/C*(T). Cooling from 60 °C to 8 °C over a seeded batch (initial load 250 g solute · kg⁻¹ solvent, seed point 56 °C, 30 min endpoint hold) generates the S-trajectory that everything else responds to [6,7,12].

Nucleation. Primary (homogeneous/heterogeneous) nucleation follows a classical power law in supersaturation with a sharp onset near the metastable limit, and secondary nucleation is driven by existing crystal surface and by agitation-induced attrition [9,10,11]. The metastable zone width is cooling-rate dependent — faster cooling pushes the metastable ratio S_meta higher before nucleation relieves it [10,12] — so a batch that cools too fast overshoots into a nucleation burst.

Growth. Crystals grow by a size-independent (McCabe ΔL) law modulated by supersaturation; growth consumes solute and thereby couples back to S. Ostwald ripening of the finest cohorts is included following Lifshitz–Slyozov–Wagner theory [20,21].

Hydrodynamics. Agitation enters twice. Below the just-suspended speed (Zwietering correlation [22]) part of the crystal surface is unavailable to the supersaturated liquor, suppressing growth; above a reference speed, power-per-volume (~rpm³) drives attrition and secondary nucleation [46]. This gives agitation a genuine interior optimum — too low starves growth, too high shatters crystals into fines.

Isolation. The isolated yield is not the equilibrium yield. Crystals finer than a filter cut (20 µm, sharpness 5) are incompletely retained on the filter cloth and are partially lost to the wash (30% loss of the sub-cut mass), so a batch that nucleates a fine tail loses yield on isolation even though it reached the same thermodynamic endpoint. This PSD-dependent isolation loss (added during verification; see §2.3) is what makes yield sensitive to how one crystallizes and not merely to where one stops cooling.

CPPs, CQAs, and the nominal recipe. The four CPPs are cooling rate, seed loading, seed-point temperature, and agitation. The four CQAs and their release specifications are d₅₀ ∈ [45, 95] µm, span ∈ [0, 1.6], yield ∈ [87, 102] %, and purity ∈ [99, 100] %. At the nominal recipe (10.4 °C·hr⁻¹, 2.0 % w/w seed, 56 °C seed point, 85 rpm) the model returns d₅₀ = 68.0 µm, span = 0.680, yield = 89.33 %, purity = 99.88 %, with peak supersaturation S_max = 1.207 comfortably below the metastable ratio S_meta = 1.351.

2.3 Model verification

Before any characterization was run, the model was subjected to a suite of 34 falsifiable assertions (tests/test_cr201_model.py) encoding physical expectations rather than fitted numbers: solute mass is conserved to tolerance across the integration; faster cooling raises peak supersaturation and, past a point, triggers a nucleation burst; more seed surface lowers d₅₀ when the batch is seed-controlled; agitation has an interior optimum for yield; the isolated yield depends on the PSD and not only on the thermodynamic endpoint; and the filter retains coarse product essentially quantitatively. The suite is part of the deliverable and is intended to be run on every change to the model.

Two verification findings changed the model and are reported for transparency. First, an early version made yield purely thermodynamic — every batch cooled to the same 8 °C endpoint, so cooling rate had exactly zero effect on yield, which is physically wrong for a process that loses fines on isolation. This motivated the PSD-dependent isolation loss described above; with it, cooling rate and agitation both act on yield through the fine tail. Second, a seed-control diagnostic was initially computed on the number fraction of nucleated crystals, which flagged the nominal recipe as uncontrolled because nuclei dominate any population by count (83 %) while carrying < 0.5 % of the mass; re-deriving the diagnostic on a mass basis resolved it. Both errors were caught by the agent during verification, not after publication — which is the point of writing the tests first.

2.4 Design-of-experiments campaign

A 54-batch campaign was generated in three blocks [28,30]: a space-filling screening block over the full studied ranges, an optimization block concentrated around the nominal recipe, and a challenge block placed deliberately at aggressive combinations. Setpoints and as-executed values (with small realization noise) were recorded, together with the "true" CQAs from the mechanism and the "measured" CQAs with analytical scatter added, and the release disposition. Of the 54 batches, 21 are in specification (block 1: 12/30, block 2: 9/18, block 3: 0/6) — the challenge block failing wholesale by construction, which is what a challenge block is for.

2.5 Global sensitivity analysis

Variance-based (Sobol) global sensitivity indices [23,24] were estimated with the Saltelli sampling scheme [26] as implemented in SALib [27], propagating uniform distributions over the four CPPs through the mechanism (6 144 model evaluations per index set) and computing first-order (S₁) and total-effect (S_T) indices with bootstrap confidence intervals [25]. A parameter is designated a CPP if S_T ≥ 0.10 on any CQA. The interaction contribution is reported as S_T − S₁.

2.6 Design space, PARs, and NOR

The design space was resolved on a dense grid over each CPP pair, classifying every point by whether all four CQAs release. PARs were derived by varying one parameter at a time from the nominal recipe until the first CQA leaves specification, with a 5 % edge margin so that no range boundary is set by analytical scatter rather than by the mechanism. The NOR was taken as the largest hyper-rectangle, centered on the nominal recipe, that is 100 % releasable at every corner. To quantify the gap between one-at-a-time ranges and the joint region [40,41,42], 4 000 points were sampled uniformly inside the PAR rectangle and the fraction failing release ("PAR-box reliability") was computed, with the failing points attributed to the CQA combinations responsible.

2.7 Empirical surrogate models

Two surrogates were fit to the 54-batch campaign for comparison against the mechanism. A partial-least-squares (PLS) regression [31,32,33] predicted all four CQAs jointly from the four (standardized) CPPs; the number of latent components was chosen by repeated cross-validation on Q². A gradient-boosted classifier (histogram GBDT [34,35,36]) predicted the binary release disposition. Both were then evaluated on a fresh mechanistic sample of 600 batches that neither model had seen. Because the mechanism generated that holdout, each surrogate's release call could be compared directly against the true disposition, and the disagreement rate was stratified by proximity to the metastable limit (supersaturation headroom S_meta − S_max below vs. above 0.05). scikit-learn [36], NumPy [37], SciPy [38], and Matplotlib [39] were used throughout.

3. Results

3.1 All four parameters are critical, but for different attributes

Figure 1 reports the Sobol decomposition. Under the S_T ≥ 0.10 convention, every studied parameter is a CPP, but no single parameter dominates the process; which one leads depends on the CQA. Cooling rate has the largest single effect anywhere in the process — S_T = 0.64 on purity — because it sets how far the batch overshoots the metastable limit and therefore how much impurity is occluded during a nucleation burst [9,48]. Agitation is the dominant driver of yield (S_T = 0.59), acting through the fine tail that isolation loses. Seed loading dominates span (0.41), and seed-point temperature dominates d₅₀ (0.34).

Between 10 % and 16 % of the variance in d₅₀, span, and yield is interaction (S_T − S₁), which is invisible to a main-effects screening design and is the reason one-factor-at-a-time development under-characterizes this process. The most important interaction is between agitation and seed loading: increasing agitation raises the critical seed loading — the smallest seed charge that keeps the batch seed-controlled — from ~1.0 % w/w at 60 rpm to ~2.75 % at 150 rpm, because more vigorous agitation generates attrition nuclei that compete with the seed surface for supersaturation. (An earlier reading of this interaction as a "slope reversal" in d₅₀ vs. seed loading was found to be an artifact of fitting a straight line across a discontinuity, and was corrected to the critical-loading description given here.) The practical consequence, flagged for automation in the dashboard companion to this study, is that a seed-loading alarm limit cannot be a constant — it must be conditioned on the current agitation setpoint.

Figure 1. Variance-based (Sobol) sensitivity. Top row: stacked first-order (S₁) and interaction (S_T − S₁) indices per CPP for each CQA, with total-effect confidence intervals; the dashed line marks the S_T = 0.10 CPP threshold. Bottom: the agitation × seed-loading interaction, showing the critical seed loading shifting with agitation. 6 144 model evaluations per index set. 100 % synthetic data.

3.2 The design space is not a rectangle, and quoting it as one costs 19 %

Figure 2 shows the resolved design space and the derived ranges. The Proven Acceptable Ranges, each with a 5 % edge margin, are: cooling rate 4–26 °C·hr⁻¹ (first failure at the high edge is yield), seed loading 0.5–6.0 % w/w (no failure within the studied range), seed-point temperature 52.5–58.0 °C (low-edge failure is d₅₀), and agitation 50–127 rpm (high-edge failure is d₅₀). The Normal Operating Range is the largest fully-releasable box, at 39 % of the PAR half-width in every dimension: cooling rate 7.9–12.9, seed 1.41–2.59 %, seed temp 55.2–56.8 °C, agitation 71–99 rpm.

The result I would lead a filing with is the PAR-box reliability of 81 %: sampling 4 000 points uniformly inside the rectangle formed by the four PARs, 19 % fail release [40,41,42]. The failures decompose as 9.1 % on d₅₀ and yield jointly, 6.9 % on d₅₀ alone, and 3.0 % on yield alone. This gap is not a modelling error; it is the geometric consequence of quoting rectangular ranges over a region whose cooling-rate/agitation boundary is diagonal, so the corner where both parameters sit at their individually-proven limits lies outside the joint space. Every filing that reports PARs without the joint region carries this gap silently; here it can be measured because the space is resolved rather than inferred. The NOR is the honest alternative — narrower, but with no unreleasable corners.

A separate check confirmed that the specifications and the mechanism agree along the PAR axes: scanning each PAR from edge to edge, no point is simultaneously in-specification yet mechanistically uncontrolled (fines mass fraction elevated). The release specs are, on these axes, a faithful proxy for the underlying physics — a reassuring result rather than a worrying one.

Figure 2. Design space and acceptable ranges. Left: releasable region over cooling rate × agitation with the diagonal boundary that makes the PAR corner unreleasable; nominal recipe starred. Right: PAR (orange, dashed) and NOR (blue, solid) as fractions of the studied ranges, with edge failure modes. Headline: PAR-box reliability 81 %. 100 % synthetic data.

3.3 The empirical surrogates are weakest exactly where the mechanism is decisive

The PLS regression selected four latent components (Q² plateauing at 0.497) and reproduced the held-out CQAs with R² of 0.70 (d₅₀), 0.69 (span), 0.67 (yield), and 0.71 (purity). Normalized by specification width, its holdout RMSE is small for purity (2.4 %) and yield (9.7 %) but larger for span (12.8 %) and d₅₀ (14.6 %) — the two attributes most affected by the nonlinear nucleation front, which a linear projection cannot represent. The classifier reached ROC-AUC 0.70 (PR-AUC 0.77) in cross-validation and 0.92 on the fresh holdout.

Figure 3 reports the analysis that a real campaign cannot perform. On the fresh 600-batch mechanistic sample, the PLS release call is wrong 22.2 % of the time overall, but 23.6 % near the metastable limit versus 20.7 % away from it; the classifier is wrong 14.8 % overall, 17.6 % near versus 12.2 % away. The degradation toward the boundary is real and consistent across both surrogates, though modest in magnitude (roughly 3–5 percentage points) — I report it at that size rather than inflating it. Its mechanism is visible in the PLS d₅₀ residuals, which fan out as the supersaturation headroom approaches zero: the region where the batch is about to tip from seed-controlled growth into a nucleation burst is exactly the region a smooth surrogate blurs. This is the empirical case for why, in this study, the mechanism and not the fitted surrogate defines the design space — the surrogate is an interpolator that is least trustworthy at the most consequential boundary [43,49].

Figure 3. PLS and classifier versus the mechanism. Panels: PLS predicted vs. actual on the fresh holdout; PLS d₅₀ residual vs. supersaturation headroom (error grows toward the metastable limit); PLS scaled coefficients; classifier P(OOS) against the true boundary; release-call error near vs. away from the limit; PLS RMSE as a fraction of spec width. 100 % synthetic data.

3.4 Scale-up: the same recipe releases or fails depending only on the agitation transfer basis

To test whether the characterization transfers across scale, the mechanistic model was extended with a vessel-volume parameter under geometric similarity (§2.2): the linear scale ratio is s = (V/V_bench)^{1/3}, the just-suspended speed falls as N_js ∝ D^{−0.85} (Zwietering [22]), and power per unit volume rises as P/V ∝ N³D² [46], so agitation's two roles diverge with scale. The bench characterization sits at 2 L; the nominal recipe was transferred to 200 L and 2000 L under three transfer bases (Figure 4, Table 1).

The result is unambiguous. Transferring on constant rpm fails at both larger scales. Holding 85 rpm from 2 L to 2000 L drives d₅₀ from 68 µm down to 26 µm (200 L) and 19 µm (2000 L), span from 0.68 to 1.48 and 1.69, yield from 89.3 % to 82.6 % and 75.4 %, and the fines mass fraction from 0.5 % to 21 % and 58 % — an uncontrolled, fines-dominated product at plant scale, because constant rpm at larger diameter means higher power/volume and therefore runaway attrition. Transferring on constant power/volume reproduces the bench batch exactly at every scale (rpm falling 85 → 30.5 → 18.3), because it holds attrition intensity constant, which is the dominant mechanism here. Constant tip speed over-corrects — it under-agitates the larger vessels (d₅₀ drifts up, span narrows) — but stays releasable.

Equivalently, the releasable agitation window shrinks with scale — 27–120 rpm at 2 L, 8–45 rpm at 200 L, 8–27 rpm at 2000 L — and the constant-power/volume setpoint stays near its center at every scale while the bench 85 rpm falls outside it above 2 L. The window, not a single setpoint, is what should be transferred. Because the mechanism is authored, these are exact statements about the model rather than extrapolations from a surrogate fit at one scale — which, being fit at 2 L, contains no information about 200 L or 2000 L at all.

Figure 4. CR-201 scale-up from 2 L to 2000 L. Left: transferring the recipe on constant rpm collapses the CQAs (shown as fraction of spec width; green band = in spec) above 2 L. Center: the releasable agitation window (shaded) shrinks with scale; the constant-power/volume setpoint (blue) tracks it while the held-constant bench 85 rpm (red dashed) leaves it. Right: the mechanism of failure — constant rpm means rising power/volume, which shatters crystals into fines. Geometric similarity; N_js ∝ D^{−0.85}, P/V ∝ N³D². 100 % synthetic data.

Table 1. Nominal recipe transferred across scale under three agitation bases.

4. Discussion

4.1 What the "answer key" buys

The single feature that distinguishes this study from a conventional characterization is that the ground truth is authored, not estimated. That changes what each downstream artifact means. The sensitivity ranking is not an estimate whose error must be argued from replicate scatter; it is a variance decomposition of the true response. The design space is resolved on a grid rather than fit and extrapolated, so the 81 % PAR-box reliability is an exact geometric statement about the region, not an inference. And the surrogate audit in §3.3 is only possible because the disposition of the 600 holdout batches is known — a real campaign, which never sees the truth, can report a surrogate's cross-validated accuracy but cannot say where the surrogate is wrong relative to the physics. The general lesson is methodological: an in-silico dry run with a known mechanism is a way to calibrate one's expectations of the empirical tools before spending them on the bench [43,49,53].

4.2 Scale-up: the transfer basis, not the recipe, decides the outcome

The scale-up analysis (§3.4) makes the discussion concrete. Because agitation enters the mechanism through both suspension (the just-suspended speed, which falls with vessel diameter [22]) and attrition (power per unit volume, which rises [46]), the two effects that are co-calibrated at bench scale diverge on scale-up, and it is that divergence — not the recipe — that determines whether the batch releases. Transferring the nominal recipe on constant rpm fails at both 200 L and 2000 L; transferring on constant power/volume reproduces the bench batch exactly at every scale. The mechanistic model already contains the variables scale-up risk depends on, whereas an empirical surrogate fit at one scale carries no information about any other. The practical reading is that a scale-up protocol must transfer on a hydrodynamic invariant and re-derive the agitation window (and the agitation-dependent critical seed loading of §3.1) at each scale, rather than hold rpm fixed.

4.3 Limitations

This is a synthetic study, and its limits should be stated plainly. The mechanism is a plausible crystallization model, not a fitted representation of any real compound; its rate constants and the isolation model were chosen to produce qualitatively correct, internally-consistent behavior, not to match experimental data. Consequently the numerical results (specific S_T values, the 81 % figure, the surrogate error rates) characterize this model, and their value is as a worked demonstration of the workflow rather than as claims about a real drug substance. No polymorphic transformation kinetics, impurity-specific rejection, or spatial inhomogeneity in the vessel are modelled beyond the lumped terms described. The surrogates were fit to a deliberately small (54-batch) campaign to mimic a realistic development budget; larger campaigns would narrow their error but would not remove the boundary-localized degradation, which is structural. Finally, the two modelling errors caught during verification (§2.3) are a reminder that a mechanistic model is only as good as the physical assertions it is tested against — the test suite is not an appendix to this work but a load-bearing part of it.

4.4 Agentic cloud-science as the instrument

The workflow reported here — author a mechanism, verify it against falsifiable assertions, run a DoE campaign, decompose the variance, resolve the design space, fit and audit surrogates, and produce publication figures — was executed as a single reproducible thread by an LLM agent driving a live compute kernel [50,51,52]. Three observations follow. First, the agent's usefulness was greatest precisely where a known ground truth allowed its output to be checked: it caught and corrected two of its own modelling errors during verification, which is the behavior one wants and which is only observable because the tests encode physics the agent could be held to. Second, the artifact-versioning discipline of the environment — every figure and table carries its generating code — makes the study reproducible in a stronger sense than a methods section usually affords. Third, and more speculatively, the same loop that made this dry run cheap (state a physical expectation, let the agent implement and test it, inspect where it fails) is a plausible template for using agentic tooling to pre-plan real campaigns: characterize the mechanism in-silico, find where the empirical tools will be weakest, and spend bench experiments where the dry run says they matter most [53]. The claim is not that the agent discovered new crystallization science; it is that the agent executed a rigorous, checkable QbD workflow end-to-end, which is itself a demonstration worth reporting.

5. Conclusion

A cooling-crystallization unit operation was characterized end-to-end in-silico, using an agentic cloud-science environment, against a mechanistic model treated as ground truth. All four studied parameters are critical, but for different quality attributes; the Proven Acceptable Range rectangle is only 81 % releasable, exposing the cost of quoting one-at-a-time ranges over a non-rectangular design space; and two empirical surrogates were shown to disagree with the generating mechanism most often near the metastable limit, the boundary that matters most and that a smooth model represents least well. Because the mechanism was authored, each of these results is checkable rather than assumed — which is the entire point of a dry run. The study is offered as a reproducible template for in-silico process understanding and as a demonstration that agentic tooling can carry the full ICH Q8/Q11 characterization workflow as one thread. All data are synthetic; the code, tests, and data are retained as versioned artifacts.

Data and code availability

All code (src/cr201_model.py, analysis scripts, and the 34-assertion test suite) and all synthetic data and figures are retained as versioned artifacts in the project repository. No proprietary or experimental data were used.

Author contributions

Z.T. conceived the study, set its framing and acceptance criteria, and reviewed all output. The modelling, analysis, and figure generation were executed by an LLM agent (Claude Science, Anthropic) under the author's direction; the author is responsible for the content.

Conflicts of interest

None declared. This is an independent, non-commercial study.

A note on the use of AI

This study was produced using an agentic large-language-model environment as its primary computational instrument, as described in §2.1 and §4.4. The author directed the work and takes responsibility for its content; the transparent account of the agent's role, including the errors it caught and corrected, is provided in the spirit of the reproducibility this work advocates.

References

[1] International Council for Harmonisation. ICH Q8(R2): Pharmaceutical Development. ICH, 2009.

[2] International Council for Harmonisation. ICH Q9: Quality Risk Management. ICH, 2005.

[3] International Council for Harmonisation. ICH Q10: Pharmaceutical Quality System. ICH, 2008.

[4] International Council for Harmonisation. ICH Q11: Development and Manufacture of Drug Substances. ICH, 2012.

[5] Yu, L. X., et al. "Understanding pharmaceutical quality by design." The AAPS Journal 16.4 (2014): 771–783.

[6] Mullin, J. W. Crystallization. 4th ed. Butterworth-Heinemann, 2001.

[7] Myerson, A. S., ed. Handbook of Industrial Crystallization. 2nd ed. Butterworth-Heinemann, 2002.

[8] Randolph, A. D., and Larson, M. A. Theory of Particulate Processes. 2nd ed. Academic Press, 1988.

[9] Kashchiev, D. Nucleation: Basic Theory with Applications. Butterworth-Heinemann, 2000.

[10] Nývlt, J. "Kinetics of nucleation in solutions." Journal of Crystal Growth 3–4 (1968): 377–383.

[11] Mersmann, A., and Bartosch, K. "How to predict the metastable zone width." Journal of Crystal Growth 183.1–2 (1998): 240–250.

[12] Sangwal, K. "Novel approach to analyze metastable zone width." Journal of Crystal Growth 311.21 (2009): 4050–4061.

[13] Fujiwara, M., et al. "First-principles and direct design approaches for the control of pharmaceutical crystallization." Journal of Process Control 15.5 (2005): 493–504.

[14] Nagy, Z. K., and Braatz, R. D. "Advances and new directions in crystallization control." Annual Review of Chemical and Biomolecular Engineering 3 (2012): 55–75.

[15] Braatz, R. D. "Advanced control of crystallization processes." Annual Reviews in Control 26.1 (2002): 87–99.

[16] Rawlings, J. B., Miller, S. M., and Witkowski, W. R. "Model identification and control of solution crystallization processes." Industrial & Engineering Chemistry Research 32.7 (1993): 1275–1296.

[17] Hulburt, H. M., and Katz, S. "Some problems in particle technology: a statistical mechanical formulation." Chemical Engineering Science 19.8 (1964): 555–574.

[18] Ramkrishna, D. Population Balances: Theory and Applications to Particulate Systems in Engineering. Academic Press, 2000.

[19] Kumar, S., and Ramkrishna, D. "On the solution of population balance equations by discretization." Chemical Engineering Science 51.8 (1996): 1311–1332.

[20] Lifshitz, I. M., and Slyozov, V. V. "The kinetics of precipitation from supersaturated solid solutions." Journal of Physics and Chemistry of Solids 19.1–2 (1961): 35–50.

[21] Wagner, C. "Theorie der Alterung von Niederschlägen durch Umlösen." Zeitschrift für Elektrochemie 65.7–8 (1961): 581–591.

[22] Zwietering, T. N. "Suspending of solid particles in liquid by agitators." Chemical Engineering Science 8.3–4 (1958): 244–253.

[23] Sobol', I. M. "Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates." Mathematics and Computers in Simulation 55.1–3 (2001): 271–280.

[24] Homma, T., and Saltelli, A. "Importance measures in global sensitivity analysis of nonlinear models." Reliability Engineering & System Safety 52.1 (1996): 1–17.

[25] Saltelli, A., et al. Global Sensitivity Analysis: The Primer. Wiley, 2008.

[26] Saltelli, A., et al. "Variance based sensitivity analysis of model output." Computer Physics Communications 181.2 (2010): 259–270.

[27] Herman, J., and Usher, W. "SALib: an open-source Python library for sensitivity analysis." Journal of Open Source Software 2.9 (2017): 97.

[28] Montgomery, D. C. Design and Analysis of Experiments. 8th ed. Wiley, 2013.

[29] Myers, R. H., Montgomery, D. C., and Anderson-Cook, C. M. Response Surface Methodology. 4th ed. Wiley, 2016.

[30] Box, G. E. P., Hunter, J. S., and Hunter, W. G. Statistics for Experimenters. 2nd ed. Wiley, 2005.

[31] Wold, S., Sjöström, M., and Eriksson, L. "PLS-regression: a basic tool of chemometrics." Chemometrics and Intelligent Laboratory Systems 58.2 (2001): 109–130.

[32] Geladi, P., and Kowalski, B. R. "Partial least-squares regression: a tutorial." Analytica Chimica Acta 185 (1986): 1–17.

[33] Höskuldsson, A. "PLS regression methods." Journal of Chemometrics 2.3 (1988): 211–228.

[34] Friedman, J. H. "Greedy function approximation: a gradient boosting machine." Annals of Statistics 29.5 (2001): 1189–1232.

[35] Ke, G., et al. "LightGBM: a highly efficient gradient boosting decision tree." Advances in Neural Information Processing Systems 30 (2017).

[36] Pedregosa, F., et al. "Scikit-learn: machine learning in Python." Journal of Machine Learning Research 12 (2011): 2825–2830.

[37] Harris, C. R., et al. "Array programming with NumPy." Nature 585 (2020): 357–362.

[38] Virtanen, P., et al. "SciPy 1.0: fundamental algorithms for scientific computing in Python." Nature Methods 17 (2020): 261–272.

[39] Hunter, J. D. "Matplotlib: a 2D graphics environment." Computing in Science & Engineering 9.3 (2007): 90–95.

[40] Peterson, J. J. "A Bayesian approach to the ICH Q8 definition of design space." Journal of Biopharmaceutical Statistics 18.5 (2008): 959–975.

[41] del Castillo, E. Process Optimization: A Statistical Approach. Springer, 2007.

[42] Lebrun, P., et al. "Design space approach in the optimization of the manufacturing process." European Journal of Pharmaceutics and Biopharmaceutics 74.3 (2010): 483–496.

[43] Kennedy, M. C., and O'Hagan, A. "Bayesian calibration of computer models." Journal of the Royal Statistical Society: Series B 63.3 (2001): 425–464.

[44] Nagy, Z. K. "Model based robust control approach for batch crystallization product design." Computers & Chemical Engineering 33.10 (2009): 1685–1691.

[45] Kitamura, M. "Controlling factors and mechanism of polymorphic crystallization." Crystal Growth & Design 4.6 (2004): 1153–1159.

[46] Gahn, C., and Mersmann, A. "Brittle fracture in crystallization processes." Chemical Engineering Science 54.9 (1999): 1273–1292.

[47] Variankaval, N., Cote, A. S., and Doherty, M. F. "From form to function: crystallization of active pharmaceutical ingredients." AIChE Journal 54.7 (2008): 1682–1688.

[48] Kim, S., et al. "Control of the particle properties of a drug substance by crystallization engineering." Organic Process Research & Development 9.6 (2005): 894–901.

[49] Sacks, J., et al. "Design and analysis of computer experiments." Statistical Science 4.4 (1989): 409–423.

[50] Vaswani, A., et al. "Attention is all you need." Advances in Neural Information Processing Systems 30 (2017).

[51] Boiko, D. A., et al. "Autonomous chemical research with large language models." Nature 624 (2023): 570–578.

[52] Bran, A. M., et al. "Augmenting large language models with chemistry tools." Nature Machine Intelligence 6 (2024): 525–535.

[53] Wang, H., et al. "Scientific discovery in the age of artificial intelligence." Nature 620 (2023): 47–60.

Manuscript prepared as an independent preprint. Figures 1–3 and all tabulated values are generated from the synthetic mechanistic model and its analysis scripts, retained as versioned artifacts alongside this document.