Predictive modeling of visible-light azo-photoswitches’ properties using structural features

Abstract

In this manuscript we present the strategy for modeling photoswitch properties (maximum absorption wavelength and thermal half-life of photoisomers) of visible-light azo-photoswitches using structural data. We compile a comprehensive data set from literature sources and perform a rigorous benchmark to select the best feature type and modeling approach. The fragment counts have demonstrated the best performance in the benchmark for both properties. We validate the models in cross-validation and on an external set. The predictions of absorption wavelengths for this set are highly accurate; on the other hand, the model for thermal half-life is less reliable, likely due to the modest size of the data set related to half-life of photoisomers, although consensus modeling approach allows to improve the predictivity. We also provide an interpretation of the modeling results using ColorAtom approach and the insights into the chemical space covered by the data set.

Scientific contribution The paper provides a machine learning approach based only on structural features to predict two important photoswitch properties. Unlike previous studies, we do not use any quantum chemical features which accelerates the modeling procedure, while the accuracy of models remains high. Moreover, the fragment counts offer unique approach to model interpretation that is useful for rational design of photoswitches with desired properties.

Introduction

Photoswitches are a category of compounds whose chemical structures and properties can be changed by light irradiation [1]. The photochemically generated species (photoisomers) are less stable and hence can be reverted to more stable species thermally or via a reverse photoisomerization. In the field of photoswitches, the wavelength with the maximum light absorption (λ_max) and the thermal half-life of a metastable photoisomer (t_1/2) decide their use in various applications in material science and biology [2]. For instance, the well-known azobenzene photoswitch has λ_max of 320 nm and hence requires UV light for light-induced changes (trans to cis photoisomerization). However, the UV light is harmful to the living beings as it damages biological cells and its penetration ability to the deep tissue is very low compared to green and yellow lights in the visible region. Thus, for biological applications such as photopharmacology, where the drug action is controlled by the light, visible-light responsive photoswitches are required [3]. Moreover, photoswitches having long λ_max and long t_1/2 of meta-stable photoisomers are useful for various applications. However, some studies [4, 5] show that introducing particular substitution patterns may lead to inverse proportionality of λ_max and t_1/2, thus, the design of visible-light photoswitches would require careful consideration of structural patterns to avoid this effect. So far, experimentalists heavily rely on the time-consuming density functional theory calculations to estimate the λ_max and t_1/2 properties of photoswitches [6,7,8]. We envisioned to develop a machine learning model to accurately predict the λ_max and t_1/2 properties focusing on the visible-light photoswitches.

Previously constructed machine learning (ML) models were focused on the broad category of azo-dyes with the prediction of single parameter such as λ_max [9]. Another important parameter of photoswitches is the photoisomerization quantum yield, which can be predicted by a neural network machine learning model [10]. Very recently, ML-based model also developed for the prediction of thermal half-lives of azobenzene derivatives, though the authors used the data obtained in quantum mechanical calculations as a training set [11]. In the category of photoswitches, prediction of multiple parameters such as λ_max at two electronic excitations (π → π* transition of trans and cis isomers) has been investigated using a multi-output Gaussian process model and molecular fingerprints (FP) [12]. Such structural descriptors derived directly from 2D representations of molecules allow for a significant acceleration of the modeling process. Indeed, the authors report to have achieved comparable predictive accuracy to time-dependent density functional theory (TD-DFT) calculation with reduced inference time on a curated data set consisting of the λ_max of trans and cis isomers of various photoswitches.

However, to the best of our knowledge, accurate prediction of two different parameters such as λ_max and t_1/2 properties in a standalone machine learning model based on structural features has not been achieved. In this paper, we present a computational study on visible-light photoswitches, employing quantitative structure–property relationship (QSPR) modeling to predict two crucial properties, λ_max and t_1/2 (Fig. 1). We compile a comprehensive data set from literature sources and perform a rigorous benchmark study in order to select the structural descriptors that would be the most appropriate to predict these properties. The accuracy of obtained ML models is also validated externally on a pre-selected test set. We also discuss the interpretability and the applicability domains of these models, as well as the challenges the endeavors of such modeling face.

