The Role of Metal Ions in the Electron Transport through Azurin-Based Junctions

Abstract

We studied the coherent electron transport through metal–protein–metal junctions based on a blue copper azurin, in which the copper ion was replaced by three different metal ions (Co, Ni and Zn). Our results show that neither the protein structure nor the transmission at the Fermi level change significantly upon metal replacement. The discrepancy with previous experimental observations suggests that the transport mechanism taking place in these types of junctions is probably not fully coherent.

1. Introduction

The study of metal–protein–metal junctions has opened up new avenues in the field of molecular electronics as it paves the way to develop new hybrid devices capable of exploiting proteins’ remarkable properties [1]. In particular, they have the inherent capability of transferring charge over long distances. Thus, one can foresee the production of sensors, flexible implants, solar cells and many other types of electrical devices incorporating proteins [2,3,4,5,6]. The significant amount of work accumulated in recent years on this topic has given rise to a new field, namely that of protein-based electronics, which has been given the name of proteotronics [7]. Research in this field has been carried out on various fronts. On the one hand, remarkable effort has been made to understand the exact nature of the charge transfer mechanism in protein-based junctions as well as to identify the components which play the dominant role in the transport process [8,9,10,11]. On the other hand, the possibility of modifying the conductance properties by the insertion of mutations [3,10,11] has also attracted a lot of interest. Research on the transport properties of the building blocks of proteins, namely amino acids and peptides, has also been carried out [12,13,14].

Recent experimental results highlight the possibility of tuning the current density through azurin-based junctions via the replacement of the central copper ion with other metals [15]. These measurements were performed on proteins monolayers lying on a conducting substrate and contacted to a macroscopic electrode at the opposite side. At low temperatures, the difference in current between the Cu(II)-azurin and the Zn-azurin was found to be nearly two orders of magnitude. This was ascribed to the fact that each metal contributes to levels lying at different energies with respect to the Fermi level. Combining organic and metal components in the central area of molecular junctions is indeed a strategy which has been used quite often in molecular electronics [16,17,18].

In the case of azurins, this possibility becomes particularly intriguing since the aforementioned results of Amdursky et al. [15] indicate that by modifying only one (metal) atom one can modify the conductance properties of a system of almost 2000 atoms. To shed light on this issue, we thus performed theoretical calculations based on density functional theory (DFT) and calculated the electron-transport properties of several metal–azurin–metal junctions based on different metallic ions. Numerous computational studies on the active site on metal-substituted azurins are reported in the literature [19,20], but none involved the study of a whole junction comprising an entire protein held between two metal electrodes. Azurins are one of the proteins which have been studied the most in the field of protein-based electronics; additionally, they have long been the object of many computational studies because of their key role in biological functions [21]. Over the years, various experimental techniques have been adopted to measure the conductance of proteins [6,8,10,22,23,24,25,26,27,28,29,30,31,32,33]. In this work, we focus on structures which are likely to be formed during the experiments based on STM (Scanning Tunneling Microscopy), which are generally designed to measure the conductance of a single protein. The exact nature of the electron transport through the blue-copper azurin has been the object of a long debate [10,34], the degree of coherence being one of the issues at the heart of the debate. More specifically, what still needs clarifying is whether the process is fully coherent or whether it is better described by a sequential-tunneling process consisting of various steps which are individually coherent. In the case of the metal-substituted proteins studied in [15], the low-temperature results showed evidence for coherent transport; it was suggested that both tunneling and hopping are possible, however, depending on the temperature and the type of junction [35]. For the present study, we are particularly interested in the possibility of fully-coherent transport, which was suggested by the temperature-independence of the conductance measured in several experiments based on azurins [8,11,23,26,32,36,37,38,39]. In [34], various types of configurations were studied, in which an STM tip was contacted to a blue-copper azurin in three different ways: in blinking mode (in which the tip is kept at a fixed distance from the surface and the protein is allowed to attach to the tip spontaneously), via top indentation (in which the tip is brought closer to the protein) and by lateral indentation (in which the tip approaches the protein sideways). It was found that, under the assumption of fully-coherent electron transport, measurable conductance values can only be obtained by slightly compressing the protein upon tip indentation from the top or by bringing the tip close to the protein from the side. In this study, we chose to focus on this latter scenario. We compared the conductance properties of the blue-copper azurin from Pseudomonas aeruginosa with that of junctions in which the Cu ion was replaced by Zn, Ni or Co. Such a substitution is well known to lead to the formation of stable compounds [40,41,42,43,44,45], which in fact were used for the aforementioned experiments [15]. The Apo (protein deprived of the metal ion) was also considered.

