Impact of Doxorubicin on Self-Organization of Congo Red: Quantum Chemical Calculations and Molecular Dynamics Simulations

Quantum-chemical calculations and molecular dynamics simulation were applied to a model self-organization process of Congo red (CR) molecules in aqueous solution and the impact of doxorubicin (DOX) molecules on such a process. It was demonstrated that both pure CR/CR and mixed CR/DOX dimers were stable. Van der Waals interactions between aromatic units were responsible for a stacked dimer formation. An important source of stabilization in the CR/CR dimer was the polarization energy. In the CR/DOX mixed dimer long range, electrostatic interactions were the main driving force leading to complexation. An implicit solvent model showed that the formation of the CR/CR dimer was favored over the CR/DOX one. Molecular dynamics simulations demonstrated rapid complexation. In the pure CR system, short sequences of ribbon-like structures were formed. Such structures might be glued by hydrogen bonds to form bigger complexes. It was shown that the aromatic part of the DOX molecule enters CR ribbons with the sugar part covering the CR ribbons. These findings demonstrated that CR may find applications as a carrier in delivering DOX molecules; however, further more extensive investigations are required.


■ INTRODUCTION
Development of supramolecular chemistry opens up new horizons for many fields of science. The chemistry of such molecular systems is ruled by weak intermolecular interactions, e.g., van der Waals forces, hydrogen bonds, and π−π and electrostatic interactions. The presence of noncovalent interactions between organic molecules allows the formation of large self-assembled supramolecular structures. 1−3 This kind of multimolecular structure is characterized by high elasticity and a large surface area. These properties enabled some supramolecular associates to form complexes with proteins, which in general are compounds that are unable to interact at sites other than the active site. 4 Congo red (CR) is an example of a self-assembling supramolecular compound creating ribbon-like structures related to the association of flat, polyaromatic parts of these molecules. It belongs to the azo dye group, and it is used especially as an indicator and biological stain. CR has been used for over a century to demonstrate the presence of amyloidal deposits in tissue, e.g., in brain tissue, or in the heart 5 or kidneys 6 in the diagnosis of Alzheimer's disease.
Amyloid fibers bind to Congo red, giving the effect of applegreen birefringence under a polarized light microscope, which remains the gold standard in recognizing amyloid deposits. 7,8 Congo red is also often used as a model dye in evaluating the effectiveness of various materials that catalyze the degradation of organic pollutants in water. 9,10 The structure of Congo red is symmetrical. Its central part consists of a biphenyl ring linked through azo bonds with naphthalene rings substituted with sulfone and amino groups. Individual molecules interact with each other to form larger supramolecular, ribbon-like structures. The polyaromatic structure of Congo red enables intercalation of various compounds with planar groups, e.g., the doxorubicin (DOX) drug.
Doxorubicin is an anthracycline antibiotic commonly used as an anticancer drug. The mechanism of doxorubicin action is associated with the formation of a stable complex with DNA, thus inhibiting further cell division ultimately causing its death. Unfortunately, chemotherapy is not an ideal treatment method. Cytostatics used in this treatment destroy all rapidly dividing cells, including normal cells, which causes numerous side effects for the human body. Currently, new solutions are sought in order to reduce the therapeutic dose of the drug while maintaining its effectiveness and reducing its toxicity to normal cells. One of these possible solutions is the use of drug carriers.
The self-assembled ribbon-like structures that form elongated Congo red systems owing to their unique properties are used in targeted drug delivery to cancer tissue. 10 Such systems are an example of a new type of ligand for proteins to which they bind using a nonclassical mode of interaction. 4 An example is binding to antigen−antibody complexes while not binding to free antibodies. 11 This interaction is the basis for the use of such systems in immunotargeting. 12 At the same time, CR systems effectively bind to various molecules, including drugs, to form co-micelles with them. 13 Research to date has demonstrated the possibility of using CR systems in vivo as potential drug carriers. These systems easily bound to the immune complexes formed in the body and then were gradually removed. This gives the potential for use in immunotargeting, and an additional advantage is the easy removal of excess unused drug from the body. 11 In this work, the formation of Congo red supramolecular assemblies is investigated and its possible use as a carrier system for doxorubicin, a model anticancer drug, is considered. We focus on the microscopic description of CR and CR/DOX clusters, including interactions responsible for dimer stabilization, the influence of solvent (water), and the mechanism of the clustering process in aqueous solution.

