Conceptio › Archive › arXiv CS

arXiv CSopen access

Computational Design and Experimental Validation of Photoactive PARP1 Inhibitors

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

Computational Design and Experimental Validation of Photoactive PARP1 Inhibitors

arXiv:2604.24634v1 [physics.chem-ph] 27 Apr 2026

Simon Axelrod,∗,†,‡,§ Miroslav Kašpar,∗,¶,§ Kristýna Jelı́nková,¶ Markéta Šmı́dková,¶ Erika Bartůňková,¶ Sille Štěpánová,¶ Eugene Shakhnovich,† Václav Kašička,¶ Martin Dračı́nský,¶ Zlatko Janeba,¶ and Rafael Gómez-Bombarelli‡ †Department of Chemistry and Chemical Biology, Harvard University, Cambridge, MA, 02138 ‡Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA, 02139 ¶Institute of Organic Chemistry and Biochemistry of the Czech Academy of Sciences, Flemingovo nám. 2, 160 00 Prague 6, Czech Republic §These authors contributed equally to this work. E-mail: [email protected]; [email protected]

Abstract

Introduction

Light-activated drugs are a promising way to treat localized diseases for which existing treatments have severe side effects. However, their development is complicated by the set of photophysical and biological properties that must be simultaneously optimized. Here we used computational techniques to find a set of promising candidates for the photoactive inhibition of the poly(ADPribose) polymerase 1 (PARP1) cancer target. Using our recently developed methods based on atomistic simulation and machine learning (ML), we screened a set of 5 million hypothetical photoactive ligands. Our workflow used protein-ligand docking to identify candidates with differential PARP1 binding under light and dark conditions; ML force fields and quantum chemistry calculations to predict pKa , absorption spectra, and thermal half-lives; graph-based surrogate models to screen additional compounds; excited-state nonadiabatic dynamics with ML force fields to estimate quantum yields; and free energy perturbation (FEP) to refine binding predictions. From these predictions, we prioritized a small set of synthetically feasible candidates expected to have red-shifted absorption spectra, thermal halflives on the order of seconds to minutes, and isomerdependent PARP1 binding under visible-light control. We synthesized 10 candidates and experimentally characterized their photobehavior and PARP1 inhibition constants. Among the validated compounds, 1 showed a 15-fold increase in inhibition of PARP1 upon greenlight irradiation at 519 nm (208.8 ± 28.3 µM vs 14.4 ± 1.9 µM). These results validate the computation-guided screening strategy for identifying red-shifted PARP1 photoinhibitors, while also underscoring current limitations such as rapid thermal relaxation in aqueous media.

Photopharmacology is an emerging field based on the optical control of drug activity 1,2 . By activating a drug at a specific location or time, one can minimize off-target effects and increase the maximum deliverable dose, thereby improving quality of life 1 . These therapeutics have been proposed as a treatment modality for cancer 3 and blindness 4 , among other conditions 1 . They are typically built around photoswitches such as azobenzene (Fig. 1(a)), which undergoes trans ↔ cis isomerism in response to light. Because bioactivity depends strongly on ligand shape, the large structural change induced by azobenzene photoisomerization can alter the drug’s biological effect. Despite its promise, photopharmacology is hindered by a difficult multi-objective optimization problem. Photoactive drugs should ideally absorb light in the near-infrared (IR) range between 700 and 900 nm 5 , since only light in this range can penetrate human tissue. Yet unsubstituted trans azobenzene absorbs in the UV and visible range 6 . The thermal lifetime of the active drug, typically the cis isomer generated by excitation of the trans form, should be long enough to allow accumulation and the onset of a biological effect. Yet it should be short enough to prevent transport through the bloodstream and hence biological activity elsewhere 5 . Lastly, the active isomer should exhibit strong biological activity, whereas the inactive isomer should exhibit weak or no activity. Further, the numerous properties that must be optimized for standard pharmaceuticals also need to be optimized for photoactive ones. These include solubility, permeability, metabolic stability, and toxicity. Balancing these objectives leads to constraints on orally absorbed drugs, most famously quantified by Lipinski’s rule of five 7 . This further complicates the design 1

b O

5 million pairs

Cl Br F

Docking O O NH S NH2 O N O H

N O

Binding free energy

combinatorial substitution

PAINS

MCF

Inspection

SC score

5 million pairs

Neural nonadiabatic MD

Absolute FEP

random sampling 4,000

explore / exploit x2

Quantum chemistry

Neural force field

Synthetic feasibility SA score

Binding free energy

30,000

Drug likeness QED

Quantum yield

48 pKa

Absorption spectrum

Thermal halflife

2D model

Figure 1: Virtual screening workflow. Molecular libraries are shown in yellow ovals, calculation types in blue rectangles, and computed properties in green rectangles. (a) Molecules were generated through combinatorial substitution and filtered based on drug likeness and synthetic feasibility. QED: quantitative estimate of drug-likeness 8 ; PAINS: pan-assay interfering compounds 9,10 ; MCFs: medicinal chemistry filters 11 ; SA score: synthetic accessibility score 12 ; SC score: synthetic complexity score 13 . (b) 5 million cis-trans pairs (10 million molecules) were docked into PARP1, and 30,000 were selected based on these docking scores. Next, 4,000 pairs were randomly chosen to receive thermal half-life and absorption spectrum calculations. These calculations involved atomistic simulations using our neural network force field (NFF). Some also received pKa and azonium half-life calculations. These results were used to train a fast 2D model, which was in turn used to select 4,000 new molecules based on favorable properties (exploitation) or uncertainty (exploration). (c) Compounds with the best chemical and photophysical properties were selected for further validation. The predicted pKa values, absorption spectra, and thermal half-lives were validated with quantum chemistry. The remaining compounds were tested for strong differential binding through absolute free energy perturbation (FEP) calculations, and non-zero quantum yields through NFF simulations of excited-state molecular dynamics (MD).

of light-activated drugs. For example, constraints on molecular mass and the number of aromatic rings 8 limit simple strategies for red-shifting azobenzene absorption. In this work, we tackle this multi-objective optimization problem in photopharmacology by focusing on photoswitchable inhibition of poly(ADP-ribose) polymerase 1 (PARP1), a cancer-relevant target that we previously predicted to be photodruggable (i.e., a protein for which cis azobenzene derivatives bind more strongly than trans) 14 . Cancer therapy is one of the most promising applications for photopharmacology, because the side effects of many existing treatments are severe 1 . To identify optimized photoactive inhibitors of PARP1, we applied computational methods of increasing accuracy and decreasing speed to a large virtual

library of azobenzene derivatives. Our methods combined computational docking, machine learning (ML), quantum chemistry, and atomistic simulation to predict properties of millions of azobenzene derivatives. We identified promising compounds with red-shifted absorption spectra, thermal half-lives in the range of seconds to minutes, and strong differential binding affinities. We then synthesized a small set of compounds and measured these properties experimentally.

Results and discussion

the most promising compounds for further calculation. Compounds with thermal half-lives between seconds and minutes, and either red-shifted trans absorption spectra or high trans pKa values, were further validated with quantum chemistry. The remaining species received the most intensive calculations. We used an NFF to perform excited-state, non-adiabatic molecular dynamics (NAMD) to check that each species had a non-zero photoisomerization quantum yield 16 . We also performed absolute free energy perturbation (FEP) calculations to compute the binding affinity of each compound to PARP1. These were the most computationally demanding calculations, taking five days per compound with four graphics processing units (GPUs).

Virtual screening To generate a large virtual library for screening, we substituted azobenzene with common literature functional groups and patterns (Fig. 1(a) and Computational Supplementary Information (CSI) Figs. S2-S4). We then filtered these molecules by drug likeness and synthetic feasibility. Filters consisted of the quantitative estimate of drug likeness (QED) 8 , medicinal chemistry filters (MCFs) 11 , removal of non-drug-like groups through inspection, and removal of pan-assay interfering compounds (PAINS) 9,10 . Synthetic feasibility was measured through the synthetic accessibility and complexity scores (SA and SC scores, respectively) 12,13 . This resulted in a library of 5 million cis-trans pairs. Further details of this process are given in CSI Sec. 2. We then filtered the library based on computed chemical, photophysical, and biological properties. We used the approach of a computational funnel, reducing the library size with calculations of increasing accuracy and decreasing speed (Fig. 1(b)). We began with protein-ligand docking to predict the binding affinity of each isomer for each compound. This was the fastest but least accurate type of calculation. We used AutoDock Vina 15 to dock all 10 million molecules into PARP1. We kept cis-trans pairs for which the cis isomer was in the best 15% of scores, and the trans isomer was in the worst 50%. This left 30,000 pairs. Next we computed chemical and photophysical properties for 12,000 cis-trans pairs. We did so using a set of ML and simulation tools that we recently developed 14,16,17 . We used a neural network force field (NFF) trained on quantum chemical data to compute cis thermal half-lives and trans absorption spectra 17 . NFFs are orders of magnitude faster than quantum chemistry, and hence can be used to screen far more candidates. We were also interested in the spontaneous protonation of azobenzene to form azonium in water. Protonation of one of the azo nitrogens leads to red-shifted absorption spectra and shorter half-lives 5,18 . Appropriate substitution of azobenzene can lead to trans azonium formation at physiological pH, and hence to the desired red-shift 5,18 . We therefore computed the azonium pKa for compounds with ortho hydrogen bond acceptors, since these can stabilize the hydrogen atom and hence elevate the pKa 5,18 . This was combined with thermal half-life and absorption spectrum calculations for the azonium species. These calculations involved molecular dynamics simulations and millions of geometry optimizations. This made them rather costly, even with an NFF. We therefore could not use them to screen all 30,000 isomer pairs. Instead we screened a portion of the library, and used the resulting data to train models that predicted these properties from the 2D molecular graph. We used these fast 2D models to screen the entire library and identify

Hits Virtual screening results are shown in Fig. 2. Panels (a) to (c) show the results of the computational funnel, starting with docking and followed by NFF calculations of absorption spectra, half-lives, and pKa values. As we found in previous work 14 , the docking scores of the cis and trans isomers are in fact highly correlated (Pearson R = 0.72, Spearman ρ = 0.69). Nevertheless, we identified 30,000 compounds with fairly strong cis binding affinities and fairly weak trans affinities (bottom 15% of scores and top 50% of scores, respectively). From these compounds we selected species with half-lives between 5 seconds and 5 hours, and absorption wavelengths above 600 nm (Fig. 2(b)). Absorption wavelengths were computed as λ5% max , defined as the maximum absorption wavelength at 5% of the maximum intensity. We also selected compounds with trans pKa values above 5 and half-lives in the desired range (Fig. 2(c)). As with docking, we see that compounds with the desired properties in the right two panels made up only a small portion of the overall library. This reinforces the difficulty of the optimization problem and the utility of large-scale screening. The top three candidates are shown in Fig. 2(d). Visualizations and computed properties are shown in Fig. 3, and further data can be found in CSI Sec. 1. The first two compounds were predicted to be moderate differential binders of PARP1, with affinities that are typical for early-stage hit finding (Fig. 3, columns (i) and (ii)). Compound 1 was predicted to form an azonium ion at physiological pH, with a predicted trans pKa of 6.8. Based on experimental results in the literature, this suggests a red-shifted absorption spectrum, though the azonium spectrum generated by TDDFT was not in fact red-shifted (Fig. 3, row (a), column (iii)). This is discussed further below. Compound 2 was predicted to have a red-shifted absorption spectrum with no azonium formation (Fig. 3, row (b), columns (iii) and (vi)). Both compounds 1 and 2 were predicted to have halflives of 1.0 min in water at 37◦ C (Fig. 3, column (iv)). Each compound had a predicted quantum yield of 0, and hence no photoswitching. Compound 3 was predicted to be a weaker binder

Strong differential binders

Redshifted, good half-lives

High pKa, good half-lives

Figure 2: Results of virtual screening. (a) We first identified strong differential binders through docking. We then filtered remaining molecules based on (b) absorption spectra and thermal half-lives (t1/2 ), or (c) high pKa values and thermal half-lives. (d) Final hits after validation with quantum chemistry and FEP. i

iii

a QY = 0%

𝐾! = 466 𝜇𝑀

𝑡!/# = 1.0 min

𝐾! = 9,800 𝜇𝑀

pKa = 6.8

QY = 0%

𝐾! = 591 𝜇𝑀

𝐾! = 4,000 𝜇𝑀

𝑡!/# = 1.0 min

pKa = 0.9

QY = 4%

𝐾! = 2,100 𝜇𝑀

𝑡!/# = 101 min

𝐾! = 6.7 ×10" 𝜇𝑀

pKa = −0.5

Figure 3: Properties, plots, and graphics for the top candidates. (a) Compound 1. (i) Cis binding pose from the last frame of FEP. (ii) As in (i), but for trans. (iii) Absorption spectrum of trans azobenzene, trans 1, and trans 1-azonium (denoted by 1[H+]). (iv) Transition state (TS) for thermal isomerization of 1. Note that the lowest free energy TS involves azonium. (v) Kernel density estimation plot of the CNNC dihedrals and NNC angles at which hopping occurs in NAMD. NAMD was performed with 1-azonium. (vi) Azonium isomer with the lowest free energy. (b) As in (a), but for compound 2. Unlike in (a), the TS does not involve azonium, and NAMD was performed with the azobenzene form. (c) As in (b), but for compound 3.

than the other two compounds, but with a more redshifted spectrum, a non-zero quantum yield, and a longer (but perhaps too long) thermal half-life (Fig. 3, row (c)). Although compound 3 exhibited a favorable predicted photophysical profile, including the only nonzero predicted quantum yield, we did not pursue its syn-

thesis or the synthesis of its analogues due to practical resource constraints. Our synthetic efforts were instead directed toward derivatives of 1 and 2, for which building blocks were commercially available, azo coupling was expected to be robust, and handling was straightforward. As part of this focus, we removed the bulky

Absorbance [AU]

2.0

Scheme 1: General synthetic strategy for the azobenzene library. Ar1 denotes the electron-poor aryl partner, and Ar2 denotes the electron-rich aryl partner.

PSS

1.5

1.0

0.5

tert-butyl group from 2 to improve synthetic accessibility, and to increase the probability of a nonzero quantum yield. The latter was motivated by the predicted zero yield, which we attributed to the bulky ortho substituent that can often inhibit photoswitching 16 . We also hypothesized, and subsequently confirmed, that the azonium form of 1 would in fact be red-shifted, despite TDDFT predictions to the contrary.

0.0 250

300

350

400

450

500

550

600

650

Wavelength [nm]

Figure 4: Absorbance spectrum change of compound 2a upon irradiation by 519 nm LED in acetonitrile at 22◦ C, at concentration 4.452e-5 M.

Synthesis The photoswitches were prepared by late-stage assembly of the azobenzene linkage. The key disconnection forms the Ar–N=N–Ar motif by diazotization of an aniline precursor, followed by azo coupling with an electron-rich arene/heteroarene nucleophile (Scheme 1). Diazotization was typically conducted with sodium nitrite at 0–4◦ C using either p-toluenesulfonic acid or tetrafluoroboric acid as a proton source. Overall, we attempted to synthesize 20 compounds and successfully created 10, with isolated yields ranging from 9–59%. Experimental procedures, characterization, and unsuccessful attempts are reported in the Experimental Supplementary Information (ESI).

Absorption spectra. To compare UV-Vis absorption profiles, we report λ5% values, defined as the wavelength at 5% of maximum intensity relative to the maximum absorbance of the most red-shifted peak. Our computational pipeline correctly prioritized compounds with red-shifted absorption spectra: All compounds exhibited red-shifted absorption relative to unsubstituted azobenzene, with λ5% values in acetonitrile ranging from 557 nm (2h) to 616 nm (2b). This red-shift arises from the push–pull electronic configuration imparted by electron-donating piperazine or pyrrolidine groups and electron-withdrawing nitrile or amide substituents. 19 The effect of the electron-withdrawing group is clearly illustrated by the ortho-dihalogenated series: decreasing halogen electronegativity (F > Cl > Br) correlates with decreasing λ5% values of 580 nm, 562 nm, and 557 nm, respectively, for compounds 2f, 2g, and 2h. Compounds generally showed modest additional red-shifts of 2–8 nm when measured in aqueous buffer compared to acetonitrile. For example, compound 2a shifted from λ5% = 576 nm in acetonitrile to 578 nm in PBS/Tween 80. Compounds 1 and 1a, which possess measurable pKa values, exhibited pH-dependent absorption spectra. As shown for compound 1 in Fig. 5, the addition of trifluoroacetic acid (5 ppm) red-shifted λ5% to 606 nm, consistent with partial azonium formation, while triethylamine (5 ppm) blue-shifted the spectrum to λ5% = 556 nm. Compound 1a displayed analogous behavior (acid: 627 nm; base: 567 nm; see ESI Figures S31–S38). The azonium red-shift agrees with our expectations, but not with TDDFT predictions (Fig. 3, row (a), column (iii)). As shown in CSI Table S1, we observed a surprising trend for asymmetrically substituted azobenzenes: the azonium isomer with the higher pKa , and thus the greater likelihood of forming, often had a blue-shifted TDDFT spectrum. The isomer with the

