Molecular simulation studies of cyanine-based chromonic mesogens: spontaneous symmetry breaking to form chiral aggregates and the formation of a novel lamellar structure

All-atom molecular dynamics simulations are performed on two chromonic mesogens in aqueous solution: 5,5’-dimethoxy-bis-(3,3’-di-sulphopropyl)-thiacyanine triethylammonuim salt (Dye A) and 5,5’-dichloro-bis-(3,3’-di-sulphopropyl)-thiacyanine triethylammonuim salt (Dye B). Simulations demonstrate the formation of self-assembled chromonic aggregates with an interlayer distance of ~0.35 nm, with neighboring molecules showing a predominantly head-to-tail anti-parallel stacking arrangement to minimize electrostatic repulsion between hydrophilic groups. Strong overlap of the aromatic rings occurs within the self-assembled columns, characteristic of H-aggregation in aqueous solution. At low concentrations, aggregates of Dye A form chiral columns, despite the presence of strictly achiral species. Chirality arises out of the minimization of steric repulsion between methoxy groups, which would otherwise disrupt the stacking of aromatic molecular cores. At higher concentrations, simulations suggest the interaction of short columns leads to the formation of an achiral-layered structure in which hydrophobic aromatic regions of the molecule are sandwiched between two layers of hydrophilic groups. This novel lamellar structure is suggested as a likely candidate for the structure of a J-aggregate. The latter are known to exhibit intense red-shifted absorption peaks in solution but their structure has not yet been characterized. Self-organization of such structures provides a route to the formation of “smectic” chromonic mesophases. demonstrate the formation of chiral chromonic stacks in dilute aqueous solution. Interaction of chromonic stacks at higher concentration leads to the formation of a novel layered structure, which is suggested as a candidate for a smectic chromonic phase.


Introduction
Chromonic liquid crystals are an unusual class of lyotropic liquid crystals (LLC), which have recently generated considerable interest. [1][2][3][4] Chromonic phases form from disc-like molecules in a solvent, usually water, through a two-stage process: involving initial formation of aggregates (stacks) and subsequent alignment of aggregates to form liquid crystal phases. [5,6] Unlike conventional lyotropics, chromonics can form ordered phases at extremely low concentrations. [4,7,8] Moreover, aggregation into stacks takes place without a critical micelle concentration, forming aggregates of variable length. Subtle changes in the balance of enthalpic and entropic (hydrophobic) interactions can easily enhance or destabilize chromonic phase behavior. [9,10] Chromonic molecules tend to have rigid disc-like aromatic cores, decorated with hydrophilic groups on the periphery, to provide water solubility. Although once thought to be rare, it is now understood that chromonic self-assembly can occur in a wide variety of molecules: including drugs [11,12] , dyes [13,14] and nucleic acids. [15] The renewed interest in chromonics arises from a range of interesting potential applications, which are enhanced by tunable aggregation and phase behavior, together with the ease of alignment of chromonic phases (provided by magnetic fields, mechanical shearing, or surface alignment layers). These applications have included controllable selfassembly of gold nanorods, [16] fabrication of highly ordered thin films, [17] formation of biosensors, [18,19] and novel experiments on the interplay between hydrodynamics and topology of active matter. [1] Chromonics have the possibility of forming a range of complex self-assembled supramolecular aggregates, which can be potentially converted into aligned thin films. Two stacking patterns predominate: H-aggregates, where molecular cores stack directly on top of each other, and J-aggregates where molecular cores adopt a staggered conformation. [7,20,21] However, within these two aggregation motifs further complexity can arise. Hence, there are intriguing literature suggestions for alternating head-to-tail assembly, [22] double columns [23] and longer ranged ordered structures, such as the hollow chimney and brickwork models ( Figure S1). [8] Consequently, self-assembled cylindrical aggregates of carbocyanine dyes have been suggested as an interesting candidate for synthetic light harvesting systems and electronic energy transport wires. [24] Atomistic simulation provides a powerful way of understanding chromonic aggregation at a structural and thermodynamic level. [9,22] For example, for the azo dye sunset yellow, [22] atomistic simulations have been used to distinguish between competing models of self-assembly, have explained many of the experimental observations of chromonics and have shown the likely structure of the chromonic nematic phase.
In this paper we simulate the behavior of two cyanine dye molecules, which demonstrate chromonic self-assembly: 5,5'-dimethoxy-bis-(3,3'-di-sulphopropyl)-thiacyanine triethylammonuim salt (Dye A) and 5,5'-dichloro-bis-(3,3'-di-sulphopropyl)-thiacyanine triethylammonuim salt (Dye B) ( Figure 1). The two mesogens share an identical molecular core, but differ in the periphery groups attached to the aromatic rings (the 4-position of the external phenyl ring). Cyanine dyes of this type have been studied extensively experimentally [7,8,25] as they undergo dramatic changes in optical properties on aggregation and have important applications including spectral sensitization of photographic materials. [26] For such systems, there are many competing ideas of how self-assembly takes place and what the structure of aggregates and mesophases may be (see Figure S1). We test these directly by means of large-scale molecular simulation. We predict a new structure for H-aggregates, including the possibility of chiral columns forming from achiral molecules and we suggest a structure for J-aggregates in solution that allows the formation of a "smectic" chromonic phase at low molecular concentrations. In the discussion, we show how simulations help rationalize the differing observations of self-assembly and mesophase formation seen for these systems.