2. Methods

2.1. Azurin Molecular Dynamics Simulations Free Solved in Water

To inspect structural modifications resulting from the ion-substitution, we performed all-atom molecular dynamics (MD) simulations. To this aim, three different azurin variants were considered: Cu-azurin or native azurin, Co-azurin and Ni-azurin. In all cases, the starting configuration is the one available in the the protein data bank [46,47] (PDBID = 4 AZU), and for each variant the Cu atom was replaced by the corresponding ion. Subsequently, protons were added accordingly to the calculated ionization states [48] of its titratable groups at a pH of 4.5—corresponding to the pH used in reference experiments [10]. In all cases, the resulting structures had a zero net charge. The amino acids inter-atomic interactions were modeled using the standard ff14SB force field [49]. As for the ion-coordination site, composed by the ion and five surrounding amino acids (GLY45, HIS46, CYS112, HIS117 and MET121), we adopted the ab-initio derived force-field parameters reported in [20]. This force-field includes all bonded and non-bonded terms in the coordination sphere for the three metal ions here considered (these parameters are detailed in the Supplementary Materials). Being all simulations performed in an explicit water medium, these solvent atoms were described using the TIP3P force-field [50], which is consistent with the choice of both aforementioned force-field parameters.

All MD simulations were performed using AMBER18 software suite [51] with NVIDIA GPU acceleration [52,53,54]. Moreover, we used periodic boundary conditions (PBC) with a rectangular box and particle mesh ewald (PME), with a real-space cutoff of 10 Å, to account for long-range electrostatic interactions. Van der Waals (vdW) contacts were truncated at the real space cutoff of 10 Å and coordinates were saved every 1000 steps. The temperature of the system was adjusted by means of a Langevin thermostat [55] with a friction coefficient of $\gamma =1$ ps${}^{-1}$. For the simulations performed in the NPT ensemble, a Berendsen barostat [56] with a relaxation time of $t_p=1$ ps was used. The SHAKE [57] algorithm was used to constrain bonds containing hydrogen atoms, thus allowing us to use an integration time-step of 2 fs. The simulation protocol is analogous to the one used in [3], which consists of the following four stages: (1) preparing the system by embedding the protein in water in such a way that the minimum solute–water distance is 1 Å, thus resulting in a system with dimensions ∼70 Å $\times 63$ Å $\times 69$ Å; (2) energy minimization of the structures using a combination of steepest descent and conjugate gradient methods; (3) heating up the system from 0 to 300 K with a 2 ns NPT simulation thus ensure an even temperature and density/presure ($T=300$ K and $P=1$ atm) on the system; and (4) performing 500 ns-long MD production runs in the NVT ensemble.

2.2. Density Functional Theory Calculations on Metal–Protein–Metal Junctions

We carried out density functional theory (DFT) calculations on metal–protein–metal junctions obtained by replacing, in the geometries obtained for the work in [34], the Cu ion with Zn, Co or Ni without re-optimizing the structure. This strategy makes it possible to identify effects in the electrical conductance due to changes in the electronic structure only, ruling out factors which may rather be related with geometrical differences. As mentioned in the Introduction, we focused on the structures that describe a gold tip approaching the protein sideways since, among the possible scenarios analyzed in [34], this was found to yield one of the highest conductance values. Further details about the protocol followed for obtaining these metal–protein–metal geometries can be found in the Supplementary Materials.