Photophysical Characterization Absorption spectra, thermal relaxation kinetics, quantum yields, and photostationary state (PSS) compositions were measured in acetonitrile and in phosphatebuffered saline (PBS) containing 0.01% Tween 80 to approximate assay-relevant conditions and prevent sample precipitation. Compounds, experimental results, and computational predictions can be found in Table 1. Protonation state. To establish compound protonation states during switching experiments, we determined thermodynamic pKa values by capillary electrophoresis. Most compounds showed no evidence of protonation at pH 2, indicating pKa < 0. Three exceptions were compounds 1 (pKa = 3.78), 1a (4.29), and 2b (4.15). Of these, compound 1 was specifically chosen because of its predicted near-neutral pKa of 6.8. The experimental value was significantly lower, indicating that our computational approach predicted an overstabilization of the azonium form. The experimental results indicate negligible azonium formation at physiological pH for all tested compounds, meaning the observed photoswitching behavior arises from the neutral azobenzene chromophore rather than the protonated azonium species. 5

140

protonated neutral

120

Response (% of Control)

Absorbance [AU]

0.8

0.6

0.4

0.2

100 80 60 40 20 0

0.0 200

300

400

500

600

700

1, PSS, 519 nm, N=2, n=3 1, dark, N=2, n=3 2a, PSS, 519 nm, N=2, n=3 2a, dark, N=2, n=3

800

Concentration (μM)

Wavelength [nm]

Figure 6: Dose–response curves for PARP1 inhibition by compounds 1 and 2a in the dark (predominantly trans isomer) and under continuous 519 nm irradiation (PSS). Compound 1 exhibits an approximately 15fold potency enhancement upon illumination (IC50 = 208.8 ± 28.3 µM in the dark vs. 14.4 ± 1.9 µM in the light), while compound 2a shows an approximately 5-fold enhancement (209.4 ± 51.9 µM in the dark vs. 39.7 ± 4.7 µM in the light). Data points represent mean ± SD of technical triplicates (n = 3) from two independent biological replicates (N = 2).

Figure 5: Titration of compound 1 with trifluoroacetic acid in acetonitrile. Each step is the addition of 10 µL of 0.01% TFA in ACN into 2.0 mL of azobenzene dissolved in ACN at starting concentration 2.609e-5 M.

lower pKa , and hence the lower likelihood of forming, usually had a red-shifted spectrum. We hypothesized that the predicted blue-shift was a TDDFT artifact, and therefore prioritized compounds with high pKa values regardless of their predicted spectra. The experimental azonium spectra for compounds 1 and 1a support this hypothesis. We had targeted compounds 1 and 1a because of their high predicted pKa values, which meant that they would undergo a red-shift upon azonium formation under physiological conditions. Although the azonium forms were indeed experimentally red-shifted, the measured pKa values were too low for this effect to be relevant under physiological conditions. This reflects a shortcoming of our computational predictions. Moreover, even if the pKa values had been high enough, the spectral red-shifts would have been insufficient to reach the near-IR absorption of 720 nm achieved by the best azonium compounds in Ref. 5 . So while our computational approach produced red-shifted azobenzene derivatives, it did not produce compounds with the desired near-IR absorption achieved in Ref. 5 . Thermal relaxation. Thermal half-lives (t1/2 ) in acetonitrile at 22◦ C ranged from < 0.01 min (2b) to 5.95 ± 0.89 min (2a), placing most compounds within the computationally predicted range of seconds to minutes. Increasing the temperature from 22◦ C to 37◦ C reduced half-lives by approximately 40–50%, as expected from Eyring theory 20 (e.g., compound 2a: t1/2 = 3.30 ± 0.19 min at 37◦ C). A critical finding was that thermal relaxation in aqueous buffer was dramatically accelerated: all compounds exhibited t1/2 < 0.01 min (<0.6 seconds) in PBS/Tween 80 at 22◦ C, faster than our instrumental time resolution (ESI Figures S39–S48). This aqueous acceleration, likely driven by hydrogen bonding and

solvent polarity effects on the rotational isomerization barrier 21–23 , was not captured by our implicit-solvent computational model and represents a significant limitation for in vivo applications. However, Herges and coworkers recently demonstrated that an azo photoswitch with millisecond aqueous half-life could nonetheless exhibit prolonged biological activity (>3 h) when bound to its receptor, presumably because the protein environment stabilizes the active isomer 24 . This precedent suggests that fast aqueous relaxation may not preclude functional photopharmacology if target binding is sufficiently strong. Quantum yields. Photoisomerization quantum yields (Φtrans→cis ) were determined at 519 nm in acetonitrile using indolyl-fulgide actinometry 25 (ESI Figures S66–S73). An example of the spectral evolution upon irradiation is shown in Fig. 4. Measured quantum yields ranged from 1.2 ± 0.2% (1) to 5.2 ± 1.7% (1a). However, no photoswitching was observed for compounds 1 or 1a under acidic conditions, where both exist in the trans azonium form (ACN + 5 ppm TFA, 22◦ C; denoted by superscript h in Table 1). This could be due to rapid relaxation of the cis isomer, leading to very little accumulation under steady-state illumination, or the complete absence of excited-state switching, and hence zero quantum yield. We cannot say for certain what the root cause is. However, the latter would be consistent with our prediction of zero quantum yield

for compound 1, which we had considered might be a simulation artifact. Derivatives of compound 2 (Fig. 2) exhibited photoswitching under all conditions, and had low but nonzero quantum yields in acetonitrile. The yields varied from 1.7% to 4.5%, which is up to 18× lower than the yield of 31% for unsubstituted trans-azobenzene in acetonitrile 26 . It is interesting that the parent compound bearing a bulky tert-butyl group had a predicted yield of zero, while its closest derivative, 2a, had a measured yield of 1.7%. The two values are not directly comparable, but it is still notable that replacing the tert-butyl group with a methyl group did not produce a high yield. Complete removal of the methyl group in 2e likewise did not alter the yield. Together, these results suggest a complex and somewhat unintuitive picture of the excited-state dynamics, one that nevertheless appears to be captured by our simulations. Photostationary state composition. PSS compositions were determined by 519 nm irradiation of the sample at −78 ◦ C in CDCl3 and subsequent measurement with low-temperature 1 H NMR. The trans:cis ratios varied substantially across the compound series. Compound 1 achieved a favorable ratio of 1:2.03, indicating that approximately 67% of molecules occupy the cis state at PSS. Compound 1a similarly showed good switching (1:1.59, ∼61% cis). In contrast, ortho-dihalogenated derivatives 2g (dichloro) and 2h (dibromo) exhibited poor PSS enrichment (1:0.15 and 1:0.14, respectively). The origin of this poor switching is not entirely clear. While the larger van der Waals radii of Cl and Br relative to F may play a role, ortho-tetrafluoroazobenzenes are known to achieve excellent PSS ratios despite bearing four ortho substituents. 6 Electronic effects, such as differences in how these halogens modulate the n→ π ∗ band separation of the trans and cis isomers, may also contribute. Photostability. Fatigue resistance was assessed through repeated irradiation–relaxation cycles in acetonitrile (ESI Figures S74–S86). Most compounds showed excellent photostability over 10 cycles with no detectable degradation. However, compounds 1 and 1a exhibited gradual baseline drift without irradiation in the presence of 5 ppm triethylamine, suggesting possible side reactions under prolonged basic conditions. After background subtraction, the irradiation–relaxation cycles revealed only minor drift, indicating that the neutral forms are likely to be photochemically robust. For biological applications under physiological pH, this base-mediated degradation pathway is unlikely to be relevant.

activity assay (BPS Bioscience). All measurements were made under green light illumination at 519 nm. This is an 80 nm red-shift relative to typical measurements under blue light, and was made possible by the compounds’ red-shifted absorption spectra. Given the rapid thermal relaxation observed in aqueous buffer (t1/2 < 0.01 min for all compounds in PBS/Tween 80), we employed continuous in situ illumination throughout the assay incubation to maintain a steady-state population of the cis isomer. A custom 3D-printed biophotoreactor was designed to uniformly irradiate half of a 384-well plate while leaving the remaining wells dark as internal controls (ESI Figures S6–S11). Solubility posed a significant challenge under assay conditions. At the screening concentration of 500 µM in PARP1 assay buffer containing 1% DMSO, several compounds exhibited visible precipitation. Addition of 0.01% Tween 80 substantially improved solubility for most compounds (ESI Figures S99–S101), and this formulation was adopted for all subsequent IC50 determinations. Among the ten synthesized compounds, 1 and 2a exhibited the most pronounced light-dependent inhibition (Fig. 6). Under dark conditions, compound 1 (predominantly the trans isomer) displayed an IC50 of 208.8 ± 28.3 µM. Upon continuous 519 nm irradiation, the IC50 decreased to 14.4 ± 1.9 µM, representing an approximately 15-fold enhancement in potency. This result qualitatively aligns with the computational prediction of stronger cis binding (Kd = 466 µM) relative to trans (Kd = 9.8 mM), though the absolute values differ. Direct numerical comparison between the computed Kd values and the experimentally measured IC50 values is not straightforward, as these are related but distinct quantities. The IC50 depends not only on the intrinsic binding affinity but also on assay-specific conditions such as substrate concentration, enzyme concentration, and the mechanism of inhibition, as formalized by the Cheng–Prusoff relationship 27 . Consequently, IC50 values are generally expected to exceed Kd values, and differences in absolute magnitude between the two metrics do not necessarily indicate poor predictive accuracy. The relevant comparison is therefore the directional agreement—namely, that the cis isomer is a substantially stronger inhibitor than the trans isomer— which is consistent between the computational predictions and experimental measurements. Compound 2a, the second computationally prioritized hit, showed more modest photoswitching: IC50 values of 209.4 ± 51.9 µM (dark) and 39.7 ± 4.7 µM (519 nm), corresponding to an approximately 5-fold potency enhancement. Compound 1a exhibited a similar trend with IC50 (PSS) = 167.8 ± 23.0 µM under illumination, although the dark IC50 could not be reliably determined due to poor curve fitting at the accessible concentration range.

Biological Evaluation To assess whether the computationally predicted differential binding translates into measurable lightdependent PARP1 inhibition, we evaluated the synthesized compounds using a chemiluminescent PARP1

Table 1: Predicted and experimental physicochemical, photophysical, and biological properties of synthesized compounds.

Structure

pKa ± SDg

trans λ5%,max (nm) pred. exp.

t1/2 (min)

pred.

exp.

6.8

3.78 ± 0.09

514

5.0

4.29 ± 0.09

513

1.3

< 0f

601

576n 578L

17.2

4.6

4.15 ± 0.10

603

616n 621L

−0.1

< 0f

613

0.4

< 0f

−0.4

PARP1 Inhibition (trans; µM) pred.r exp.s

PARP1 Inhibition (PSS; µM) pred.r exp.s

exp.

0±0

1.2 ± 0.2j

1:2.03

9.8 208.8 ×103 ± 28.3

466

–

5.2 ± 1.7j

1:1.59

5.0 ×103

NDe

1.1 167.8 ×103 ± 23.0

5.95 ± 0.89n 3.30 ± 0.19o < 0.01a,L

–

1.7 ± 0.1n

1:0.78

44.6

209.4 ± 51.9

2.9

39.7 ± 4.7

0.2

< 0.01a,m < 0.01a,L

–

NDc

–

NDe

1.8 ×105

NDe

577n 579L

20.7

1.80 ± 0.22n 1.37 ± 0.18o < 0.01a,L

–

4.5 ± 0.1n

NDb

356

NDe

1.5 ×104

NDe

582

562n 562L

9.3

1.62 ± 0.23n 0.93 ± 0.39o < 0.01a,L

–

NDd

1:0.48

–

232.7 ± 20.3

–

209.2 ± 7.9

< 0f

584

563n 570L

36.3

3.32 ± 0.42n 1.29 ± 0.24o < 0.01a,L

–

1.7 ± 0.3n

1:0.74

297

NDe

−1.8

< 0f

603

580n 614L

1.0

0.45 ± 0.01n 0.16 ± 0.01o < 0.01a,L

–

2.6 ± 0.2n

1:0.47

2.4 ×105

NDe

4.5

93.2 ± 11.2

−1.5

< 0f

644

562n 570L

–

0.06 ± 0.01n < 0.01a,L

–

4.3 ± 0.3n

1:0.15

3.5 ×104

NDe

2.0 ×104

NDe

−1.3

< 0f

640

557n 560L

–

0.07 ± 0.01n < 0.01a,L

–

4.0 ± 0.1n

1:0.14

3.0 ×104

NDe

13.8

144.3 ± 28.4

627h 567j 572L

1.0

1.8

exp.

PSS (trans: cis)q exp.

pred.

606h 556j 594L

pred.

QY (%)p

N/Ab,h 0.03 ± 0.01i 0.02 ± 0.01j < 0.01a,L N/Ab,h 0.17 ± 0.08j 0.08 ± 0.01k < 0.01a,L

14.4 ± 1.9

6.4 78.6 ± ×10−7 9.7

Cannot be determined precisely with our experimental setup. b No switching observed. c Not determined; switching too fast to reliably measure with our setup. d Not determined due to observed precipitation. e Not determined due to unreliable fit. f Estimated value. g Thermodynamic acidity constant (related to zero ionic strength) determined by capillary electrophoresis; see ESI for full details. h ACN + 5 ppm TFA, 22 ◦ C. i ACN + 5 ppm Et3 N, 10 ◦ C. j ACN + 5 ppm Et3 N, 22 ◦ C. k ACN + 5 ppm Et3 N, 37 ◦ C. L PBS/Tween 80, 22 ◦ C. m ACN, 4 ◦ C. n ACN, 22 ◦ C. o ACN, 37 ◦ C. p Quantum yields were determined at 519 nm spectrophotometrically against fulgide actinometer. q PSS (trans/cis) ratios were determined by irradiation of 1.5 mM CDCl3 solutions with 519 nm diode at -78 ◦ C, and subsequent measurement in low-temperature 1 H NMR. r Kd values. s IC50 values were determined against recombinant human PARP1 using a chemiluminescent activity assay (BPS Bioscience PARP1 Chemiluminescent Assay Kit, cat. no. 80551) according to the manufacturer’s protocol (see ESI). IC50 (trans) corresponds to dark-adapted samples (predominantly trans), and IC50 (PSS) corresponds to samples irradiated in situ at 519 nm. IC50 values represent the mean ± SD of technical triplicates (n = 3) from two independent biological replicates (N = 2). Irradiation was performed in a custom-built 3D-printed photoreactor; see ESI for reactor details and irradiation parameters. Photophysical values are mean ± SD of technical triplicates (n = 3) from a single synthesized batch (N = 1). Abbreviations: exp, experimental; pred, predicted; ACN, acetonitrile; TFA, trifluoroacetic acid; Et3 N, triethylamine; PBS, phosphate-buffered saline; PSS, photostationary state; QY, quantum yield. Predictions: The quantum yield was only predicted for compound 1, since predicting it for the other compounds would require training a new set of models with new quantum chemical data generated through active learning. Dashes for other predictions indicate failed calculations. PSS prediction was not part of the computational pipeline, and predicted PSS PARP1 inhibition is simply the predicted cis Kd .