Self-assembly in dilute solution
Starting from a random arrangement of dyes and counterions, both dyes show spontaneous selfassembly in solution. At low concentrations, dimer aggregates form within ~1 ns for both dyes, and continued simulations show the formation of two stacks in solution over a 200 ns timespan. Locally, association occurs through face-to-face interactions, with the aromatic parts of the molecule stacking on top of each other with a mean spacing of ~0.35 nm (Figure 2), which can typically be seen in WAXS studies of many chromonic systems. The spread of distances in figure 2 shows the liquid-crystallinenature of the columns with molecular centers fluctuating with respect to each other as a function of time.
There is a very strong preference for antiparallel stacking, as seen in Figure 3. This allows the charged sulfonate groups from a dye molecule to lie on the opposite side of the stack to those from the molecule below, and allows them to be well-solvated by water molecules, away from the aromatic core.
The aggregation behavior of dye B has been studied experimentally alongside several other cyanine dyes. [26,8,7,25,27,28] At extremely low concentrations, 6 1 10   M, [26] the dye exists as a monomer in solution with an absorbance maximum at 428 nm. However, increasing the concentration leads to the growth of a more intense, and wider, absorption peak at lower wavelength, 406 nm, which overlaps the monomer absorption signal. [7,26] The solution remains isotropic in this regime. This hypsochromic shift (blue shift) is believed to be characteristic of the growth of H-aggregates, with face-to-face overlap of molecules in dimers and possibly larger molecular aggregates. At slightly higher concentrations ~4 1.4 10   M, an additional absorption is evident, corresponding to an intense narrow J-band, which occurs at longer wavelengths. The latter is believed to indicate the presence of J-aggregates. According to Moll et al., [26] a first order transition to a mesophase containing J-aggregates, occurs via a heterogeneous two-phase regime consisting of dye solution and dye liquid crystal. During this regime, there is an equilibrium between monomer/H-aggregates/J-aggregates, M H J , which shifts towards J-aggregation with increasing concentration.
In Figure 3, we see H-aggregate behaviour, due to the strong overlap of aromatic parts of our molecule within an aggregate. Simple ZINDO calculations, carried out on snaps hots of two and three molecule aggregates taken from the molecular dynamics runs, confirm a blue shift in the UV-Vis spectra in moving from a monomer to a dimer and trimer. These H-aggregates are similar to those expected in several other chromonic systems, [13,29] where there is a strong enthalpic driving force to "hide" the aromatic rings from water. In our systems, the driving force is sufficient for parallel aggregation to occur in addition to the preferred anti-parallel association. Within multiple repeated simulations of the self-assembly process, we observed also several cases where parallel stacking occurred (including some continuous stacks of up to five parallel molecules). In most cases, this generates strained-bent aggregates that eventually fragment (see Figure S2 and Figure S3), allowing future recombination in a preferred anti-parallel stacking arrangement.
To quantify the strength of aggregation, we carried out potential of mean force, PMF, calculations for Dye A, as shown in Figure 4. From the depth of the well, the binding energy is estimated to be ~35 kJ mol -1 (~14 B kT). This is towards the higher end of values found for chromonic systems by extraction from experiments (TP6EO2M: [30] 14 B kT, Sunset Yellow: [13,31] ~7 B kT, disodium cromoglycate: [32] 7 B kT, Blue 27: [32] 12 B kT, Benzopurpurin 4B: [33] 10.3 B kT) or from molecular simulation (TP6EO2M: [9] 14.7 B kT, Sunset Yellow: [22] 7 B kT), but consistent with the formation of Haggregates at very low concentrations. [7] For the two dyes, we calculated solvent accessible computational cross sectional areas (CSA) for the aromatic regions of the aggregates. These are given in Table 1 and correspond to typical singlemolecule cross-sections. Hence, the unusual packing arrangement seen cannot be used to explain the large cross-sections deduced for some other chromonic systems. However, calculating the solvent accessible surface around a number of long stacks of Dye A (using a solvent probe radius of 0.14 nm, Figure S4) and calculating the internal volume within, gives a volume of 30 12457 10   m 3 with an error of <10%, which corresponds to a density within the stack of 1.37 g cm 3 , rather higher than that typically used to calculate cross-sectional areas from X-ray data.

