A computational study of ligand binding a ffi nities in iron ( III ) porphine and protoporphyrin IX complexes †

The search for novel anti-malarial drugs that can disrupt biomineralization of ferriprotoporphyrin IX to haemozoin requires an understanding of the fundamental chemistry of the porphyrin’s iron(III) centre at the water–lipid interface. Towards this end, the binding affinities for a diverse set of 31 small ligands with iron(III) porphine have been calculated using density functional theory, in the gas phase and also with implicit solvent corrections for both water and n-octanol. In addition, the binding of hydroxide, chloride, acetate, methylamine and water to ferriprotoporphyrin IX has been studied, and very similar trends are observed for the smaller and larger models. Anionic ligands generally give stronger binding than neutral ones; the strongest binding is observed for RO and OH ligands, whilst acetate binds relatively weakly among the anions studied. Electron-rich nitrogen donors tend to bind more strongly than electrondeficient ones, and the weakest binding is found for neutral O and S donors such as oxazole and thiophene. In all cases, ligand binding is stronger in n-octanol than in water, and the differences in binding energies for the two solvents are greater for ionic ligands than for neutrals. Finally, dimerization of ferriprotoporphyrin IX by means of iron(III)–carboxylate bond formation has been modelled. The results are discussed in terms of haemozoin crystal growth and its disruption by known anti-malarial drugs.


Introduction
Although malaria is both treatable and preventable, it remains one of the most prevalent infectious diseases worldwide, imposing an intolerable burden of disease and death across endemic countries.An effective vaccine has yet to be developed, and strategies for preventive chemotherapy are complicated by the emergence of parasite strains that are resistant to currently available drugs, including artemisinin. 1 Tackling this situation requires a broad approach, including the targeted use of a range of chemotherapeutics such as quinine (QN), quinidine (QD) and chloroquine (CQ).The disease stage of malaria infection arises when the Plasmodium parasite is growing in the host's red blood cells, and a key aspect of the mechanism of action of both QN, CQ and related antimalarials is accepted to be their disruption of the parasite's ability to detoxify the haem released from digestion of the host's haemoglobin. 2Normally, at least 95% of the iron(III) ferriproto-porphyrin IX (FePPIX, ‡ Scheme 1) released from haemoglobin is biomineralized by the parasite into insoluble haemozoin crystals; reducing the efficiency of this process kills the parasite. 3,4Crystallization of haemozoin does not occur spontaneously under normal physiological conditions, but is promoted in the digestive vacuole of the parasite by a combination of low pH (ca.4.8) and the presence of lipids; in particular, the water-lipid interface is heavily implicated as the site of crystal growth. 4,5he powder X-ray diffraction structure of β-haematin, which is the synthetic equivalent of haemozoin, shows that the structure is composed of FePPIX dimers, in which each iron centre is bonded to a propionate group from its partner's porphyrin (Scheme 1).The dimers are in turn linked together by cyclic hydrogen bond pairs between the remaining neutral propionic acid groups. 6It has been suggested from quantitative estimates that CQ and related drugs inhibit haemozoin crystallization by surface adsorption onto the growing crystal faces, 6 and the possible binding modes of a wide range of antimalarials to the various crystal faces of haemozoin have been reviewed. 3lthough our understanding of the structure-function relationships of drugs that inhibit haemozoin crystal growth is still incomplete, the X-ray crystal structures of the iron(III) species of FePPIX complexed with halofantrine, 7 QN and QD 8 have provided evidence that these drugs can bind directly to the iron(III) centre.In all three structures, the antimalarial drug is bound to the metal site through an alkoxide oxygen; there are also hydrogen bonds between protonated amine groups on the drugs and the FePPIX propionate groups, and π-π stacking interactions.In contrast, there is currently no published crystal structure of a FePPIX-CQ complex.X-ray absorption spectroscopy studies have suggested that in DMSO and acetic acid solutions, CQ forms weak complexes with the FePPIX dimer and monomer respectively; in DMSO, the Fe atom is closest to (but not coordinated by) the quinoline nitrogen. 9s a potential drug target, haemozoin is very unusual in that it has a relatively simple crystalline structure, in marked contrast to the much more elaborate structures of more common drug targets such as proteins and nucleic acids.The inhibition of haemozoin crystal growth can be expected to arise primarily from some combination of a few specific interactions between the drug and FePPIX (whether in monomeric, dimeric, or crystalline form), as follows.First, the strongest interactions would arise from covalent attachment of the drug to the haem.This has been demonstrated for artemisinin, which can alkylate FePPIX at the meso-positions. 10Second, a drug may coordinate to the iron atom, as observed in the crystal structures of halofantrine, 7 QN and QD. 8 Third, a drug might form specific hydrogen bonds with FePPIX, either in its molecular form or on the surface of haemozoin crystals.The options for classical hydrogen bonding are however limited to the propionate groups; non-classical hydrogen bonds to the hydrocarbon parts of the molecule should be much weaker.Fourth, the planar haem lends itself to π-π stacking interactions; these are likely to be important for drugs containing aromatic groups, such as CQ.The energies of π-π stacking interactions cover a similar range of energies as for hydrogen bonds; for example the benzene dimer interaction energy is approximately 2.6 kcal mol −1 , 11 whilst the enthalpy of π-π stacking in a zinc porphyrin system has been measured as 11.5 kcal mol −1 . 12Experiments on substituted porphyrins gave ΔG values for their π-π interactions with aromatic systems of 6, 10 and 14 electrons of ca.1.7, 3.8 and 4.4 kcal mol −1 respectively. 13For non-alkylating drugs (i.e.most drugs apart from artemisinin and its derivatives), some combination of coordination to iron, hydrogen bonding and/or π-π stacking interactions must be sufficiently favourable to disrupt the normal process of haemozoin crystal growth at the water-lipid interface.
In principle, quantum calculations can provide an estimate of the strength of binding of different drugs to FePPIX.However, although haem-based molecules have been extensively studied by means of quantum calculations, the large size of the porphyrin ring and the complications arising from accurate description of the iron atom and its spin state have made the calculation of reliable ligand binding energies a challenging task.5][16][17] In particular, papers by Wondimagegn and Rauk have described the properties of iron(II) and iron(III) porphine complexes in the context of modelling the possible roles of FePPIX in Alzheimer's disease. 14They considered a range of small ligands inspired by the naturally available amino acid side chains, including implicit corrections for water and benzene as solvents.
In view of the pressing need to understand the mode of action of known inhibitors of haemozoin crystallization, and to use that knowledge to develop novel antimalarials, the present paper describes the results of density functional theory (DFT) calculations on the binding of a diverse set of 31 different ligands with iron(III) porphine.The importance of the water-lipid interface in haemozoin crystallization has been taken into account by the use of implicit solvent calculations for both water and n-octanol, using the solute electron density (SMD) implicit solvent model. 18The calculations have also been extended to the binding of a few ligands to the complete iron(III) FePPIX structure, including estimates of the pK a 's of the two propionate ligands and coordinated water.The results are discussed in terms of our current understanding of haemozoin crystallization, and its inhibition by drugs such as QN and QD.

