Theoretical Characterization of Structural Disorder in the Tetramer Model Structure of Eumelanin

TL;DR

This paper develops a piecewise Ising Model to predict formation enthalpies of eumelanin tetramers, offering a new theoretical approach to understanding its complex disorder and properties.

Contribution

The work introduces a modular Ising Model for eumelanin, enabling detailed energetic landscape analysis and extending to larger polymer structures.

Findings

Superior predictive accuracy of the Ising Model for tetramer formation enthalpies

Potential to analyze larger eumelanin structures like hexamers and octomers

Framework for integrating absorbance data with polymer abundance predictions

Abstract

Eumelanin is regarded to be an attractive candidate material for biomedical applications. Despite many theoretical studies exploring the structure of eumelanin, an exact mapping of the energetic landscape of the very large phase space of eumelanin is still elusive. In this work, we implement a piecewise Ising Model to predict formation enthalpies of Eumelanin single and double tetramers, and demonstrate its superior predictive and generalizable capabilities. We believe this model will prove very useful in theoretically characterizing the many unique properties attributed to its disorder. The modular nature of the predictive Ising model built up in this work is well-suited for analysis and characterization of a larger phase space of eumelanin polymers, including hexamers and octomers, as well as larger stacked structures, such as potential triple and quadruple eumelanin tetramers. Absorbance data can be incorporated with population-wide predictions of polymer abundance to produce weighted-average predictions of broadband absorbance of bulk eumelanin.

Figures

Tables

Table 1. Table 1: Statistics on calculated Ising Model coefficient values, in eV, summarized across each general type of coefficient. Mean of means and standard deviation of means data provided for coefficients calculated using the 1-Step and the 3-Step calculation schemes.

Coeff	Grand $\mu$ (1S)	Grand $\sigma$ (1S)	Grand $\mu$ (3S)	Grand $\sigma$ (3S)
h	-0.28	0.08	6.49	1.11
jin	-0.26	0.23	-4.03	0.24
jid	-0.46	0.30	-7.40	0.22
jou	-0.46	0.37	-0.14	0.17
jon	-0.20	0.23	-0.05	0.11
jod	-0.44	0.04	-0.12	0.06

Table 2. Table 2: Statistics on calculated coefficient values, in eV, in the eumelanin Ising Model. Mean and standard deviation data provided for individual coefficients across the 2000 values exported with the BEEF-vdW ensemble, for Ising Models computed using the 3-Step calculation scheme.

Coeff.	$\mu$ (3S)	$\sigma$ (3S)
hH	5.23	1.21
hI	6.92	1.25
hM	7.33	1.29
jinHH	-3.82	0.56
jinHI	-3.89	0.56
jinHM	-4.32	0.57
jinII	-3.88	0.56
jinIM	-4.37	0.57
jinMM	-3.92	0.56
jidHH	-7.16	1.03
jidHI	-7.23	1.03
jidHM	-7.55	1.04
jidII	-7.24	1.03
jidIM	-7.68	1.04
jidMM	-7.56	1.03
jouHH	-0.18	0.1
jouHI	-0.43	0.14
jouHM	-0.17	0.14
jouII	-0.07	0.07
jouIM	-0.01	0.11
jouMM	0.04	0.09
jonHH	-0.11	0.06
jonHI	-0.17	0.05
jonHM	-0.06	0.06
jonII	0.11	0.09
jonIM	-0.1	0.05
jonMM	0.04	0.07
jodHH	-0.09	0.07
jodHI	-0.15	0.07
jodHM	-0.06	0.07
jodII	-0.08	0.09
jodIM	-0.11	0.07
jodMM	-0.23	0.05

Equations

$$\Delta H_{T-T}^a=U_{T-T}^{DFT}-\sum _C{E}_C^{Ref}-\sum _H{E}_H^{Ref}-\sum _N{E}_N^{Ref}-\sum _O{E}_O^{Ref}+\Delta pV$$

$$E_O^{Ref}=U_{H_{2}O}^{DFT}-U_{H_2}^{DFT}-\Delta H_{H_{2}O}^{Exp}$$

$$E_N^{Ref}=U_{N{H}_3}^{DFT}-3/2U_{H_2}^{DFT}-\Delta H_{N{H}_3}^{Exp}$$

