Next Article in Journal
Development of Azo Dye Immobilized Sulfonated Poly (Glycidyl Methacrylate) Polymer Composite as Novel Adsorbents for Water Treatment Applications: Methylene Blue Immobilization Isotherm, Kinetic, Thermodynamic, and Simulations Studies
Next Article in Special Issue
Atoms-In-Molecules’ Faces of Chemical Hardness by Conceptual Density Functional Theory
Previous Article in Journal
Discovery and Mechanistic Investigation of Piperazinone Phenylalanine Derivatives with Terminal Indole or Benzene Ring as Novel HIV-1 Capsid Modulators
Previous Article in Special Issue
Synthesis, Characterization, Thermal Analysis, DFT, and Cytotoxicity of Palladium Complexes with Nitrogen-Donor Ligands
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Molecules
Volume 27
Issue 23
10.3390/molecules27238417
molecules-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Structural Investigation of DHICA Eumelanin Using Density Functional Theory and Classical Molecular Dynamics Simulations

by
Sepideh Soltani
1,
Shahin Sowlati-Hashjin
2,
Conrard Giresse Tetsassi Feugmo
3 and
Mikko Karttunen
1,4,*
1
Department of Physics and Astronomy, The University of Western Ontario, 1151 Richmond Street, London, ON N6A 3K7, Canada
2
Institute of Biomedical Engineering, University of Toronto, Toronto, ON M5S 3G9, Canada
3
Department of Chemistry, University of Waterloo, 200 University Ave. West, Waterloo, ON N2L 3G1, Canada
4
Department of Chemistry, The University of Western Ontario, 1151 Richmond Street, London, ON N6A 5B7, Canada
*
Author to whom correspondence should be addressed.
Molecules 2022, 27(23), 8417; https://doi.org/10.3390/molecules27238417
Submission received: 28 October 2022 / Revised: 21 November 2022 / Accepted: 22 November 2022 / Published: 1 December 2022
(This article belongs to the Special Issue Density Functional Theory in the Age of Chemical Intelligence)
Browse Figures
Versions Notes

Abstract

:
Eumelanin is an important pigment, for example, in skin, hair, eyes, and the inner ear. It is a highly heterogeneous polymer with 5,6-dihydroxyindole-2-carboxylic acid (DHICA) and 5,6-dihydroxyindole (DHI) building blocks, of which DHICA is reported as the more abundant in natural eumelanin. The DHICA-eumelanin protomolecule consists of three building blocks, indole-2-carboxylic acid-5,6-quinone (ICAQ), DHICA and pyrrole-2,3,5-tricarboxylic acid (PTCA). Here, we focus on the self-assembly of DHICA-eumelanin using multi-microsecond molecular dynamics (MD) simulations at various concentrations in aqueous solutions. The molecule was first parameterized using density functional theory (DFT) calculations. Three types of systems were studied: (1) uncharged DHICA-eumelanin, (2) charged DHICA-eumelanin corresponding to physiological pH, and (3) a binary mixture of both of the above protomolecules. In the case of uncharged DHICA-eumelanin, spontaneous aggregation occurred and water molecules were present inside the aggregates. In the systems corresponding to physiological pH, all the carboxyl groups are negatively charged and the DHICA-eumelanin model has a net charge of − 4 . The effect of K + ions as counterions was investigated. The results show high probability of binding to the deprotonated oxygens of the carboxylate anions in the PTCA moiety. Furthermore, the K + counterions increased the solubility of DHICA-eumelanin in its charged form. A possible explanation is that the charged protomolecules favor binding to the K + ions rather than aggregating and binding to other protomolecules. The binary mixtures show aggregation of uncharged DHICA-eumelanins; unlike the charged systems with no aggregation, a few charged DHICA-eumelanins are present on the surface of the uncharged aggregation, binding to the K + ions.

1. Introduction

Eumelanins are black-brown insoluble pigments [1] that are found, for example, in mammalian skin, hair, and eyes [2,3]. Eumelanins have unique properties, such as UV protection [4,5,6,7], antioxidant activity, and free-radical scavenging [8,9]. They have electronic charge carrying properties [10,11], and their amorphous semiconductivity differs from typical polymer semiconductors [12,13].
Eumelanins are highly heterogeneous polymers with no fully established molecular structure [14,15]. Experiments have, however, confirmed that the building blocks of eumelanins are approximately 15–20 Å long, as reported by X-ray measurements [16,17], scanning tunneling microscopy (STM) [18,19], and ultra-high resolution scanning electron microscopy (UHR)-SEM [20]. Furthermore, STM and transmission electron microscopy (TEM) measurements have shown that synthetic crystallized eumelanins can be graphite-like [21]. Their characteristic features have been identified to originate from the aromatic rings with π − π interactions [22], hydrogen bonding, van der Waals interactions, and the high number of carboxylic acid residues with negative charges [23,24].
Even though the exact chemical structure of eumelanin remains elusive, it is well-established that eumelanins are composed of two main building blocks, 5,6-dihydroxyindole-2-carboxylic acid (DHICA) and 5,6-dihydroxyindole (DHI), at various levels of oxidation and linked through covalent C-C bonds into small oligomers [1,12]. The experimentally observed stacking of DHI [25,26,27,28] and bundling of DHICA [3,29] in eumelanin is due to structural disorder; this structural complexity is further emphasized by the fact that the shapes and proportions of monomers differ depending on the source of eumelanin extraction [20,30]. For instance, Pezzella et al. [31] showed that eumelanin obtained from sepia melanin is composed of approximately 20 % DHI and 75 % DHICA, and Magarella et al. [32] demonstrated that the DHICA:DHI ratio for eumelanin obtained from sepia melanin depends on the purification procedure.
TEM images of eumelanins with DHICA moiety show aggregation with rod-like granular shapes [3,33,34] with dimers larger than 100 nm. Similarly, SEM analysis of bovine melanosomes [26] shows elongated-shapes. DHI eumelanin shows different behaviour; experiments have reported black-colored onion-like aggregates of 50 nm [3,28] and stacking of planar protomolecules. Featureless UV/Vis absorption, which could be due to non-covalent π − π interaction between the rings in heteroaromatic systems, has been reported in simulations of DHI-eumelanins [28,35,36]. It is plausible that the carboxylate group in DHICA may be responsible for twisting the structure and leading to weak internal interactions and aggregation. In experiments [3], DHICA-eumelanins have shown a distinct band around 310 nm with only little contribution in the visible region of UV/Vis absorption. Based on computer simulations, it has been suggested the that broadening of eumelanin spectra may be due to DHICA monomers with π − π interactions [35]. Compared to DHI, DHICA-eumelanins are lighter in color, have weaker aggregation due to their twisted structure, and have stronger redox and photo-protective properties than DHI-eumelanin [3,37]. It has been reported that DHICA shows potent hydroxyl radical scavenging activity [38] and higher antioxidant activity compared to DHI-eumelanin [9,39].
While several density functional theory (DFT) calculations [28,40,41,42,43,44] and molecular dynamics (MD) simulations [28,35,45] have been reported for DHI-eumelanins, only a few DFT calculations [11,46,47] and one MD simulation [35] of DHICA-eumelanins could be found in the literature review. Although it is not immediately clear why DHICA has received significantly less attention, it can be speculated that it is due to its chemical and structural complexity [48]. DHICA-eumelanin’s complex twisted structure together with the fact that the aggregates it forms are more amorphous and much less well-defined compared to the stacks and onion-like aggregates that DHI-eumelanin forms are most likely the main reasons. In addition, although DHI is less complex, it poses many challenges of its own [28,40,45,49]. From the physical perspective, DHICA is a very interesting molecule that even has potential in bioelectronics [11].

2. Materials and Methods

