The self-association equilibria of doxorubicin at high concentration and ionic strength characterized by fluorescence spectroscopy and molecular dynamics simulations

Abstract The self-association equilibria of doxorubicin hydrochloride (DX), at high drug and NaCl concentrations, are studied by temperature scan fluorescence spectroscopy, with the support of molecular dynamics (MD) calculations. Even though all anthracyclines show dimerization equilibria, DX only can further associate into long polymeric chains according to DXmon ⇄ DXdim ⇄ DXpol. This is reflected not only in the mechanical properties of DXpol solutions (behaving as thixotropic gels) but also in their spectroscopic behaviour. Fluorescence, in particular, is the technique of election to study this complex set of equilibria. Upon increasing the temperature, DXpol melts into DXdim, which in turn is in equilibrium with DXmon. Since DXdim is non fluorescent, with a fluorescence temperature scan experiment the DXpol⇄ DXmon equilibrium is probed. However, also information on the DX dimerization equilibrium can be derived together with the relevant thermodynamic parameters ruling the dimerization process ( Δ H d i m ° = - 56 k J m o l - 1 ; Δ S d i m ° = - 97 J m o l - 1 K - 1 ). The residence time of DX molecules in the dimer (74.7 μs), as well as the monomers mutual orientation in the dimer, are characterized by means of theoretical and computational modelling.


Introduction
The anticancer antibiotic doxorubicin hydrochloride (DX -Chart 1) is a powerful drug used for the treatment of several malignancies. The first trials with doxorubicin (at that time named adriamycin) started back in 1969; thereafter the drug was approved for clinical use with exceptional results in spite of the many toxic side-effects [1] and the development of multidrug resistance. [2] Very active research started very soon to try to overcome these drawbacks, and it is still alive, trying to find new formulations or new administration strategies. [1][3] [4] A C C E P T E D M A N U S C R I P T From a physico-chemical point of view, DX is an extremely intriguing molecule. Despite its relatively low molar mass (580 g/mol), it exhibits features typical of high molecular weight species, such as the capability to act as a water gelator. This feature has been described in several papers back in the nineties [5] [6][7] [8] and has recently been the subject of a deep investigation that revealed the peculiar association of DX molecules into supramolecular aggregates made of hundreds of units (distorted dimers). [9] Interestingly, the closely related molecule epirubicin hydrochloride, that differs from DX only for the stereochemistry at the C4' carbon atom in the aminosugar moiety of the drug, does not form gels in the same conditions that induce DX gelling. Molecular dynamics simulations revealed that this was due to the critical network of hydrogen bonds holding together the supramolecular assembly of DX in the gels, possible only for the specific stereochemistry of donors and acceptors fulfilled in the DX molecule. [9] The stereochemically driven supramolecular polymerization of DX is a consequence of the capacity of DX to form dimers in solution, even though the gelification is not just the extension to many molecules of the dimerization process. In fact, dimerization is a characteristic shared by all anthracyclines, [10] while gelification, as already said, is a DX peculiarity. Self-association of anthracyclines has been the subject of several investigations, [11][12] [13][14] [15] however an analysis of the all set of DX association equilibria by means of fluorescence spectroscopy, is not present, to our knowledge, in the literature. Among anthracyclines, DX is the only molecule, so far, to show a complex pattern of selfassociation equilibria that goes beyond the classical monomer ⇄ dimer equilibrium. When the analytical concentration exceeds the Critical Gel Formation Concentration (CGFC), [8] even in the absence of salts, as we recently demonstrated, [9] the complete set of equilibria DXmon ⇄ DXdim ⇄ DXpol has to be taken into account to properly describe the DX solutions behaviour. In the present study an extension of the analysis of the association equilibria of DX leading to its polymerization into gels will be presented exploiting the joint use of fluorescence spectroscopy, theoretical modelling and Molecular Dynamics (MD) simulations. In addition, the thermodynamic properties of these equilibria have been characterized. It has to be pointed out, finally, that only at high DX concentration and in the presence of salt the supramolecular DX polymerization leading to gels occurs, making thus possible to probe the whole set of association equilibria above introduced.