Methodology
The initial goal of this work was to use DFT calculations to provide a reasonably accurate description of the known chemistry of FePPIX.Here, it is helpful to reduce the system down to the porphine complex (3) in Scheme 1.Not only does this speed up the calculations, but it also avoids the complications arising from the chemistry of the propionate groups; in particular, these are capable of three different ionization states, as well as a range of conformations and intramolecular hydrogen bonding (these matters are dealt with in more detail below).We have previously found that the B3LYP functional and LanL2DZ basis set (denoted below as method 1) is generally sufficient to provide useful semi-quantitative models for a range of transition metal complexes. 19However, although this level of theory is adequate for the comparison of relative bond energies among related species, it provides rather inaccurate absolute values of bond energies.Moreover, correct prediction of the spin ground state of iron complexes is a difficult problem; in particular, the B3LYP functional is generally considered to be biased toward high spin configurations of transition metal complexes, and consequently may perform less well than more recent functionals such as OPBE for spin state prediction. 20Therefore, additional single point calculations have been carried out with the OPBE functional and triple-ζ 6-311+G(d,p) basis set (method 2), which preliminary calculations indicated was the largest basis set to give tractable calculations for the largest molecules considered in this work.Finally, single point calculations have also been undertaken at the B3LYP/6-311+G(d,p) level of theory (method 3).The relative strengths and weaknesses of different DFT approaches in the particular context of bioinorganic chemistry have recently been reviewed by de Visser et al., 16 who found that the combination of B3LYP with a basis set of triple-ζ quality and solvent corrections can give reliable energies for biologically relevant iron complexes, provided that the calculations are calibrated with experimental data (see below).
For each iron(III) complex, initial geometry optimizations were carried out in vacuo using method 1 for all three possible spin states (S = 1/2, 3/2 or 5/2).These were followed by frequency calculations at the optimized geometries to check that the optimized structures were true minima and to obtain zero point energy (ZPE) values, plus additional single point calculations using the SMD implicit solvent model 18 to obtain energies for the molecules in both water and n-octanol.All reported energies include ZPE corrections (see the ESI † for raw energy and ZPE values).According to these preliminary method 1 calculations, the S = 1/2 spin state was never the ground state; the closest approach of the S = 1/2 state to the ground state was found for the in vacuo calculations on the parent iron(III) porphine (0.1 kcal mol −1 higher in energy), and its imidazolide and methanethiolate complexes (both 0.8 kcal mol −1 higher in energy).Therefore, these species were also included in subsequent OPBE (method 2) calculations.However, the method 2 calculations agreed that the S = 1/2 state was not the ground state for these molecules, being at least 3.8 kcal mol −1 above the ground state.Therefore, the S = 1/2 spin state was not considered further in the method 2 and method 3 calculations.For most species, the method 2 OPBE calculations predicted the same ground state as the method 1 B3LYP calculations, but for seven of the ligands there were some disagreements; in these cases, the method 2 calculations invariably favoured the higher S = 5/2 spin state when compared to the method 1 calculations, which predicted an S = 3/2 spin state.This outcome is contrary to the normally stated view that B3LYP is biased toward higher spin states; however, the choice of basis set is also known to be an important factor, 20 and it would seem that the solvent correction is also important, since almost all of the disagreements were found for the calculations that included solvent corrections.When methods 1 and 2 both predicted the same spin state, only this state was included in the method 3 calculations; if there was a disagreement, both the S = 3/2 and 5/2 calculations were included.In these cases, the ground state predictions for method 3 were evenly spread between agreement with methods 1 and 2.Although the well-known difficulties associated with prediction of the relative energies of different spin states in haem complexes by DFT and other quantum methods should always be kept in mind, this does not negate their value for the prediction of reaction energies, especially when used in a comparative mode. 17omparing the complete sets of binding energy data for all the ligands in this study, including both the porphine and porphyrin models, with and without solvent corrections, it is apparent that method 1 predicts the most favourable binding energies, method 2 the least favourable, and method 3 gives intermediate values (for example, see Table 1).Although the absolute calculated binding energy values are contingent on the method, there are good linear correlations between the data sets; plots of all data for each pair of methods gave correlation coefficients R 2 of 0.987 for method 1 versus method 2, 0.985 for method 1 versus method 3, and 0.994 for method 2 versus method 3, indicating that all three methods agree quite well in predicting the relative binding energies of different ligands, notwithstanding occasional disagreements over the ground spin state.
Next, the theoretical results obtained in this study have been compared with available experimental data.Although a number of quantitative studies of the binding of exogenous ligands to iron(III) haem sites have been reported, 21,22 it has proved difficult to obtain experimental data for monomeric FePPIX species with only a single axial ligand.Nevertheless, this has been achieved in a recent paper by Kuter et al., 21 who have reported pH-independent association constants for FePPIX with a number of ligands, specifically identified as coordinating to iron rather than π-stacking.In their report, these authors found that 4-dimethylaminopyridine appeared to be an outlier in a plot of log K obs versus pK a of the free ligand, and initial comparison with the calculated results from the present work also suggested that this compound is anomalous.An explanation is provided by the fact that the experiments were carried out in 40% aqueous DMSO, whereas the pK a values used were those pertaining to water.Other workers have measured the pK a values of pyridine, 4-methylpyridine and 4-dimethylaminopyridine, 23 of butylamine and morpholine, 24 and of imidazole 25 in water-DMSO mixtures.Using these data to calculate corrected pK a values and so replot the graph of experimentally obtained log K obs versus pK a resolved the discrepancy, and also improved the agreement between the observed and calculated data.Therefore, the experimentally derived ΔG values used in the following analysis have been calculated using the solvent-corrected pK a values.These are compared with calculations on the equivalent porphine complexes in Table 1.As usual, the calculated ΔE values include ZPE's; for all three computational methods, omitting the ZPE's gave somewhat poorer correlations with the experimental data, as did the use of full thermal corrections to provide calculated values of ΔH.Interestingly, method 2 using the OPBE functional gave the poorest correlation with the experimental data (R 2 = 0.699); the best correlation was obtained with method 3 (R 2 = 0.865).Moreover, the values of ΔE obtained with method 3 were closest to the experimental ΔG values.In all three cases, including the calculated entropies to give a theoretical estimate of ΔG resulted in a poorer correlation with the experimental data; the calculated entropy terms for the five ligands were all unfavourable for ligand binding, and ranged from 11.7 to 12.6 kcal mol −1 at 298 K.It should be noted, however, that these values do not include the contribution from the surrounding solvent molecules, which would probably offset much of the entropic cost of ligand coordination.If it is assumed for the sake of analysis that the ΔE values obtained with method 3 are close to the true enthalpy values, empirical fitting of the data suggests that the overall entropy term for coordination of these ligands is in the region of +1 to +3 kcal mol −1 at 298 K.