In all cases, the substrate and tip employed in these calculations consist of 969 and 252 atoms, respectively. Note that the corresponding number of gold atoms employed in the MD simulations was instead much larger (2538 and 835, respectively): this is because, when importing the MD output structures for the Cu-based junctions as an input for the DFT calculations, the outermost layers in both the surface and the tip as well as the surrounding water molecules were removed in order to comply with the computational limitations. The OpenMX code was used, which is based on highly optimized numerical pseudoatomic orbitals (PAOs) [58,59] and which we already used for our previous work on proteins [8,34,60,61]. We used the Perdew–Burke–Ernzerhof (PBE) exchange and correlation functional [62] and norm-conserving pseudopotentials [63]. Single-$\zeta$ basis sets were used for all species involved in the calculations (including the gold tip and substrate) except for the metallic ion for which a double-$\zeta$ basis set was employed instead. The cutoff radii used for the PAOs were 8 Bohr for the metal ion, 7 Bohr for S and 5 Bohr for C, N, O and H. As for the gold electrodes, the Au atoms belonging to the tip and the substrate were described by the cutoff radii of 9 Bohr in those regions closer to the protein, whereas a 7-Bohr cutoff radius was chosen for the rest of the atoms. Due to computational reasons, these calculations were performed without spin polarization. The electronic self-consistency is achieved by a Pulay mixing scheme based on the residual minimization method in the direct inversion iterative subspace (RMM-DIIS) [64] with a Kerker metric [65], using an energy cutoff of ${10}^{-8}$ Ha as convergence criterion. Further details about the strategies used to reach the self-consistent-cycle convergence can be found in [60,61]. The overlap and Hamiltonian matrices obtained by the DFT calculations were subsequently used to compute the zero-bias transmission as a function of energy in the spirit of the Landauer formalism and making use of Green’s function techniques. Within this method, the structure is divided in three regions (left (L), center (C) and right (R)) and the energy dependent transmission is given by

$$\tau \left(E\right)=\mathrm{Tr}\left[{\Gamma }_{\mathrm{L}}{\mathbf{G}}_{\mathrm{CC}}^{\mathrm{r}}{\Gamma }_{\mathrm{R}}{\mathbf{G}}_{\mathrm{CC}}^{\mathrm{a}}\right],$$

(1)

where ${\Gamma }_{\mathrm{L}}/\mathrm{R}$ are the scattering rate matrices and ${\mathbf{G}}_{\mathrm{CC}}^{\mathrm{r}\left(\mathrm{a}\right)}$ is the retarded (advanced) Green’s function (see [34] for further details).

3. Results and Discussion

First, we inspect the influence of ion replacement on the protein structure. Henceforth, we refer to X-azurin (X = Ni, Cu, Co and Zn) as the azurin containing the X-ion. The four metals considered all belong to the first row of the d block of the periodic table. In Figure 1a–c, we represent the time averaged configuration obtained from a 0.5 μs long MD simulation of the three different proteins, i.e., Cu-azurin (black), Co-azurin (red) and Ni-azurin (green). From the structure overlap between the different proteins and native azurin (i.e., Cu-azurin), it becomes apparent that the ion replacement has little influence on the protein conformation. This is further corroborated and quantified, by both small root-mean-square-deviation (RMSD) between average structures (shown in Figure 1a–c) as well as the preservation of the most prevalent secondary structure on the protein, i.e., $\beta$-sheets shown in Figure 1e. The $\beta$-sheets refer to a particular spacial organization of the protein’s amino acids regulated/characterized by a specific hydrogen bond network among them. In the particular case of azurin, the most prevalent sort of secondary structure is the $\beta$-sheet, which are further arranged in a circular manner, thus forming a characteristic $\beta$-barrel. In Figure 1e, we quantify the content of this secondary-structure for the different variants to realize that the ion replacement does not alter/disrupt the H-bond network sustaining the characteristic shape and structure of the wild-type azurin. Moreover, from the RMSD, we observed that all proteins exhibit a structural deviation from the native protein smaller than 1.5 Å. Such RMSD values are within the structural variations experienced by the protein resulting from thermal fluctuations. Moreover, the process of protein adsorption and subsequent tip indentation results in a structural deformation which is at least three times larger than the one observed from altering the protein ion (see Figure 2). All considered, one may state that, by altering the coordination ion, the structural modifications thus introduced are substantially smaller than the ones the protein experiences during the contact formation process. Note that the similarity of the parameters governing the bond motif of the different metal ions (e.g., charge and bond distances provided in Figures S2–S5) precluded the negligible structural modifications here observed. All considered, by the aforementioned reasons, the subsequent simulations were performed using a tip-protein contact equilibrated.

