β-Hematin Crystal Formation: New Insights from Molecular Dynamics

Mar 15, 2016 - Understanding the driving forces observed during the first steps of β-hematin crystal formation is vital in view of developing new the...
0 downloads 14 Views 4MB Size
Download PDF
Article pubs.acs.org/crystal

β‑Hematin Crystal Formation: New Insights from Molecular Dynamics Simulations of Small Clusters in Condensed Phase Jean-Paul Becker, Fan Wang, Pascal Sonnet, and François-Yves Dupradeau* UMR CNRS 7378, UFR de Pharmacie, Université de Picardie Jules Verne, 1, rue des Louvels, F-80037 Amiens, France S Supporting Information *

ABSTRACT: Understanding the driving forces observed during the first steps of β-hematin crystal formation is vital in view of developing new therapeutic treatments against malaria. In this context a new additive Amber force field specific to the cyclic dimer of ferriprotoporphyrin IX coordinated via Fe−O propionate bonds (FPD) was developed and validated in the context of a crystal model. The structure and dynamics of small clusters or oligomers of the FPD entity were studied by molecular dynamics in condensed phase. New hydrogen bond patterns were identified in the dynamical process of the simulations: the characteristic pair of coplanar hydrogen bonds reported in crystallographic structures between FPD sheets (Pagola et al. Nature 2000, 404, 307) is found to be in a subtle equilibrium with new hydrogen bonds between FPD located within the same sheet. Hydrophobic interactions appeared to play a key role in oligomer cohesion: clusters with intersheet hydrogen bonds are found unstable compared to clusters with embedded-type molecular systems. We propose that oligomers enclosing this kind of embedded structure could serve as elementary building motifs for crystal growth in the digestive vacuole of Plasmodium, and could represent an attractive target for impairing hemozoin formation in a therapeutic approach.

■

consists of “head-to-tail” hematin dimers or cyclic dimers of ferriprotoporphyrin IX coordinated via Fe−O propionate bonds (FPD). In the crystal FPD entities are stacked allowing π−π interactions between the porphyrin rings of two successive FPD. This sequence of porphyrin rings alternatively associated by ferri-propionate coordinate bonds and π−π interactions forms arrangements comparable to parallel sheets. The relative orientation of FPD belonging to two neighboring sheets allows the formation of pairs of intersheet hydrogen bonds (HB) between free carboxylic acid groups. Hence, π−π interactions between the porphyrin rings, the HB between the carboxylic acid groups, and the ferri-propionate bonds are considered as the main interactions participating in the crystal cohesion.29−31 Although β-hematin and hemozoin structures have been described, the relative importance of the different interactions stabilizing these tridimensional structures is not fully understood, and the role of π−π interactions in the crystal stabilization is poorly described.31−34 Interactions known to stabilize experimental tridimensional structures were deduced from static observations obtained mainly from crystallographic methods. So far no information deals with the stability and the dynamics of these interactions during atomic and molecular motions in condensed phase over time. In order to identify the driving forces observed during the early steps of the β-hematin crystal formation, we report on the study of the structure and the dynamics of small FPD clusters or oligomers by molecular

INTRODUCTION Malaria remains one of the deadliest diseases: 214 million new cases, 438 000 deaths, and 3.2 billion people at risk were estimated in 2015. The death toll is particularly high in children under five years old in Sub-Saharan Africa.1,2 The pathogen of malaria is a blood-feeding parasite of the genus Plasmodium transmitted to human beings by the bite of the mosquito Anopheles. While the proportion of Plasmodium vivax cases still rises in certain regions of the world, the majority of malaria infections and deaths are due to Plasmodium falciparum. These species have developed resistance mechanisms against available drugs challenging even artemisinin derivatives in Southeast Asia and South America.3−5 The growing emergence of resistance to existing drugs coupled with the appearance of more virulent species of Plasmodium call for the development of new treatments.6−8 One common chemotherapeutic strategy against Plasmodium is to interfere with the hemozoin detoxification pathway. During the erythrocytic stage of its life cycle, Plasmodium digests the hemoglobin of its human host to obtain the nutrients essential for its growth.9,10 As a waste byproduct the hemoglobin digestion produces heme,11 which crystallizes into hemozoin in the food vacuole of the parasite.12,13 While hemozoin is innocuous,14 free heme can disrupt membranes, thus being lethal for the parasite.15−17 In a therapeutic approach several drugs active against malaria, especially chloroquine or quinoline and related 4-aminoquinoline, are known to interact with hemozoin impairing its formation.18−24 Comparative observations of hemozoin and of β-hematin synthetic crystals reveal important similarities in terms of geometry, shape, and proportions.18,25−29 β-Hematin © 2016 American Chemical Society

Received: January 12, 2016 Revised: March 15, 2016 Published: March 15, 2016 2249

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

Figure 1. RESP charge derivation, atom typing, FF library building, and FF parameter generation for FPD stereoisomers and analogs were carried out by using the PyRED program available at R.E.D. Server Development. (A) Description of the elementary building blocks (numbered from 1 to 7) involved in FF generation; blue circles: capping groups involved in intramolecular charge constraints applied during the charge fitting step. (B) Listing of the molecular fragments (numbered 1−7) issued from these building blocks; fragments 1, 3−5 (or 2−5) allowed the construction of FPD stereoisomers; fragments 6 and 7 within the gray box were involved in the construction of FPD analogs; covalently and ionically bound Fe(III) porphyrin analogs were differentiated (fragments 1 and 2, respectively); open valencies are underlined with blue arrows; FPD derivatives and oligomers were constructed by connecting these molecular fragments. empirical molecular fragments, which were issued from the monomeric heminic molecule, that was duplicated and from the organic molecules: after removing a methyl or an ethyl capping group, ethane, propene, pentanoic acid, pentanoate, and pentane were used as donors of methyl, vinyl, 2-carboxyethyl, 2-carboxylatoethyl, and propyl groups for the porphyrin ring of the heme derivative, respectively. The porphyrin ring was used as an acceptor of these molecular fragments removing its eight capping methyl groups and the ethyl group of the pentanoate ligand (Figure 1B). Quantum mechanical computations were achieved by using the Gaussian 2009 program (v D.01).39 For the heme building block the high spin state (S = 5/2) of Fe(III) was considered in density functional theory (DFT) calculations.40 The geometries of the different input molecules considered in this work were optimized by using the 6-31G(d) basis set and the B3LYP method.41−44 For each building block the lowest minimum found after a conformational search was selected. To account for the dimeric structure of FPD the geometry of the monomeric heme building block was optimized with constraints: the positions of the α-carbon atom of the pentanoate ligand and of the two closest heminic nitrogen atoms to this α-carbon atom were frozen during geometry optimization.45 These constraints were chosen because they did not lead to any imaginary frequency, and yielded an optimized structure, which was close to that observed in the diheminic crystal structure. The wave function of each minimum was checked: if found unstable it was optimized and the geometry of the corresponding minimum was optimized again using the newly determined wave function. Molecular electrostatic potentials (MEP) were also computed at the B3LYP/6-31G(d) theory level, and involved the Connolly surface algorithm with the default density of points and a radius of 1.8 Å for the iron atom.46 For each optimized geometry a pair of molecular orientations was automatically selected based on a rigid-body reorientation algorithm ensuring the control of the orientation of the geometry involved in MEP computation, and consequently the reproducibility of the atomic charge values.35,47 To