The model that is used for the classical atomistic MD simulations in this paper is named DHICA-eumelanin (Figure 1), as it is based on 5,6-dihydroxyindole-2-carboxylic acid (DHICA) moieties. DHICA-eumelanin has three building block units: (1) DHICA’s tautomer, the indole-2-carboxylic acid-5,6-quinone (ICAQ) moiety; (2) the DHICA moiety; and (3) pyrrole-2,3,5-tricarboxylic acid (PTCA), derived from photo-induced oxidative degradation of the DHICA moiety [50,51]. Figure 1a shows the chemical stucture of the model: ICAQ and DHICA are bonded at atoms C7-C4 ′ , PTCA is bonded to the DHICA moiety at C-C7 ′ , and one of the three carboxylic acids converts to ketone in covalent binding. This model contains the 4-7 ′ DHICA dimer, which is the most stable dimer of DHICA [11]. These units generate a non-planar structure [36]; unlike DHI eumelanin, which has a planar shape [28,40,42,45], they bundle and are not able to fully stack. Hydrogen bonding between the carboxylic groups and π − π interactions between the aromatic rings are the key interactions between the DHICA-eumelanin molecules [14,35].

2.1. Geometry Optimization and Validation of the DHICA-Eumelanin Model

First, a model for the DHICA monomer was generated, and its geometry was optimized and validated using DFT calculations. The DFT results are collected in Tables S1–S4, and the atomic coordinates are provided in Table S1. The optimized structure for the DHICA monomer is shown in Figure S1. After geometry optimization (Tables S1–S3), the dipole moment as well as the molecular vibrational modes and their infrared intensities (Table S4) were evaluated; only the vibrational modes of the main functional groups (those larger than 1000 cm − 1 ) are tabulated. Comparison of the bond lengths and angles with the DFT calculations of Powell et al. and Okuda et al. [48,54] in Tables S2 and S3, respectively, show good agreement. The wavenumbers and infrared intensities from Okuda’s results are reported for comparison (Table S5). The vibrational modes are marginally higher than in Okuda’s calculations, which used the Becke three-parameter Lee-Yang-Parr (B3LYP) functional and 6-31G(d,p) basis set and the homogeneous scale factor of 0.975 [54]; it is worth noticing that the DFT calculations and experiments of Okuda et al. [54] showed differences. Potential explanations for this are: (1) the free DHICA monomer in the DFT calculations was in the gas phase, while the experiments were performed using solid samples; (2) the vibrational frequencies were calculated using the harmonic approximation: and (3) the calculated absorption bands were not pure but were a combination of vibrational modes, and as such the intermolecular hydrogen bonds had an influence on the vibrational modes.
The second verification of the model was the calculation of the dipole moment of an individual DHICA monomer. The result, 3.0 Debye, is in good agreement with the DFT calculation of Matta et al., who reported 2.9 Debye [11].
The next step towards classical MD simulations is optimizing the geometry of the DHICA-eumelanin model (ICAQ-DHICA-PTCA) using DFT calculations. The optimized structure is shown in Figure S2. Geometry optimization was carried out at the B3LYP level of theory in the gas phase [55,56,57,58,59,60] combined with the split valence 6-311G(d,p) basis set [61,62,63]. Quantum mechanical calculations were performed using the Gaussian 09 [64] software. The last verification of the optimized geometry was the C6-C7-C4 ′ -C5 ′ dihedral angle (see Figure 1a for the structure). The result of 59.2 ∘ is in good agreement with Pezzella et al. [36,65], who reported twisted conformations of tautomeric forms of the two-electron oxidation products of DHICA-dimers 4,7 ′ -biindolyls (negatively charged carboxylate). They performed geometry optimization at the PBE0/6-31+G(d,p) level of theory and with the dihedral angle varying in the range of 47–63.1 ∘ or 112–120 ∘ . The lowest energy was at the angle of 47.4 ∘ . On the other hand, our model contains uncharged ICAQ-DHICA-PTCA moieties, and the closest tautomer to all of their structures has the angle of 63.1 ∘ (ICAQ-DHICA and the same C5,OH and C6,OH angles). The variation could be due to our system being uncharged, or because in our system the PTCA moiety is bonded to the ICAQ-DHICA moieties, meaning that PTCA could affect the dihedral angle of DHICA-ICAQ.
The LigParGen server [66] was used to obtain the initial parameters of the bonded and non-bonded interactions; the parameterization is compatible with the Optimized Potentials for Liquid Simulations for All Atoms (OPLS-AA) [67,68] force field. The partial charges were computed using electrostatic potential (ESP) [69,70] fitting over the van der Waals surface grid with the Merz-Singh-Kollman (MK) scheme [70,71] at DFT level.

2.2. Classical MD Simulations

The MD simulations were performed using the Gromacs 2019.3 package [72] with the OPLS-AA force field [67,68]. This force field has been successfully used for similar systems [45]. The positions of the atoms are shown in Table S6. The bond lengths and angles from the DFT calculations and MD after energy minimization are compared in Tables S7 and S8.
The simple point charge (SPC) model [73] was used for water. An integration time step of 0.5  fs was applied. The Lennard–Jones interactions and the real-space part of threlectrostatic interactions were cut off at 1.0  nm. The particle-mesh Ewald (PME) [74,75] method was used for long-range electrostatics, with the reciprocal-space interactions evaluated on a 0.12 nm grid. The DHICA-eumelanin and water molecules were separately coupled to a heat bath using the V-rescale algorithm [76] at 300 K and 0.1  ps coupling constant. Periodic boundary conditions were used in all directions and hydrogens were constrained using the P-LINCS algorithm [77].
After energy minimization, a pre-equilibration step was performed in the NVT (constant particle number, volume and temperature) ensemble for 2 ns. This was followed by a second pre-equilibration step in the constant particle number, pressure, and temperature (NPT) ensemble using a Parrinello–Rahman barostat [78] at 1 bar with compressibility of 4.5 × 10 − 5 bar − 1 and time constant of 2.0  ps for 2.0  ns.
To investigate size dependence, four different system sizes and nine different systems were simulated: four (with 2, 4, 27 and 50 DHICA-eumelanins) using uncharged DHICA-eumelanin, four (with 2, 4, 27 and 50 DHICA-eumelanins) with fully deprotonated carboxylic groups, and one mixed system with 25 uncharged and 25 deprotonated DHICA-eumelanins. Potassium counterions were added to maintain overall charge neutrality when necessary.
The uncharged systems were considered as the reference systems. In the range of acidic to physiological pH, all the carboxylic groups of DHICA-eumelanin are ionized [24,52,53] (Figure 1b); system details are provided in Table 1. The sizes of the simulation boxes were chosen to be large enough to ensure that the DHICA-eumelanins were well surrounded by water molecules and the size dependent artificial hydrophobic effect was absent [79,80]. Visualizations were performed using Visual Molecular Dynamics (VMD) [81] and PyMol [82].

3. Results and Discussion

Eumelanins are known to be insoluble in water and to form aggregates. Because DHICA-eumelanin has three twisted moieties/planes, the aggregation process and the resulting structures appear amorphous. When two eumelanins are close enough to attract each other, two of their aromatic rings can stack, and the rest of moieties are free to interact with water molecules and/or neighboring DHICA-eumelanins, Figure 2.

3.1. Radial Distribution Functions

