Data-Driven and Multiscale Modeling of DNA-Templated Dye Aggregates

Abstract

Dye aggregates are of interest for excitonic applications, including biomedical imaging, organic photovoltaics, and quantum information systems. Dyes with large transition dipole moments ($\mathsf{\mu }$) are necessary to optimize coupling within dye aggregates. Extinction coefficients ($\mathsf{\epsilon }$) can be used to determine the $\mathsf{\mu }$ of dyes, and so dyes with a large $\mathsf{\epsilon }$ (>150,000 M⁻¹cm⁻¹) should be engineered or identified. However, dye properties leading to a large $\mathsf{\epsilon }$ are not fully understood, and low-throughput methods of dye screening, such as experimental measurements or density functional theory (DFT) calculations, can be time-consuming. In order to screen large datasets of molecules for desirable properties (i.e., large $\mathsf{\epsilon }$ and $\mathsf{\mu }$), a computational workflow was established using machine learning (ML), DFT, time-dependent (TD-) DFT, and molecular dynamics (MD). ML models were developed through training and validation on a dataset of 8802 dyes using structural features. A Classifier was developed with an accuracy of 97% and a Regressor was constructed with an ${\mathrm{R}}^2$ of above 0.9, comparing between experiment and ML prediction. Using the Regressor, the $\mathsf{\epsilon }$ values of over 18,000 dyes were predicted. The top 100 dyes were further screened using DFT and TD-DFT to identify 15 dyes with a $\mathsf{\mu }$ relative to a reference dye, pentamethine indocyanine dye Cy5. Two benchmark MD simulations were performed on Cy5 and Cy5.5 dimers, and it was found that MD could accurately capture experimental results. The results of this study exhibit that our computational workflow for identifying dyes with a large $\mathsf{\mu }$ for excitonic applications is effective and can be used as a tool to develop new dyes for excitonic applications.

1. Introduction

Organic molecules, which absorb and emit light, also known as dyes, are useful for many applications, such as biomedical imaging [1,2], organic photovoltaics [3,4], non-linear optics [5], and quantum information systems [6,7,8,9]. Key parameters that determine the performance of the dyes in those applications include the extinction coefficient ($\mathsf{\epsilon }$) and transition dipole moment ($\mathsf{\mu }$), as well as aggregation ability. Thus, optimizing the key electronic (e.g., $\mathsf{\mu }$) and molecular (e.g., aggregate) features is crucial for the desired applications of dye molecules. The interaction of dyes with light can be quantified via their extinction coefficient, ε. The value of $\mathsf{\epsilon }$, resulting from the absorption of light by the dye, can be measured using optical spectroscopy. From the measured value of $\mathsf{\epsilon }$, the transition dipole moment μ can be extracted [10,11,12]. This relationship not only allows for the measurement of the μ of dyes, but also helps select the dye candidates with optimal electronic properties for excitonic applications (i.e., large $\mathsf{\mu }$) from numerous dyes. The value of $\mathsf{\epsilon }$ strongly depends on the molecular structure of the dye. Some efforts have been made to augment the $\mathsf{\epsilon }$ of dyes, such as adding donor or acceptor groups to the π-conjugation network [13], extending the π-conjugation network [14], and making the dye structures more planar [15,16,17]. However, the relationship between dye structure and $\mathsf{\epsilon }$ remains unclear. In addition, it is time-consuming to conduct either experimental measurements or computational modeling to screen many dye candidates for desirable properties (e.g., high $\mathsf{\epsilon }$ and μ).

Another key feature of dyes for excitonic applications is dye aggregation. Dye aggregation has been observed in natural systems [18,19] as well as artificial systems [20,21]. Dye aggregates feature exciton delocalization, which facilitates energy transfer through the aggregate [22]. The dynamics of excitons residing on a dye aggregate can be described using the Frenkel Hamiltonian [23], where the exchange of excitons is largely dependent on the transition dipole coupling of the dyes [22,23,24]. The dipole coupling strength, exciton delocalization, and corresponding dynamics depend on the electronic properties of individual dyes, or monomers, as well as the orientations of the dyes in the aggregate [22,23,24,25,26]. One method of facilitating dye aggregation in a controlled and predictable manner is using DNA scaffolds. Dyes attached to DNA scaffolds, such as duplexes, Holliday junctions, and origami, have been shown to aggregate into dimers, trimers, and tetramers [27,28,29,30,31,32,33,34,35,36,37,38,39]. Characterization of the optical properties of the aggregates reveals that the dyes can adopt various ideal orientations [22,24,26,40]. One orientation, called an H-aggregate, occurs when the dyes are stacked, and results in a blue-shifted absorption spectrum. When the dyes are oriented head-to-tail, called a J-aggregate, the absorption spectrum is red-shifted. Another aggregate, termed oblique, occurs when the dyes are at 90° to one another, and results in Davydov splitting of the absorption spectrum. To maximize or fine-tune the coupling between dye molecules, the aggregation should be predictable and controlled.