Materials and Methods
Doxorubicin hydrochloride was a kind gift of Farmitalia Carlo-Erba Chem. Co. All the other chemicals used, of the highest available purity, were from Sigma-Aldrich. Ultrapure water (18.2 Chart 1 Doxorubicin hydrochloride (DX) chemical structure. The red arrow indicates the long axis of the molecule.
MΩ/cm) from Arium pro UV system (Sartorius Stedim Biotech) was used in the preparation of the solutions. Samples were prepared by diluting a 0.02M solution of DX with suitable amounts of 5M NaCl stock solution and water directly into a 1.5 mL Eppendorf vial to obtain the final DX concentration of 0.01M and the desired final NaCl concentration. Samples were vortex mixed for several minutes and then centrifuged at 1500 r.p.m. (Hettich Mikro 20). They were left to equilibrate in the dark at room temperature for 30 minutes before transferring them into the cuvette for the spectroscopic measurements. Steady-state fluorescence measurements were performed on a Fluoromax 2 instrument (Horiba-Jobin Yvon), equipped with a Peltier thermostatted LUMA 40/JY1 sample holder (Quantum Northwest) controlled by the TC1 unit via the T-App software from the same manufacturer. Sample temperature was measured through a QNW-1 thermistor probe located in close proximity of the 0.01 cm sandwich cells (Hellma) used, directly connected to the TC-1 unit (accuracy: ± 0.01°C). The cell was oriented at 45° respect to the exciting beam (angular geometry) paying attention to avoid direct reflection reaching the emission monochromator. Temperature was varied using a ramp temperature of 0.5 °C/min and once reached the target temperature, the sample was left for further 4 min to ensure thermal equilibration. Samples were excited at 410 nm and the excitation and emission slits were both set at 2.5 nm; integration time was 0.25 s and spectra were acquired in the 500-800 nm region, with a step of 0.5 nm. The temperature range explored along fluorescence T-scan experiments broadly covered the interval 293-328 K, depending on samples. The all-atoms MD simulations were performed using the Gromacs software package version 2018.1 [16] on two DX molecules solvated in a cubic box using the SPC water model. [17] The force field for DX, successfully used in our previous work [9] and obtained from the ATB server, is derived using a multistep process in which results from quantum mechanical calculations are combined with a knowledge-based approach to ensure compatibility with a specific parameter set of well-known force fields. [18] After thermalization, an equilibration phase of 10 ns has been performed to let the systems reach their local minimum. Then, two MD simulations of two DX molecules in water (at 300K and 450K lasting 200 ns each) were performed in the canonical ensemble (NVT) with periodic boundary conditions and using the v-rescale temperature coupling. [19] The temperature of 450 K has been chosen after a sensitivity analysis obtained by short MD simulations in the interval between 300 K and 500 K. Such data indicate that the temperature of 450 K guarantees the sampling of a significant number of monomer-dimer transitions in a reasonable amount of time (vide infra). Following our previous works, [20][21] the box volume (~ 91 nm 3 ) at each temperature has been tuned to reproduce the pressure of 560 bar, the value at which the SPC model has a density corresponding to liquid water, i.e. 55.32 mol/L. [20] In this way the solute-solvent systems could be considered as inserting a solute into a solvent box at the same temperature and pressure of the reference solvent box, thus mimicking the experimental conditions of solvating a solute into liquid water. The bonds were constrained by means of the LINCS algorithm [22] and the particle mesh Ewald method [23] was used to compute long range interactions. A cut- off of 1.1 nm was used. The atomic coordinates were written every 2 ps and all the analyses have been performed on 50000 frames.