$$\Delta H_T^a=U_T^{DFT}-\sum _C{E}_C^{Ref}-\sum _H{E}_H^{Ref}-\sum _N{E}_N^{Ref}-\sum _O{E}_O^{Ref}$$

$$\Delta H_M^a=U_M^{DFT}-\sum _C{E}_C^{Ref}-\sum _H{E}_H^{Ref}-\sum _N{E}_N^{Ref}-\sum _O{E}_O^{Ref}$$

$$\begin{gathered} \Delta H_{T-T}^a=\sum _{⟨i⟩}{h}_i{\sigma }_i+\sum _{⟨ik⟩}{j}_{i{n}_{i,k}}{\sigma }_i{\sigma }_k+\sum _{⟨ik⟩}{j}_{i{d}_{i,k}}{\sigma }_i{\sigma }_k+ \\ +\sum _{⟨ik⟩}{j}_{o{u}_{i,k}}{\sigma }_i{\sigma }_k+\sum _{⟨ik⟩}{j}_{o{n}_{i,k}}{\sigma }_i{\sigma }_k+\sum _{⟨ik⟩}{j}_{o{d}_{i,k}}{\sigma }_i{\sigma }_k \end{gathered}$$

$$\Delta H_{M_i}^a-2E_H^{Ref}=h_i{\sigma }_i$$

$$\Delta H_T^a=\sum _{⟨i⟩}{h}_i{\sigma }_i+\sum _{⟨ik⟩}{j}_{i{n}_{i,k}}{\sigma }_i{\sigma }_k+\sum _{⟨ik⟩}{j}_{i{d}_{i,k}}{\sigma }_i{\sigma }_k$$

$$\Delta H_T^a-\sum _{⟨i⟩}{h}_i{\sigma }_i=\Delta H_T^M$$

$$\Delta H_T^M=\sum _{⟨ik⟩}{j}_{i{n}_{i,k}}{\sigma }_i{\sigma }_k+\sum _{⟨ik⟩}{j}_{i{d}_{i,k}}{\sigma }_i{\sigma }_k$$

$$\Delta H_{T-T}^a-(\sum _{⟨i⟩}{h}_i{\sigma }_i+\sum _{⟨ik⟩}{j}_{i{n}_{i,k}}{\sigma }_i{\sigma }_k+\sum _{⟨ik⟩}{j}_{i{d}_{i,k}}{\sigma }_i{\sigma }_k)=\Delta H_{T-T}^T$$

$$\Delta H_{T-T}^T=\sum _{⟨ik⟩}{j}_{o{u}_{i,k}}{\sigma }_i{\sigma }_k+\sum _{⟨ik⟩}{j}_{o{n}_{i,k}}{\sigma }_i{\sigma }_k+\sum _{⟨ik⟩}{j}_{o{d}_{i,k}}{\sigma }_i{\sigma }_k$$

Full text

I Introduction

Eumelanin is a subgroup of melanin pigments found in living organisms that plays an important role in skin coloration and UV protection.Friedmann and Gilchrest (1987) The precise chemical structure of eumelanin is still not completely known because it is highly cross-linked and insoluble in available solvents.Ito et al. (2013) Among organic polymers, eumelanin occupies a unique position because of: (i) its widespread occurrence in nature, from people and mammals to fish, birds and molluscs;Edwards et al. (2014) (ii) the variety of biological roles, from photoprotection to scavenging of reactive oxygen speciesNofsinger et al. (2002); Bustamante et al. (1993) and metal chelation;Im et al. (2017) and (iii) distinct physical and chemical properties, including broadband photoabsorption throughout the visible range,Tran et al. (2006) water-dependent ionic-electronic semiconductor-like behavior,Meredith and Sarna (2006) stable free radical character and efficient nonradiative energy dissipation,Ju et al. (2011) making eumelanin an attractive candidate for biomedical and technological applications. Despite growing interest in eumelanin-type functional materials and systems, the exact structural underpinnings due to the highly insoluble and heterogeneous character of these polymers, has proved challenging.d’Ischia et al. (2009)

