Intermolecular interactions in eumelanins: a computational bottom-up approach. I. small building blocks

Abstract

The non-covalent interactions between pairs of the smallest eumelanins building blocks, 5,6-dihydroxy-indole (DHI) and its redox derivatives, are subjected to a systematic theoretical investigation, elucidating their nature and commenting on some of their possible effects on the layered structure of eumelanin. An accurate yet feasible protocol, based on second order perturbation theory, was set up and validated herein, and thereafter used to sample the intermolecular potential energy surfaces of several DHI related dimers. From the analysis of the resulting local minima, the crucial role of stacking interactions is assessed, evidencing strong effects on the geometrical arrangement of the dimer. Furthermore, the absorption spectra of the considered dimers in their most stable arrangements are computed and discussed in relation to the well known eumelanin broadband features. The present findings may help in elucidating several eumelanin features, supporting the recently proposed geometrical order/disorder model (Chen et al., Nat. Commun. 2014, 5, 3859).

---

1 Introduction

As a consequence of their unique features, eumelanins have attracted increasing attention, which has stimulated continuous research in many fields, as testified in several recent reviews dedicated to these materials.^1–7 Besides their peculiar physicochemical polymer-like properties,² the two more well known characteristics of both natural and synthetic eumelanins are the intrinsic free radical character and the broadband UV/vis absorption, fundamental in its photo-protective function.⁸ Thanks to this fascinating ensemble of different properties, eumelanins have found an astonishing wealth of applications as bio-compatible multifunctional platforms in several fields, such as bio-interfaces,⁹ polymer stabilization,¹⁰ organic electronics^4,11 or, more recently, hybrid materials.¹² A further boost has come from the multiple applications of a eumelanin like material, namely polydopamine, whose features have been recently reviewed as being strongly related to those of eumelanins.⁷ In fact, since its first application in 2007 as an efficient coating material,¹³ the properties of polydopamine have been exploited in many technological fields, ranging from nanomedicine to bio-sensors and energy production.^14–16

Notwithstanding the many theoretical^17–38 and experimental^{2,5,13,33,39–50} studies, the structure–property relationships underpinning the above mentioned multipurpose applications are still far from clearly understood, and the patterns of supramolecular assembly or even the chemical nature of eumelanin constituting molecules are still under debate.^7,35,36 The reason for this has been clearly pointed out in a very recent review by D’Ischia et al.,⁷ who traced back the difficulties in establishing a detailed picture of eumelanin structure to the fact that “unlike the vast majority of natural pigments, eumelanins cannot be described in terms of a single well-defined structure, and it is not possible to provide an accurate picture beyond a statistical description of main units and functional groups”. More specifically, since the original amorphous semiconductor model³⁹ was definitively abandoned in light of experimental evidence of a certain degree of supramolecular structure,^{40,41,43,49,51,52} many of eumelanin’s peculiar characteristics (and in particular the broadband absorption spectra) were explained through the so-called chemical disorder model,^3,4,27 where eumelanins are constituted by a plethora of distinct species, which may differ in structure, redox state or geometrical arrangement. In this scenario, most of the resulting eumelanin properties can be rationalized as a consequence of such disorder: the broadband absorption, for instance, might be interpreted as a superposition of the spectra of the multiple constituents of eumelanin. The main drawback of the chemical disorder model is undoubtedly the neglect of both the non-covalent interactions occurring among the various eumelanin constituents and the effect that the supramolecular motifs that originate from these interactions might have on the optical and mechanical properties of eumelanin, as reported in many investigations.^{18,23,36–38,48} Very recently, Chen et al.^36,38 have complemented the chemical disorder model with the geometrical order/disorder picture, where the existence of stacked layers containing a number of eumelanin forming units (protomolecules) and the conformational and configurational disorder within and among them is also taken into account (see also Scheme 3 in ref. 7). Despite the elusive nature of this description, several points seem now to be shared by the scientific community:^1–7 (i) the basic² unit components (see also Fig. 1) of natural and synthetic eumelanins are 5,6-dihydroxy-indole-2-carboxylic acid (DHICA), 5,6-dihydroxy-indole (DHI), its redox form indole-quinone (IQ) and its tautomer quinone-imine (MQ); (ii) the above units are linked through covalent C–C bonds into small oligomers (protomolecules), ranging from tetramers to octamers, although a monomeric-based picture has also been recently proposed;³⁴ (iii) to date, the most probable and investigated oligomeric structures are those of a cyclic, porphyrin-like planar tetramer and linear or cross linked oligomers. It might be interesting to note that the former structures were mainly deduced from computational results,²⁶ whereas the latter were hypothesized, based on experimental measurements;^46,47 (iv) whatever their structure, the protomolecules are arranged into stacked layers, with a typical separation of 3–4 Å; (v) the dimensions of these stacked aggregates are nanometric and no long-range orientational or positional order was found among the layers. Additionally, it is most probable that, within each layer, protomolecules can adopt a variety of geometrical conformations and/or redox states.

Fig. 1 Monomeric species involved in eumelanins^2,7 that are considered in the present work: 5,6-dihydroxy-indole (DHI), its redox forms semi-quinone (SQ) and indole-quinone (IQ) and its tautomer quinone-imine (MQ).

Based on these few points, it appears that it is now of fundamental importance to deeply understand and characterize the non-covalent interactions among protomolecules that are responsible for both the formation of the stacked layers and the tertiary structure arising from their spatial disposition. In this context, the potential of a computational approach can be exploited to unravel the inherent complexity of the supramolecular structure of eumelanin. Surprisingly, to our knowledge, only a few attempts to thoroughly investigate non-covalent interactions in eumelanins have been reported to date in the literature. Indeed, apart from two recent works, where quantum mechanical (QM) methods, such as density functional theory (DFT)³¹ and its tight binding extension combined with the self-consistent charge technique (SCC-DFTB),³⁵ were employed, most of the computational studies dealing with interactions between protomolecules were based on empirical, transferable force-fields (FFs),^23,36,38 whose accuracy in handling non-covalent forces has been recently questioned.^53–55 Moreover, all previous studies dealt only with a few spatial arrangements of the protomolecules. On the contrary, considering the crucial role of the non-covalent interactions that may occur between eumelanin building blocks, the effect of their relative orientation as well as the possible balance between van der Waals (stacked) and hydrogen bond (HB) forces should be also taken into account. Therefore, a more systematic though accurate investigation appears to be necessary.

The first step towards a deeper comprehension of the structural motifs of protomolecules consists of devising a computational method able to accurately investigate the interactions between pairs of the smallest eumelanin building blocks, i.e. DHI and its redox (IQ) and tautomeric (MQ) forms. This is the main objective of the present work. Considering the aromatic character of all species considered, and the presence of one or more –OH functional groups, it is not unlikely that the most favorable arrangements of such pairs will be driven by a delicate balance between stacked (π–π) and HB interactions. The most accurate and reliable QM method to handle such non-covalent forces is probably the coupled cluster technique, evaluated at complete basis set (CCSD(T)@cbs), often referred to as the gold standard of quantum chemistry.⁵⁶ Unfortunately, the prohibitive computational cost of this technique rules out its adoption in sampling a large number of non-covalent dimer conformations. Computational convenience would suggest to turn to DFT dispersion corrected calculations, but the accuracy of the results is known to depend on the choice of the functional and its empirical correction.^57,58 In a rather time-consuming procedure, some purposely computed CCSD(T)@cbs data could be used as a reference to select the most accurate empirical dispersion correction,^57,59–61 although it has been recently reported that the performances of a given functional may vary dramatically with the target dimer and even with its spatial arrangement.^62,63