Several compounds displayed minimal or no measurable differential activity. For compound 2g, reliable IC50 values could not be determined under either condition. For 2h, the dark IC50 could not be reliably determined, although an illuminated IC50 of 144.3 ± 28.4 µM was obtained, precluding quantitative assessment of light-dependent potency shifts. Moreover, their poor PSS enrichment (trans:cis ≈ 1:0.15) means that irradiation generated only a small increase in the proportion of the cis isomer. It is therefore difficult to determine whether both isomers are weak inhibitors, or whether the cis isomer is a moderate inhibitor whose activity is masked by its low proportion in the PSS. It is important to note that the measured IC50 values for all compounds remain in the micromolar range, substantially weaker than clinically used PARP1 inhibitors such as olaparib (IC50 ≈ 5 nM under our assay conditions; ESI Figure S102). This modest potency is typical of early-stage hits identified through computational screening and would require further optimization in a hit-to-lead campaign. Nevertheless, the demonstration of photoswitchable PARP1 inhibition under green light—particularly the 15-fold potency enhancement observed for compound 1—validates the core photopharmacological concept and the utility of the computational screening funnel for identifying differentially binding photoswitches under red-shifted illumination.

More significantly, thermal relaxation in aqueous buffer was dramatically accelerated (t1/2 < 0.6 s for all compounds in PBS/Tween 80), a limitation not captured by implicit-solvent calculations. Despite this rapid aqueous relaxation, compound 1 exhibited a ∼15fold enhancement in PARP1 inhibition upon green light irradiation at 519 nm (IC50 = 208.8 ± 28.3 µM in the dark vs. 14.4±1.9 µM under illumination), demonstrating that functional photopharmacology remains achievable under continuous irradiation conditions. Together, these results validate computation-guided multi-property optimization as a viable strategy for photopharmacology, while highlighting the need to explicitly model aqueous-phase dynamics and to develop scaffolds with improved thermal stability in physiological media. Future efforts should focus on extending cis-isomer lifetimes in aqueous environments, improving target potency into the nanomolar range, and exploring additional red-shifted chromophores compatible with deeper tissue penetration.