In parallel, interest in miniaturized medical implants and edible biometric sensors has led to a need to develop novel biodegradeable batteries, based on eumelanin extracted from the common cuttlefish ($Sepia~{}officinalis$)Kim et al. (2016, 2013). Electrochemical characterization indicates the possibility of electrochemical intercalation of up to two sodium ions per eumelanin unit.Kim et al. (2016) The theoretical analysis in the work of Kim et al.Kim et al. (2016) utilizes a stacked tetramer model to rationalize these findings, originally proposed by Kaxiras et al.Kaxiras et al. (2006) Recently, another study has explored the geometric complexity possible within the Kaxiras model.Di Mauro et al. (2017); Vahidzadeh et al. (2018); Crescenzi et al. (2017) While these studies represent important strides, a complete energetic landscape of double tetramer is still elusive given the large phase space.

Systematic exploration of large phase spaces for crystalline materials is enabled through the Ising Model (or cluster expansion).Ceder (1993); De Fontaine et al. (1992); Tepesch et al. (1995) The general methodology involves calculating the system energy for a subset of structures and training a model that can be used to then subsequently predict the rest of the phase space with very high accuracy. In this work, we develop an Ising model to describe the interactions in-plane and out-of-plane for a double tetramer. Utilizing 108 density functional theory calculations, we train an Ising model that predicts on an test set with an accuracy of $0.15$ eV. The model is generalizable and allows for an accurate mapping of the energetic landscape of a double tetramer. We believe this analysis can form the basis for further characterization of broadband absorbance, electrochemical ion interacalation, etc.

II Methods

II.1 DFT Simulations and Structure

In this work, the primary molecular structure considered for this study is a double eumelanin tetramer, building on the work of Kaxiras et al.Kaxiras et al. (2006) While other possible structures such as hexamerRandhawa et al. (2009); Stark et al. (2005) and octomerSlominski et al. (2004) are important, we have chosen to utilize the double tetramer structural model to demonstrate our methodology. As will be discussed later, this model can be extended to other possible structures of eumelanin.

The specific structures simulated for this study includes 3 types of eumelanin monomers: hydroquinone (HQ), indolequinone (IQ), and quinone-methide (MQ).Meng and Kaxiras (2008) The fourth and final monomer type described by Kaxiras, quinone-imine (NQ), was not used in this study, since it is an isomer of quinone-methide (MQ). In the context of planar tetramer assembly and ion intercalation, its structural difference is negligible, and thus it is expected to produce structures with similar formation enthalpies. The monomer structures are represented in Fig. 1.

We introduce the following naming scheme to describe the geometric arrangement of simulated structures. The double eumelanin tetramers are planar in the X-Y plane, stacked one above the other along the Z axis. This is one of the stable stacked structures when vdW interactions are taken into account.Kim et al. (2016) Looking down the Z axis onto the X-Y plane, the monomers of the lower tetramer are named first, starting with the monomer in Quadrant 2 and proceeding clockwise. Following the 4 monomers of the lower tetramer, the 4 monomers of the upper tetramer are named, in the same order. All monomers are designated by the first letter of their acronym (H for hydroquinone, I for indolequinone, M for quinone-methide). Thus, for instance, in an example tetramer HHHM-HHHM, the 2 sets of 3 hydroquinones, in Quadrants 2, 1, and 4, are situated one above the other, and the 2 sets of single quinone-methides, in Quadrant 3, are likewise situated one above the other. Fig. 2 illustrates the structures of an example planar tetramer and an example double tetramer.

In addition to simulations of the double eumelanin tetramer, single eumelanin planar tetramers were also simulated, to compute formation enthalpies of the double tetramers from single tetramers. As will be discussed later, the purpose of these calculations was to ensure the out-of-plane interactions between the simulated pairs of planar double tetramers were properly quantified and accounted for.

Finally, single eumelanin monomers were also simulated, to compute formation enthalpies of both the single and double eumelanin tetramers. These single monomers were simulated in their stable standalone configuration, with an additional 2 hydrogen atoms, which are lost in the polymerization process into tetramers, retained on the monomers.

Self-consistent DFT calculations were performed using the Projector Augmented Wave Method as implemented in GPAW,Enkovaara et al. (2010) with the Bayesian Error Estimation FunctionalMortensen et al. (2005) with van-der-Waals correlationWellendorff et al. (2012) (BEEF-vdW) exchange-correlation functional. All molecules were simulated in standalone bulk cells measuring 20 x 20 x 17 Å, with all periodic boundary conditions disabled. Calculations were run with a real-space grid of 0.18 Å The Conjugate Gradient eigensolver was used to facilitate convergence of the simulations. Fermi-Dirac occupation smearing of 0.01 eV was used to expedite simulation convergence while providing accurate results. All structures were geometrically relaxed until the net force in the simulated molecule was decreased to 0.05 eV.