Computational studies of dyes aim to identify optimal candidates for excitonic applications. For example, density functional theory (DFT) and time-dependent (TD-) DFT can be used to screen the effects of functional groups on dye electronic properties [41,42,43,44]. Our prior studies indicated that functional group substitution can affect the solvation free energy ${\mathsf{\Delta }\mathrm{G}}_{\mathrm{solv}}$, transition dipole moment $\mathsf{\mu }$, and absorption wavelength ${\mathsf{\lambda }}_{\max }$ of a dye [45,46]. This effect is correlated with the empirically derived Hammett constant (${\mathsf{\sigma }}_{\mathrm{p}}$), which demonstrates the electron-donating or electron-withdrawing strength of a substituent. The DFT and TD-DFT methods are applicable to dye monomers, but not to dye aggregates attached to DNA scaffolds due to structural size and complexity. An alternative method to further screen dyes with favorable electronic properties attached to DNA scaffolds is molecular dynamics (MD), which has been used to study dye–DNA interactions [47,48,49,50]. In a recent study, Mathur et al. used MD to study the orientations of cyanine dyes attached to DNA bundles and found that they were able to accurately capture dye dynamics and orientations using MD [47]. Nicoli et al. utilized MD to study the aggregation of Cy3 dimers attached to DNA duplexes, and found that MD could accurately capture the stacking of dyes leading to H-aggregation [50]. However, both DFT and MD are time-consuming for high-throughput screening of dye candidates with desired properties. Recently, machine learning (ML) has been shown to be a viable method of screening thousands of molecules to identify structure–property relationships based on both computational (e.g., DFT) [51] and experimental data [52,53]. The problem of searching through chemical-based datasets that contain labeled data for optimal molecules is a common task for pharmaceuticals and dye-sensitized solar cells, but there has not yet been work specifically targeting optimization of dye candidates for dye aggregate–DNA constructs. In particular, our group is interested in near-IR dye molecules exhibiting a large $\mathsf{\epsilon }$ (>150,000 M⁻¹cm⁻¹) and hydrophobic properties. The same photophysical data used to create the chemical-based datasets for dyes have been of interest for organic photovoltaics, and there are several public datasets available [52,53,54,55]. ML techniques applied to chemical space exploration is a rich field with a variety of methods from which to choose. The methods can be hierarchical with the size of the dataset, spanning from well-established supervised learning to more complex artificial neural networks [56].

In this work, a systematic approach, combining ML, DFT, TD-DFT, and MD methods, was used to screen dye monomers from an expansive dataset and provide insight into dye aggregate–DNA duplex interactions. We first used ML to identify ideal dye candidates with high extinction coefficients (ε) from a dataset of around 18,000 molecules. Then, for the 100 ML-selected dye candidates with desirable structural features and high values of ε, DFT and TD-DFT calculations were performed to predict their ground and excited state properties. Finally, benchmark MD simulations were conducted to reveal the interactions between the selected dye dimers and the DNA duplexes.

2. Methods

2.1. Machine Learning

Classifier and Regressor models were trained to identify ideal dye candidates with high extinction coefficients (ε) based on dye structure features. The Classifier model could quickly classify the dyes with either high or low ε, where we set a threshold of 150,000 M⁻¹cm⁻¹ for strong exciton coupling in dye aggregates. As it learned from the Classifier model, the Regressor model could further estimate the values of ε for the dyes. Three data sources, including Deep4Chem [57], PhotoChem CAD 3 [54], and Dyomics GmbH [58] (8802 molecule datapoints in total), were used to train and validate the models. We utilized SMILES format for the molecule. We also utilized RDKit [59] to calculate 284 different features, such as the maximum carbon chain length and aromatic, amide, and ester group counts. The development dataset, containing 90% of data points, was used to train and test various model hyperparameters. Based on the accuracy for the Classifier and R² for the Regressor, a model was selected for further analysis. We then used the validation dataset, containing 10% of data points, to validate the selected model’s effectiveness. The molecules with ε of above 800,000 M⁻¹cm⁻¹ were excluded because the ε values deviated too greatly from the threshold of 150,000 M⁻¹cm⁻¹. Figure 1 shows the dataset breakdown with high and low ε values. The three data sources all have a data imbalance, where the number of molecules with low ε values is larger than that of molecules with high ε values. The effect of imbalanced data is discussed in Section 3.1.

2.2. Density Functional Theory

Density functional theory (DFT) and time-dependent (TD-) DFT calculations were performed to optimize dye structures in the ground state and calculate solvation energies and transition dipole moments. Similarly to our previous work [46], the dyes were optimized with the M06-2X [60] functional and 6-31+G(d,p) basis set, to a residual force of 4.5 × 10⁻⁴ Hartree/Bohr. The M06-2X functional with 6-31+G(d,p) basis set was validated in our prior studies of pristine and substituted cyanine and squaraine dyes [45,46] and has been used successfully for the calculations of the excited state properties of similar systems [61,62,63]. Frequency calculations were conducted to confirm the ground state structures were true minima. To determine $\mathsf{\mu }$, single point, vertical excited state calculations using the M06-2X functional were performed on the ground state structures to obtain transitions to the first 30 excited singlet states and identify the state with the largest oscillator strength. Calculations of the ground and excited state properties were conducted with implicit water solvation using the integral equation formalism polarizable continuum model (IEFPCM) [64,65], which was successfully used for the excited state property calculations of similar systems [66,67,68]. Excited state calculations were conducted assuming nonequilibrium solvent conditions.

To approximate the relative hydrophobicity of the dyes, the partitioning coefficient between n-octanol and water, $\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)$, was calculated according to [69,70]

$$\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)=-\frac{{\mathsf{\Delta }\mathrm{G}}_{\mathrm{o}}-{ \mathsf{\Delta }\mathrm{G}}_{\mathrm{w}}}{2.3\mathrm{RT}}$$

(1)

where ${\mathsf{\Delta }\mathrm{G}}_{\mathrm{o}}$ and ${\mathsf{\Delta }\mathrm{G}}_{\mathrm{w}}$ are the Gibbs free energy of solvation for a dye in n-octanol and water, respectively; $\mathrm{R} =8.31\frac{\mathrm{J}}{\mathrm{mol}·\mathrm{K}}$; and $\mathrm{T} =273.15 \mathrm{K}$. A more positive value of $\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)$ means a molecule is more hydrophobic, and a more negative value means a molecule is more hydrophilic. In general, the Gibbs free energy of solvation for a molecule (${\mathsf{\Delta }\mathrm{G}}_{\mathrm{solv}}$) is a measure of the amount of energy required to dissolve the dye in solvent, and was calculated according to [45,46,68,71]