dynamics (MD) simulations in an aqueous medium. A new Amber-type additive force field (FF), named the Becker et al. FF, was developed for FPD and FPD clusters, and validated based on a 27 FPD crystal model. Oligomers were built using the β-hematin crystal structure as a template, and can be considered as crystal synthons. An ensemble of 15 representative dimers, trimers, and tetramers, which can be grouped into two types of oligomers was studied: (i) clusters with fully characterized molecular contacts between FPD entities observed within the crystal, and (ii) clusters, that are unlikely to be stable, but provide diverse input structures for MD. The concern here is more to enclose in this ensemble interactions and structures that are significant for the crystal formation rather than exhaustiveness.

■

COMPUTATIONAL SECTION

The Becker et al. FF was developed as follows: atomic charge derivation, atom typing, FF library building, and FF parameter generation were carried out by using the PyRED program available at the R.E.D. Server Development Internet site http://q4mdforcefieldtools.org/REDServer-Development/.35−37 This empirical work was achieved on the basis of six different input molecules or elementary building blocks (Figure 1A):38 a well-designed monomeric derivative of heme and a list of five organic molecules: ethane, propene, pentanoic acid, pentanoate, and pentane. The porphyrin ring of the monomeric heme derivative presents methyl capping groups in place of the characteristic vinyl and propionic acid groups of heme B, and contains a ferric iron, which is coordinated by the four heminic nitrogen atoms and a fifth pentanoate ligand. This heme derivative was repeated twice in the list of input molecules allowing the consideration of both a covalently bound ferric ion bearing a noninteger charge value and an ionically bound ferric ion with an integer charge value of +3 e. FPD stereoisomers and analogs of interest can be constructed from 2250

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

Figure 2. (A) Stick representation of the input structure of the FPD entity involved in MD: hydrogen, carbon, nitrogen, oxygen, and iron atoms are colored in white, cyan, blue, red, and white, respectively. For clarity, only the hydrogen atoms of the carboxylic acids are represented. (B) Representative structure extracted from a MD trajectory of FPD: two intra-FPD HB are formed, and are circled in black. (C) MEP map of FPD computed by using the Jmol program: the Connolly surface of FPD generated with a probe of radius 1.4 Å is colored according to the MEP values. The latter range from −0.19 to 0.17 au, while the color scheme spans from −0.10 to 0.10 au. The most negative regions appear in red, the most positive in blue, and the neutral ones in green. FPD conformation and orientation are the same as in A. parameters were combined into frcmod files.55 This project was submitted in the RESP ESP Charge DDataBase http://q4mdforcefieldtools.org/REDDB/,56 and is available for download under the “F-94” R.E.DD.B. code. The Becker et al. FF was exported within the LEaP program using a dedicated script to generate parameter/ topology and Cartesian coordinate files,57 and was validated by MD simulations in condensed phase (see the Discussion section for justification). The force field parameter file, force field libraries, and the leaprc file are also provided in the Supporting Information. MD simulations were performed using the pmemd module of the Amber 12 distribution.54,58 FPD molecular systems were solvated in a truncated octahedral box with a buffer distance of 10.0 Å. An aqueous condensed phase was chosen for MD, and the parameters used for water were taken from the TIP3P model.59 Molecular systems were equilibrated as follows: the procedure involved a minimization step followed by one canonical ensemble (NVT) and four isothermal− isobaric ensemble (NPT) MD periods of 100 ps each, where the temperature was ramped up to the target value. Positional restraints were applied to the FPD oligomers for relaxing the molecular systems during the equilibration procedure. Productive MD simulations were performed using the NPT ensemble at a pressure of 1 atm and a temperature of 300 K. The weak coupling algorithm was used to regulate the temperature and the pressure.60 The temperature was maintained close to the intended value by weak coupling to an external bath with a relaxation time of 2 ps, and the pressure to an external bath of 1 atm with a coupling constant of 2 ps. The Shake algorithm was used to constrain C−H bonds, and a time step of 2 fs was used to integrate the equations of motion.61 Periodic boundary conditions were imposed during simulation. A distance cutoff of 9.0 Å was applied to nonbonded interactions and the particle mesh Ewald method was used to compute long-range interactions.62 Configurations of the systems were stored at intervals of 1 ps. Trajectory outputs were saved in the NetCDF binary file format developed by Unidata (http://www. unidata.ucar.edu/software/netcdf/). To ensure reliable results MD simulations involved up to four different input conditions and between 100 and 500 ns for each cluster depending on the size of the molecular system. Moreover, the stability of each cluster characterized at 300 K was systematically checked in another series of MD at 325 K.63 MD analyses as well as structure and MEP visualization were carried out using the ptraj/cpptraj (Amber 12 distribution), VMD, and Jmol programs.64−66 Free energy values were calculated from the MMGBSA approach.67,68 The AYILIE crystallographic structure was used as reference in this work.30,52 The stability of the different clusters was mainly followed in term of root-mean-square deviation (RMSD), radius of gyration (RG) calculations, and presence of HB during MD.

further check each optimized geometry frequencies were also computed at the B3LYP/6-31G(d) theory level. RESP atomic charge fitting was carried out using a standalone version of the RESP program, and following the standard two-stage procedure.48 The molecular fragments required for the construction of the FPD stereoisomers and analogs were designed by setting specific intramolecular charge constraints during the charge fitting step for the capping groups. For each elementary building block a single intramolecular charge constraint with a value of zero was applied to the alkyl capping groups. More particularly the atomic charges of the eight methyl capping groups of the porphyrin ring and of the ethyl capping group of the ethanoate ligand were constrained to have a total charge of zero.49 The atoms involved in these intramolecular charge constraints were removed from the molecule to yield the required molecular fragments and the corresponding FF libraries to the Mol3 file format.50 This approach was found to be the one that led to the lowest relative root-mean-square (RRMS) value and to the highest correlation coefficient (r2) value for the charge fitting step. An additional intramolecular charge constraint of +3.0 e was used to define the ionically bound iron atom for a heme derivative. RRMS and r2 values of 0.057 and 0.996 between the MEP calculated by quantum chemistry and that generated using the derived charge values were obtained for the charge fitting step, respectively. Highly similar values were obtained in the absence of intramolecular charge constraints. The small RRMS value and high r2 value obtained for this two-stage charge fitting step as well as the small differences of RRMS and r2 values between the charge fitting steps carried out with and without intramolecular charge constraints demonstrate the accuracy of the fitting step performed in this study, and the weak effect of these constraints.48 Atom types and FF parameters for the input molecules were determined by the PyRED program based on an internal dictionary of atom types and a database of FF parameters corresponding to the Amber99 FF.51 Missing FF parameters for the organic parts of the input molecules were determined by analogy to existing ones. FF parameters involving the covalently bound ferric ion were determined as follows: (i) bond and angle equilibrium values were taken from the crystal structure (Cambridge structural database code: AYILIE),30,52 (ii) bond and angle FF constants were calculated from the Hessian matrix (using the Gaussian formchk file obtained from frequency computation) and by applying the Seminario’s method,53 (iii) the energy barrier of the dihedral FF parameters were set to zero, (iv) improper FF parameters were not considered, and (v) van der Waals parameters for the Fe(III) atom were taken from the Amber12 program contribution (R* = 1.2 Å, ε = 0.05 kcal/mol).54 FF 2251

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