Results and Discussion
Fluorescence Spectroscopy All anthracyclines and DX in particular, are fluorescent molecules (the fluorescence quantum yield of the DX monomer is Φ DX mon = 3.9 • 10 −2 )[9] but their fluorescence dramatically drops upon dimerization (Φ DX dim ≈ 10 −5 ). [24] When above a critical concentration of both the drug and NaCl (the so called CGFC [8]), the DX solutions turn into gels and we have shown that this process is due to a supramolecular polymerization of DX in a two-state process between the gel and dimers, leading to the formation of DXpol. [9] The fluorescence spectrum of DXpol is markedly different from that of DXmon and since the DX supramolecular polymerization is a T-controlled process, upon heating a DXpol sample, due to the non-fluorescent nature of DXdim, the formation of only DXmon was probed (even though CD spectroscopy revealed that DXpol melts into DXdim solution). [ [25] In the case of DXpol samples, upon increasing the temperature, since they melt into non fluorescent dimers, a decrease in the fluorescence intensity is firstly observed followed by a shift in the maximum emission wavelength (λ max DX pol = 610 nm; λ max DX mon = 495 nm) and an increase of the fluorescence intensity. This feature can be appreciated in Figure 1 -Panel A where the spectra evolution with temperature of a DXpol sample is shown. Along a fluorescence temperature scan experiment, therefore, it is the equilibrium DXpol ⇄ DXmon that is probed. By considering the DXpol emission spectrum acquired at low temperature as due only to polymerized DX molecules, and that at the highest temperature as due only to DXmon, we can fit the intermediate emission spectra as the weighted sum of these two limiting spectra, i.e. Epol and Emon, according to E (T) = χ pol E pol + ϕ mon E mon An example of the results obtained is in Figure 1 -Panel B. While the weight factor for DXpol is the true molar fraction of polymerized DX (χpol; at low temperature and in the presence of salt all the DX molecules are in the form of polymer, as evidenced by time resolved fluorescence data [9]), that of DXmon is only the fraction (ϕmon) with which the spectral shape of DXmon contributes to the total fluorescence emission spectrum at a certain temperature (i.e. ϕmon∝ χmon). In the case of DX samples at the very same concentration but without NaCl, upon increasing the temperature, the fluorescence intensity linearly increases since, as already said, in this situation only the DXdim ⇄ DXmon equilibrium is active and it moves to the right upon heating the solution. By assuming that at the highest probed temperature (52-53 °C, depending on samples) all the drug is in monomeric form, a normalization of all the spectra for this one allowed the weight factor of the DXmon in the absence of salt to be calculated. In Figure 2 the result of these data treatment is summarized and it is evident that once DXpol turned into a DXmon+DXdim solution, the systems behaviour is salt independent (all the ϕmon data coming from fluorescence emission spectra calculated at [NaCl]≠ 0.00 M lie on the data obtained at [NaCl]= 0.00 M). From the data in the absence of salt, the proportionality constant between ϕmon and χmon can be calculated since Kdim is known from the literature. By taking one of the literature data for the DX dimerization constant ( 298 = 6.3 • 10 4 −1 ), [25] the true value of DXmon molar fraction at 298K can be calculated thus allowing to solve the relation ϕ mon = κ • χ mon . The proportionality constant κ will depend only on the different fluorescence quantum yields of the species present in the system and on instrumental parameters. As long as the DX concentration and the acquisition parameters remain the same, it is possible to compare fluorescence measurements obtained at different salt concentrations thus calculating the  Since both χpol and χmon are known thank to the above described procedure, it is straightforward to calculate the χdim as 1-( χmon + χpol). Even though the presence of NaCl is fundamental for the DX polymerization into gels, once they are melted there is no dependence of χdim on NaCl concentration (Figure 3 -panel A). At this point the composition of the investigated systems along temperature scan experiments is known and this allows the relevant dimerization constant of DX to be calculated according to where χdim and χmon are the molar fractions of dimer and monomer DX, respectively and [DX]anal is the analytical concentration of the drug. At [DX]anal= 1•10 -2 M used throughout this work, at variance with literature reports, [10] [12] even in the presence of NaCl there is no dependence of the dimerization constant on the ionic strength ( Figure 3 -Panel B). This last point requires a few words of comment. It is evident that the Kdim values calculated from DXpol samples lie on the same straight line of the sample without NaCl if two conditions are satisfied: i) the temperature is above Tpol (Tpol is the temperature at which the DX polymerization process starts upon cooling [9]); ii) the χdim values are not too small due to the high value of χpol (high uncertainty in the determination of the system composition due to the prevalence of DXpol). Due to these reasons, at high NaCl concentrations and at low temperatures, deviations from the Van't Hoff prediction were observed (Figure 3 - [9]).