$${\mathsf{\Delta }\mathrm{G}}_{\mathrm{solv}}={ \mathrm{E}}_{\mathrm{solvated}}-{ \mathrm{E}}_{\mathrm{vacuum}}$$

(2)

where ${\mathrm{E}}_{\mathrm{solvated}}$ is the total energy of the dye in implicit solvent and ${\mathrm{E}}_{\mathrm{vacuum}}$ is the total energy of the dye in vacuum. Calculations for the solvation energy were conducted using the universal solvation model based on density (SMD) variation of IEFPCM [72], which was useful for predicting the solvation energies of organic molecules [73] and calculating the relative hydrophobicity of modified squaraine dyes [70]. All DFT and TD-DFT calculations were conducted using the Gaussian16 software package [74].

2.3. Molecular Dynamics

Molecular dynamics (MD) simulations were performed with the GROMACS 2020.3 software package [75]. Dye–DNA structures were built using the UCSF ChimeraX software [76] with the dyes initialized on the outside of the DNA backbone. The OL15 force-field [77] with non-bonded modifications [78] was used for DNA parameters and the generalized amber forcefield (GAFF) [79] was used for dye parameters. Atomic charges for the dyes were calculated using the HF/6-31G* theory level [80]. The 26-basepair dsDNA duplex sequence and dye locations from Huff et al. [29] was used. The dye–DNA structures were solvated in TIP3P water [81] in a truncated octahedron box with 1.2 nm between the dye–DNA structure and the box edge. Mg²⁺ ions were used to neutralize the system. Cannon et al. and Huff et al. showed, experimentally, that by adding excess MgCl₂ to solutions containing DNA duplexes with two pentamethine indocyanine Cy5 dyes, DNA Holliday junctions could be formed [28,29]. Because of this, no excess MgCl₂ was used in the MD simulations apart from replacing a necessary number of water molecules with Mg²⁺ to achieve neutral charge. Neighbor-searching was used with a cutoff of 1.2 nm. Van der Waals interactions were limited to 1.2 nm, and the particle mesh Ewald (PME) was used with a real-space coulomb cutoff of 1.2 nm. Bonds to hydrogen atoms were constrained using the LINCS algorithm [82]. A timestep of 2 fs was used.

The initial systems were energy-minimized with the steepest descent method for 1000 steps. Then, to achieve a well-relaxed starting structure, two subsequent 10 ns equilibration steps were performed with harmonic constraints, the first with 1000 $\frac{\mathrm{kJ}}{\mathrm{mol}·{\mathrm{nm}}^2}$ spring constants applied to non-hydrogen atoms, and the second with 100 $\frac{\mathrm{kJ}}{\mathrm{mol}·{\mathrm{nm}}^2}$ spring constants applied to non-hydrogen atoms, keeping the number of atoms, volume, and temperature constant. A final 10 ns equilibration was performed with no restraints. Following equilibration, 1 $\mathsf{\mu }\mathrm{s}$ production simulations were carried out, keeping the number of atoms, pressure, and temperature constant. The velocity-rescale thermostat [83] was used to maintain a constant temperature of 300 K with a coupling time of 0.1 ps, with the DNA–dye and solvent being coupled separately. The Parrinello–Rahman barostat [84] was used to keep the pressure at 1 atm with a coupling time of 1.0 ps. Coordinates were written every 10 ps and the first 100 ns of the production simulations were treated as equilibration periods, and so were not used for analysis.

To determine the transition dipole coupling strength between two dyes, the orientation of the dyes with respect to one another and the transition dipole moments $\mathsf{\mu }$ are needed. The dye–dye center-to-center distances $({\mathrm{R}}_{\mathrm{m},\mathrm{n}})$ and orientation factors ($\mathsf{\kappa }$) were determined every 10 ps. The values of $\mathsf{\kappa }$ were determined using [48]

$$\mathsf{\kappa } ={\widehat{\mathsf{\mu }}}_{\mathrm{m}}·{\widehat{\mathsf{\mu }}}_{\mathrm{n}}-3\left({\widehat{\mathrm{R}}}_{\mathrm{m},\mathrm{n}}·{\widehat{\mathsf{\mu }}}_{\mathrm{m}}\right)\left({\widehat{\mathrm{R}}}_{\mathrm{m},\mathrm{n}}·{\widehat{\mathsf{\mu }}}_{\mathrm{n}}\right)$$

(3)

where ${\widehat{\mathsf{\mu }}}_{\mathrm{i}}$ is the transition dipole moment unit vector of dye $\mathrm{m}$ or $\mathrm{n}$ (taken along the long axis of the dye), and ${\widehat{\mathrm{R}}}_{\mathrm{m},\mathrm{n}}$ is the unit vector between the centers of dyes $\mathrm{m}$ and $\mathrm{n}$. When $\left|\mathsf{\kappa }\right|=0$ or $\left|\mathsf{\kappa }\right|=1.5$, the dyes are in a stacked oblique orientation or tail-to-tail oblique orientation, respectively. When $\left|\mathsf{\kappa }\right|=1$ or $\left|\mathsf{\kappa }\right|=2$, the dyes are in a stacked (H-aggregate) or head-to-tail (J-aggregate) orientation, respectively.

The exciton exchange energy (${\mathrm{J}}_{\mathrm{m},\mathrm{n}}$), which is a measure of the strength of the transition dipole coupling between two dyes, depends on the transition dipole moment $\mathsf{\mu }$, which is related to $\mathsf{\epsilon }$ [10,11,12] and can be obtained experimentally or using TD-DFT. The ${\mathrm{J}}_{\mathrm{m},\mathrm{n}}$ of a dimer was approximated using the extended dipole model according to [85]

