Magnetic interactions in the catalyst used by nature to split water: a DFT + U multiscale study on the Mn₄CaO₅ core in photosystem II

Photocatalytic energy production by direct conversion of solar energy into a chemical fuel, such as hydrogen or hydrocarbons, is a major challenge in the research of new materials for energy. In particular, the discovery of new catalytic materials based on earth-abundant elements that are capable to perform the oxidation of water with high efficiency is a key ingredient for the future development of devices that are able to capture and store solar energy on a large scale [1]. Natural photosynthetic processes efficiently perform this task, capturing solar photons and converting them into chemical energy that can be stored in stable organic molecules. One of the most important steps in this process is represented by the oxidation of water, which is carried out by the oxygen evolving complex (OEC) of photosystem II (PSII), which is a protein/pigments complex that has remained almost evolutionary unchanged since its first appearance about 2.5 billion years ago [2]. The active site of the enzyme is composed of a core of four Mn ions and one Ca ion that are connected together through μ-oxo-bridges in an open cubane fashion [3]. In the last years, an inorganic cobalt-based catalyst film, having favorable working and stability conditions, was shown to be able to catalyze electrochemical water oxidation with small overpotential [4, 5]. The film is an amorphous solid that is self assembled from low cost cobalt salts at the anode. It has been recently combined with a silicon junction and with a hydrogen evolution catalyst to form a prototype 'artificial leaf' [6]. Interestingly, x-ray scattering measurements have revealed that the structure of this inorganic catalyst shares some common features with its natural counterpart [7–9].

Recently, Dau and co-workers [10, 11] have reported an Mn oxide that has also been characterized as water-oxidizing. The design of new materials based on Mn would be desirable since this element, which is that chosen by nature, is earth-abundant, inexpensive and largely non-toxic. Electronic structure calculations and ab initio molecular dynamics (MD) have been successfully used to unravel the structural and functional properties of molecular models of materials for water oxidation based on transition metals [12–14]. These studies are important to understand the similarities and differences of the electronic properties of the catalysts and they might also help in the design of new biomimetic materials, when compared to parallel studies on natural enzymes.

A high-resolution crystallographic structure of the PSII complex has recently revealed the atomistic details of the catalytic Mn₄CaO₅ core [3]. This structure provides the grounds for parallel computational investigations of the catalytic mechanism of the OEC [15–18]. It also helps the study of natural and synthetic analogous for water splitting reactions [9]. During the catalytic pathways toward oxygen evolution, the natural system cycles through different states, which are the S₀–S₄ states of the so-called Kok's Cycle [19]. Despite a large amount of data being available for the functional and spectroscopic properties of such states in different organisms and experimental conditions [20–23], the exact reaction mechanism and the strategy adopted by the enzyme to carry on its task are still elusive. An important aspect of the OEC, as well as of their artificial analogues, is represented by the magnetic interactions between the transition metals connected by oxo-bridges. In this respect, electronic paramagnetic resonance is successfully and routinely used to probe the oxidation state and the chemical environment of the PSII core along the Kok's cycle and its intermediates [21]. From the computational point of view, the description of such magnetic states can be conveniently described through broken-symmetry (BS) approach based on density functional theory (DFT) [24–29]. In a recent work, Pantazis et al [18] have proposed two distinct gas-phase (GP) molecular models of the Mn₄CaO₅ cluster in PSII, differing by the magnetic interactions between Mn ions. These two structures have two different overall spin ground states and their interconversion can explain the presence of two distinct signals observed by EPR spectroscopy in the S₂ state [20, 21, 30, 31]. Two important factors should be considered to fully describe the interplay of such magnetic interactions in biologically relevant conditions, which are: the effect of temperature and the role of the protein field. To achieve this goal, ab initio MD simulations coupled to quantum mechanics/molecular mechanics (QM/MM) scheme can be used to unravel both the effect of dynamics and the effect of the environment on the average values and fluctuations of the magnetic properties [32]. To properly describe the electron correlation of the transition metal centers, it is necessary to go beyond the approximations used by the standard DFT, either with the use of hybrid functionals or with the DFT + U scheme [33–37]. Since the latter approach is particularly computationally efficient for ab initio MD calculations, it has been adopted for the present study.