Figure 3. (A) Stick representation of the input structure of the 27 FPD crystal model used in MD: the color scheme is the same as in Figure 2A. For clarity, the hydrogen atoms are omitted. (B) Stick representation of the input structure of the crystal model viewed along the arrow direction of Figure A. (C) Representative structure extracted from a MD trajectory of the crystal model: only two FPD are represented: the central entity among the 27 ones is in blue, and a neighboring entity is in orange. Two inter-FPD HB are formed between these FPD entities, and are circled in black. The oxygen and nitrogen atoms participating in these inter-FPD HB are represented by red and white spheres, respectively. RMSD and RG value calculations for the different clusters only involve heavy atoms, and were carried out versus the input structure, that was also the crystallographic structure. HB criteria defining its occurrence were a distance of 3.3 Å and an angle of 120° between the donor and the acceptor. Computations were carried out using cluster nodes running CentOS 6.4 and featuring 20 Intel Xeon cores (CPU E5− 2690 v 2 at 3.0 GHz) and 256 GB RAM.

HB induces the concomitant observation of gauche+ populations for the C−C−C−C dihedral angle of the carboxylate-iron bond and for the C−C−C−O dihedral angle of the carboxylic acid group in one heminic region, and of gauche− populations for the two reported dihedral angles in the second heminic counterpart. On the contrary, the geometry of the iron chelation site, including distances and angles between atoms, is in good agreement with the crystallographic reference. The iron atom remains outside the plane defined by the four nitrogen atoms in the direction of the carboxylate group as enlightened by the selected N−N−Fe−N improper dihedral angle, whose values stay close to the one measured in the crystallographic reference (Figure S6). Crystal Model Composed of 27 FPD Entities. The structure and dynamics of a 27 FPD crystal model constituted by a central FPD entity surrounded by 26 others (Figure 3A,B) was studied by MD allowing comparison with results obtained for a single FPD entity (Figure 2A,B). Simulations involving this crystal model were also used to validate the Becker et al. FF developed in this work (see the Discussion section for justification). RMSD values of the whole crystal model fluctuate around 1.3 Å with an amplitude of 0.4 Å (Figure S7). This indicates good overall agreement between the geometry of the FPD entities composing the crystal model and the crystallographic reference. The structure and dynamics of the central FPD entity within the crystal model were further studied. By contrast with simulations involving a single FPD the characteristic C−C−C−C dihedral angle of each propionate bond exhibits an anti conformation with values around ±180° (Figure S8) as in the crystallographic structure. In the meantime, the observed Fe−Fe distance values are slightly higher than 9.0 Å (Figure S9), in agreement with the one observed in the reference structure, and the geometry of the iron binding site is also well preserved (Figure S10a-c). Thus, these results, together with the RMSD values, which stay under 0.7 Å (Figure S1), indicate that the geometry of the central FPD entity is closer to the crystallographic reference compared to the single FPD molecular system. Meanwhile a noticeable difference appears in the HB network of the crystal model during MD. In the static crystallographic structure, the extended conformation of each

■

RESULTS FPD Entity. The structure and dynamics of the FPD entity were studied by MD. RMSD values oscillate around 1.6 Å with a magnitude of 0.4 Å (Supporting Information Figure S1). This reveals the evolution of the FPD structure (Figure 2A) toward a stable conformation slightly different from the one depicted in the crystallographic reference. Key features observed along the simulations include the torsion of the C−C−C−C dihedral angle of the carboxylate-iron bonds connecting the two heminic rings, and the decrease of the distance between the two ferric cations. In the crystallographic reference the values for this dihedral angle are found around ±180° corresponding to an anti population, and the Fe−Fe distance is measured at 8.9 Å.69 During MD the C−C−C−C dihedral angle shows gauche populations with values fluctuating around ±80°, and the Fe− Fe distance oscillates between 6.0 and 7.5 Å (Figures S2 and S3).69 A distinctive aspect in MD for this FPD entity is the emergence of a new type of interaction. Indeed in the crystallographic structure the two free carboxylic acid groups present a conformation allowing interactions between two sheets. In contrast, in the dynamic process of MD these groups undergo a conformational change resulting in the formation of new intra-FPD HB with the carbonyl oxygen atoms of the carboxylate-iron bonds belonging to the same heminic ring (Figure 2B). The evidence of the formation of these two HB is brought by the distances observed between the associated oxygen atoms, which stay close to 3.5 Å as well as the C−C− C−O dihedral angle values of the corresponding carboxylic acid groups, which fluctuate around gauche values (±100°; Figures S4 and S5). This is further sustained by the FPD MEP map reported in Figure 2C, where the most negative and positive MEP values (−0.19 and 0.17 au, respectively) are found for the atoms involved in these intra-FPD HB. The formation of these 2252

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

Figure 4. Stick representation of the 15 studied clusters of FPD molecular systems: the color scheme is the same as in Figure 2A. For clarity, the hydrogen atoms are omitted. The acronym of each cluster appears under its representation. Dimers, trimers, and tetramers appear on first, second, and third rows, respectively.

free carboxylic acid moiety allows the formation of the characteristic pair of coplanar HB between two sheets of FPD. In the crystal model two different carboxylic acid moieties can be found: a first one buried in the structure between two FPD sheets, and a second one exposed to the solvent. During the dynamical process of MD two-thirds of the pairs of intersheet HB are broken, leading to a new HB pattern with a significant lifetime (Figure S11): this HB, called inter-FPD HB in the remainder of the manuscript, is formed between a hydroxyl donor of the carboxylic acid group in a first FPD entity and a carbonyl oxygen acceptor of a carboxylate group in a second FPD located in the crystal direction [001] or [001̅] (Figure 3C). The formation of this HB is further evidenced by the distances observed between the associated oxygen atoms, which stay below 3.5 Å, as well as the C−C−C−O dihedral angle values of the corresponding carboxylic acid group, which fluctuate around gauche values (±100°; Figures S12 and S13). Moreover, among the 18 carboxylic acid moieties, which are exposed to the solvent within the crystal model, 12 can form a new HB with another FPD as just defined: thus, these carboxylic acid moieties prefer to form such an inter-FPD HB, rather than to be fully exposed to the solvent. Indeed occupancy rate values above 80% are found for these interFPD HB. All these observations are in agreement with the MEP map reported in Figure 2C, where the most negative and positive MEP values are also observed for the atoms involved in these inter-FPD and intersheet HB. Definition of the Studied FPD Oligomers. A list of 15 FPD clusters or crystal synthons made of two, three, and four FPD entities were studied by MD. These oligomers were defined as follows. Two different types of molecular systems were considered: clusters with representative interactions observed within the crystal, and clusters, which are unlikely to be stable, and which provide diverse input structures for MD. The dimers are constituted by a reference FPD entity and a neighbor entity involved in a representative assemblage with specific interactions in the crystal. Thus, an ensemble of five dimers named D1 to D5 was set up. In D1 the second FPD is localized in the crystal direction [001] to the reference. This results in two embedded FPD entities, whose large surface of