Iron(III) porphine complexes
Having identified method 3 as giving the best correlation with experiment, the following discussion refers to this set of calculations unless mentioned otherwise.The next goal is to understand the coordination chemistry of the iron(III) centre.In general, methods 2 and 3 both tended to favour the high spin S = 5/2 state for anionic ligands, and the S = 3/2 spin state for neutral ligands; the only exceptions were chloride and imidazolide (S = 3/2 ground state predicted by method 3 in water and octanol), and Me 3 NO and Me 3 PO (for which some of the solvent jobs predicted the S = 5/2 ground state).In all cases, the vacuum calculations identified the more typical spin states as the ground states.According to method 2, the largest energy difference between the S = 3/2 and 5/2 spin states was found for the parent porphine in a vacuum (S = 3/2 more stable by 11.5 kcal mol −1 ); the difference was reduced by both the presence of an exogenous ligand and by the solvent correction, such that for several species, the higher and lower spin states were separated by less than 1 kcal mol −1 .Hence, the energy differences between the S = 3/2 and 5/2 spin states are often within the error margin of the calculations.Taking the data from all both methods 2 and 3 into account, the Mulliken electron spin densities on iron for the vacuum jobs covered discrete ranges of 2.6-3.0 and 4.0-4.3 for the S = 3/2 and 5/2 spin states respectively; the remaining spin was largely divided among the atoms directly coordinated to iron.
Analysis of the optimized geometries shows that the iron atom is displaced out of the plane of the four haem N atoms upon coordination to a fifth ligand; the largest displacements are generally seen for the S = 5/2 spin state, amounting to 0.47-0.50Å for OR − anionic ligands (including hydroxide).As expected, smaller displacements are seen for the ground states of weakly bound ligands such as O-coordinated oxazole, water and thiophene (0.13-0.15 Å).The Fe-N( porphine) bond lengths for the S = 5/2 structures are in the range of 2.052-2.115Å, longer than the range observed for the 3/2 spin state (1.977-2.036Å).Among the O donor ligands, the Fe-O bond lengths are 1.797-1.928Å for the anionic ligands and 1.883-2.359for the neutrals; the N donor ligands have Fe-N bond lengths of 1.857-2.280Å, with the anionic ligands tending to give shorter bonds.It is worth mentioning that the orientation of the acetate ligand in the porphine model complex is different to that observed in the X-ray structure of the FePPIX coordination dimer, 6 in that the positions of the non-coordinated O and C are transposed; however this output geometry was obtained regardless of the initial orientation of the acetate in the porphine model, indicating that it is preferred in the absence of other constraints.Further details of the calculated geometries are given in the ESI.]26 The method 2 and 3 binding energies of the exogenous ligands for the ground state structures in vacuum, water, and n-octanol are given in Table 2, and selected data are also presented graphically in Fig. 1.The calculated ground state struc-tures of selected species are shown in Fig. 2. Seven of these species (imidazole, AcO − , MeS − , MeNH 2 , H 2 O, OH − and Cl − ) were also modelled by Wondimagegn and Rauk, 14 using a somewhat different model chemistry; in addition, they included MeC 6 H 4 O − , whilst in the present study PhO − was used.On comparing data for these eight species between the two studies, plots of the binding energies in vacuo and the reaction solvation energies both gave excellent straight line correlations, with R 2 values of 0.995 and 0.996 respectively, although the absolute values of the in vacuo binding energies obtained with method 3 in the present study were systematically lower.
A number of trends are evident from the data in Table 2. Thus, the binding of all of the ligands is stronger in n-octanol than in water, particularly for the anions; the stabilization energies on transferring from water to octanol are 10.6 to 16.6 kcal mol −1 for the anions, compared to 0.7 to 6.5 kcal  a Values not in parentheses were obtained with method 3; values in parentheses were obtained with method 2; all values include ZPE corrections obtained with method 1. b Abbreviations; AcOH = acetic acid, AcO − = acetate, Me 2 NPy = 4-dimethylaminopyridine, MePy = 4-methylpyridine, PyO − = 2-hydroxypyridyl.c S = 3/2 ground state for method 3 calculations in octanol and water; otherwise S = 5/2 ground state.d S = 5/2 ground state for method 2 calculations in water and octanol, and method 3 calculation in water; otherwise S = 3/2 ground state.e S = 5/2 ground spin state for method 2 calculations in water and octanol; otherwise S = 3/2 ground state.f DMSO ligand dissociates for both the S = 1/2 and 5/2 spin states.
mol −1 for the neutral ligands.In general, the iron(III) centre prefers anionic over neutral ligands, although the two ranges of binding energies in water overlap.Iron(III) is considered to be a hard metal in Pearson's classification, and this is reflected in the binding energies of several ligand pairs (MeO − > MeS − ; F − > Cl − ; O-DMSO > S-DMSO).On the other hand, Me 3 P binds more strongly than Me 3 N, but this can be explained by steric considerations arising from the markedly different Fe-N and Fe-P bond lengths (2.280 and 2.651 Å respectively; and see below).Focusing on those ligands most relevant to haematin formation, hydroxide and methoxide are among the strongest binding, whilst acetate is among the weakest binding of the anionic ligands, although it still benefits from a notable stabilization of 12.7 kcal mol −1 on transferring from water to octanol.Comparing heterocyclic ligands, binding is stronger for more electron-rich species (imidazolide > imidazole, pyridine > triazine).Morpholine and oxazole are both preferentially coordinated through nitrogen rather than oxygen.Steric effects are also evident; these are most clearly seen for the weaker binding of quinoline compared to pyridine, where the porphine ring of the quinoline complex is noticeably warped (see Fig. 2).Experimental studies indicate that quinoline interacts with FePPIX by π-π stacking rather than coordination. 21Steric effects may also come into play for the 2-hydroxypyridyl anion, which preferentially binds through oxygen, and the two amines MeNH The pK a of the aqua complex is a point of particular interest.Wondimagegn and Rauk predicted a very low pK a value of  −13 for [Fe( porphine)(H 2 O)] + . 14Experimental determination of the pK a of [Fe(PPIX)(H 2 O)] is complicated by the dimerization of this species in solution, which is thought to involve π-π stacking interactions; nevertheless a pK a of 7.3 has been reported for the monomer, 27 and values of (6.2 and 7.0), and (8.5 and 8.06), have been reported for ionization of the first and second aqua ligands respectively in the [Fe(PPIX)(H 2 O)] 2 dimer. 27,28We have previously reported a method for estimation of pK a values for aqua complexes, using DFT calculations with PCM solvent corrections; 19a applying this method to [Fe( porphine)(H 2 O)] + gave a calculated deprotonation energy of +31.1 kcal mol −1 and pK a of 7.5, in excellent agreement with the experimental values for [Fe(PPIX)(H 2 O)].

