Dynamic language model representations for multi-objective reaction optimisation

main_text.pdf

Introduction

Chemical reaction optimisation is essential to chemical synthesis, and optimising reaction conditions requires navigating vast combinatorial spaces of reaction components such as ligands, catalysts, and bases, while balancing multiple objectives such as yield, selectivity, and cost [1, ?, ?, ?, ?]. This challenge is particularly acute in pharmaceutical process and preclinical chemistry, where identifying robust conditions for active pharmaceutical ingredient (API) synthesis entails additional environmental, health, and safety considerations [?, ?, ?, ?], and in academic settings where enabling novel transformations often requires extensive exploration of many diverse reaction parameters [?, ?].

Machine learning approaches have proved remarkably effective, with Bayesian optimisation emerging as the leading framework for data-efficient reaction optimisation. It has been applied in low-data regimes [?, ?, ?, ?, ?, ?, ?, ?], and more recently, in parallel high-throughput settings [?, ?]. Despite these algorithmic advances, how best to represent reactions for model-driven optimisation remains an open and consequential challenge.

Molecular descriptors have been the dominant featurisation strategy in reaction optimisation and prediction (Figure 1a). Pre-computed descriptor libraries such as Kraken [?] and the COSMO-RS [?] database provide readily available steric and electronic properties for established compound classes and have been widely adopted [?, ?, ?, ?, ?, ?, ?, ?]. Alternatively, tailored and often laborious feature en-

Preliminary work. Under review. gineering workflows, involving computationally expensive density functional theory (DFT) or semi-empirical calculations, conformer ensemble analysis, and reaction-specific mechanistic interrogation, can yield informative descriptors [?, ?, ?, ?, ?, ?, ?, ?, ?].

Despite their widespread use, descriptor-based approaches face fundamental limitations. Selecting which descriptors to compute relies on a priori chemical intuition about which molecular properties govern reactivity, a process that is biased by existing mechanistic understanding and must often be revisited for each new reaction system. More critically, descriptors are often not transferable across chemically distinct reaction components. Steric and electronic

parameters for monodentate phosphines, for example, are only partially applicable to bidentate ligands, and reactions involving chemically heterogeneous components such as inorganic bases and organometallic additives may lack any meaningful shared descriptor space, making feature engineering increasingly prohibitive as chemical diversity grows. One-hot encoding offers a simpler alternative that is universally applicable and has shown competitive performance in some settings [?, ?, ?]. However, it produces sparse, high-dimensional representations that encode no chemical similarity, and its dimensionality scales unfavourably with the number of reaction components [?, ?] (Figure 1).

Advances in large language models (LLMs) have shown strong performance across chemical tasks [?, ?, ?, ?, ?, ?, ?, ?, ?] and demonstrated that their learned embeddings can serve as dense, informative representations for predictive modelling [?, ?, ?, ?]. Recent work has further shown that LLM representations can be adapted to specific optimisation tasks through joint fine-tuning with Gaussian process surrogates, demonstrated for single-objective, sequential optimisation [?]. Crucially, LLMs can encode arbitrarily complex chemical information into fixed-length representations directly from text, offering a naturally unified featurisation method across component classes. These important advances notwithstanding, a general framework for multi-objective reaction optimisation, applicable across chemically heterogeneous reaction components and validated across experimental regimes from low-data to high-throughput settings, has yet to be realised.

Here, we introduce Alice, a multi-objective chemical reaction optimisation framework that uses trainable LLM embeddings as a unified reaction representation (fig:1b). Textual descriptions of reaction conditions are encoded by a low-rank adaptation (LoRA)-fine-tuned language model and passed through reaction objective-specific projection heads into independent Gaussian process surrogates, with the architecture trained end-to-end via joint marginal log-likelihood maximisation. This jointly trained architecture, hereafter referred to as the LLM-GP framework, produces dynamic embeddings that are refined at each optimisation iteration, adapting the reaction representation to the optimisation task at hand. We first benchmark the LLM-GP framework retrospectively against established descriptor libraries and one-hot encoding baselines across sequential low-data optimisation of nickel- and palladium-catalysed cross-couplings and large-batch 96-well plate high-throughput experimentation campaigns, converging on high-performing conditions in fewer experiments across all settings.

We then apply the framework prospectively to two wet-lab campaigns executed on an automated high-throughput platform, each chosen for a form of chemical heterogeneity for which descriptor libraries are not readily available. In a palladium-catalysed cyanation, a single representation spans mono- and bidentate phosphine catalysts alongside additives ranging from elemental zinc to organic and inorganic bases; within this space the framework identified conditions using a lower-hazard cyanide source that gave $>99\%$ conversion and $>99\%$ selectivity, and 94% isolated yield on gram scale. In an asymmetric dynamic kinetic hydrogenation, the framework balanced conversion, diastereomeric excess (de), and enantiomeric excess (ee) across chiral catalysts spanning different metals, donor sets, and coordination geometries, delivering the target syn alcohol in 84% isolated yield and 99.6% ee on scale-up. Both campaigns used only two rounds of 96 experiments, sampling under 3% of their respective

design spaces. Neither required descriptor computation or feature engineering, indicating that effective multi-objective optimisation is achievable in reaction spaces where constructing a shared descriptor representation would itself be a substantial undertaking.

Computational results

We benchmarked the LLM-GP framework against established and readily available featurisation methods in reaction optimisation: DFT descriptors from the Kraken [?] and COSMO RS [?] descriptor libraries, and one-hot encoding baselines. Descriptor coverage is uneven across component classes: Kraken supplies descriptors for the monophosphine ligands and COSMO-RS for the solvents, but no comparable library exists for the bases and precursors, which were therefore represented by one-hot encoding within the descriptor baseline. All three representations were embedded in the same multi-objective Bayesian optimisation loop, using the same acquisition function (qLogNParEGO) and an identical experimental budget (see Methods).

