Structure of supramolecular astaxanthin aggregates revealed by molecular dynamics and electronic circular dichroism spectroscopy

Grzegorz Zajac,1 Ewa Machalska,1 Agnieszka Kaczor,1,2 Jiří Kessler,3 Petr Bouř, *,3 Malgorzata Baranska*,1,2 1 Faculty of Chemistry, Jagiellonian University, Gronostajowa 2, Krakow 30-387, Poland 2 Jagiellonian Centre for Experimental Therapeutics (JCET), Jagiellonian University, Bobrzynskiego 14, Krakow 30-348, Poland 3 Institute of Organic Chemistry and Biochemistry, Academy of Sciences, Flemingovo náměstí 2, Prague, 16610, Czech Republic


Introduction
We concentrate on the astaxanthin dye (3,3 0 -dihydroxy-b,bcarotene-4,4 0 -dione, AXT, Fig. 1), because it is a convenient system to study, important in biology and potentially useful in human medicine.2][3][4][5][6][7][8] Humans and animals cannot produce it and because AXT exhibits similar metabolic and physiological functions to b-carotene, zeaxanthin and lutein 5 it must be supplied by food. 2,4,9XT chemistry is important for biological functions.1][12] The relatively long polyene chain is responsible for strong absorption of green and blue light, and it also largely determines the aggregation properties via p-p stacking interactions.][15][16][17][18] The conjugated bonds can adopt both the trans and cis conformations, 19,20 although the trans-form prevails. 14,21Each molecule has two chiral centres at the C3 and C3 0 atoms giving rise to two enantiomers (3R,3 0 R) and (3S,3 0 S) and one meso form (3R,3 0 S). 17 Six conformers are generated due to mutual orientation of the ionone rings (gauche and trans conformations). 19,22hese variations of structure influence the bioavailability of AXT. 16,20,21,23][26][27][28][29][30][31][32][33][34][35] Typically, aggregation occurs when water is added to an AXT solution in an organic solvent. 36,37Aggregation almost always results in significant spectroscopic changes, often sensitive to aggregate type.AXT forms two types: H-aggregates, characterized by a blue-shifted strong electronic absorption band around 400 nm, and J-aggregates exhibiting a red-shift.The H-aggregates are characterised by tight ''card-pack'' stacking where the polyene chains are more or less parallel to each other, while J-aggregates with ''head-to-tail'' or ''herring bone'' orientation of the conjugated chains are looser. 26,33,34,38,39Formation of a particular assembly depends on the balance of intermolecular forces between molecular constituents, such as hydrogen bonds, electrostatic forces, and van der Waals interactions. 30,34,36In spite of gross phenomenological models of the H and J aggregates, however, detailed mutual arrangements of the monomers in them is not known.
However, while (non-resonance) ROA intensities can be computed relatively easily, a straightforward simulation methodology for AIRROA does not exist.Simulations are needed in structure-function studies of biomolecules, to provide a link between the spectral shapes and system geometry.As a first step, in the present study, we concentrate on the dependence of the CD spectra on the aggregation, because CD intensities were shown to be directly related to the resonance ROA effects. 42n particular, according to a one electronic resonant state model, the ROA sign is determined by the CD intensity at the excitation wavelength.
Although we currently cannot model larger-scale ordering of the aggregates, molecular dynamics simulations on limited models (the AXT dimer and decamer) are used to understand formation of different supramolecular structures in different solvent mixtures, and the results are compared to experimental data.The molecular dynamics methodology has been suggested for AXT also in the past, 43 but according to our current knowledge, the dependence of the aggregation geometry and its spectral features on the solvent type has not been theoretically studied yet.][49][50][51][52] The results indicate that the combination of experimental and theoretical approaches provides a solid basis for interpretation of the experimental data, including information about the structure and chirality of the aggregates, and the aggregation mechanism.

Experimental
The aggregate preparation and experimental spectra were reported previously. 25,26Briefly, the aggregates were prepared by mixing 0.3 mL 1.0 Â 10 À4 M AXT acetone stock solution with 2.7 mL distilled water (H-aggregate), or with 0.6 mL acetone and 2.1 mL distilled water (J-aggregate), thereby obtaining solutions with 1.0 Â 10 À5 M concentration.A monomer AXT reference solution was prepared by mixing 0.3 mL stock solution with 2.7 mL acetone.
Electronic absorption (EA) and electronic circular dichroism spectra were measured in a quartz cell (optical path length 1 cm), in the 330-900 nm spectral range, using Jasco J-815 and UV-vis PerkinElmer Lambda-35 spectrometers, respectively.