$${\mathrm{J}}_{\mathrm{m},\mathrm{n}}={ \mathrm{J}}_0\left(\frac{1}{\left|{\mathrm{r}}_{\mathrm{m}}-{ \mathrm{r}}_{\mathrm{n}}\right|}-\frac{1}{\left|{\mathrm{r}}_{\mathrm{m}}-{ \mathrm{s}}_{\mathrm{n}}\right|}-\frac{1}{\left|{\mathrm{s}}_{\mathrm{m}}-{ \mathrm{r}}_{\mathrm{n}}\right|}+\frac{1}{\left|{\mathrm{s}}_{\mathrm{m}}-{ \mathrm{s}}_{\mathrm{n}}\right|}\right)$$

(4)

where the values ${\mathrm{r}}_{\mathrm{i}}$ and ${\mathrm{s}}_{\mathrm{i}}$ correspond to either end of dye $\mathrm{m}$ or $\mathrm{n}$ along the dye’s long axis. The pre-factor term, ${\mathrm{J}}_0$, is defined as [85]

$${\mathrm{J}}_0=\frac{{\mathsf{\mu }}_{\mathrm{m}}{\mathsf{\mu }}_{\mathrm{n}}}{4\mathsf{\pi }{\mathsf{ϵ}}_{\mathrm{o}}{\mathrm{n}}^2{\mathrm{l}}_{\mathrm{m}}{\mathrm{l}}_{\mathrm{n}}}$$

(5)

where ${\mathsf{\mu }}_{\mathrm{i}}$ is the transition dipole moment magnitude (calculated using TD-DFT) of dyes $\mathrm{m}$ or $\mathrm{n}$, ${\mathsf{ϵ}}_0$ is the vacuum permittivity constant, $\mathrm{n}$ is the refractive index of water (1.33), and ${\mathrm{l}}_{\mathrm{m}}$ and ${\mathrm{l}}_{\mathrm{n}}$ are the lengths of dyes $\mathrm{m}$ and $\mathrm{n}$, respectively (such that ${\mathrm{l}}_{\mathrm{i}}=\left|{\mathrm{r}}_{\mathrm{i}}-{ \mathrm{s}}_{\mathrm{i}}\right|$).

3. Results

3.1. Dye Screening Using Machine Learning and Density Functional Theory

A Random Forest Classifier was trained on the development dataset utilizing a five k-fold to determine the best ε threshold; the max feature, which refers to the maximum number of features to consider; the max depth, which refers to the maximum depth of the tree; and the class weight for the model. Every model we created has a high accuracy of 97% or above. Figure 2a shows the accuracy of the Random Forest Classifier with various ε thresholds in comparison with a model (labeled Always Low ε) that always classifies a molecule with low ε. The Always Low ε model has a high accuracy of around 92% or above, indicating that the high accuracy of 97% or above for the Random Forest Classifier is reasonable. In addition, the accuracy of the Random Forest Classifier starts to converge at 150,000 M⁻¹cm⁻¹, which supports our selection of the threshold of 150,000 M⁻¹cm⁻¹. To address a concern of the effect of imbalanced data on the accuracy of the models, we also trained and validated the models with different datasets that included various percentages of molecules with low ε values. Figure 2b shows the accuracy vs. percentage of low ε in the dataset. The value of 0.5 represents the balanced dataset containing 50% molecules with high ε values and 50% molecules with low ε values. The number of molecules with low ε values gradually increases, so the dataset becomes imbalanced. However, the accuracy of the Random Forest Classifier remains at 97% or above. These results indicate that the effect of the imbalanced data we used to train and validate the models is negligible.

Based on the performance of the Random Classifier model, we developed a Random Forest Regressor model to further predict the precise ε value of a molecule. We trained the Regressor to determine its best hyperparameters of class weight, max features, max depth, and criterion as we did with the Random Forest Classifier. Figure 3 compares the predicted and actual ε values for the development and validation datasets, which have 0.95 and 0.91 for R² (i.e., the coefficient of determination), respectively. The datasets and codes developed in this study will be published online.

To identify dyes with high $\mathsf{\epsilon }$ (which is indicative of high $\mathsf{\mu }$) the optimized Regressor model was further applied to a dataset containing around 18,000 potential dye candidates, including the commercially available cyanine dyes such as Cy3, Cy5, Cy5.5, and Cy7, shown in Figure 4. These dyes are known to exhibit large $\mathsf{\epsilon }$ (and thus, large $\mathsf{\mu }$ and strong dye coupling) [86], and are of interest to our research group and collaborators [27]. Four modified Cy5 dyes with hydrophobic substituents were also considered in this study. These dyes, labeled as Cy5-Cl, Cy5-Peg, Cy5-hex, and Cy5-tBu, are also shown in Figure 4 and are hypothesized to exhibit stronger dye coupling due to being more hydrophobic compared to Cy5, which may result in shorter inter-dye distances [87]. Two other modified dyes, Cy5-CN and Cy5-NMe₂, were considered since CN and NMe₂ were shown to have a large effect on excess dipole moments (the difference in the dipole moments of the ground and excited states) of similar dyes [45,46]. The rest of the dataset consisted of dyes obtained from PubChem [88], including dyes with similar structures to cyanine, porphyrin, and methyl violet molecules. Those three classes of dyes were chosen for their prominent $\mathsf{\pi }$-conjugation [89,90,91], absorption in the visible light range [90,92,93], and excitonic applications [90,92,94,95,96].