contact favors interactions between the porphyrin rings and between the vinyl and methyl groups. In D2 the second FPD is placed in the crystal direction [011̅] with respect to the reference. This corresponds to the situation where two FPD entities are superposed and allows interactions between the porphyrin rings and the vinyl and methyl groups. In D3 the second FPD is translated in the crystal direction [100]. In D4 the second FPD is located in the crystal direction [01̅2]. These two molecular systems permit weak interactions between the methyl and vinyl groups. In D5 the second FPD is situated in the crystal direction [101̅]. The characteristic interaction of this dimer is the pair of coplanar HB between the carboxylic acid groups. Trimers, T1 to T5, as well as tetramers, Q1 to Q5, were also considered, and were constructed based on dimers. In T1 a second FPD is positioned in the crystal direction [001], and a third is placed in the direction [001̅] with respect to the reference. This cluster corresponds to the dimer D1 with one additional embedded FPD. In T2 a second FPD is positioned in the crystal direction [001], and a third is placed in the direction [01̅1]. In T3 a second FPD is positioned in the crystal direction [001̅], and a third is placed in the crystal direction [01̅1]. In T4 a second FPD is positioned in the crystal direction [001], and a third is placed in the crystal direction [101]̅ . In T5 a second FPD is positioned in the crystal direction [101]̅ , and a third is placed in the crystal direction [1̅01]. This results in three FPD entities, which interact via two pairs of coplanar HB found between sheets. In Q1 three FPD entities are positioned according to the crystal directions [001], [002], and [003] with respect to the reference. This results in a sequence of four embedded FPD entities. In Q2 three FPD entities are positioned according to the crystal directions [001], [010], and [011]̅ with respect to the reference. In Q3 three FPD entities are positioned according to the crystal directions [001], [101̅], and [100]. In Q4 three FPD entities are positioned according to the crystal direction [001], [100], and [101]. In Q5 three FPD entities are positioned according to the crystal directions [101̅], [202̅], and [303̅]. This results in four FPD entities, which interact via three pairs of coplanar HB found between sheets. This whole list of clusters is described in Figure 2253

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

carboxylate-iron bond are also consistent with a structure close to the experimental one (Figure S31 and S32). The most noticeable event is the formation inter-FPD HB as already described for D1 with an occupancy of at least 60% (Figure S33). As in D1 the structure of dimer T1 is stabilized by interactions between porphyrin rings in an embedded molecular system (Figure 4). T2 RMSD values oscillate between 2.0 and 8.0 Å demonstrating a structure reorganization (Figure S34). More specifically the RMSD values computed for the substructure constituted by two FPD entities associated as in D1 remain close to 1.0 Å, while the ones computed for the substructure corresponding to D2 fluctuate between 2.0 and 8.0 Å (Figures S35 and S36). The agreement of the Fe−Fe distance values with the reference structure of D1 as well as the extended conformation of the C−C−C−C dihedral angle of each carboxylate-iron bond are observations in favor of the stability of this D1 substructure in T2 (data not shown). The inter-FPD HB similar to those described in D1 participate in this stability (data not shown). The family of substructure D2 in T2 displays similar structural and dynamical features to dimer D2 with a decrease of the solvent accessible surface area, and a structure stabilization by interactions between porphyrin rings as previously exposed (Figure S37). In the early stages of MD T3 RMSD values rise up to 4.0 Å demonstrating the instability of T3 (Figure S38). The RMSD computed for each pair of FPD is consistent with the formation of the D2 family of substructures. The third FPD moves and interacts with D2 to limit solvent accessibility (Figure S39). None of the structures generated during MD are in agreement with the crystal arrangement. T4 RMSD values increase up to values around 3.0 Å demonstrating a structure reorganization (Figure S40). The RMSD computed for each pair of FPD shows the preservation of the D1 substructure found in the input structure. The third FPD moves and interacts with D1 to limit solvent accessibility in agreement with the crystal structure (Figure S41). It is important to underline that the pair of HB found between sheets and observed in T4 is rapidly broken (Figure S42), and is overtaken by interactions involving porphyrin rings. T5 RMSD values as well as these computed on each pair of FPD demonstrate a structure reorganization. The two pairs of HB found between sheets and observed in T5 (Figure 4) are quickly broken (Figure S43), and are once again overtaken by interactions involving porphyrin rings. Representative structures of T5 involve D1 and D2 substructures depending on the considered simulation (Figure S44). Tetramers of FPD Entities. Tetramers Q1 and Q2 (Figure 4) are constituted by D1, and D1 and D2 substructures, respectively. RMSD values of Q1 and Q2 indicate that these clusters are stable (Figures S45 and S46): each input structure is conserved during MD. Q1 and Q2 are stabilized by interactions between porphyrin rings as shown in D1 and D2, and by pairs of inter-FPD HB as observed in D1, T1, and in the crystal model composed of 27 FPD entities. Q3 and Q4 (Figure 4) are composed of two D1 dimers connected by one and two pairs of intersheet HB, respectively. A brusque increase in the Q3 and Q4 RMSD values is observed after 60 ns indicating a reorganization of the Q3 and Q4 structures after a significant period of time, independently of the considered simulation (Figure S47). RMSD values computed for pairs of FPD reveal the stability of substructure D1 over the simulation. Moreover, study of HB patterns shows

4, and the Cartesian coordinates of each MD input cluster are provided in the Supporting Information. Dimers of FPD Entities. RMSD values of dimer D1 fluctuate between 0.6 and 1.6 Å (Figure S14), while RG values stay close to 7.5 Å (Figure S15), the value computed for the experimental structure. The distances between the two iron atoms in each FPD remain in good agreement with the reference (Figure S16). These results indicate that the two FPD entities are maintained in the embedded position (Figure 4). It is noteworthy to mention that the flexibility of the C−C−C−O dihedral angle of the carboxylic acid moieties of D1 permits the formation of inter-FPD HB as previously observed on the crystal model: among the four carboxylic acid groups only two have a partner to form this HB. Such HB are stable during MD, and have occupancy rates above 60% (Figure S17). Moreover, no solvent molecule is observed between the two embedded FPD entities: D1 remains surrounded by solvent molecules as demonstrated by the small variations of the RMSD and RG values. Thus, the structure of D1 is mainly stabilized by interactions between porphyrin rings and by inter-FPD HB. D2 RMSD values oscillate between 3.0 and 8.5 Å (Figure S18). RG values are found between 7.0 and 8.5 Å (Figure S19), while the experimental value is 10.4 Å. The distance between the two closest iron atoms of the two FPD entities decreases from 7.8 Å to values around 5.0 Å (Figure S20). These observations are manifestations of conformational changes of FPD resulting in molecular systems with a more important surface of contact between porphyrin rings. This phenomenon results in a decrease of the solvent accessible surface area, and is stabilized by interactions between porphyrin rings (Figure S21): this process is rapid, and achieved in less than 100 ps. The important oscillations of the RMSD values are related to complex movements of translation and rotation of the porphyrin rings between each others, and yield an ensemble of relatively heterogeneous conformations. The evolution of the Fe−Fe−Fe−Fe dihedral angle values characterizes these movements (Figure S22). Thus, the family of D2 structures does not lead to a stable and well-defined molecular arrangement found in the crystal. D3 and D4 RMSD values quickly reach relatively high values just after the equilibration period (Figures S23 and S24) revealing the reorganization of the D3 and D4 structures. As observed for cluster D2, the D3 and D4 oligomers reorganize to reduce solvent accessible surface area (data not shown). Representative structures of D3 and D4 strongly differ from the crystal structure demonstrating the instability of the two FPD entities in these two input molecular systems (Figures S25 and S26). The pair of HB observed in D5 are exposed to the solvent (Figure 4), and quickly break during the equilibration period. This is clearly demonstrated by the corresponding RMSD values, which increase up to 5.0 Å (Figure S27). RG values also decrease toward 7.0 Å (Figure S28), while the experimental value is 9.9 Å. The structure of D5 reorganizes into more compact conformations, limiting the accessibility to solvent, and results from the rupture of the HB found between sheets. Representative structures of D5 are in agreement with the D1 or D2 cluster depending on the considered simulation (Figure S29). Trimers of FPD Entities. RMSD values of trimer T1 fluctuate around 1.2 Å in agreement with a stable structure (Figure S30). The values of the characteristic Fe−Fe distances as well as these of the C−C−C−C dihedral angle of each 2254

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