Arbitrary models
Arbitrary models of H and J-aggregates of (3S,3 0 S)-AXT were created from 10 and 8 molecules in the most stable conformation (-GG). 16,19,41The monomer was optimized at the CAM-B3LYP/ Def2SVP theoretical level, using the Gaussian 09 E.01 package. 535][56] The monomers in the H-aggregate have been arranged into a helical, parallel orientation, with 5 Å distance and 201 angle between the polyene chains of neighbouring molecules (Fig. 2).For the J-aggregation the monomers were arranged into two twisted layers (401) and separated by 25 Å (Fig. S1, ESI †).

Molecular dynamics
Two solutions of acetone and water mixed in 1 : 9 and 3 : 7 ratios were selected based on experimental knowledge, 25,26,57 as these mixtures favour formation of the H and J-aggregates of AXT, respectively.Two models of both the H and J aggregates were investigated (dimers and decamers), differing by the numbers of included AXT molecules.To compare obtained electronic spectra of aggregates with non-aggregated AXT, we performed a monomer simulation as well, where a single AXT molecule was placed in acetone.MD simulations were performed using the Amber 14, 58 Packmol 59,60 and VMD 61 programs.The general amber force field (GAFF) has been used.To model the monomer, one AXT molecule was placed into the centre of a cubic box (40 Â 40 Â 40 Å 3 ) and the remaining space was filled with 511 acetone molecules.In the dimer models, two AXT molecules parallel to each other separated by 20 Å were put in the same box and the remaining space was filled with a mixture of 50 : 1869 and 150 : 1454 acetone : water molecules, corresponding to the 1 : 9 or 3 : 7 volume ratio, respectively, needed for H and J-aggregation.Finally, the decamer models consisted of ten AXT molecules randomly distributed in a larger box (80 Â 80 Â 80 Å 3 ), and additional 407 : 15 128 and 1221 : 11 767 acetone : water molecules (for the 1 : 9 or 3 : 7 volume ratios).
All the simulations have been performed using nVT ensembles (T = 300 K), Langevin dynamics with the collision frequency 2 ps À1 , a non-covalent interaction cutoff at 8 Å, the particle mesh Ewald (PME) method, and a 1 fs integration time step.To account for the electron conjugation presumably not well represented in GAFF, CQC-CQC torsion angles in the polyene chain (except for C5-C6-C7-C8 and C5 0 -C6 0 -C7 0 -C8 0 ) were restrained to 1801, using a harmonic force constant of 1.8 Â 10 À2 kcal mol À1 rad 2 .After 0.5 ns equilibration, production runs of 20 ns (the monomer), 205 ns (dimers) and 110 ns (decamers) were performed.Geometry frames were saved every 10 ps.From hundred monomer frames, evenly distributed in time within the 20 ns simulation, the solvent molecules were deleted and CD and EA spectra were simulated.The dimers and decamers were allowed to aggregate for 5 ns and 30 ns, respectively.Thousand dimer frames evenly distributed over the remaining 200 ns and 400 decamer frames evenly distributed over the remaining 80 ns were used to simulate the spectra, as for the monomer.The simulation conditions are summarized in Table 1.

Electronic absorption and electronic circular dichroism
Transition energies, oscillator and rotation strengths for the first 100 (arbitrary models), 40, 80 and 100 (monomer, dimer and decamer) singlet electronic transitions were calculated using the ZIndo/S method.For the monomers and dimers analogous calculations were performed using TD-DFT, with the CAM-B3LYP functional and the Def2SVP basis set.The EA and ECD curves were obtained by convolution with Gaussian functions of 0.2 eV half width at full maximum, using the CD Spec Tech program 62,63 and averaged over the MD snapshots.

Main features in experimental EA and ECD of AXT aggregates
The electronic absorption is dominated by the strong transition between the S 0 (1 1 Ag À ) -S 2 (1 1 Bu + ) states coming from the polyene chain chromophore. 64The S 0 -S 2 excitation does not exhibit any resolvable vibrational structure, which is sometimes observed for carotenoids with conjugated carbonyl groups.The H-and J-aggregates are characterized by the blue and red shift of the main electronic absorption band and the appearance of corresponding strong Cotton effects in CD.In the case of the J-aggregate in acetone-water mixed solutions, a red shift (from 478 to 515 and 551 nm) of the S 0 -S 2 excitation is observed.Note that two major types of AXT H-aggregates have been discussed in the literature, 25,26 one with a blue shift of the S 0 -S 2 excitation from 492 to 386 nm in DMSO-water mixed solutions (H1-aggregate), and one with a minor blue shift from 478 to 469 nm in acetone-water mixed solutions (H2-aggregate).The monomers do not possess prominent rotatory strength at this transition.Supposedly, prepared samples with prevalent H-and J-aggregates are always mixtures of more aggregate subtypes and monomers.