While $\mathsf{\mu }$ can be extracted from experimental absorption spectra [10,11,12], it is not possible to determine $\mathsf{\mu }$ based on the peak $\mathsf{\epsilon }$ alone. Because of this, TD-DFT was used to calculate $\mathsf{\mu }$, to compare with experimentally available $\mathsf{\epsilon }$. Figure 5 shows a comparison of ML-predicted $\mathsf{\epsilon }$ and TD-DFT-calculated $\mathsf{\mu }$ with commercially available Cy3, Cy5, Cy5.5, and Cy7 dyes, as well as modified Cy5 dyes for which experimentally measured values of $\mathsf{\epsilon }$ are available. The $\mathsf{\epsilon }$ values for the commercially available dyes were obtained from their respective commercial websites, including AAT Bioquest [97], Lumiprobe [98], Glen Research [99], and Interchim [100], and from Huff et al. [29]. Because multiple vendors advertise slightly different $\mathsf{\epsilon }$ values, a range is given in Figure 5. Meares et al. synthesized Cy5-hex, Cy5-Peg, Cy5-tBu, and Cy5-Cl for incorporation into DNA and measured the $\mathsf{\epsilon }$ of the dyes incorporated into DNA strands at their peak wavelengths [87]. The ranges of $\mathsf{\epsilon }$ for Cy5, Cy5-hex, Cy5-Peg, Cy5-tBu, and Cy5-Cl are likely due to small differences in local environments when the dyes are attached to DNA sequences and relative purities [87]. In general, ML-predicted $\mathsf{\epsilon }$ values agree with the trend of experimental $\mathsf{\epsilon }$ obtained from literature. Notably, ML-predicted $\mathsf{\epsilon }$ values for Cy5, Cy5.5, and Cy7 are within the experimental $\mathsf{\epsilon }$ range indicated by the shaded region in Figure 5. Similarly, the ML-predicted $\mathsf{\epsilon }$ trend agrees with the trend of TD-DFT-calculated $\mathsf{\mu }$ values, which are assumed to be correlated (i.e., larger $\mathsf{\epsilon }$ leads to larger $\mathsf{\mu }$). The dyes that do not fall into the range of experimental $\mathsf{\epsilon }$ consist of the Cy5 derivatives developed by Meares et al. [87]. The percent errors from the experiments include 45% for Cy5-hex, 49% for Cy5-Peg, 10% for Cy5-tBu, and 20% for Cy5-Cl. Such differences could be caused by solvent, DNA-dye interactions, and dye purities [87]. Furthermore, the specific functional groups for Cy5-hex and Cy5-Peg might not be well represented by the ML training dataset, which could lead to inaccuracies when predicting $\mathsf{\epsilon }$. In general, the Regressor is able to predict the overall trend of $\mathsf{\epsilon }$, which is necessary for the screening of numerous new dyes for $\mathsf{\epsilon }$.

The Regressor model was also applied to the dataset obtained from PubChem [88] to identify additional potential dye candidates. The top 100 dyes that were predicted to have $\mathsf{\epsilon }$ above 150,000 M⁻¹cm⁻¹ were then screened using DFT and TD-DFT to determine their $\mathsf{\mu }$ values by calculating vertical excited state transitions to the lowest 30 excited states. Of the 100 dye candidates, the 15 dye candidates with desirable properties, such as absorption wavelength in the visible region, large $\mathsf{\pi }$-conjugated networks, and $\mathsf{\mu }$ comparable to that of Cy5 (within 50%), are shown in Figure 5 and are labeled 1–15. Their corresponding ML-predicted $\mathsf{\epsilon }$ and TD-DFT-calculated $\mathsf{\mu }$ values are listed in

Table 1. ML-Regressor predicted extinction coefficient ($\mathsf{\epsilon }$) and TD-DFT-calculated transition dipole moments ($\mathsf{\mu }$) of dyes shown in Figure 4. Cartesian $\mathsf{\mu }$ vector components along x, y, and z axis are provided in Table S1.

Dye	ML-Predicted $\mathsf{\epsilon }$ $\times 1000 {\mathbf{M}}^{-1}\mathbf{c}{\mathbf{m}}^{-1}$	TD-DFT $\mathsf{\mu }$, Debye
Cy3	126	12.25
Cy5	222	15.35
Cy5-CN	138	15.66
Cy5-NMe2	216	15.83
Cy5-Cl	214	15.74
Cy5-hex	110	16.22
Cy5-Peg	107	16.19
Cy5-tBu	224	15.90
Cy5.5	205	15.57
Cy7	235	17.62
1	307	9.08
2	288	7.84
3	265	20.25
4	240	7.99
5	235	20.17
6	229	14.00
7	227	10.66
8	227	15.49
9	226	14.68
10	223	16.19
11	222	14.00
12	218	11.26
13	216	9.52
14	212	9.34
15	210	10.48

.

Comparing the TD-DFT calculated values of $\mathsf{\mu }$ for Cy3, Cy5.5, Cy7, and the Cy5 derivatives, Cy3 has the lowest $\mathsf{\mu }$. Disregarding Cy5-hex and Cy5-Peg, the values of $\mathsf{\mu }$ and the ML-predicted $\mathsf{\epsilon }$ generally follow the same trend, with Cy3 having the smallest $\mathsf{\mu }$ and Cy7 having the largest. Comparing the 15 selected dyes from the Regressor model, all dyes are predicted to have an $\mathsf{\epsilon }$ above 210,000 M⁻¹cm⁻¹. Dye 1 has the largest overall ML-predicted $\mathsf{\epsilon }$ of 309,000 M⁻¹cm⁻¹, with a TD-DFT $\mathsf{\mu }$ of 9.08 D. Dye 3 has the largest overall TD-DFT $\mathsf{\mu }$ of 20.25 D and an ML-predicted $\mathsf{\epsilon }$ of 265,000 M⁻¹cm⁻¹.