Table 1. Free Energy Values of the Q1−Q4 Tetramers Computed Using the MM-GBSA Approach (in kcal/mol)a Q1 Average TS EGB GGb

420.42 10.68 −409.74

Q2 SD

Average 0.35 3.75 3.72

Q3 SD

424.86 10.99 −413.87

Average 0.07 1.61 1.55

420.37 33.09 −387.28

Q4 SD

Average 0.42 2.37 1.99

420.52 34.35 −386.17

SD 0.15 1.03 0.97

a

TS: entropic contribution at 298.15 K computed by normal-mode analysis. EGB: energy component computed using the Generalized Born model. GGB: free energy of the considered tetramers. For each tetramer four different MM-GBSA computations were performed on four independent trajectories sampling 4000 frames over a period of 40 ns after equilibration of the system. The values reported in this table are averaged over four MM-GBSA computations and presented with the associated standard deviation (SD).

presence of carboxylic acid groups accounts for mimicking the protonation state of FPD entities in the range of the pH values required for intersheet HB formation. Force Field Development and Validation. It is important to acknowledge that the FPD entity was experimentally characterized only in the β-hematin crystal: no evidence of an elementary FPD building block was reported during crystal growth.32 On the contrary different theoretical works involving DFT methods have considered FPD as a preformed entity.30,31,77 Indeed defining a topology for a bioinorganic complex with coordinate bonds between the metal center and its different ligands is a strategy widely adopted in structural studies involving MD simulations. Examples of such an approach are legion in the literature.78−80 This empirical strategy is the one we have chosen in this work. We have designed the Becker et al. FF: it is constructed using the building block approach,38 and is compatible with the known stereoisomers of FPD and with several analogs of interest.31 This FF involves the RESP atomic charge model, which is recognized to correctly handle intermolecular properties.35,38,46 These are essential in condensed phase simulations, where the solute−solvent and solvent−solvent interactions have to be well balanced. Indeed, nonbonded interactions play a key role in FPD cluster cohesion, and are modeled by electrostatic and van der Waals potentials in the additive Amber FF model.51 The topology of the pentacoordinated iron III complex in FPD is accounted by the definition of coordinate bonds between the iron atom and its different ligand atoms, and uses a noninteger partial charge value for the iron cation. MEP values computed from RESP charges present four main hot spots related to the chemical groups involved in intra-FPD and inter-FPD as well as in intersheet hydrogen bonding. Other regions of FPD present smoother MEP values, which is consistent with the establishment of interactions between the porphyrin rings belonging to different FPD entities ensuring crystal synthon cohesion. Considering the goals of this study, i.e., studying the driving forces involved in β-hematin crystal formation as well as the structure and dynamics of FPD oligomers, we have decided to validate the Becker et al. FF based on a crystal environment rather than on a single FPD entity. This choice was reinforced by the lack of experimental data for a single FPD entity and for the studied FPD oligomers in particular in a dynamical view. The ensemble of 3*3*3 FPD entities within the crystal model represents in our opinion a mimic of the crystallographic environment of FPD, and is the minimal oligomeric structure required for FF validation. Indeed in this crystal model the central FPD entity is surrounded by 26 others imitating the crystal environment and packing. Structural features obtained from MD demonstrate the good agreement between the central

a rupture of intersheet HB after 60 ns. This phenomenon is mediated by the solvent accessibility of the oxygen atoms, which are involved in intersheet HB: during MD intersheet HB are replaced by HB between the carboxylic acid groups and solvent molecules and by inter-FPD HB (Figure S48a−d). The two D1 substructures move toward FPD entity arrangements limiting the solvent accessibility of the global structure. Simulations lead to the Q1 type structure (data not shown). Q5 (Figure 4) is composed by four FPD entities connected by three pairs of HB found between sheets. The study of Q5 RMSD values reveals that the structure of Q5 breaks down in the early stages of MD. This is further demonstrated by the study of HB patterns: the three pairs of HB quickly break, and are replaced by HB between carboxylic acid groups and solvent molecules and by inter-FPD HB. FPD entities rearrange into more compact FPD molecular systems. Representative structures of Q5 involve D1 or D2 substructures depending on the considered simulation (Figure S49). MM-GBSA Study. Free energy computations for selected clusters were carried out using the AMBER MM-GBSA approach for tetramers Q1−Q4 (Table 1), which exhibit long stability in condensed phase allowing the collection of a significant number of representative snapshots. Independently of the considered simulation these tetramers can be separated into two groups according to free energy values. Tetramers Q1 and Q2 display the lowest values, while tetramers Q3 and Q4 have the highest ones: free energy differences between the two groups are greater than 20 kcal/mol. The entropic contribution is similar for each tetramer: 420 kcal/mol for Q1, Q3, and Q4, and slightly higher for Q2. Thus, most of the free energy differences between tetramers originate from the EGB energy component. The relatively higher free energy values of Q3 and Q4 compared to Q1 and Q2 are in agreement with weaker stability observed during MD.

■

DISCUSSION Condensed Phase Simulations. β-Hematin crystal formation is carried out in vitro in polar protic environments such as water−alcohol mixtures, while hemozoin synthesis occurs in vivo at the interface of water−lipid systems in the Plasmodium vacuole. Moreover, a pH around five is required for the formation of characteristic pairs of intersheet HB.70−76 In this work we have studied in silico the structure and the dynamics of small FPD clusters by MD simulations in condensed phase using the TIP3P explicit solvent model: to our knowledge this work represents the first attempt reported in the literature to evaluate the driving forces observed during the first steps of the β-hematin crystal growth using MD. In this context this solvent model leads to new insights obtained in condensed phase, is intended to model experimental conditions, and appears as an opening representation. The 2255

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