■ COMPUTATIONAL DETAILS
All quantum-chemical calculations were carried out using the Gaussian 09 suite of programs. 14 The structures of DOX and CR (see Figure 1) were located at the density functional level of theory, i.e., the hybrid B3LYP functional with a 6-311G(d,p) basis set. The Grimme dispersion correction was taken into account. 15 The located optimum structures of the monomers were used to prepare the initial structures of CR/CR and CR/ DOX dimers. These structures were further optimized at the B3LYP/6-311G(d,p) level of theory. All located minima on the energy surface were verified by vibration frequency analysis. The nature of interactions was characterized by a selfconsistent charge and configuration method for subsystems (SCCCMS). 16 The SCCCMS scheme has been described elsewhere. 16  η BB − 2η AB > 0. The reactant chemical potentials (μ A and μ B ) and hardness values (η AA , η BB , and η AB ) are given by the following partial derivatives: , is the external potential due to nuclei. Here, we do not consider derivatives with respect to v(r⃗ ). In other words, the positions of nuclei were fixed. Notice that the first (electron density) and second-order energy derivatives (linear response function) with respect to the external potential as well as the mixed second-order derivatives (Fukui function) are of special importance in the theory of chemical reactivity. For details, see the excellent works in the field of conceptual density functional theory. 18 The information from the energy surface can be utilized to compute for the charge reorganization accompanying P and CT processes as well as the CT Fukui function and its components. 21 Additional information regarding affinity of CR and DOX toward self-organization was derived from molecular dynamics (MD) simulations. The CR/DOX systems were previously discussed by Jagusiak and Panczyk; 22 however, the AMBER force-field parameterization of CR did not support selforganization of these molecules, which was observed experimentally. 23 This divergence can be attributed to the charge distribution that in the AMBER force field is derived by fitting to the molecular electrostatic potential of the isolated molecule. For this reason, CR and DOX molecules were parameterized in the CHARMM force field using a force field toolkit of VMD graphical software. 24 In contrast to the AMBER charges, the CHARMM charges are derived by reproducing the interaction energies of hydrogen donor and hydrogen acceptor atoms of the CR or DOX molecule with water molecules. It can be expected that this will improve the description of CR and CR/DOX aggregates.
The MD calculations were carried out using the NAMD package 25 with periodic boundary conditions imposed on the system. An all-atom CHARMM force field 26,27 was used. The TIP3P model 28 was adopted for water. The VMD software was applied to visualize the trajectories and compute selected quantities. Density-based cluster analysis of the trajectories was performed with the dbcscan function of the R package. 29 In all calculations, the number of atoms (N) in the simulation box was greater than 200,000. The initial configuration was obtained by placing the target molecules on a cube grid. A random orientation of each molecule was used. The CR anions were neutralized by sodium cations. Simulations at constant temperature (T = 293 K) and pressure (p n = 1 bar) were employed. These simulations were preceded by energy minimization and 10 ns canonical (NVT, number of particles, volume, and temperature) ensemble simulations. The (Np n T) ensemble was applied for 120 ns. The ratio of the x and y box sizes was constant. Five independent trajectories were computed for each molecular system, and their average is reported as the final result.