Fig. 1

General workflow of the study. Literature data set of azobenzenes and azoheteroarenes is used for the benchmark of models built on a variety of molecular fingerprints and fragment count descriptors. The best performing model is selected via cross-validation and validated on an external test set

Data and methods

Data collection and curation

The basis for the predictive models in this work is literature data. The data on the absorption properties is the most abundant; on the other hand, there are no large publicly available datasets that contain curated experimental data on the isomerization of azo-photoswitches. Following entries make part of the integrated data set: (1) 191 azobenzene dyes reported in [13, 14]; (2) 212 various heteroaryl dyes with spectra measured in ethanol, reported in [9]; (3) 91 heteroaryl dyes with spectra measured in acetonitrile, reported in [15,16,17,18,19,20,21]; (4) 80 molecules with both λ_max and half-life time measured, reported in [5, 22, 23]; (5) photoswitch data set reported in [12], containing azo-photoswitches with a variety of experimentally measured properties, including absorption maxima and isomerization rate. All data in the sources is annotated by room temperature or explicitly 25 ℃. Moreover, while most data sources explicitly indicate the stable isomer, some entries lack this information. For modeling purposes, we assume that the E-isomer is the more stable form for all azo-photoswitches in the dataset. The structures where this is not explicitly stated are labeled as “E*”. More detailed description of sources and data is available in SI.

It should be noted that not all properties are reported for each compound in these sets. We have performed an additional curation to merge all available data into a single data set. For merging, chemical structures were compared between all sets. In the case of duplicates, the most recent or the most complete (containing more reported properties) data was kept. In total, there are 798 unique compounds with measured λ_max, and 134 with measured t_1/2. We have divided our data based on their properties (wavelength and half-life) into two sets: 90% for cross-validation (CV), hyperparameter optimization and training (718 and 120 data points for λ_max and t_1/2, respectively), and 10% for external validation (80 and 14 data points for λ_max and t_1/2, respectively). The separation of data sets was performed randomly, with stratification by property values, so that the full range of values would be present in both training and test sets (see Fig. 2 for the property distribution in each set). The full data set is available in Additional file 2.

Fig. 2

The histograms of the distribution of λ_max and logt_1/2 in the training and test sets

Model hyperparameter optimization and modeling

To build a robust QSPR model, compounds need to be represented as relevant numerical parameters—molecular descriptors. It is impossible to know in advance which descriptors will provide the best predictivity. Most often, the choice of the best descriptors is done via a rigorous benchmark study. While prediction of photoswitch properties, especially the isomerization rate and half-life, would intuitively benefit from the knowledge of 3D structures, the data that is used in this work does not reflect any conformational information. Therefore, we only consider several types of 2D structural descriptors, including CircuS and ChyLine fragments [24], Morgan [25], Atom pairs [26], Avalon [27], and native RDKit FP [28]. Molecular FP are binary features (they only take values 0 or 1) commonly used in QSPR modeling, especially for biological activities [29]. Morgan, or extended circular, fingerprints encode a structure by presence/absence of substructures of circular topology (atom and neigboring environment of a certain radius). It is also possible to generalize such representation further by using chemical features (aromaticity, H-bond donor/acceptors, etc.) instead of atom symbols, herein called Morgan features. Native RDKit FP similarly encode the substructures of linear and branched topologies. Like Morgan FP, aditional degree of generalization may be introduced by “layering” the substructures with different information (topology, aromaticity, bond orders, ring sizes, etc.), further referenced as “RDKit FP layered”. Atom pairs generally contain the least information: only pairs of atoms with distances (number of bonds) between them are recorded. Avalon FP account for a variety of structural information, although this information is unavailable externally via the RDKit library. Fragment descriptors, like Circular Substructures (CircuS) and Chython Linear (ChyLine) fragments that we use here, improve on the information content of the features by recording the actual count of substructures that each descriptor represents. We have previously shown that CircuS descriptors outperform FP in other tasks [24]. Details on the parameters for each type of descriptors are available in Additional file 1. All models also include a preprocessing step that removes features with zero variance and scales the features to the range of 0 to 1.