Ferriprotoporphyrin IX complexes
The present study has been extended to consider the binding of OH − , Cl − , OAc − , MeNH 2 and H 2 O ligands to the iron(III) centre of FePPIX, using the same methodology as for the porphine models.As well as the possibility of various spin states, an additional complication in modelling these systems compared to iron porphines is the presence of two propionate side chains, with their associated protonation equilibria.Therefore, all three protonation states have been considered for both the S = 3/2 and 5/2 spin states, giving six different permutations for each ligand.In general, the ground spin state for each of these complexes was found to be the same as for the corresponding porphine model; the only exception was for [Fe(H 2 O)-(PPIXH)], for which the S = 5/2 spin state was preferred for the method 2 calculation in water.In this case, the higher spin state was 1.3 kcal mol −1 lower in energy, whereas the method 3 calculation predicted it to be 1.3 kcal mol −1 higher in energy.Hence, for such species, the difference in spin state energies is within the error margin of the calculations.For the vacuum calculations by methods 2 and 3, the Mulliken spin densities on iron covered ranges of 2.1 to 3.0 and 3.8 to 4.2 for the S = 3/2 and 5/2 spin states respectively; spin distributions were generally similar to the porphine models, except that the deprotonated carboxylate oxygen atoms provided an additional site for significant spin density in some cases (up to 0.40); these spins were diminished to near-zero when either solvent was included, with concomitant increase in the spin on iron.Checks for spin contamination revealed nothing untoward; rather, such cases suggest an internal redox reaction in which the negatively charged carboxylates are partially oxidized by the iron(III) when modelled in a vacuum.
The calculated geometries of the FePPIX models are very similar to those of the porphine complexes discussed above.The binding energies of the exogenous ligands for the FePPIX species are given in Table 3.For all three methods, there is very good agreement between the bond energy calculations for the neutral [Fe(PPIXH 2 )] species and the porphine models; this is reasonable since the charges on the two sets of models are the same for this series.In general, the binding energies of all the exogenous ligands decrease successively in the order: porphine species > [Fe(PPIXH 2 )] species > [Fe(PPIXH)] − species > [Fe(PPIX)] 2− species.As expected, for the vacuum calculations the binding energies of the anionic ligands are very sensitive to the overall charge on the FePPIX species, but the solvent corrections eliminate most of this variation.Binding of all the exogenous ligands to all three FePPIX species is stronger in n-octanol than in water, as was also observed for the porphine models.Overall, it may be concluded that the simpler porphine models retain the essential features of the larger FePPIX system, particularly the form for which both propionates are in the neutral -CO 2 H state. Considering the results for the H 2 O ligand, OPBE method 2 again appears to somewhat under-esti-Table 3 Calculated binding energies a and ground spin states (in water and n-octanol) for coordination of ligands to iron(III) ferriprotoporphyrin IX