■ RESULTS
In Table 1, we have listed the interaction (INT) energies and their components for dimer systems. To the best of our knowledge, CR/CR and CR/DOX dimers were not described previously in the literature. The interaction energy was computed using a supermolecular approach: where the energies E A and E B are calculated for monomers having the same geometry as in the dimer. The basis set superposition error was eliminated using the counterpoise method of Boys and Bernardi. 30 The total binding (b) energy should be supplemented by the deformation energy (shown in the last row of Table 1) This term is positive, and the energies of the reactants in the minimum structures E A 0 and E B 0 are more negative than in dimer geometries. As one would expect, P and CT energies are stabilizing (negative) while DEF is destabilizing (positive). The highest destabilization contribution to the INT energy in the CR/CR system is the ES term. This behavior was not surprising because the CR is negatively charged. This energy term in a mixed CR/DOX system is negative; therefore, the ES term should enhance the self-organization of molecules or at least promote the incorporation of DOX molecules into ribbon-like supramolecular structures of CR. The polarization energy is the most stabilizing effect in the CR/CR dimer. The polarization contribution in CR/CR is almost three times higher than in the CR/DOX dimer. The second highest is VdW energy. Except DEF energy, the remaining contribution to CR/DOX stabilizes the dimer formation. Dimerization enthalpies and Gibbs free energies are summarized in Table 2. Data corresponding to vacuum and the implicit solvent model, namely, the polarizable continuum model (PCM), 31 are reported. The temperature and entropy effects did not change the trends observed in INT (b) energies. Binding enthalpy and Gibbs free energies collected in Table 2 favor the formation of a CR/DOX dimer. This picture changes for the PCM water model. The formation of a CR/CR dimer is  Table S1. The same conclusions follow from reported data, namely, the reverse of stability of CR/DOX and CR/CR dimers when going from vacuum to the implicit water model. Additional support for the observed reverse in stability of CR/ DOX and CR/CR dimers when going from vacuum to the implicit water model results from polarization ( Table 1). The stabilization due to polarization energy in CR/CR dimers is three times greater than in the CR/DOX dimer. The selfconsistency in PCM calculations, namely, the polarization of solute and backpolarization of dielectric medium, will depend on this particular feature of charge distribution.
To understand the nature of stabilization in dimers, quantum theory of atoms in molecules (QTAIM) was applied. Topological properties of electron density have been shown to correlate with the strength of intermolecular interactions. 32,33 Here, bond critical points (BCPs) were used to detect the bonding pattern. The existence of bond paths between the atoms of different monomers is a clear indicator of bonding contacts. Further classification can be done by the value of the electron density at the BCP, ρ BCP . In the CR/CR system, there are 28 BCPs between CR monomers. The value of the electron density at BCPs ranges from 0.004 to 0.031 a.u. In this variety of BCPs, there are four hydrogen bonds of the N−HO type (ρ BCP ≈ 0.03 a.u.) and four hydrogen-like bonds of the C− HO type involving aromatic carbon atoms (ρ BCP ≈ 0.01 a.u.). The Laplacian of the electron density at the hydrogen bonds is within the 0.038−0.105 a.u. range. Both sets of data at bond critical points are in agreement with topological criteria of the existence of hydrogen bonding as proposed by Koch and Popelier. 33 In accordance with these rules, there is a loss of electron density in the area of hydrogen atoms. It is reflected by negative isosurfaces Δρ P (r⃗ ) = ρ P (r⃗ ) − ρ 0 (r⃗ ) shown in Figure S1 of the Supporting Information. There is no charge flow due to charge transfer, so the change in density is dominated by polarization process: Δρ M (r⃗ ) = ρ M (r⃗ ) − ρ 0 (r⃗ ) ≈ Δρ P (r⃗ ). The charge distribution in the CR/CR dimer is softer than in the CR/DOX dimer. Namely, the changes in Δρ M (r⃗ ) (CR/CR dimer) are spread over the whole system, while in the mixed dimer, they have a more local character. This observation also supports the reverse of stability of CR/DOX and CR/CR dimers in the implicit water model compared to vacuum. The remaining contacts are of van der Waals nature: 4 NC (ρ BCP ≈ 0.007 a.u.), 13 CC, 2 HC, and 1 NN contacts with the Laplacian of the electron density positive and lower than 0.02 a.u. Majority of detected contacts are responsible for stacking interactions. The number of bonding paths in the CR/DOX dimer is equal to 22: one relatively strong O−HO hydrogen bond (ρ BCP ≈ 0.04 a.u., Δρ BCP ≈ 0.13 a.u.), six HO, five HN, three CC, three CH, one NO, one NC, one CO, and one HH contacts. The value of the electron density at the remaining contacts ranges from 0.003 to 0.014 a.u. Molecular graphs obtained from QTAIM analysis for both dimers are depicted in Figure 3.
Molecular Dynamics Simulations. Vacuum and PCM calculations performed for both dimers indicate the increasing ability of CR to form organized structures in water. Such findings are qualitative since no explicit water molecules were taken into account. This was the main reason for using molecular dynamics simulations. For both systems (CR/CR and CR/DOX), the initial configuration was random without any self-organizing pattern. Practically from the beginning of simulations, the formation of dimers could be observed, which subsequently aggregated into larger clusters. Clustering patterns for the last frames of selected CR/CR and CR/ DOX trajectories are shown in Figure 4. Panel a corresponds to the pure CR system, while panel b corresponds to the mixed equimolar CR/DOX molecular system. Clusters of different sizes can be observed. Sample clusters, corresponding to the sizes of 8 and 13 molecules, are highlighted with red surfaces and drawn in detail in panels c and d. In the case of the CR system, the molecules form two ribbon-like substructures that are glued by hydrogen bonds. A qualitatively similar picture is observed in a mixed CR/DOX system. Ribbon-like structures are also present with the sugar moiety of DOX molecules protruding to the surface and the aromatic moiety built into CR ribbons. The clustering patterns for the remaining CR and CR/DOX systems are shown in Figures S2 and S3, respectively.
Changes in the average number of monomers forming a cluster during the courses of both types of simulations are presented in Figure 5. The resultant curves were averaged over five trajectories. Deviations from the average values are shown in the Supporting Information (Figures S4 and S5). Two approaches were employed to define the center of mass (CM) of DOX molecules in the CR/DOX system. In the first one (red line in Figure 5), the CM was computed as the center of mass of the whole DOX molecule. In the second case (green line), only the aromatic part of a DOX molecule was taken into account. The differences are associated with the mechanism of the clustering process; the aromatic part of the DOX molecule is incorporated into the ribbon while the sugar part protrudes