A set of reaction conditions is Pareto optimal when no other condition improves one objective without degrading another; these points define the Pareto front (fig:2a). Optimisation performance was measured by hypervolume, the volume of objective space enclosed by the Pareto-optimal conditions identified so far. A larger hypervolume means the conditions found are both closer to the best achievable trade-offs and better spread across them, and we report it as a percentage of the value attained by the true Pareto front of each dataset. We compared methods by the mean number of iterations required to reach 90% and 95% of the maximum hypervolume, with standard error across 20 random seeds, as these thresholds represent practically relevant levels of convergence. For each dataset the optimisation budget was set by the number of iterations needed for at least one method to reach 95% of the maximum hypervolume on average. We additionally report the percentage of optimisation runs reaching each hypervolume threshold within the allocated experimental budget as a measure of robustness; these results are provided in Extended Data Section A. Among the pre-trained language models and pooling strategies evaluated, T5-base [?] with mean pooling consistently yielded the best optimisation performance and was used throughout (full ablation studies across model architectures and pooling strategies are provided in Extended Data Section B).

Sequential low-data reaction optimisation

We first evaluated our approach in sequential low-data optimisation regimes, where the primary goal is sample efficiency: minimising the number of experiments required to identify high-performing conditions (fig:2b). We used two open-source multi-objective reaction datasets, each comprising experimental yield and selectivity as the optimisation objectives and presenting large combinatorial search spaces with diverse categorical reaction components: a nickel-catalysed Suzuki coupling (384 experiments from a space of 88,000 combinations, varying ligand, Ni-precursor, base, solvent, co-solvent, and temperature) [?], and a

palladium-catalysed Suzuki coupling (404 experiments from 59,040 combinations, varying ligand, Pd-precursor, base, solvent, and co-solvent) [?]. Each campaign was initialised with 5 experiments selected using Sobol sampling to ensure diverse coverage of initial points across the search space (see Methods), after which the Bayesian optimisation loop selected one candidate per iteration. Iteration counts below refer to these optimisation iterations and exclude the initial experiments. On the nickel-catalysed Suzuki coupling, the LLM-GP framework reached the 90% hypervolume threshold in approximately 13 optimisation iterations on average, roughly half the number required by both the DFT descriptor-based ($\sim26$ iterations) and one-hot encoding ($\sim24$ iterations) baselines (Figure 2b). At the more stringent 95% threshold, the gap narrowed, though the framework retained a consistent advantage, converging in approximately 24 iterations compared to $\sim28$ for DFT descriptors and $\sim30$ for one-hot encoding. On the palladium-catalysed Suzuki coupling benchmark, the framework achieved both the 90% and 95% thresholds in approximately 65 iterations on average, while the descriptor-based and one-hot encoding baselines required approximately 80 iterations each to reach the same performance levels (Figure 2b). Across both case studies, the LLM-GP framework consistently reached practical convergence thresholds in fewer experiments than either baseline, improving sample efficiency. Over 20 independent runs initialised with different random seeds, the framework also achieved a higher proportion of successful optimisations reaching both convergence thresholds within the allocated experimental budget, indicating more robust convergence (see Extended Data Section A).

Highly parallel HTE reaction optimisation

Modern high-throughput experimentation (HTE) platforms have combined parallel screening with data-driven optimisation [?, ?], enabling the exploration of large reaction spaces in compressed experimental timescales. We next sought to evaluate our approach in this batched optimisation setting, where reducing the number of HTE plate iterations directly translates to savings in experimental time. We used an open-source palladium-catalysed sulfonamide coupling dataset (a Buchwald–Hartwig-type C–N coupling with 21,648 virtual experiments, varying ligand, Pd-precursor, base, and solvent) [?] as a benchmark, optimising yield and selectivity in parallel batches of 96 experiments to simulate 96-well plate HTE campaigns (Figure 2c). Each campaign was initialised with a Sobol-sampled plate of 96 experiments, after which each iteration selected 96 conditions at once, so one optimisation iteration corresponds to one physical plate. The LLM-GP framework reached practical convergence thresholds in approximately 1.2 plate optimisation iterations beyond initialisation, whereas both the DFT descriptor-based and one-hot encoding baselines required approximately double ($\sim2.4$) to achieve the same thresholds (Figure 2c). The framework also achieved a higher proportion of successful optimisation runs reaching both convergence thresholds, and results for 24- and 48-well plate batch sizes show the same trend (Extended Data Section A). Given that each 96-well plate campaign typically requires approximately one week of experimental and analytical time, this reduction translates directly into meaningful savings in time and resources. These results demonstrate that our approach offers consistent advantages across both sequential low-data regimes, where sample efficiency is crucial, and parallel large-batch settings increasingly adopted in pharmaceutical process development and academic high-throughput screening [?, ?, ?].

Prospective experimental case studies

Building on our retrospective benchmarks across conventional Buchwald–Hartwig and Suzuki–Miyaura datasets, we next sought to apply our framework prospectively to wet-lab reaction optimisation campaigns, moving beyond standard reaction spaces. Here, we targeted transformations containing chemically heterogeneous components that traditional reaction representations struggle to capture. In each campaign, we used our approach to identify high-performing conditions within the search space, iteratively generating ML predictions and conducting the suggested experiments with our HTE robotic platform (see Supplementary Information Section 1 for more details). All experimental data and characterisation of any isolated products are reported in the Supplementary Information.

Case study 1: Palladium-catalysed cyanation

Aryl nitriles are prevalent across pharmaceuticals, agrochemicals, and functional materials [?, ?], serving as compact, metabolically stable hydrogen-bond acceptors and hydroxyl and carboxyl isosteres [?]. Among synthetic routes accessing these motifs, the transition-metal-catalysed cyanation of aryl halides remains one of the most widely adopted and functionally tolerant strategies [?]. While substantial methodology development has produced a broad repertoire of candidate catalysts, additives, and cyanide sources, this very diversity introduces three intersecting sources of chemical heterogeneity that challenge traditional modelling approaches. First, both monodentate and bidentate phosphine ligands are competitive across cyanation substrate classes, yet descriptor libraries developed for monophosphines are not directly transferable to bisphosphines (and vice versa) thus fragmenting the feature space. Second, productive cyanation conditions often require additives spanning chemically disparate classes such as elemental reductants (e.g., Zn), organic bases (e.g., NEt₃), and inorganic salts (e.g., KOAc), which resist straightforward unification within standard chemical feature spaces. Finally, the choice of cyanide source carries practical implications beyond reactivity. Classical sources such as NaCN and Zn(CN)₂ are highly toxic, whereas safer but less reactive alternatives such as K₄[Fe(CN)₆] are increasingly preferred in industrial process settings despite their poor organic solubility, which necessitates biphasic reaction media [?, ?]. Together, these intersecting dimensions of heterogeneity in cyanation reactions pose a distinct challenge for conventional descriptor frameworks, highlighting the value of more flexible, unified representation strategies.