In this scenario, an alterative approach was adopted in our group, exploiting the combination of the second order Möller–Plesset perturbation theory (MP2) with purposely modified basis sets (therefore named⁶⁴ MP2mod). Very recently, the MP2mod method was successfully applied by us^64,65 to a non-covalent dimer, similar to those encountered in eumelanins, namely quinhydrone. As a matter of fact, MP2mod results obtained for quinhydrone agreed well with reference high-level CCSD(T)@cbs data, whereas, thanks to the reduced basis set dimensions, the observed computational times were remarkably reduced with respect to the other investigated methods.⁶⁴ Given the similarities (also recently pointed out in ref. 34) between quinhydrone’s component monomers and the eumelanin building blocks considered here, MP2mod is expected to give good performances in this case also. Nonetheless, for a more robust assessment of the accuracy of the method, it is here again validated against CCSD(T)@cbs data, purposely computed on a few selected arrangements. Subsequently, MP2mod is employed to calculate the interaction energy of different non-covalent pairs arising from the combination of the considered eumelanin minimal building blocks. Finally, the effect of the resulting dimer conformations on the shape of the absorption spectrum is also considered and discussed in some detail.

2 Computational details

The initial geometry of each monomer was obtained by optimizing its energy at DFT level, using the B3LYP functional with the Dunning’s correlation consistent cc-pVTZ basis set.

The interaction potential energy surface (IPES) of each investigated pair was first sampled by computing the interaction energy in several dimer arrangements, constructed without altering the internal (B3LYP/cc-pVTZ) geometry of each monomer. In this case, the interaction (ΔE^inter) and binding (ΔE^bind) energies between monomers A and B have the same expression:

ΔE^inter = E_AB − (E⁰_A + E⁰_B) = ΔE^bind (1)

where E⁰_A, E⁰_B are the energies of the two isolated monomers in their equilibrium geometry, whereas E_AB is the energy of the whole dimer AB in the chosen conformation.

Next, several local minima of each considered IPES were further explored by optimizing the dimer geometries corresponding to the most representative arrangements found in the previous sampling. In this case, since the monomer geometries may change during the full optimization, ΔE^inter and ΔE^bind may differ, i.e.

ΔE^inter = E_AB − (E_A + E_B) (2)

where E_A and E_B are the energies of the two monomers, computed for each monomer in the geometries assumed within the AB dimer. In this case, the binding energy is

ΔE^bind = E_AB − (E⁰_A + E⁰_B) = ΔE^inter + (ΔE_A + ΔE_B) (3)

where

ΔE_M = E_M − E⁰_M; M = A, B (4)

is the relaxation/distortion energy of monomer A (or B), due to the AB complex formation. To evaluate this term, for each investigated complex, E⁰_M was obtained after a geometry optimization performed on the isolated monomer M at the same level of theory (MP2mod, vide infra) employed in the calculation of ΔE^inter for dimer AB.

For each investigated AB dimer geometry, the computation of the interaction energies was performed at two different levels of theory.

(i) For selected arrangements (either obtained by assembling the two rigid monomers or through dimer optimization) ΔE^inter was computed according to eqn (2), at the CCSD(T)@cbs⁵⁶ level, following the procedure employed by us for quinhydrone,⁶⁴ which is reported in detail in the ESI.†

(ii) All considered dimer pairs were investigated through the MP2mod method,⁶⁴ computing both ΔE^inter and ΔE^bind according to eqn (2) and (3), respectively. All calculations were performed with 6-31G*(α_C, α_O, α_N), where α_C, α_O and α_N are the modified polarization exponents for carbon, oxygen, and nitrogen atoms, respectively. The first two exponents, namely 0.25 and 0.44, were transferred from the ones optimized for quinhydrone,⁶⁴ whereas α_N was optimized in this work through the EXOPT^64,66,67 procedure against the reference CCSD(T)@cbs data obtained for pyrrole, as detailed in the ESI† and in the next section.

Finally, the absorption vertical transition energies for the first 25 excited states were obtained for eumelanin related monomers and dimers at the TD-DFT level, using the B3LYP functional with the aug-cc-pVDZ basis set. Absorption spectra were subsequently obtained by broadening each transition with a Gaussian function with an half-width at half maximum of 0.33 eV. All calculations were performed with the Gaussian09 package.⁶⁸

3 Results and discussion

3.1 MP2mod set up

The eumelanin related monomers considered in this work, and displayed in Fig. 1, all consist of two condensed aromatic moieties, namely an o-hydroquinone derivative and one pyrrole ring. Therefore, a first validation of the reliability of MP2mod can be obtained by testing the method performances on these two smaller molecules. To obtain the α polarization exponents for the modified 6-31G* basis set required by the MP2mod method, two different strategies were followed.

As far as o-hydroquinone derivatives are concerned, the polarization exponents for carbon (0.25) and oxygen (0.44) atoms were transferred from previous work performed on a very similar dimer (quinhydrone).⁶⁴ To test this assumption, o-benzoquinone/o-hydroquinone pairs were prepared in three different geometries (see Fig. A of the ESI†): face-to-face (FF), antiparallel face-to-face (AFF) and hydrogen bonded (HB). The intermolecular distance between the centers of mass of the monomers was varied stepwise in the 2.8–7.0 Å range, and the interaction energy was computed for each of the resulting dimer arrangements at the MP2mod/6-31G*(0.25, 0.44) level. From the computed curves, reported in detail in Fig. B of the ESI,† two main observations may be made. First, as could be expected, it is clear from the difference between the MP2mod energies and their HF contributions that the dispersion forces routed in the dynamic electron correlation dominate the stacked wells (FF and AFF), and play only a minor (but not negligible) role in the H-bonded (HB) dimers. As a consequence, the FF and AFF arrangements are much more favorable than the HB one.

Starting from the most stable AFF geometry, a complete energy optimization was carried out at the MP2mod level. The final arrangement (dimer I) is shown in the left panel of Fig. 2. The interaction energy ΔE^inter of this geometry has been computed according to eqn (2) at both the MP2mod and CCSD(T) levels, as reported in Table 1. For a more extended validation, a second optimization was also carried out at the MP2mod level, starting with the HB minimum conformation. As could be expected from the comparison of the AFF and HB curves (ESI, Fig. B†), the HB dimer optimization ends up in a stacked conformation (dimer II), shown in the right panel of Fig. 2. In contrast with dimer I, in dimer II one hydroquinone hydrogen is found out of the aromatic plane, allowing for the formation of two weak hydrogen bonds, which concur to further stabilize the stacked geometry. MP2mod and CCSD(T) interaction energies for dimer II were then computed, and are compared to each other in Table 1. For both dimers, the agreement between MP2mod and CCSD(T) reference values is rather good, with a ΔE^inter difference less than 0.5 kcal mol⁻¹, and the adoption of the α_C and α_O exponents to handle o-hydroquinone derivatives appears reliable.

Fig. 2 Optimized geometries of the o-benzoquinone/o-hydroquinone pair computed at the MP2mod level, starting from the AFF (dimer I, left) and HB (dimer II, right) arrangements, displayed in Fig. A of the ESI.† Hydrogen bonds in dimer II are highlighted with green dashed lines.

Table 1 Interaction (ΔE^inter) and binding (ΔE^bind) energies computed with different methods for quinone dimer I and dimer II. The ΔE^bind value at the CCSD(T)/cbs level is not reported as the calculation of the relaxation/distortion energy at cbs is not strictly defined

Dimer	Method	Basis set	ΔE^inter (kcal mol^−1)	ΔE^bind (kcal mol^−1)
Dimer I	CCSD(T)	cbs	−5.34	—
Dimer I	MP2mod	6-31G*(0.25, 0.44)	−5.80	−5.35
Dimer II	CCSD(T)	cbs	−9.54	—
Dimer II	MP2mod	6-31G*(0.25, 0.44)	−9.91	−7.04

---