Figure 3 and Figure 4 show the radial distribution functions (RDF) between the centers of masses (COMs) of selected planes (see Figure 1a for the structures) from the last 100 ns. RDFs for systems of two, four, and fifty uncharged DHICA-eumelanins were computed. The RDFs for systems of two DHICA-eumelanins are shown in black. As Figure 3a shows, in the system of two DHICA-eumelanins the ICAQ moieties are at distances of 6 and 10.6 Å. The DHICA moieties of two DHICA-eumelanins, however, are in the range of π − π interactions, with the first peak at 3.8 Å followed by peaks at further distances of 6 Å and 6.8 Å; see Figure 3b. The PTCA moieties have the closest peak in the range of π − π interactions, at 3.6 Å, Figure 3c. Classical MD simulations do not contain electrons, and hence π − π interactions have only an operational definition. It is typically based on the intermolecular distance and sometimes supplemented by an angle defined between the centroids (here, COMs) of the two rings [83]. The typical distance criterion is 3.2–4.0 Å and the centroid–centroid distance is then slightly longer than the intermolecular distance, as the rings are slightly shifted [83] (parallel-displaced configuration [84]); see the ring orientations in Figure 2.
The RDFs for the systems of four DHICA-eumelanins are shown in red in Figure 3. Because these systems are larger and have more structural flexibility than the systems with two DHICA-eumelanins, the RDFs show more structure. The first and dominant distances are generally further than in the systems of two DHICA-eumelanins, and only the PTCA moieties, shown in Figure 3c, have their first peak in the range of π − π interactions at 3.7 Å. The dominant peak, however, is at 7.2 Å. The shift to larger distances continues with the largest system of 50 neutral DHICA-eumelanins; similarly to the case of four DHICA-eumelanins, only the case of PTCA shows a peak in the range of π − π interactions at 3.7 Å, with the dominant peak at 4.3 Å.
We investigated the RDFs between the COMs of the ICAQ and the PTCA moieties (Figure 4a), the DHICA and the PTCA moieties (Figure 4b), and the ICAQ and the DHICA moieties (Figure 4c). All the systems show multiple peaks with the primary peak being mostly in the range of about 4.5–5.2 Å. In all cases, there is a minor structure around 4 Å which is in the range of π − π interactions.
In summary, the RDFs revealed the possibility of π − π interactions. The other mechanism that promotes aggregation is hydrogen bonding between the carboxylic acids of two neighboring uncharged DHICA-eumelanins. The DHICA-eumelanin structure is not rigid. Apart from being able to rotate around each of the bonds linking the planes, each planar moiety may bend slightly. Due to the twisted structure, at most only two aromatic rings are able to interact via π − π interactions, while the remaining planes cannot form π − π interactions with the same neighboring eumelanin; these may interact via H-bonds or van der Waals interactions with other DHICA-eumelanins or water molecules. These situations are shown in Figure 2a,b with two and four eumelanins.

3.2. Aggregation

Aggregation of eumelanins in aqueous solutions was studied under conditions corresponding to pH = 7.4 using systems of two, four, and 50 charged DHICA-eumelanins. The details of the respective systems are provided in Table 1. The charged (protonated) eumelanin has four carboxylic acid groups (Figure 1a), which have been indicated as potential binding sites for drugs [85,86]. No aggregation occurred in the charged systems, as shown in Figure 5. At higher concentrations, however, pairing of eumelanins via π − π interactions was observed, as shown in Figure 5c.
Next, the distances between the counterions (K + ) and the oxygen atoms in the system of four charged DHICA-eumelanins were examined in the time window of 1000–1100 ns. The selected distance distributions (the probability density functions) are shown in Figure 6. The results show that the K + ions favour interactions with the oxygen atoms of PTCA. In addition, the first shell between the K + ions and the oxygen atoms of the PTCA moieties is at 2.7–2.8 Å, which is in the same range (2.73–2.79 Å) as the first shell of the water oxygens and the K + ions reported in experiments [87]. Thus, K + ions interact with water molecules as well as with DHICA-eumelanins.
Snapshots from a mixture of 25 uncharged and 25 charged DHICA eumelanins after about 1900 ns are shown in Figure 7. As the figure shows, the uncharged eumelanins (green) become shielded by the charged ones (pink). As Figure 7a shows, a number of the charged molecules and counterions remain dispersed in the solution. Figure 7b shows a zoomed-in view displaying aggregated uncharged eumelanins and charged eumelanins that are within 3 Å of uncharged eumelanins as well as the K + ions that are within 3 Å of charged eumelanins. The K + ions that bind with carboxylic acids are shown as dots in order to identify them better. During the aggregation process, one of the charged eumelanins became trapped inside the cluster. The figure shows that due to hydrophobic sites in the uncharged molecules, they prefer to interact with themselves rather than with water molecules. The charged eumelanins, however, tend to form hydrogen bonds with water molecules or bind with the K + ions. In addition to the above, a few of the charged eumelanins are in contact with uncharged eumelanins via π − π interactions.
The RDFs for the mixed systems can now be re-examined with the above behaviour in mind. The green lines in Figure 3 and Figure 4 show the mixed systems. Comparison of the RDFs from the system of 50 neutral DHICA-eumelanins (blue lines in Figure 3 and Figure 4) and the snapshots in Figure 7 makes it clear that the uncharged DHICA-eumelanins are compressed by the charged DHICA-eumelanins in the mixed systems.

3.3. Dihedral Angles

The DHICA-eumelanin protomolecule has a twisted structure, as is clear from Figure 4. Previous quantum computations [11,36,65] have shown that the dihedral angle between the DHICA (5,6-dihydroxyindole-2-carboxylic acid) and ICAQ moieties (see Figure 1a for the definitions of the groups) is in the range of 47–63.1 ∘ . Our DFT calculations in the gas phase resulted in 60.3 ∘ , which is in good agreement with Pezzella’s results [36], though higher than Matta’s [11] DFT result of 50 ∘ for the DHICA dimer (see Section 2.1 Geometry optimization and validation for more details).
The dihedral angle distributions from the classical MD simulations in explicit water are shown in Figure 8. The dihedral angle between the DHICA group and its tautomer, ICAQ, is defined by the C6-C7-C4 ′ -C5 ′ atoms (Figure 1a). The results are shown in Figure 8a. In the uncharged systems (solid lines) of the two protomolecules the most probable angle is around 60 ∘ , with a second peak present around 120 ∘ . The amplitude of the first peak decreases as the number of molecules increases. The second peak displays opposite behaviour. The main peak is around 60 ∘ , while the second peak is shifted slightly to around 110 ∘ . The result of 60 ∘ is in very good agreement with the experiments of Corani et al. [88].
Considering that there are three twisted planes, the next dihedral angle to investigate is the one between the DHICA and PTCA moieties (Figure 1a). The results are shown in Figure 8b. The geometries of individual uncharged and charged DHICA-eumelanins are presented as aligned structures in Figure 9a. The dihedral angles for the uncharged systems vary between − 50 ∘ to − 145 ∘ , with a peak at − 110 ∘ . For the charged systems, the dihedral angles are in the range of 50–130 ∘ , with a peak at around 100 ∘ .
To better understand the results in Figure 8, Figure 9a shows the ICAQ-DHICA (C6-C7-C4 ′ -C5 ′ atoms) dihedral angle in blue. For the DHICA-PTCA angle (C8 ′ -C7 ′ -C-C3 ″ atoms), black is used for the uncharged and red for the charged DHICA-eumelanin. Snapshots of 100 superimposed structures from the system with four uncharged DHICA-eumelanins are shown in Figure 9b, and Figure 9c shows the same for the corresponding charged system. The main difference in their behaviour is that the charged DHICA-eumelanins are soluble in water and that the PTCA moiety is in favor of interacting with the K + ions. This is likely the the reason for the DHICA-PTCA dihedral angle distribution being narrower in the charged systems.

3.4. Hydrogen Bonding

One of the most key molecular interactions is hydrogen bonding. We computed H-bonds based on a donor–acceptor distance of <3.5 Å and hydrogen donor–acceptor angle of < 25 ∘ . Figure 10a shows the average number of H-bonds between the DHICA-eumelanins and water per DHICA-eumelanin molecule. In the small systems of two and four molecules (the upper panels), the number fluctuates, while in the two larger systems (lower panels) it settles to about 11.5. The number of H-bonds between the DHICA-eumelanins settles to approximately two after 60 ns (Figure 10b).
The H-bonding properties are a consequence of the complex twisted structure, as is apparent from the analysis of the dihedral angles and the orientations of the different moieties with respect to each other. Certain potentially favourable interactions are not possible due to the steric constraints and the orientations of the different groups. The complexity of the situation is clear from Figure 2. Figure 2a,b shows the non-planar structures; as Figure 2d shows, water is trapped within the complex.
We analyzed all the H-bonds in the system; one special feature is that the ICAQ group forms on average only about 0.07 intermolecular H-bonds with the neighboring molecules’ ICAQ groups, while all the other intermolecular combinations have either around 0.3 (PTCA-PTCA, DHICA-DHICA and PTCA-ICAQ) or even higher (>0.5) H-bonds per moiety (ICAQ-DHICA and PTCA-DHICA). That is, the DHICA is the key part for intermolecular H-bonding. Regarding the about two H-bonds/molecule (intermolecular), that number is enough to provide the aggregates with enough stability against thermal fluctuations to allow them to form.
As for intramolecular bonding, the dominant bond is between PTCA and ICAQ, although its contribution is very small; over the simulation time, the number of bonds fluctuates between 0.02 and 0.06.
The presence of water is critical in these systems. For example, when compared to the DFT calculations of Matta et al. [11], which were carried out in the absence of water, intramolecular H-bond networks form in either a zig-zag or helical fashion. In addition, their systems consisted of DHICA moieties only. The case here is very different; due to the presence of water, no such networks appear.