Formation of chiral aggregates
To further test the observed preference for anti-parallel association, we pre-prepared columns of 30 molecules in both anti-parallel and parallel configurations for both dyes, with random arrangements of counterions, solvated to give 9 wt% solutions. For the parallel cases, we see fragmentation of columns and eventual growth of anti-parallel order ( Figure S3). For dye A, starting from an anti-parallel initial configuration, the system anneals over ~200 ns to give a chiral column ( Figure 3). Chirality appears to arise through the system finding the most efficient packing arrangement for the bulky methoxy group.
Adoption of a chiral column allows the achiral dye molecules to minimize the strain caused by the methoxy group wanting to be further apart than the optimal 0.35 nm spacing for aromatic rings. We note that similar chiral packing arrangements have been found in coarse-grained models for achiral discshaped molecules. [34,35] In the absence of a methoxy group, chirality is not seen for Dye B. Instead, slow fragmentation of the long pre-assembled columns occurs over hundreds of nanoseconds, (see below).
For Dye A, we calculate a 2d density map (  (Figure S5). At this concentration the stack starts to interact within itself via the periodic boundary conditions to form a layer (see below), and the twist is unwound to accommodate in-layer dye interactions with the charged arms expelled into the water above and below the layer. For Dye B, in the absence of bulky methoxy groups, nearest neighbor molecules favor a smaller twist of ~160° ( Figure S6), next nearest neighbors favor a twist of ~20°, and the long-range column twist seen for Dye A is completely absent.

Formation of layered structures
For Dye B at 30 wt%, fragmentation of an anti-parallel seeded column occurs over 200 ns to give smaller aggregates with a range of aggregation numbers from 2 upwards (see Figure S7). At this concentration, aggregates start to interact strongly within the space confines of the periodic box. To look at this further we simulated a system of 100 molecules, starting from a random initial configuration at this concentration. The results are shown in Figure 6.
Here, we see a layer structure forming, from the coming together of small aggregates. The layer is one molecule thick and the aggregates organize to allow sulfonate groups to lie above and below the plane of the layer. The time scales for this process are extremely slow, and there are insufficient molecules to form a full layer across the entire simulation box, or to form a multilayer structure (hence we cannot guarantee that we see the full equilibrium structure of a smectic chromonic). However, the structure seen in Figure 6 matches experimental observation quite well. [8,25] X-ray results suggests that the layers are one molecule thick, [25] that the molecular long axis lies in the plane of the layer, [8] that configurations occur where neighboring molecules are shifted in terms of the overlap of aromatic ring structures (J-aggregates) and the layers formed are uniaxial (not biaxial). [25] This also shows a mechanism for the conversion of H-aggregates to J-aggregates, which occurs continually as concentration is increased across the two-phase region, where dye mesophase and dye solution are in dynamic equilibrium. Hence, in many ways this provides a better picture than the previously suggested biaxial brickwork arrangement of dyes [8,25] (Figure S1).
For Dye A, the binding of the H-aggregates appears to be stronger than for Dye B at 30 wt%. We see columns coming together to form a layered structure (Figure 7) but, on the (200 ns) time-scale of the simulation, we do not observe the fragmentation into short columns seen with Dye B. We note that for a structurally related compound, pinacyanol chloride, Berlepsch et al., [36] have reported H-aggregates that transform into J-aggregates on the time-scale of weeks for a macroscopic sample. So the time-scale to see these changes, even at the nanoscale, is expected to be long.
We note in passing that throughout this study we found no evidence for the stability of a double width column that has been suggested as a possible chromonic dye structure [23] and was hypothesized as a possible structure for dye A by Bottrill. [28] To test this hypothesis we created a pre-seeded assembly of 32 dye A molecules (as in Figure S1) with each 'layer' composed of two dye molecules with the flexible arms pointing away from the center of the column, with the next layer rotated through 90° to avoid unfavorable interactions with the sulfonate groups of the neighboring layer. This was then solvated with TIP3P water and randomly placed counterions to give a single double stack at 26 wt%. With an initial aromatic-ring separation of 0.4 nm the structure is free from both unfavorable steric and electrostatic interactions but, at 326 K, spontaneously fragments over a few ns of simulation time and then reforms a single width column over a period of 200 ns (as shown in Figure S8).

Discussion
It is clear from both experiments, and the current simulation study, that self-assembly and mesophase formation can be very complex for cyanine dye systems. There are competing modes of association: Haggregates and J-aggregates. [7,36] The latter can lead to the formation of smectic chromonic mesophases through the formation of semi-infinite sheets of dyes, which can be in equilibrium with a dye solution containing other aggregates. There is also potential for aggregation to change, due to small changes in chemical structure, or due to changes in pH or electrolyte concentration, through changes in counterion, [36] or simply through the presence of impurities. There is also strong evidence that, for some dyes, sheets of J-aggregates can curl up to form water-containing cylinders, which can in turn selforganize into a nematic chromonic phase. [24,36,37] This also helps to explain why the smectic chromonic phase nucleates from the isotropic phase, rather than a nematic phase. The long time scales seen in terms of some experimentally observed structural changes [36] have contributed to making molecular order in these systems challenging to understand.
The atomistic simulation approach in the current study provides useful pointers to the structure of H-aggregates that are seen to self-assemble spontaneously within a few hundred nanoseconds. They provide also clues as to possible structures of more complex J-aggregates. We note that the time scale and potential length scales associated with mesophase formation remains very challenging to study at an atomistic level. Hence, the possible key to a fuller understanding of these complex systems may be through development of new coarse-grained models that can be studied on time and length scales that are at least three orders of magnitude longer than atomistic simulation can provide currently. Here, it will be necessary to capture molecular structure and dimensions correctly, together with the correct free energy of association (as in Figure 3). Candidates for these models include recent work on SAFT- and MARTINI models for chromonics, [10] or dissipative particle dynamics (DPD) models, [38,39] which may allow the study of competing modes of aggregation, counterion effects and hierarchical self-assembly.

Conclusion
Atomistic molecular dynamics studies, for two cyanine dyes in aqueous solution have provided, for the first time, molecular simulation insights into chromonic aggregation in solution for these complex dye systems. We see the spontaneous formation of H-aggregates, which show a ~0.35 nm stacking of aromatic discs and a strong preference for a head-to-tail stacking arrangement; thus minimizing electrostatic repulsion between charged sulfonate moieties on neighboring dye molecules. Spontaneous chirality is observed in aggregates for Dye A at low concentrations despite the presence of strictly achiral species. Chirality arises from the need to minimize steric repulsion from the methoxy groups, while maintaining efficient - stacking.
At high concentrations, we observe the self-assembly of molecular monolayers. This process is most easily seen in Dye B, where the pair binding energy is smaller than for Dye A. Here, we see fragmentation of H-aggregates and formation of a novel one-molecule thick uniaxial-layered structure in which the sulfonate groups point above and below the plane of the aromatic layer. These non-chiral planar monolayers are candidate structures for J-aggregates, which are known to occur from experimental absorption measurements and X-ray diffraction studies of aqueous solutions of cyanine dyes. The self-organization of these structures provides an easy mechanism for the formation of a "smectic" chromonic phase.

Atomistic Simulations
The atomistic simulations were carried out using molecular dynamics and a classical force field. The simulations used the GROMACS 4.5.5 package [40] and make use of the GAFF force field parameter set. [41] Here, the interaction potential is given by eq eq bonds angles dihedrals impropers  , are respectively natural bond lengths and angles,  are dihedral angles,  and d  are phases angles, r K , K  , n V and d k are respectively bond, angle, and torsional force constants, ij  and ij  are the usual Lennard-Jones parameters, and i q , j q are partial electronic charges. The standard Lorentz-Berthelot mixing rules have been applied throughout this work. The Antechamber software from AmberTools was used to generate GAFF topologies, with partial atomic charges calculated from BCC partitioning [42,43] with the AMBER sqm program. [44] Charge is delocalized charge across the center linking group of the dye molecules, and hence AmberTools fails to assign good torsional angle potentials for the central fragment. Consequently, torsional energies for the central fragment were calculated using DFT with the hybrid density functional B3LYP [45] together with a 6-31G** basis set, [46] employing the Gaussian 03 [47] program. The GAFF parameters were then modified accordingly to match the DFT energies. This leads to a planar configuration for the molecular core for the dye molecules, as drawn in Figure 1 and a very high barrier height to rotation of 70 kJ mol -1 . The TIP3P water model [48] was chosen, as most applicable to use with the GAFF force field. The GAFF topologies and coordinate files were converted into the GROMACS format using the ACPYPE script. [49] Simulations were initially equilibrated within the canonical (constant-NVT) ensemble, followed by extensive runs in the isobaric-isothermal (constant-NPT) ensemble. All simulations, unless otherwise stated, were performed at a temperature of 300 K, with a Nosé-Hoover thermostat [50,51] being employed.
Pressure was controlled with a Parrinello-Rahman barostat [52,53] at atmospheric pressure and an isotropic pressure coupling. Bond constraints were applied using the LINCS algorithm [54] with a 2 fs time step.
Interaction cut-offs were applied for Lennard-Jones (1.1 nm) and Coulombic interactions (1.2 nm). The long range part of the Coulomb potential was accounted for by employing a Particle Mesh Ewald (PME) summation. [55,56] .

System details
Initial simulations consisted of a system of 30 dye molecules and 30 counterion molecules randoml y oriented and positioned in a cubic box and solvated to give a 9 wt% dye concentration. The same process was applied for concentrations of 20 and 30 wt%. For the 30 wt% solutions, additional larger simulations were performed consisting of 100 dye-counterion pairs in solution. Additional simulations, starting from seeded structures, were carried out as described in the results section. After equilibration, production runs were carried out for, typically, 200 ns, with sampling of coordinate sets carried out at 400 fs increments. Some runs were extended to longer time scales (300 ns).

PMF calculation
A PMF curve was calculated for a dimer of Dye A molecules using the pull code implemented within GROMACS. Distances were constrained for the centers of mass of the aromatic units, with spacing between neighboring points varying between 0.01 nm and 0.1 nm, with smaller spacing used near the global minimum and larger spacing for long distances. Each point was equilibrated for 1 ns, followed by a 50 ns production run.

Analysis work
The interlayer distance, cos The cross sectional area (CSA) for a column was defined as the cross section normal to the column axis. CSAs were calculated by first running a 0.14 nm spherical probe (corresponding to the radius of a water molecule) over the hydrophobic core region only of the aggregate structure (excluding the flexible side arms) to obtain a solvent accessible surface area. The aggregate was then split into two and the procedure was repeated for each segment. The CSA was defined as half the difference of the additional surface area created in this process. This allows a direct comparison to typical experimental CSAs, where the experimental assumption is usually that the hydrophilic arms do not contribute to the core of a chromonic column. The calculations used the molecular graphics program PYMOL. [57] Supporting Information Supporting Information is available from the Wiley Online Library or from the author. Following the experimental assumption, CSA calculations for simulation results have been performed using the hydrophobic core of the dye only, excluding the side arms.       Interaction of chromonic stacks at higher concentration leads to the formation of a novel layered structure, which is suggested as a candidate for a smectic chromonic phase.

Conversion between parallel and anti-parallel associated molecules.
The molecular strain in an aggregate caused by parallel stacking can be relieved by one of two possible mechanisms (Fig S2). One mechanism occurs where one of the dyes in a parallel configuration travels towards the end of the column where it can be eliminated. The other requires a large rotation of the sections of the aggregate either side of this molecule. It is this latter mechanism, which has been observed for both dye systems ( Figure S3). Parallel dye pairs appending a column are far easier to re-orient, with the dye molecule having a far greater degree of rotational freedom than a trapped molecule in the center. Substantial continuous parallel stacking (~4-5 adjacent dye molecules) has also been observed in some simulations and leads to a localized bend in the column ( Figure S3). A strained-bent aggregate ultimately fragments and produces two smaller aggregates. Diffusion of the aggregates within solution eventually leads to recombination, in an attempt to form a fully anti-parallel stacking arrangement between neighboring dye molecules in a column.      Figure S7. Fragmentation of a 30-dye molecule system for Dye B at 30 wt%, and initial steps in the formation of a layered structure. Three snapshots have been selected at different time intervals: a) 40 ns, b) 80 ns, c) 120 ns. Counterions and water have been removed for visual clarity. Figure S8. A seeded double-column stack structure of Dye A (as in Figure S1a,b) is transformed into a single column stack (shown right) over a 200 ns time period. Figure S9. Molecular structure for the core of Dye A and Dye B, highlighting the atoms used to measure the plane of the molecule highlighted in red and the atoms used to describe the direction vector in blue.