II.2 Formation Enthalpies from Bulk Simulations

Several sets of formation enthalpies were calculated for the assemblage of simulated structures, including relevant monomers, planar tetramers, and double tetramers.

First, formation enthalpies of only the double tetramers (T-T) were calculated, using DFT-calculated internal energies of the double tetramers and DFT-calculated reference energies for the constituent elemental species: carbon, hydrogen, nitrogen, and oxygen.

[TABLE]

The pressure-volume work term, $\mathrm{\Delta pV}$, can be disregarded, as it is typically about 5 orders of magnitude smaller than internal energy contributions in formation enthalpy calculations.Aydinol et al. (1997); Obrovac and Chevrier (2014)

Elemental reference energy for hydrogen was calculated using simply the DFT-calculated internal energy of hydrogen gas, H2. Elemental reference energy for carbon was calculated using bulk graphite.Pande and Viswanathan (2018) The reference energies for oxygen and nitrogen, however, required additional correction, since both are well-known to be poorly described within DFT.Jacob (2006); Rod et al. (1999) Oxygen reference energy was calculated using the standard water-reference scheme,Baran et al. (2014) using DFT-computed internal energies of water and hydrogen gas, as well as the experimental formation enthalpy of water:

[TABLE]

Nitrogen reference energy was calculated using a similar, ammonia-reference scheme,Wang and Pantelides (2011) using DFT-computed internal energies of ammonia and hydrogen gas, as well as the experimental formation enthalpy of ammonia:

[TABLE]

In addition to calculating formation enthalpies of the double tetramers from the elemental basis, formation enthalpies of the single planar tetramers (T) and the single monomers (M) were also computed, using the same elemental basis.

[TABLE]

[TABLE]

II.3 Ising Model Coefficients

In order to map out the energetic interactions between the eumelanin monomers, both within the planar tetramers and across the planes of the double tetramers, we utilize a modified Ising Model for the lattice Hamiltonian. As implemented, the model consists of a lattice of $N$ sites $i$, whose filling is described by occupation terms, $\sigma_{i}$. All terms within the Ising Model are calculated using formation enthalpies of the monomers and polymers, referenced to their constituent atomic species, as described above. The occupation energies, $h$, in this implementation of the Ising Model, represent the formation enthalpies of the monomers directly. As there are 3 types of monomers used in the explored configuration space of eumelanin double tetramers, 3 types of occupation terms were used: $h_{H}$, $h_{I}$ and $h_{M}$, each term corresponding to the formation enthalpy of its corresponding monomer. In addition to the occupation terms $h$, energy contributions to the full polymer due to monomer-monomer interactions were captured by the interaction terms $j$. Several types of interaction terms are identified, separated into in-plane interactions between monomers situated in the same planar tetramer, $j_{i}$, and out-of-plane interactions between monomers situated in opposing tetramers, $j_{o}$. Within each category, interactions between monomers occupying nearby quadrants are accounted for, $j_{in}$ and $j_{on}$, as well as interactions between monomers occupying diametrically opposed quadrants, $j_{id}$ and $j_{od}$. Finally, interactions between monomers situated in opposing dimers but in the same geometric quadrant, one underneath another, are accounted for as $j_{ou}$. Thus, a total of 5 major types of interaction terms are used in the Ising Model: $j_{in}$, $j_{id}$, $j_{ou}$, $j_{on}$, $j_{od}$. Within each type of interaction term, there is further distinction as to the types of monomers involved, with the full set of permutations between H, I, and M accounted for: HH, HI, HM, II, IM, and MM.

Two schemes were used to calculate the above occupation and coupling terms for the Ising Model. In the first scheme, the formation enthalpy data of eumelanin double tetramers, calculated through DFT simulations, was used to derive the full set of corresponding Ising Model coefficients in a single step, though a least-squares regression fit using Wolfram Mathematica. Each double tetramer formation enthalpy was described as:

[TABLE]