3.5. Dipole Moment

The dipole moment for the uncharged DHICA-eumelanin model was calculated in the gas phase with the B3LYP/6-311G(d) basis sets and ESP fitting [69,70] following the Merz-Singh-Kollman (MK) scheme [70,71] and a 500 self-consistent field cycles. The calculation yielded 11.0 Debye. For comparison, the same system was optimized with the ω B97XD/6-311G(d,p) basis sets without ESP fitting and 500 self-consistent field cycles. The result was 11.5 Debye. Both methods provided higher values than the results for dimeric DHICA reported by Matta et al. [11], who obtained 5.8 Debye. However, they obtained 12.7 Debye for tetrameric DHICA (four carboxylic acids). Our result for the trimeric form is between these two values. In addition, differences originate from the carbonyl groups in ICAQ and from the PTCA moiety due to polarity of the carboxylic acids.

4. Conclusions

We studied the aggregation of DHICA-eumelanin using a combination of multi-microsecond MD simulations and DFT calculations. The latter were performed to parameterize and verify the model and to determine the dipole moments. Classical MD simulations were performed at different concentrations of uncharged DHICA-eumelanins, charged DHICA-eumelanins (the protonation state corresponding to about pH 7.4), and a mixture of the two. The DFT calculation of the dihedral angle of uncharged ICAQ-DHICA moieties in the gas phase resulted in approximately 60 ∘ , which is in good agreement with the prior DFT studies of Pezzella et al., who used negatively charged ICAQ-DHICA moieties [36]. The dipole moment of uncharged DHICA-eumelanin in the gas phase was carried out at B3LYP/6-311G(d) and ω B97XD/6-311G(d,p) levels. The calculations resulted in 11–11.5 Debye. This value is above the DHICA-dimer dipole moment of 5.8 Debye and below the tetramer-DHICA value of 12.7 Debye reported by Matta et al. [11] using DFT calculations at ω B97XD-D/6-31G+ level of theory.
In the classical MD simulation with explicit water, the uncharged DHICA-eumelanins were aggregated to bundles. Due to their twisted structure, only minor or no π − π interactions were present. Instead, hydrogen bonds formed. Because the interactions between the molecules were relatively weak, water was trapped inside the aggregates at higher concentrations. This is very different from the aggregation of DHI-eumelanins, which form well-defined tight stacks with no water inside the aggregates [28,45].
The simulation results show that, upon aggregation, the average number of H-bonds between the individual DHICA-eumelanins settles to about two H-bonds/molecule. The number of H-bonds between the DHICA-eumelanins and water plateaus at around 11/DHICA-eumelanin. The dihedral angle distribution between the DHICA and ICAQ moieties in water solution is semi-bimodal, with a peak around 60 ∘ . This is in very good agreement with the experiments of Corani et al. [88] (≈ 60 ∘ ), who studied 4-7 ′ DHICA moiety dimers in water, though higher than Matta’s [11] results for oligomers of DHICA moieties 50 ∘ using ω B97X-D/cc-pVTZ level of theory.
All carboxylic acids in DHICA-protomolecules become ionized at physiological pH, and the protomolecule attains a net charge of − 4 . The charged protomolecules favour interactions with water and potassium ions rather than with themselves. Distributions of K + ions around the DHICA-eumelanins show the first shell at ≈2.7 Å, which is in the same range as the experimentally measured potassium–water (using water’s oxygen) first shell at 2.73–2.79 Å [87]. The reason that PTCA moieties contribute to ion binding more than DHICA or ICAQ moieties could arise from the fact that the dihedral angle distribution between PTCA and DHICA moieties for charged DHICA-eumelanins is around 110 ∘ . However, the DHICA and ICAQ moieties are closer together, and the dihedral angles are ≈ 65 ∘ to minimize electrostatic interactions, which is slightly higher than in the case of uncharged DHICA-eumelanin.
The systems of mixed charged and uncharged DHICA-eumelanins show aggregation of uncharged eumelanin, with charged DHICA-eumelanins H-bonding to them. Potassium ions bind to the carboxylic acids of the charged molecules. The results of this study show that DHICA-eumelanins have structural disorder, as well as that charged DHICA-eumelanins, unlike the uncharged ones, are soluble in an aqueous solution with potassium ions. Potassium ions prevent fast aggregation. This could potentially allow for investigations of absorption profiles of natural eumelanin without scattering effects. In addition, the experimentally observed aggregates in naturally occurring eumelanin are most likely composed of uncharged DHICA, with charged molecules H-bonding to them at the surfaces of the aggregates. Finally, the bundled and disorderd shapes of the aggregates suggest that structural disorder should be considered in models for drug–eumelanin binding.

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/molecules27238417/s1. Figure S1: DHICA monomer structure. Figure S2: Geometry optimized DHICA-eumelanin chemical structure (ICAQ-DHICA-PTCA). Table S1: Atom types and coordinates of monomer DHICA’s Dreyer structure (in  Å). Table S2: Bond lengths of monomer DHICA’s Dreyer structure (in  Å) compared with Okuda’s [54] and Powell’s [48] results. Table S3: Comparison of the angles of the DHICA monomer with Okuda’s [54] and Powell’s [48] results. Table S4: Calculated molecular vibrational intensity for DHICA monomer. Table S5: Calculated molecular vibrational intensity for DHICA monomer of Okuda’s results for comparison with Table S4. Table S6: Atom types and coordinates of optimized DHICA-eumelanin (in  Å). Table S7: Comparison of bond lengths of DHICA-eumelanin calculated using DFT and MD simulations with explicit water model. Table S8: Comparison of angle degrees of DHICA-eumelanin calculated using DFT and MD simulations with explicit water model.

Author Contributions

Conceptualization, M.K.; methodology, validation, formal analysis, investigation, and data curation, S.S. and M.K.; resources, M.K.; writing—original draft preparation, S.S.; writing—review and editing, S.S., S.S.-H., C.G.T.F. and M.K.; visualization, S.S.; supervision, M.K.; funding acquisition, M.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC) and Canada Research Chairs Program.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Parameters for the DHICA-eumelanin protomolecule are available at https://github.com/SoftSimu/melanin (accessed on 15 September 2022) including Gromacs compatible gro and itp files.

Acknowledgments

MK acknowledges financial support by the Natural Sciences and Engineering Research Council of Canada (NSERC) and Canada Research Chairs Program. Computational resources were provided by Compute Canada and SharcNet.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
DHICA5,6-DiHydroxyIndole-2-Carboxylic Acid
DHI5,6-DiHydroxyIndole
PTCAPyrrole-2,3,5-TriCarboxylic Acid
DFTDensity Functional Theory
STMscanning tunneling microscopy
UHR-SEMultra-high resolution scanning electron microscopy
MDMolecular Dynamics
OPLS-AAOptimized Potentials for Liquid Simulations for All Atoms
SPCSimple Point Charge
PMEParticle-Mesh Ewald
P-LINCSParallel LINear Constraint Solver
NVTconstant particle Number, Volume, and Temperature
NPTconstant particle Number, Pressure, and Temperature
VMDVisual Molecular Dynamics
RDFRadial Distribution Functions
COMCenter of Mass
ESPElectroStatic Potential
B3LYPBecke three-parameter Lee–Yang–Parr
MKMerz–Singh–Kollman

