Efficient discovery of visible light-activated azoarene photoswitches with long half-lives using active search
Photoswitches are molecules that undergo a reversible, structural isomerization after exposure to different wavelengths of light. The dynamic control offered by molecular photoswitches is favorable for applications in materials chemistry, photopharmacology, and catalysis. Ideal photoswitches absorb visible light and have long-lived metastable isomers. We used high throughput virtual screening to predict the absorption maxima (λ max ) of the E-isomer and half-lives (t 1/2 ) of the Z-isomer. However, computing the photophysical and kinetic properties of each entry of a virtual molecular library containing 10 3 -10 6 entries with density functional theory is prohibitively time-consuming. We applied active search, a machine learning technique to intelligently search a chemical search space of 255 991 photoswitches based on 29 known azoarenes and their derivatives. We iteratively trained the active search algorithm based on whether a candidate absorbed visible light (λ max > 450 nm). Active search was found to triple the discovery 1 rate compared to random search. Further, we projected 1 962 photoswitches to 2D using the Uniform Manifold Approximation and Projection (umap) algorithm and found that λ max depends on the core, which is tunable with substituents. We then incorporated a second stage of screening with to predict the stabilities of the Z-isomers for the top 1% of candidates. We identified four ideal photoswitches that concurrently satisfy λ max > 450 nm and t 1/2 > 2 hours; the range of λ max and t 1/2 range from 465 to 531 nm and hours to days, respectively.
Introduction<!>Methods<!>Results and discussion
<p>Light is an ideal external stimulus to promote organic reactions. Photoswitches are a class of molecules that absorb light and reversibly interconvert between their thermodynamically stable and meta-stable forms to create photostationary states. Azobenzenes are a class of well-studied photoswitches that undergo efficient isomerization from their thermodynamically stable form (i.e., E) to their metastable form (i.e., Z) using ultraviolet light (314 nm). 1 The Z → E isomerization is promoted with 365 nm light. 1 This relatively high-energy light (e.g., ultraviolet) may promote undesired side reactions that compete with the isomerization pathway (e.g., electrocyclic ring-closing reactions). UV light can also promote [2+2]-dimerizations that alter the structure and function of nucleotides and has a limited (epidermal depth, 0.1 mm) 2 tissue penetration depth, thus limiting the therapeutic potential of photoswitches in photopharmacology. The Z-isomer of azobenzene has a thermal half-life (t 1/2 ) of 4.7 hours, which prevents the establishment of photostationary states. Ideal photoswitches feature long absorption wavelengths and long t 1/2 ; unfortunately, the simultaneous optimization of these parameters is challenging and has been empirically observed to compete. Functionalizing the phenyl rings has been shown to shift the λ max of azobenzene-based photoswitches into the visible range. Konrad et al. 3 recently demonstrated that functionalizing the phenyl rings with halogens at the ortho positions led to a substantial red shift to 410 nm. This functionalization strategy also increased the (t 1/2 ) to 16 hours. Another strategy involves replacing one or both phenyl rings with heteroaryl ring(s), thus creating a more general class of photoswitches, azoarenes. Azoarenes are substantially more diverse than azobenzenes, and multiple examples show λ max in the visible range and t 1/2 exceeding 1.5 hours. Figure 1 highlights some of the most promising synthesized azoarenes with respect to λ max and t 1/2 . 3-9 While this relatively new class of azoarene photoswitches is attractive, the full enumeration of the chemical space approaches 10 6 . Density functional theory (dft) calculations are used to predict structures and photophysical properties at a relatively low computational cost. 10,11 Thus, dft has been previously used in high throughput virtual screening (htvs) [12][13][14][15] for virtual libraries containing 500-500 000 molecules. The vastness of the chemical space cannot be understated; conservative estimates suggest that 10 23 organic molecules are theoretically possible. 16 This figure can be narrowed to roughly 10 6 for azoarenes by focusing on those already experimentally realized. Abreha et al. 17  Even with high-performance computing and efficient quantum chemistry codes, computing the photophysical properties and stabilities of 10 5 photoswitches is a substantial undertaking. We have employed the machine learning algorithm 'active search' 19 to intelligently search the vast chemical space (255 991 candidates) of azoarene photoswitches. Active search (as) was created to discover as many target molecules as possible while balancing computational resources. as uses the data observed at any given point throughout a search and adaptively makes decisions informed by the latest observations. The prediction accuracy of our predictive model improves as we frequently query from quantum chemical calculations.</p><p>We now combine these existing tools (the verde materials db, 17 EZ-TS, 18 and active search 19 ) to automatically identify top photoswitch candidates featuring visible-light λ max and long t 1/2 . Scheme 1 shows an illustration of the iterative processes used to identify ideal photoswitches Scheme 1: The multipronged iterative procedure used to update the active search algorithm with dft results.</p><p>Phase 1 : An initial screen of 50-100 molecules is processed through an automated computational workflow developed by Abreha et al. 17 . RDKit 20 is used to generate 3-D coordinates from a simplified molecular-input line-entry system (smiles) 21 string, followed by a low-mode conformational search where each conformer (4 total) is minimized with the Universal Force Field. 22 The lowest energy conformer is determined through semi-empirical optimizations and a single-point energy calculation. The lowest energy structure is optimized with M06 23 /6-31+G(d,p) 24,25 and IEFPCM MeCN , 26 and a vibrational analysis confirms the stationary point as the true minimum if it has only positive frequencies. The λ max is calculated with a single point energy calculation using ωB97XD 27 /6-31+G(d,p)//M06 23 /6-31+G(d,p).</p><p>Figure 2 shows the automated workflow of quantum chemical calculations used to compute the excitation energies and corresponding λ max for selected molecules from our virtual library.</p><p>Figure 2: Quantum chemical workflow for computing the λ max for all molecules considered in this study.</p><p>Phase 2 : An in-house Python script assigns a "core id" (1-29) to each computed structure.</p><p>Cores are determined using a substructure analysis included in RDKit. True or False labels are assigned to each smiles string based on the pre-determined threshold, λ max > 450nm.</p><p>Phase 3 : A machine learning model is trained on the set of molecules that are labeled to guide the search algorithm. First, we generate the Morgan fingerprint 28 of each molecule and compute the Tanimoto similarity coefficient 29 between each pair of molecules. We then build a k-nearest neighbors (k-nn) 30 predictive model that computes the probability of a given unlabeled molecule having a positive label, given the data we have observed thus far. This k-nn model is then utilized by the search algorithm. Note that the Morgan fingerprints and Tanimoto similarity coefficients only need to be computed once, while the k-nn is updated with newly labeled data at each iteration Phase 4 : The active search algorithm builds the set of 50 recommendations, selecting among all unlabeled molecules. These recommendations are then sent to Phase 1 to be computed and labeled. This procedure repeats for a total of 40 iterations, sampling 1 962 molecules from the space. We include a more detailed description of our methods in the following section.</p><!><p>We adapted the active search method, which has shown impressive performance in molecular discovery in previous studies. [31][32][33][34][35] The method was first introduced by Garnett et al. 19 and extended to the batch setting by Jiang et al. 35 . Formally, suppose we have a large set of elements X = {x i }, among which there is a small subset R ⊂ X of valuable elements that we wish to search for (i.e., molecules exhibiting a desired property). We do not know which members of X belong to R a priori, but whether a specific element x belongs to R can be determined by querying an oracle, requesting for the binary label y = 1{x ∈ R}, where 1{•} is the indicator function. In this work, the binary label denotes whether a molecule exceeds the λ max threshold of 450nm. Further, we assume that at each iteration of the search, b elements are inspected simultaneously, requiring that queries to the oracle be made in batches of size b. This models experimental settings in which multiple experiments may be run in parallel to maximize throughput, contrasting with the fully sequential setting where queries are made one after another; here, b = 50. The goal is to design a sequence of queries limited by a predetermined budget, such that the number of target elements uncovered by querying the oracle is maximized. As such, we naturally define the utility of a given set of observations D = {(x i , y i )} to be the total number of targets found:</p><p>We aim to determine the sequence of queries that maximizes our definition of utility in the expected case using Bayesian decision theory. This framework first requires a classification model that computes the posterior probability that an unlabeled point x belongs to R, given the elements we have inspected thus far in D, Pr (y = 1 | x, D). The active search method is model-agnostic and does not make any further assumptions about this predictive model.</p><p>In the next section, we describe the k-nearest neighbors model we use for this classification task.</p><p>We denote T = t b to be the total number of queries allowed to be made given our budget, where t is the number of search iterations). We further denote by D i the observations collected at the end of iteration i. At iteration i + 1 ≤ t, the best batch of queries (of size b) we can make, denoted as X i+1 , maximizes the expected value of the utility of the dataset at termination D t :</p><p>Although this expected utility can be derived using the standard procedure of backward induction, 36 it involves t−i nested steps of sampling over unknown labels of candidate queries and maximizing the future expected utility. This computation is prohibitively expensive for horizons t − i ≥ 3, rendering the optimal query infeasible to calculate in practice.</p><p>We adopt the sequential simulation strategy proposed by Jiang et al. 35 as an efficient approximation to the optimal batch of queries. First, the strategy builds on the efficient nonmyopic search algorithm ens 34 in the sequential setting where only one query is made at each iteration. ens itself approximates the optimal sequential strategy by assuming that all future queries after the current iteration are made at the same time. Jiang et al. 34 demonstrated that ens actively explores the search space when the remaining budget is large, recommends increasingly promising molecules as the search progresses, and achieves significant improvements in performance over greedy strategies Our sequential simulation active search algorithm under the batch setting builds its recommendations by iteratively adding elements to an initially empty set using the ens algorithm until the desired size (b = 50) is reached. As a new element is added, we assume that this element will return a negative label (i.e., the element is assumed to lack the desired property). Jiang et al. 35 showed that by taking on this pessimistic view, the algorithm encourages the elements within the same batch to be diverse, which helps explore the search space more effectively.</p><p>Finally, we aim to distribute our queries equally across the 29 cores. Our sequential simulation strategy may be naturally modified in service of this goal as follows. As a new element is added to the running batch in the iterative procedure described above, we temporarily remove other candidates having the same core id as the newest batch member from the search space. When no candidate remains, we add all removed molecules back to our search space. This simple procedure effectively forces each batch of queries to be constructed to span the available cores equally.</p><p>As previously described, our active search algorithm requires a probabilistic model that computes the probability that an unlabeled element has a positive label (i.e., exhibiting the desired property), given the current set of observations we have made so far. We first generate the Morgan fingerprint 28 of each molecule in our search space and compute the Tanimoto similarity coefficient 29 between each pair of elements x and x , denoted as t(x, x ).</p><p>We then implement a k-nearest neighbor (k-nn) 30 predictive model, which computes the probability of an uninspected molecule being an active compound as:</p><p>where nn(x) is the labeled subset of the k nearest neighbors of x in X. γ is a parameter of the model that acts as a "pseudo count" to define the prior probabilities for molecules that do not have any labeled neighbor; we set γ = 0.1. This k-nn performs well in previous work, 19,31,32,34,35 as well as in our experiment. It can further be rapidly updated in light of new observations, allowing for efficient lookahead computations that are central in active search.</p><!><p>We generated a relatively small virtual molecular library of 1 636 azobenzene, bisazopyrrole, bisazothiophene, and bisazofuran photoswitches (Figure 3). The substituent sites (red circles) were replaced with the disubstituted alkenyl, alkynyl, or aryl (spacer) groups. The unfunctionalized end of the growing molecule (-R) was substituted with functional (terminal) groups. Figure 3 shows the sites where a set of 4 azoarene cores were substituted with spacer and terminal groups to generate the initial training set. From the 1 636 initial azoarenes, 198 were selected to train the active search algorithm.</p><p>A histogram of the λ max of these 198 azoarenes is shown in Figure 4.</p><p>Figure 4 shows that the λ max ranges from 301 to 541 nm for the selected 198 azoarenes.</p><p>To train the as algorithm, we assigned each candidate a label of True or False, depending on whether the following expression is satisfied, λ max > 450 nm. 62 of the 198 azoarenes were assigned True and 136 were assigned False. We designed a virtual molecular library with 29 bisdiazoarene cores (Figure 5) to apply the trained algorithm. Each of these has at least one functionalization site substituted with functional groups (i.e., terminals).</p><p>The cores were selected based on a literature search of previously synthesized azoarenes.</p><p>1-29 range from symmetric bisazoarenes to azoheteroarenes and known functionalization strategies inspire the substitution sites. Figure 6 describes these positions for a smaller subset of cores.</p><p>We then iteratively applied the algorithm 40 times on our new molecular dataset. Each molecular batch featured 50 as-suggested candidates that would enter our computational workflow. The first 20 iterations used an "equidistributed" policy, which equally sampled molecules belonging to each core family of the 29. Since the as selected 50 molecules for each iteration, we sampled the 29 cores by constraining the algorithm to select at least one molecule per core. The remaining 21 slots for each batch were selected in a similar fashion where no more than two molecules were selected for each core. The remaining iterations (21-40) used a "targeted" policy that only selected molecules from a subset of 15 cores that had derivatives where the λ max > 450 nm. Cores that did not show derivatives that fit the criteria were excluded from the subset. After each iteration, we added a binary label to each  molecule based on whether λ max > 450 nm. Figure 3 summarizes this iterative procedure.</p><p>We compared the as strategy to the performance of a random search strategy by sampling three molecules selected at random from each of the 29 cores. Figure 7 shows the distribution of the λ max values from as and the random search. The random search showed that 11 out of the 87 molecules (13%) had λ max > 450 nm.</p><p>Figure 8 shows how the proportion of hits changes with respect to the first 20 iterations using the equidistributed policy. We define the hit rate as the percentage of molecules with a λ max > 450 nm from the current batch.</p><p>The dotted orange line indicates a random search hit rate of 13%. The black data points indicate the hit rate as the active search is iteratively applied. The equidistributed search shows a range of hit rates from [12% to 35% (batch 3 and 18, respectively)]. The slope is +0.82; the hit rate is improved relative to the random search in nearly all iterations. We then turned our attention to the targeted as policy to maximize the number of hits corresponding to the subset of cores with molecules that had a λ max > 450 nm, shown in Figure 9.  For iterations 21-40, the as algorithm selected three derivatives corresponding to each of the 15 cores for a total of 45 selected molecules. To keep the batch size consistent to 50, as chooses five more from the top-ranked derivatives of the 15 core subset. Figure 10 shows the hit rate for iterations 21-40 with the targeted policy. In the targeted policy, the hit rate varied from 44% to 56%; the average hit rate was 49%. Unlike the equidistributed policy, Figure 10 does not show an increase in hit rate as a function of the batch number. The relatively high hit rate led to the rapid discovery of 485 candidates with λ max > 450 nm in batches 21-40.</p><p>Overall, we identified a total of 717 photoswitches with λ max > 450 nm after the 40 batches (1 962 molecules) of as-assisted virtual screening. The resulting hit rate is 37%, corresponding to a tripling of the 13% hit rate from the random search. A two-sample z-test rejects the null hypothesis that the two strategies result in equal hit rates with overwhelming confidence, yielding a p-value of 5 × 10 −6 .</p><p>We represented the complex molecular data with a Uniform Manifold Approximation (umap) 37 to visualize the molecular motifs responsible for candidates with λ max > 450 nm.</p><p>Each of the 1 962 structures was plotted based on the Tanimoto similarity 29 in Figure 11.</p><p>The clusters are grouped based on structural similarity and color-coded based on computed λ max results.</p><p>Figure 11a shows the umap results with each azoarene candidate overlaid with the color corresponding to the λ max . The data points shown in grey correspond to the ultraviolet range of the electromagnetic spectrum (λ max < 400 nm). Cores 1-5, 17, 24, and 25 formed distinct clusters, indicated by the dotted lines in the umap plot. These cores also had considerably more derivatives with a λ max in the visible range, suggesting that these cores have especially tunable λ max values and should be explored experimentally in the future.</p><p>We examined the influence of substituents on each core by plotting the distribution of λ max . Figure 11b shows the range of λ max for 1 962 azoarenes. Spacings within each box represent the degree of dispersion and skewness within the data. Cores with larger boxes indicate a higher variation in absorbance due to the substitution pattern. We compared unsubstituted cores 1-5, 17, 24, and 25 to the derivative with the highest λ max . These values are summarized in Table S2 of the supporting information.  with the longest λ max , identified with active search. Figure 12 illustrates the candidate from each family of cores subjected to transition state calculations with PBE036-D3/6-31+G(d,p) to optimize the transition states. This was reported to give activation free energies that approach chemical accuracy. Scheme 2 shows the Z → E isomerization transition state. The λ max for these top 29 candidates ranges from 382 to 602 nm. The range of activation free energies is 8.1 to 30.0 kcal mol −1 . We plotted the activation free energies (∆G ‡ ) against Figure 13 shows no linear relationship between the λ max and activation free energy (R 2 of 0.0002). However, we divided the plot into four quadrants to highlight those candidates that meet both, one, or none of the λ max and t 1/2 optimization criteria. Quadrants A (red tint) and B (green tint) contain molecules that have λ max > 450 nm or 2.6 eV. Quadrants</p>