Figure 6 shows the log(P_o/w) values calculated with vacuum and implicit solvent DFT for the ML-selected dyes. A more positive value of log(P_o/w) indicates a more hydrophobic dye, and a more negative value of log(P_o/w) indicates a more hydrophilic dye. It is hypothesized that by increasing hydrophobicity, dyes may aggregate closer, thus improving coupling. This has been demonstrated in a set of squaraine dyes modified with hydrophobic substituents [70], and the values of log(P_o/w) for the hydrophobic squaraine dyes are similar to those for the Cy5 derivatives. Comparing the Cy5 derivatives, Cy3, Cy5.5, and Cy7, Cy5-hex is the most hydrophobic and Cy3 is the least hydrophobic, closely followed by Cy5-CN. Most of the 15 dyes chosen from the ML Regressor model predictions exhibit hydrophobicity similar to that of Cy5. Three dyes are hydrophilic (dyes 2, 4, and 13), with dye 13 being the most hydrophilic. These three dyes also exhibit relatively low $\mathsf{\mu }$ values compared to the rest of the dataset, indicating that they may not be suitable for excitonic applications that require close inter-dye separations and large transition dipole moment couplings. Conversely, dyes 7, 12, and 15 exhibit the most positive log(P_o/w) values, meaning they are estimated to be the most hydrophobic. Furthermore, dyes 7, 12, and 15 have $\mathsf{\mu }$ values within 25% of that for Cy5, making those dyes suitable for excitonic applications. Dyes 3 and 5 have log(P_o/w) values slightly larger than that of Cy5, and $\mathsf{\mu }$ values about 5 D larger than Cy5. Overall, based on our criteria—a large $\mathsf{\epsilon }$ (indicating large $\mathsf{\mu }$) and a large positive log(P_o/w)—dyes 3 and 5 are the most promising candidates in the dataset for excitonic applications.

3.2. Molecular Dynamics Simulations of Dye Aggregate–DNA Duplex Interactions

To study the effects of DNA on the dye orientations, 1 $\mathsf{\mu }\mathrm{s}$ MD simulations were performed with two dyes covalently bound to the backbone of DNA duplexes via dual phosphoramidite linkers. In our study, we started with commercially available Cy5 and Cy5.5 as reference dyes to guide other dye candidate selection. Our research group experimentally demonstrated that Cy5 can exhibit aggregation, strong absorption, and excitonic coupling when attached to DNA duplexes and Holliday Junctions [27,28,29]. Due to its similar structure to Cy5, Cy5.5 should exhibit similar properties [38]. The extra aryl groups on Cy5.5 extend the conjugation and add to the size of the molecule, which could affect dye-packing and make ${\mathsf{\mu }}_{\mathrm{Cy}5.5}$ slightly larger than ${\mathsf{\mu }}_{\mathrm{Cy}5}$, as shown in Table 1. Furthermore, as shown in Figure 6, log(P_o/w) for Cy5.5 is 40% larger than that for Cy5, indicating that Cy5.5 might pack closer than Cy5 when aggregated. Thus, we chose a Cy5 dimer and a Cy5.5 dimer for MD simulations as a benchmark for comparison with future simulations of other selected dyes.

The simulations were performed in water at a 1 atm pressure and 300 K. Dye orientations were quantified using the orientation factor, $\mathsf{\kappa }$, calculated using Equation (3). Dimer exciton exchange energies, $\left|\mathrm{J}\right|$, were quantified using Equation (4) with inputs from TD-DFT for the transition dipole moments, where ${\mathsf{\mu }}_{\mathrm{Cy}5}=15.35 \mathrm{D}$ and ${\mathsf{\mu }}_{\mathrm{Cy}5.5}=15.57 \mathrm{D}$, as shown in Table 1. The vectors corresponding to $\mathsf{\mu }$ were found to primarily reside on the long axis of the dyes, and thus, the values of ${\mathrm{r}}_{\mathrm{i}}$ and ${\mathrm{s}}_{\mathrm{i}}$ in Equation (4) were chosen as the centers of the terminal aryl groups of the dyes.

The MD results for the Cy5 dimer attached to a DNA duplex are presented in Figure 7 and Figure 8. Figure 7a shows a heatmap plot of $\left|\mathsf{\kappa }\right|$ versus $\mathrm{R}$, and Figure 7b shows a heatmap plot of $\left|\mathrm{J}\right|$ versus $\mathrm{R}$ for the 900 ns of data collection. Based on the $\left|\mathsf{\kappa }\right|$ values shown in Figure 7a, there are two distinct dimer orientations, labeled “O1” and “O2”. The approximate $\left|\mathrm{J}\right|$ regions corresponding to O1 and O2 are also shown in Figure 7b. As shown in Figure 7c, O1 corresponds to where the dyes are located outside of the duplex (i.e., non-intercalated). This orientation has a $\left|\mathsf{\kappa }\right|$ value ranging approximately from 0–1 (oblique and H-like aggregate) and an $\mathrm{R}$ of approximately 2.5–3.0 nm, which results in a relatively low $\left|\mathrm{J}\right|$ of less than 20 meV. Examining Figure 8, a shift in the orientation of the dyes occurs after around 150 ns of simulation, where the dyes re-orient and intercalate into the base-stack region of the DNA. This change in orientation results in a mostly head-to-tail (J-like) configuration with some obliqueness, corresponding to a $\left|\mathsf{\kappa }\right|$ of about 1.25–1.5 and an $\mathrm{R}$ of about 0.9–1.5 nm, as indicated by O2 and represented in Figure 7d. This more closely spaced orientation results in a larger $\left|\mathrm{J}\right|$ of roughly 40–80 meV. Since |J| was relatively stable after 200 ns, all dimer orientations beyond that time were averaged, yielding $\left|\mathsf{\kappa }\right|=1.35 \pm 0.22$ and $\mathrm{R} =1.26 \pm 0.21$ nm. Thus, averaging $\left|\mathrm{J}\right|$ over this period of time results in $\left|\mathrm{J}\right|=58.28\pm 12.74$ meV. The post-200 ns average values of $\mathrm{R}$ and $\left|\mathrm{J}\right|$ agree well with the experimentally derived values of $\mathrm{R} =1.32$ nm and $\left|\mathrm{J}\right|=48$ meV, as obtained by Huff et al. [29]. Furthermore, they determined the dimer to have a similar orientation, with a red-shift in the main absorption peak observed for the dimer relative to the monomer, indicating a mostly J-like orientation [29]. The $~10$ meV larger $\left|\mathrm{J}\right|$ value obtained from MD might be caused by a slight overestimation of ${\mathsf{\mu }}_{\mathrm{Cy}5}$ using TD-DFT, which is shown to be $~2$ D larger than the ${\mathsf{\mu }}_{\mathrm{Cy}5}$ obtained from experimental measurements [29,46].