Since no α exponent was available for the nitrogen atom in pyrrole, α_N was optimized following the procedure adopted for quinhydrone in ref. 64, and discussed in some detail in the ESI.† The optimized polarization exponent for the nitrogen atom resulting from the EXOPT minimization is 0.37, a value not far from that found for the oxygen atom in hydroquinone (0.44).⁶⁴ The pyrrole IPES was scanned according to selected cross sections, described in Fig. C of the ESI.† From looking at Fig. D in the ESI,† the MP2mod results match well with the scenario that emerged from previous calculations performed on the pyrrole dimer and reported in the literature (see ref. 69 and references therein). Indeed, the stacked face-to-face (FF) arrangements are found not to be stable, whereas the antiparallel displaced one shows a significant well (∼3.5 kcal mol⁻¹) at rather short ring–ring distances (∼3 Å). More importantly, the most favorable conformations (>5 kcal mol⁻¹) among the four investigated arrangements are those where a hydrogen bond between the –NH group of one monomer and the π cloud of the other is established. To further check the reliability of the MP2mod protocol, two additional dimer optimizations were performed, starting from the two most favorable geometries (T-shaped (TS) and APD, see Fig. C and D in ESI†). The final structures, labeled dimer I and II, are reported in Fig. 3. In the first case, the two monomers are found in a T-shaped like arrangement, where the hydrogen of on the –NH group strongly interacts with the π cloud of the neighboring monomer, the distance with the aromatic plane being around 2 Å. However, because the HF contribution is only slightly attractive (−0.5 kcal mol⁻¹), dispersion forces still play the major role in dimer stabilization. As far as dimer II is concerned, the planar rings are found to be in a stacked conformation, though displaced by ∼2 Å and with the NH bonds pointing in opposite directions. The interaction energies were computed for both dimers at the MP2mod level with the optimized 6-31G*(0.25, 0.37) basis set, finding ΔE^inter = −5.41 and −5.97 for dimer I and II, respectively. Despite these values being well within the interaction energy range (−4.0 to 6.3 kcal mol⁻¹) found with previous high level calculations,⁶⁹ the many local minima expected for the pyrrole IPES do not allow for a direct comparison. To better assess the reliability of MP2mod, the ΔE^inter was computed for dimer I and II also at the CCSD(T)@cbs level, and the results are compared with their MP2mod counterparts in Table 2. It appears that the agreement with the reference values is remarkable, with the differences between the reference and MP2mod values being 0.1 kcal mol⁻¹ and 0.5 kcal mol⁻¹for dimer I and II geometries, respectively.

Fig. 3 Optimized geometries of the pyrrole dimer computed at the MP2mod level, starting from the TS (dimer I, left) and APD (dimer II, right) arrangements, displayed in Fig. C and D of the ESI.†

Table 2 Interaction (ΔE^inter) and binding energies (ΔE^bind) computed with different methods for pyrrole dimer I and dimer II. The ΔE^bind value at the CCSD(T)/cbs level is not reported as the calculation of the relaxation/distortion energy at cbs is not strictly defined

Dimer	Method	Basis set	ΔE^inter (kcal mol^−1)	ΔE^bind (kcal mol^−1)
Dimer I	CCSD(T)	cbs	−5.51	—
Dimer I	MP2mod	6-31G*(0.25, 0.37)	−5.41	−5.35
Dimer II	CCSD(T)	cbs	−6.52	—
Dimer II	MP2mod	6-31G*(0.25, 0.37)	−5.97	−4.77

---

Exploiting the exponent reliability found for both the o-quinone and pyrrole dimer cases, the 6-31G*(0.25, 0.44, 0.37) basis set was implemented, and the MP2mod procedure was applied in the calculations of the interaction energy between pairs of the eumelanin-related monomers displayed in Fig. 1.

3.2 DHI/IQ IPES

The most probable protomolecules that may concur to form the eumelanin structure have been recently reviewed in ref. 7, 36 and 38 . Despite several models being proposed (ranging from monomers to octamers), all candidates share the same kind of building blocks,² i.e. the DHI derivatives shown in Fig. 1. Among all possible combinations, we started the IPES investigation with the DHI/IQ pair. Indeed, beside being one of the most probable combinations that are generated when two of the aforementioned protomolecule models are considered in a stacked arrangement, it is expected that the most stable DHI/IQ non-covalent dimers will originate from dispersion (stacking) interactions and hydrogen bonds, whose delicate interplay is a challenging benchmark for the MP2mod protocol. Finally, many features of the DHI/IQ dimer closely resemble those of quinhydrone, recently⁶⁴ studied by us with the same computational technique.

DHI/IQ pairs were initially placed in eight different arrangements, devised to consider all types of intermolecular interactions (stacking, H–π, hydrogen bond, etc.) that could arise between such non-covalent dimers, i.e. stacked face-to-face (FF), anti-parallel face-to-face (AFF), side-to-face (SF), hydrogen bonded (HB), parallel displaced (PD), anti-parallel displaced (APD), stacked rotated (ST) and not stacked rotated (NST). For each arrangement, a number of dimer geometries were built by moving one monomer with respect to another along one IPES coordinate (either a distance or an angle), without changing the internal geometries of the monomer. All considered configurations and the chosen displacement coordinates are displayed in the insets of Fig. 4. Finally, for each arrangement, the IPES cross sections, also displayed in Fig. 4, were obtained at the MP2mod/6-31G*(0.25, 0.44, 0.37) level by computing the interaction energy of the proper set of configurations.

Fig. 4 MP2mod interaction energy curves, computed for the DHI/IQ pair in different arrangements. From left to right: stacked face-to-face (FF, top) and parallel displaced (PD, bottom), antiparallel face-to-face (AFF, top) and antiparallel displaced (APD, bottom), side-to-face (SF, top) and stacked rotated (ST, bottom), hydrogen bonded (HB, top) and not stacked rotated (NST). The scanned IPES coordinates in abscissa is highlighted in blue in the inset of each panel.

The first general observation that arises by looking at the resulting curves is that the DHI and IQ monomers show rather strong non-covalent interactions, which remarkably depend on their relative orientation, and span a ΔE^inter range between −10 and −3 kcal mol⁻¹. Notwithstanding the dominant stacking interactions, the most stable conformers are found when the molecular planes are twisted (ST, −9.51 kcal mol⁻¹) or displaced along the plane direction (APD, −9.28 kcal mol⁻¹). In all cases, an antiparallel arrangement of the CO dipoles is preferred to a parallel one. It is also worth noticing that in all the stacked configurations considered, the inter-ring distance corresponding to the curve minimum is found between 3.4 to 3.6 Å, a value in good agreement with the range (3–4 Å) proposed^7,36,40,41 for the eumelanin inter-layer distance. Interestingly, while the lack of stacking interactions causes the HB configurations to be the least stable, the side-to-face like pairs (SF, −7.95 kcal mol⁻¹ and NST, −8.10 kcal mol⁻¹) also show a non-negligible interaction, although the distance between the centers of mass of the monomers is around 5 Å. As a consequence, even in the (most likely) hypothesis that eumelanin protomolecules consist of DHI related oligomers, the latter side-to-face arrangements between the building blocks of two neighboring oligomers cannot be ruled out, as they could be realized by a rotation around the C–C bond connecting two protomolecules moieties, if the gain in non-covalent energy is more than the one needed for the dihedral torsion. However, a full investigation of the subtle balance of the interaction energy between oligomer building blocks and the torsional barriers driving the internal structure of the protomolecules is beyond the aims of the present work, but will be the main topic of a future investigation currently in progress in our group.

3.3 MP2mod validation on DHI/IQ dimers

Before attempting to study other DHI related pairs possibly involved in eumelanin, for instance those involving the MQ tautomer, a final test was performed on the optimized basis set to validate the assumption of transferability of polarization exponents. For this purpose, four dimer optimizations were performed at the MP2mod level starting from selected DHI/IQ arrangements, namely ST, APD, NST and HB. The final geometries, labeled A–D, are shown in Fig. 5 and their interaction energies, computed both at the MP2mod and CCSD(T)@cbs levels, are reported in Table 3.