Species
Ground spin state ΔE (vacuum)/kcal mol −1 ΔE (water)/kcal mol −1 ΔE (octanol)/kcal mol −1 mate the ligand binding energies.Given the experimental observation that water acts as a labile ligand in such species, 15 the small positive values obtained with B3LYP method 3 seem more plausible than the larger positive values obtained with OPBE.
A notable feature of the optimized FePPIX geometries is a hydrogen bond between the two propionate carboxylic groups, which is observed whenever either or both carboxylates are protonated.This is illustrated for the chloride complexes in Fig. 3. Calculations on [Fe(H 2 O)(PPIXH)] in which the geometry was re-optimized with the non-protonated propionate rotated away from the protonated propionate, gave method 3 energies that were 11.5, 6.2 and 7.4 kcal mol −1 higher in vacuum, water and n-octanol respectively.Such intramolecular hydrogen bonds are well known to modulate the pK a values of simple diacids.Therefore, in order to estimate the pK a 's of the propionate groups in these complexes, the aqueous protonation energies ΔE p of a series of hydrocarbon mono-and diacids were calculated using method 1 and the SMD implicit solvent model, according to eqn (1); The input geometries of the monoprotonated diacids were constructed so as to include the internal hydrogen bond.A plot of experimental pK a against calculated ΔE p values obtained with method 1 gave a straight line, with correlation coefficient 0.862 for 45 pK a values ranging from 2.83 to 6.74 (see ESI †).The pK a values of the FePPIX propionic acid groups were then obtained from their calculated first and second protonation energies by applying the equation of the straight line for the reference acids, and are reported in Table 4.The first pK a is essentially invariant at ca. 4.3, whilst the second is increased by ca.0.8 when an anionic axial ligand is present, such that the doubly deprotonated [Fe(X)(PPIX)] 2− species is a dianion.The pK a 's of the FePPIX propionate groups have proved difficult to measure experimentally, although it has been reported that the compound precipitates as β-haematin between pH 3-6, but is soluble outside of this pH range. 28hese observations seem consistent with the estimates in Table 4; this is the range in which the FePPIX should have one protonated and one deprotonated propionate, which is the form required for β-haematin crystallization.
Estimation of the pK a of the water ligand in the FePPIX aqua species using our published method 19a is hindered by the inability of the PCM solvent method to cope with tight hydrogen bonding geometries.Therefore, PCM calculations were only possible for the [Fe(H 2 O)(PPIXH 2 )] + and [Fe(H 2 O)-(PPIX)] − forms.These gave deprotonation energies of +32.4 and +35.1 kcal mol −1 , leading to pK a values of 7.7 and 8.3 respectively.These again are in good agreement with the experimental values discussed above.
Next, formation of the propionate-bonded FePPIX dimer was modelled, including final single point calculations using method 3. The structures of these models are shown in Fig. 4. On geometry optimization starting from the X-ray crystal structure, 6 the non-coordinated propionates realigned to form a hydrogen bond with the coordinated propionates [structure (a)]; this behaviour is similar to the monomeric systems discussed above.By re-optimizing the geometry in a slightly modified form to remove these hydrogen bonds [structure (b)], the final energies were increased by 6.8, 3.0 and 2.3 kcal mol −1 in vacuum, water and octanol respectively.Hence, the internal propionate hydrogen bonding remains marginally favourable even when the anionic propionate is coordinated to iron.In the structure of β-haematin, these intramolecular hydrogen bonds are replaced by intermolecular hydrogen bonds between the dimer units.The intermediate structure (c), containing only one Fe-O bond, optimized to a geometry in which the two porphyrin rings are essentially perpendicular.The successive method 3 energies for the first and second steps of the dimerization [corresponding to formation of structures (c) and (a) respectively in Fig. 4] were as follows; −50.1 then −33.9 kcal  , and indeed this species is not observed in aqueous solution, but rather the facesharing π-π dimer. 27In octanol, formation of the Fe-O bonded dimer is more favourable, at −17.0 kcal mol −1 , whilst π-π interactions are known to be attenuated in less polar solvents. 21ence, the lipid environment of the parasite's digestive vacuole serves to promote the Fe-O bonded form of FePPIX dimer required for crystallization of haematin, over the π-π form observed in aqueous solution.It is worth emphasizing that although the present calculations cannot give a reliable description of the π-π interactions in this system, these energies are reasonably well known from experiment, as summarized in the Introduction.
The effect of solvent polarity on the strength of the hydrogen bond between the propionic acid groups in β-haematin is an additional point of interest, since this is an important stabilizing interaction in the crystals.In order to probe this, a series of calculations on the acetic acid dimer were carried out.Using B3LYP with a large basis set, as described in the Computational methods, gave an in vacuo dimerization energy of −14.5 kcal mol −1 , in excellent agreement with the experimental heat of dimerization 29 of −14.4 kcal mol −1 .In water and octanol, the calculated dimerization energies were −6.9 and −7.6 kcal mol −1 respectively.Extending this approach to 3-(2,4-dimethyl-1H-pyrrol-3-yl)-propionic acid [Scheme 1, structure (4)] as a better model for FePPIX gave dimerization energies for this compound of −14.8, −6.9 and −7.3 kcal mol −1 in vacuum, water and octanol respectively.Hence, the lipid environment provided by the parasite for FePPIX crystallization should give only a marginal increase in the strength of the hydrogen bonds between the propionate groups, much less than the enhancement in Fe-O bond strength.