The machine learning algorithm of choice here was Support Vector Machines (SVM) [30], as it is capable of handling smaller data sets with high precision due to high customizability of the fitting function (for results of other methods in the benchmarks, see Additional file 1). However, such adaptability also requires rigorous optimization of method hyperparameters, including the descriptor space. To tackle this, we employ the in-house optimization library—DOPtools [31]—that uses well-established optimization strategies to identify the optimal hyperparameters for the ML model. The construction and evaluation of these models adheres to best practices in cheminformatics [32]. A repeated 10-fold CV strategy was carefully applied during hyperparameter optimization, serving as a robust means of assessing the model's predictive performance. The performance is statistically evaluated by Root Mean Square Error (RMSE) and determination coefficient R² on the test set predictions:

$$RMSE= \frac{1}{N}\sum_{N}{\left({y}_{obs,i}-{y}_{pred,i}\right)}^{2}$$

$${R}^{2}= 1-\frac{\sum_{N}{\left({y}_{obs,i}-{y}_{pred,i}\right)}^{2}}{\sum_{N}{\left({y}_{obs,i}-{\widehat{y}}_{obs,i}\right)}^{2}}$$

where N is the number of points in the set, y_obs,i is the experimentally observed value of the ith data point, y_pred,i is the predicted value of the ith data point, and ŷ_obs,i is the average observed value across the set.

The analysis of chemical space is performed using UMAP method [33] in Python 3.10. The interpretation of the model is done using ColorAtom method [34] as implemented in the DOPtools library. ColorAtom calculates the weights of fragments as partial derivatives of the model predictions. The atoms accumulate the weights of the fragments they participate in, which eventually are summed up to calculate the atomic contribution. The visualization of these atomic contributions allows to see which parts of the molecule have the highest contribution to the prediction. The code for the analysis of benchmark and modeling results is available in GitHub repository.

Results and discussion

Prediction of λ _max

The model for the prediction of λ_max was based on the ensemble of data with experimentally measured spectra, 798 compounds in total. The benchmark study on molecular descriptors was performed along with the hyperparameter optimization, so that every type of descriptors achieves the best possible predictivity in CV. The results can be seen in the Fig. 3. Here, the CircuS and ChyLine fragments show the best performance. Both types are based on the fragment enumeration, with the difference in the substructure topologies. CircuS fragments count the number of occurrences of circular fragments, and ChyLine—of linear fragments, improving the information content over molecular fingerprints. The accuracy of predictions here is on par with the state-of-the-art approaches reported elsewhere (e.g., Griffiths et al. [12], which is the largest model to-date, report RMSE of 20.9 nm vs 21.6 nm in cross-validation for our model. Note that the scores are not directly comparable, as they report random train/test splits instead of CV).

Fig. 3

Descriptor benchmark results (left) for the model predicting λ_max. Each boxplot represents the distribution of scores (RMSE, in nm) for 10 repeats of CV on the training set with random shuffling (white square for the mean score, the box for the interquartile range IQR, whiskers for 1.5 IQR, other points are outliers). The best model was obtained on CircuS fragments, achieving R² = 0.907 in CV and R²(test) = 0.905 on an external test set, the observed vs predicted plot (right) shows the CV and external results for this model. Error bar on the points indicated the standard deviation of the prediction across CV repeats

Note that this model does not consider the physico-chemical parameters of the solvents, although the solvents may have a considerable effect on the absorption of photoswitches [35]. Since not all data points in the data set presented here are annotated by the solvent in which the measurements are performed, building a model with solvent parameters included would reduce the size of the data set significantly. Still, we have performed such modeling, and the results are very similar to the model above (R² = 0.907, RMSE = 21.7 nm in CV); the details are available in Additional file 1.

Prediction of t _1/2