result observed during MD is the equilibrium observed between dimers D5 and D1, which shifts toward D1. D5 with its pair of coplanar HB observed between sheets is found unstable. This type of intersheet HB is overtaken by hydrophobic interactions with the emergence of new inter-FPD HB, and the decrease of the solvent accessible surface area. MD simulations as well as MM-GBSA computations related to tetramers Q3 and Q4 also show the relative instability of the clusters enclosing one and two pairs of coplanar HB. Indeed when exposed to solvent these HB break leading to new molecular organizations with a minimized solvent accessible surface area. In view of simulations in aqueous media this phenomenon, and more particularly the D5 to D1, T5 to T1, and Q5 to Q1 equilibria, could recall the hydrophobic collapse observed in proteins.82−84 A second series of results gleaned from the structural and dynamical study of the FPD oligomers is the following: (i) embedded structures D1, T1, and Q1 are in agreement with the experimental references, and are found to be among the most stable structures from MD and MM-GBSA studies; (ii) D2 structures and substructures in trimers and tetramers correspond to a heterogeneous family of structures, are found roughly as stable as the D1 structure, but do not correspond to a motif found in the crystal; (iii) D5, T5, and Q5 structures, which present HB patterns as those observed between FPD sheets, are rapidly broken, and are overtaken by interactions between porphyrin rings and between methyl and vinyl groups; and (iv) other clusters are not stable and lead to the D1 structure and substructure, or to nonrepresentative structures, i.e., molecular systems not observed within the crystal. These results demonstrate that these stabilizing interactions, classified here under the generic “hydrophobic interaction” term, are essential and ensure the cohesion of the clusters. In the literature, different hypotheses have been proposed about βhematin and hemozoin crystal formation. Thus, Strasso et al. have proposed that the β-hematin crystal growth is envisaged through the assembly of FPD dimers via hydrogen bonding, which subsequently form larger structures via π−π interactions.30,31 Klonis et al. have described a different mechanism, where the initial interactions responsible for crystal growth involve π−π interactions of monomeric heminic structure.32 Buller et al. have described that the fastest face of crystal growth is that corresponding to the [001] crystal face, i.e., to the embedded molecular systems reported in this study.18 From the results we have obtained here it is unlikely that the assembly of FPD dimers only hold via the pairs of coplanar HB found between HB sheets. Indeed the lifetime of the D5, T5, and Q5 clusters is limited, and the corresponding structures reorganize at the beginning of productive MD. Thus, the presence of two HB is clearly not sufficient to stabilize the dimeric FPD structure in a dynamic process. On the contrary, D1, T1, and Q1 clusters are found highly stable. Considering that the Becker et al. FF was designed for a preassembled FPD entity, we propose that the clusters presenting embedded-type structures could serve as elementary building motifs for crystal growth.18 However, extreme care should be taken considering the different media used to crystallize β-hematin in vitro, to form hemozoin in vivo, and the medium involved in this in silico study. In this context we have also studied the structure and stability of the D1 and D5 FPD dimers in explicit organic solvents such as chloroform and methanol by MD simulations: preliminary results obtained for D1 and D5 in these organic solvents confirm those reported using TIP3P, i.e., the great stability of the D1 cluster and the instability of the D5 one.

FPD entity within the crystal model structure and the crystallographic reference validating the Becker et al. FF. Hydrophobic Interactions in the β-Hematin Crystal. Different authors report on the importance of π−π interactions between porphyrin rings in the crystal stabilization.25,26,29−31 Debates have also appeared about the criteria defining these π−π interactions in β-hematin.34 Indeed from the measures of the distances between the porphyrin rings of two FPD entities π−π criteria are not fully satisfied, resulting in weak π−π interactions.32,33 Authors also describe that pronounced π−π overlap yields a low density crystal.31 On the contrary and in the context of an aqueous medium, π−π interactions between porphyrin rings as well as weak interactions between vinyl and methyl groups can be classified as hydrophobic interactions resulting from nonbonded interactions. From the MD simulations carried out here in condensed phase these hydrophobic interactions clearly play a major role in the cohesion of the studied clusters. Dynamical Features of the FPD Entity. MD involving a single FPD entity and MD performed on the crystal model made of 27 FPD entities allow comparing the conformation of this single FPD entity and that of the central FPD entity within the crystal model. Striking conformational differences are observed for these two types of FPD. When fully solvated a FPD entity undergoes key conformational reorganization: the two characteristic porphyrin rings come closer together. This is evidenced by the torsion of the C−C−C−C dihedral angle of the carboxylate-iron bonds connecting the two heminic rings (transition between anti to gauche populations), and the decrease of the distance between the two ferric cations. This conformational change can be explained by the absence of surrounding FPD entities and crystal packing forces, which maintain the characteristic structure observed for a FPD buried within the crystal. Dynamical Features of the FPD Oligomers. A first series of results obtained in this study is related to the HB present between FPD sheets in the crystallographic reference. The reported crystal structures display a single type of HB, i.e., intersheet HB.25,29,30 By contrast, three equilibria between different HB patterns were characterized during MD simulations: (i) intra-FPD HB within the single FPD entity versus HB between the FPD carboxylic acid groups and solvent molecules, (ii) inter-FPD HB versus HB between the FPD carboxylic acid groups and solvent molecules in oligomers enclosing a D1 substructure, and (iii) inter-FPD HB versus intersheet HB in the crystal model composed of 27 FPD entities. First it is interesting to note that the instability of the HB between the carboxylic acid groups and bulk water tends to demonstrate the relative hydrophobicity of the carboxylic acid groups in an aqueous medium. In addition such an intra-FPD HB has already been reported in a computational study involving the FPD entity, and was observed for a hexacoordinated gallium complex.77,81 Furthermore, the equilibria between the different HB reported above highlights the importance of MD in the presence of an explicit solvent to account for the dynamical feature of these HB compared to the static view found in a crystallographic structure. In this context it is relevant to underline that the authors of the AYILIE crystallographic structure acknowledge that during the refinement procedure the conformation of the carboxylic acid moiety was manually adjusted so that the pair of intersheet HB is formed.30 This remark clearly reinforces our observation of a dynamical process between different HB patterns. Finally a key 2256

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

■

Article