The second scheme broke down the calculation of Ising Model coefficients into three distinct steps, using the formation enthalpies of the monomers, single tetramers, and double tetramers for each step, respectively. First, the formation enthalpies of the monomers were used to calculate the occupation terms $h$ directly. A correction had to be implemented for hydrogen, since each monomer lost two hydrogen atoms when it was polymerized into a tetramer:

[TABLE]

Next, the formation enthalpies of the single planar tetramers were used to calculate the in-plane coupling coefficients $j_{i}$. The formation enthalpy of a single tetramer contains contributions both from the presence of individual monomers and their in-plane interactions:

[TABLE]

To isolate the energetic contributions of in-plane interactions, we subtracted the occupation-energy contributions of the constituent monomers, which were quantified in the preceding step, from the formation enthalpy of the tetramer. In effect, the formation enthalpy of a tetramer from the monomer basis was calculated. The resultant energy was used to calculate the in-plane coupling coefficients, though a least-squares regression fit using Wolfram Mathematica:

[TABLE]

[TABLE]

Finally, the formation enthalpies of the double tetramers were used to calculate the out-of-plane coupling coefficients $j_{o}$. As discussed in the one-step scheme, the formation enthalpy of the double tetramer contains information about all types occupation and coupling terms, so we subtracted the occupation-energy contributions, calculated using the monomers, and the in-plane interaction-energy contributions, calculated using the single tetramers. Effectively, the remaining energy represented the formation enthalpy of a double tetramer from its constituent single tetramers. This formation enthalpy was used to fit the out-of-plane interaction coefficients, though a least-squares regression fit using Wolfram Mathematica:

[TABLE]

[TABLE]

II.4 Double Tetramer Phase Space

Following the derivation of the interaction coefficients, these coefficients were applied to predict the formation enthalpies of the entire phase space of relevant eumelanin double tetramers. If every possible permutation of HQ, IQ, and MQ within a double tetramer is considered, there exist a total of 6561 distinctly named double tetramers. However, many of these tetramers are simply rotated or flipped equivalents of other tetramers within the same phase space. To account for this phase space degeneracy, the full set of Ising Model coupling coefficients for each of the 6561 configurations was calculated analytically, and families of configurations which matched the set of coupling coefficients had all but one member removed from the phase space. This left a total of 1032 distinct, non-degenerate double tetramer configurations in the phase space.

To provide training data for the Ising Model, a total of 85 double tetramers composed of HQ, IQ, and MQ monomers were simulated in DFT using the BEEF-vdW exchange-correlation functional. In addition, the 3 monomers themselves were simulated in DFT, as well as 20 intermediate single planar tetramers.

This study yielded two sets of interaction coefficients, as described above: one derived from the formation enthalpies of the double tetramers calculated directly from the atomic basis in one step, the other derived from the formation enthalpies of the double tetramers, single tetramers, and monomers, all calculated from the atomic basis, in three steps. The use of the BEEF-vdW exchange correlation functional generates a non-self consistent ensemble of energies, which was used to train an ensemble of Ising models. It has been shown that the ensemble of energies reliably reproduces trends in energies at the GGA-level.Houchins and Viswanathan (2017); Vinogradova et al. (2018); Krishnamurthy et al. (2016, 2018)

To most adequately represent the distribution of predicted formation enthalpies of the melanin double tetramers, it was decided to use a histogram plot of the calculated formation enthalpies predicted for every double tetramer configuration. To produce these plots, the double tetramers were sorted in order of increasing average formation enthalpy, as calculated across the full ensemble of 2000 values predicted by the BEEF-vdW functional. For every double tetramer, the standard deviation of formation enthalpy was calculated across all 2000 predicted BEEF-vdW values, to indicate uncertainty in the calculation. This standard deviation was represented as a band on the histogram of formation enthalpies, centered on each individual column of the histogram. The band provided for each column represents the standard deviation of formation enthalpies relevant only to the double tetramer structures present in that column. Since the standard deviation slightly differed for individual double tetramers, the resultant bands are slightly but negligibly different for each column.

III Results $\&$ Discussion