Figure 2 shows the initial and final output geometries of the MD simulation performed on the Cu-azurin from which the structures used for the present work were obtained (see more details in [34]). In this simulation, the tip is brought close to the protein sideways and the contact is established via the $\alpha$-helix arm and hidrophobic patch.

In the left column of Figure 3, we show the zero-bias transmission as a function of energy for the four proteins considered, at four different values of the tip-ion distance (d${}_{t-X}$). Such a distance was evaluated as the distance between the center of the lowest tip layer and the metal ion. In all curves, all resonances appear very sharp, indicating very low hybridization between the corresponding orbitals and the metallic leads. As expected, for most energy values the transmission decreases for increasing distance.

In the right column of the same figure, we report the projected density of states (PDOS) on the metal ion for the shortest tip-protein distance. The curve for Zn clearly differs from the other three in that the d-level peaks lie at lower energies with respect to the Fermi level as compared to the other metal ions (this, in turn, derives from the closed-shell electronic Zn configuration). This difference is reflected in the corresponding transmission curve, in which no peaks appear in a broad range immediately below the Fermi level as opposed to the other three proteins. Indeed, the lower density currents measured for the Zn-protein junctions and reported by Amdursky et al. [15] were ascribed to the different energy alignment of the metal-derived levels. A common feature to the transmission curves of all proteins is the energetic position of the lowest occupied orbital (localized on HIS35), which shows minimal changes from one curve to another. This is not surprising because, in our simulations, the position of the amino acids (including HIS35) within the protein is the same for all cases. Consequently, they contribute to the transmission with similar features. Note that, although the lack of re-optimization might in principle lead to some artifacts, the MD simulations discussed above revealed minimal changes upon ion replacement. This suggests that they would result into energy shifts of individual orbitals the effect of which would most likely be hidden by the rest of the contributions of the whole electronic structure of the protein (given the extreme sharpness of all peaks in the transmission curves of Figure 3, it is hard to believe that any shift of one of these peaks would affect the transmission at the Fermi level). The same applies to possible shifts of individual peaks arising from the spin polarization of the metallic ion which could not be described in our closed-shell calculations.