Fig. 5 Optimized geometries of the DHI/IQ pair, starting from selected minima of the APD (dimer A), HB (dimer B), NST (dimer C) and ST (dimer D) curves (see Fig. 4). The geometrical parameters reported in Table 3 are highlighted with blue and green dashed lines.

Table 3 Geometrical parameters, interaction (ΔE^inter, computed at MP2mod/6-31G*(0.25, 0.37, 0.44) and CCSD(T)@cbs levels) and binding (ΔE^bind, last column, only at MP2mod level) energies computed for DHI/IQ dimers A–D displayed in Fig. 5; a and b (or b′) are highlighted in the same Figure, and correspond to the distance between the two aromatic planes and to the minimum distance between a DHI hydrogen and one IQ heavy atom, respectively. In the last row, CCSD(T) and MP2mod computational times on a Xeon E5-2670 2.60 GHz processor are also reported for comparison

Dimer	a (Å)	b(b′) (Å)	ΔE^interCCSD(T) (kcal mol^−1)	ΔE^interMP2mod (kcal mol^−1)	ΔE^bindMP2mod (kcal mol^−1)
A	2.82	2.11	−11.4	−12.4	−10.0
B	2.80	1.92	−14.4	−15.6	−10.7
C	2.95	2.25 (2.54)	−11.0	−12.4	−10.1
D	2.95	1.92	−14.3	−14.8	−12.1
CPU time (hours)			∼15 000	∼1

---

The comparison between the MP2mod values and those obtained using the CCSD(T)@cbs method is rather good, with a standard deviation of ∼1 kcal mol⁻¹, i.e. less than 10% of the total ΔE^inter value. It is worth noticing that the accuracy obtained with the MP2mod protocol is even more appreciable if one considers the computational times, reported in the last row of Table 3, that underline the feasibility of MP2mod calculations, even on a large number of dimer arrangements.

Once the reliability of MP2mod has been finally assessed, a closer look at the four minimum energy conformations found using MP2mod optimizations can help in understanding the type and strength of the non-covalent interactions established between DHI and IQ units. First, a general feature shared by the four pairs is the stacked arrangement of the two monomers, which is a clear indication of the fundamental role of dispersion forces. Indeed, the two molecular planes are in all cases very close, the inter-plane distance being slightly less than 3 Å. Another specific characteristic of the DHI/IQ pair is the formation of a hydrogen bond between the two units, in addition to the intramolecular one found in the DHI monomer for all A–D conformations. The intermolecular H-bond is formed through an out-of-plane rotation of one hydrogen of the –OH group around the C–OH bond. For dimers B and D, this rotation is very effective, as the intermolecular H-bond distance is found to be below 2 Å. It is worth noticing that the same kind of torsion also takes place in the quinhydrone dimer,⁶⁴ and is involved in both dimer stability and electron and proton transfers between the redox pair in solution.⁶⁵ Finally, by comparing the ΔE^inter_MP2mod and ΔE^bind_MP2mod values reported in Table 3 for all investigated dimers, it appears that the distortion energy [ΔE_DHI + ΔE_IQ] reported in eqn (4) ranges from 2 to 5 kcal mol⁻¹, depending on dimer conformation. However, the former contribution (ΔE_DHI) is more than 75% of the total distortion energy, indicating that larger distortions from the monomer equilibrium geometries are located in the DHI species, and essentially consist of the aforementioned dihedral rotation.

3.4 Additional pairs of eumelanin building blocks

In the hypothesis that the structure of eumelanins can be described by groups of stacked layers, in turn composed of small-sized oligomers of covalently bound DHI related units,⁷ it can be expected that the intermolecular interactions between two neighboring building blocks belonging to different layers are very similar to those settled in the above discussed DHI/IQ arrangements. In this scenario, it is clear how the relatively high interaction energies found between DHI and IQ moieties can strongly influence both the relative orientation of two facing oligomers (through dispersion) and the distance between the stacked layers (through H-bonds) to which they belong. Nonetheless, these considerations rely on the results obtained for only one of the possible redox/tautomeric pairs that can be thought to take place in eumelanins; to draw more sound conclusions, other DHI related non-covalent dimers should be investigated. Among the many possible combinations between the monomers shown in Fig. 1, to restrict computational burden here the only pairs that will be only considered are those that can arise when two of the most reliable and investigated⁷ oligomer models proposed for eumelanin (i.e. the planar cyclic tetramer first proposed by Kaxiras et al.²⁶ from computational findings, and the cross-linked, non-cyclic tetramers hypothesized on experimental evidences by the D’Ischia group^46,47) are disposed in stacked arrangements. As far as the non-cyclic models are concerned, the most recurring monomeric building block is DHI, although IQ, MQ and SQ units can be also contained in the cross-linked tetramer.^7,46,47 Hence, besides the DHI/IQ pair, considering the higher occurrence of DHI units within the oligomer, only the DHI/DHI, DHI/MQ and DHI/SQ combinations have been taken into account for additional MP2mod calculations. Moreover, since the cyclic tetrameric model is made up of IQ and MQ building blocks,^26,36 the interactions between these two moieties (IQ/MQ) have been also investigated.

The starting conformations employed in the previous MP2mod optimizations, carried out on the DHI/IQ pair (that eventually led to the A–D dimers displayed in Fig. 5), were employed as starting points also for the other considered dimers, after proper addition/removal of hydrogen atoms. The final dimer arrangements (labeled E to T) are displayed in Fig. 6. By looking at these optimized conformers, and at the geometrical and energy data reported in Table 4 for all pairs, some general conclusions may be drawn.

Fig. 6 From top to bottom, MP2mod-optimized geometries of the DHI/DHI (E–H), DHI/MQ (I–L), IQ/MQ (M–P) and DHI/SQ (Q–T) pairs, starting from the same initial conformations (i.e. antiparallel displaced (APD), hydrogen bonded (HB), not stacked (NST) and stacked (ST)) used in the optimization of DHI/IQ shown in Fig. 5, as indicated in the top panel. The geometrical parameters reported in Table 4 are highlighted with blue and green dashed lines.

Table 4 Geometrical parameters, interaction (ΔE^inter) and binding (ΔE^bind) energies, computed for the DHI related pairs shown in Fig. 6; a and b (or b′) are highlighted in the same Figure, and respectively correspond to the distance between the two aromatic planes and to the minimum distance between a hydrogen of one monomer and one heavy atom of the other

Dimer	a (Å)	b(b′) (Å)	ΔE^inter (kcal mol^−1)	ΔE^bind (kcal mol^−1)
DHI/DHI
E	3.00	2.21	−11.4	−7.8
F	3.10	2.04	−6.9	−4.5
G	2.95	1.95	−15.7	−10.7
H	2.92	2.16 (2.17)	−15.8	−11.3
DHI/MQ
I	2.72	2.00	−21.9	−14.3
J	2.90	1.85 (2.90)	−10.5	−8.0
K	2.70	1.66	−11.4	−6.4
L	2.70	2.02	−18.5	−13.8
IQ/MQ
M	3.20	—	−14.0	−6.3
N	2.80	2.51	−11.8	−7.1
O	2.85	—	−6.7	−4.8
P	—	2.56	−5.9	−5.7
DHI/SQ
Q	3.14	2.56	−6.8	−4.9
R	—	1.91	−3.3	−3.1
S	2.87	2.92	−11.4	−7.2
T	—	2.30	−3.9	−3.5

---