Compared to the Cy5 dimer, the Cy5.5 dimer exhibits a similar trajectory. The Cy5.5 dyes were initialized, outside of the DNA duplex (i.e., non-intercalated). After about 400 ns, the dyes intercalated, reducing $\mathrm{R}$ and increasing $\left|\mathsf{\kappa }\right|$ (more J-like) and $\left|\mathrm{J}\right|$. The two regions corresponding to the non-intercalated and intercalated states of the Cy5.5 dimer are labeled as “O1” and “O2” in Figure 9a,b, which are represented in the snapshots in Figure 9c,d. The O1 region of the Cy5.5 dimer has an $\mathrm{R}$ that ranges from about 2.2–3.2 nm. However, $\left|\mathsf{\kappa }\right|$ for the Cy5.5 dimer has a range of roughly 0–1.5, which is larger than that for the Cy5 dimer. However, $\left|\mathrm{J}\right|$ for the O1 region of the Cy5.5 dimer is between 0–15 meV, comparable to that of the Cy5 dimer. The ranges of $\mathrm{R}$ and $\left|\mathsf{\kappa }\right|$ of the O2 region of Cy5.5 are roughly 0.8–1.7 nm and 1.2–1.75, respectively, similar to that of the Cy5 dimer. This range of orientations results in a $\left|\mathrm{J}\right|$ range of about 30–60 meV, slightly smaller than that of the Cy5 dimer.

Averaging the dimer orientations past 400 ns (after which $\left|\mathsf{\kappa }\right|$ and $\mathrm{R}$ are relatively stable, as shown in Figure 10) results in $\left|\mathsf{\kappa }\right|=1.33\pm 0.15$ and $\mathrm{R} =1.17\pm 0.17$ nm. Similarly, averaging $\left|\mathrm{J}\right|$ over this time period results in $\left|\mathrm{J}\right|=44.98 \pm 11.6$ meV. Even though the two dimers exhibit similar orientations, the smaller $\left|\mathrm{J}\right|$ for the Cy5.5 dimer might be caused by a small structure difference between Cy5.5 and Cy5. Cy5.5 is slightly larger and longer than Cy5 (by about 0.1 nm) due to the extra aryl groups at the two ends of the dye. Despite the smaller average $\mathrm{R}$ and similar orientations, the larger dye length of Cy5.5 compared to that of Cy5 results in a smaller pre-factor term ${\mathrm{J}}_0$, leading to the smaller $\left|\mathrm{J}\right|$. In a similar study, it was found that Cy5.5 homodimers attached to transverse strands on a DNA Holliday junction exhibited closer dye distances and larger $\left|\mathrm{J}\right|$ values than Cy5 homodimers [38]. A potential reason for this difference might be that the DNA Holliday junctions are more flexible than the duplexes, which allow dye orientations to promote larger J values.

4. Discussion

Based on the results shown in Figure 5, the developed Regressor is useful for the prediction of experimental $\mathsf{\epsilon }$ trends. This was shown by predicting the $\mathsf{\epsilon }$ of a dataset of over 18,000 molecules, from which 15 were identified using a combination of Regressor predictions and TD-DFT calculations. In general, according to desired criteria the developed ML model can quickly and accurately screen a large dataset of dyes for optimum $\mathsf{\epsilon }$ for further study, which could save computation and experimental time.

There are improvements that can be made to the developed ML models. Primarily, more data would improve the chemical space on which the ML models would be trained, such as more chemical structures, as well as the inclusion of the solvents used. This is highlighted in Figure 1, which shows a disproportion in large versus small $\mathsf{\epsilon }$ values in the training set. Furthermore, the present model could benefit from optimization and feature engineering, which could further reduce the computation time and help to improve the model by identifying more features that may better describe dye properties. Other types of ML models such as Support Vector Machine (SVM) algorithms could also be explored, which could potentially lead to better predictions. Finally, an improvement could be implemented by including human knowledge in the workflow so as to identify desirable structural features; this could be useful in interpreting data, improving efficiency, and enhancing model performance.