Conclusion
Notwithstanding some minor disagreements, particularly over ground spin states for some complexes, all three computational methods used in this study provide a common overall picture of the relative ligand binding preferences of the iron(III) centre in FePPIX.In particular, this site prefers anionic ligands, particularly those of weak conjugate acids such as alkoxide and hydroxide.Among neutral ligands, electron-rich, sterically undemanding donors such as primary amines, and phosphine or amine oxides, are most favoured.Binding of all ligands is more favourable in n-octanol than in water, and the difference is greater for anionic ligands than for neutrals.Both specific and more general comparisons with the available experimental data indicate that in terms of bond energy calculations, the B3LYP functional performs better than OPBE, and in particular the combination of B3LYP with the 6-311+G(d,p) basis set gives bond energies that appear to agree quite well with experiment.The results of this theoretical study tally well with the known chemistry of β-haematin formation.The ligand binding energy and pK a calculations indicate that in aqueous solution, at the usual physiological pH of 7.2, FePPIX should exist predominantly as [Fe(PPIX)] − , in which both propionates are deprotonated and the fifth coordination site is either unoccupied, or transiently occupied by water.Ionization of an aqua ligand to the much more strongly bound hydroxide requires a somewhat higher pH (ca.8).The unfavourable energy for coordination of an acetate ligand to the [Fe(PPIX)] − state in water (see Table 3) suggests that displacement of the aqua ligand to give the iron propionate dimer will be unfavourable.This agrees with experimental observations that the form of dimer under these conditions is not the iron propionate, but rather the π-π stacked dimer. 27Hence, the digestive vacuole of the parasite provides an acidic environment to aid haemozoin crystallization; the pK a estimates in Table 4 suggest that a pH of ca. 5 is optimal for the protonation of a single propionate to give the required [Fe(PPIXH)] species.This also moves the possibility of interference from a strongly bound hydroxide ligand firmly out of range, and makes carboxylate more competitive against water as a ligand.Propionate coordination is further enhanced by the lipid-rich environment of the vacuole, which again selectively promotes binding of carboxylate over water (cf.Fig. 1), as well as lowering the local concentration of water.The advantages provided by the acidic lipid environment for β-haematin formation can be summarised by considering reactions (2) and (3) below; Eqn (2), which pertains to the situation in aqueous solution at pH 7.2, has an overall reaction energy in water of +4.1 kcal mol −1 (see Table 3); whilst eqn (3), which is more appropriate to the parasite's digestive vacuole, has an energy in n-octanol of −13.7 kcal mol −1 .Inspection of the data in Table 3 indicates that this promotion of carboxylate binding to the iron centre arises more or less equally from the changes in protonation state of the FePPIX species (worth 8.4-10.0kcal mol −1 ) and the switch from water to octanol (worth 7.8-9.4kcal mol −1 ).Similarly, considering the Fe-O bonded [Fe(PPIXH)] dimer (Fig. 4), formation of this species is promoted by 13.6 kcal mol −1 on switching from water to octanol.Therefore, the digestive vacuole provides a suitable environment for the formation of the propionate-bridged FePPIX dimer that is the prerequisite for β-haematin crystallization.
Finally and most importantly, the results presented in this paper offer some insights for antimalarial drug development.The X-ray crystal structures of FePPIX complexes with QN, QD and halofantrine 7,8 reveal that these antimalarial drugs act as zwitterions, coordinating to the iron atom by alkoxide groups.As revealed by the data in Table 2, alkoxides are among the strongest binding ligands for the iron(III) centre of FePPIX, and compare very favourably with the carboxylate binding required to form β-haematin. Therefore, the cost of tautomerization to the zwitterionic form is likely to be more than offset by formation of the Fe-O bond for these drugs, particularly in a lipid-rich environment.With a view to developing novel drugs that can bind to the iron centre, none of the neutral ligands in Table 2 are very competitive against carboxylate, particularly in octanol solution (cf.Fig. 1).Furthermore, some of the anionic ligands that can bind more strongly than carboxylate have potential shortcomings as possible starting points for novel drugs.For example, imidazolide would be even more difficult to access in the acidic vacuole than the alkoxide groups of known drugs, whilst thiolates would probably be unstable with respect to oxidation.Phenoxide and O-coordinated 2-hydroxypyridyl are both competitive against carboxylate, but their planar geometries might prove difficult to elaborate into structures that can also provide secondary interactions with FePPIX, such as hydrogen bonds to the propionates and/or π-π stacking interactions.The most promising novel drug leads are probably dimethyl phosphinic acid and trimethylsilanol.Both are predicted to bind more strongly than acetate, and are also more acidic than simple alcohols, suggesting that zwitterion formation would be easier.The possibilities of these two functional groups in antimalarial drug development do not appear to have been investigated so far, and is a focus of our ongoing research.Finally, it is interesting to note that CQ is an effective inhibitor of haemozoin crystallization, even though it appears to lack any functional groups capable of binding strongly to the iron centre.The possible molecular basis of CQ's action will be discussed in a forthcoming paper.