Molecular Dynamics calculations
The ΔH dim° value was estimated using ΔH dim°= < H dim°> +< H mon°> , where <H dim°> and <H mon°> are the system energy averages in the presence of the dimer and of the monomer, respectively. Within the Van't Hoff approximation, i.e. the enthalpy does not change along the  isobar, our estimate at 450K can be compared with the experimental value in the explored temperature range. The association/dissociation kinetic constants, kass/kdiss, were calculated assuming a diffusive model for the monomer motion. The rate constant at the distance r0 where the monomers freely diffuse in the solvent is given by: [26] ass = 4π(D M,1 + D M,2 ) 0 where is the Avogadro number and DM,1 and DM,2 are the DXmon diffusion coefficients. From kass and the experimental value of ΔG dim°, kdiss was estimated by diss (T) = ass (T)e ΔG dim°( T) T The MD simulation at 300K was used to describe the conformational behaviour of the dimer, whereas the simulation at 450K to estimate the thermodynamics of the dimer-monomer equilibrium. The high temperature MD allows to sample a significant number of monomer/dimer transitions in a reasonable amount of time (200 ns), thus giving the possibility to estimate ΔH dim° and the associated statistical error. The dimer conformational behaviour at 300K was analysed in terms of the doxorubicin long axis as defined by the vector connecting the midpoints of the C2/C3 and C8/C9 bonds (see Chart 1). For each MD frame the scalar product between the long axes of the two molecules was calculated to obtain the main orientations between the two molecules, i.e. antiparallel (-1 to -0.5), orthogonal (-0.5 to 0.5) and parallel (0.5 to 1). The value of the monomer diffusion constant (DM=2.89•10 -11 m 2 /s) needed to estimate the kass was taken from the literature, [27] whereas the distance at which the monomers start to freely diffuse was calculated (as 1.5 nm) from the probability distribution of the distance between the monomer centres of mass as provided by the MD simulation at 450K (for further details, see [26]). The same 1.5 nm value was also used in the evaluation of ΔH dim° to define the monomeric state. As expected, the MD simulation at 300K shows the presence of the dimer only.
Using the long DX axis as a measure of the mutual orientation of two DX molecules and the boundary definitions of the three main conformational states (parallel, antiparallel and orthogonal), we found that the antiparallel arrangement is present in more than half of the trajectories (Figure 4). The parallel and the orthogonal orientations were observed in the 15% and 30% of the trajectories, respectively. The mean distance between the centre of mass of the two DX molecules in the dimer is 0.37 nm, in excellent agreement with the literature value of ∼0.34 nm based on NMR spectroscopy and molecular mechanics modelling. [14] Making reference to the mode of the distributions obtained by the MD simulations of Figure 4, for the antiparallel DX arrangement a value of -0.94 is obtained while, for the parallel orientation the Figure 4 The normalized histogram of the scalar product between the long DX axes (red arrow in Chart 1) as obtained from the 450 K MD simulation. Representative geometries for the antiparallel (-1 to -0.5), orthogonal (cross, between -0.5 and 0.5) and parallel (0.5 to 1) DX arrangements are included as stick representation.
mode is 0.93. These values allow to calculate the most probable angle between the long axes of two adjacent DX molecules. For the two orientations, angles of 160° and 21°, respectively, were obtained. These values are in excellent agreement with the literature data (168° and 15° for the antiparallel and parallel orientations, respectively). [14] From the MD simulation at 450K, the dimer-monomer thermodynamics was modelled. The value of ΔH dim° at 450 K obtained from the MD simulation (-45±3 kJ/mol), indicating a favourable interaction between the two molecules in water, well compares with our experimental measurements (Figure 3 -Panel B) obtained in the temperature interval 293K-328K.
The kinetic constants associated to the monomer-dimer equilibrium -estimated from the experimental value of Kdim and assuming a diffusive model when the distance between the two monomers is greater than 1.5 nm -are kass = 6.56•10 8 M -1 s -1 and kdiss = 1.338•10 4 s -1 yielding a dimer lifetime estimate of 74.7 μs (a direct estimate of the dimer lifetime from the MD simulation at 450K provides a value of ~ 11 ns, compatible with the 150 K temperature increase; see ESI - Figure S1). Finally, the analysis of the DX structures as sampled by the 450 K MD simulations indicates that its internal degrees of freedom are not significantly affected by the dimer formation, i.e. the distributions of selected dihedral angles of the monomeric form -in line with a previous literature report [28] -do not change upon dimerization (see ESI - Figure S2).

Conclusions
The intrinsic and peculiar fluorescent nature of DX has been exploited in this work to improve our knowledge of the drug association processes, even from a thermodynamic point of view. Also in this instance, the joint use of experimental techniques and MD simulations proved to be the winning strategy to unravel the complex association processes ruling the DX behaviour in solution. Despite its low molar mass, in fact, DX has a chemical complexity that is at the basis of its capability not only to form dimers but also to undergo a T-controlled polymerization that leads to the formation of long chains of DX molecules responsible of the gel properties of the system. [9] The interesting point of the DX dimerization process is that, when the DX concentration is above the CGFC, a substantial independency on ionic strength was found ( Figure 3). The all-atom MD simulations here presented, show mutual orientations of the DX molecules in the dimer closely matching those reported in the literature obtained using ad-hoc restrained MD simulations. [14] In addition, theoretical-computational modelling not only lends support to the experimental findings but also offers an estimate of the residence time of DX molecules in the dimer hardly achievable experimentally even with T-jump apparatus of the last generation.
The thermodynamic parameters derived for the dimerization process both from the fluorescence measurements and from MD calculations here presented are in fairly good agreement with literature data. [14][29] However, at variance with literature data, [10] [12] our results point to an independence of the dimerization constant on the NaCl concentration ( Figure 3 -Panel B). This point will be the subject of a forthcoming paper, where the dimerization process of DX will be addressed to by means of both fluorescence and MD A C C E P T E D M A N U S C R I P T simulations by varying in a systematic way the drug concentration below 1•10 -2 M and the NaCl concentration. The resulting systems will be studied as a function of temperature in order to have access to the thermodynamics of the dimerization process as a function of ionic strength.