In this work, we present a QM/MM study of the Mn₄CaO₅ cluster of PSII carried on at the DFT + U level for both the low-spin and the high-spin proposed structures. Our extensive analysis will reveal the role of the temperature and of the protein field in the modulation of the magnetic interactions in the PSII catalytic core.

The portion of the system close to the active site was treated by DFT as an isolated cluster (GP calculations) or immersed into its protein environment (QM/MM calculations). As suggested by Pantazis et al [18], for both calculations we have considered two different configurations of the cluster, which are hereafter referred to as model-A and model-B, differing by the oxidation and spin pattern of the Mn ions as detailed below.

2.1. Gas-phase model

The GP model was defined consistently with the DFT calculations that were previously reported for the S₂ state [18, 31, 38]. The system is shown in left panel of figure 1, and it is composed by 238 atoms: the metal core (Mn₄CaO₅), the side chains of its ligands (Asp170, Glu189, His332, Glu333, Asp342, Ala344, and CP43-Glu354), their neighbors residues (Asp61, Tyr161, His190, CP43-Arg357 and His337), the two water molecules coordinated to the Ca²⁺, a H₂O plus an OH⁻ bound to the terminal Mn4 (see also figure 2), and the closest eight crystallographic water molecules (not shown for clarity in figure 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In agreement with a previous study [18] we defined and investigated two different conformers of the catalytic core in S₂ state (see figure 2), the GP model-A (GP-A) and the GP model-B (GP-B). The geometries of these two models differ from each other by the position of an oxo-bridge (O5) which is bound in one case to Mn4 (model-A) and in the other case to Mn1 (model-B). Consequently, the oxidation states of the Mn1 and Mn4 in the two models differ according to the respective coordination pattern (i.e. Mn1(III) Mn4(IV) in model-A and Mn1(IV) Mn4(III) in model-B).

From the electronic structure point of view, the electron open-shell Mn₄CaO₅ cluster is composed by one Mn^III and three Mn^IV, with four and three unpaired electrons, respectively. These electrons occupy the quasi-degenerate d-like orbitals, localized close to the metallic center, and they give rise to a local high-spin configuration on each metal center. However, due to the exchange and super-exchange interactions, the alignment (parallel or antiparallel) between the resulting spins of each Mn ion is strongly sensitive to the structural parameters of the Mn₄CaO₅, thus changing from model-A to model-B. The exchange interactions between the electrons of the Mn ions and the energy spin ladder can be described at a phenomenological level by the Heisenberg Hamiltonian

$$\begin{equation} \skew3\hat{H}=-2\sum_{i,j(i>j)}J_{ij} \skew3\hat{S}_i\skew3\hat{S}_j, \end{equation} \tag{ 1 }$$

where $\skew3\hat{S}_i$ and $\skew3\hat{S}_j$ are the spin operators of the resulting magnetic momentum on the metal center i and j, whereas J_ij represents their exchange coupling constant. In order to describe the electronic states using a single Slater-determinant in DFT, we use the BS approach that leads to a good description of the electron density and allows an estimation of the exchange coupling constants J_ij [39–41]. In the four-metal core Mn₄CaO₅, one total high-spin state and seven BS states with different multiplicities have to be calculated in order to evaluate the exchange coupling constants between the four manganese atoms (see the supplementary data, available from stacks.iop.org/NJP/16/015020/mmedia).

The DFT calculations were performed using the ORCA [42] and CP2K [43] packages. In the ORCA calculations, the all-electron Ahlrichs TZVP basis set [44] and def2-TZVP/J auxiliary basis sets [45] are employed for all atoms, except for the C and H atoms that are described with SVP and SVP/J. We used the B3LYP hybrid functional [46–48], which has already been used to reproduce accurate energetics and magnetic properties for manganese complexes [29]. To investigate the effect of a polarizable environment, we also used an implicit solvent model (with dielectric constant  = 8), as implemented in [49]. In the calculations performed with CP2K, all of the atoms are described by the GTH pseudo-potentials [50, 51]. The DZVP-MOLOPT-SR-GTH basis set, optimized for molecular systems, is employed for all atoms. For the B3LYP calculations, the SZV-MOLOPT-SR-GTH auxiliary basis set is used for the Mn and Ca atoms and the cFIT3 is used for the other atoms [43].

In order to reduce the computational cost in ab initio MD simulations, we have adopted the DFT + U scheme [33, 52, 53]. This approach uses a basic density functional (i.e. the generalized gradient corrected PBE functional [54] in our case), including in the KS orbitals calculation an Hubbard correction term (inspired to the Hubbard model) that acts on the d orbitals of the Mn atoms. This Hubbard correction is system dependent and it is defined by the effective parameter U as detailed in the supplementary data. The value of the parameter U was chosen in order to reproduce the J_ij values obtained by the B3LYP functional on two reference geometries. The final value U = 1.16 eV was derived from the average of the linear interpolations on the exchange coupling constants J_ij, as explained in detail in the supplementary data.

2.2. Quantum mechanics/molecular mechanics model

On the basis of the models proposed by Pantazis et al [18], we built GP models of the two conformers that were representative of the S₂ state of the Kok's cycle proposed to be interconvertible in the PSII structure. The two geometries differ mainly by the position of the oxygen O5, which is bound to Mn4 in model-A, whereas in model-B it is bound to Mn1, as shown in figure 2. The two models have different magnetic properties, which are: an overall low-spin configuration for model-A, which is represented by a BS state of multiplicity 2, and an overall high-spin configuration for model-B, which is described by a BS state with multiplicity 6 (see the right panel of figure 2).

Upon geometry optimization, we calculated all the J_ij coupling constants using different functionals (table S2 of the supplementary data (see the supplementary data, available from stacks.iop.org/NJP/16/015020/mmedia)). The B3LYP Hybrid functional, which has been shown to fairly reproduce coupling constants in Mn dimers [29], provides a trend that is similar to that obtained following the procedure of [18]. This has been taken as reference for the choice of the parameter U for the Hubbard correction in PBE + U calculations. The optimal choice of the U value has been taken by fitting the J_ij parameters on the B3LYP reference for both model-A and model-B, following the procedure detailed in the supplementary data (see also table S3). The resulting value U = 1.16 eV is, therefore, used in all of the following DFT + U calculations.

For both the high-spin(A) and the low-spin(B) states, we report in table 1 the geometrical parameters of the optimized GP models, the optimized QM/MM models, and the averaged values along the QM/MM dynamics. Although it has been pointed out that the cluster geometry is affected by x-ray radiation damage [3, 38], in the last column distances and angles measured in the x-ray diffraction (XRD) structure [3] are shown for comparison. As expected, the main difference between the distances found in the x-ray structure and those of the A and B models arises from the Mn1–O5 and the Mn4–O5 distances. The O5 atom, which is equidistant from the two manganese ions in the crystallographic structure, in our models is clearly bound to only one manganese ion (Mn4 for model-A and to Mn1 for model-B).

With respect to the GP optimized geometries, the QM/MM relaxed structures display two main differences: a reduced distance between Mn4 and the Ca ion (about 0.4 Å in both models, and a slightly shorter distance between Mn1 and Mn4, which is reflected in a shorter distance between O5 and the Mn that is not bound to it. Whereas the first effect is attenuated by the dynamics, the second remains basically unchanged at finite temperature. As a general remark, we observe that the overall structures remain fairly stable along the simulated time. In particular, the bond distances between the Mn atoms of the cluster remain close to the optimized structure values with fluctuations about 0.1 Å, with the exception of the Mn1–Mn4 distance that experienced rather larger fluctuations due to a reduced rigidity of the Mn4 connection with respect to the cubane structure. Both the enhanced mobility of O5 and its delocalization in QM/MM simulations might suggest a possible role of the protein field to promote the transition from model-A to model-B. This issue can have important implications since the O5 oxo-bridge atom is thought to be a substrate for the molecular oxygen formation catalyzed by the PSII [66, 67]. An improved protein-assisted mobility might, therefore, be relevant to understand the molecular mechanism behind water oxidation.

The magnetic coupling constants have been calculated on the QM/MM optimized structures as well as along the MD trajectories. As shown previously [18], the calculated J coupling constants for the two investigated models are consistent with a large amount of EPR experiments. The comparison between the GP and QM/MM single point calculations are reported in table 2. Interestingly, in model-A the physically meaningful J_ij values (i.e. J_1,2 J_1,3 J_2,3 J_3,4, since the respective atom pairs are connected by oxo bridges) calculated on the QM/MM optimized structure show weaker ferromagnetic and antiferromagnetic interactions with respect to the GP model. Thermal fluctuations act on the magnetic coupling constants, increasing both the ferromagnetic and the antiferromagnetic interactions when compared to the values calculated on the optimized QM/MM structure. In the case of model-B, the differences between the J_ij values evaluated on the single QM/MM optimized structure and averaged over the QM/MM trajectory are smaller when compared to the trend found in model-A. In particular, only the J_1,3 value differs substantially when the effect of the dynamics at finite temperature is explicitly taken in account. Moreover, in the case of model-B, the QM/MM calculation, independently on the inclusion of the dynamics, decreases the magnetic interactions between Mn2 and Mn3 and between Mn3 and Mn4 with respect to the GP calculations. In particular, as shown in figure 3 and table 2, we found that J_3,4 in model-B largely depends on whether the GP models or QM/MM simulations are considered. The reasons for this effect can be attributed to the high mobility of Mn4, which is external to the cubane structure. The small variations of its distance and angle parameters can easily largely affect the value of the J_3,4 coupling constant.

Zoom In Zoom Out Reset image size

Table 2. A comparison of the exchange coupling constants calculated on the optimized GP structure (GP-OPT), on the optimized QM/MM structure (QM/MM-OPT) and averaged over 50 snapshots extracted from the QM/MM MD trajectory (QM/MM-MD). Optimizations of model-A and model-B were carried out in the respective spin ground states. J_1,4 and J_2,4 (marked in gray) should be considered physically unmeaningful because they refer to Mn ions not linked by a μ-oxo-bridge.

Apart from the average J_ij values, another important information available from a dynamical description of the system concerns the fluctuations of the coupling constants due to the thermal motions. An histogram of the coupling parameters calculated along the QM/MM trajectories is shown in figure 3. It is interesting to point out that the fluctuations of the coupling constants are remarkably large, often being about half of the average signal. The points on the tail of the J_ij are often rather off their average values, even instantaneously inverting the sign of the magnetic coupling: from antiferromagnetic to ferromagnetic (as in the case of J_3,4 in both model-B and model-A) and from ferromagnetic to antiferromagnetic (e.g. J_2,3 in model-A, J_1,2 in model-B).

These fluctuations, changing instantaneously the kind of the magnetic interaction between manganese atoms, might have important consequences on the interconversion between the two structures and on the catalytic mechanism of the Mn₄CaO₅ cluster. In particular, it is thought that a specific spin configuration for the manganese atoms, as well as for the oxygen atoms involved in the formation of the O₂ molecule, is required in order to allow the water splitting reaction [66]. A proper (dynamical) description of the spin electron density of the Mn₄CaO₅ cluster in the QM/MM scheme, which takes into account both the biological environment and the thermal fluctuations, is crucial in order to fully characterize the two structures model-A and model-B and to understand their role in the S₂–S₃ transition. On the basis of similar calculations, in a recent work [31] we have provided a framework in which it is possible to interpret several puzzling EPR experiments that identified metastable states between S₂ and S₃.

In the present work we have investigated the structural and magnetic properties of the Mn₄CaO₅ catalytic core of PSII. Thanks to the use of DFT + U framework within BS approach, we were able to describe the electronic structure and the magnetic properties of both model-A (low spin) and model-B (high spin) of the PSII in S₂ state. The two systems, simulated at room temperature, were structurally stable during the simulated dynamics. The magnetic coupling constants calculated on snapshots along dynamics show large fluctuations in both models, especially for J_3,4, which is the most sensitive to structural changes. The thermal fluctuations are so large that they induce instantaneous swapping between antiferromagnetic and ferromagnetic interactions, even in the case of largest $| J_{ij} |$ values. These effects might suggest that the relative stabilities of high-spin and low-spin surfaces in both model-A and model-B can be instantaneously inverted; therefore, facilitating surface hopping between the two spin states. Despite this broad distribution, their average values were in general close to that calculated on QM/MM optimized structures. In addition, the present work demonstrates the possibility of performing QM/MM MD simulations with an accurate treatment of the Mn–Mn interactions. The use of such methodologies opens the way to the study of the water splitting mechanisms in PSII within a framework that takes into account the effects of temperature fluctuations and of the protein environment.

The authors acknowledge funding provided by the European Research Council project no. 240624 within the VII Framework Program of the European Union. Computational resources were supplied by CINECA, PRACE infrastructure and the Caliban-HPC centre at the University of L'Aquila.