We use the half-life (in s) measurements in logarithmic scale in predictive models due to the extremely wide range of values in a non-logarithmic scale. The photoswitch data set provides the values of the isomerization rate, which were transformed into t_1/2 according to the first-order kinetic equation. The total number of data points for the model training was 134. Modeling followed the same protocol as described above. The benchmark was performed on models built only using structural parameters of the compounds, and the results are shown in Fig. 4 below. As the figure shows, the best model was built on linear ChyLine fragments, although the performance is overall quite poor (R² = 0.584 in CV and 0.482 in external validation).

Fig. 4

Descriptor benchmark results (left) for the model predicting logt_1/2. Each boxplot represents the distribution of scores (in log units of t_1/2) for 10 repeats of CV on the training set with random shuffling (the boxplot elements follow the same designation as in Fig. 2). The best model was obtained on ChyLine fragments, the observed vs predicted plot (right) shows the CV results for this model. Error bar on the points indicated the standard deviation of the prediction across CV repeats

One way to overcome weak predictions by an individual model is using consensus modeling [32], i.e., averaging the predictions made by several models. This allows to complement the initial model’s bad predictions by the insights of another model, thus, hopefully, reducing the errors of predictions for outliers. In this case, we have started by calculating consensus predictions of the best model (ChyLine) with the second best (RDkit layered FP) and have observed a great improvement of the results, especially for two notable outliers by the middle of the property range (Fig. 5, center). Adding a model built on RDkit linear FP further improved the predictions for these molecules (Fig. 5, on the right), however, a notable outlier by the lower end of the property had worse predictions here. It should be noted that this molecule is the only thiazole derivative in the test set, and all other thiazole-based azo-dyes in the training set have a much higher logt_1/2 value, indicating that this compound is most likely outside of the applicability domain of these models.

Fig. 5

Comparison of performance of the individual best model and its consensus with other models. Observed vs predicted are shown in each plot, with the vertical lines on the consensus plots showing improvement (in green) or deterioration (in red) of prediction compared to the individual model. Consensus models improve the predictions for most outliers of the test set, except the molecule with the lowest logt_1/2 value which is most likely outside of the applicability domain

Interpretation of predictions and analysis of chemical space

The validation of the model’s accuracy on an external data set is supplemented by the interpretation of the model and the analysis of the chemical spaces of training and test sets for each property. Note that this is performed for the individual best model, as the descriptor space or a single hyperparameter setup for the consensus model cannot be defined. We have studied the chemical spaces of training and test set for both properties with UMAP method by projecting the descriptor spaces selected in the benchmark by the best model into a 2-dimensional space (see Fig. 6, on the left for both sets). The Fig. 6A shows that the large and diverse data set related to λ_max covers the chemical space quite well, and the molecules of the test set are close to the clusters of the training set. On the other hand, the data set related to logt_1/2 is relatively small and, more importantly, less diverse. Most molecules of the test set are found outside of the initial clusters, thus, which increases the likelihood of them being outside of the applicability domain of the model.

Fig. 6

Interpretation of models for prediction of (a) λ_max and (b) logt_1/2 by ColorAtom (on the right) and visualization of chemical space of each model (on the left). In ColorAtom, the contributions of atoms are coded blue for negative contributions and red for positive, with the intensity of color indicating the scale of the effect. White-coded atoms have virtually no contribution to the prediction. The contributions are scaled to the maximum in the test set for each property, to allow the comparison of the effect; the scale is indicated as a colorbar. Visualization of the chemical space is made by UMAP on the descriptor shown to be the best in the benchmark; molecules are colored red for training set and blue for test set. IDs are given for some molecules used for model interpretation; they follow the indexing in the Additional file 2

The interpretation of the predictions here is done using ColorAtom methodology [34]. In ColorAtom, the weights of fragment descriptors are calculated as partial derivatives of prediction and define the atomic contributions of all atoms in a molecule. This allows to see how the presence of different functional groups affects the predicted values and guide the rational synthesis of compounds with desired properties. Generally speaking, the ColorAtom is applicable to one molecule at a time; however, here (Fig. 6) we have scaled the atomic contributions to the maximum in the subset, so that all molecules can be compared (full tables with the whole test set are available in SI, Additional files 1 and 3). It should be noted, however, that some limitations should be considered for ColorAtom. First, it is limited to fragments descriptors (CircuS, ChyLine, etc.) and would not work with fingerprints; luckily, the best models in both cases were based on fragments. Second, the algorithm that calculates the weights of descriptors assumes linear, or at least monotonous, relationship, as a partial derivative is used, so the interpretations should be used with caution. Finally, although it is true for any ML model, interpretations of the predictions are limited to the model and, thus, to the training data, meaning that if some effects are not accounted for in the data, they would not be relevant in the atomic contributions.