EA and ECD of arbitrary models
The theoretical spectra of the arbitrary aggregate models are compared to experimental data in Fig. 3.One can see that the positions of the main absorption bands both for the H-and J-aggregates are very close to the experimental ones.For J the main experimental band exhibits a substructure that is not reproduced.The calculated characteristic blue and red shifts of this transition for the H-and J-aggregates with respect to the monomer can be seen in Fig. S3 in the ESI, † and agree with the observations.ECD of both models exhibits Cotton effects (positive/negative at longer/shorter wavelengths) in the S 0 -S 2 excitation region, similar to the experimental spectra, however theoretical ECD couplets (zero-line crossings) are shifted to longer wavelengths, more than the absorption band, which we attribute to an effect of inhomogeneous line broadening and splitting caused by longer-range coupling, impossible to simulate at this stage.

MD simulations
The aggregation process and basic structural information (parallel or head to tail organization) of the four MD models are apparent from the AXT-AXT distances and the C6 AXT1 -C6 0 AXT1 -C6 0 AXT2 -C6 AXT2 torsional angles.Their histograms and dependencies on the simulation time are plotted in Fig. 4 and ESI † (Fig. S4-S7).
In the case of the 1 : 9-dimer system (Fig. S4(a), blue curve, ESI †), the aggregation starts almost immediately (o3 ns) when the initial distance of 20 Å shortens to about 5 Å, more precisely it fluctuates between 3 and 10 Å during the remaining 202 ns.This corresponds to a tightly organized structure of the H-aggregate.The parallel orientation of aggregated AXT is confirmed by the AXT-AXT dihedral angle histogram (Fig. 4(a), blue).The highest probability is observed between À301 and 301, with two maxima at 131 and À131.
In the 3 : 7-dimer model the molecules move more freely as can be seen in Fig. S4(a) in the ESI † (red curve).The AXT-AXT distance varies between 3 and 25 Å during the entire 205 ns simulation.The looser structure of the 3 : 7-dimer is also confirmed by the dihedral angle histogram (Fig. 4(a), red).Contrary to the 1 : 9-dimer, all angles can be adopted, albeit with a low probability density.Probability maxima (À181 and 141), however, are nearly at the same angles as for the 1 : 9-dimer.In terms of partial integral probabilities, 8.3% of AXT pairs in the 1 : 9-dimer and 59.5% of pairs in the 3 : 7-dimer adopt the dihedral angle between À1801 and À301 or 301 and 1801 characteristic of weaker binding.
For the 1 : 9-decamer the starting averaged distance between all AXT molecules of approximately 33 Å shrank to about 15 Å after 20 ns; it oscillates between 12 and 18 Å, thus indicating a tightly packed structure and similar behaviour to the 1 : 9-dimer model.However, the aggregation time is considerably longer (Fig. S6(a), blue curve, ESI †).The angular distribution is broader, but the two maxima (À181 and 171) corresponding to the two preferred stacking angles are the same (Fig. 4(b), blue).
The simulation thus agrees well with the observations, indicating that the 1 : 9 solvent mixture supports tight and nearly parallel organization of AXT molecules, with two mostpreferred values of the dihedral angle.On the other hand, the 3 : 7 mixture leads to a more weakly associated supramolecular structure, allowing considerable freedom to AXT-AXT orientations and distances.

