In recent years, the design of artificial metalloenzymes obtained by the insertion of homogeneous catalysts into biological macromolecules has become a major field of research. These hybrids, and the corresponding X-ray structures of several of them, are offering opportunities to better understand the synergy between organometallic and biological subsystems. In this work, we investigate the resting state and activation process of a hybrid inspired by an oxidative haemoenzyme but presenting an unexpected reactivity and structural features. An extensive series of quantum mechanics/molecular mechanics calculations show that the resting state and the activation processes of the novel enzyme differ from naturally occurring haemoenzymes in terms of the electronic state of the metal, participation of the first coordination sphere of the metal and the dynamic process. This study presents novel insights into the sensitivity of the association between organometallic and biological partners and illustrates the molecular challenge that represents the design of efficient enzymes based on this strategy.
The interaction of organometallic compounds with biological macromolecules has become a major focus of attention in several areas of chemistry and its interfaces . One of them is the design of artificial metalloenzymes obtained by the insertion of homogeneous catalysts into protein cavities; an approach with a tremendous potential in biocatalysis [2–4]. Numerous systems built on this concept have already been reported that include enzymes able to perform Diels–Alder reactions [5–6], transfer hydrogenation [7,8], sulfoxidation [9–11] and hydration  among others. In these hybrids, the molecular partnership between homogeneous and biological subsystems is reminiscent of those happening in natural haemoenzymes: the organometallic group affords the catalytic functionality, whereas the biomolecular one mainly dictates the substrate specificity and controls catalytic features such as enantioselectivity and/or regiospecificity .
Nowadays, the engineering of natural metalloenzymes is a widespread strategy. It mainly stands on the genetic alterations of the metal-binding site by a reduced number of amino acid substitutions, which leads to the modulation of their catalytic profiles [14–16]. However, the control of the structural and electronic properties of the metal centre using non-natural cofactors is more challenging because it involves a different chemobiological space from the one provided by the natural framework of evolution [1,17,18]. Unprecedented complementarities between organometallic and proteic moieties can be observed which could lead to binding and catalytic features different from the objectives pursued by enzyme designers. In other words, each partner may suffer molecular stresses which are absent in their conventional media such as unusual electronic states of the metal, metal–ligand exchanges upon binding of the catalyst to its receptor or even metal-mediated conformational changes.
One of the few artificial metalloenzymes for which the structure is reported in the Protein Data Bank (PDB code 1WZD) clearly illustrates the previous statement. This system, developed by Ueno et al. , was obtained by the substitution of haem by a Fe(Schiff base) (salophen) in the Corynebacterium diphtheriae haem oxygenase (cdHO). Natural and artificial systems share important aspects of their catalytic mechanisms, i.e. both proceed throughout an oxido-reduction mechanism with the cytochrome P450 reductase as a partner, leading to the reduction of the iron atom of the cofactor from Fe(III) to Fe(II) and the subsequent coordination of a molecule of dioxygen (scheme 1). However, the natural system performs an entire oxidative process ending at the cleavage of the haem, whereas the artificial enzyme stops at the activation of the oxygen to form superoxide ions. No subsequent oxidation of the cofactor or any external substrate proceeds. This superoxidase activity differs drastically from the natural function of cdHO; a surprising result considering the high chemical similarity between salophen Schiff bases and porphyrinic cofactors in terms of structure and reactivity.
The X-ray structure of the system provides some molecular hints on those differences. In this structure, the Fe(Schiff base) is highly distorted and has two residues of the receptor, His20 and Glu24, coordinating the metal (figure 1). Hypothesized as the resting state of the system, the iron atom is afforded a +3 oxidation state of unknown spin state and presents a distorted octahedral configuration with no room for the binding of the oxygen. This structure highly differs from those known for natural haemoenzymes including holo-cdHO. Indeed, these species generally display resting states with either a square pyramidal high-spin (i.e. catalase ) or an octahedral low-spin configuration (i.e. peroxidases , cytochromes P450  or haem oxygenase ) with only one residue of their apo-protein bound to the metal in the axial position and, eventually, one labile water molecule coordinated on the distal site of the haem to complete the octahedral configuration of the metal (figure 1). The strength of the coordination of this water molecule can be influenced by its environment and, in turn, conditions the transition from hexacoordinated resting states to active pentacoordinated species .
The experimental geometry of Fe(Schiff base) · cdHO therefore suggests that both the resting state and the activation process could be different from those of naturally occurring Fe(III) haemoenzymes. Its divergence from the resting states of natural ones also questions the possibility of a crystallographic artefact; the known flexibility of apo-HO and salophens could be accentuated in the crystallographic media. A better understanding of the complementarities between the organometallic and the proteic moieties in Fe(Schiff base) · cdHO could represent a major step forward in bioinorganics and further help in the design of artificial metalloenzymes and other bioinorganic architectures.
Computational chemistry is now widely used for the study of metalloenzymes. Most of its application focuses on unravelling the molecular grounds of their catalytic mechanisms . Little has yet been done to decode the relative contributions of inorganic and biological partners and less still for enzyme design. This is due to the complexity of dealing with the binding of metal-containing systems to proteins by standard approaches . Simulations involving large conformational samplings as well as accurate representation of the metallic environment are indeed necessary. We recently showed that a combination of structural bioinformatics, quantum-based methods (quantum mechanics (QM) and quantum mechanics/molecular mechanics (QM/MM)) and protein–ligand dockings represents an interesting solution in characterizing low-energy complexes of organometallic systems bound to proteins [27,28]. Here, we further investigate the molecular complementarity between artificial organometallic cofactors and proteic partners by performing an extensive QM/MM study on the experimental Fe(Schiff base) · cdHO system, focusing on (i) the nature of its Fe(III) resting state, its comparison with those of natural haemoenzymes and its possible artefactual nature, (ii) the activation process, with the simulation of the full transition occurring between hexa- and pentacoordinated forms and considering the reduction of the metal centre, and (iii) the relative contributions of the inorganic and biological parts to the structural and energetic properties of the biocomposite. Additionally, this study illustrates the potential of QM/MM approaches to understand novel chemobiological processes.
2. Material and methods
Calculations have been performed on the entire Fe(Schiff base) · cdHO system using the QM/MM ONIOM (DFT : AMBER) method implemented in Gaussian 09 . The charges and protonation states of all titrable amino acids were automatically assigned using the interface provided by the UCSF Chimera package  with the exception of the iron-chelating histidine His20, which was manually set to be consistent with the coordination rules of the metal. Visual inspection was subsequently performed. The total charge of the system is −8 or −7 depending on the oxidation state of the iron. It can be divided by −10 for the isolated protein and +2 or +3 for the inorganic complex. The QM part in the QM/MM partition has a charge +1 in the Fe(II) species and +2 for the Fe(III) ones.
One of the main objectives of this work is to determine whether the experimental structure corresponds to a plausible electronic structure of the metal centre in a resting state configuration or not. Thus, we are interested in keeping the initial structure for our study as close as possible to the crystallographic one. With this purpose in mind, and considering that the absence of optimum force field parameters of the iron in the different coordination spheres of the salophen under consideration may lead molecular dynamics towards less realistic equilibrium structures, no such experiments were performed on the initial X-ray structure . Moreover, as we were also interested in analysing the transition between hexa- and pentacoordinated forms of the Fe(Schiff base) · cdHO complex, an accurate electronic representation of the first and second coordination spheres of the metal was mandatory, and QM/MM calculations were chosen as the best method for this purpose.
For the hexacoordinated system, calculations were carried out on the original crystallographic X-ray coordinates corresponding to this state (referred to as 6, PDB reference: 1WZD). The equatorial environment of the metal results from an N2O chelation provided by the Schiff base and the Oγ of the Glu24 of the host; the remaining axial positions are occupied by the second oxygen of the cofactor on the distal side and the Nε of residue His20 on the proximal one (figure 1). For the pentacoordinated systems with or without an axial water (referred to as 5′ and 5, respectively), the equatorial environment is entirely provided by the cofactor bound in an N2O2 fashion, and the axial proximal position is occupied by the His20. In these structures, the side chain of Glu24 has been removed from the first coordination sphere of the metal and moved towards the most probable rotameric state for this residue using the Dunbrack rotamer library of the UCSF Chimera package [3,2]. A fourth system, referred to as 6′, was built to allow energetic comparison with 5′. It consists of a similar system to 6 but with a water molecule added inside the binding site of the biocomposite in the nearest hole in the vicinity of the cofactor (between Lys13 and His25). In all these systems, the cofactor remains in a hydrophobic environment and packed between helices A and F of apo-cdHO. The water molecules that could interfere with the mechanism would therefore be on the superficial region of the protein and would tremendously reorganize during the simulation of transition between penta- and hexacoordinated structures. This would lead to a highly unsmoothed potential energy landscape that is difficult to study and with relatively little improvement in the description of the chemical process at the metal site. Therefore, other molecules than those involved in the binding of the iron were removed from the calculations [3,1].
Two QM/MM partitions were used. The larger partition includes the metal, the side chain of the coordinating residues His20 and Glu24 up to their Cα atoms and the entire aromatic part of the Schiff base (electronic supplementary material, scheme S1). The second partition contains the metal, the same coordinating residues of the protein up to the Cα and a model of the first coordination sphere provided by the salophen where the aromatic moieties are replaced by simpler ethylenic groups (electronic supplementary material, scheme S2). In both cases, the rest of the cofactor was included in the MM part of the system and treated using the AMBER force field  for the amino acids and the GAFF  one for the atoms of the cofactor. In the first partition, 55 of the 3347 total atoms are included in the QM region, and in the second one, only 37. The smaller partition was used in the post-analysis and focused on the contribution of the first coordination sphere of the metal to the mechanism of the enzyme. To maintain the geometry of the system close to the experimental system, only atoms (QM and MM) within a sphere of about 15 Å from the iron where allowed to relax during the optimization process with the entire helix A included in this flexible scheme. In both partitions, the number of flexible atoms was 192 atoms and the rest remained fixed.
For each structure, QM/MM minimizations were undertaken for low-, intermediate- and high-spin Fe(II) (multiplicity values: 1, 3 and 5) and Fe(III) (multiplicity values: 2, 4 and 6) species. Optimizations have been performed with the B3LYP functional [35,36] using the 6–311+g* basis set [37,38] for the main group elements and the aug-cc-VTZ [39,40] one for the iron. Single point calculations on the optimized structures have been performed with the M06L , PBE  and B97D  functionals to identify possible drawbacks on the relative stability between different spin states [44–46]. For the same reason, calculations using electronic embedding have been also tested.
3. Results and discussion
3.1. Analysis of the Fe(III)(Schiff base)·HO resting state
Calculations have been carried out on the crystallographic system (referred to as 6) and on configurations generally observed in natural haemoenzymes. The latter correspond to, on one side, the structure obtained by the removal of the glutamate from the first coordination sphere of the metal (leading to five coordinated iron and referred to as 5) and, on the other, to the configuration obtained when substituting the glutamate by a water molecule (leading to an alternative six coordinated iron and referred to as 5′) (scheme 2). All the optimizations were performed with the larger QM/MM partition described in the previous section.
Independently of their spin states, all optimized structures of the ferric systems 6 are in good agreement with their crystallographic counterpart (electronic supplementary material, table S1). The iron remains in an octahedral configuration with a distorted cofactor and coordinated to the protein by His20 and Glu24 (figure 2; electronic supplementary material, figures S1 and S2). In all cases, bond lengths and angles of the first coordination sphere of the metal are in good agreement between theory and experiment, but closer values are obtained for the Fe(III) high-spin system (generally lower by 0.1 Å with the exception of interaction Fe–Oglu, which increases by 0.2 Å for distances and 0.3° for angle values).
When removing the glutamate from the first coordination sphere (5), the system adopts a square pyramidal geometry (electronic supplementary material, table S1). The cofactor becomes more planar than in the experimental system with deviations from an ideal plane (OOP Fe values in the electronic supplementary material, table S1) of 0.2, 0.4 and 0.6 Å for low, intermediate and high spin, respectively, instead of the 0.7 Å observed for the experimental system. Overall, the same tendencies are also obtained for the substitution of the glutamate by a water molecule, being the cofactor slightly more planar in 5′. Nonetheless, it must be noted that the Fe–(Owater) bond is almost cleaved in the intermediate- and high-spin complexes, which is consistent with previous observations in cytochromes P450 (electronic supplementary material, figure S2 and table S2) .
In addition to the general changes in the first coordination sphere of the metal and the reorganization of the cofactor, the removal of the glutamate is also correlated with a reorganization of helix A. Optimized structures of 5 and 5′ show the Cter end of the helix pointing towards the solvent in a conformation similar to those observed in the haem-bound haem oxygenase structure . This contrasts with the structures of haem-bound systems in which the macrocyclic nature of the cofactor prevents the Glu24 from getting close enough to the Fe(III).
Relative stabilities between the different optimized geometries of the Fe(III)(Schiff base) · cdHO show a clear preference for the experimental configuration (6) with respect to the ferric 5 and 5′ structures (figure 4; electronic supplementary material, figure S7). Moreover, calculations performed with the (B3LYP : AMBER)-mechanical embedding (ME), (B3LYP : AMBER)-electronic embedding (EE), (M06L : AMBER)-ME and (B97D : AMBER)-ME level are consistent in predicting 6 to be about 5–10 kcal mol−1 more stable in its high-spin configuration. Only the (PBE : AMBER)-ME scheme diverges, giving isoenergetic low- and high-spin wave functions in the case of configurations 6 and 5′, a behaviour that has been associated with the weakness of this approach to discriminate between spin states of different multiplicity . Finally, the differences in energy between 6 and 5 is comparable to those observed between 6 and 5′ but larger for the spin state of higher multiplicity. This has been associated with the weakening of the Fe–(Owater) bond that could lead to a dissociation (electronic supplementary material, figure S2).
To further investigate the ferric system, the full transition paths between configurations 6 and 5 have been simulated. Because the use of electronic embedding does not present major improvement and taking into account previous considerations regarding the quality of the different functionals, the following part of the work was carried out under the (B3LYP : AMBER)-ME scheme.
The transition state structures have been characterized at Fe–Oglu distances of 3.01, 2.55 and 3.37 Å for the low-, intermediate- and high-spin species, respectively (see figure 3 and electronic supplementary material, table S3, for more details). No major differences are observed in the overall structures of the different transition states and a unique negative vibrational mode has been identified. In the first coordination sphere of the metal, this mode shows major components in the breaking of the Fe–Glu24 bond as well as in the displacement of the His20 towards a perpendicular position relative to the cofactor and passing through the metal and the displacement of the glutamate towards the solvent (figure 3). Major contributions are also observed for the cofactor, which bends to get close to planarity, as well as for helix A, whose Nter region gets closer to the Schiff base and whose Cter region gets more exposed to the solvent. This result clearly shows how the distortion of the cofactor, the modification of the first coordination sphere of the metal and the rearrangement of the structure of the receptor are intrinsically correlated and not individual phenomena.
The energy of the transition state is 20.7, 9.8 and 18.6 kcal mol−1, above the absolute minimum of the system, for low-, intermediate- and high-spin systems, respectively (figure 4). They are sufficiently large to discard the possibility that the atomic motions at room temperature naturally allow the transition from glutamate bound (6) to glutamate unbound configurations (5 and 5′). Interestingly, the energetic breakdown in QM and MM contributions of the total QM/MM energy shows that the electronic (QM) term dominates in the relative stability of 6, 5 and 5′ (electronic supplementary material, figure S3). Therefore, the changes in the structure of the receptor at helix A and the conformation of the salophen are mainly driven by electronic and geometric properties of the first coordination sphere of the metal.
Our calculations on this part of the study support that the experimental configuration of the Fe(Schiff base) · cdHO observed in the crystallographic structure is the real resting state of the enzyme and that its electronic state corresponds to a high-spin configuration. The large differences in energy between the experimental configuration and other possible resting state geometries imply that dynamical equilibrium between them does not naturally take place in solution. In this, the constraint on the system is manly associated with the first coordination sphere of the metal.
As, in this configuration, neither labile ligand nor sufficient room is available on the top of the iron to coordinate the oxygen, the activation process of the enzyme most likely would be triggered by the reduction of the metal and would differ from the chemistry of haemoenzymes.
3.2. Electronic transition and activation mechanism
The previous part of the study reveals that the transition from glutamate bound (6) to glutamate unbound (5 and 5′) metal is unreachable for Fe(III) species. In this configuration, the oxygen cannot reach the iron and bind to it because of the absence of room in the upper position of the cavity placed in the axial position of the salophen. Therefore, the reduction of iron is necessary for the system to be activated and the reaction to proceed (scheme 1). To understand the activation mechanism of the enzyme, the entire transition paths from 6 to 5 have been simulated for Fe(II). The impact of the solvent molecule (5′) has only been estimated on the minima.
Optimizations of the Fe(II) systems lead to stable 6 and 5 conformations for all spin states. The general structural profiles of these species are similar to those observed in the ferric structures. Metal ligand distances of the first coordination sphere are slightly longer for the ferrous system. The position of the glutamate with respect to the metal is particularly affected, being approximately 0.1 Å closer for Fe(III) than Fe(II) systems at identical spin states. This is consistent with stronger electrostatic interactions in the former. Importantly, for 5′, the water molecule remains coordinated to the iron in the axial position in the low-spin state Fe(II) only (see electronic supplementary material, table S3). For intermediate- and high-spin species, the water molecule leaves the first coordination sphere of the metal and lies at 3.33 and 3.56 Å, respectively. Such spontaneous removal of the solvent molecule is compatible with the creation of a vacant site on the top of the iron, as observed in haemoenzyme systems . Therefore, simulation of transitions from 6 to 5 should be indicative enough of the activation mechanism.
For (B3LYP : AMBER)-ME calculations, the transition between 6 and 5 is exothermic for intermediate and high spin and slightly endothermic for the low spin (figure 5; electronic supplementary material, figure S6). The high-spin state is substantially more stable than the two others, which suggests that the process should take place entirely in this spin state. These results are qualitatively supported by (M06L : AMBER)-ME, (B97D : AMBER)-ME and (B3LYP : AMBER)-EE (electronic supplementary material, figure S7). The (PBE : AMBER)-ME approach slightly diverges from the other functional by suggesting possible spin crossing between low- and intermediate-spin wave functions. Probably due to the weakness of this functional in dealing with the difference in energy between spin states, these results still have no significant impact on the chemistry of the system. For the transition from 6 to 5′, the high-spin wave function remains the most probable spin state for the transition but the differences in energy observed are lower than for the transition between 6 and 5 (electronic supplementary material, figure S6). This is probably due to the presence of a water molecule in 5′ that can stabilize the system. Taken together, the stable geometries of penta- and hexacoordinated structures of the ferrous systems show that a conformation adequate for the binding of the oxygen is favoured after reduction of the iron and produce a stable high-spin configuration. This is strongly different from what is known of naturally occurring haemoenzymes. For example, in cytochromes P450, the entrance of the substrate in the binding site leads to the removal of the labile water ligand and the electron transfer to generate the ferrous species.
The transition path was therefore investigated. As the water chelating the iron in 5′ is dissociated in the high- and intermediate-spin function, further calculations were carried out only on the transition between 6 and 5. Transition states between these two structures occur at Fe–Oglu distances at 3.44, 2.73 and 3.02 Å for low-, intermediate- and high-spin species, respectively. The general features of the lowest vibrational mode are similar to those of the ferric system as reported in the previous section. Additionally, the transition states for these different species are located at 13.4, 2.5 and 10.1 kcal mol−1 above the reactant structure for low, intermediate and high spin, respectively. For the high-spin species, which represents the ground state of the system, this means a diminution of about 8 kcal mol−1 and a value relatively accessible in standard conditions. It appears therefore that, when the system has been reduced, the transition between 6 and 5 has a low enough barrier to take place.
Enthalpic and kinetic considerations confirm that, once the electron has been transferred to the ferric ion, the removal of the glutamate from the first coordination sphere of the metal becomes favourable. The mechanism should take place via a high-spin state throughout the process and lead to stable square pyramidal structures compatible with oxygen binding at the vacant site of the metal. This would take place through a transition mechanism where a mixed contribution of the conformational change of the cofactor, the rearrangement of helix A and the change in coordination state of the metal occurs. Which of these features dominates in the process is yet to be determined but is fundamental to decoding the real extent of non-natural bioinorganic complementarities.
3.3. First coordination sphere of metal versus protein–cofactor complementarities
The final question we aim to answer is to identify the molecular variables that condition the unexpected geometrical features of the Fe(Schiff base) · cdHO system and so shed light on what are the energetic contributions that dominate in the formation of bio-organometallic systems. We therefore embarked on analysing the relative energy of the system, partitioning it in terms of (i) first coordination sphere of the metal versus the rest of the system and (ii) ligand versus protein (electronic supplementary material, scheme S2). To do so, we performed additional calculations where the QM/MM partition was constructed in such a way that the first coordination sphere of the metal was modelled in the quantum mechanics region, while the rest of the system was included in the MM one (electornic supplementary material, scheme S3). We first benchmark this novel partition by carrying on calculations on the transition from 6 to 5 and comparing these results with the larger one. Despite slight geometric and energetic nuances, the overall profiles are very similar (figure 6; electronic supplementary material, figure S4) and further analyses were therefore performed on these calculations.
Decomposition between QM and MM terms of the total QM/MM energy unambiguously shows that the QM term dictates the shape of the overall profile in all oxidation and spin states. The QM curves are only a few kcal mol−1 higher than the QM/MM ones and both of them can almost overlap. The MM terms have a minor role and only modulate the shape of the QM/MM profile. However, it is still interesting to note that, in most cases, this part of the system remains almost unaffected (Fe(II) cases) or even stabilized (Fe(III) cases) for Fe–Oglu going from reactant to a few steps after the transition state. From then on, a destabilization of the system is observed, reaching up to ca 15 kcal mol−1 (electronic supplementary material, figure S4). Structural analysis showed us that the first part of the transition implies a slight relaxation of helix A while the glutamate is removed. However, at longer distances some clashes at the hinge between the flexible region of our partition and the rigid one are observed, hence leading to higher energies. Therefore, our partition does not allow us to deal with the full extent of the molecular flexibility engendered by the relaxation of the first coordination sphere of the metal; something that definitely represents a major tour de force in molecular modelling. However, our results are still quite indicative. Once the Glu24 goes out of the first coordination sphere of the metal, the rearrangement of helix A can be energetically meaningful.
The larger partition allows us to consider the relative contribution of the cofactor and receptor in the total energy of the complex. In this case, a clear different pattern is observed. For all spin and oxidation states, the energy of both isolated cofactor and receptor decrease as a function of the distance between the glutamate and the metal (figure 7; electronic supplementary material, figure S5). The protein tends to stabilize to ca 10 kcal mol−1 and the cofactor up to 14 kcal mol−1. This means that a planar conformation of the Fe(Schiff base), the presence of the glutamate out of the first coordination sphere of the metal as well as a conformation of helix A with its Cter end pointing towards the solvent would be preferred by considering the energetic properties of both individual species. All these statements are consistent with current knowledge on porphyrinic, salen and salophen moieties  as well as the crystallographic structure of haem-bound haem oxygenase .
This part of the study clearly shows that the structure observed experimentally corresponds to a stable Fe(III) distorted conformation of the cofactor with the Glu24 coordinating the metal. This structure can only be understood if considering that the coordination rules of the metal drive the tuning between the inorganic and proteic moieties. The Fe(III)Schiff base moiety, by being smaller and more flexible than the haem in the natural system, affords an additional coordination site that can be reached by a residue acting as a Lewis basis. Hence, it appears that all possible conformational changes in the cofactor or the proteic receptor that could satisfy the best coordination of the metal would take place even if both chemical and proteic systems should be penalized.
In this study, we embarked on the analysis of a novel resting state and investigated the nature of the activation process of a haemoenzyme to finally discuss the results in terms of physico-chemical properties governing the interaction between inorganic moieties and proteins. This study clearly shows that, despite the structural similarity of the artificial cofactor to the naturally occurring haem, the structure reported by Ueno et al. represents indeed the conformation of the resting state of the enzyme and corresponds to a high-spin system. Despite the well-known tendency of salophen to adopt planar geometries and the apo-HO system to coordinate the haem through a unique histidine, the molecular rules that govern the structure of the resting state of the novel enzyme and its activation process are imposed by the first coordination sphere of the metal. The removal of the side chain of the Glu24 needs though a synergic mechanism driven by the first coordination sphere of the metal but implying an overall structural change in the relative position of the proximal helix and the conformation of the cofactor. These results show that foreseeing the degree of conformational variability of the inorganic cofactor and the proteic receptor and any possible ligand exchange that could take place in the coordination sphere of the metal is a significant variable for the design of novel biometallic hybrids. This work sheds light on the major impact of the metal in dictating the structure of bioinorganic composites and provides key information for the rational design of novel members of this family of catalysts.
We are particularly thankful to the Spanish ‘Ministerio de Economía y Competividad’ for financial support through projects CTQ2011-23336 and ORFEO Consolider-Ingenio 2010 Programme (grant no. CSD2007-00006), the Generalitat de Catalunya through project 2009SGR68, and E.O-C. thanks the Universitat Autònoma de Barcelona for her scholarship (UAB-PIF).
- Received January 27, 2014.
- Accepted April 17, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.