The first observation is that, as already noticed for the DHI/IQ pair, most of the considered optimizations end up in a stacked arrangement, independently of the starting conformation. More precisely, after optimization, only P, R and T dimers are found to not be in a stacked-like arrangement, but it should be mentioned that in all three cases the starting conformation was also not stacked. However, by looking at the ΔE^bind values reported in Table 4, it appears that they are among the less stable conformers. On the contrary, dimers B, F, J and N all start from a side-by-side, H-bonded conformer (see Fig. 4, top left panel), but, during optimization, the molecular planes rotate around the intermolecular H-bond(s) and the two moieties are eventually found in an almost aligned, stacked geometry. By looking at the IQ/MQ pair, where only one or no H-bonds can be formed, it is clear that the stabilization of the stacked arrangements is ensured essentially by dispersion forces, which may contribute up to ∼8 kcal mol⁻¹ to the pair stability. The formation of one or more H-bonds, of course, favors the stability of the dimer with an additional gain in binding energy that can be quantified, depending on the case considered, to be between 2 and 4 kcal mol⁻¹. Furthermore, when H-bonds are formed between the two moieties, the two monomers tend to orient their long axes in an almost parallel (or antiparallel) disposition, whereas, in their absence, the molecular planes are found to be mostly twisted with respect to each other. However, by looking at the binding energies reported in Tables 3 and 4, it seems that there is no systematic trend in the strength of the non-covalent interaction exhibited by different types of monomers pairs, although DHI/IQ, DHI/DHI and DHI/MQ dimers are on average more stable than those involving IQ/MQ and, in particular, DHI/SQ.

3.5 The effect of stacking on the absorption spectrum

Another feature shared by all the investigated pairs concerns the inter-monomer separation between the stacked planes, which falls, in almost all cases, between 2.7 and 3.1 Å. At such small separations, the π clouds of the two monomers show a non-negligible overlap, hence affecting the frontier orbitals in their shapes and energies. For instance, in the quinhydrone dimer, recent TD-DFT calculations⁶⁴ have shown that the stacking distance strongly influences the absorption spectrum, causing the appearance of a broad band in the visible region as the separation between the two monomers forming quinhydrone is reduced. As a matter of fact, concerning the broadband absorption of eumelanin, the important consequences of the excitonic coupling that takes place in stacked arrangements was first pointed out by Stark and coworkers,²³ and has been very recently invoked to complement the chemical disorder model with the latest geometrical order/disorder one.³⁸ However, a possible limitation of these previous computational works may stand in the method employed to originate the stacked conformers (empirical force-fields) as well as in the low level of theory employed for computing the spectrum (ZINDO or semi-empirical Hamiltonian). In this work, thanks to the reduced dimensions of the investigated targets, a more reliable approach was employed both to select the involved geometries (MP2mod) and to compute the optical properties (TD-DFT).

To investigate the effects of the afore-discussed dimers on the absorption properties of eumelanin, TD-DFT calculations were performed on both the isolated monomers and their optimized combinations. The accuracy of the chosen functional (B3LYP) was validated (see Fig. E in the ESI†) by comparing the computed absorption spectrum of the DHI monomer with the available experimental data.⁴² The absorption spectra of all monomers considered in this work have been computed and reported in the top panel of Fig. 7. By looking at the resulting spectral lines, it appears that the differences in the redox form (DHI vs. SQ, vs. IQ) and in the geometrical arrangement of their constituent atoms (MQ vs. IQ) may cause significant shifts in the position of the strong absorption peak in the UV region and/or changes in intensity in the broad band in the visible range. In agreement with the chemical disorder model, when the separate monomer signals are averaged (see light green line in the bottom panel of Fig. 7), a broad band is formed in the UV region, whose intensity gradually increases with decreasing wavelength, as expected for eumelanin.^2,3 On the contrary, not all of the visible region is covered by the average absorption band, and small or null absorption intensities are found in the 350–550 nm range, confirming that the eumelanin broadband absorption, which covers monotonically (see top panel of Fig. 10) the whole UV/visible spectrum, cannot simply originate from the superposition of the absorption lines of the non-interacting monomers.

Fig. 7 Top panel: TD-DFT absorption spectra of the monomer species displayed in Fig. 1: DHI (green line), IQ (blue), SQ (orange) and MQ (magenta). Bottom panel: average of all monomer spectra (light green line) and average (red line) of the absorption spectra computed for DHI/IQ dimers A–D, displayed in Fig. 5.

In Fig. 8 and 9 the absorption spectra obtained for the considered non-covalent dimer arrangements A to T are displayed. From looking at Fig. 8, where only the spectra relative to the DHI/IQ pair are displayed, it is evident that the formation of the stacked complex, and the consequent excitonic coupling, causes the appearance of a much more diffuse absorption band, which in some cases (dimer C) almost covers the whole visible region or presents maxima (dimer A) at wavelengths around 400–450 nm, where the separate monomers do not absorb at all. This last observation is also evidenced in Fig. 7, where the absorption band of the DHI/IQ dimer, averaged over the four considered minima (A–D), is directly compared to the one arising from the monomeric species.

Fig. 8 TD-DFT absorption spectra of the DHI/IQ dimers, in the optimized conformations displayed in Fig. 5: A (black line), B (cyan), C (magenta) and D (blue). The spectra resulting from the superimposition of the former band shapes is reported with a solid red line.

Fig. 9 TD-DFT absorption spectra of the DHI related pairs in the optimized conformations displayed in Fig. 6.

A more complete picture can be achieved by also taking into account the absorption due to the other dimer configurations E–T (see Fig. 9), built by assembling different monomer pairs. Similarly to what was previously found for the DHI/IQ dimer, all the DHI/MQ and IQ/MQ complexes (central panels of Fig. 9) show rather diffuse bands, whose position and intensity strongly depend upon the dimer arrangement. In particular, dimer M, formed from the IQ and MQ monomers, shows an intense band centered around 400 nm, a region not covered when considering all the isolated monomers. Conversely, dimers E–H and Q–T are characterized by a much reduced dependency upon the relative orientation of their constituent monomers, all bands being centered in the high energy region (∼100 nm).

Finally, it is interesting to superimpose all the computed signals, either originating from monomers or non-covalent dimers, and compare the resulting curve with the experimental³ absorption band shape of real eumelanin. The outcomes of both these operations are displayed in Fig. 10. As shown in the upper panel, the computed curve matches well with most of the characteristics of the experimental signal, i.e. the monotonic increase towards lower wavelengths and the broad band covering a great part of the investigated spectrum. Nonetheless, it is important to stress that the computed average displayed in Fig. 10 neglects the contributions to the spectral signal due to chemical heterogeneity, which comes into play when covalent combinations of different building blocks (e.g. tetramers) and their relevant redox states are also taken into account. This however goes beyond the aims of the present study, and is the object of investigations currently in progress in our group.

Fig. 10 Top: experimental³ absorption spectrum of eumelanin. Bottom: absorption spectrum computed in this work by averaging the TD spectra obtained for monomers and their combinations in the discussed optimized conformations.

4 Conclusions

A key point in the investigation of the properties of eumelanin is certainly a thorough comprehension of the intermolecular interactions that determine eumelanin’s peculiar structural properties. Because of the large dimensions of the protomolecules that are thought to compose eumelanin layers, a mandatory step in investigating their interactions is to set up an affordable yet reliable computational protocol. In this context, the Mp2mod method was adopted here, and the heavy atom polarization exponents of the small 6-31G* basis set were purposely tuned to reproduce the intermolecular interactions occurring between the smallest building blocks of eumelanin.