As far as the interpretations go here, we can observe some trends that are reported in literature. Functional groups (e.g., NO₂, CN, OMe) in their role of chromophore or auxochromes on the azo-dye affect λ_max and logt_1/2 of a photoisomer [17, 23]. For example, interpreting the model for λ_max, we observed high positive contribution of functional groups (e.g., NO₂ and NEt₂ in entries 541 and 317, Fig. 6A) for the prediction, which is consistent with previous reports [36, 37]. Similarly, introducing heteroarenes (entries 738, 805, Fig. 6A) into a dye has a positive, albeit smaller, effect [38, 39]. It should be noted that these effects are often associated with the push/pull character of interaction between electron-donating and withdrawing groups or other electronic effects; our models, however, only consider the structural descriptors, which may only implicitly contain such information. Similarly, logt_1/2 can be increased by the presence of heteroarenes [40] (entries 60, 68, Fig. 6B) or by substitution patterns (e.g., F in ortho position has been reported to increase half-life [38, 41], as in entry 50). Sometimes, however, ColorAtom may omit the expected effect from functional groups (e.g., OMe in entry 780, Fig. 6B), likely due to relatively small contribution compared to other fragments in a molecule.

An important limitation of the presented models is the omission of the solvent effects. While these are important for both λ_max and t_1/2, the currently collected dataset does not contain sufficient data to account for them. Specifically, some of the data sources for λ_max lack solvent annotations. On the other hand, while all data points in the set related to half-life have those, their distribution is highly uneven, with the majority (95 entries) measured in acetonitrile. We report the models following different strategies—either isolating the dataset of the majority solvent or using solvent properties as descriptors—in SI (see Additional file 1). Neveretheless, netiher approach led to a significant improvement of predictions. We believe that a more systematic dataset is required to assess solvent effects properly and encourage the community to contribute to its creation and collection.

Conclusions

In this manuscript we outline the modeling strategies for prediction of both λ_max and t_1/2 of photoisomers of visible-light azo-photoswitches. The structural features (molecular fingerprints and fragment counts) that we investigate provide a solid base for the accurate and fast modeling and virtual screening. Moreover, model interpretation by the ColorAtom approach allows for the rational design of new compounds, by pinpointing the structural motifs leading to the desired property values. The predictions of the model were validated on an external set of 80 and 14 for λ_max and t_1/2, respectivelly. The data for λ_max of azobenzene and azoheteroarene families of photoswitches is abundant, and the accuracy of predictive models, even built on purely 2D representations of molecules without any electronic information, is quite high, as we demonstrate here. The external test set of photoswitches has demonstrated relatively small errors of prediction of λ_max, despite the property not being directly correlated to the structural features we employ. On the other hand, the lack of data on the t_1/2, as well as the general complexity of the mechanism of isomerization in azo-photoswitches lead to not quite satisfactory performance of such models. While generally we observe quite good predictivity in the external set, several outliers with quite large prediction errors (up to 2 log units for t_1/2) are found. We would like to note that the chemical space covered by the data set for t_1/2 is quite narrow, which may lead to a limited reliability of predictions in the higher range of the property. Moreover, the application of consensus modeling allows to significantly reduce the error of prediction for most outliers. We provide the full integrated data set to encourage further investigations using other machine learning methods or features, including quantum chemical calculations, and the expansion of the chemical space of photoswitches. A more systematic data set for the half-life, incorporating a variety of both structural features and measurement conditions such as solvents, would benefit the future efforts in modeling this property with machine learning methods.