Computational Methods Library Generation The initial molecular library was generated using inhouse scripts. As described in CSI Sec. 2, we used a set of common literature substitution patterns with substituents from both the literature and the Enamine REAL catalogue 28 . We also expanded the library using a set of random fragment additions, and a set of modifications of promising motifs. QED and SA scores were computed with RDKit 29 , PAINS and MCF calculations were performed with MOSES (https://github.com/ molecularsets/moses), and SC scores were computed with code from https://github.com/connorcoley/ scscore.

Conclusions In this work, we combined large-scale computational screening with experimental synthesis and validation to identify photoresponsive PARP1 inhibitors. From an initial library of 5 million azobenzene derivatives, our computational funnel—comprising docking, neural network force field simulations, and free energy perturbation calculations—prioritized candidates predicted to exhibit red-shifted absorption, appropriate thermal half-lives, and isomer-dependent binding. We synthesized 10 compounds and experimentally characterized their photophysical and biological properties. The photophysical measurements validated several computational predictions: all compounds exhibited redshifted absorption relative to unsubstituted azobenzene (λ5% = 557–616 nm), and thermal half-lives in acetonitrile fell within the computationally predicted range of seconds to minutes. However, we also observed critical discrepancies. The experimental pKa value for compound 1 was substantially lower than predicted (3.78 vs. 6.8). This indicates that our computational model can overestimate azonium stability for critical compounds, even when results for other derivatives are quite accurate (Table 1). This highlights an inherent risk of computational screening: by screening for exceptional properties, one can unintentionally prioritize artifacts generated by model error over truly exceptional candidates.

Binding Free Energy Calculation Docking was performed with Dockstring version 0.2.0 30 . We removed all compounds that did not maintain their cis-trans isomerism during docking. Absolute FEP was performed with Yank version 0.25.2 31 . For FEP we used the ff14SB force field for the protein 32 , the TIP4P-Ew model for water 33 , and the GAFF2 force field 34 with re-parameterized azo terms for the ligands (CSI Sec. 5.3.3). To generate starting structures for FEP, we began with the docked pose and then performed a minimization in Gromacs 35 , followed by 250 ps of MD in the NVT ensemble and 500 ps in the NPT ensemble with restrained protein heavy atoms. We then performed FEP for 10 ns. An equilibration time period was determined by Yank and removed from the analysis 36 . λ windows for FEP were also determined automatically; on average about 70 windows were used for uncharged species and 100 windows for charged ones. Replica exchange among the windows was performed automatically 37 . The total binding free energy was computed with MBAR 38 . FEP 9

calculations took an average of 4.8 days per uncharged molecule using four GeForce RTX 2080 Ti GPUs. CSI Table S6 provides a benchmark of FEP and docking against experimental data. Both perform rather well, with FEP outperforming docking (e.g. R2 and Spearman ρ are 0.32 and 0.52 for docking, and 0.69 and 0.81 for FEP). FEP predictions were shifted by +4.81 kcal/mol to match the mean of the experimental data. This shift was also applied to virtual screening results, such that Kd = exp((∆GFEP + shift)/kB T ), where kB is Boltzmann’s constant and T is the temperature. Further details of docking can be found in CSI Sec. 5.2, and details of FEP can be found in CSI Sec. 5.3.

wavelengths. These geometries were added to the training set and the model was re-trained. 83,481 geometries were used in total. Performance statistics for species outside the training set are given in CSI Table S10. The model is quite accurate: the MAE for absorption wavelengths is 3.51 nm for azobenzene compounds and 7.03 nm for azonium derivatives (R2 = 0.98 and 0.97, respectively). The MAE of the transition dipole moment is 0.05 Debye for azobenzene and 0.16 Debye for azonium (R2 = 0.98 for both). The predicted spectrum is in good agreement with experiment for unsubstituted azobenzene. The predicted peak and extent are approximately 450 and 550 nm, respectively, which agree well with the spectrum measured in Ref. 6 in acetonitrile. However, the performance is significantly worse for azonium derivatives. For example, the predicted peak and extent are approximately 509 nm and 604 nm for the azonium compound 22 in Ref. 5 ; the experimental values in water are 597 nm and 750 nm 5 . By contrast, the predicted pKa for this compound is in good agreement with experiment (predicted and experimental values of 6.6 and 6.7, respectively). This suggests that the problem comes from TDDFT, not from incorrect sampling of ground state geometries. A final NFF was trained on ground and excited state energies and forces to perform NAMD simulations. This model was trained on MRSF-TDDFT data generated through active learning. We used a similar active learning setup to the one shown in CSI Fig. S5, with slight differences in the sampling strategies (see CSI Sec. 7.2). Surface hopping was performed with the Zhu-Nakamura method 42 , which does not require nonadiabatic coupling vectors. We performed 15 ps of ground state MD, followed by 200 independent excited state simulations. Excited state simulations were initialized from randomly sampled MD geometries and propagated for 4 ps. The quantum yield was computed as the proportion of trajectories that ended in a different isomer. The uncertainty in the yield was computed through statistical bootstrapping with 1,000 samples. The predicted yield of trans azobenzene was 43%, in good agreement with the experimental value of 35% in an 80-20 ethanol-water mixture 26 .

Chemical Property Calculations The ground state NFF was trained through active learning (CSI Fig. S5). To make use of previously generated data, we first pre-trained the models with 680,736 gasphase spin-flip (SF) TDDFT 39 calculations from Ref. 16 . We then refined the model with implicit solvent SFTDDFT data from Ref. 17 , plus newly generated data from active learning. Calculations were performed with Q-Chem 5.3 40 . This procedure yielded 90,460 calculations in total. However, we found that spin contamination was severe for azonium TSs, and so we finetuned the model with data from spin-complete MRSFTDDFT. We performed MRSF-TDDFT calculations for 17,027 geometries using GAMESS 41 . Activation free energies and pKa values were computed as described in CSI Secs. 6.1 and 6.2. The pKa predictions are compared to experimental data in CSI Fig. S13(a). The predictions are quite good, with an R2 value of 0.81, and a mean absolute error (MAE) of only 0.62. Half-life performance statistics for species outside the training set are given in CSI Table S9. The model is quite accurate for azobenzene derivatives, with an MAE of 0.77 kcal/mol relative to MRSF-TDDFT for the activation energy (R2 = 0.88). It is less accurate for azonium, with an MAE of 1.92 kcal/mol. This is mainly due to a subset of geometries with very large errors. Removing errors above 5 kcal/mol, which make up 8% of the test set, yields an MAE of 1.20 kcal/mol (R2 = 0.90). Note also that we performed a benchmark against experiment in Ref. 17 .

Graph-to-Property Model Training Graph-based models were trained to predict isomerization barriers, absorption wavelengths, and pKa values from NFF simulations. We used the Attentive FP architecture 43 implemented in the DeepChem library 44 . Model performance is given in CSI Table S12. The models are reasonably accurate, with R2 values between 0.54 and 0.74 depending on the property and training generation.

Photophysical Property Calculations The NFF for absorption spectra was trained on excitation energies, gradients, and oscillator strengths from TDDFT with the CAM-B3LYP functional. Geometries were selected from metadynamics simulations of non-TS geometries, which were performed during active learning for half-life predictions. The initial training set was selected randomly. A model was then trained and used to select geometries with high predicted absorption

Experimental Methods

with a Teflon-coated stirring bar, capped, and placed into a spectrometer-photoreactor. Spectra were accumulated, and the cuvette was irradiated with an LED until at least 10 full scans were obtained (∼0.5 minute). LEDs were then turned off and the sample was left to thermally re-equilibrate (at least 10 × t1/2 ). Values are reported as mean ± SD of technical triplicates (n = 3) from a single synthesized batch (N = 1). For mathematical details and representative graphs, see the ESI, section Calculation of Experimental Photophysical Properties, pages S183–S185.

Photostationary State Ratio Measurement The sample (∼1 mg, weighed on a microbalance) was weighed into a 1.5 mL HPLC screw cap vial and diluted with CDCl3 through a dilution series, resulting in a final analyte concentration of 1.5 mM in approximately 500 µL of CDCl3 . The tube was placed in an NMRphotoreactor, cooled to −78 ◦ C, and irradiated with an appropriate LED. After a defined time, the tube was quickly transferred to a pre-cooled NMR at −20 ◦ C. The collected spectra were immediately evaluated, and the irradiation/measurement procedure was repeated until equilibrium was achieved. The custom-made, 3Dprinted NMR photoreactor is described in more detail in the ESI, section NMR-photoreactor, pages S38–S39.

Chemiluminescent Assay for the Determination of PARP1 Inhibition Constants The PARP1 chemiluminescent assay was optimized for screening in 384-well plates. White 384-well plates with transparent bottoms were used for this purpose. Prior to performing the assay, plates were prepared by charging the wells with Ni2+ ions to capture 6His-tagged histones. Briefly, 20 µL of 10 mM BCML solution in 0.1 M sodium phosphate buffer, pH 8.0, was added to each well and incubated overnight at room temperature. After incubation, each well was washed three times with 200 µL of 0.05% Tween 20. Subsequently, the plates were blocked by incubating with 3% BSA in 50 mM Tris-HCl buffer (pH 7.5, containing 150 mM NaCl and 0.05% Tween 20) for 2 hours at room temperature. After blocking, the plates were sequentially washed with a series of buffers in the following order: first with 50 mM Tris-HCl (pH 7.5, containing 500 mM imidazole and 0.05% Tween 20), then with 0.05% Tween 20, followed by 100 mM EDTA (pH 8.0), and finally again with 0.05% Tween 20. Next, the plates were incubated with 10 mM NiSO4 solution for 20 minutes at room temperature, followed by washing first with 0.05% Tween 20 and then with 50 mM Tris-HCl buffer containing 500 mM NaCl, pH 7.5. Finally, plates were washed three times with PBS and then immediately used for the PARP1 assay. PARP1 activity was evaluated using the PARP1 Chemiluminescent Assay Kit according to the manufacturer’s instructions, with certain modifications. First, 12.5 µL of 1X histone mixture was added to each well, and the plate was incubated overnight at 4 ◦ C. After incubation, the plate was washed three times with PBS-T. Next, wells were blocked by adding 50 µL of blocking buffer per well and incubating at room temperature for 90 minutes. After blocking, 6.25 µL of master mix (containing PARP1 assay buffer, PARP substrate mixture 1, activated DNA, 0.5 mM DTT solution, and 0.01% Tween 80) was added to each well. At this step, the plate was placed in a custom-made, 3D-printed biophotoreactor (described in more detail in the ESI, section Bio-photoreactor, pages S35–S38), which illuminated half of the plate at either 519 nm or 590 nm, while the other half of the plate remained in the dark. After a 10-minute illumination period, 1.25 µL of com-

UV-Vis Spectrophotometry and Measurement of Thermal Isomerization A stock solution (∼3 × 10−4 M) was prepared as follows: the sample (∼1 mg, weighed on a microbalance) was placed into a 10 mL volumetric flask and diluted with dry spectroscopy-grade acetonitrile. The stock solution was transferred into a 10 mL amber screw-cap vial, sealed with vinyl duct tape, and stored under argon in a freezer at −20 ◦ C. Samples were prepared by diluting the stock solution directly into a fluorescence quartz cuvette just before measurement. The cuvette, equipped with a Teflon-coated stirring bar, was capped and placed into a spectrometer-photoreactor. Spectra were recorded while irradiating the sample with an LED until a photostationary equilibrium was reached. Subsequently, the LED was turned off, and the thermal relaxation process was monitored spectroscopically until thermal equilibrium was attained. Absorption coefficients for the pure (trans) isomer were calculated from the initial spectral scan. Thermal relaxation half-lives (t1/2 ) and rate constants (k) were determined by fitting the spectral data obtained after LED irradiation was ceased to first-order kinetic models. Values are reported as mean ± SD of technical triplicates (n = 3) from a single synthesized batch (N = 1). Detailed mathematical procedures, schematics of spectrometer-photoreactor, and representative graphs illustrating the thermal isomerization behavior for each compound are available in the ESI, section spectrophotometer-photoreactor, pages S39–S40.

Spectrophotometric Quantum Yield

Determination

Material (∼1–2 mg) was placed into a 4 mL amber screw-cap vial and dissolved in spectroscopy-grade acetonitrile (4 mL). The sample was transferred into a fluorescence quartz cuvette. The cuvette was equipped 11

sulfoxide (DMSO) as an electroneutral osmotic flow marker. Compounds exhibiting positive charge under these conditions were selected for further analysis. Selected compounds were subsequently analyzed in ten individually prepared background electrolytes: covering a pH range from 1.71 to 12. All buffers, except the neat NaOH solution (pH = 12), were adjusted to an ionic strength of 25 mM and equilibrated at 25 ◦ C before use. From the pH-dependent effective electrophoretic mobilities (µeff ), mixed acidity constants (pKa,mix ) and actual ionic mobilities of the cationic forms (µ1+ ) were determined. These mixed constants, which relate to the activity of hydronium ions at the experimental ionic strength of 25 mM, were subsequently converted to thermodynamic acidity constants (pKa,th ) using the Debye–Hückel theory as described. 45 Values are reported as mean ± SD of technical triplicates (n = 3) from a single synthesized batch (N = 1). Electropherograms, pH dependence of effective mobilities, pKa1,th , pKa2,th , pKa1,mix , pKa2,mix , and µ1+ are presented in the ESI, Table S2, Figures S87–S98.

pound solution at the desired concentrations (ranging from 2 µM to 500 µM) was added to wells containing the master mix. The concentration of DMSO was adjusted to 1% in all samples, including the blank and positive controls. The reaction was initiated by adding 5 µL of PARP1 enzyme solution (final reaction concentration: 0.132 ng/mL). Positive controls corresponded to reactions containing enzyme but no tested compound, while blanks corresponded to reactions without enzyme or compound. After incubation for 1 hour at room temperature under either illuminated or dark conditions, the plate was washed three times with PBS-T. Detection was performed using streptavidin-HRP and an ELISA ECL substrate. Chemiluminescence was measured using a plate reader (Tecan Spark, Switzerland). Blank values were subtracted from all measured signals, and resulting values were normalized to the positive control, defined as 100%. IC50 values were calculated using a custom Python script by fitting data to a four-parameter sigmoid curve using the curve fit function from the SciPy module, and plots were generated using the Matplotlib library. Data are reported as mean ± SD of technical triplicates (n = 3) from two or one independent biological replicates (N = 2 for measurements at 519 nm and N = 1 for measurements at 590 nm). All dose response graphs can be found in the ESI, Figures S102–S112.

Acknowledgments This work was financially supported by DARPA (Award HR00111920025) and the Institute of Organic Chemistry and Biochemistry of The Czech Academy of Sciences (RVO 61388963). Harvard Cannon cluster, MIT Engaging cluster, and MIT Lincoln Lab Supercloud cluster at MGHPCC are gratefully acknowledged for computational resources and support. We also thank Dr. Jiřı́ Kaleta, Dr. Tomáš Slanina, and Dr. Milan Mašát for their valuable photochemical consultations. We gratefully acknowledge Dr. Lucie Bednárová for her consultations regarding vibrational spectra; Dr. Tomáš Mašek for his assistance with retrosynthetic design; Dr. Robert Reiberger for his deep Buchwald–Hartwig amination expertise; Dr. Karel Kudláček for assisting with solubility studies. Dr. Veronika Vetýšková and Dr. Jan Voldřich for their consultations regarding the PARP1 assay; Dr. Ema Chaloupecká for assistance during NMR measurements; Dr. Ondřej Ticháček for his expert consultation on mathematical implementation in Python scripts; Dr. Jakub Copko and Dr. Eva Bednářová for their assistance with UV-Vis spectroscopy; Dr. Jaroslava Hniličková for measuring Karl Fischer titrations; and Jakub Brož for his help with CAD modeling and 3D printing.

Determination of Solubility in PARP1 Buffer Compounds (500 µM final concentration; 12.5 µL per well in a 384-well plate) were prepared in 1× PARP1 assay buffer supplied with the BPS Bioscience chemiluminescent PARP1 activity kit (cat. #80551). According to the manufacturer, the exact buffer formulation is proprietary; however, from our experience such buffers typically contain 50 mM Tris-HCl, 150 mM NaCl, 5 mM MgCl2 , and 1 mM DTT at pH ≈ 7.5. All samples contained 1% (v/v) DMSO; where indicated, Tween 80 was added to a final concentration of 0.01% (v/v). Samples were incubated for 1 hour at room temperature, irradiated at 519 nm, and imaged with a Zeiss Axio microscope (50× magnification) using AxioVision SE64 Rel. 4.9.1. The resulting photomicrographs are shown in the ESI, Figures S99–S101.

Capillary Electrophoretic Determination of Acidity Constants Compounds were dissolved in methanol to achieve final concentrations of 0.16–0.25 mg/100 µL. Samples were introduced into the capillary via hydrodynamic injection using pressures of 13.8–20.7 mbar for 10–20 seconds, resulting in injection volumes of 5.1–15.4 nL. Detection was performed using a PDA (200 nm). Initial screening of all compounds was conducted in acidic background buffer (pH 2.01) using dimethyl 12

Computational Supplementary Information 1

Properties of top candidates

Here we provide FEP results and quantum chemistry data for the top candidates. CSI Fig. S1 shows the candidates that reached the final stage of screening, before being filtered based on NAMD and FEP results. Compounds 14 and 20 were generated manually based on 1. CSI Table S2 shows the FEP results for molecules 1 to 20. CSI Table S1 contains chemical and photophysical data for the candidates with predicted Kd < 5 mM.

Molecule generation

2.1

Base library

Following our previous work 17 , we combined common literature functional groups and substitution patterns to create a virtual library of azobenzene derivatives. We also expanded the library and filtered the molecules based on druglikeness and synthetic accessibility (Fig. 1(a)). To expand the library, we fused five- and six-member rings to the ortho and meta positions of benzene in azobenzene (CSI Fig. S2). We used the conjugated rings in Ref. 46 , plus the six-member di-methoxy ring in Ref. 5 . Fusing conjugated rings to benzene can red-shift the absorption spectrum by expanding the π conjugated system. The di-methoxy ring can red-shift absorption by stabilizing the azonium species 5 . The oxygens act as hydrogen bond acceptors, which stabilize azonium and hence increases its presence in solution 5 . We further expanded our library by adding drug-like fragments from the Enamine REAL catalogue 28 . We used all compounds with 11 heavy atoms or fewer in the Essential Fragment Library and the High Fidelity Fragment Library. We also used the basic fragments from Fig. 2-4 in Ref. 47 . 1

N N

N O

F F F

O H 2N S O

N N

N N N

N N

O S NH2 O

8 N

H 2N

N S

10 N

N N

N F N

NH2

O N

H 2N

NH2+ HN

N N

N O

S O

O H 2N S O

N N

N H 2N

N N N

N N

O N

N HO

HN F

F N O

19 N N

H 2N

O N

O H 2N S O

NH2

N+

N N

O N

NH2

O N N

N O

Supplementary Figure S1: Candidates that reached the final stage of screening (NAMD and FEP).

O O

Supplementary Table S1: Chemical and photophysical predictions for compounds with predicted cis binding affinities below 5 mM. ∆G† is the activation free energy for thermal isomerization, and ∆G†calib is the calibrated activation free energy, shifted so that the value for azobenzene matches experiment. ∆Gprot is the free energy cost of protonating the eff † molecule, and ∆G†, calib is the effective value of ∆Gcalib after adding the protonation cost. Compounds are separated into azobenzene and their two azonium forms. We have underlined the lowest value of ∆G†,eff calib and its associated half-life. An asterisk indicates that TS optimization did not converge, and so the azonium barrier was simply taken as the barrier of the other azonium form. Units are kcal/mol for free energies, minutes for half-lives, and nm for absorption wavelengths. Calculations were performed using MRSF-TDDFT with BHHLYP applied to NFF critical points for all quantities other than λ. For λ we used TDDFT with the CAM-B3LYP functional. Further details can be found in CSI Secs. 6.1, 6.2, and 7.1. cis Compound 1, azobenzene 1, azonium a 1, azonium b 2, azobenzene 2, azonium a 2, azonium b 3, azobenzene 3, azonium a 3, azonium b 7, azobenzene 7, azonium a 7, azonium b 8, azobenzene 8, azonium a 8, azonium b

†

∆G 26.2 15.3 26.9 23.5 23.7∗ 23.7 26.3 20.0 19.7 23.8 21.7∗ 21.7 26.2 20.1∗ 20.1

∆G†calib 23.6 12.7 24.2 20.9 21.1 21.1 23.7 17.4 17.1 21.2 19.1 19.1 23.6 17.5 17.5

pKa – 1.4 3.9 – 3.4 −3.7 – −5.6 −7.4 – 3.4 3.5 – 4.4 −0.8

∆Gprot – 8.2 4.8 – 5.5 15.2 – 17.7 20.2 – 5.4 5.3 – 4.1 11.2

trans eff ∆G†, calib

23.6 20.9 29.1 20.9 24.3 34.2 23.7 35.1 37.3 21.2 24.5 24.3 23.6 21.6 28.7

t1/2 84 1.0 5.9 × 105 1.0 2.8 × 102 2.5 × 109 101 1.0 × 1010 3.9 × 1011 1.7 3.3 × 102 2.8 × 102 87 3.4 3.2 × 105

pKa – 2.1 6.8 – −5.6 0.9 – −4.8 −0.5 – 0.8 1.5 – 0.3 4.8

λ5% max 534 667 514 617 657 530 649 534 510 619 531 712 600 643 460

Supplementary Table S2: FEP binding predictions for candidates that reached the final stage of screening (CSI Fig. S1). ∆GFEP is the binding free energy computed with FEP. ∆GFEP calib is the calibrated binding free energy, computed by adding a shift of 4.81 kcal/mol to match experiment (see CSI Table S6). Kd is the dissociation constant computed from ∆GFEP calib . Binding affinities below 5 mM are underlined. Trans binding scores are only shown for cis compounds with binding affinities below 5 mM. Compound 1 cis 1 trans 2 cis 2 trans 3 cis 3 trans 4 cis 5 cis 6 cis 7 cis 7 trans 8 cis 8 trans 9 cis 10 cis 11 cis 12 cis 13 cis 14 cis 15 cis 16 cis 17 cis 18 cis 19 cis 19 trans 20 cis

∆GFEP (kcal/mol) −9.4 −7.6 −9.2 −8.1 −8.5 −5.0 −5.1 −6.9 −3.8 −9.3 −10.0 −10.0 −10.5 −6.2 −3.6 −7.5 −4.8 −4.5 −2.6 −0.6 −7.8 −3.9 −7.1 −8.8 −7.9 0.4

∆GFEP calib (kcal/mol) −4.5 −2.7 −4.4 −3.3 −3.7 −0.2 −0.3 −2.1 1.0 −4.5 −5.2 −5.2 −5.7 −1.4 1.2 −2.7 0.0 0.3 2.2 4.2 −3.0 0.9 −2.3 −4.0 −3.1 5.2

Kd (µM ) 466 9.8 × 103 591 4.0 × 103 2.1 × 103 6.7 × 105 6.2 × 105 3.0 × 104 5.3 × 106 534 163 154 70 1.0 × 105 7.6 × 106 9.9 × 103 1.1 × 106 1.6 × 106 4.2 × 107 1.2 × 109 6.6 × 103 4.9 × 106 2.0 × 104 1.2 × 103 5.0 × 103 6.7 × 109

The substitution patterns are shown in CSI Fig. S3. We enumerated all possible combinations of patterns and fragments. For each generated molecule, we optionally added a small substituent with ≤ 4 atoms to a random carbon on the azobenzene scaffold. When adjacent ortho and meta positions were unsubstituted, we fused one of the rings in CSI Fig. S2 with benzene at up to two different sites. More fusions would bring the number of aromatic rings above four, which is already quite high for drug-like molecules (see Fig. 1g in Ref. 8 ). We then narrowed our search space to synthetically accessible drug-like molecules with a reasonable chance of having good pharmacokinetic properties. To do so we used several different filters. First we used the quantitative estimate of drug-likeness (QED) 8 , a continuous generalization of Lipinski’s rule of 5 7 that estimates drug-likeness on a scale of 0 to 1 (higher is better). We filtered out all molecules with QED < 0.4. We then removed compounds that failed any of the medicinal chemistry filters (MCFs) in Ref. 11 , apart from the N=N filter and the limitation of ≤ 3 halogens. The first filter was removed because N=N is present in azobenzene. This bond can be cleaved by azoreductase 48 , but may be metabolically stable if the molecule is appropriately substituted 49 . The halogen filter was removed because tetra-halogenation is a good red-shifting strategy 50 . We also removed groups containing panassay interfering compounds (PAINS) 9,10 . These groups often lead to false positives through promiscuous binding. Code for filtering by PAINS and MCFs was accessed from Ref. 11 . We also inspected a few hundred compounds by hand. We found a small number of non drug-like groups, such as aldehydes and nitroso groups, and removed all molecules that contained them. Lastly, we filtered molecules by synthetic accessibility. We used the synthetic accessibility (SA) score 12 and synthetic complexity (SC) score 13 . The SA score uses a number of chemistry-informed rules to estimate accessibility, and its scores are well-correlated with those of synthetic organic chemists 12 . The SC score uses a neural network

O N

H N

N N

Supplementary Figure S2: Groups fused to the azobenzene core. The fused bond is shown with a wavy line.

N R1

R2 N

N R2

R2 R2 R2

R2 N

N R1

R2 N

Supplementary Figure S3: Substitution patterns used for combinatorial molecule generation.

trained on published reactions. It is trained to implicitly recognize the number of steps needed to synthesize a compound 13 . While less accurate than the SA score 13 , it is nevertheless a useful second filter for reducing the number of compounds. We removed molecules with SA scores above 5 (scale of 1 to 10; lower is better) and SC scores above 4 (scale of 1 to 5; lower is better). After generating the candidate pool and filtering by drug-likeness and synthetic accessibility, we were left with 9.5 million molecules in total. Half of these were cis and half were trans. We then generated two azonium compounds for each species, one for each nitrogen protonation site. In some cases the two azonium compounds were equivalent by symmetry. Adding the azonium species yielded 26.8 million molecules in total. These molecules were sampled randomly in each round of active learning, while azobenzenes were sampled together with their two azonium species during screening.

2.2

Library expansion

After docking all 10 million compounds, we expanded our virtual library in three ways. First, we identified docking hits with at least one fused dioxane ring. These scaffolds are of interest because they were used in Ref. 5 to increase the trans pKa , and hence red-shift the absorption spectrum. We generated new species by converting some dioxane rings to methoxy groups in the ortho position (CSI Fig. S4(a)). This strategy was also used in Ref. 5 . We performed this step to increase the pool of hits with the potential for achieving near-neutral pKa . Second, for hits with one or two dioxanes fused to the same ring, we added two halogens to the ortho position of the other benzene ring (CSI Fig. S4(b)). We also generated all combinations of methoxy and dioxane for these compounds. We used fluorine, chlorine, and bromine as the halogen substituents. We added halogens because azonium compounds tend to have low barriers 5 , while ortho halogenation tends to increase barriers 6 . Ortho halogenation red-shifts the absorption spectrum, with the size of the effect increasing when moving down the periodic table 50 . Third, for hits with two ortho acetylpiperazine substituents, we added two ortho halogens on the opposite benzene ring. In our first round of screening we found some hits with this piperazine motif, similarly to our previous work 14 , and so we generated more candidates with this scaffold. We used ortho halogenation for the same reasons described above. We filtered the new molecules by drug-likeness and synthetic accessibility, as described above, and docked the remaining compounds. After filtering by docking score, we were left with 898 cis-trans pairs in total. All 898 molecules were then screened for chemical and photophysical properties, independent of the predictions of the

N O

O X N

N O

O O

X O

O O

O N

O O

N N N N

N N

N O

O O

OO O

X N

OO O

O O

X N

N O

X X X N N N N N O O X X X O O

N O

X X N N N O X X

X O

O X

O O

X X X N N N NO N O X X X

Supplementary Figure S4: Strategies used to expand the library after docking. (a) For all docking hits with at least one fused dioxane ring, all permutations of (dioxane → methoxy) and (dioxane → dioxane) were generated. (b) For a hit with one or two dioxanes fused to the same benzene, two halogens were added in the ortho positions on the other benzene. All permutations of (dioxane → methoxy) and (dioxane → dioxane) were also generated.

graph-to-property model.

Virtual screening

In the first NFF stage of screening, we began by randomly sampling 4,000 molecules. 1,000 were chosen to have hydrogen bond acceptors in the ortho position, or one atom away from it. These compounds received pKa and azonium barrier calculations. After performing NFF calculations on these molecules, we trained 2D models to predict the pKa , the maximum absorption wavelength at 5% of the spectrum peak, and the thermal half-life. We trained an ensemble of three models for each property. The models were then used to select 4,000 new species. A maximum of 50% were chosen based on optimal chemical and photophysical properties or high pKa values. This was called exploitation. 25% were chosen to maximize the model uncertainty, as measured by the variance of the ensemble predictions. 25% were chosen randomly. The latter two were called exploration; their purpose was to explore areas of chemical space in which the model was not well-trained. The selected molecules received new NFF calculations. Molecules chosen because of high trans pKa values, and molecules with ortho hydrogen bond acceptors, received pKa and azonium barrier calculations. We performed three rounds of screening and 2D model training. Only exploitation was used in the final round. We then validated the chemical and photophysical predictions with quantum chemistry. We kept all molecules with predicted half-lives between 5 seconds and 5 hours (raw ∆G† between 22 and 27 kcal/mol), and λ5% max above 600 nm.

ML-accelerated simulation

Properties were computed through ML-accelerated simulation, using methods that we developed in previous work 14,16,17 . To compute pKa values, we used an NFF to compute the free energies of the azonium and azo forms of the molecules. We then used a linear regression to fit free energy differences to experimental pKa values, as described in CSI Sec. 6.2. To compute absorption spectra, we performed MD with the same NFF, and computed oscillator strengths and absorption energies for each frame with a second NFF. The histogram of the absorption energies weighed by the oscillator strengths gave the spectrum, as described in CSI Sec. 7.1. For the quantum yield, we used a third NFF to perform NAMD on the top candidates. This model was trained on both S0 and S1 energies and forces, and was used with the Zhu-Nakamura surface-hopping method to predict isomerization quantum yields 42 . Thermal half-lives were computed with the ground state NFF using the workflow developed in Ref. 17 and shown in CSI Fig. S5(b). In Ref. 17 we argued that thermal isomerization in azobenzene occurs through rotation mediated by intersystem crossing (ISC). However, we also found a strong correlation between rates from ISC and rates from Eyring transition state (TS) theory using a rotational TS. Therefore, to avoid computing triplet forces to train a separate NFF, we simply used the ground-state singlet NFF and computed rates from Eyring TS theory. We

message

NAMD pre-training

&!"

message

…

update

&!"

TS generation

update

neural network committee

×2 mechanisms

TS generation

database X

trans SMILES

RDKit confgen

cis

trans

NN confgen

relaxed scan

cis ensemble

trans ensemble

TS guess

N Y

! = 𝐸𝜓 𝐻𝜓 ⃗ 𝑓 = −∇! 𝐸 quantum chemistry

cis SMILES

Y Y

NN confgen

? energy, uncertaintybased selection

optimized TS

EVF

TS ensemble

Supplementary Figure S5: Approach to active learning for TS generation. (a) Active learning loop for training the NN. (b) Workflow for generating equilibrium and TS geometries. “NN confgen” means conformer generation with the NFF using metadynamics. Initial conformations were generated with RDKit 29 and refined with NN conformer generation. TSs were generated through relaxed scans followed by conformer generation and eigenvector following (EVF). report half-lives below by subtracting a constant shift in ∆G† , equal to the difference between the computed and experimental values for azobenzene (28.0 kcal/mol and 25.4 kcal/mol in benzene 51 , respectively). The NFFs used for ground-state properties and NAMD were trained on mixed-reference spin-flip time-dependent density functional theory (MRSF-TDDFT) data 52,53 . We used the BHHLYP functional 54 , the 6-31G* basis 55 , D3BJ dispersion 56,57 , and a C-PCM description of water 58–60 . MRSF-TDDFT can account for some multi-reference effects, which are critical for NAMD 16 and for TS searches in azobenzene 17 . The excited-state NFF used for absorption spectra was trained with TDDFT using the CAM-B3LYP functional 61 , the def2-SVP basis 62 , and a C-PCM model of water 58–60 . Further quantum chemistry details can be found in CSI Secs. 6.1, 7.1, and 7.2. All NFF models used the PaiNN architecture 63 , which predicts molecular properties through equivariant messagepassing. This approach generates a feature vector for each atom that incorporates information from its surrounding environment. The initial feature vector is generated from the atomic number alone, and is then updated through a set of “messages” that incorporate the distance, orientation, and features of atoms within a cutoff distance. This process is repeated several times to gather information from increasing distances. Geometries for training the models were generated through active learning (CSI Fig. S5(a)). Three ground state models were used to generate TSs for a set of sampled reactions, and geometries from the TS workflow were sampled based on energy and uncertainty to receive new quantum chemistry calculations. The new data was added to the training set and the models were retrained. A similar process was used for NAMD. The excited-state model used for absorption spectra was trained on metadynamics samples of cis and trans equilibrium geometries. Further details of the training process and model accuracy can be found in Methods and CSI Secs. 6.1, 7.1, and 7.2. All models performed well relative to the quantum chemistry methods used for training (see Methods and CSI Tables S9, S10, and S11). With the exception of azonium absorption wavelengths, the models also performed well relative to experiment (see Methods, CSI Table S6, and CSI Fig. S13(a)).

Protein-ligand binding

5.1

PARP1

Previously we docked 9,000 azobenzene derivatives into 58 targets using Dockstring 14 . We showed that the docking scores of cis and trans isomers were highly correlated for most targets. We also found that trans ligands bind more strongly than cis to most proteins. Based on this work, we identified PARP1 as a promising photo-druggable protein

CNNC

CCNN CNN

CCN

Supplementary Figure S6: Internal coordinates of interest in azobenzene. Dihedrals are shown in blue and angles in red. The NN and CN bonds are those within the CNNC dihedral. Supplementary Table S3: New bond parameters for cis and trans azobenzene. Units are Å for distances and kcal/mol/Å2 for force constants.

cis trans

Atoms NN CN

r0 1.225 1.420

k 1792.543 669.216

NN CN

1.234 1.405

1732.791 717.017

Supplementary Table S4: New angle parameters for cis and trans azobenzene. Units are degrees for angles and kcal/mol/rad2 for force constants.

cis trans

Atoms CCN CNN

θ0 111.0 116.0

k 59.751 143.403

CCN CNN

116.0 110.5

95.602 167.304

for several reasons. The first is that it is a cancer target. The second is that cis isomers bind more strongly to PARP1 than trans isomers, which is uncommon in the proteins that we studied 14 . This is desirable because one usually wants the inactive drug to be thermodynamically stable (trans). If the inactive drug were thermodynamically unstable (cis), then it would eventually become active throughout the whole body due to thermal isomerization. The last reason is that docking scores for PARP1 are fairly reliable 30 . This makes it a good target for computational screening. It is especially important given that FEP is quite slow, and so it is important for a reasonable proportion of the docking hits to be actual hits.

5.2

Docking

Docking was performed with Dockstring 30 . This open-source package performs docking into one of 58 different proteins, and requires only the ligand SMILES string and protein name as input. It performs (de)protonation at physiological pH using OpenBabel 64 , generates initial ligand conformers using RDKit 29 , and performs docking with AutoDock Vina 15 . The protein structures come from the DUD-E database 65 . The search boxes were specified by hand by the authors of Ref. 30 for each protein. We docked all 9.5 million azobenzene derivatives into PARP1.

5.3

FEP

5.3.1

Theory

Docking scores were improved through absolute FEP calculations 66 . These calculations use a thermodynamic cycle with two stages. In the first stage, the ligand is decoupled from the protein and the solvent. A restraint is usually added to keep the uncoupled ligand close to the protein. In the second stage, the ligand is re-coupled to the solvent in the absence of the protein. The free energy change of the cycle, plus a correction for the restraint, is then the free energy difference between the solvated ligand and the bound ligand. The free energy differences are computed with MD runs from a series of windows, called λ windows. The Hamiltonians of these windows interpolate between the fully coupled system and the fully uncoupled system. An estimator such as the multi-state Bennett acceptance ratio (MBAR) 38 is used to determine the free energy change between the first and last window, using MD samples from all windows.

Supplementary Table S5: New dihedral parameters for cis and trans azobenzene. Units are degrees for dihedral angles and kcal/mol for force constants.

cis trans

Atoms CNNC CCNN

ϕ1 0.0 180.0

ϕ2 0.0 180.0

k1 0.0 0.0

k2 −9.560 1.076

CNNC CCNN

0.0 0.0

−8.485 −1.649

Supplementary Figure S7: Distribution of internal coordinates for cis azobenzene, computed with our NFF (solid blue) and our re-fit GAFF2 classical force field (hatched white).

5.3.2

Implementation

Absolute FEP calculations were performed with Yank (version 0.25.2) 31 . Yank is an open-source package that uses OpenMM to run MD 67 , determines λ windows and MD equilibration times automatically 36 , performs replica exchange among windows to accelerate sampling 37 , and computes the binding free energy with MBAR. It also provides convergence analysis and uncertainty estimates. We ran Yank for 10 ns per window, using an average of 70 windows for uncharged species (further details below). This took an average of 4.8 days per molecule using four GeForce RTX 2080 Ti GPUs. Each run was divided into segments of 1 ps of MD, followed by displacement and rotation of the ligand with Monte Carlo moves. Replica exchange was then performed among all windows using Gibbs sampling 37 . This process was repeated 10,000 times to generate equilibrium samples for all windows. For MD we used the NPT ensemble with the LangevinSplittingDynamicsMove function, which performs Langevin dynamics using the BAOAB split integrator 68 . We used a time step of 2 fs and a collision rate of 1 ps−1 . Velocities were sampled from the MaxwellBoltzmann distribution at the start of each segment. The pressure was enforced with the default Monte Carlo barostat in Yank. The temperature was set to 298.15 K and the pressure to 1 atmosphere. Windows were chosen automatically by Yank. Yank determines a schedule for first turning off electrostatic interactions and then turning off van der Waals interactions. In the protein-ligand stage, it also determines a schedule for turning on the restraint, chosen here to be harmonic. The electrostatic interactions were turned off by linearly scaling all ligand charges to zero. The van der Waals interactions were turned off by scaling Lennard–Jones terms to zero using a soft-core schedule 69 , with (a, b, c, α) in Ref. 69 set to (1, 1, 6, 0.5). All windows were chosen such that the standard deviation of energy differences between states were close to equal. These states were generated through an iterative algorithm that uses samples from short MD runs in each window. On average, about 70 windows were used in the protein-ligand stage for neutral species, and 100 were used for charged species. For the solvent stage, 60 and 90 windows were used for neutral and charged molecules, respectively. Long-range electrostatics were computed with the particle-mesh Ewald method 70 . We used a non-bonded cutoff of 9 Å and a box clearance of 10 Å. We used Na+ and Cl− ions to form a salt buffer of 150 mM, and to neutralize

Supplementary Figure S8: Distribution of internal coordinates for trans azobenzene, computed with our NFF (solid blue) and our re-fit GAFF2 classical force field (hatched white).

Supplementary Figure S9: Distribution of internal coordinates for cis azobenzene, computed with our NFF (solid blue) and the unmodified GAFF2 classical force field (hatched white).

any net charges in the system. 5.3.3

Force fields

We used the ff14SB force field 32 for the protein and the TIP4P-Ew model 33 for water. In our benchmark against experimental data, we used the GAFF2 force field 34 with AM1-BCC partial charges 71 . We used GAFF atom types instead of GAFF2 atom types; surprisingly, we found that this significantly improved the agreement with retrospective experimental data. Charges were computed on the lowest-energy structure generated by a conformer search. We performed conformer searches with CREST 72 for the benchmark and our NN version of CREST 17 for azobenzene screening. For prospective calculations of azobenzene derivatives, we used GAFF2 with a re-parameterization of the azo terms (see below), again with GAFF atom types. We also re-ran FEP with GAFF2 atom types and compared to the GAFF/GAFF2 predictions (not shown). We found that the top GAFF/GAFF2 binders were also

Supplementary Figure S10: Distribution of internal coordinates for trans azobenzene, computed with our NFF (solid blue) and the unmodified GAFF2 classical force field (hatched white).

good GAFF2/GAFF2 binders. However, there were some strong GAFF2/GAFF2 binders that were not good GAFF/GAFF2 binders (e.g. compounds 11 and 13). We chose not to focus on these compounds because of the lack of FEP consensus, and because of synthetic accessibility concerns. We re-fit the azo terms because GAFF2 gave poor results for azobenzene. The internal coordinates of interest for azobenzene are shown in CSI Fig. S6. We modified parameters for the NN and CN bonds, the CNN and CCN angles, and the CCNN and CNNC dihedrals. We used harmonic functions for bonds and angles, and periodic functions for dihedrals: 1 k(r − r0 )2 2 1 Vangle (θ) = k(θ − θ0 )2 2 2 X kn (1 + cos(ϕ − ϕn )). Vdihedral (ϕ) = Vbond (r) =

(1) (2) (3)

n=1

Here V is the potential energy, k is the force constant, r is the bond length, θ is the angle, and ϕ is the dihedral angle. The new parameters are given in CSI Tables S3-S5; they were determined by hand to reproduce distributions from NFF dynamics. NFF distributions are compared with re-fitted GAFF2 distributions in CSI Figs. S7 and S8. Statistics were collected from 1 ns simulations in vacuum using the Nosé-Hoover thermostat 73,74 at 298.15 K. Classical FF simulations were performed with Gromacs 35 , and NFF simulations were performed with ASE 75 . No bond constraints were used. Distributions from equivalent internal coordinates were aggregated in each plot (there are two equivalent CNN angles, four equivalent CCN angles, and four equivalent CCNN dihedrals). All distributions agree quite well. The optimized geometries also agree well: the cis and trans geometries have RMSDs of only 0.048 and 0.041 Å with respect to NFF geometries. By contrast, the agreement is poor for standard GAFF2 (CSI Figs. S9 and S10). Cis-azobenzene isomerizes to trans, while the CCNN dihedrals in trans azobenzene are centered at the wrong angles. The latter leads to non-planarity of the optimized trans molecule. The CCN distributions are also too narrow, and the CNNC distributions are too wide. Modified ligand force field parameters were supplied to Yank through frcmod files. To check that we rendered these files correctly, we used Yank functions to load their parameters, and used the results to create OpenMM simulations of azobenzene. We then ran 1 ns simulations in OpenMM in the NVT ensemble, and confirmed that the distributions matched the NFF results.

Supplementary Figure S11: Docking and FEP results for the PARP1 benchmark. 1 kcal/mol and 2 kcal/mol ranges are shown with dashed and solid gray lines, respectively. All results are shifted by a method-dependent constant, ∆Gshift , equal to the mean of the experimental results minus the mean of the predictions. The constants and the prediction statistics are given in CSI Table S6. (a) Docking and (b) FEP run for 10 ns. Supplementary Table S6: Statistics of docking and FEP predictions for the PARP1 benchmark. R is the Pearson R value, R2 is the coefficient of determination, and ρ is the Spearman rank correlation coefficient. 95% confidence intervals are given in square brackets below; they were computed with statistical bootstrapping using 1,000 samples. ∆Gshift , MAE, and the root-mean-squared error (RMSE) are given in kcal/mol.

Docking FEP (10 ns)

5.3.4

R 0.63 [0.18, 0.90] 0.86 [0.68, 0.97]

R2 0.32 [−0.21, 0.73] 0.69 [0.33, 0.89]

ρ 0.52 [−0.12, 0.92] 0.81 [0.49, 1.0]

MAE 1.68 [1.09, 2.27] 1.05 [0.74, 1.34]

RMSE 1.98 [1.39, 2.52] 1.19 [0.92, 1.44]

∆Gshift 1.64 4.81

Initial structures

To generate initial structures for FEP, we started with docked protein-ligand complexes, and then relaxed them with restrained MD. This was necessary to ensure stable Yank simulations, as starting with docked complexes often led to divergences. Relaxations were performed in Gromacs. We started with the docked structure of the ligand and a PARP1 structure from the PDB database (PDB ID 3L3M). This is the same structure used in DUD-E and hence in docking, apart from a few extra atoms that had been removed by the DUD-E curators. The structure was prepared using PDBFixer 76 . We performed 50,000 minimization steps with steepest descent, followed by 50,000 steps with conjugate gradient descent, each with a force convergence tolerance of 5 kJ/mol/nm. We then performed 250 ps of MD in the NVT ensemble at 298.15 K, followed by 500 ps of MD in the NPT ensemble at 1 atm. Protein heavy atoms were restrained with a force constant of 1,000 kJ/mol/nm2 . Gromacs simulations were performed with the ff14SB protein force field and the TIP3P water model 77 . The re-fitted GAFF2 force field was used with AM1-BCC charges for the ligand. Ligand files were converted from PDB to mol2 using OpenBabel 64 followed by Antechamber 78 . The ligand mol2 files were supplied to acpype 79 . The acpype script was locally modified to allow the GAFF2 force field to be used. We used acpype to produce Gromacs input files, which were then modified to account for changes in the GAFF2 force field. Long-range interactions were treated with PME, using a 10 Å cutoff for long-range electrostatics and a 10 Å cutoff for Lennard-Jones terms. All bonds were constrained with the LINCS algorithm 80,81 . The Nosé-Hoover thermostat 73,74 was used for NVT and NPT simulations, and the Berendsen barostat 82 was used for NPT simulations. Separate heat baths were used for the solute and the solvent. Frames from the NPT MD run were clustered using the single linkage method, with a clustering cutoff of 1 Å. The centroid of the largest cluster was used as the starting point for FEP. The structure was split into protein and ligand PDB files; the protein PDB file was supplied to Yank, while the ligand PDB file was converted to mol2 using OpenBabel and Antechamber. A frcmod file was supplied with the custom azobenzene parameters. 5.3.5

Benchmark

To evaluate the performance of FEP and docking, we performed a benchmark of 10 ligands with experimental PARP1 binding data. IC50 values for human PARP1 were retrieved for 1,542 ligands from the ExCape database 83,84 . We

𝛼′

N N

𝜔

𝛼

rotation

inversion

N N

cis

trans N N

N N

inversion-assisted rotation

concerted inversion

Supplementary Figure S12: Possible mechanisms of thermal isomerization. The inversion TS has α or α′ ≈ 180◦ , while the rotational TS has ω ≈ 90◦ . Inversion-assisted rotation combines inversion and rotation. Concerted inversion combines α inversion and α′ inversion.

converted IC50 values to binding free energies through ∆Gbind = kB T · log Kd , where we approximated the binding dissociation constant as Kd ≈ IC50. We selected the ten ligands with binding free energies that evenly interpolated between the highest and lowest values of ∆G. All selected ligands were neutral. Docking was performed with Dockstring, using the same setup as in azobenzene screening. FEP was performed as described above. Dockstring (de)protonates each ligand with OpenBabel, using a default pH of 7.4. However, this brings two complications when comparing to experiment. The first is that different experiments use different pHs, so care must be taken to choose the right pH for each data point. The second is that the ligand protonation state can change upon binding. In particular, for a ligand that is protonated at a given pH, the binding free energy is approximately the lower of the following two options: 1. The binding free energy of the protonated ligand. 2. The binding free energy of the unprotonated ligand, plus the free energy cost of deprotonation. The latter is equal to (kB T log 10) · (pKa − pH), which follows from the Henderson-Hasselbalch equation. Therefore, for the most accurate results, one should compute the binding free energy for both protonation states. The lower of (protonated binding free energy, unprotonated binding free energy + deprotonation cost) should then be compared to experiment. We used this approach in FEP for the two ligands in the benchmark that were protonated by OpenBabel. We also tried using it for docking, but found that performance was unchanged. For docking we therefore only used the protonated forms generated by OpenBabel. For azobenzene screening, we only ran FEP with the protonated forms from OpenBabel to save computational time. FEP and docking scores were adjusted by a constant shift to match the mean of the experimental data (using a linear regression did not improve results for either). Results are shown in CSI Fig. S11 and CSI Table S6. The docking results are fairly good, which is consistent with the performance reported in Ref. 30 . The Spearman rank correlation coefficient, which is the most important metric for ordering hits, is above 0.5. The Pearson correlation coefficient is above 0.6, while the R2 value is only 0.32. The FEP results are better, with a Spearman rank correlation of 0.81, R2 of 0.69, and R of 0.86. The RMSE is under 1.2 kcal/mol.

Chemical properties

6.1

Thermal half-lives

6.1.1

Mechanisms

The possible mechanisms for thermal azobenzene isomerization are given in CSI Fig. S12. Previously we argued that isomerization occurs through partial rotation to ω ≈ 70◦ , followed by intersystem crossing to a lower-energy triplet state, and then re-crossing to the singlet state at ω ≈ 105◦ 17 . We also showed that the isomerization rate from singlet rotation is highly correlated with that of singlet-triplet-singlet rotation. Therefore, to avoid the extra cost of calculating triplet gradients to train a triplet model, we simply computed the barrier of singlet rotation during screening. Since rotation can occur in the clockwise or counter-clockwise direction, this led to two mechanisms per species.

Supplementary Table S7: Training parameters used for all models. LR is the initial learning rate, LRmin is the minimum learning rate, and patience is the number of epochs without improvement in the validation loss before dropping the learning rate. Factor is the factor by which the learning rate is dropped. Train, validation, and test refer to the proportion of the total dataset used for the corresponding split.

All models

Optimizer Adam 85

LR 1e−4

LRmin 5e−6

Patience 50

Factor 2

Train 94%

Validation 3%

Test 3%

Supplementary Table S8: Loss coefficients for different models. All models use an MAE loss, with different coefficients for different quantities. Ei is the energy of state i, F⃗i are the forces in state i, and µ ⃗ ij is the transition dipole moment between states i and j. Training was performed in units of kcal/mol for energies, kcal/mol/Å for forces, and Debye for transition dipoles.

Thermal barrier Absorption spectrum Quantum yield

E0 0.4 – 0.4

F⃗0 1 – 1.0

E1 – – 0.4

F⃗1 – – 1.0

E1 − E 0 – 1.0 1.0

F⃗1 − F⃗0 – 0.2 1.0

µ ⃗ 01 – 4.0 –

Supplementary Table S9: TS model performance for 158 species outside the training set. We report the accuracy of ∆E † = ETS − Ecis , and of the forces F⃗ for both TS and cis geometries. Only rotational TSs are considered. Units are kcal/mol for energies and kcal/mol/Å for forces. “Remove ≥ 5 kcal/mol errors” means that we remove all species for which the prediction error in ∆E exceeds 5 kcal/mol. R2 is not provided for forces, because, by construction, all predicted forces are zero at critical points. Base compound Azobenzene Azonium All

Dataset All data Remove ≥ 5 kcal/mol errors All data Remove ≥ 5 kcal/mol errors All data Remove ≥ 5 kcal/mol errors

Data points 57 57 101 93 158 150

∆E † MAE (↓) 0.77 0.77 1.92 1.20 1.51 1.04

∆E † R2 (↑) 0.88 0.88 0.59 0.90 0.73 0.93

F⃗ MAE (↓) 0.61 0.61 1.13 1.12 0.94 0.93

While inversion might be favored over singlet- or triplet-mediated rotation for certain derivatives, preliminary work suggests that this is rare 17 . Moreover, affordable quantum chemical methods are not accurate enough to determine when inversion is favored 17 . Nevertheless, we also generated inversion TSs during training (but not screening), so that the model could be fine-tuned with more accurate methods in future work. 6.1.2

Data generation and training

Data was generated with the active learning loop shown in CSI Fig. S5. In each round we trained a committee of three NFFs, and then generated geometries for reactants, products, and transition states (TSs) for 1,000 cis-trans reactions. The isomer pairs were randomly sampled from our virtual library, which included both the azo and azonium forms of each compound. The TSs were generated with one of the three NFFs using the workflow in CSI Fig. S5(b). The workflow involves NN conformer generation for the reactants and products, together with a relaxed scan, conformer generation, and eigenvector following (EVF) to optimize TSs. Activation free energies ∆G† were computed as GTS − Gcis . Each free energy was computed with a modified rigid rotor-harmonic oscillator (mRRHO) approximation 86 applied to Hessian vibrational frequencies, plus a conformational entropy correction 87 . Further details can be found in Ref. 17 ; our only change here was to flip all negative frequencies above −30 cm−1 , and to discard equilibrium geometries with frequencies below −30 cm−1 , rather than zero. Geometries from the TS workflow were then sampled for quantum chemistry calculations. We selected 10,000 geometries in each round (15,000 in the final one), from a random sample of 100,000 configurations. As in Ref. 17 , 50% of the final geometries were selected based on uncertainty, 30% by high energy, and 20% randomly. Uncertainty was estimated as the standard deviation of the forces computed by the three models. These calculations were added to the training set, and the models were re-trained. This was repeated five times in total. Removing outliers and spin-contaminated geometries as in Ref. 17 yielded 47,313 data points in total. Combining these with 43,147 calculations of Ref. 17 led to a final training set of 90,460 geometries. The models were based on the PaiNN architecture 63 , with an analytical D3-BJ dispersion correction 56,57 added 13

to the output energies. They were trained with the settings in CSI Tables S7 and S8. We first pre-trained the models with 680,736 gas-phase SF-TDDFT 39 calculations from Ref. 16 . These calculations were performed for non-adiabatic molecular dynamics (NAMD) simulations. 95% of the data was used for training, 4% was used for validation, and 1% for testing. The models were then trained with SF-TDDFT calculations using an implicit solvent model of water, without freezing any parameters or changing the learning rate. We used SF-TDDFT with the BHHLYP functional 54 , the 6-31G* basis 55 , D3-BJ dispersion 56,57 , and a C-PCM description of water 58–60 using Q-Chem 5.3 40 . 6.1.3

MRSF-TDDFT refinement

We chose SF-TDDFT because it can account for some multi-reference effects 88 , which are important for the rotation mechanism in azobenzene 17 . However, we found that spin contamination was often severe for the rotational TSs of azonium compounds. About 50% of these azonium TSs had square spin ⟨S 2 ⟩ > 1. Moreover, azonium NN-EVF calculations only converged 49% of the time, while azo calculations converged 98% of the time. To fix this problem, we refined the final models with MRSF-TDDFT calculations 52,53 . MRSF-TDDFT is a spincomplete variant of SF-TDDFT, so it does not suffer from spin contamination. While the method was developed several years ago, it only recently became publicly available, and so we had not used it previously. We ran MRSFTDDFT on 17,027 geometries using GAMESS 41 (version release Sept. 30, 2022), with the same basis, functional, solvation model, and dispersion correction as in SF-TDDFT. The models were refined using the same learning parameters in CSI Table S7, and without freezing any parameters. EVF with this final model converged for 99.3% of azo species and 86.5% of azonium species. 6.1.4

Model accuracy

To evaluate the model’s performance, we used it to generate equilibrium and TS geometries for 158 species that were not in the training set. If a cis azobenzene derivative was in the training set, then neither the trans form nor any azonium forms could be in the evaluation set. Since a maximum of 5 TSs per species were generated from EVF (see Ref. 17 ), we randomly selected one TS per species to receive quantum chemical calculations. We performed MRSF-TDDFT calculations to evaluate the predicted activation energy, ∆E † , and the forces on each geometry. The results are given in CSI Table S9. The performance is excellent for azobenzene derivatives: the MAEs of ∆E † and F⃗ are 0.77 kcal/mol and 0.61 kcal/mol/Å, respectively. The R2 value for the activation energies is 0.88, which is also quite good. The performance is much worse for azonium derivatives, with MAEs of 1.92 kcal/mol and 1.13 kcal/mol/Å for energies and forces, respectively. Interestingly, however, the high azonium error for ∆E † mainly comes from a small number of compounds. Indeed, 8% of all azonium species have errors in ∆E † that exceed 5 kcal/mol. Removing these compounds leads to a 38% reduction in error, with the energy MAE dropping from 1.92 to 1.20 kcal/mol. This new MAE is rather good when one considers the spread in azonium activation energies. As shown in CSI Fig. S14(a), the distribution of ∆G† is much wider for azonium than for azobenzene (standard deviations of 4.30 and 2.95 kcal/mol, respectively). This leads to an increase in R2 , since the R2 metric measures the error relative to the width of the underlying distribution. Hence this new azonium subset has a higher R2 than azobenzene (0.90 vs. 0.88), despite also having a higher MAE. We therefore conclude that the predictions are quite good for 92% of azonium species and quite bad for 8%. Further, combining this 92% subset with all azobenzene species leads to a very good R2 value of 0.93. We inspected the eight species with high prediction errors, but could not find any obvious trends in the molecules or geometries. There were no common functional groups or abnormal configurations. However, we did find that the model systematically overestimated the barriers, with a mean signed error of 7.45 kcal/mol. We also found that re-training the model with the outlier data did not remove all outlier predictions. Therefore, it appears that there is an effect that is quite difficult to predict. 6.1.5

Quantum chemical validation

Hits were validated through single-point energy calculations with MRSF-TDDFT. For each species we used the TS geometry and cis geometry with the lowest NFF free energy. The resulting DFT energies were combined with the NFF-computed mRRHO and entropic corrections to obtain the free energy.

6.2

Protonation

6.2.1

Theory

We computed the pKa of azonium formation for both cis and trans isomers, as well as the free energy barrier to azonium isomerization. Trans azonium is of interest because it is red-shifted relative to azobenzene. Cis azonium

Supplementary Figure S13: Accuracy of pKa predictions and distribution of pKa values during virtual screening. (a) Predicted vs. experimental pKa values. (b) Distribution of computed pKa values during virtual screening. The pKa was only computed for compounds with hydrogen bond acceptors at ortho positions.

is of interest because it is required to compute the free energy of azonium-mediated isomerization. To compute this activation barrier, we require both the azonium barrier and the cis azonium pKa : (4) ∆G† = Gazonium TS + G(H2 O)n−1 (OH− ) − Gazo cis + G(H2 O)n = ∆G†azonium + ∆Gcis protonation

(5)

= ∆G†azonium + kB T log[1 + 10pH−cis pKa ].

(6)

The first bracketed term in the first line is the free energy of the azonium TS with n − 1 water molecules and one hydroxide. Subtracted from this term is the free energy of the cis isomer with n water molecules. The second term in the last line follows from the Henderson-Hasselbalch equation. 6.2.2

pKa implementation

The pKa for the m-fold protonation of a species A is given by 89 ∆G kB T log10 ∆G = Gproducts − Greactants pKa = −

(7)

reactants = (H2 O)n + A

(9)

products = (H2 O)n−m (OH− )m + A(H+ )m ,

(10)

(8)

where n is the number of water molecules. One can compute ∆G in two steps. First, one can compute the free energy of n water molecules in a box, followed by the free energy of n − m water molecules and m hydroxide ions. One can then compute the free energy of the protonated and unprotonated forms of A. Supplying the free energy difference to Eq. (7) then yields the pKa . However, quantum chemical free energy differences are typically not accurate enough to reproduce experimental data. Therefore, Eq. (7) is usually passed through a linear regression that is fit to experiment 89 . Computing free energies with double-hybrid DFT applied to geometries from semi-empirical DFT, and fitting to experiment with a linear regression, has yielded R2 values above 0.93 89 . In this work the calculation is further simplified because m = 1 for all compounds. This means that the change in free energy of water is constant for all species. Hence it can be absorbed into the intercept during linear regression, and does not need to be calculated. We simply calculated ∆G = Gazonium − Gazo and applied a linear regression to experiment to get the pKa . The free energies of the azo and azonium species were computed with our NN. We used the mRRHO approximation 86 applied to Hessian vibrational frequencies, plus a conformational entropy correction 87 , as described in Ref. 17 . Free energy differences were fit to nine experimental pKa measurements of azobenzene derivatives. Five came from Ref. 5 and four from Ref. 18 . We excluded compound 20 from Ref. 5 , since it was predicted by OpenBabel to have protonated substituents. We did not want this to cause problems in the pKa fit if the OpenBabel prediction was incorrect. We also removed the glutathione groups attached to the ends of the functional groups in Ref. 18 . This is

Supplementary Figure S14: Distribution of computed azonium properties. These calculations were only performed for azonium compounds with ortho hydrogen bond acceptors. (a) Comparison of the activation free energies of azo and azonium compounds. (b) Activation free energy vs. cis pKa of azonium species.

because glutathione is very flexible, which can lead to unconverged free energy estimates. It is also not representative of molecules in the training set. Both of these issues could have led to inaccuracies in the ML pKa predictions. The pKa predictions are compared to experiment in CSI Fig. S13(a). The results are quite good. R2 is greater than 0.8, the Spearman rank correlation coefficient is above 0.9, and the MAE is only 0.62. Note that, while there is a fairly even distribution of pKa values from 3 to 9 in the experimental data, this is unusual for the compounds screened here. Indeed, as shown in CSI Fig. S13(b), only 13% of compounds screened had a pKa value above 3. 6.2.3

Barrier implementation

Azonium activation free energies were computed as in Ref. 17 , but with one notable change. In Ref. 17 , we generated rotational TSs by performing a relaxed scan that forced the CNNC dihedral to ±90◦ . We also forced the CNN angles to 122◦ , so that we would not generate an inversion TS. However, when we did this for azonium derivatives, we found that only 10% of EVF jobs converged. In successful optimizations, we found that the azonium nitrogen formed a planar improper dihedral with its neighbors. In cis and trans equilibrium geometries, this improper dihedral was non-planar. We therefore added a third constraint to the scans, forcing the N[NH+ ]C improper dihedral to 180◦ . Further, we held the azonium hydrogen atom fixed during TS conformer generation, in addition to the CNNC atoms held fixed in Ref. 17 . The total activation free energy of proton-mediated isomerization was computed with Eq. (6). 6.2.4

Virtual screening

We computed azonium barriers and pKa values for all species with hydrogen bond acceptors at the ortho positions, or one atom away. We also did so for any species predicted to have pKa > 3 by the graph-to-property model. To compute the pKa for a given species, we first checked that free energies had been calculated for all azonium forms (i.e. the N[NH+ ] and [NH+ ]N forms for asymmetrically substituted compounds, which are equivalent for symmetric substitutions). We also checked that all vibrational frequencies were above −30 cm−1 , and discarded the species if they were not. We flipped the signs of negative frequencies above this value. We also removed any geometries with broken bonds, as measured through the D3 coordination number described in Ref. 17 . The azonium free energy was approximated as the lower of the free energies of the two azonium isomers. Compounds with simulated trans pKa > 5 and thermal half-lives in the target range were flagged as potential hits. 6.2.5

Quantum chemical validation

Hits were validated through single-point energy calculations with MRSF-TDDFT. For each species we used the geometry with the lowest NFF free energy. The resulting DFT energies were combined with the NFF-computed mRRHO and entropic corrections to obtain the free energy. This was in turn used to compute the pKa . For some species, the EVF calculation failed to converge for one of the two azonium forms. Typically EVF failed for the azonium form with the higher pKa . To estimate the free energy cost of this mechanism [Eq. (6)], we used the pKa of this azonium isomer with the barrier of the other isomer (i.e., the one for which EVF succeeded) in CSI Table S1.

To handle the failed EVF calculations, we also tried replacing NFF conformer generation with CREST 72 using GFN2-xTB 90 . As in NFF conformer generation, we fixed all CN[NH+ ]C atoms. We inspected the TS geometries generated by NFF and CREST, and found that both looked qualitatively wrong. EVF with the NFF failed for all TS guesses from CREST. Therefore, it seems that our general approach to TS generation caused the failure, not the NFF itself.

Photophysical properties

7.1

Absorption spectrum

7.1.1

Simulation

We 14 and others 50 have shown that it is important to consider the full absorption spectrum of azobenzenes, not simply the absorption wavelength of an optimized geometry. Azobenzene derivatives that absorb in the near-IR have absorption peaks in the visible 5 ; it is the high-wavelength tail of their spectra that allows for near-IR absorption. Moreover, dynamical effects can lead to asymmetric absorption spectra with long tails that could not be anticipated by single-point calculations 14 . We therefore simulated the full absorption spectrum of the azo and azonium forms of each trans compound. To do so we performed ground-state MD, and computed the absorption spectrum I as the weighted histogram of absorption frequencies at each frame 91 : I(ω) ∝ histogram({fi ωi }i∼MD ), 2

fi ∝ ωi |⃗ µi | .

(11) (12)

Here ω is the frequency, ωi is the absorption frequency of frame i sampled from MD, and fi is its oscillator strength. The latter is proportional to the absorption frequency and the square of the transition dipole moment µ ⃗ i . We discarded absorption energies outside of (1.3, 6.2) eV, and created histograms with 30 bins. Frequency spectra were then converted to wavelength spectra. This was done by first converting frequencies to wavelengths, and then multiplying by the Jacobian of the transformation, which is proportional to ω 2 . MD was initialized with conformers from NN conformer generation 17 . We performed 5 ps of MD per conformer and used 20 conformers per molecule, yielding 100 ps in total. Conformer i was sampled from the conformer ensemble with probability pi ∝ exp(−Ei /kB T ), where E is the energy. We performed MD with the Nosé-Hoover thermostat 73,74 and our trained NFF using the same MD parameters as in Ref. 17 . 7.1.2

Model

Absorption energies and transition dipole moment vectors were predicted with a single model. The model used the same PaiNN architecture as the ground state model. However, instead of outputting a single energy, it outputted the transition energy and a vector representing the transition dipole moment. The D3 correction was not applied to the output energy, since the D3 corrections of the ground- and excited-state energies cancel. The vector was produced by applying equivariant operations to the vector features ⃗vi generated in the message-passing stage, as described in Ref. 63 . 7.1.3

Training data

Training data was produced with TDDFT using the CAM-B3LYP functional 61 , the def2-SVP basis 62 , a C-PCM model of water 58–60 , and D3 dispersion 56 . As described below, we trained on transition dipole moments, transition energies, and transition energy gradients. Since the gradient of the transition energy is the difference of two energy gradients, we computed both S0 and S1 energy gradients with CAM-B3LYP. Along with the PBE0 92 functional, CAM-B3LYP is among the most reliable functionals for TDDFT vertical excitation energies of organic molecules 93 . Previously we used the PBE0 functional, and found that it reproduced the absorption spectrum of azobenzene quite well 14 . However, we used CAM-B3LYP here because, as a rangeseparated hybrid, it is presumably more accurate at describing excitations with charge-transfer character. It is therefore more likely to be reliable across the space of azobenzene derivatives. We have also found that it reproduces the absorption spectrum of azobenzene well (see Methods). Training geometries were sampled randomly from metadynamics simulations of equilibrium cis and trans configurations. These were performed during NN conformer generation in active learning. To increase the number of red-shifted geometries in the training set, we first trained an initial model on these random geometries, and then used it to predict the absorption energies of geometries without TDDFT labels. We then selected the geometries

Supplementary Table S10: Performance of the excited state absorption spectrum model for species outside the training set. λ is the absorption wavelength and µ ⃗ is the transition dipole vector. µ ⃗ was multiplied by ±1 for each geometry to minimize the error relative to DFT. Units are nm for λ and Debye for µ ⃗. Base compound Azobenzene Azonium All

Data points 170 329 499

λ MAE (↓) 3.51 7.03 5.83

λ R2 (↑) 0.98 0.97 0.98

µ ⃗ MAE (↓) 0.05 0.16 0.12

µ ⃗ R2 (↑) 0.98 0.98 0.98

with the highest predicted absorption wavelengths and performed TDDFT on them. We used 83,481 geometries in total. 7.1.4

Training

Parameters used for training are given in CSI Tables S7 and S8. We trained on both transition energies and their gradients to improve transition energy predictions. Previously we found that training on transition energies alone led to disappointing performance for red-shifted species 17 . One training complication is that quantum chemical transition dipole vectors have arbitrary signs. This is because the wavefunction of each electronic state has an arbitrary sign. These cancel for properties of a given state, since matrix elements are taken with respect to the same state, but do not cancel for transition properties. Therefore, as we did previously for non-adiabatic couplings 16 , we multiplied the model output by ±1 for each geometry in each training batch. The sign was chosen separately for each geometry to minimize the prediction error for that configuration. 7.1.5

Model performance

CSI Table S10 shows the model accuracy for 50 species outside of the training set. 10 geometries were sampled from NN-driven metadynamics for each compound. The performance is excellent for both azobenzene and azonium derivatives. The error in the absorption wavelength is 3.5 nm for azobenzenes and 7.0 nm for azonium species. Both are significantly lower than typical errors in TDDFT itself 94 . The R2 value is above 0.97 for both the absorption wavelength and the transition dipole moment for all groups of compounds. 7.1.6

Quantum chemical validation

Hits were validated with TDDFT by computing the absorption wavelength of the lowest energy conformer. The difference between the TDDFT result and the model result was used to shift the NFF-computed spectrum.

7.2

Quantum yield

Certain substituents can inhibit photoisomerization, e.g. through hydrogen bond-induced locking 95 or other mechanisms 96 . We therefore performed NN excited-state dynamics on our top candidates to predict whether they would isomerize under light. Previously we trained an excited-state NFF across the chemical space of azobenzene derivatives 16 . Combining it with surface-hopping MD 97 , we showed that it can predict quantum yields that are in moderate agreement with experiment. We developed an architecture based on diabatic electronic states 98 , and used it to improve the accuracy of NN fitting and dynamics for unseen species. Here we applied a similar workflow to our top candidates. We started with our pre-trained model from Ref. 16 , which was trained on over 600,000 calculations from over 8,000 species. We then refined it with quantum chemistry data from the top molecules. We performed active learning similar to CSI Fig. S5. We ran non-adiabatic excited state MD with the NFF, and then chose new geometries for quantum chemistry calculations either randomly or based on uncertainty, the S0 /S1 energy gap, or excited state barriers. Further details on these sampling techniques can be found in Ref. 16 . We began with active learning for unsubstituted azobenzene while we waited for screening results. We then focused on our top 20 candidates, performing calculations on geometries from equilibrium and TS metadynamics for each of the hits. We then sampled 250 geometries per species in each round of active learning, and removed non-hits from the training set as FEP results came back. The final dataset had 6,985 geometries. Analytic gradients are implemented for MRSF-TDDFT, but analytic non-adiabatic coupling vectors (NACVs) are not. We therefore could not use our diabatic architecture, since it requires NACVs for training. Instead we used a simple adiabatic model that predicts ground- and excited-state energies as separate outputs. We used the same dispersion-corrected PaiNN architecture as for thermal barriers, but with two energy outputs instead of one. We generated training data using MRSF-TDDFT with the 6-31G* basis 99 , the BHHLYP functional 54 , and Grimme’s 18

Supplementary Table S11: Performance of the quantum yield model. Geometries were sampled from NAMD performed with the final trained model. Results are divided by the method used to select samples. “ZN” refers to Zhu-Nakamura sampling to select low-gap geometries, “barrier” refers to sampling of excited-state barriers, and “uncertainty” refers to sampling based on model uncertainty. Details of these sampling methods can be found in Ref. 16 . Energies are given in kcal/mol, and forces are given in kcal/mol/Å. Sampled by ZN

Data points 150

Barrier

149

Random

214

Uncertainty

208

Metric MAE (↓) R2 (↑) MAE (↓) R2 (↑) MAE (↓) R2 (↑) MAE (↓) R2 (↑)

E0 0.57 1.00 0.63 1.00 0.60 1.00 1.01 1.00

E1 0.80 0.99 0.69 1.00 1.46 0.99 2.94 0.98

∆E01 0.71 0.90 0.62 0.99 1.59 0.95 2.92 0.91

F⃗0 1.55 0.98 1.02 1.00 1.17 0.99 2.65 0.79

F⃗1 1.60 0.97 1.08 0.99 2.00 0.95 7.04 −0.48

Supplementary Table S12: Performance of the graph-to-property models after each round of exploration and exploitation. We report the mean and standard deviation of each metric from three models. Units are kcal/mol for ∆G† and nm for λ. Iteration 1 2

∆G† MAE (↓) 2.4 ± 0.2 2.3 ± 0.1

∆G† R2 (↑) 0.54 ± 0.05 0.58 ± 0.04

λ MAE (↓) 21 ± 1 16 ± 1

λ R2 (↑) 0.60 ± 0.04 0.74 ± 0.03

pKa MAE (↓) 1.38 ± 0.03 1.49 ± 0.07

pKa R2 (↑) 0.62 ± 0.03 0.70 ± 0.03

Supplementary Table S13: Combined size of the training, validation, and test set for graph-to-property models in each iteration. We used a train/validation/test split of 0.9/0.05/0.05. Iteration 1 2

∆G† 4,855 9,207

λ 6,109 11,825

pKa 4,500 9,947

D3-BJ dispersion correction 56,57 . We used a C-PCM model for water 58–60 . Details of the training process are given in CSI Tables S7 and S8. Performance statistics are given in CSI Table S11. Surface hopping was performed with the Zhu-Nakamura method 42 , which requires only energies and gradients, not NACVs. Initial geometries were sampled from 15 ps of MD at 298.15 K in the NVT ensemble. 200 surface hopping trajectories were then performed for 4 ps each. Hops were restricted to geometries with S1 /S0 gaps below 0.5 eV. The quantum yield was computed as the proportion of trajectories that isomerized, excluding those that stayed in the excited state. Further details of our implementation can be found in Ref. 16 , and the code is available at https://github.com/learningmatter-mit/NeuralForceField.

Graph-to-property models

We trained graph-to-property models to predict isomerization barriers, absorption wavelengths, and pKa values generated by our NFF. We used the Attentive FP architecture for each model 43 . This architecture combines an attention mechanism 100 with neural network message-passing 101 to generate an embedding for each atom. These node embeddings are then converted to graph embeddings with attention-based pooling, before being mapped to the final property with a neural network. This model substantially improves on previous state-of-the-art results for nearly all properties in standard benchmarks 43 . We trained an ensemble of three models for each property. Barrier models were trained on activation free energies of both azo and azonium compounds. Absorption models were trained on the maximum absorption wavelength at 5% of the peak intensity. We converted these wavelengths to frequencies in eV before training, which we found to significantly improve model performance. The pKa models were trained on per-azonium pKa values. That is, they were trained to map the azonium SMILES to the pKa associated with that azonium isomer. There were usually two azonium isomers per azo compound, for the two nitrogen protonation sites, unless they were equivalent by symmetry. Models were trained with an initial learning rate of 10−3 for a maximum of 1,000 epochs. The learning rate was

dropped by a factor of two if the validation loss had not improved in 10 epochs. Training was stopped if the learning rate fell below 10−5 . The best model was chosen as the one with the lowest validation loss. We used a batch size of 16 and train/validation/test splits of 0.9/0.05/0.05. Models were created and trained using the DeepChem library 44 . Model performance is shown in CSI Table S12. The models are reasonably accurate, with R2 values between 0.54 and 0.74 depending on the property and generation. The R2 value increased for each model going from iteration 1 to 2. MAEs were around 2 kcal/mol for ∆G† , 15 to 20 nm for λ, and 1.5 pKa units for pKa . Dataset sizes are given in CSI Table S13.

References (1) Lerch, M. M.; Hansen, M. J.; van Dam, G. M.; Szymanski, W.; Feringa, B. L. Emerging targets in photopharmacology. Angewandte Chemie International Edition 2016, 55, 10978–10999. (2) Broichhagen, J.; Frank, J. A.; Trauner, D. A roadmap to success in photopharmacology. Accounts of chemical research 2015, 48, 1947–1960. (3) Szymanski, W.; Ourailidou, M. E.; Velema, W. A.; Dekker, F. J.; Feringa, B. L. Light-controlled histone deacetylase (HDAC) inhibitors: towards photopharmacological chemotherapy. Chemistry–A European Journal 2015, 21, 16517–16524. (4) Bonardi, F.; London, G.; Nouwen, N.; Feringa, B. L.; Driessen, A. J. Light-Induced Control of Protein Translocation by the SecYEG Complex. Angewandte Chemie International Edition 2010, 49, 7234–7238. (5) Dong, M.; Babalhavaeji, A.; Collins, C. V.; Jarrah, K.; Sadovski, O.; Dai, Q.; Woolley, G. A. Near-infrared photoswitching of azobenzenes under physiological conditions. Journal of the American Chemical Society 2017, 139, 13483–13486. (6) Knie, C.; Utecht, M.; Zhao, F.; Kulla, H.; Kovalenko, S.; Brouwer, A. M.; Saalfrank, P.; Hecht, S.; Bléger, D. ortho-Fluoroazobenzenes: Visible light switches with very long-lived Z isomers. Chemistry–A European Journal 2014, 20, 16492–16501. (7) Lipinski, C. A.; Lombardo, F.; Dominy, B. W.; Feeney, P. J. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Advanced drug delivery reviews 2012, 64, 4–17. (8) Bickerton, G. R.; Paolini, G. V.; Besnard, J.; Muresan, S.; Hopkins, A. L. Quantifying the chemical beauty of drugs. Nature chemistry 2012, 4, 90–98. (9) Capuzzi, S. J.; Muratov, E. N.; Tropsha, A. Phantom PAINS: Problems with the Utility of Alerts for P an-A ssay IN terference Compound S. Journal of chemical information and modeling 2017, 57, 417–427. (10) Senger, M. R.; Fraga, C. A.; Dantas, R. F.; Silva Jr, F. P. Filtering promiscuous compounds in early drug discovery: is it a good idea? Drug Discovery Today 2016, 21, 868–872. (11) Polykovskiy, D.; Zhebrak, A.; Sanchez-Lengeling, B.; Golovanov, S.; Tatanov, O.; Belyaev, S.; Kurbanov, R.; Artamonov, A.; Aladinskiy, V.; Veselov, M.; others Molecular sets (MOSES): a benchmarking platform for molecular generation models. Frontiers in pharmacology 2020, 11, 565644. (12) Ertl, P.; Schuffenhauer, A. Estimation of synthetic accessibility score of drug-like molecules based on molecular complexity and fragment contributions. Journal of cheminformatics 2009, 1, 1–11. (13) Coley, C. W.; Rogers, L.; Green, W. H.; Jensen, K. F. SCScore: synthetic complexity learned from a reaction corpus. Journal of chemical information and modeling 2018, 58, 252–261. (14) Axelrod, S.; Shakhnovich, E.; Gómez-Bombarelli, R. Mapping the space of photoswitchable ligands and photodruggable proteins with computational modeling. arXiv preprint arXiv:2302.11490 2023, (15) Trott, O.; Olson, A. J. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. Journal of computational chemistry 2010, 31, 455–461.

(16) Axelrod, S.; Shakhnovich, E.; Gómez-Bombarelli, R. Excited state non-adiabatic dynamics of large photoswitchable molecules using a chemically transferable machine learning potential. Nature Communications 2022, 13, 3440. (17) Axelrod, S.; Shakhnovich, E.; Gómez-Bombarelli, R. Thermal half-lives of azobenzene derivatives: Virtual screening based on intersystem crossing using a machine learning potential. ACS Central Science 2023, 9, 166–176. (18) Dong, M.; Babalhavaeji, A.; Hansen, M.; Kalman, L.; Woolley, G. Red, far-red, and near infrared photoswitches based on azonium ions. Chemical Communications 2015, 51, 12981–12984. (19) Aleotti, F.; Nenov, A.; Salvigni, L.; Bonfanti, M.; El-Tahawy, M. M.; Giunchi, A.; Gentile, M.; Spallacci, C.; Ventimiglia, A.; Cirillo, G.; Montali, L.; Scurti, S.; Garavelli, M.; Conti, I. Spectral Tuning and Photoisomerization Efficiency in Push–Pull Azobenzenes: Designing Principles. The Journal of Physical Chemistry A 2020, 124, 9513–9523, DOI: 10.1021/acs.jpca.0c08672. (20) Eyring, H. The activated complex and the absolute rate of chemical reactions. Chemical Reviews 1935, 17, 65–77. (21) Dokic, J.; Gothe, M.; Wirth, J.; Peters, M. V.; Schwarz, J.; Hecht, S.; Saalfrank, P. Quantum chemical investigation of thermal cis-to-trans isomerization of azobenzene derivatives: substituent effects, solvent effects, and comparison to experimental data. The Journal of Physical Chemistry A 2009, 113, 6763–6773. (22) Garcia-Amorós, J.; Sánchez-Ferrer, A.; Massad, W. A.; Nonell, S.; Velasco, D. Kinetic Study of the Fast Thermal Cis-to-Trans Isomerisation of Para-, Ortho- and Polyhydroxyazobenzenes. Physical Chemistry Chemical Physics 2010, 12, 13238–13242, DOI: 10.1039/C004340K. (23) Dunn, N. J.; Humphries, W. H., IV; Offenbacher, A. R.; King, T. L.; Gray, J. A. pH-Dependent Cis → Trans Isomerization Rates for Azobenzene Dyes in Aqueous Solution. The Journal of Physical Chemistry A 2009, 113, 13144–13151, DOI: 10.1021/jp903102u. (24) Ewert, J.; Heintze, L.; Jordà-Redondo, M.; von Glasenapp, J.-S.; Nonell, S.; Bucher, G.; Peifer, C.; Herges, R. Photoswitchable Diazocine-Based Estrogen Receptor Agonists: Stabilization of the Active Form inside the Receptor. Journal of the American Chemical Society 2022, 144, 15059–15071, DOI: 10.1021/jacs.2c03649. (25) Reinfelds, M.; Hermanns, V.; Halbritter, T.; Wachtveitl, J.; Braun, M.; Slanina, T.; Heckel, A. A Robust, Broadly Absorbing Fulgide Derivative as a Universal Chemical Actinometer for the UV to NIR Region. ChemPhotoChem 2019, 3, 441–449, DOI: 10.1002/cptc.201900010. (26) Bortolus, P.; Monti, S. Cis-trans photoisomerization of azobenzene. Solvent and triplet donors effects. Journal of Physical Chemistry 1979, 83, 648–652. (27) Cheng, Y.-C.; Prusoff, W. H. Relationship between the inhibition constant (Ki) and the concentration of inhibitor which causes 50 per cent inhibition (I50) of an enzymatic reaction. Biochemical Pharmacology 1973, 22, 3099–3108. (28) Fragment Libraries. https://enamine.net/compound-libraries/fragment-libraries, Accessed: 2023-0114. (29) RDKit: Open-source cheminformatics. http://www.rdkit.org. (30) Garcı́a-Ortegón, M.; Simm, G. N.; Tripp, A. J.; Hernández-Lobato, J. M.; Bender, A.; Bacallado, S. DOCKSTRING: easy molecular docking yields better benchmarks for ligand design. Journal of chemical information and modeling 2022, 62, 3486–3502. (31) YANK. https://github.com/choderalab/yank, Accessed: 2023-03-02. (32) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. Journal of chemical theory and computation 2015, 11, 3696–3713.

(33) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. The Journal of chemical physics 2004, 120, 9665–9678. (34) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. Journal of computational chemistry 2004, 25, 1157–1174. (35) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1, 19–25. (36) Chodera, J. D. A simple method for automated equilibration detection in molecular simulations. Journal of chemical theory and computation 2016, 12, 1799–1805. (37) Chodera, J. D.; Shirts, M. R. Replica exchange and expanded ensemble simulations as Gibbs sampling: Simple improvements for enhanced mixing. The Journal of chemical physics 2011, 135, 194110. (38) Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. The Journal of chemical physics 2008, 129, 124105. (39) Shao, Y.; Head-Gordon, M.; Krylov, A. I. The spin–flip approach within time-dependent density functional theory: Theory and applications to diradicals. The Journal of chemical physics 2003, 118, 4807–4818. (40) Epifanovsky, E.; Gilbert, A. T.; Feng, X.; Lee, J.; Mao, Y.; Mardirossian, N.; Pokhilko, P.; White, A. F.; Coons, M. P.; Dempwolff, A. L.; others Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package. The Journal of chemical physics 2021, 155, 084801. (41) Barca, G. M. J. et al. Recent developments in the general atomic and molecular electronic structure system. The Journal of Chemical Physics 2020, 152, 154102, DOI: 10.1063/5.0005188. (42) Yu, L.; Xu, C.; Lei, Y.; Zhu, C.; Wen, Z. Trajectory-based nonadiabatic molecular dynamics without calculating nonadiabatic coupling in the avoided crossing case: Trans⇀ ↽cis photoisomerization in azobenzene. Physical Chemistry Chemical Physics 2014, 16, 25883–25895. (43) Xiong, Z.; Wang, D.; Liu, X.; Zhong, F.; Wan, X.; Li, X.; Li, Z.; Luo, X.; Chen, K.; Jiang, H.; others Pushing the boundaries of molecular representation for drug discovery with the graph attention mechanism. Journal of medicinal chemistry 2019, 63, 8749–8760. (44) Ramsundar, B.; Eastman, P.; Walters, P.; Pande, V.; Leswing, K.; Wu, Z. Deep Learning for the Life Sciences; O’Reilly Media, 2019; https://www.amazon.com/Deep-Learning-Life-Sciences-Microscopy/dp/ 1492039837. (45) Šolı́nová, V.; Kašička, V. Determination of acidity constants and ionic mobilities of polyprotic peptide hormones by CZE. Electrophoresis 2013, 34, 2655–2665, DOI: 10.1002/elps.201300119. (46) Mukadum, F.; Nguyen, Q.; Adrion, D. M.; Appleby, G.; Chen, R.; Dang, H.; Chang, R.; Garnett, R.; Lopez, S. A. Efficient discovery of visible light-activated azoarene photoswitches with long half-lives using active search. Journal of Chemical Information and Modeling 2021, 61, 5524–5534. (47) Cheron, N.; Jasty, N.; Shakhnovich, E. I. OpenGrowth: an automated and rational algorithm for finding new protein ligands. Journal of medicinal chemistry 2016, 59, 4171–4188. (48) Misal, S. A.; Gawai, K. R. Azoreductase: a key player of xenobiotic metabolism. Bioresources and Bioprocessing 2018, 5, 1–9. (49) Brown, M. A.; De Vito, S. C. Predicting azo dye toxicity. Critical reviews in environmental science and technology 1993, 23, 249–324. (50) Gelabert, R.; Moreno, M.; Lluch, J. M. Predicting the Electronic Absorption Band Shape of Azobenzene Photoswitches. International Journal of Molecular Sciences 2023, 24, 25. (51) Talaty, E. R.; Fargo, J. C. Thermal cis–trans-isomerization of substituted azobenzenes: a correction of the literature. Chemical Communications (London) 1967, 65–66. 22

(52) Lee, S.; Filatov, M.; Lee, S.; Choi, C. H. Eliminating spin-contamination of spin-flip time dependent density functional theory within linear response formalism by the use of zeroth-order mixed-reference (MR) reduced density matrix. The Journal of chemical physics 2018, 149, 104101. (53) Lee, S.; Kim, E. E.; Nakata, H.; Lee, S.; Choi, C. H. Efficient implementations of analytic energy gradient for mixed-reference spin-flip time-dependent density functional theory (MRSF-TDDFT). The Journal of Chemical Physics 2019, 150, 184111. (54) Becke, A. D. A new mixing of Hartree–Fock and local density-functional theories. The Journal of chemical physics 1993, 98, 1372–1377. (55) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. The Journal of Chemical Physics 1982, 77, 3654–3665. (56) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. The Journal of chemical physics 2010, 132, 154104. (57) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. Journal of computational chemistry 2011, 32, 1456–1465. (58) Truong, T. N.; Stefanovich, E. V. A new method for incorporating solvent effect into the classical, ab initio molecular orbital and density functional theory frameworks for arbitrary shape cavity. Chemical Physics Letters 1995, 240, 253–260. (59) Barone, V.; Cossi, M. Quantum calculation of molecular energies and energy gradients in solution by a conductor solvent model. The Journal of Physical Chemistry A 1998, 102, 1995–2001. (60) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Energies, structures, and electronic properties of molecules in solution with the C-PCM solvation model. Journal of computational chemistry 2003, 24, 669–681. (61) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange–correlation functional using the Coulombattenuating method (CAM-B3LYP). Chemical physics letters 2004, 393, 51–57. (62) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Physical Chemistry Chemical Physics 2005, 7, 3297–3305. (63) Schütt, K.; Unke, O.; Gastegger, M. Equivariant message passing for the prediction of tensorial properties and molecular spectra. International Conference on Machine Learning. 2021; pp 9377–9388. (64) O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An open chemical toolbox. Journal of cheminformatics 2011, 3, 1–14. (65) Mysinger, M. M.; Carchia, M.; Irwin, J. J.; Shoichet, B. K. Directory of useful decoys, enhanced (DUD-E): better ligands and decoys for better benchmarking. Journal of medicinal chemistry 2012, 55, 6582–6594. (66) Chodera, J. D.; Mobley, D. L.; Shirts, M. R.; Dixon, R. W.; Branson, K.; Pande, V. S. Alchemical free energy methods for drug discovery: progress and challenges. Current opinion in structural biology 2011, 21, 150–160. (67) Eastman, P.; Friedrichs, M. S.; Chodera, J. D.; Radmer, R. J.; Bruns, C. M.; Ku, J. P.; Beauchamp, K. A.; Lane, T. J.; Wang, L.-P.; Shukla, D.; others OpenMM 4: a reusable, extensible, hardware independent library for high performance molecular simulation. Journal of chemical theory and computation 2013, 9, 461–469. (68) Leimkuhler, B.; Matthews, C. Rational construction of stochastic numerical methods for molecular sampling. Applied Mathematics Research eXpress 2013, 2013, 34–56. (69) Pham, T. T.; Shirts, M. R. Identifying low variance pathways for free energy calculations of molecular transformations in solution phase. The Journal of chemical physics 2011, 135, 034114. (70) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. The Journal of chemical physics 1995, 103, 8577–8593. 23

(71) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method. Journal of computational chemistry 2000, 21, 132–146. (72) Pracht, P.; Bohle, F.; Grimme, S. Automated exploration of the low-energy chemical space with fast quantum chemical methods. Physical Chemistry Chemical Physics 2020, 22, 7169–7192. (73) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. The Journal of chemical physics 1984, 81, 511–519. (74) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Physical review A 1985, 31, 1695. (75) Larsen, A. H. et al. The atomic simulation environment—a Python library for working with atoms. Journal of Physics: Condensed Matter 2017, 29, 273002. (76) PDBFixer. https://github.com/openmm/pdbfixer, Accessed: 2023-04-19. (77) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. The Journal of chemical physics 1983, 79, 926–935. (78) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Antechamber: an accessory software package for molecular mechanical calculations. J. Am. Chem. Soc 2001, 222 . (79) Sousa da Silva, A. W.; Vranken, W. F. ACPYPE-Antechamber python parser interface. BMC research notes 2012, 5, 1–8. (80) Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. LINCS: A linear constraint solver for molecular simulations. Journal of computational chemistry 1997, 18, 1463–1472. (81) Hess, B. P-LINCS: A parallel linear constraint solver for molecular simulation. Journal of chemical theory and computation 2008, 4, 116–122. (82) Berendsen, H. J.; Postma, J. v.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. The Journal of chemical physics 1984, 81, 3684–3690. (83) Sun, J.; Jeliazkova, N.; Chupakhin, V.; Golib-Dzib, J.-F.; Engkvist, O.; Carlsson, L.; Wegner, J.; Ceulemans, H.; Georgiev, I.; Jeliazkov, V.; others ExCAPE-DB: an integrated large scale dataset facilitating Big Data analysis in chemogenomics. Journal of cheminformatics 2017, 9, 1–9. (84) ExCAPE-DB: ExCAPE chemogenomics database. https://solr.ideaconsult.net/search/excape/, Accessed: 2023-02-07. (85) Kingma, D. P.; Ba, J. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 2014, (86) Grimme, S. Supramolecular binding thermodynamics by dispersion-corrected density functional theory. Chemistry–A European Journal 2012, 18, 9955–9964. (87) Chang, C.-e. A.; Chen, W.; Gilson, M. K. Ligand configurational entropy and protein binding. Proceedings of the National Academy of Sciences 2007, 104, 1534–1539. (88) Lee, S.; Shostak, S.; Filatov, M.; Choi, C. H. Conical intersections in organic molecules: benchmarking mixedreference spin–flip time-dependent DFT (MRSF-TD-DFT) vs spin–flip TD-DFT. The Journal of Physical Chemistry A 2019, 123, 6455–6462. (89) Pracht, P.; Wilcken, R.; Udvarhelyi, A.; Rodde, S.; Grimme, S. High accuracy quantum-chemistry-based calculation and blind prediction of macroscopic pKa values in the context of the SAMPL6 challenge. Journal of Computer-Aided Molecular Design 2018, 32, 1139–1149. (90) Bannwarth, C.; Ehlert, S.; Grimme, S. GFN2-xTB–An accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions. Journal of chemical theory and computation 2019, 15, 1652–1671.

(91) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Electronic absorption spectra and solvatochromic shifts by the vertical excitation model: Solvated clusters and molecular dynamics sampling. The Journal of Physical Chemistry B 2015, 119, 958–967. (92) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. The Journal of chemical physics 1999, 110, 6158–6170. (93) Jacquemin, D.; Wathelet, V.; Perpete, E. A.; Adamo, C. Extensive TD-DFT benchmark: singlet-excited states of organic molecules. Journal of Chemical Theory and Computation 2009, 5, 2420–2435. (94) Greenman, K. P.; Green, W. H.; Gómez-Bombarelli, R. Multi-fidelity prediction of molecular optical peaks with deep learning. Chemical science 2022, 13, 1152–1162. (95) Bandara, H. D.; Friss, T. R.; Enriquez, M. M.; Isley, W.; Incarvito, C.; Frank, H. A.; Gascon, J.; Burdette, S. C. Proof for the concerted inversion mechanism in the trans→cis isomerization of azobenzene using hydrogen bonding to induce isomer locking. The Journal of organic chemistry 2010, 75, 4817–4827. (96) Bandara, H. D.; Cawley, S.; Gascón, J. A.; Burdette, S. C. Short-circuiting azobenzene photoisomerization with electron-donating substituents and reactivating the photochemistry with chemical modification. Chemistry – A European Journal 2011, 2011, 2916–2919. (97) Tully, J. C. Molecular dynamics with electronic transitions. The Journal of Chemical Physics 1990, 93, 1061–1071. (98) Mead, C. A.; Truhlar, D. G. Conditions for the definition of a strictly diabatic electronic basis for molecular systems. The Journal of Chemical Physics 1982, 77, 6090–6098. (99) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. Ab initio Molecular Orbital Theory. 1986. (100) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, L.; Polosukhin, I. Attention is all you need. Advances in neural information processing systems. 2017; pp 5998–6008. (101) Gilmer, J.; Schoenholz, S. S.; Riley, P. F.; Vinyals, O.; Dahl, G. E. Neural message passing for quantum chemistry. International conference on machine learning. 2017; pp 1263–1272.

Record · ID 138922 · SHA-256 8f879ee92ae68fe6
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.