(3) Worldwide Antimalarial Resistance Network: Malaria Drug Resistance: http://www.wwarn.org/about-us/malaria-drug-resistance (accessed December 16, 2015). (4) Dondorp, A. M.; Nosten, F.; Yi, P.; Das, D.; Phyo, A. P.; Tarning, J.; Lwin, K. M.; Ariey, F.; Hanpithakpong, W.; Lee, S. J.; Ringwald, P.; Silamut, K.; Imwong, M.; Chotivanich, K.; Lim, P.; Herdman, T.; An, S. S.; Yeung, S.; Singhasivanon, P.; Day, N. P.; Lindegardh, N.; Socheat, D.; White, N. J. N. Engl. J. Med. 2009, 361, 455−46. (5) Pribluda, V. S.; Evans, L., 3rd; Barillas, E.; Marmion, J.; Lukulay, P.; Chang, J. Malar. J. 2014, 13, 77. (6) Eastman, R. T.; Fidock, D. A. Nat. Rev. Microbiol. 2009, 7, 864− 874. (7) Sanchez, C. P.; Dave, A.; Stein, W. D.; Lanzer, M. Int. J. Parasitol. 2010, 40, 1109−1118. (8) Phyo, A. P.; Nkhoma, S.; Stepniewska, K.; Ashley, E. A.; Nair, S.; McGready, R.; ler Moo, C.; Al-Saai, S.; Dondorp, A. M.; Lwin, K. M.; Singhasivanon, P.; Day, N. P. J.; White, N. J.; Anderson, T. J. C.; Nosten, F. Lancet 2012, 379, 1960−1966. (9) Bannister, L. H.; Hopkins, J. M.; Fowler, R. E.; Krishna, S.; Mitchell, G. H. Parasitol. Today 2000, 16, 427−433. (10) Nagaraj, V. A.; Sundaram, B.; Varadarajan, N. M.; Subramani, P. A.; Kalappa, D. M.; Ghosh, S. K.; Padmanaban, G. PLoS Pathog. 2013, 9, e103522. (11) Weissbuch, I.; Leiserowitz, L. Chem. Rev. 2008, 108, 4899− 4914. (12) Ursos, L. M. B.; DuBay, K. F.; Roepe, P. D. Mol. Biochem. Parasitol. 2001, 112, 11−17. (13) Egan, T. J.; Combrinck, J. M.; Egan, J.; Hearne, G. R.; Marques, H. M.; Ntenteni, S.; Sewell, B. T.; Smith, P. J.; Taylor, D.; van Schalkwyk, D. A.; Walden, J. C. Biochem. J. 2002, 365, 343−347. (14) Oliveira, M. F.; Timm, B. L.; Machado, E. A.; Miranda, K.; Attias, M.; Silva, J. R.; Dansa-Petretski, M.; de Oliveira, M. A.; de Souza, W.; Pinhal, N. M.; Sousa, J. J. F.; Vugman, N. V.; Oliveira, P. L. FEBS Lett. 2002, 512, 139−144. (15) Fitch, C. D.; Chevli, R.; Kanjananggulpan, P.; Dutta, P.; Chevli, K.; Chou, A. C. Blood 1983, 62, 1165−1168. (16) Becker, K.; Tilley, L.; Vennerstrom, J. L.; Roberts, D.; Rogerson, S.; Ginsburg, H. Int. J. Parasitol. 2004, 34, 163−189. (17) Balla, J.; Vercellotti, G. M.; Jeney, V.; Yachie, A.; Varga, Z.; Jacob, H. S.; Eaton, J. W.; Balla, G. Antioxid. Redox Signaling 2007, 9, 2119−2137. (18) Buller, R.; Peterson, M. L.; Almarsson, O.; Leiserowitz, L. Cryst. Growth Des. 2002, 2, 553−562. (19) Sullivan, D. J. Int. J. Parasitol. 2002, 32, 1645−1653. (20) Egan, T. J. Mini-Rev. Med. Chem. 2001, 1, 113−123. (21) Tilley, L.; Loria, P.; Foley, M. in Antimalarial Chemotherapy, Rosenthal, P. J., Ed.; Humana Press, Totowa, NJ, 2001; pp 87−122. (22) Slater, A. F. G. Pharmacol. Ther. 1993, 57, 203−235. (23) Hanscheid, T.; Egan, T. J.; Grobusch, M. P. Lancet Infect. Dis. 2007, 7, 675−685. (24) Ziegler, J.; Linck, R.; Wright, D. W. Curr. Med. Chem. 2001, 8, 171−189. (25) Pagola, S.; Stephens, P. W.; Bohle, D. S.; Kosar, A. D.; Madsen, S. K. Nature 2000, 404, 307−310. (26) Slater, A. F. G.; Swiggard, W. J.; Orton, B. R.; Flitter, W. D.; Goldberg, D. E.; Cerami, A.; Henderson, G. B. Proc. Natl. Acad. Sci. U. S. A. 1991, 88, 325−329. (27) Bohle, D. S.; Dinnebier, R. E.; Madsen, S. K.; Stephens, P. W. J. Biol. Chem. 1997, 272, 713−716. (28) Oliveira, M. F.; Kycia, S. W.; Gomez, A.; Kosar, A. J.; Bohle, D. S.; Hempelmann, E.; Menezes, D.; Vannier-Santos, M. A.; Oliveira, P. L.; Ferreira, S. T. FEBS Lett. 2005, 579, 6010−6016. (29) Gildenhuys, J.; le Roex, T.; Egan, T. J.; de Villiers, K. A. J. Am. Chem. Soc. 2013, 135, 1037−1047. (30) Straaso, T.; Kapishnikov, S.; Kato, K.; Takata, M.; Als-Nielsen, J.; Leiserowitz, L. Cryst. Growth Des. 2011, 11, 3342−3350. (31) Straaso, T.; Marom, N.; Solomonov, I.; Barfod, L. K.; Burghammer, M.; Feidenhans’l, R.; Als-Nielsen, J.; Leiserowitz, L. Cryst. Growth Des. 2014, 14, 1543−1554.