Using the 1-step calculation for formation enthalpy of the double tetramers, it was found that the average formation enthalpy across all simulated double tetramers, calculated from the atomic reference basis, was $-14.7$ eV, with a standard deviation of $3.1$ eV. Double tetramers containing a higher proportion of HQ (hydroquinone) were found to be the most stable, down to $-19.3$ eV, while double tetramers containing higher proportions of MQ (quinone-methide) exhibited the least negative formation enthalpies, up to $-5.22$ eV. The standard deviation of calculated formation enthalpy for each individual double tetramer, throughout the phase space of the BEEF-vdW ensemble was calculated as well. The average of these standard deviations for the set of simulated structures was $2.9$ eV, with an overall standard deviation of just $0.1$ eV.

The average formation enthalpy of an individual planar tetramer, calculated from the atomic basis, was $-5.09$ eV, with a standard deviation of $1.92$ eV. Thus, it was expected the formation enthalpy of the double tetramer from single planar tetramers would be of comparable order-of-magnitude to the formation enthalpy of the said tetramers. As calculated, the average formation enthalpy of a double tetramer from its constituent single tetramers was $-1.88$ eV, with a standard deviation of $0.58$ eV. From this observation, it was apparent that the out-of-plane coupling coefficients derived through the Ising Model should be lower than the in-plane coupling coefficients.

Finally, formation enthalpies of the individual monomers were calculated from the atomic basis. The average magnitude was found to be $-1.46$ eV, with a standard deviation of $1.11$ eV. Since each corresponding $h$ term, derived from the formation enthalpies of the monomers, had to account for the loss of 2 hydrogen atoms, the average magnitude of the occupation energies used in the 3-step Ising Model was $6.49$ eV.

In both least-squares fits of the coupling coefficients used with the Ising Model, as fit to the training data, the reported $r^{2}$ was very high, in excess of $0.99$, for all ensemble-specific sets of coefficients. However, the values derived for the coupling coefficients differed strongly. Tab. 1 outlines the results summarized across each type of coupling coefficient, while Table 2 provides a more complete database of coefficient values, with statistics available for each individual coefficient present in the system.

It is evident that although both methods are able to match the formation enthalpies of the double tetramers well, the 1-Step calculation method provides a less physically interpretable model. In addition, this set of coupling coefficients fails to predict the formation enthalpies of single tetramers well. The predicted coefficient values using the 1-Step calculation are not generalizable as no data regarding the monomers and single-tetramers data were used in the training. It is worth highlighting that this 1-Step analysis would constitute a black-box application of the Ising model for double tetramer model.

On the other hand, the results provided by the 3-Step calculation method are far more systematic, generalizable and physically meaningful. In the full set of fit Ising Model coefficients, it was found that the occupation coefficients and the in-plane coupling coefficients contribute the bulk of the formation enthalpy of a double tetramer, as was expected since out-of-plane interactions are dominantly vdW interactions and do not involve any chemical bonding. Of the in-plane interaction coupling coefficients, the coefficient across the diagonal of a single tetramer was found to be almost twice the magnitude of the nearest-neighbor coupling coefficient. We interpret this as indicative of the stability of the tetramer as a standalone polymer: if the nearest-neighbor coupling coefficient was much greater in magnitude than the diagonal coupling coefficient, the monomers may have preferred to form a ribbon that could stretch indefinitely, rather than forming a self-contained tetramer unit. Finally, all out-of-plane coupling coefficients were found to be much lower than the in-plane coefficients, as expected from the formation enthalpy calculations.

The effect of this disparity between the two methods can be seen in the predicted statistics regarding the expected formation enthalpies of the full double tetramer phase space. Fig. 3 shows the histogram of formation enthalpies of all double tetramers, as calculated using the 3-step computational method described above. In the plot, it is easy to see that there are some independent peaks, which correspond to double tetramer phases with a specific number of hydroquinone monomers, but the histogram overall follows a standard Bell curve. For comparison, Fig. S1 (provided with the Supporting Information document) shows the same histogram produced for the formation enthalpies calculated using the 1-step computational method. That histogram demonstrates much more significant peak formation, indicating a much stronger dependence of formation enthalpy on the presence of hydroquinone monomers, and, by extension, the out-of-plane coupling interactions, as these are magnified in the 1-step computation.