We assembled a design space for the cyanation of methyl 3-bromo-5-fluorobenzoate, comprising 40 commercially available palladium catalysts spanning both mono- and bidentate ligand classes (e.g., [Pd(tBuXPhos)(allyl)]OTf, [Pd(DPPF)Cl$_2$]), three cyanide sources of varying toxicity (NaCN: median lethal dose (mg/kg) $LD_{50} = 3.6$, Zn(CN)$_2$: $LD_{50} = 54$, K$_4$[Fe(CN)$_6$]: $LD_{50} = 3613$), five additive options (Zn, KOAc, NEt$_3$, K$_2$CO$_3$, or none), ten solvents, two co-solvent options, and three temperatures ([]{#fig:3}a). Conditions in which the reaction temperature exceeded the solvent boiling point were removed, giving a final design space of 30,000 conditions. While constructing a bespoke descriptor representation across these heterogeneous components is in principle feasible, doing so would require laborious, trial-and-error feature engineering with no guarantee of success; our approach bypasses this bottleneck by providing a unified, out-of-the-box featurisation directly from textual descriptions, enabling immediate campaign execution without manual descriptor curation. The full search space is provided in Extended Data Figure 7.

30,000 reaction condition combinations

Pd catalyst (0.01 eq)

Cyanide source (1.5 - 3.0 CN− eq)

Additive (0.5 - 1.5 eq)

Solvent (240 µL) / Co-solvent (120 µL)

T (°C), 20 h

19.0 µL

Mono- and bidentate Palladium catalysts (40)

[PdCl₂(PPh₃)₂] [PdCl₂(DPPF)] [PdCl₂(DPPP)] [PdCl₂(AmPhos)] [Pd(allyl)(tBuBrettPhos)]OTf XPhos Pd G3 AmPhos Pd G3 MorDalphos Pd G3

BrettPhos Pd G3 [Pd(allyl)(Xantphos)]Cl [PdCl₂((rac)-BINAP)] [Pd(cinnamyl)Cl(IPr)] [PdCl₂(XantPhos)] [Pd(crotyl)Cl(RuPhos)] Xantphos Pd G3 RuPhos Pd G3

cataCXium A Pd G3 Pd PEPPSI IPent [PdCl₂(DTBPF)] AdBrettPhos Pd G3 [Pd(allyl)Cl((R)-BINAP)] [Pd(allyl)(tBuXPhos)]OTf PtBu₃ Pd G4 MeCgPPh Pd G3

[PdCl₂(DPEPhos)] EPhos Pd G3 BrettPhos Pd G4 CPhos Pd G4 PhCPhos Pd G3 DTBPF Pd G3 PCy₃ Pd G4 TrixiePhos Pd G3

N-Xantphos Pd G3 [HtBuXPhos]₂[Pd₂Cl₆] PAd₃ Pd G3 QPhos Pd G3 sSPhos Pd G2 GPhos Pd G3 Me₄tBuXPhos Pd G3 tBuXPhos Pd G4

Solvents (10)

MeTHF

DMF

tBuOH

NBP

DMC

PhMe

MeCN

CPME

EtOAc

γ-valerolactone

Cyanide sources (3)

NaCN

Zn(CN)₂

K₄[Fe(CN)₆]·3H₂O

Co-solvents (2)

none

water

Additives (5)

none

Zn

KOAc

K₂CO₃

diethylamine

Temperatures (3)

50

75

100

Optimisation was initialised with a Sobol-sampled first plate of 96 experiments to probe the global design space. While the first round identified several high-performing conditions, the top hits relied exclusively on toxic cyanide sources (NaCN and Zn(CN)$_2$), while K$_4$[Fe(CN)$_6$] yielded only mediocre outcomes with either poor conversion or selectivity ([]{#fig:3}b). To test whether high-performing conditions could be identified using a lower-hazard cyanide source, we fixed the cyanide source for round 2 to K$_4$[Fe(CN)$_6$] and re-optimised within this constrained domain. The second optimisation round successfully met this target, identifying multiple conditions achieving $>99\%$ conversion and $>99\%$ selectivity ([]{#fig:3}b). Notably, high-performing hits spanned both monodentate ([Pd(tBuXPhos)(allyl)]OTf) and bidentate ([Pd(Xantphos)(allyl)]Cl) ligands, exemplifying the benefit of jointly searching over both ligand families within a unified representation. To confirm that these milligram-scale HTE hits translate to preparative scale-up, we scaled the lead conditions ([Pd(Xantphos)(allyl)]Cl with K$_4$[Fe(CN)$_6$] · 3H$_2$O in DMC/H$_2$O) to gram scale ([]{#fig:3}c). The reaction proceeded cleanly, delivering the target aryl nitrile with 94% isolated yield (see Supplementary Information Section 2). This scalable, low-hazard protocol was identified in just two iterative campaign rounds (one week per round) totalling 192 experiments, evaluating less than 1% of the 30,000-condition design space.

Case study 2: Asymmetric ketone hydrogenation

Having established the approach on reaction systems with multi-component heterogeneity, we next sought to evaluate its capacity to navigate stereoselective transformations. Catalytic asymmetric hydrogenation of prochiral ketones is among the most widely deployed enantioselective methods in pharmaceutical and agrochemical synthesis [?]. In particular, the dynamic kinetic hydrogenation of α-substituted aryl ketones, in which two contiguous stereocentres are set in a single step to yield up to four potential stereoisomers [?, ?], exemplifies a class of multi-objective optimisation problems that require balancing trade-offs across a Pareto front defined by conversion, diastereomeric excess (de), and enantiomeric excess (ee). Productive catalysts include chiral iridium and ruthenium complexes drawn from a structurally diverse pool of ligand motifs, spanning bidentate phosphines, mixed phosphine–nitrogen donors (P, N and P, N, N), and a variety of associated ancillary ligands and counter-ions. Representing this coordination diversity in a unified descriptor framework is non-trivial; physical featurisation schemes engineered to quantify the 3D chiral pocket of one ligand geometry or metal centre fail to generalise across fundamentally disparate coordination spheres. Jointly balancing three reaction objectives across this heterogeneous chiral space compounds the modelling challenge. We assembled a design space of 8,064 conditions for the asymmetric dynamic kinetic hydrogenation of 1,2-diphenylpropan-1-one towards the syn-(1S, 2R)-1,2-diphenylpropan-1-ol diastereomer, comprising 32 chiral iridium and ruthenium catalysts, 9 organic and inorganic bases, 7 solvents, varied base loading (mol%), and temperature (fig:4a). The full search space is provided in Extended Data Figure 8.

Optimisation was initialised with a Sobol-sampled first HTE plate of 96 experiments spanning the full design space (fig:4b). The first round of experiments revealed that the undesired anti diastereomer was favoured across the majority of conditions evaluated, yielding multiple hits with high anti selectivity (−99% de (syn)). In contrast, accessing the target syn diastereomer was more challenging: only 5 out of 96 experiments were syn selective, and out of the 32 chiral catalysts tested, only a single catalyst pair, [IrClH₂((R/S) – DTB – SpiroPAP – 3-Me)], demonstrated activity towards the syn isomer. The top hit for the desired isomer in Plate 1 achieved quantitative conversion but moderate stereoselectivity (73.5% de (syn) and 86.7% ee (syn)) (fig:4b). In achiral media, enantiomeric catalyst pairs yield mirror-image products, so each training observation was augmented with its reflected counterpart under a sign inversion of the enantiomeric excess (see Supplementary Information Section 3.3).

Guided by the updated LLM representations, the framework refocused optimisation towards syn-producing regions, increasing syn-selective outcomes from 5/96 to 59/96 conditions in Plate 2 (fig:4b). The Pareto front was expanded towards higher stereoselectivity, with the most stereoselective conditions identified reaching 87.1% de (syn) and 93.2% ee (syn). The most balanced lead hit achieved >99% conversion, 80.5% de (syn), and 89.7% ee (syn). Identifying the productive catalyst family did not by itself resolve the optimisation problem: across the reaction conditions employing [IrClH₂((R/S) – DTB – SpiroPAP – 3-Me)], diastereoselectivity spanned −67% to +87.1% de (syn), and in 12 instances its sign inverted upon changing temperature or base loading alone, consistent with the dynamic kinetic nature of the transformation. Systematically removing the highest-performing conditions from the Plate 1 training data did not prevent the model from recovering high-performing syn-selective conditions, indicating that the outcome did not depend on a small number of fortunate initial hits (see Supplementary Information Section 3.4). Translating the lead conditions identified from optimisation to preparative gram scale proceeded smoothly (fig:4c): the lead conditions, [IrClH₂((S) – DTB – SpiroPAP – 3-Me)] with DBU (50 mol%) in tAmOH at 60 °C, delivered the target syn alcohol with the desired enantiomer in 84% isolated yield with 99.6% enantiomeric excess (see Supplementary Information Section 3). This scalable protocol was identified in two iterative campaign rounds (one week per round) using 192 experiments out of 8,064, corresponding to approximately 2.4% of the design space.

Outlook

This work demonstrates that dynamic language model representations provide a general route to multi-objective reaction optimisation, applicable across chemically diverse reaction systems. By encoding reaction components directly from textual descriptions, the framework searches jointly over chemical entities of different classes (e.g., ligand families, additives, and catalysts) within a single representation, without descriptors being selected or computed for each new system. Applied prospectively to two wet-lab campaigns, it identified conditions that translated directly to preparative scale within two rounds of high-throughput experimentation. Engineered descriptors are chemically interpretable, but their selection presupposes an understanding of which molecular properties govern reactivity in the system at hand. Learning the representation dynamically from text removes this requirement, allowing optimisation to begin before such understanding is established, where mechanistic insight could emerge from the conditions identified rather than being needed to find them. By removing expert-guided featurisation as a prerequisite, we anticipate this approach will extend model-guided optimisation to a wider range of chemistry, including biocatalysis and enzymatic reactions, and ultimately enable more efficient chemical processes.

Methods

Metrics for multi-objective optimisation problems

We consider the simultaneous maximisation of $M$ objective functions $y_m:\mathcal{X} \to\mathbb{R}$, $m = 1,\ldots,M$ (e.g., reaction yield and selectivity), over a finite chemical design space $\mathcal{X} = \{x_1,\ldots,x_N\}$, where each $x_i$ is a vector of reaction conditions. The aim is to identify the Pareto optimal set, also referred to as the Pareto front. We say that $x'$ dominates $x$ (written $x' > x$) if $y_m(x') \ge y_m(x)$ for all $m = 1,\ldots,M$ and $y_j(x') > y_j(x)$ for some $j \in\{1,\ldots,M\}$. The Pareto optimal set is then defined as:

We seek to approximate $\mathcal{P}^{*}$ using as few experimental iterations or evaluations as possible with multi-objective Bayesian optimisation. The quality of a candidate Pareto front $\mathcal{P} \subseteq\mathbb{R}^{M}$, obtained from the current set of observa-sequence tions, is measured by its dominated hypervolume (fig:2a) relative to a reference point $\mathbf{r} \in\mathbb{R}^{M}$, chosen such that $r_m < p_m$ for all $\mathbf{p} \in\mathcal{P}$ and $m = 1, \ldots, M$:

where $\operatorname{Vol}(\cdot)$ denotes the Lebesgue measure. For all benchmark datasets, the optimisation objectives are percentages bounded on $[0, 100]$, and we set $\mathbf{r} = (0, 0)$. The same reference point was used for all featurisation methods. Hypervolume is a widely used quality indicator for multi-objective optimisation, rewarding both convergence to the true Pareto front and diversity along it [?, ?, ?, ?]. Importantly, it is monotonic and Pareto-compliant [?, ?], meaning that if one candidate Pareto front strictly dominates another, it will achieve a strictly higher hypervolume. In this work, we report the normalised hypervolume percentage, defined as $\operatorname{HV}(\mathcal{P},\mathbf{r})/\operatorname{HV}(\mathcal{P}^{*},\mathbf{r}) \times100\%$, to evaluate the quality of candidate Pareto fronts $\mathcal{P}$ identified by optimisation algorithms relative to the true Pareto front $\mathcal{P}^{*}$.

Representing reaction conditions

We benchmarked three featurisation approaches for representing reaction conditions.

One-hot encoding. As a universally applicable and inexpensive baseline that has shown competitive performance in prior work [?, ?, ?], each unique categorical variable (e.g., XPhos, dioxane, PhMe) was represented as a binary indicator vector, and all component vectors were concatenated to form the input representation (fig:1a).

Molecular descriptors. Widely adopted as the primary featurisation strategy in reaction optimisation, molecular descriptors from established descriptor libraries were included as a chemically informative baseline [?, ?, ?, ?, ?, ?, ?, ?]. Monophosphine ligands were represented using 190 DFT descriptors from the Kraken library [?], with Principal Component Analysis (PCA) [?] applied to retain components explaining 99% of the variance, reducing dimensionality while preserving information. Solvents were parameterised using four COSMOtherm-derived DFT descriptors from the COSMO-RS database [?]. The remaining categorical variables (e.g., bases and precursors), for which comparable descriptor libraries are not readily available, were represented using one-hot encoding.

Large language model representations. Natural language provides a flexible modality for encoding chemically heterogeneous reaction components without requiring component-class-specific featurisation. Each reaction condition $x_i \in\mathcal{X}$ was represented as a textual prompt of component-value pairs, for example:

Reaction condition: ligand: {ligand name} solvent: {solvent name} precursor: {precursor name} base: {base name}

with one pair per parameter varied in the design space.

A pre-trained language model $h_{\phi}$ maps the tokenised prompt to a sequence of embedding vectors $\mathbf{H}_i = h_{\phi}(x_i) \in\mathbb{R}^{L \times d_{\mathrm{emb}}}$, where $L$ is the sequence length and $d_{\mathrm{emb}}$ is the model’s hidden dimension. A pooling operator $\operatorname{POOL}: \mathbb{R}^{L \times d_{\mathrm{emb}}} \rightarrow\mathbb{R}^{d_{\mathrm{emb}}}$ aggregates the token-level representations into a fixed-dimensional embedding $\mathbf{e}_i = \operatorname{POOL}(\mathbf{H}_i) \in\mathbb{R}^{d_{\mathrm{emb}}}$ to produce a sequence-level representation. We adopted pooling strategies based on each model’s architecture. For encoder-decoder models, only the encoder stack was used, and pooling was applied over its hidden states. The T5 encoder-decoder models (T5-base [?], T5-small, and T5-chem [?]) used mean pooling, which averages hidden states over non-padded tokens. We used last-token pooling for the decoder-only model Qwen2.5 [?], which takes the hidden state at the last non-padding position, as causal attention ensures that only this token has attended to the full input sequence [?]. Mean pooling was also evaluated for Qwen2.5. For the BART-base encoder-decoder model [?], we used mean pooling following the T5 models, and additionally evaluated CLS pooling as BART inherits a <s> classification token. All models were accessed via the Hugging Face Transformers library [?]. Ablation studies over pre-trained language model architectures and pooling strategies across all benchmark datasets are presented in Extended Data Section B.

Overview of optimisation loop

We initialised our Bayesian optimisation workflows using low-discrepancy Sobol sequences [?] to provide broad initial coverage of the design space, which were evaluated to form the initial training data. At each subsequent iteration of the optimisation loop, the surrogate model was fitted to all currently observed data, a batch of candidates was selected by optimising the acquisition function over the remaining design space, and the selected experiments were evaluated and added to the training set.

We used the $\mathrm{qLogNEI}$ $\mathrm{qEGO}$ acquisition function [?, ?, ?], a log acquisition function that conditions on noisy baseline observations and applies Chebyshev scalarisation to reduce the multi-objective problem into single-objective subproblems. For batch acquisition, candidates were selected in a greedy sequential fashion, with each selection conditioned on previously chosen pending points.

For surrogate models using fixed input representations (one-hot encoding and molecular descriptors), we used an optimised Gaussian process (GP) configuration adapted from EDBO+ [?] and Minerva [?] following standard marginal likelihood optimisation. For our approach using an LLM-based deep kernel surrogate with learned input representations, the GP hyperparameters and language model parameters were optimised jointly as described below. All computational experiments were repeated over 20 random seeds to report statistical variation.

LLM-based deep kernel Gaussian process surrogate

The LLM first maps tokenised reaction condition prompts to pooled embeddings $e_i \in\mathbb{R}^{d_{\mathrm{emb}}}$. For each objective $m = 1,\ldots,M$, a trainable projection head $g_{\phi_m}$ transforms the language model embedding $e_i$ to an objective-specific representation $z_i^{(m)} = g_{\phi_m}(e_i)$. A Gaussian process (GP) then models the mapping from each projected representation to its corresponding objective value:

where $f_m$ is the latent function for objective $m$, modelled as a Gaussian process with constant mean function $\mu_m$ and kernel function $k_{\theta_m}$, and $\epsilon_m$ is Gaussian observation noise with variance $\sigma_m^2$. We use the Matérn-5/2 kernel:

where $\sigma_{f,m}^2$ is the signal variance, and $\ell_m$ is the length-scale. Together with the constant mean function $\mu_m$ and noise variance $\sigma_m^2$, these constitute the GP hyperparameters $\theta_m=\{\mu_m,\sigma_{f,m}^2,\ell_m,\sigma_m^2\}$. The kernel matrix for objective $m$ is given by $K_m=k_{\theta_m}\left(Z^{(m)},Z^{(m)}\right)+\sigma_m^2I$, where $Z^{(m)}$ collects the projected embeddings of all observed reaction conditions. The $M$ GPs are queried jointly by the multi-objective acquisition function. The language model $h_\phi$, projection heads $\{g_{\phi_m}\}_{m=1}^{M}$, and GP hyperparameters $\{\theta_m\}_{m=1}^{M}$ are jointly optimised by minimising the sum of negative marginal log-likelihoods:

where the marginal log-likelihood for objective $m$ is:

where $n$ is the number of observations and $\tilde{y}_m=y_m-\mu_m\mathbf{1}$ denotes the targets centred by the constant mean function. GP targets are standardised to zero mean and unit variance for training. Gradients from all $M$ objectives are backpropagated through the kernel and projection heads to the shared LLM parameters $\phi$, where each projection head $g_{\phi_m}$ receives gradients only from its corresponding objective. This enables the LLM to learn representations that are jointly optimised for GP predictive performance across all objectives, while each projection head specialises for its target objective.

Parameter-efficient fine-tuning of language models

The language model embeddings of reaction conditions evolve across BO iterations as the LLM is fine-tuned on accumulating observations. To do so efficiently, we apply parameter-efficient fine-tuning (PEFT) with Low-Rank Adaptation (LoRA) [?]. Rather than updating all of the pre-trained language model’s parameters, we use LoRA to update only a small subset. For a pre-trained weight matrix $\mathbf{W}_0 \in\mathbb{R}^{d \times k}$, LoRA injects trainable low-rank decompositions such that the weight update takes the form:

where $r\ll\min(d,k)$ is the rank and $\alpha$ is a scaling hyperparameter. We use LoRA to target a subset of linear projection matrices within the pre-trained language models, using rank $r=4$ and $\alpha=16$.

Computational implementation

All models, Bayesian optimisation workflows, and computational analyses were implemented in Python 3.10. We used PyTorch [?] (v2.8.0), GPyTorch [?] (v1.14), and BoTorch [?] (v0.13.0) to construct our multi-objective Bayesian optimisation algorithms. Pre-trained language models were loaded using Hugging Face transformers [?] (v4.51.2), with Low-Rank Adaptation (LoRA) implemented using peft (v0.15.1) for parameter-efficient fine-tuning. We used PyTorch Lightning [?] (v2.5.4) for reproducible random seed management and workflow control. Scikit-learn [?] (v1.7.1) and NumPy [?] (v1.26.4) were utilised for machine learning pipelines and numerical data processing, including Principal Component Analysis (PCA) for descriptor dimensionality reduction. Experiment tracking and logging were conducted using Weights & Biases (wandb, v0.21.0) [?]. Matplotlib [?] (v3.10.1) and seaborn [?] (v0.13.2) were used for plotting and visualising all computational results presented in this work. All computations were run on the Roche HPC cluster, using NVIDIA A100 Tensor Core and Blackwell GPUs.

Experimental procedures

High-throughput experimentation procedures, analytical methods, and scale-up syntheses and characterisation for the palladium-catalysed cyanation and asymmetric ketone hydrogenation reactions are provided in the Supplementary Information.

Data availability

All benchmark datasets, reaction condition search spaces, and high-throughput experimentation (HTE) experimental data generated in this study are included in the manuscript, Supplementary Information, and on the accompanying public GitHub repository.

Code availability

The custom code developed for this work is implemented in Python and is made available in a public GitHub repository under the Apache 2.0 licence: https://github.com/schwallergroup/alice

Acknowledgements

J.W.S., R.P.B., K.P. and R.B. thank Roche and its Technology Innovation and Science (TIS) initiative for financial support. We thank the Global Internship Program in Innovation & Sustainability (IP2TIS) 2025 (https://careers.roche.com/global/en/ip2tis-program) for financial support of this project. D.M.S., B.R., and P.S. acknowledge support from NCCR Catalysis (grant no. 225147), a National Centre of Competence in Research funded by the Swiss National Science Foundation. D.M.S. was funded by the Swiss National Science Foundation (SNSF) [226509].

Author contributions

J.W.S., B.R., R.B. and P.S. conceived the project. J.W.S., D.M.S. and B.R. developed the computational workflow with input from S.L.C., A.A., R.P.B. and P.S. J.W.S. and D.M.S. performed the computations and trained the models. J.W.S., M.J.N. and R.B. carried out the experimental work. J.W.S., K.P., M.J.N. and R.B. designed and planned the experimental case studies, with input from M.D.R.L. J.W.S., D.M.S. and B.R. analysed the results with help from M.D.R.L., K.P., M.J.N., R.B. and P.S. S.L.C., A.A., R.P.B., K.P., M.J.N., R.B. and P.S. provided supervision. J.W.S. wrote the paper with input from all authors.

Competing interests

J.W.S., D.M.S., M.D.R.L., A.A., R.P.B., K.P., M.J.N. and R.B. declare potential financial and non-financial conflicts of interest as full employees of F. Hoffmann-La Roche Ltd. The other authors declare no competing interests.

References

Supplemental main text benchmarking studies

Influence of different pre-trained models and pooling strategies

Design spaces for prospective optimisation campaigns

SI.pdf

†These authors contributed equally to this work.

*Corresponding authors. Email: [email protected]

High-throughput experimentation (HTE) platform

We conducted all high-throughput experiments using a parallel experimentation platform customdesigned by UnchainedLabs (Figure 1). This system comprises two interconnected Big Kahuna platforms, one of which is further integrated with a LiCONiC LiCotel system. The entire setup is encased within an LC Technology Solutions glove box featuring dual circulation systems, with one dedicated to solid dispensing and another to reaction execution. LC-MS analysis was performed using an ACQUITY UPLC I-Class system with QDa from Waters. For Chiral HPLC analyses, a Chiralcel OJ-3 column from Daicel was used, with ethanol and heptane as the mobile phase.

Solid components were dispensed with a Mettler balance in vial dispense mode, using SV hoppers for precursors and ligands and SV hoppers ($\le15$ mg) or 10 mL classic hoppers ($>15$ mg) for solid additives. For target dispense quantities $<0.4$ mg, materials were dispensed as coated ChemBeads. Following automated solid dispensing, liquid components (substrates, liquid additives, and solvents) were manually added using an Eppendorf Multipette E3 single channel pipette (4987000010). Reactions were performed in standard 96-position parallel synthesis reaction blocks (Analytical Sales and Services, SKU: 96960), with V&P Scientific super tumble stir discs (VP 721F-1) for stirring, or in a Screening Pressure Reactor (SPR) from UnchainedLabs. Data analysis was performed using the HTE OS workflow described by [1] HTE OS is integrated with a SpotFire application that enables tagging of LC-MS and HPLC signals into categories (e.g., “limiting SM”, “other SM”, “solvent”, “ignore peak”). During analysis, all peaks were tagged accordingly, with peaks corresponding to ligands and precatalyst components designated as “ignore peak”. For LCAP calculations, signals tagged as “other SM”, “solvent”, or “ignore peak” were excluded from consideration.

Prospective application: Palladium-catalysed cyanation

HTE experimental procedure

For all HTE plates, the following experimental procedure was performed. Palladium catalysts ($1$ mol%), cyanide sources, and solid additives were dispensed into 1 mL vials with stirring disks in a 96-well plate. Methyl 3-bromo-5-fluorobenzoate (19.0 μL, 129 μmol), solvents (240 μL), and cosolvent (120 μL) were added to the vials, followed by the addition of liquid additive NEt$_3$ (26.9 μL, 1.50 eq) where applicable. The plate was sealed and stirred (400 rpm) at the specified temperature for 20 h. 4/1 MeCN/H$_2$O (400 μL) was added to each well and shaken for 20 min at room temperature. Samples were taken and analysed by Liquid Chromatography-Mass Spectrometry (LC-MS).

Scale-up and synthesis of methyl 3-cyano-5-fluorobenzoate

A nitrogen-flushed 100 mL Easymax flask with overhead stirring was charged with methyl 3-bromo-5-fluorobenzoate (3.17 mL, 21.5 mmol) and DMC (40 mL). In a glovebox, a solution of K$_4$Fe(CN)$_6$ · 3 H$_2$O (4.53 g, 10.7 mmol, 0.50 eq) in water (20 mL) was prepared and added to the Easymax flask under argon. The mixture was heated to 75 °C. Once the temperature was reached, [Pd(allyl)(Xantphos)]Cl (163.4 mg, 215 μmol, 1.0 mol%) was added and the mixture was stirred at 75 °C for 4 h. The reaction was transferred with degassed 25% aq. NaCl solution (100 mL) to a separation funnel and extracted three times with degassed EtOAc (3 × 100 mL). The combined organic phases were washed with sat. aq. NaCl solution, dried over Na$_2$SO$_4$, filtered, and the solvent was removed using a rotary evaporator. The crude product was purified by flash column chromatography on silica gel (EtOAc:heptane = 0:100 to 30:70) over 20 min, 80 g, only UV active fractions collected), affording the product as a white solid. Yield: 3.61 g (94%).

$^{1}$H NMR (400 MHz, CDCl$_3$): δ 8.18 – 8.13 (m, 1H, ArH), 8.00 – 7.95 (m, 1H, ArH), 7.58 – 7.53 (m, 1H, ArH), 3.98 (s, 3H, OCH$_3$).

$^{13}$C{$^{1}$H} NMR (101 MHz, CDCl$_3$): δ 164.0 (d, $^{4}J_{\mathrm{F,C}}$ = 2.9 Hz, CO$_2$Me), 162.1 (d, $^{1}J_{\mathrm{F,C}}$ = 252.4 Hz, arom.), 133.8 (d, $^{3}J_{\mathrm{F,C}}$ = 7.7 Hz, arom.), 129.2 (d, $^{4}J_{\mathrm{F,C}}$ = 3.7 Hz, arom.), 123.1 (d, $^{2}J_{\mathrm{F,C}}$ = 24.6 Hz, arom.), 121.4 (d, $^{2}J_{\mathrm{F,C}}$ = 23.1 Hz, arom.), 116.7 (d, $^{4}J_{\mathrm{F,C}}$ = 2.9 Hz, CN), 114.3 (d, $^{3}J_{\mathrm{F,C}}$ = 8.8 Hz, arom.), 53.0 (OCH$_3$).

${}^{1}\mathrm{H}$ NMR (400 MHz, $\mathrm{CDCl}_{3}$): $\delta$ 8.18–8.13 (m, 1H, ArH), 8.00–7.95 (m, 1H, ArH), 7.58–7.53 (m, 1H, ArH), 3.98 (s, 3H, $\mathrm{OCH}_{3}$).

${}^{13}\mathrm{C}\{{}^{1}\mathrm{H}\}$ NMR (101 MHz, $\mathrm{CDCl}_{3}$): $\delta$ 164.0 (d, ${}^{4}J_{\mathrm{F,C}} = 2.9$ Hz, $\mathrm{CO}_{2}\mathrm{Me}$), 162.1 (d, ${}^{1}J_{\mathrm{F,C}} = 252.4$ Hz, arom.), 133.8 (d, ${}^{3}J_{\mathrm{F,C}} = 7.7$ Hz, arom.), 129.2 (d, ${}^{4}J_{\mathrm{F,C}} = 3.7$ Hz, arom.), 123.1 (d, ${}^{2}J_{\mathrm{F,C}} = 24.6$ Hz, arom.), 121.4 (d, ${}^{2}J_{\mathrm{F,C}} = 23.1$ Hz, arom.), 116.7 (d, ${}^{4}J_{\mathrm{F,C}} = 2.9$ Hz, CN), 114.3 (d, ${}^{3}J_{\mathrm{F,C}} = 8.8$ Hz, arom.), 53.0 ($\mathrm{OCH}_{3}$).

Prospective application: Asymmetric ketone hydrogenation

HTE experimental procedure

For all HTE plates, the following experimental procedure was performed inside a nitrogen filled glovebox. Iridium and ruthenium catalysts (1 mol%) and solid bases were dispensed into 1 mL vials in a 96-well plate. Stock solutions of 1,2-diphenylpropan-1-one (20 mg, 95.1 µmol in solvent (200 µL)) were prepared and added to the vials, followed by addition of liquid bases, where applicable. The plate was sealed inside the glovebox, placed in a screening pressure reactor (SPR) and shaken overnight at the indicated temperature under 20 bar H₂. The solvent was removed using a GeneVac and the residue redissolved in EtOH (200 µL). Samples were taken and analysed by chiral HPLC.

Scale-up and synthesis of (1S, 2R)-1,2-diphenylpropan-1-ol

In an argon filled glovebox, [Ir(Cl)H₂((S)-DTB-SpiroPAP-3-Me)] (93.1 mg, 95.1 µmol, 0.01 eq) and 1,2-diphenylpropan-1-one (2000 mg, 9.51 mmol, 1.0 eq) were added to a 50 mL autoclave. Degassed tAmOH (16.1 g, 20 mL, 19.2 eq) and DBU (724.0 mg, 711.2 µL, 4.76 mmol, 0.5 eq) were subsequently added to the reaction mixture. The autoclave was then sealed, pressurised with Ar (5 bar) and removed from the glove box. The autoclave was connected to a H₂ line, purged five times with H₂ (10 bar) and stirred at 60 °C and 500 rpm under 20 bar of H₂. Reaction progress was monitored via H₂ uptake, reaching completion within 4 h. The autoclave was allowed to cool to room temperature and carefully depressurised. The deep yellow reaction mixture was transferred with EtOAc (130 mL) to a separation funnel and washed sequentially with 1 M aqueous HCl solution (2 * 50 mL), water (50 mL) and saturated aq. NaCl solution (50 mL). The organic phase was dried over Na₂SO₄, filtered, and concentrated under reduced pressure to afford a lightly brown oil.

This crude oil was purified by flash column chromatography, eluting with EtOAc/heptane (18%, v/v). Solvent removal under reduced pressure yielded 1,2-diphenylpropan-1-ol as a transparent, oily solid (1.99 g). The crude solid was then dissolved in pentane (10 mL) and allowed to crystallise by slow evaporation over 48 h at room temperature. The resulting crystal bearing oily residue on the surface was transferred into a 15 mL glass vial. Ice-cold pentane (3 mL) was added to submerge the crystal. The vial was gently swirled for 3–5 seconds, after which the supernatant was drawn off. The crystal was rinsed a second time with ice-cold pentane (3 mL) and immediately decanted. The purified crystal was transferred onto lint-free filter paper and was allowed to dry to constant mass at room temperature, affording the product. Yield: 1687.89 mg (83.6%).

¹H NMR (400 MHz, DMSO-d₆): δ 7.22 – 7.07 (m, 10H, ArH), 5.31 (d, ³JH,H′ = 4.8 Hz, 1H, (H)COH), 4.64 (dd, ³JH,H′ = 4.6 Hz, ³JH,H′ = 6.4 Hz, 1H, (H)COH), 2.96 (p, ³JH,H′ = 6.9 Hz, 1H, (H)CCH₃), 1.24 (d, ³JH,H′ = 7.0 Hz, 3H, (H)CCH₃).

${}^{13}\mathrm{C}\{{}^{1}\mathrm{H}\}$ NMR (101 MHz, DMSO–$\mathrm{d}_{6}$): $\delta$ 145.41 (arom.), 145.18 (arom.), 128.52 (arom.), 128.23 (arom.), 127.95 (arom.), 126.99 (arom.), 126.90 (arom.), 126.25 (arom.), 77.49 (COH), 47.53 (CCH$_{3}$), 16.97 (CCH$_{3}$)

Enantiomeric data augmentation

All solvents and bases in the design space are achiral, and all 32 chiral catalysts are present as enantiomeric pairs. A catalyst and its enantiomer are therefore expected to yield mirror-image product distributions under otherwise identical conditions. Conversion and diastereomeric excess are unchanged, while the enantiomeric excess changes sign. We exploited this symmetry as a data augmentation strategy. For each observed reaction, a reflected counterpart was added to the training data in which the catalyst was replaced by its enantiomer, the sign of ee ($syn$) inverted, and conversion and de ($syn$) retained. Augmented observations were used to fit the surrogate models only. The acquisition function was evaluated over the full design space excluding conditions already run experimentally.

Model sensitivity to initial training data

The first initialisation plate in the HTE reaction optimisation campaign contained only five conditions producing the target $syn$ diastereomer. To assess how strongly the round-2 optimisation outcome depended on these observations, we repeated the optimisation with the top $n$ conditions removed from the Plate 1 training data, for $n = 1$, 3, and 5, ranking conditions by de ($syn$).

In each case the surrogate models were refitted on the reduced dataset and the acquisition function was evaluated over the full design space, excluding conditions already queried. No new experiments were performed: suggested conditions for which experimental data were already available were compared against the best condition remaining in the ablated training set. Because only suggestions already measured in round 2 could be evaluated, the values below are lower bounds on the performance recoverable at each ablation level.

With the single most diastereoselective condition removed ($n = 1$), the best remaining training observation reached 47.9% de ($syn$). The model recovered conditions delivering 99.4% conversion, 75.4% de ($syn$), and 89.5% ee ($syn$), improving stereoselectivity. With the top three removed ($n = 3$), the best remaining training observation reached 23.5% de ($syn$), and the model recovered conditions at 64.2% de ($syn$), again with higher enantioselectivity. Removing all five conditions ($n = 5$) that were $syn$-selective left no positive diastereoselectivity examples in the training data, yet the model still recovered conditions reaching 83.4% de ($syn$) and 93.3% ee ($syn$). Even with no $syn$-selective examples to learn from, the model identified conditions delivering high $syn$ selectivity.