CONCLUSION In summary we have developed a new additive Amber FF compatible with FPD, a list of related small clusters and a crystal model. This FF was involved in MD simulations carried out in condensed phase to study the structure and dynamics of small FPD clusters and to identify the driving forces involved in the crystal cohesion. The dynamical aspect of MD reveals new HB patterns within FPD and between FPD entities, and demonstrates the importance of hydrophobic interactions in the cluster stability in condensed phase. Preventing crystal growth of hemozoin is an established therapeutic strategy against Plasmodium: developing new drugs, which would inhibit the formation of the embedded D1 motif, i.e., the smallest stable cluster synthon found in this work could be an attractive strategy for impairing hemozoin formation. This theoretical study provides new insights on FPD oligomer structure, dynamics, and stability. Simulation conditions involving aqueous conditions with alcohols or lipids, or at the interface of two phases, could provide additional pieces of information.85−87 Work is in progress in our group to address these different issues.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.cgd.6b00052. Force field related files (libraries, parameters, and leaprc files), Cartesian coordinates of the 17 input structures of MD simulations in the PDB file format, and 49 figures related to MD results (S1 to S49) (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: fyd@q4md-forcefieldtools.org. Tel: +33 (0)3-22827498. Fax: +33 (0)3-2282-7469. Funding

This study was funded by DGA (Direction Générale de l’Armement, Ministère de la Défense, France) and ANR Astrid (Agence Nationale de la Recherche, project ANR-12-ASTR003). The Ph.D. thesis of F.W. was supported by a grant from Conseil Régional de Picardie and European Regional Development Fund. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful to Dr. de Villiers (Stellenbosch University, South Africa) for helpful discussions, and to Dr. Cézard (UMR CNRS 7378, UPJV) for determining bonded FF constants related to iron in the Becker et al. FF.

■

ABBREVIATIONS DFT, density functional theory; FF, force field; FPD, cyclic dimer of ferriprotoporphyrin IX; HB, hydrogen bond; MD, molecular dynamics; RG, radius of gyration; RMSD, rootmean-square deviation; SASA, solvent accessible surface area

■

REFERENCES

(1) Sachs, J.; Malaney, P. Nature 2002, 415, 680−685. (2) World Health Organization. World Malaria Report 2015; WHO Press: Geneva, Switzerland, 2015. 2257

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

(32) Klonis, N.; Dilanian, R.; Hanssen, E.; Darmanin, C.; Streltsov, V.; Deed, S.; Quiney, H.; Tilley, L. Biochemistry 2010, 49, 6804−6811. (33) Bohle, D. S.; Dodd, E. L.; Stephens, P. W. Chem. Biodiversity 2012, 9, 1891−1902. (34) Martinez, C. R.; Iverson, B. L. Chem. Sci. 2012, 3, 2191−2201. (35) Dupradeau, F.-Y.; Pigache, A.; Zaffran, T.; Savineau, C.; Lelong, R.; Grivel, N.; Lelong, D.; Rosanski, W.; Cieplak, P. Phys. Chem. Chem. Phys. 2010, 12, 7821−7839. (36) Vanquelef, E.; Simon, S.; Marquant, G.; Garcia, E.; Klimerak, G.; Delepine, J. C.; Cieplak, P.; Dupradeau, F.-Y. Nucleic Acids Res. 2011, 39, W511−W517. (37) Wang, F.; Becker, J.-P.; Cieplak, P.; Dupradeau, F.-Y. RED Python: Object oriented programming for Amber force fields; 247th ACS National Meeting; Dallas, TX; March 16−20, 2014. (38) Cieplak, P.; Cornell, W. D.; Bayly, C. I.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1357−1377. (39) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision D.01; Gaussian, Inc., Wallingford, CT, 2009. (40) Bohle, D. S.; Debrunner, P.; Jordan, P. A.; Madsen, S. K.; Schulz, C. E. J. Am. Chem. Soc. 1998, 120, 8255−8256. (41) Rassolov, V. A.; Pople, J. A.; Ratner, M. A.; Windus, T. L. J. Chem. Phys. 1998, 109, 1223−1229. (42) Rassolov, V. A.; Ratner, M. A.; Pople, J. A.; Redfern, P. C.; Curtiss, L. A. J. Comput. Chem. 2001, 22, 976−984. (43) Kim, K.; Jordan, K. D. J. Phys. Chem. 1994, 98, 10089−10094. (44) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (45) This corresponds to the atoms with the 60, 61, and 65 indexes in the atom order of the input molecule. The description of the constraints applied during geometry optimization is described in the PyRED ’Project.config’ configuration file. The latter can be downloaded from the “F-94” R.E.DD.B. project. (46) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129− 145. (47) The atoms involved in the rigid-body reorientation algorithm (RBRA) are the following: the atoms with the 63, 1 (i.e., the Fe atom), 60 and 60, 1, 63 indexes for the monomeric heme derivative; with the 1, 5, 4 and 4, 5, 1 indexes for ethane; with the 6, 4, 1 and 1, 4, 6 indexes for propene; with the 5, 8, 11 and 11, 8, 5 indexes for pentanoic acid; with the 10, 7, 4 and 4, 7, 10 indexes for pentanoate and with the 5, 8, 11 and 11, 8, 5 indexes for pentane. The atoms involved in the RBRA procedure are automatically determined by the PyRED program, and their description is available in the “F-94” R.E.DD.B. project. (48) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269−10280, and http://upjv.q4md-forcefieldtools. org/REDServer-Development/resp/ accessed December 16, 2015. (49) 39 atoms are involved in a single intramolecular charge constraint. The atoms involved in intramolecular charge constraints are described in the PyRED ’Project.config’ configuration file. The latter can be downloaded from the “F-94” R.E.DD.B. project. (50) The Mol3 file format: http://q4md-forcefieldtools.org/ Tutorial/leap-mol3.php (accessed December 16, 2015). This file format contains RESP atomic charge values, force field atom types, the notion of head/tail for a molecular fragment, the topology (atom

connectivities) of the molecule/molecular fragment as well as the Cartesian coordinates, which were optimized by quantum mechanical computations. (51) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049−1074. (52) Allen, F. H. Acta Crystallogr., Sect. B: Struct. Sci. 2002, B58, 380− 388. (53) Seminario, J. M. Int. J. Quantum Chem. 1996, 60, 1271−1277. (54) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. Comput. Chem. 2005, 26, 1668−1688. (55) The format for the force field parameter modification file is described at http://ambermd.org/formats.html#frcmod (accessed December 16, 2015). The corresponding file generated for the FPD stereoisomers and analogs can be downloaded from the “F-94” R.E.DD.B. project. (56) Dupradeau, F.-Y.; Cézard, C.; Lelong, R.; Stanislawiak, E.; Pecher, J.; Delepine, J. C.; Cieplak, P. Nucleic Acids Res. 2008, 36, D360−D367. (57) PyRED automatically generates a script for the LEaP program (AmberTools). This script has been extended to construct and study numerous FPD-based molecular systems in condensed phase, and is available for download from the “F-94” R.E.DD.B. project. (58) Case, D. A.; Darden, T. A.; Cheatham, III, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Götz, A. W; Kolossvary, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H; Seetin, M. G.; Salomon-Ferrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12; University of California, San Francisco, 2012. (59) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (60) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. (61) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327−341. (62) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577−8593. (63) Simmerling, C.; Strockbine, B.; Roitberg, A. E. J. Am. Chem. Soc. 2002, 124, 11258−11259. (64) Roe, D. R.; Cheatham, T. E., III J. Chem. Theory Comput. 2013, 9, 3084−3095. (65) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33−38. (66) Jmol: an open-source Java viewer for chemical structures in 3D; http://www.jmol.org/ (accessed December 16, 2015). (67) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srinivasan, J.; Case, D. A.; Cheatham, T. E., III Acc. Chem. Res. 2000, 33, 889−897. (68) Miller, B. R., III; McGee, T. D., Jr.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E. J. Chem. Theory Comput. 2012, 8, 3314− 3321. (69) Nasipuri, D. Stereochemistry of organic compounds: principles and applications; New Age International Ltd., Publishers, New Delhi, India, 1991. (70) Pisciotta, J. M.; Sullivan, D. Parasitol. Int. 2008, 57, 89−96. (71) Egan, T. J.; Chen, J. Y.-J.; de Villiers, K. A.; Mabotha, T. E.; Naidoo, K. J.; Ncokazi, K. K.; Langford, S. J.; McNaughton, D.; Pandiancherri, S.; Wood, B. R. FEBS Lett. 2006, 580, 5105−5110. (72) Kapishnikov, S.; Berthing, T.; Hviid, L.; Dierolf, M.; Menzel, A.; Pfeiffer, F.; Als-Nielsen, J.; Leiserowitz, L. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 11184−11187. (73) Kapishnikov, S.; Weiner, A.; Shimoni, E.; Guttmann, P.; Schneider, G.; Dahan-Pasternak, N.; Dzikowski, R.; Leiserowitz, L.; Elbaum, M. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 11188−11193. 2258

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259

Crystal Growth & Design

Article

(74) Kapishnikov, S.; Weiner, A.; Shimoni, E.; Schneider, G.; Elbaum, M.; Leiserowitz, L. Langmuir 2013, 29, 14595−14602. (75) Olafson, K. N.; Rimer, J. D.; Vekilov, P. G. Cryst. Growth Des. 2014, 14, 2123−2127. (76) Ketchum, M. A.; Olafson, K. N.; Petrova, E. V.; Rimer, J. D.; Vekilov, P. G. J. Chem. Phys. 2013, 139, 121911−121919. (77) Durrant, M. C. Dalton Trans. 2014, 43, 9754−9765. (78) Rouge, P.; Dassonville-Klimpt, A.; Cézard, C.; Boudesocque, S.; Ourouda, R.; Amant, C.; Gaboriau, F.; Forfar, I.; Guillon, J.; Guillon, E.; Vanquelef, E.; Cieplak, P.; Dupradeau, F.-Y.; Dupont, L.; Sonnet, P. ChemPlusChem 2012, 77, 1001−1016. (79) Shahrokh, K.; Orendt, A.; Yost, G. S.; Cheatham, T. E., III J. Comput. Chem. 2012, 33, 119−133. (80) Autenrieth, F.; Tajkhorshid, E.; Baudry, J.; Luthey-Schulten, Z. J. Comput. Chem. 2004, 25, 1613−1622. (81) Bohle, D. S.; Dodd, E. L. Inorg. Chem. 2012, 51, 4411−4413. (82) Kauzmann, W. Adv. Protein Chem. 1959, 14, 1−63. (83) Dill, K. A.; Bromberg, S.; Yue, K.; Ftebig, K. M.; Yee, D. P.; Thomas, P. D.; Chan, H. S. Protein Sci. 1995, 4, 561−602. (84) Chandler, D. Nature 2005, 437, 640−647. (85) Olafson, K. N.; Ketchum, M. A.; Rimer, J. D.; Vekilov, P. G. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 4946−4951. (86) Vekilov, P. G.; Rimer, J. D.; Olafson, K. N.; Ketchum, M. A. CrystEngComm 2015, 17, 7790−7800. (87) Olafson, K. N.; Ketchum, M. A.; Rimer, J. D.; Vekilov, P. G. Cryst. Growth Des. 2015, 15, 5535−5542.

2259

DOI: 10.1021/acs.cgd.6b00052 Cryst. Growth Des. 2016, 16, 2249−2259