ACS Omega
http://pubs.acs.org/journal/acsodf Article from the ribbon surface. For CR molecules, both approaches give the same result. In both systems, a saturation effect can be observed. After saturation, in the CR system, there are on average five CR molecules in the ribbon. For a mixed CR/ DOX system, the average number of molecules is higher and amounts to 7 or 8 molecules depending on the approach employed to calculate the center of mass. For this system, larger fluctuations in the average cluster size can be observed.
They are connected with greater mobility of DOX molecules associated with the surface of larger clusters whose aromatic parts are not incorporated into the ribbons. Fluctuations observed between the successive points indicate that the monomer can detach from the aggregate. However, eye examination of the trajectories reveals that this detachment is only temporary and the displaced molecules quickly return to the clusters. Such temporary dissociation is more frequent for the mixed CR/DOX system than for the CR one due to the presence of a hydrophilic sugar group in DOX and a more disordered structure of CR/DOX clusters. The clustering process is summarized in Figure 6. The histogram illustrates the formation of clusters of various sizes.
Data is acquired from the terminal parts (after saturation, Figure 5 and Figures S4 and S5) of all five trajectories obtained for each system. In summary, ca. 1500 frames were gathered for both CR and CR/DOX systems. Histograms for individual runs are shown in Figures S6 and S7 of the Supporting Information. Again, two approaches were used to calculate CM of DOX molecules. In the case of the CR system, clusters composed of 2 to 17 molecules can be observed. The distribution of cluster sizes is relatively narrow and exhibits a clear maximum at the value of four molecules in a cluster. In this system, free monomers are practically not observed. In contrast, the cluster size distribution in the CR/DOX system is broader and extends toward larger clusters. The largest cluster observed in the simulations consisted of as many as 31 molecules. Surprisingly, free monomers also occur in this system. In the calculations employing CMs of whole DOX molecules, they are the most abundant entity in the system. However, eye examination of the trajectories reveals that these are not free monomers but molecules (mostly DOX) loosely associated with the surface of larger clusters. Among aggregates, the most populated number of molecules in a cluster is 4. However, the maximum is less pronounced than in the CR system and larger clusters are also present. Two approaches used in cluster analysis give qualitatively similar predictions. A higher green bar than a blue one indicates that whole DOX molecules cover the ribbon and the aromatic part . Visualization of the last MD simulation frames for (a) pure CR and (b) mixed CR/DOX molecular systems. Example clusters, marked in red color, are depicted in panels (c) and (d) together with all water molecules (cyan color) and Na + cations (green balls) within a 10 Å distance from the cluster. In panels (a) and (b), CR molecules are drawn in gray and DOX molecules in black. In panels (c) and (d), CR molecules are colored according to atom names (C -black, Ored, N -blue, S -yellow, and H -gray). In DOX molecules, the aromatic part is drawn in orange and the sugar part in magenta.