Examining Figure 5 and Table 1, the trends of the TD-DFT-calculated $\mathsf{\mu }$ and the available experimentally determined $\mathsf{\epsilon }$ mostly agree when considering the upper range of $\mathsf{\epsilon }$ values. The deviations in $\mathsf{\epsilon }$ compared to $\mathsf{\mu }$ could be due to the solvent chosen for TD-DFT calculations (in our case, water) or the exchange–correlation functional used for the calculations. For example, our prior studies showed that the CAM-B3LYP functional produced slightly different $\mathsf{\mu }$ values for similar dyes, however, the overall trend remained the same [46]. In the case of the Cy5 derivatives, the $\mathsf{\epsilon }$ values might be affected by the specific DNA strand to which the dyes are attached, leading to the range in $\mathsf{\epsilon }$ shown in Figure 5 for those dyes. This could also lead to some degree of disagreement with the TD-DFT calculations of $\mathsf{\mu }$ (only considering free dyes in water) and the ML-predicted $\mathsf{\epsilon }$. Furthermore, in general, there is no strong correlation between the ML-predicted $\mathsf{\epsilon }$ and TD-DFT $\mathsf{\mu }$ for the 15 dyes shown in Table 1. This highlights the need for TD-DFT calculations as a validation for the initial dye screening using the ML Regressor, as the ML model is not 100% accurate. Studying the values of $\mathsf{\mu }$ in Table 1, the main dye feature leading to a large $\mathsf{\mu }$ is long conjugated chains, as is the case for dyes 3 and 5. The same trend is observed for Cy5 and Cy7, which have two and four more carbons in the conjugated chain compared to Cy3, respectively. Despite having similar structures, dye 10 has a ~1 D larger $\mathsf{\mu }$ compared to dye 8, which might be caused by the Cl atom bonded to the polymethine chain in dye 8.

Examining Figure 6, the Cy5 derivatives modified with hydrophobic groups exhibit an enhanced hydrophobicity (as exemplified by more positive $\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)$ values), and the trend of $\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)$ follows the same trend obtained by Meares et al. who used the Percepta Platform to predict $\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)$ based on chemical structure [87]. Since these dyes are derived from the Cy5 dye, which has been shown to form aggregates when attached to DNA scaffolds [27,28,29], they are promising candidates for further studies with an aim of enhancing $\left|\mathrm{J}\right|$. Examining the $\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)$ for the 15 ML-selected dyes, most exhibit values similar to Cy5 and its derivatives, with the notable exceptions of dyes 2, 4, 7, 12, 13, and 15. Dyes 2, 4, and 13 have a negative $\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)$ and are, therefore, expected to be hydrophilic. Dye 2 is structurally similar to dyes 1, 7, and 12, which have positive $\log \left({\mathrm{P}}_{\mathrm{o}/\mathrm{w}}\right)$ values. The hydrophilicity of dye 2 might be attributed to the Se atoms in the ring, which dyes 1, 7, and 12 do not contain. Dyes 4 and 13 contain S groups, which are known to increase hydrophilicity [101]. Dyes 7, 12, and 15 stand out as the most hydrophobic of the dyes tested and all belong to the same family of porphyrins, which are known to be hydrophobic [102,103]. Based on their relatively high $\mathsf{\mu }$ and enhanced hydrophobicity, dyes 7, 12, and 15 (porphyrin-based dyes) are promising candidates for dye–DNA applications. Furthermore, it has been shown that porphyrin dyes are suitable for bonding to DNA and are able to intercalate between bases, which may be beneficial for dye aggregation in DNA [102].

Our MD simulations show that the Cy5 dimer and Cy5.5 dimer orient similarly when attached to the same positions on a DNA duplex. Comparing the simulation results for the Cy5 dimer to Huff et al., our predicted $\left|\mathrm{J}\right|$ and $\mathrm{R}$ are within 20% and 5%, respectively [29]. It should be noted that the results shown by Huff et al. were obtained using a program developed by our group based on the Kuhn–Renger–May (KRM) theory, which only considers absorbance and circular dichroism spectra to obtain dye orientations [29]. However, the KRM-based program does not provide information on the details of the DNA–dye orientations or the impact of the DNA on the dye orientation as our MD simulations do. Based on the MD results, one possible mechanism is dye intercalation into the DNA duplex from outside of the DNA backbone to form a J-dimer, which occurs for both Cy5 and Cy5.5 simulations. Our results highlight the effectiveness of MD simulations for predicting dye orientations and the excitonic properties of dyes bound to DNA duplexes. For future studies, we will continue applying MD to additional ML-selected dye aggregates in DNA scaffolds.

5. Conclusions

ML models were developed through training and validation on a set of 8802 molecules to predict $\mathsf{\epsilon }$ based on dye structures. A Random Forest Classifier was developed with an accuracy of 97%, and based on the performance of the Classifier, a Random Forest Regressor was constructed. Comparing ML-predicted $\mathsf{\epsilon }$ from the Regressor and experimental $\mathsf{\epsilon }$, the Regressor was found to have a maximum ${\mathrm{R}}^2$ of 0.95. Using the Regressor, the $\mathsf{\epsilon }$ values of molecules in a dataset of around 18,000 were predicted, and the top 100 dyes were used in TD-DFT calculations to calculate $\mathsf{\mu }$. Overall, 15 dyes were identified to have a relatively large $\mathsf{\mu }$ comparable to Cy5. MD simulations were conducted on reference dyes (Cy5 and Cy5.5) to determine the dye dimer orientations and transition dipole moment couplings, $\left|\mathrm{J}\right|$. For Cy5, the MD simulations were able to predict dye orientations and $\left|\mathrm{J}\right|$ within 20% of the experiment. The Cy5.5 dimer yielded similar results to the Cy5 dimer. The successful use of the combined ML and DFT/TD-DFT screening to identify dyes with a large $\mathsf{\epsilon }$ and $\mathsf{\mu }$ highlights the effectiveness of our workflow to screen numerous dyes for desired properties. The agreement of our MD simulations with the experiment show that we can accurately detect dye dimer properties when attached to DNA scaffolds, which is crucial to guide dye design and synthesis for excitonic applications.