Computational methods
All DFT calculations were performed with Gaussian09. 30Input geometries were created using Hyperchem 8.0.10, 31 and processed for Gaussian input with MolDraw 2.0. 32Figures of molecular structures were prepared using Ortep-3 for Windows. 33The standard procedure for DFT calculations on the iron complexes was as follows.For method 1, an initial gas phase calculation and wavefunction stability check using the B3LYP/LanL2DZ level of theory was followed by geometry optimization and frequency calculations, also at the B3LYP/ LanL2DZ level, and finally single point implicit solvent calculations at the optimized geometry using the same level of theory with SMD solvent corrections and additional wavefunction stability checks (method 1 results are given in the ESI †).
To assist convergence, quadratic convergence was routinely used for the initial step, but used only as required for subsequent steps.Two of the optimized porphine structures had persistent imaginary frequencies, which could not be eliminated by the usual methods.For the Me 3 PO complex in the S = 3/2 spin state, a persistent imaginary frequency at −6 cm −1 was associated with rotation about the P-O bond, whilst the Me 3 SiO − complex in the S = 5/2 spin state had a persistent imaginary frequency at −10 cm −1 , associated with rotation about the Si-O bond.Both of these were eliminated by 5°r otations about the relevant bonds, with negligible increases in energy (<0.1 kcal mol −1 ).Method 2 and 3 calculations consisted of single point calculations at the method 1 optimized geometries, using SMD solvent corrections and the OPBE/ 6-311+G(d,p) and B3LYP/6-311+G(d,p) levels of theory respectively, and including wavefunction stability checks.Ligand binding energies reported in Tables 1-3 have been corrected for zero-point energies (ZPE's), which were scaled by a factor of 0.981. 34Calculations on the FePPIX dimer structures were carried out for the high spin (S = 6/2) states; a trial calculation on the antiferromagnetically coupled S = 0 state made a negligible difference to the geometry and energies.
For the organic molecules, the B3LYP/6-311++G(3df,3pd) level of theory was employed.Structures were optimized separately in vacuo and using SMD solvent corrections for water and n-octanol, and ZPE corrections (scale factor 0.970 34 ) are included.Preliminary calculations on the acetic acid dimer using the OPBE/6-311++G(3df,3pd) level of theory gave a dimerization energy in vacuo of −8.7 kcal mol −1 , in much poorer agreement with the experimental heat of dimerization than the equivalent result with the B3LYP functional.
2 and Me 3 N.The former shows stronger binding, a shorter Fe-N bond (2.185 versus 2.280 Å), a greater Fe-N-C angle (117.5 versus a mean of 109.7°) and a smaller displacement of Fe out of the plane of the porphine nitrogens (0.19 versus 0.24 Å).

Fig. 4
Fig. 4 Calculated structures of the [Fe(PPIXH)] dimer.(a) Initial optimization of the X-ray structure results in the formation of intramolecular hydrogen bonds.(b) Modified structure, without internal hydrogen bonds.(c) Intermediate structure incorporating only one Fe-O bond.Atom colouring as in Fig. 3.The intramolecular hydrogen bonds are shown by dashed bonds; non-propionate hydrogen atoms are omitted for clarity.

Table 1
Comparison of experimentally obtained free energies (ΔG) a and calculated bond energies (ΔE) for coordination of N donor ligands in water Values calculated from the data given in ref.21, using the solvent corrected pK a values shown in the table.bValuescalculatedby applying a correction, calculated from the data in ref.23-25,to the aqueous pK a values given in ref.21. a

Table 2
Binding energies and ground spin states (in water and n-octanol) for coordination of exogenous ligands to iron(III) porphine, calculated by methods 2 and 3, a and ordered by binding energies in water b

Table 4
Calculated pK a values (in water) for the propionate groups in iron(III) ferriprotoporphyrin IX species