EA and ECD of MD clusters
A few examples of typical MD clusters of AXT dimers and decamers (without solvent molecules) are shown in Fig. 5 and 6.In almost all extracted 1 : 9-dimer clusters, the AXT molecules are oriented in parallel to each other and tightly packed.About three aggregate types can be recognized (Fig. 5(a-c)), characterized by  different values of the dihedral angles, distances, and ring orientations.In the a and c clusters, the terminal rings are connected by the hydrogen bonding between the hydroxyl and carbonyl group and bridging water molecules (Fig. 7).
The 3 : 7-dimers, contrary to 1 : 9, are less associated, and prefer different orientations.Four major types of the 3 : 7-dimers are exemplified as the serial dimer (Fig. 5d), X-dimer (Fig. 5e), T-dimer (Fig. 5f) and L-dimer (Fig. 5g).Parallel dimers similar to these found for the 1 : 9 system are also encountered, but they are less frequent.UV-vis and ECD spectra of selected structures are presented together with experimental spectra in Fig. S8 and S9, ESI.† Some of the structures reproduce well the experimental shapes (1 : 9-c and 3 : 7-d,e,f), but others rather give opposite signs (1 : 9-b and 3 : 7-g) or the wrong intensity (1 : 9-a).It is worth noting that the presented ''typical'' structures cannot substitute for the variety of all dimers in the MD clusters.Spectral averaging over typical structures would be problematic due to the large degree of freedom in the real solution. 65he decamers show similar behaviour to the dimers.The computed geometries are consistent with structural models proposed previously. 43(Note that a helical longer-range arrangement was also proposed for carotenoid aggregates, 34 this however, is not possible to simulate reliably with our computational means).
The simulation thus provides useful insight into the underlying processes of aggregation.It confirms that the low concentration of acetone (in the 1 : 9 model) is favourable for hydrophobic attraction of AXT molecules, although all weak interactions (electrostatic, p-p stacking, van der Waals forces, and hydrogen bonding) participate.At higher acetone concentrations the hydrophobic attraction is weaker.A similar situation was described for AXT aggregation in mixed water-ethanol solutions. 43XT dimers.The computed UV-vis and ECD spectra of the dimer MD models, obtained as averages over 1000 clusters, are presented together with the experimental spectra in Fig. 8.The convergence of selected ECD intensities during the MD simulation can be seen in Fig. 9 and in the ESI † (Fig. S10 and S11).We can see that the spectra converge more slowly for the ''3 : 7'' system, nevertheless they stabilize at about the 700 cluster average.Note that the equilibration required rather long MD simulation times, to allow an average of a sufficient number of independent (uncorrelated) geometry frames, i.e., separated by long time intervals.
For both dimers, the maximum of the main UV-vis band is predicted at shorter wavelengths compared to the experimental spectra.For the 3 : 7-dimer model this shift is considerably more prominent.Both dimer models and both theoretical methods lead to the same positive/negative (at high/low wavenumber) ECD couplet, which is also observed in the experimental spectra of AXT aggregates.This couplet results from splitting of the AXT excited energy levels during the interactions.The sign of the couplet corresponds to the ''positive'' chirality in the dimers, i.e. to the slight prevalence of the positive torsional angles as indicated by the histograms, both for the 1 : 9 and 3 : 7-dimers (Fig. 4(a)).
An alternate representation of the results may be a comparison of the monomer, J and H aggregate EA and ECD spectra, as done for the experiment and computation in the ESI: † Fig. S2 (experiment), S15 and S16 (ZIndo/TDDFT theories).The detailed line structure of the simulated spectra can also be seen in the ESI † (Fig. S18-S21).Some experimental trends are not well-reproduced, such as the blue and red shifts of the absorption maxima of the H and J aggregates, if compared to the monomer.The theoretical maxima remain almost constant.The aggregation also broadens the experimental signal, leading to a distinct signal at 551 nm for the J aggregate attributed to its two different forms, 57,66 which is not seen in theory (see Fig. S2,  ESI †).Supposedly the dimer model is too crude to reproduce all  aspects of the structure of the loose J-aggregation; the H-aggregate experiment is reproduced much better, which is consistent with the ''card-packing'' concept where the shortrange interactions reasonably-well reproduced by the model dominate.
The TDDFT and ZIndo/S methods gave quite similar results in the region of the main absorption band, while the calculation times were significantly (B50 times) shorter for the former one.TDDFT gives only slightly higher transition energies and ECD intensities.
AXT decamers.The theoretical decamer results can be seen in Fig. 10.The spectra calculations had to be restricted to the ZIndo/S method; the TD-DFT approach for such big systems (400 clusters times 960 atoms for each decamer model!) is not feasible with the available computational infrastructure.The convergence of ECD intensities with respect to the number of clusters (Fig. S22, ESI †) is similar to that for the dimers.We suppose that neither can the decamer model reproduce some features of the J-aggregate geometry, which could explain the too narrow simulated absorption band.

Fig. 2
Fig.2Arbitrary models of H and J-aggregates.

Fig. 3
Fig. 3 Calculated and experimental EA and ECD spectra of the H and J-aggregates' arbitrary models.

Fig. 7
Fig. 7 Solvent distribution and hydrogen bond network around the AXT aggregate geometry, taken from a single MD snapshot of the 1 : 9-dimer model, and the radial distribution function (RDF) of oxygen atoms from water (blue) or acetone (green) around oxygen atoms from AXT hydroxyl groups (solid lines) or AXT carbons C15 (dotted lines).

Fig. 8
Fig. 8 Calculated EA and ECD spectra of the 1 : 9 and 3 : 7 AXT dimers, averaged over 1000 MD clusters, combined with the experimental data for H and J-aggregates of AXT.

Fig. 10
Fig. 10 Calculated EA and ECD spectra of 1 : 9 and 3 : 7 AXT decamers, averaged over 400 MD clusters, combined with the experimental data of H-and J-aggregates of AXT.

Table 1
Systems investigated in molecular dynamics simulations