References

  1. d’Ischia, M.; Napolitano, A.; Pezzella, A.; Meredith, P.; Buehler, M. Melanin Biopolymers: Tailoring Chemical Complexity for Materials Design. Angew. Chem. Int. Ed. 2020, 59, 11196–11205. [Google Scholar] [CrossRef] [PubMed]
  2. d’Ischia, M.; Napolitano, A.; Pezzella, A.; Meredith, P.; Sarna, T. Chemical and Structural Diversity in Eumelanins: Unexplored Bio-Optoelectronic Materials. Angew. Chem. Int. Ed. 2009, 48, 3914–3921. [Google Scholar] [CrossRef] [PubMed]
  3. Panzella, L.; Gentile, G.; D’Errico, G.; Della Vecchia, N.F.; Errico, M.E.; Napolitano, A.; Carfagna, C.; d’Ischia, M. Atypical Structural and π-Electron Features of a Melanin Polymer That Lead to Superior Free-Radical-Scavenging Properties. Angew. Chem. Int. Ed. 2013, 52, 12684–12687. [Google Scholar] [CrossRef] [PubMed]
  4. Schraermeyer, U.; Heimann, K. Current understanding on the role of retinal pigment epithelium and its pigmentation. Pigment Cell Res. 1999, 12, 219–236. [Google Scholar] [CrossRef] [PubMed]
  5. Ortonne, J.P. Photoprotective properties of skin melanin. Br. J. Dermatol. 2002, s61, 7–10. [Google Scholar] [CrossRef]
  6. Sarna, F. Properties and function of the ocular melanin–a photobiophysical view. J. Photochem. Photobiol. B 1992, 12, 215–258. [Google Scholar] [CrossRef]
  7. Ju, K.Y.; Kang, J.; Chang, J.H.; Lee, J.K. Clue to Understanding the Janus Behavior of Eumelanin: Investigating the Relationship between Hierarchical Assembly Structure of Eumelanin and Its Photophysical Properties. Biomacromolecules 2016, 17, 2860–2872. [Google Scholar] [CrossRef]
  8. Micillo, R.; Iacomino, M.; Perfetti, M.; Panzella, L.; Koike, K.; D’Errico, G.; d’Ischia, M.; Napolitano, A. Unexpected impact of esterification on the antioxidant activity and (photo)stability of a eumelanin from 5,6-dihydroxyindole- 2-carboxylic acid. Pigment Cell Melanoma Res. 2018, 31, 475–483. [Google Scholar] [CrossRef]
  9. Liberti, D.; Alfieri, M.L.; Monti, D.M.; Panzella, L.; Napolitano, A. A Melanin-Related Phenolic Polymer with Potent Photoprotective and Antioxidant Activities for Dermo-Cosmetic Applications. Antioxidants 2018, 9, 270. [Google Scholar] [CrossRef] [Green Version]
  10. Wünsche, J.; Deng, Y.; Kumar, P.; Mauro, E.D.; Josberger, E.; Sayago, J.; Pezzella, A.; Soavi, F.; Cicoira, F.; Rolandi, M.; et al. Protonic and Electronic Transport in Hydrated Thin Films of the Pigment Eumelanin. Chem. Mater. 2015, 27, 436–442. [Google Scholar] [CrossRef]
  11. Matta, M.; Pezzella, A.; Troisi, A. Relation between Local Structure, Electric Dipole, and Charge Carrier Dynamics in DHICA Melanin: A Model for Biocompatible Semiconductors. J. Phys. Chem. Lett. 2020, 11, 1045–1051. [Google Scholar] [CrossRef]
  12. Tran, M.L.; Powell, B.J.; Meredith, P. Chemical and Structural Disorder in Eumelanins: A Possible Explanation for Broadband Absorbance. Biophys. J. 2006, 90, 743–752. [Google Scholar] [CrossRef] [Green Version]
  13. Abbas, M.; D’Amico, F.; Morresi, L.; Pinto, N.; Ficcadenti, M.; Natali, R.; Ottaviano, L.; Passacantando, M.; Cuccioloni, M.; Angeletti, M.; et al. Structural, electrical, electronic and optical properties of melanin films. Eur. Phys. J. E 2009, 28, 285–291. [Google Scholar] [CrossRef]
  14. d’Ischia, M.; Napolitano, A.; Ball, V.; Chen, C.T.; Buehler, M.J. Polydopamine and Eumelanin: From Structure–Property Relationships to a Unified Tailoring Strategy. Acc. Chem. Res. 2014, 47, 3541–3550. [Google Scholar] [CrossRef]
  15. Mostert, A.B. Melanin, the What, the Why and the How: An Introductory Review for Materials Scientists Interested in Flexible and Versatile Polymers. Polymers 2021, 13, 1670. [Google Scholar] [CrossRef]
  16. Cheng, J.; Moss, S.; Eisner, M.; Zschack, P. X-ray Characterization of Melanins—I. Pigment Cell Res. 1994, 7, 255–262. [Google Scholar] [CrossRef]
  17. Cheng, J.; Moss, S.; Eisner, M. X-ray Characterization of Melanins—II. Pigment Cell Res. 1994, 7, 263–273. [Google Scholar] [CrossRef]
  18. Zajac, G.; Gallas, J.; Cheng, J.; Eisner, M.; Moss, S.; Alvarado-Swaisgood, A. The fundamental unit of synthetic melanina verification by tunneling microscopy of X-ray-scattering results. Biochim. Biophys. Acta 1994, 1199, 271–278. [Google Scholar] [CrossRef]
  19. Zajac, G.; Gallas, J.; Alvarado-Swaisgood, A. Tunneling microscopy verification of an x-ray scattering-derived molecular model of tyrosine-based melanin. J. Vac. Sci. Technol. B. 1994, 12, 1512–1516. [Google Scholar] [CrossRef]
  20. Liu, Y.; Simon, J. The effect of preparation procedures on the morphology of melanin from the ink sac of Sepia officinalis. Pigment Cell Res. 2003, 16, 72–80. [Google Scholar] [CrossRef]
  21. Chen, C.T.; Martin-Martinez, F.J.; Jung, G.S.; Buehler, M.J. Polydopamine and eumelanin molecular structures investigated with ab initio calculations. Chem. Sci. 2017, 8, 1631–1641. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Reilly, J.; Williams, S.L.; Forster, C.J.; Kansara, V.; End, P.; Serrano-Wu, M.H. High-throughput melanin-binding affinity and in silico methods to aid in the prediction of drug exposure in ocular tissue. J. Pharm. Sci. 2015, 104, 3997–4001. [Google Scholar] [CrossRef] [PubMed]
  23. Karlsson, O.; Lindquist, N.G. Melanin affinity and its possible role in neurodegeneration. J. Neural Transm. 2013, 120, 1623–1630. [Google Scholar] [CrossRef] [PubMed]
  24. Karlsson, O.; Lindquist, N.G. Melanin and neuromelanin binding of drugs and chemicals: Toxicological implications. Arch. Toxicol. 2016, 90, 1883–1891. [Google Scholar] [CrossRef] [PubMed]
  25. Nofsinger, J.; Forest, S.; Eibest, L.; Gold, K.; Simon, J. Probing the building blocks of eumelanins using scanning electron microscopy. Pigment Cell Res. 2000, 13, 179–184. [Google Scholar] [CrossRef]
  26. Watt, A.A.R.; Bothma, J.P.; Meredith, P. The supramolecular structure of melanin. Soft Matter. 2009, 5, 3754. [Google Scholar] [CrossRef]
  27. Pezzella, A.; Panzella, L.; Crescenzi, O.; Napolitano, A.; Navaratman, S.; Edge, R.; Land, E.J.; Barone, V.; d’Ischia, M. 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. [Google Scholar] [CrossRef]
  28. Chen, C.T.; Ball, V.; de Almeida Gracio, J.J.; Singh, M.K.; Toniazzo, V.; Ruch, D.; Buehler, M.J. 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. [Google Scholar] [CrossRef]
  29. Micillo, R.; Panzella, L.; Koike, K.; Monfrecola, G.; Napolitano, A.; d’Ischia, M. “Fifty Shades” of Black and Red or How Carboxyl Groups Fine Tune Eumelanin and Pheomelanin Properties. Int. J. Mol. Sci. 2016, 17, 746. [Google Scholar] [CrossRef]
  30. Jakubiak, P.; Lack, F.; Thun, J.; Urtti, A.; Alvarez-Sánchez, R. Influence of Melanin Characteristics on Drug Binding Properties. Mol. Pharm. 2019, 16, 2549–2556. [Google Scholar] [CrossRef]
  31. Pezzella, A.; d’Ischia, M.; Napolitano, A.; Palumbo, A.; Prota, G. An integrated approach to the structure of Sepia melanin. Evidence for a high proportion of degraded 5,6-dihydroxyindole-2-carboxylic acid units in the pigment backbone. Tetrahedron 1997, 53, 8281–8286. [Google Scholar] [CrossRef]
  32. Magarelli, M.; Passamonti, P.; Renieri, C. Purification, characterization and analysis of sepia melanin from commercial sepia ink (Sepia Officinalis). Rev. CES Med. Vet. Zootec. 2010, 5, 18–28. [Google Scholar] [CrossRef]
  33. Liu, Y.; Hong, L.; Bowers, C.; Wakamatsu, K.; Ito, S.; Simon, J. chemical and spectroscopic properties of human black and red hair melanosomes. Photochem. Photobiol. 2005, 81, 135–144. [Google Scholar] [CrossRef]
  34. Terranovaa, M.L.; Tamburri, E. Understanding the way eumelanin works: A unique example of properties and skills driven by molecular heterogeneity. Polymer 2021, 229, 123952. [Google Scholar] [CrossRef]
  35. Supakar, S.; Banerjee, A.; Jha, T. Intermolecular association of some selected melanin monomers and their optical absorption. Comput. Theor. Chem. 2019, 1151, 43–49. [Google Scholar] [CrossRef]
  36. Pezzella, A.; Panzella, L.; Crescenzi, O.; Napolitano, A.; Navaratnam, S.; Edge, R.; Land, E.J.; Barone, V.; d’Ischia, M. 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. [Google Scholar] [CrossRef]
  37. Panzella, L.; Napolitano, A.; d’Ischia, M. Is DHICA the key to dopachrome tautomerase and melanocyte functions? Pigment Cell Melanoma Res. 2011, 24, 248–249. [Google Scholar] [CrossRef]
  38. Jiang, S.; Liu, X.M.; Dai, X.; Zhou, Q.; Lei, T.C.; Beermann, F.; Wakamatsu, K.; Xu, S. Regulation of DHICA-mediated antioxidation by dopachrome tautomerase: Implication for skin photoprotection against UVA radiation. Free Radic. Biol. Med. 2010, 48, 1144–1151. [Google Scholar] [CrossRef]
  39. Cecchi, T.; Pezzella, A.; Di Mauro, E.; Cestola, S.; Ginsburg, D.; Luzi, M.; Rigucci, A.; Santato, C. On the antioxidant activity of eumelanin biopigments: A quantitative comparison between free radical scavenging and redox properties. Nat. Prod. Res. 2020, 34, 2465–2473. [Google Scholar] [CrossRef]
  40. Kaxiras, E.; Tsolakidis, A.; Zonios, G.; Meng, S. Structural model of eumelanin. Phys. Rev. Lett. 2006, 97, 218102. [Google Scholar] [CrossRef]
  41. Pezzella, A.; Iadonisi, A.; Valerio, S.; Panzella, L.; Napolitano, A.; Adinolfi, M.; d’Ischia, M. 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. [Google Scholar] [CrossRef] [PubMed]
  42. Chen, C.T.; Chuang, C.; Cao, J.; Ball, V.; Ruch, D.; Buehler, M.J. Excitonic effects from geometric order and disorder explain broadband optical absorption in eumelanin. Nat. Commun. 2014, 5, 1. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Assis Oliveira, L.B.; L Fonseca, T.; Costa Cabral, B.J.; Coutinho, K.; Canuto, S. Hydration Effects on the Electronic Properties of Eumelanin Building Blocks. J. Chem. Phys. 2016, 145, 084501. [Google Scholar] [CrossRef] [PubMed]
  44. Chen, C.T.; Buehler, M.J. Polydopamine and eumelanin models in various oxidation states. Phys. Chem. Chem. Phys. 2018, 20, 28135–28143. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Soltani, S.; Sowlati-Hashjin, S.; Feugmo, C.G.T.; Karttunen, M. Free energy and stacking of eumelanin nanoaggregates. J. Phys. Chem. B. 2022, 22, 563–579. [Google Scholar] [CrossRef] [PubMed]
  46. Choudhury, A.; Ghosh, D. Charge transfer in DHICA eumelanin-like oligomers: Role of hydrogen bonds. Chem. Commun. 2020, 56, 10481–10484. [Google Scholar] [CrossRef]
  47. Micillo, R.; Panzella, L.; Iacomino, M.; Prampolini, G.; Cacelli, I.; Ferretti, A.; Crescenzi, O.; Koike, K.; Napolitano, A.; d’Ischia, M. Eumelanin broadband absorption develops from aggregation-modulated chromophore interactions under structural and redox control. Sci. Rep. 2017, 7, 1. [Google Scholar] [CrossRef] [Green Version]
  48. Powell, B. 5,6-dihydroxyindole-2-carboxylic acid (DHICA): A First Principles Density-Functional Study. Chem. Phys. Lett. 2005, 402, 111–115. [Google Scholar] [CrossRef] [Green Version]
  49. Meng, S.; Kaxiras, E. Theoretical Models of Eumelanin Protomolecules and Their Optical Properties. Biophys. J. 2008, 94, 2095–2105. [Google Scholar] [CrossRef] [Green Version]
  50. Wakamatsu, K.; Ito, S. Advanced Chemical Methods in Melanin Determination. Pigment Cell Res. 2002, 15, 174–183. [Google Scholar] [CrossRef]
  51. Ito, S. Reexamination of the structure of eumelanin. J. Biol. Chem. 1986, 883, 155–161. [Google Scholar] [CrossRef]
  52. Schroeder, R.L.; Pendleton, P.; Gerber, J. Physical factors affecting chloroquine binding to melanin. Colloids Surf. B Biointerfaces 2015, 134, 8–16. [Google Scholar] [CrossRef]
  53. Kirla, K.T.; Groh, K.J.; Poetzsch, M.; Banote, R.K.; Stadnicka-Michalak, J.; Eggen, R.I.L.; Schirmer, K.; Kraemer, T. Importance of Toxicokinetics to Assess the Utility of Zebrafish Larvae as Model for Psychoactive Drug Screening Using Meta-Chlorophenylpiperazine (mCPP) as Example. Front. Pharmacol. 2018, 9, 414. [Google Scholar] [CrossRef] [Green Version]
  54. Okuda, H.; Nakamura, A.; Wakamatsu, K.; Hanson, G.; Ito, S.; Sota, T. Mid-infrared absorption spectrum of 5,6-dihydroxyindole-2-carboxylic acid. Chem. Phys. Lett. 2007, 433, 355–359. [Google Scholar] [CrossRef]
  55. Becke, A.D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A. 1988, 38, 3098–3100. [Google Scholar] [CrossRef]
  56. Lee, C.; Yang, W.; Parr, R.G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. A. 1988, 37, 785–789. [Google Scholar] [CrossRef] [Green Version]
  57. Becke, A.D. A new mixing of Hartree–Fock and local density-functional theories. J. Chem. Phys. 1993, 98, 1372–1377. [Google Scholar] [CrossRef]
  58. Becke, A.D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. [Google Scholar] [CrossRef] [Green Version]
  59. Stephens, P.J.; Devlin, F.J.; Chabalowski, C.F.; Frisch, M.J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623–11627. [Google Scholar] [CrossRef]
  60. Kim, K.; Jordan, K.D. Comparison of Density Functional and MP2 Calculations on the Water Monomer and Dimer. J. Phys. Chem. 1994, 98, 10089–10094. [Google Scholar] [CrossRef]
  61. Hehre, W.J.; Ditchfield, R.; Pople, J.A. Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257–5261. [Google Scholar] [CrossRef]
  62. Francl, M.M.; Pietro, W.J.; Hehre, W.; Binkley, J.S.; Gordon, M.S.; DeFrees, D.J.; Pople, J.A. Selfconsistent molecular orbital methods. XXIII. A polarizationtype basis set for secondrow elements. J. Chem. Phys. 1982, 77, 3654–3665. [Google Scholar] [CrossRef] [Green Version]
  63. Davidson, E.R.; Feller, D. Basis set selection for molecular calculations. Chem. Rev. 1989, 86, 681–696. [Google Scholar] [CrossRef]
  64. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Petersson, G.A.; Nakatsuji, H.; et al. Gaussian 09 Revision E.01; Gaussian Inc.: Wallingford, CT, USA, 2016. [Google Scholar]
  65. Pezzella, A.; Crescenzi, O.; Panzella, L.; Napolitano, A.; Land, E.J.; Barone, V.; d’Ischia, M. Free Radical Coupling of o-Semiquinones Uncovered. J. Am. Chem. Soc. 2013, 135, 12142–12149. [Google Scholar] [CrossRef] [PubMed]
  66. Dodda, L.; Cabeza de Vaca, I.; Tirado-Rives, J.; Jorgensen, W. LigParGen web server: An automatic OPLS-AA parameter generator for organic ligands. Nucleic Acids Res. 2017, 45, W331–W336. [Google Scholar] [CrossRef] [Green Version]
  67. Jorgensen, W.; Tirado-Rives, J. The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin? J. Am. Chem. Soc. 1988, 110, 1657–1666. [Google Scholar] [CrossRef]
  68. Robertson, M.; Tirado-Rives, J.; Jorgensen, W. Improved Peptide and Protein Torsional Energetics with the OPLS-AA Force Field. J. Chem. Theory Comput. 2015, 11, 3499–3509. [Google Scholar] [CrossRef] [Green Version]
  69. Cox, S.R.; Williams, D.E. Representation of the molecular electrostatic potential by a net atomic charge model. J. Comput. Chem. 1981, 2, 304–323. [Google Scholar] [CrossRef]
  70. Singh, U.C.; Kollman, P.A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 1984, 5, 129–145. [Google Scholar] [CrossRef]
  71. Besler, B.; KM, M.J.; Kollman, P. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11, 431–439. [Google Scholar] [CrossRef]
  72. Kutzner, C.; Páll, S.; Fechner, M.; Esztermann, A.; de Groot, B.L.; Grubmüller, H. More Bang for Your Buck: Improved use of GPU Nodes for GROMACS 2018. J. Comp. Chem. 2019, 40, 2418–2431. [Google Scholar] [CrossRef] [Green Version]
  73. Berendsen, H.; Postma, J.; van Gunsteren, W.; Hermans, J. Intermolecular Forces; Pullman, B., Ed.; Springer: Dordrecht, The Netherlands, 1981. [Google Scholar]
  74. Essmann, U.; Perera, L.; Berkowitz, M.; Darden, T.; Lee, H.; Pedersen, L. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577–8593. [Google Scholar] [CrossRef] [Green Version]
  75. Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N.Log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089–10092. [Google Scholar] [CrossRef] [Green Version]
  76. Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. [Google Scholar] [CrossRef] [Green Version]
  77. Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116–122. [Google Scholar] [CrossRef]
  78. Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182–7190. [Google Scholar] [CrossRef]
  79. Gapsys, V.; de Groot, B.L. Comment on ‘Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size’. eLife 2019, 8, e44718. [Google Scholar] [CrossRef]
  80. Gapsys, V.; de Groot, B.L. On the Importance of Statistics in Molecular Simulations for Thermodynamics, Kinetics and Simulation Box Size. eLife 2020, 8, e45318. [Google Scholar] [CrossRef]
  81. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33–38. [Google Scholar] [CrossRef]
  82. DeLano, W.L. The PyMOL Molecular Graphics System; Version 1.8; DeLano Scientific LLC: South San Francisco, CA, USA, 2015. [Google Scholar]
  83. Liu, Y.; Ao, X.; Wang, Q.; Wang, J.; Ge, H. PiViewer: An open-source tool for automated detection and display of π-π interactions. Chem. Biol. Drug Des. 2018, 92, 1809–1814. [Google Scholar] [CrossRef]
  84. Carter-Fenk, K.; Herbert, J.M. Reinterpreting π-stacking. Phys. Chem. Chem. Phys. 2020, 22, 24870–24886. [Google Scholar] [CrossRef] [PubMed]
  85. Hong, L.; Simon, J.D. Physical and chemical characterization of iris and choroid melanosomes isolated from newborn and mature cows. Photochem. Photobiol. 2005, 81, 517–523. [Google Scholar] [CrossRef] [PubMed]
  86. Hong, L.; Simon, J.D. Current Understanding of the Binding Sites, Capacity, Affinity, and Biological Significance of Metals in Melanin. J. Phys. Chem. B. 2007, 111, 7938–7947. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. Dang, L.X.; Schenter, G.K.; Glezakou, V.; Fulton, J.L. Molecular Simulation Analysis and X-ray Absorption Measurement of Ca2+, K+ and Cl- Ions in Solution. J. Phys. Chem. B 2006, 110, 23644–23654. [Google Scholar] [CrossRef]
  88. Corani, A.; Huijser, A.; Gustavsson, T.; Markovitsi, D.; Malmqvist, P.; Pezzella, A.; d’Ischia, M.; Sundström, V. Superior Photoprotective Motifs and Mechanisms in Eumelanins Uncovered. J. Am. Chem. Soc. 2014, 136, 11626–11635. [Google Scholar] [CrossRef]