Exploiting the computational feasibility of the tuned basis set, the IPES of several dimer pairs was investigated. From the information retrieved, some general conclusions on the type and extent of the most likely non-covalent interactions characterizing the structure of eumelanin may be drawn. First of all, the stacking interactions are dominant, even if some hydrogen bonds can occur in the considered pairs. The latter indeed reinforce dimer stability but are not able to alter the preferred stacked disposition. More specifically, most of the investigated building block pairs are found with their molecular planes in parallel arrangements, though often twisted or slightly displaced. The computed interaction energies are found to be remarkably dependent on the relative orientation of the monomer, spanning a rather wide range of values from 5 to 12 kcal mol⁻¹. It might be thus speculated that such strong stacked interactions, sometimes complemented by non-negligible hydrogen bonds, should be able to promote the layered structure hypothesized for eumelanin by many authors. As a matter of fact, the extended IPES sampling carried out herein revealed that for the most stable local minima the distance between two stacked rings of the considered building blocks is found to be between 3.0 and 3.5 Å, in agreement with the 3–4 Å range proposed^7,36,40,41 for eumelanin inter-layer distance.

Once the presence of closely stacked aromatic planes was confirmed by calculations, its consequences in the optical properties of the resulting non-covalent dimers were investigated. As a matter of fact, considering the rather small values found for the inter-ring separation, it can be hypothesized that the π clouds that characterize each of the two aromatic moieties will realistically overlap, hence affecting the frontier orbitals involved in the electronic transition. Such supposition was confirmed by observing, for most of the considered pairs, the appearance of a broad band upon dimer formation, whose position and intensity remarkably depend on dimer type and spatial arrangement. On the contrary, when the intermolecular interactions are completely neglected and the absorption spectrum is built as a superimposition of all the signals arising from the isolated monomers, a negligible absorption intensity is found in the 350–550 nm range, in contrast with the well known broad band features of real eumelanin. On the other hand, when the contribution of the stacked non-covalent dimers conformations is included in the calculation of the spectrum, the averaged computed signal shows two specific features that characterize the experimental absorption of eumelanin: a broad band that covers practically the whole investigated spectrum and a monotonic increase towards higher energies.

Despite the fact that these results seem to suggest that many of the considered configurations do effectively take place in eumelanin structure, the covalent bonds that may connect the investigated building blocks to form small protomolecules (as tetrameric structures) should be taken into account. Indeed, work in this direction is currently in progress in our group. In this context, the validation of a reliable and effective method to compute the non-covalent interactions in eumelanin, taking place either between the simple building blocks investigated herein or between larger DHI based oligomers, paves the way to a further step in the comprehension of eumelanin supramolecular motifs.

It should be thus pointed out that this is the first part of a multi-step work, which will hopefully lead us toward an accurate large scale simulation of a model eumelanin material. Indeed, in the present work the interaction between the minimal building blocks was investigated and the proposed computational method reliably validated: a second step will consist of extending the investigation to the interaction between larger oligomers (protomolecules). This could be done, for instance, by resorting to the idea that the interaction energy between two large protomolecules (for instance two DHI tetramers) can be expressed as a sum of the interactions between all pairs of their building blocks.⁷⁰ In this hypothesis, the interaction energy between two protomolecules would strongly depend on the relative distances and orientations of their building blocks. However, unlike the case investigated here, the intramolecular flexibility of the involved oligomers should be also accounted for, since the constraints imposed by the structure of each protomolecule could realistically alter the distribution of geometrical conformations of all possible building blocks pairs. Beside the dimensions of the protomolecules, a further obstacle to such an approach is that the same level of accuracy is required to handle intramolecular flexibility and intermolecular interactions, if a correct balance is to be achieved. Therefore, the next step of this work will consist of applying to different kinds of protomolecules a number computational protocols, developed in our group^70–73 to reliably investigate van der Waals dimers of large dimensions.

Finally, an accurate sampling of the six dimensional interaction potential energy surface (IPES) of the building block pairs considered here could be also exploited to validate the performances of the FFs adopted in previous studies or employed as a reference to parameterize a specifically tailored FF.^66,67,74,75 Once both the achieved intra- and inter-molecular FF descriptions have been validated, the last step of this project will consist of a large-scale simulation on different eumelanin models, based on different combinations of DHI building blocks.

Acknowledgements

The Italian Ministry of Instruction, University and Research (MIUR) is acknowledged for financial funding through the PRIN 2010-11 projects number 2010PFLRJR (Oxidative and Radical Processes, Innovative Aspects and Applications for the Development of Melanin Biopolymers and Antioxidants of Biomedical and Technological Relevance, PROxi) and number 2010FM738P (Photophysical and Photochemical Properties of Organic Compounds in Solution and in Supramolecular Systems).

References