From these two formation enthalpy plots, it can be seen that the 1-Step method and the 3-Step methods generally predict the ensemble of formation enthalpies within the same range of values, both throughout the phase space of eumelanin double tetramer configuration, as well as the phase space of the BEEF ensemble. However, the characteristics of this phase space as predicted by the 1-Step and 3-Step methods differ. The 1-step method over-predicts an excessive reliance on the specific geometric configuration of the double tetramer, most likely due to its excessive emphasis on out-of-plane interactions within the phase space of tetramer geometry. The 3-step method reduces the enthalpic contribution from the out-of-plane interactions, and thereby produces a much smoother, more continuous plot. It is observed that the 3-Step method is more accurate at matching the formation enthalpies of calculated test data.

Formation enthalpy calculations showed that the single and double tetramers with the highest fraction of HQ monomers possess the lowest, most stable formation enthalpies, as compared to eumelanin structures with a higher fraction of IQ and MQ monomers. This observation was replicated throughout the entire phase space of double tetramer structures predicted by the Ising Model. This observation implies that within a mixture of eumelanin monomers polymerizing into double planar eumelanin tetramers, in the thermodynamic limit, double tetramers rich in HQ will be preferentially formed out of the available HQ population first, with only a small fraction of HQ monomers ultimately bound within IQ- and MQ-rich double tetramers.

The coefficients calculated from the Ising Model indicate that in-plane interactions are strongest between an MQ monomer and a dissimilar monomer; both the nearest-neighbor and the diagonal next-nearest-neighbor interactions between HQ and MQ or IQ and MQ are stabilized by 0.3 - 0.5 eV compared to analogous in-plane interactions not involving MQ. This indicates that single planar tetramers formed with a small amount of MQ monomers are more stable than monomers with a mixture of HQ and IQ monomers alone. Out-of-plane coefficients calculated from the Ising Model indicate that the strongest interactions involve HQ monomers, either in combination with fellow HQ monomers or other monomers. This is expected, since the HQ-rich double tetramers were found to be the most stable structures investigated, stabilized both by the presence of HQ monomers and their out-of-plane interactions.

These observation has significant implications for the anticipated intercalation potentials of ions into these double tetramers, as one of the preferred locations for intercalation of ions is the double tetramer’s inner ring. The HHHH-HHHH-type double tetramer carries 8 hydrogen atoms protruding into the inner ring, which may hinder ion intercalation, while the IQ-rich and MQ-rich double tetramers carry fewer hydrogen atoms protruding into the inner ring, and may possess strong binding centers for ion intercalation.

Recent experimental and $ab-initio$ simulation work indicates that the eumelanin tetramer and the eumelanin stacked double tetramer play an important role in metal ion binding and broadband absorption of bulk eumelanin.Kim et al. (2016); Kaxiras et al. (2006) Our results arre consistent with earlier work by Kaxiras et al.Kaxiras et al. (2006) where they show that hydroquinone-rich phases of eumelanin display are consistently more stable of the eumelanin polymers, ranging from dimers to tetramers, as compared to polymers poor in HQ. The modular nature of the predictive Ising model built up in this work is well-suited for analysis and characterization of a larger phase space of eumelanin polymers, including hexamers and octomers, as well as larger stacked structures, such as potential triple and quadruple eumelanin tetramers. Furthermore, as the model allows for phase-space-wide prediction of relative stability of specific polymers, it can be incorporated in a predictive model used to study the broadband absorbance of bulk eumelanin, by utilizing a weighted average of bulk polymer composition and calculated absorbance data.

IV Conclusions

In this study, we have built a generalizable, physically meaningful Ising model to describe the energetic interactions of eumelanin. This model allows a precise mapping of the very large phase of eumelanin structures within the double tetramer model. We highlight the importance of carrying out a step-wise training of the model coefficients using first the monomer enthalpies, followed by the single tetramer enthalpies and finally the double tetramer. The developed model will prove to be extremely useful in rapidly exploring the absorbance, ion interacalation of eumelanin within the double tetramer model. We also believe the developed methodology can be extended to other possible structural models for eumelanin such as the hexamer, octamer, etc.

V Acknowledgements

The computational portion of this work was performed on the Arjuna computer cluster, which was funded through Carnegie Mellon College of Engineering and the departments of Mechanical Engineering, Electrical and Computer Engineering and Chemical Engineering. The computational portion of this work was performed on Hercules computer cluster, which was funded through a Carnegie Mellon College of Engineering Equipment grant. O.B.S. acknowledges support from the DoD SMART Fellowship. The authors thank Professor Christopher Bettinger of Carnegie Mellon University for helpful discussions on Eumelanin structure and function.