Molecules 27 08417 g001 550
Figure 1. (a) The chemical structure of the uncharged eumelanin model (DHICA-eumelanin): the ICAQ (DHICA’s tautomer), 5,6-dihydroxyindole-2-carboxylic acid (DHICA), and pyrrole-2,3,5-tricarboxylic acid (PTCA) moeities. In addition to uncharged DHICA-eumelanin, deprotonated DHICA-eumelanin corresponding to physiological pH was simulated. Due to the PKa value of the DHICA moiety (b), all four carboxylic groups deprotonate, and the molecule obtains a net charge of − 4 [24,52,53].
Figure 1. (a) The chemical structure of the uncharged eumelanin model (DHICA-eumelanin): the ICAQ (DHICA’s tautomer), 5,6-dihydroxyindole-2-carboxylic acid (DHICA), and pyrrole-2,3,5-tricarboxylic acid (PTCA) moeities. In addition to uncharged DHICA-eumelanin, deprotonated DHICA-eumelanin corresponding to physiological pH was simulated. Due to the PKa value of the DHICA moiety (b), all four carboxylic groups deprotonate, and the molecule obtains a net charge of − 4 [24,52,53].
Molecules 27 08417 g001
Molecules 27 08417 g002 550
Figure 2. Size dependence of aggregation of uncharged DHICA eumelanin in an aqueous solution. For clarity, water molecules are not shown. (a) Two and (b) four uncharged DHICA molecules after 1073 ns and 983 ns, respectively. Colours in (a,b): blue: nitrogen, red: oxygen, white: hydrogen, and green: carbon. (c) 50 uncharged eumelanins after 2050 ns, carbon atoms are shown in rainbow to distinguish the residues. (d) Same as in (c), except the DHICA-eumelanins are shown in gray and water molecules within the aggregate as spheres.
Figure 2. Size dependence of aggregation of uncharged DHICA eumelanin in an aqueous solution. For clarity, water molecules are not shown. (a) Two and (b) four uncharged DHICA molecules after 1073 ns and 983 ns, respectively. Colours in (a,b): blue: nitrogen, red: oxygen, white: hydrogen, and green: carbon. (c) 50 uncharged eumelanins after 2050 ns, carbon atoms are shown in rainbow to distinguish the residues. (d) Same as in (c), except the DHICA-eumelanins are shown in gray and water molecules within the aggregate as spheres.
Molecules 27 08417 g002
Molecules 27 08417 g003 550
Figure 3. Size dependence of the radial distribution function (RDF). The planes of eumelanin are shown in Figure 1a. RDFs between (a) the COMs of the ICAQ moieties, (b) the DHICA moieties, and (c) the PTCA moieties. (d) RDFs between COMs of the DHICA eumelanin molecules. The systems are two uncharged, four uncharged, 25 uncharged and 25 charged, and fifty uncharged eumelanins; see details in Table 1.
Figure 3. Size dependence of the radial distribution function (RDF). The planes of eumelanin are shown in Figure 1a. RDFs between (a) the COMs of the ICAQ moieties, (b) the DHICA moieties, and (c) the PTCA moieties. (d) RDFs between COMs of the DHICA eumelanin molecules. The systems are two uncharged, four uncharged, 25 uncharged and 25 charged, and fifty uncharged eumelanins; see details in Table 1.
Molecules 27 08417 g003
Molecules 27 08417 g004 550
Figure 4. Size dependence of the radial distribution function (RDF). The planes of eumelanin are shown in Figure 1a. RDFs between the COMs of (a) the ICAQ and PTCA moieties, (b) the DHICA and PTCA moieties, and (c) the ICAQ and DHICA moieties. The systems are two uncharged, four uncharged, 25 uncharged and 25 charged, and fifty uncharged eumelanins; see details in Table 1.
Figure 4. Size dependence of the radial distribution function (RDF). The planes of eumelanin are shown in Figure 1a. RDFs between the COMs of (a) the ICAQ and PTCA moieties, (b) the DHICA and PTCA moieties, and (c) the ICAQ and DHICA moieties. The systems are two uncharged, four uncharged, 25 uncharged and 25 charged, and fifty uncharged eumelanins; see details in Table 1.
Molecules 27 08417 g004
Molecules 27 08417 g005 550
Figure 5. Aggregation of charged DHICA eumelanin molecules. Water molecules are not shown for clarity. Colours in (a,b): blue: nitrogen, red: oxygen, white: hydrogen, green: carbon. The potassium ions are shown in purple. (a) Two charged DHICA eumelanins after 1000 ns, (b) four charged DHICA eumelanins after 1100 ns, and (c) 50 charged DHICA eumelanins after 1361 ns. In (c). rainbow colouring is used to distinguish the residues.
Figure 5. Aggregation of charged DHICA eumelanin molecules. Water molecules are not shown for clarity. Colours in (a,b): blue: nitrogen, red: oxygen, white: hydrogen, green: carbon. The potassium ions are shown in purple. (a) Two charged DHICA eumelanins after 1000 ns, (b) four charged DHICA eumelanins after 1100 ns, and (c) 50 charged DHICA eumelanins after 1361 ns. In (c). rainbow colouring is used to distinguish the residues.
Molecules 27 08417 g005
Molecules 27 08417 g006 550
Figure 6. Distribution of K + - charged DHICA’s oxygen distances for the system of four charged DHICA-eumelanins. O21, O31, and O33 stand for the oxygen atoms of the hydroxyl group in carboxylic anions, and O22, O32, and O34 are the oxygen atoms of the carboxyl group in the carboxylic ions group. O is the oxygen atom of the ketone in PTCA, and O11 was chosen as reference.
Figure 6. Distribution of K + - charged DHICA’s oxygen distances for the system of four charged DHICA-eumelanins. O21, O31, and O33 stand for the oxygen atoms of the hydroxyl group in carboxylic anions, and O22, O32, and O34 are the oxygen atoms of the carboxyl group in the carboxylic ions group. O is the oxygen atom of the ketone in PTCA, and O11 was chosen as reference.
Molecules 27 08417 g006
Molecules 27 08417 g007 550
Figure 7. Mixture of 25 uncharged (green) and 25 charged (pink) DHICA eumelanins. K + ions (purple spheres) were added to the system to neutralize it. (a) Snapshot of the system at 1971 ns and (b) a zoomed-in view showing charged DHICA eumelanins within 3 Å of uncharged ones. Water molecules are not shown for better clarity.
Figure 7. Mixture of 25 uncharged (green) and 25 charged (pink) DHICA eumelanins. K + ions (purple spheres) were added to the system to neutralize it. (a) Snapshot of the system at 1971 ns and (b) a zoomed-in view showing charged DHICA eumelanins within 3 Å of uncharged ones. Water molecules are not shown for better clarity.
Molecules 27 08417 g007
Molecules 27 08417 g008 550
Figure 8. Dihedral angle distributions for all the systems over the last 100 ns: (a) C6-C7-C4 ′ -C5 ′ atoms and (b) C8 ′ -C7 ′ -C-C3 ″ atoms. Solid lines: uncharged molecules. Dotted lines: charged molecules.
Figure 8. Dihedral angle distributions for all the systems over the last 100 ns: (a) C6-C7-C4 ′ -C5 ′ atoms and (b) C8 ′ -C7 ′ -C-C3 ″ atoms. Solid lines: uncharged molecules. Dotted lines: charged molecules.
Molecules 27 08417 g008
Molecules 27 08417 g009 550
Figure 9. (a) Snapshots of two single-aligned charged and uncharged DHICA-eumelanins: (b,c) 100 superimposed structures over the last 100 ns from the systems of four uncharged (b) and four charged (c) DHICA-eumelanins, respectively. The angle C6-C7-C4 ′ -C5 ′ is shown in blue, C8 ′ -C7 ′ -C-C3 ″ for the charged DHICA-eumelanin is in red, and for uncharged is in black.
Figure 9. (a) Snapshots of two single-aligned charged and uncharged DHICA-eumelanins: (b,c) 100 superimposed structures over the last 100 ns from the systems of four uncharged (b) and four charged (c) DHICA-eumelanins, respectively. The angle C6-C7-C4 ′ -C5 ′ is shown in blue, C8 ′ -C7 ′ -C-C3 ″ for the charged DHICA-eumelanin is in red, and for uncharged is in black.
Molecules 27 08417 g009
Molecules 27 08417 g010 550
Figure 10. Time evolution of the number of hydrogen bonds per monomer for two, four, 27, and fifty DHICA-eumelanin systems (a) between DHICA eumelanin and water and (b) between the eumelanin molecules.
Figure 10. Time evolution of the number of hydrogen bonds per monomer for two, four, 27, and fifty DHICA-eumelanin systems (a) between DHICA eumelanin and water and (b) between the eumelanin molecules.
Molecules 27 08417 g010
Table
Table 1. Details of the simulated systems.
Table 1. Details of the simulated systems.
Number of Eumelanins (nm 3 )Number of Water MoleculesBox Size (nm 3 )Counterions (K + )Duration ( μ s )
2 uncharged4218 ( 5.07 ) 3 01.073
4 uncharged5640 ( 5.59 ) 3 00.983
27 uncharged6341 7.16 × 7.30 × 7.14 02.237
50 uncharged20,000 ( 8.60 ) 3 02.050
2 charges3487 ( 4.75 ) 3 81.000
4 charges4501 ( 5.18 ) 3 161.100
27 charges6272 ( 5.87 ) 3 1083.802
50 charges19,800 ( 8.57 ) 3 2001.361
25 uncharged-25 charged19,900 11.74 × 7.39 × 7.32 1001.974
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Soltani, S.; Sowlati-Hashjin, S.; Tetsassi Feugmo, C.G.; Karttunen, M. Structural Investigation of DHICA Eumelanin Using Density Functional Theory and Classical Molecular Dynamics Simulations. Molecules 2022, 27, 8417. https://doi.org/10.3390/molecules27238417

AMA Style

Soltani S, Sowlati-Hashjin S, Tetsassi Feugmo CG, Karttunen M. Structural Investigation of DHICA Eumelanin Using Density Functional Theory and Classical Molecular Dynamics Simulations. Molecules. 2022; 27(23):8417. https://doi.org/10.3390/molecules27238417

Chicago/Turabian Style

Soltani, Sepideh, Shahin Sowlati-Hashjin, Conrard Giresse Tetsassi Feugmo, and Mikko Karttunen. 2022. "Structural Investigation of DHICA Eumelanin Using Density Functional Theory and Classical Molecular Dynamics Simulations" Molecules 27, no. 23: 8417. https://doi.org/10.3390/molecules27238417

APA Style

Soltani, S., Sowlati-Hashjin, S., Tetsassi Feugmo, C. G., & Karttunen, M. (2022). Structural Investigation of DHICA Eumelanin Using Density Functional Theory and Classical Molecular Dynamics Simulations. Molecules, 27(23), 8417. https://doi.org/10.3390/molecules27238417

Article Metrics

Back to TopTop