ACS Omega
http://pubs.acs.org/journal/acsodf Article is completely immersed in water. This can be considered as the first step in incorporating DOX molecules into CR ribbons. In cluster analysis, information about the source of stabilization is lost. To get some information regarding the nature of stabilization (bonding contacts) in clusters, the radial distribution function (RDF) was plotted. The RDFs between the center of mass in pure and mixed systems are shown in Figure 7a. In this graph, the CM of the aromatic part of DOX molecules was taken into account as the plots based on the CM of the whole molecule do not provide an accurate picture of the interactions between the aromatic parts of the molecules. Panel b in Figure 7 presents the RDF describing the ability of the formation of hydrogen bonds between the oxygen atoms from the SO 3 group and DOX oxygen and nitrogen atoms. In Figure 7a, an oscillatory pattern is observed. This indicates the formation of ribbon-like structures. Such stacked structures are stabilized by van der Waals interactions as demonstrated by the bond critical paths shown in Figure 3 for CR/CR and CR/DOX dimers. This type of interaction occurs in pure CR (dashed red line) as well as in a mixed CR/ DOX system. In a mixed system, the RDF for CMs of CR/CR, CR/DOX, and DOX/DOX molecules is plotted. Hydrogen bonds are another important source of stability between CR and DOX molecules. Panel b in Figure 7 indicates the formation of hydrogen bonds between O(SO 3 )H−O and O(SO 3 )H−N. The former hydrogen bonds are particularly strong due to the flexibility of the DOX sugar unit. The BCP in the CR/DOX dimer has the highest value of electron density (0.04). The later hydrogen bonds are weaker with a ρ BCP close to 0.01.

■ CONCLUSIONS
It was demonstrated that Congo red molecules self-assemble in water in accordance with experimental data. This indicates that proper force-field parameters for CR were obtained. Ribbonlike structures of varying lengths are formed. These flat stacked structures are stabilized by van der Waals contacts. Hydrogen bonds are responsible for gluing short ribbons of CR. In the mixed CR/DOX system, doxorubicin is easily incorporated into the CR ribbons. The aromatic parts of DOX penetrate the CR ribbon, while the sugar part covers the surface of the ribbon. Hydrogen bonds are responsible for covering the CR ribbon. In the mixed CR/DOX system, larger clusters are formed than in a pure CR system. Weak intermolecular interactions are responsible for complexation. Self-assembled CR structures can be used as carriers for DOX molecules.

* sı Supporting Information
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acsomega.0c01095. Table S1: binding enthalpies and Gibbs free energies computed at the M062X/6-311G(d,p) level of theory; Figure S1: isosurfaces illustrating the polarization process in CR/CR and CR/DOX dimers; Figure S2: clustering patterns for the remaining four CR trajectories; Figure S3: clustering patterns for the remaining four CR/DOX trajectories; Figure S4: average number of monomers in CR clusters; Figure S5: the average number of monomers in CR/DOX clusters; Figure S6: cluster size distribution in CR systems; and Figure S7: cluster size distribution in CR/DOX systems (PDF)