1. M. D’Ischia, A. Napolitano, A. Pezzella, E. J. Land, C. A. Ramsden and P. A. Riley, 5,6-Dihydroxyindoles and Indole-5,6-diones, Adv. Heterocycl. Chem., 2005, 89, 1–63.
2. P. Meredith and T. Sarna, The Physical and Chemical Properties of Eumelanin, Pigm. Cell Res., 2006, 19, 572–594.
3. P. Meredith, B. J. Powell, J. Riesz, S. P. Nighswander-Rempel, M. R. Pederson and E. G. Moore, Towards Structure-Property-Function Relationships for Eumelanin, Soft Matter, 2006, 2, 37–44.
4. M. D’Ischia, A. Napolitano, A. Pezzella, P. Meredith and T. Sarna, Chemical and Structural Diversity in Eumelanins: Unexplored Bio-Optoelectronic Materials, Angew. Chem., Int. Ed., 2009, 48, 3914–3921.
5. M. E. Lynge, R. van der Westen, A. Postma and B. Städler, Polydopamine – a Nature-Inspired Polymer Coating for Biomedical Science, Nanoscale, 2011, 3, 4916–4928.
6. W. E. Bentley and G. F. Payne, Materials science. Nature’s Other Self-Assemblers, Science, 2013, 341, 136–137.
7. M. D’Ischia, A. Napolitano, V. Ball, C.-T. Chen and M. J. Buehler, Polydopamine and Eumelanin: From Structure–Property Relationships to a Unified Tailoring Strategy, Acc. Chem. Res., 2014, 47, 3541–3550.
8. A. Corani, A. Huijser, T. Gustavsson, D. Markovitsi, P.-A. K. Malmqvist, A. Pezzella, M. D’Ischia and V. Sundström, Superior Photoprotective Motifs and Mechanisms in Eumelanins Uncovered, J. Am. Chem. Soc., 2014, 136, 11626–11635.
9. C. J. Bettinger, J. P. Bruggeman, A. Misra, J. T. Borenstein and R. Langer, Biocompatibility of Biodegradable Semiconducting Melanin Films for Nerve Tissue Engineering, Biomaterials, 2009, 30, 3050–3057.
10. K. Shanmuganathan, J. H. Cho, P. Iyer, S. Baranowitz and C. J. Ellison, Thermooxidative Stabilization of Polymers Using Natural and Synthetic Melanins, Macromolecules, 2011, 44, 9499–9507.
11. P. Meredith, C. J. Bettinger, M. Irimia-Vladu, A. B. Mostert and P. E. Schwenn, Electronic and Optoelectronic Materials and Devices Inspired by Nature, Rep. Prog. Phys., 2013, 76, 034501.
12. E. A. Prasetyanto, P. Manini, A. Napolitano, O. Crescenzi, M. D’Ischia and L. De Cola, Towards Eumelanin@zeolite Hybrids: Pore-Size-Controlled 5,6-Dihydroxyindole Polymerization, Chem.–Eur. J., 2014, 20, 1597–1601.
13. H. Lee, S. M. Dellatore, W. M. Miller and P. B. Messersmith, Mussel-Inspired Surface Chemistry for Multifunctional Coatings, Science, 2007, 318, 426–430.
14. D. R. Dreyer, D. J. Miller, B. D. Freeman, D. R. Paul and C. W. Bielawski, Perspectives on Poly(dopamine), Chem. Sci., 2013, 4, 3796–3802.
15. J. Sedó, J. Saiz-Poseu, F. Busqué and D. Ruiz-Molina, Biomimetics: Catechol-Based Biomimetic Functional Materials, Adv. Mater., 2013, 25, 792.
16. Y. Liu, K. Ai and L. Lu, Polydopamine and its Derivative Materials: Synthesis and Promising Applications in Energy, Environmental, and Biomedical Fields, Chem. Rev., 2014, 114, 5057–5115.
17. L. E. Bolívar-Marinez, D. S. Galvão and M. J. Caldas, Geometric and Spectroscopic Study of Some Molecules Related to Eumelanins. 1. Monomers, J. Phys. Chem. B, 1999, 103, 2993–3000.
18. K. Bochenek and E. Gudowska-Nowak, Fundamental Building Blocks of Eumelanins: Electronic Properties of Indolequinone-Dimers, Chem. Phys. Lett., 2003, 373, 532–538.
19. K. Bochenek and E. Gudowska-Nowak, Electronic Properties of Random Polymers: Modelling Optical Spectra of Melanins, Acta Phys. Pol., B, 2003, 34, 2775–2790.
20. Y. V. Il’ichev and J. D. Simon, Building Blocks of Eumelanin: Relative Stability and Excitation Energies of Tautomers of 5,6-Dihydroxyindole and 5,6-Indolequinone, J. Phys. Chem. B, 2003, 107, 7162–7171.
21. K. B. Stark, J. M. Gallas, G. W. Zajac, M. Eisner and J. T. Golab, Spectroscopic Study and Simulation from Recent Structural Models for Eumelanin: I. Monomer, Dimers, J. Phys. Chem. B, 2003, 107, 3061–3067.
22. K. B. Stark, J. M. Gallas, G. W. Zajac, M. Eisner and J. T. Golab, Spectroscopic Study and Simulation from Recent Structural Models for Eumelanin: II. Oligomers, J. Phys. Chem. B, 2003, 107, 11558–11562.
23. K. B. Stark, J. M. Gallas, G. W. Zajac, J. T. Golab, S. Gidanian, T. McIntire and P. J. Farmer, Effect of Stacking and Redox State on Optical Absorption Spectra of Melanins – Comparison of Theoretical and Experimental Results, J. Phys. Chem. B, 2005, 109, 1970–1977.
24. B. J. Powell, T. Baruah, N. Bernstein, K. Brake, R. H. McKenzie, P. Meredith and M. R. Pederson, A First-Principles Density-Functional Calculation of the Electronic and Vibrational Structure of the Key Melanin Monomers, J. Chem. Phys., 2004, 120, 8608–8615.
25. B. Powell, 5,6-Dihydroxyindole-2-carboxylic Acid: a First Principles Density Functional Study, Chem. Phys. Lett., 2005, 402, 111–115.
26. E. Kaxiras, A. Tsolakidis, G. Zonios and S. Meng, Structural Model of Eumelanin, Phys. Rev. Lett., 2006, 97, 218102.
27. M. L. Tran, B. J. Powell and P. Meredith, Chemical and Structural Disorder in Eumelanins: a Possible Explanation for Broadband Absorbance, Biophys. J., 2006, 90, 743–752.
28. A. Pezzella, L. Panzella, O. Crescenzi, A. Napolitano, S. Navaratman, R. Edge, E. J. Land, V. Barone and M. D’Ischia, Short-Lived Quinonoid Species from 5,6-Dihydroxyindole Dimers en Route to Eumelanin Polymers: Integrated Chemical, Pulse Radiolytic, and Quantum Mechanical Investigation, J. Am. Chem. Soc., 2006, 128, 15490–15498.
29. A. L. Sobolewski and W. Domcke, Photophysics of Eumelanin: Ab Initio Studies on the Electronic Spectroscopy and Photochemistry of 5,6-Dihydroxyindole, ChemPhysChem, 2007, 8, 756–762.
30. H. Okuda, K. Wakamatsu, S. Ito and T. Sota, Possible Oxidative Polymerization Mechanism of 5,6-Dihydroxyindole from Ab Initio Calculations, J. Phys. Chem. A, 2008, 112, 11213–11222.
31. S. Meng and E. Kaxiras, Theoretical Models of Eumelanin Protomolecules and Their Optical Properties, Biophys. J., 2008, 94, 2095–2105.
32. S. Meng and E. Kaxiras, Theoretical Models of Eumelanin Protomolecule and Its Photoprotection Mechanism, Biophys. J., 2009, 96, 300a.
33. A. Pezzella, L. Panzella, O. Crescenzi, A. Napolitano, S. Navaratnam, R. Edge, E. J. Land, V. Barone and M. D’Ischia, Lack of Visible Chromophore Development in the Pulse Radiolysis Oxidation of 5,6-Dihydroxyindole-2-Carboxylic Acid Oligomers: DFT Investigation and Implications for Eumelanin Absorption Properties, J. Org. Chem., 2009, 74, 3727–3734.
34. D. R. Dreyer, D. J. Miller, B. D. Freeman, D. R. Paul and C. W. Bielawski, Elucidating the Structure of Poly(dopamine), Langmuir, 2012, 28, 6428–6435.
35. J. Liebscher, R. Mrówczyski, H. A. Scheidt, C. Filip, N. D. Hdade, R. Turcu, A. Bende and S. Beck, Structure of Polydopamine: a Never-Ending Story?, Langmuir, 2013, 29, 10539–10548.
36. C.-T. Chen, V. Ball, J. J. de Almeida Gracio, M. K. Singh, V. Toniazzo, D. Ruch and M. J. Buehler, Self-Assembly of Tetramers of 5,6-Dihydroxyindole Explains the Primary Physical Properties of Eumelanin: Experiment, Simulation, and Design, ACS Nano, 2013, 7, 1524–1532.
37. S. Lin, C.-T. Chen, I. Bdikin, V. Ball, J. Grácio and M. J. Buehler, Tuning Heterogeneous Poly(dopamine) Structures and Mechanics: in Silico Covalent Cross-Linking and Thin Film Nanoindentation, Small, 2014, 10, 457–464.
38. C.-T. Chen, C. Chuang, J. Cao, V. Ball, D. Ruch and M. J. Buehler, Excitonic Effects from Geometric Order and Disorder Explain Broadband Optical Absorption in Eumelanin, Nat. Commun., 2014, 5, 3859.
39. J. McGinness, P. Corry and P. Proctor, Amorphous Semiconductor Switching in Melanins, Science, 1974, 183, 853–855.
40. J. Cheng, S. C. Moss, M. Eisner and P. Zschack, X-Ray Characterization of Melanins – I, Pigm. Cell Res., 1994, 7, 255–262.
41. J. Cheng, S. C. Moss and M. Eisner, X-Ray Characterization of Melanins – II, Pigm. Cell Res., 1994, 7, 263–273.
42. X. Zhang, C. Erb, J. Flammer and W. M. Nau, Absolute Rate Constants for the Quenching of Reactive Excited States by Melanin and Related 5,6-Dihydroxyindole Metabolites: Implications for Their Antioxidant Activity, Photochem. Photobiol., 2000, 71, 524–533.
43. K. C. Littrell, J. M. Gallas, G. W. Zajac and P. Thiyagarajan, Structural Studies of Bleached Melanin by Synchrotron Small-angle X-ray Scattering, Photochem. Photobiol., 2003, 77, 115–120.
44. S. P. Nighswander-Rempel, J. Riesz, J. Gilmore, J. P. Bothma and P. Meredith, Quantitative fluorescence Excitation Spectra of Synthetic Eumelanin, J. Phys. Chem. B, 2005, 109, 20629–20635.
45. S. P. Nighswander-Rempel, J. Riesz, J. Gilmore and P. Meredith, A Quantum Yield Map for Synthetic Eumelanin, J. Chem. Phys., 2005, 123, 194901.
46. L. Panzella, A. Pezzella, A. Napolitano and M. D’Ischia, The First 5,6-dihydroxyindole Tetramer by Oxidation of 5,5′,6,6′-Tetrahydroxy-2,4′-Biindolyl and an Unexpected Issue of Positional Reactivity en Route to Eumelanin-Related Polymers, Org. Lett., 2007, 9, 1411–1414.
47. A. Pezzella, L. Panzella, A. Natangelo, M. Arzillo, A. Napolitano and M. D’Ischia, 5,6-Dihydroxyindole Tetramers with “Anomalous” Interunit Bonding Patterns by Oxidative Coupling of 5,5′,6,6′-Tetrahydroxy-2,7′-Biindolyl: Emerging Complexities on the Way Toward an Improved Model of Eumelanin Buildup, J. Org. Chem., 2007, 72, 9225–9230.
48. A. Pezzella, A. Iadonisi, S. Valerio, L. Panzella, A. Napolitano, M. Adinolfi and M. D’Ischia, Disentangling Eumelanin “Black Chromophore”: Visible Absorption Changes as Signatures of Oxidation State- and Aggregation-Dependent Dynamic Interactions in a Model Water-Soluble 5,6-Dihydroxyindole Polymer, J. Am. Chem. Soc., 2009, 131, 15270–15275.
49. A. A. R. Watt, J. P. Bothma and P. Meredith, The Supramolecular Structure of Melanin, Soft Matter, 2009, 5, 3754–3760.
50. A. B. Mostert, G. R. Hanson, T. Sarna, I. R. Gentle, B. J. Powell and P. Meredith, Hydration-Controlled X-band EPR Spectroscopy: a Tool for Unravelling the Complexities of the Solid-state Free Radical in Eumelanin, J. Phys. Chem. B, 2013, 117, 4965–4972.
51. W. L. Cheun, The Chemical Structure of Melanin, Pigm. Cell Res., 2004, 17, 422–423.
52. J. D. Simon and S. Ito, Reply to Letter: The chemical Structure of Melanin, Pigm. Cell Res., 2004, 17, 423–424.
53. M. Kolár, K. Berka, P. Jurecka and P. Hobza, On the Reliability of the AMBER Force Field and its Empirical Dispersion Contribution for the Description of Noncovalent Complexes, ChemPhysChem, 2010, 11, 2399–2408.
54. M. Zgarbová, M. Otyepka, J. Sponer, P. Hobza and P. Jurecka, Large-scale Compensation of Errors in Pairwise-Additive Empirical Force Fields: Comparison of AMBER Intermolecular Terms with Rigorous DFT-SAPT Calculations, Phys. Chem. Chem. Phys., 2010, 12, 10476–10493.
55. S. Grimme, A General Quantum Mechanically Derived Force Field (QMDFF) for Molecules and Condensed Phase Simulations, J. Chem. Theory Comput., 2014, 4497–4514.
56. J. Rezáč and P. Hobza, Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the Gold Standard, CCSD(T) at the Complete Basis Set Limit?, J. Chem. Theory Comput., 2013, 9, 2151–2155.
57. S. Grimme, Density functional theory with London dispersion corrections, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 211–228.
58. G. Paytakov, T. Dinadayalane and J. Leszczynski, Toward Selection of Efficient Density Functionals for van der Waals Molecular Complexes: Comparative Study of C–H⋯π and N–H⋯π Interactions, J. Phys. Chem. A, 2015, 119(7), 1190–1200.
59. L. A. Burns, A. Vázquez-Mayagoitia, B. G. Sumpter and C. D. Sherrill, Density-functional approaches to noncovalent interactions: a comparison of dispersion corrections (DFT-D), exchange-hole dipole moment (XDM) theory, and specialized functionals, J. Chem. Phys., 2011, 134, 084107.
60. A. Vazquez-Mayagoitia, C. D. Sherrill, E. Apra and B. G. Sumpter, An Assessment of Density Functional Methods for Potential Energy Curves of Nonbonded Interactions: The XYG3 and B97-D Approximations, J. Chem. Theory Comput., 2010, 6, 727–734.
61. S. N. Steinmann, C. Piemontesi, A. Delachat and C. Corminboeuf, Why are the Interaction Energies of Charge-Transfer Complexes Challenging for DFT?, J. Chem. Theory Comput., 2012, 8, 1629–1640.
62. S. S. Leang, F. Zahariev and M. S. Gordon, Benchmarking the performance of time-dependent density functional methods, J. Chem. Phys., 2012, 136, 104101.
63. M. Marianski, A. Oliva and J. J. Dannenberg, A Reinvestigation of the Dimer of para-Benzoquinone and Pyrimidine with MP2, CCSD(T), and DFT Using Functionals Including Those Designed to Describe Dispersion, J. Phys. Chem. A, 2012, 116, 8100–8105.
64. V. Barone, I. Cacelli, O. Crescenzi, M. D’Ischia, A. Ferretti, G. Prampolini and G. Villani, Unraveling the interplay of different contributions to the stability of the quinhydrone dimer, RSC Adv., 2014, 4, 876–885.
65. V. Barone, I. Cacelli, A. Ferretti, G. Prampolini and G. Villani, Proton and Electron Transfer Mechanisms in the Formation of Neutral and Charged Quinhydrone-Like Complexes: A Multilayered Computational Study, J. Chem. Theory Comput., 2014, 10, 4883–4895.
66. I. Cacelli, C. F. Lami and G. Prampolini, Force-field modeling through quantum mechanical calculations: molecular dynamics simulations of a nematogenic molecule in its condensed phases, J. Comput. Chem., 2009, 30, 366–378.
67. I. Cacelli, A. Cimoli, P. R. Livotto and G. Prampolini, An automated approach for the parameterization of accurate intermolecular force-fields: pyridine as a case study, J. Comput. Chem., 2012, 33, 1055–1067.
68. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, P. V. Parandekar, N. J. Mayhall, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fo, Gaussian09, Revision C.01, Gaussian, Inc., Wallingford CT, 2009.
69. M. Kabelác, P. Hobza and V. Spirko, The Structure and Vibrational Dynamics of the Pyrrole Dimer, Phys. Chem. Chem. Phys., 2009, 11, 3885–3891.
70. C. Amovilli, I. Cacelli, S. Campanile and G. Prampolini, Calculation of the intermolecular energy of large molecules by a fragmentation scheme: application to the 4-n-pentyl-4-cyanobiphenyl (5CB) dimer, J. Chem. Phys., 2002, 117, 3003–3012.
71. I. Cacelli and G. Prampolini, Parametrization and Validation of Intramolecular Force Fields Derived from DFT Calculations, J. Chem. Theory Comput., 2007, 3, 1803–1817.
72. I. Cacelli, A. Cimoli and G. Prampolini, Geometry Optimization of Large and Flexible van der Waals Dimers: A FragmentationReconstruction Approach, J. Chem. Theory Comput., 2010, 6, 2536–2546.
73. V. Barone, I. Cacelli, N. De Mitri, D. Licari, S. Monti and G. Prampolini, Joyce and Ulysses: integrated and user-friendly tools for the parameterization of intramolecular force fields from quantum mechanical data, Phys. Chem. Chem. Phys., 2013, 15, 3736–3751.
74. I. Cacelli, G. Cinacchi, G. Prampolini and A. Tani, Computer simulation of solid and liquid benzene with an atomistic interaction potential derived from ab initio calculations, J. Am. Chem. Soc., 2004, 126, 14278–14286.
75. C. Amovilli, I. Cacelli, G. Cinacchi, L. Gaetani, G. Prampolini and A. Tani, Structure and dynamics of mesogens using intermolecular potentials derived from ab initio calculations, Theor. Chem. Acc., 2007, 117, 885–901.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra03773e

---