Figure 4 (left) shows a comparison for the transmission curves obtained for the four proteins at the shortest tip-ion distance in a narrow energy range around the Fermi level. The transmission curve for the Apo (the protein deprived of the metal ion) is also shown in the same graph. It mainly differs from the other curves in the obvious absence of the peaks which derive, in the case of the other four proteins, from the d states of the metal ions. Quite strikingly, the transmission at the Fermi level for the metal-substituted azurins and the Apo is almost identical. This is due to the extreme sharpness of the HOMO-related resonances deriving from the metal ions. Similar conclusions can be drawn for the other three tip-protein distance values (see Figure 4 (right) for a comparison of the transmission values at the Fermi level at all distances for all proteins). In the case of Co, Ni and Cu, these results suggest that, despite these metal ions contributing to levels close to the Fermi level, their role in a fully-coherent electron transport is certainly important but not so dominant as one might naively expect. A large part of the contribution seems instead to be given by the rest of the protein, which provides the same background transmission for the Apo and for the other proteins analyzed in this study. For the short-distance structure corresponding to the curves shown in Figure 4, the main role seems to be played by the LUMO (this is less visible in the case of the curves of Figure 3 corresponding to larger distances, which seem to rather indicate an equal contribution between HOMO and LUMO). In any case, one has to bear in mind that the occupied energy range immediately below the Fermi level includes the contribution of several residues (mostly glutammic and aspartic acids: for a detailed analysis in the case of the Cu-azurin, we refer the reader to the work of Romero-Muñiz et al. [34]). In the case of the Zn-azurin and for the Apo, part of the contribution to the transmission at the Fermi level seems indeed to be provided by these residues as well. As an aside, it is worth mentioning that recent experiments on protoporphyrins have revealed that the metal redox center has no role in the electron transfer, which was instead found to be mediated solely by the conjugated backbone of the molecule [68]. Our results are based on DFT, which is known to be affected by uncertainties concerning the size of the HOMO–LUMO gap because of self-interaction errors in the standard exchange-correlation functionals and image charge effects [69,70,71,72,73,74]. Nevertheless, it is hard to foresee that the main conclusions would change upon use of a more sophisticated exchange-correlation functional. Although in this study we focused on sideways approach of the tip, other mechanisms would be possible such as indentation from the top. Geometries obtained via this kind of simulation, however, would probably lead to similar conclusions as those drawn for the geometries discussed above. Such a prediction is based on the fact that, for the Cu-azurin for instance, the main features of the transmission curves (concerning coupling and energetic alignment) were found to be quite similar [34].

We now turn to comparing our results with the experimental ones reported in [15]. There, the current density was found to follow the trend Cu > Ni > Co > Zn. This is at odds with our results which, as explained above, did not show significant changes for the transmission at the Fermi level. However, our findings, based on single-protein junctions, cannot be directly compared with those experiments, which were carried out on self-assembled monolayers; the latter are affected by dipole–dipole interaction and other effects [75] which are not present in our systems. Moreover, the electrodes which were there used (a highly-doped Si substrate and, in certain cases, a top electrode made of Hg) are very different from the gold electrodes employed in our simulations. Amdursky et al. [15] also reported a decrease of two orders of magnitude for the current through the Apo with respect to the Cu-azurin for a large range of temperatures. However, we also note that, in [11], the Apo was actually found to show similar conductance to that of Cu-azurin, although this was attributed to different tunneling distances. In any case, the discrepancies between our results and those reported in [15] suggest that probably the electron-transport mechanism taking place in these systems is not as completely coherent as was claimed in multiple works [8,11,23,26,32,36,37,38,39]. It also highlights the fact that probably a proper description of the transport process needs other effects to be taken into account such as the time dependence of the fluctuations of the surrounding ligands. Similar conclusions were actually drawn by Romero-Muñiz et al. [34] concerning the Cu-azurin on the basis of the very low transmission values obtained for a broad range of geometries. We believe that these latest results of ours further confirm the need of further investigation regarding the nature the charge transfer in these kinds of junctions.

4. Conclusions

In summary, we studied metal–protein–metal junctions based on modified blue-copper azurins in which the central copper ion was replaced by Zn, Ni or Co, while leaving the rest of the structure unmodified. Comparison with the Apo protein was also performed. The striking similarity of the transmission values at the Fermi level among the five proteins suggest that, within a fully-coherent-transport picture, the dominant role would not actually be played by the metal ions, opposite to what common wisdom would suggest. The discrepancy with recent experimental observations suggest that a different kind of electron-transport mechanism is probably involved, although simulations employing more similar structures to those formed in the experiments of Amdursky et al. [15] should be carried out in order to draw definitive conclusions. Molecular dynamics simulations performed on the Cu-, Ni- and Co-azurin indicate that the ion substitutions do not alter the secondary structure significantly. The structural variation induced is much smaller than that arising from contact with the top electrode. Understanding whether such a difference can result in major changes in the conductance will be pursued in future studies.