Binding of Fidarestat Stereoisomers with Aldose Reductase

The stereospecificity in binding to aldose reductase (ALR2) of two fidarestat {6fluoro-2',5'-dioxospiro[chroman-4,4'-imidazolidine]-2-carboxamide} stereoisomers [(2S,4S) and (2R,4S)] has been investigated by means of molecular dynamics simulations using free energy integration techniques. The difference in the free energy of binding was found to be 2.0 ± 1.7 kJ/mol in favour of the (2S,4S)-form, in agreement with the experimental inhibition data. The relative mobilities of the fidarestats complexed with ALR2 indicate a larger entropic penalty for hydrophobic binding of (2R,4S)-fidarestat compared to (2S,4S)fidarestat, partially explaining its lower binding affinity. The two stereoisomers differ mainly in the orientation of the carbamoyl moiety with respect to the active site and rotation of the bond joining the carbamoyl substituent to the ring. The detailed structural and energetic insights obtained from out simulations allow for a better understanding of the factors determining stereospecific inhibitor-ALR2 binding in the EPF charges model.


Introduction
The hyperglycemia observed in diabetes mellitus is the primary instigator of long-term diabetic complications such as retinopathy, neuropathy, nephropathy and cataracts [1][2].The elevated glucose concentration in blood causes increased flux through the polyol pathway and the activation of this pathway is believed to underlie various diabetic complications.Aldose reductase (ALR2), the first enzyme of the polyol pathway reduces glucose into sorbitol with concomitant conversion of NADPH to NADP+ [3][4].Excessive accumulation of intracellular sorbitol through the polyol pathway is linked to the pathogenesis of diabetic complications and prevention of sorbitol accumulation by inhibition of aldose reductase activity is an effective treatment for these complications [5][6].A large number of structurally diverse aldose reductase inhibitors have been synthesized [7][8][9][10][11][12][13][14][15][16][17], and many molecular modeling studies have been performed to understand, on a structural basis, the interactions between inhibitors and ALR2 [18][19].
Aldose reductase contains a (β/α)8 barrel motif with a large hydrophobic pocket, approximately 60 Å wide and 15 Å deep, which greatly favors aromatic and apolar substrates over highly polar monosaccharides.Crystal structures of ALR2 complexed with NADPH and diverse inhibitors reveal that this hydrophobic pocket is the only possible active site.The cofactor NADPH is bound in an unusual extended conformation across the barrel with the nicotinamide ring centered at the bottom of this cavity.Upon binding of NADPH, the enzyme undergoes a large conformational change involving the shifting of loop 7 to an orientation which appears to lock the coenzyme into place [20][21][22].
The active site can be divided into two regions.The first includes three proton-donating side chains, tyrosine, histidine, and tryptophan, along with the 4-pro-R hydrogen of the nicotinamide ring of NADPH, and provides an anion binding site at the bottom of the active site pocket.The second region forms a hydrophobic wall, composed of the hydrophobic residues phenylalanine, tryptophan, cysteine, alalanine, leucine, and offers the possibility of a π-π stacking interaction which could stabilize the aromatic ring system of a substrate and/or an inhibitor [23][24].A large number of aldose reductase inhibitors with significant activity in vitro and in animal models have been discovered but relatively few have been tested on diabetic patients and four of these, tolrestat, zopolrestat, ponalrestat and zenarestat were withdrawn from clinical trials due to the lack of efficacy or adverse side effects [25][26][27].Fidarestat, however, recently showed encouraging results in that it normalized erythrocytic sorbitol contents in neuropathy patients with no significant side effects [28].

Table 1. Functional residues around (2S,4S)-and (2R,4S)-forms of Fidarestat complexed with Human
Aldose reductase (ALR2).In the analysis of interactions, the following geometrical criterion for hydrogen bonds was used: a donor-acceptor distance 2.5 Å and a donor-hydrogen-acceptor angle of at least 135°.The IC50 value is taken from ref.In this study, we report the calculation of the relative free energies of binding to ALR2 of fidarestat and its (2R,4S) stereoisomer using molecular dynamics (MD) simulations in an explicit aqueous environment.X-ray structural and mass spectrometric studies of binding of these and other ligands to the same template have been reported [32][33][34].In our calculations, we have employed the thermodynamic integration (TI) approach [35] utilizing either Mulliken [36] or electrostatic potential fitted (EPF) [37][38] charges under CHARMM/MMFF [39][40][41][42] force fields for the perturbed part of the substrate.MD trajectories of a total of 8 ns were generated for each of the two fidarestat stereoisomers bound to ALR2.Structural characteristics of the fidarestat complexed with ALR2 were analyzed and the results were compared with the experimental inhibition data.

Partial charge and force constant calculation for fidarestat
The partial charges and force constants for the bonded energy terms were determined from a Hartree-Fock/6-31G* ab initio calculation.A full geometry optimization was carried out at the Hartree-Fock/6-31G* level using GAMESS [43].Optimization was started from standard CHARMM [39][40][41] values for the internal coordinates and the crystal structure geometry for each of the isomers of fidarestat [29][30][31][32][33]. Geometrical convergence was obtained after 43 cycles of optimization with a final gradient of 0.214×10 -3 Hartree/Bohr.Partial charges centered on the nuclei were determined by fitting the electrostatic potential derived from the ab initio calculation in the presence of special constraints.
Force constants for the internal coordinates bonded energy terms were obtained from the second derivative of the Hartree-Fock/6-31G* energy surface, the Hessian matrix, with respect to the internal coordinates, assuming that the nonbonded energy contributions to the force constants are small compared to the bonded energy terms.

Molecular dynamics simulations
MD simulations extending for over 8 ns in total were carried out with CHARMM [39].Force field parameters for the protein were taken from the CHARMM22 all-atom parameter set [41].The MMFF parameter set [42] in CHARMM was used to parameterize the charges of the NADP + atoms.The van der Waals parameters were taken from the MMFF/CHARMM22 parameters of analogous atoms.The interactions were cut off sharply with no switching function at 12 Å.SHAKE [44] module in CHARMM was used to restrain the bond lengths and a time step of 2 fs was used.To keep the temperature at 300 K, the system was coupled to a heat bath with a time constant of 0.1 ps.
Structures of (2S,4S)-Fidarestat and (2R,4S)-Fidarestat were taken from the crystal structures [31][32][33], in which fidarestat setereoisomers are complexed to ALR2.The complex was solvated in a periodic truncated octahedron using a modified TIP3P water model [45].Water outside a sphere of 20 Å radius around the center of the active site was cut off.This center was determined as the middle point between fidarestat and the C α of the side chains in the active site.The coordinates of water were optimized by steepest descent energy minimization in which ALR2 and each fidarestat was positionally restrained using a harmonic interaction and equilibrated for 100 ps with the protein atoms constrained to the crystal structure positions.All the protein atoms and waters inside a sphere of 12 Å radius around the center of the active site were allowed to move while the remaining atoms were constrained to the positions in the equilibrated structure using a harmonic potential.Prior to the free energy calculations, a final 100 ps equilibration run was then performed on the entire system using these constraints.
With these potential parameters under electrostatic potential fitted (EPF, see note in Table 2) charges, we got rather large fluctuations in the free energy change as well as in the position and orientation of the fidarestat between different simulations.In an effort to make these fluctuations smaller, we made another series of runs with all the Mulliken charges but keeping a zero net charge on all residues.
We consider the free energies in the reaction scheme in Figure 2. The free energy integration method has evolved as an efficient method with which to estimate differences in the free energy of binding between relatively similar molecules [46][47][48].It depends upon molecular dynamics simulations where one molecule is successively converted into the other form.This is accomplished [49][50][51] by adding into the simulations a perturbing force which depends upon a parameter λ that is 0 in the initial state and 1 in the final state.Perturbation is imposed by use of an improper dihedral potential: Here, Ψ is the improper dihedral that determines the angle between the C 1 -C 2 bond and the plane containing the N 1 , C 1 and O 1 atoms of carbamoyl moiety of fidarestat, K Ψ is the associated force constant and Ψ 0 = 36.2°,the result of the Hartree-Fock/6-31G* level optimization.The minimum of the potential varies from + Ψ 0 to -Ψ 0 during the conversion, the extreme values corresponding to the (2S,4S)F and (2R,4S)F stereoisomers.Use of the united atom method i.e., inclusion of hydrogen in the C 1 atom simplifies the conversion.The free energy difference G(λ) -G(0) can be calculated as a function of λ from the expression.
We are interested in G(λ) -G(0).We used the slow growth method to do the conversion and changed λ linearly in time from 0 to 1 in 20 ps.We also tried windowing techniques [48] instead of integration.The system was sampled for 1 ps at 10 discrete values of λ.The system was equilibrated between successive λ values.Half the computing time was then spent on sampling and half on changing λ and equilibrating the system.This gave results that were consistent with those obtained from the slow growth method.
For the (2R,4S)F to (2S,4S)F conversion in water we set the free energy change ∆G aq to zero.To check the integration procedure, we ran conversions of a system consisting of a periodic box with a modified TIP3P water model and one dissolved fidarestat molecule.During the equilibration, the density was adjusted to give a pressure of 1 atm.

Partial charges and force constant for fidarestat atoms
The calculated partial charges for the fidarestat atoms are shown in Table 2. Mulliken charges are given for comparison.These partial charges reproduce the electrostatic potential with a high level of accuracy, as measured by the χ 2 value of 2.85×10 -3 (A.U.), or by the χ 2 /N G of 1.60×10 -6 (A.U), where N G is the number of grid points to which the potential was fitted.The computed charges led to a dipole vector (0.43825, 1.21534, -2.31056×10 -6 ) (A.U.) in excellent agreement with the QM dipole (0.43925, 1.26534, -2.36220×10 -6 ).The correlation between the Mulliken charges and the electrostatic potential fitted (EPF) charges is excellent, with a correlation coefficient of 0.94 and a slope of 0.90.give rise to a very similar potentials, additional constraints were introduced for numerical stability and to select the most reasonable solution.The quantity to be minimized is: is the matrix of the inverse distances between the grid point r k and the ith atom position R i with a partial charge q i , and N A is the number of atoms.
For molecular dynamics simulation and energy minimization, the accuracy of the electrostatic potential at some distance from the molecule is more important than the actual values of the partial charges.Thus, the EPF charges were used instead of the Mulliken charges.Another reason for using EPF charges was that they are more consistent with the CHARMM charges [41].
Force constants were determined from the Hessian matrix.The correlation of the ab initio values with the corresponding CHARMM values is shown in Figure 3.The slope of the linear regression curve is approximately 2 and because of the form of the harmonic energy function [39], the force constants were divided by two before implementation in the CHARMM force field.It might be noted that off-diagonal terms can be of the same order of magnitude as the diagonal terms.Those terms are not accounted for by the present functional form of the potential energy function.

Molecular dynamics and Free Energy calculation
The values of ∆GALR2 from 16 simulations extending over 8 ns in total are shown in Table 3 and Figure 4.In Figure 4, the atom-positional RMSD of the fidarestat atoms in the trajectory structures from the initial structure is shown for the ALR2-(2R,4S)-Fidarestat and ALR2-(2S,4S)-Fidarestat complexes.3, respectively.
Figure 5 shows the free energy averaged over successive nanoseconds as a function of the integration parameter λ during various 500 ps runs.It can be seen that there is a systematic difference between the runs in the two directions and this hysteresis is a measure of the accuracy of the method.Additional errors in free energy integrations are discussed in [44][45][46].Consequently, the (2R,4S)F → (2S,4S)F and (2S,4S)F → (2R,4S)F integrations are averaged separately.The error in the final estimate is found by dividing the standard deviation between the four runs in each direction with the square root of the number of independent runs.
As a final estimate we take the averages over the (2R,4S)F ↔ (2S,4S)F runs with their appropriate signs.We then get a ∆G ALR2 value of -2.0 kJ/mol for the S series using EPF charges and -6.0 kJ/mol for the S* series using Mulliken charges with the signs for the (2R,4S)F to (2S,4S)F transition.The non-systematic errors in these numbers are estimated to ± 1.7 kJ/mol and ± 0.8 kJ/mol.In addition there is a systematic (hysteresis) difference between the (2R,4S)F → (2S,4S)F and (2S,4S)F → (2R,4S)F conversions of 2.9 and 2.0 kJ/mol, respectively.
The previous calculated free energy value [34] falls within the error limits of the final estimate for the S series, while it is clearly different from the result of the S* series simulations.The simulations (S* series) with the Mulliken charges show relatively small fluctuations between different runs while those with EPF charges (S series) differ more.This suggests that the EPF charges (S series) give a better description of this system.When the Mulliken charges are increased, we get a more precise binding with less fluctuations, but this occurs at the cost of more free energy to accommodate the (2R,4S)-form compared to the (2S,4S)-form of fidarestat.
We now turn to a description of the structural changes that occur during the molecular dynamics.To describe these, we use distances between some important residues in ALR2 and the atoms of fidarestat, as shown in Table 4.The orientation of the carbamoyl and the cyclic imide moiety, and the N1-C1-C2-O2 dihedral angle θ are also important quantities which are discussed below.In the structures from the runs with the Mulliken charges (S* series), the carbamoyl moiety of the fidarestats stays fixed with C1 close to the hydrophobic network cavity formed by Phe122, Cys298, Ala299, and Leu300.The variation in these distances among the different S* structures is only about 0.5 Å.In contrast, the cyclic imide moiety of both fidarestats changes its position during the simulations (Figure 6).In the first two structures it remains in the hydrophilic pocket close to Trp20 and Tyr48 but in subsequent simulations, it moves out of the pocket and over to the other side of the active site, near to His110 and Trp111.Once there, it moves only slightly back and forth during the remaining integration runs.
With EPF charges (S series), the carbamoyl group of the fidarestat appears to be attached less rigidly to the nearby region of Phe122.This is seen in Table 4 as the distances to these residues are similar to the corresponding distances in the S* series of runs.Further, as a consequence of the weaker electrostatic interactions in the system with EPF charges, the fluctuations in these distances between different structures are now 0.5 -1.7 Å compared to 0.2 -0.7 Å in the S* series.In this case, the cyclic imide moiety of the fidarestats remains in the hydrophilic pocket close to NADP+ (Figure 6).

Order parameter calculation
To describe the orientation of the fidarestat carbamoyl moiety in the active site we can use two vectors, one connecting the carboxyl-oxygen to the nitrogen and the other one perpendicular to the O 1 -C 1 -N 1 plane of the carbamoyl group (Figure 7).To compare the orientation in a series of structures we define two order parameters, S A and S B , associated with these vectors in the following way.The vectors are first normalized and then averaged over the series of structures.The average lengths of the vectors obtained in this way are taken as the order parameters, S A and S B .If the orientation of the fidarestat is exactly the same in all structures, the order parameters will be 1 and the more different the orientations of the Fidarestat ring are in the different structures, the closer the order parameters will be to 0.   5 it is seen that order parameter S A , is fairly close to 1 both for the S and S* structures, indicating that the O 1 -N 1 vector does not change orientation much during the integration runs.The S B order parameter is much smaller, especially for the S structures.This indicates that the carbamoyl plane has a quite different orientation in the different structures.However, if S B is calculated separately for the (2R,4S) and (2S,4S) conformations, we get values that are closer to one.This suggests that stereoisomerization occurs by rotation of the carbamoyl plane through 70°, just as the (2R,4S) → (2S,4S) conversion is accomplished by changing the improper dihedral angle around the O 1 -N 1 vector from +35° to -35°.The conversion occurs therefore with the ring fixed in the hydrophilic pocket and the free energy difference is small.The (N 1 -C 1 -C 2 -O 2 ) dihedral angle θ for the different structures is shown in Table 6.It lies close to one of the minima of the dihedral potential at ±60° or 180° for all the structures.If the changes of this dihedral angle during the isomerization are considered, it may be seen that the mechanism suggested previously can not explain the phenomenon.If the ring is kept fixed and the carbamoyl moiety is rotated 70° around the O 1 -N 1 vector to achieve the isomerization, the dihedral θ will change by +60° for the (2R,4S)F into (2S,4S)F conversion and -60° for the opposite transition.This can be calculated by the three-dimensional geometry and seen readily from a molecular model.That rotation will bring the dihedral from a potential minimum up on the top of a barrier, which is clearly unfavorable.The dihedral will therefore either remain in the initial well or make an additional 60° rotation over to reach the next well.This is also what is seen from Table 6.For the (2R,4S)F → (2S,4S)F conversions we have ∆θ approximately 0° or +120° while the change is 0° or -120° for the opposite transition.A single exception is that in which the (2R,4S)F → (2S,4S)F conversion involves a change in dihedral angle of +126°.In this case, the fidarestat ring changes orientation considerably.

Interaction energy between carbamoyl moiety and hydrophobic cavity
To interpret the differences between these free energy profiles, we calculated the average total interaction energy of carbamoyl group in fidarestat with its surrounding residues Phe122, Cys298, Ala299 and Leu300, using the trajectories of the free energy profile simulations.In the ALR2-Fidarestat complexes, the carbamoyl group is expected to be most often associated with Leu300 or Ala299.In Figure 8, the energy profile, and its decomposition into van der Waals and electrostatic contributions are plotted as functions of the O 1 (Fidarestat)-N(Leu300) distance.In both complexes, the total energy profile has a shape determined by the electrostatic interaction energy and similar to the corresponding free energy profile, with a barrier at intermediate distances and low values at the distances corresponding to Leu300 and Ala299.The interactions of the carbamoyl moiety with Leu300 and Ala299 stabilize the profile at the extreme positions, whereas the loss of electrostatic interactions at intermediate distances creates the large energy barrier.In the (2S,4S)F-ALR2 complex, inspection of the contributions from individual residues shows that the interaction with Leu300 is more stable energetically than that with Ala299 due to better interactions with residue Cys298, a result that is consistent with published observations [34].In the (2R,4S)F-ARL2 complex, Ala299 is more stable energetically and this is attributable to the electrostatic energy term.The main contribution to this stabilization is from interactions with residue Phe122 and Tyr123.The Cys298 residue favors the interaction of the carbamoyl moiety with Leu300, but to a smaller degree than in the (2R,4S)F-ALR2 complex.This analysis shows that the superior stabilization of the (2S,4S)F-ALR2 complex as compared to the (2R,4S)F-ALR2 complex stems from electrostatic interactions, particularly with Leu300.Inspection of the structures from the restrained-distance simulations, in which the carbamoyl moiety was restrained at Leu300, shows that the average distance between the carbamoyl moiety hydrogens and the backbone nitrogen of Leu300 is smaller in the (2S,4S)F-ALR2 complex.At the same time, the distance between atoms Leu300 N and Ala299 N is approximately 0.18 Å larger in the (2S,4S)F-ALR2 complex as a result of backbone dihedral angle distortion; this could contribute to the improved interactions between carbamoyl moiety and Leu300 in the (2S,4S)F-ALR2 complex.

Conclusions
We investigated two models with different fractional charges.The intent was to bind the substrate spirohydantoins better in the active sites and so reduce positional fluctuation as well as free energies.It was concluded that the rather weak charges in the EPF charges model for the fidarestat molecule, noncovalently linked to the active site, supports correct estimates of the relative free energy of conversion of the inhibitors.This model allows loose binding with fairly big variations in positions and energies, a result which is necessary to accommodate both the (2S,4S) and (2R,4S) isomers of fidarestat into the same position in the active site and still have a fairly small difference in free energy as required by the experimental inhibition constants [33][34].
This change of cyclic imide group orientation does not alter the interactions greatly, especially since the fidarestats enjoy reduced flexibility in the active site of ALR2.Consequently there is a very small free energy difference between the (2S,4S) and (2R,4S) stereoisomers bound in the active site.However, when more distorted molecules such as (2S,4R)F and (2R,4R)F bind, the difference in the orientation of cyclic imide and carbamoyl moieties has a more dramatic effect and these molecules have a greater difference in binding enthalpy than the (2S,4S) and (2R,4S) stereoisomers [33][34].
When the thermodynamic integration method was applied in a study of binding free energy differences of (2S,4S)F and (2R,4S)F to ALR2, the preferential binding of (2S,4S)F over (2R,4S)F was reproduced in agreement with the experimental inhibition data.Slight differences between the fidarestats were seen in the loss of mobility and orientation of the inhibitors upon binding to ALR2.This loss was more pronounced for (2R,4S)F, leading to a more unfavorable entropic contribution to the free energy of binding for this molecule, which partially explains its lower binding affinity.

Figure 4 .
Figure 4. Atom-positional RMSD of the fidarestat in the trajectory structures from the initial structure under EPF and Mulliken charges, for the simulated ALR2-Fidarestat complex.

8 Figure 5 .
Figure 5. Average free energy variation with λ during the thermodynamic integration of 500 ps run in the case of (A) EPF charges and (B) Mulliken charges during (2R,4S)F ↔ (2S,4S)F conversions i. e. Sand S* series simulations described in Table3, respectively.

Figure 7 .
Figure 7.A schematic representation of the vectors and the order parameters to define the orientation of the inhibitor at the active site.S A = |ΣV 1 /N|, S B = |ΣV 2 /N|, where N is the number of structures and V 1 and V 2 are the vectors.

Figure 8 .
Figure 8. Change in the average energy of interaction between the carbamoyl moiety and surrounding residues, in the translocation of carbamoyl from the region of Leu300 to that of Ala299 in the (A) ALR2-(2R,4S)F and (B) ALR2-(2S,4S)F complex.The energy curves have been shifted to 0.0 at the point close to Leu300.
CA b HB c AR d PH e DC f Dist a CA b HB c AR d PH e DC f a Nearest distance (Å) between atoms of the Fidarestat and the residue of Human ALR2.b Contact surface area (Å 2 ) between the substrate and the residue of Human ALR2.c Hydrophilic-hydrophilic contact (hydrogen bond).d Aromatic-aromatic contact.e Hydrophobic-hydrophobic contact, f Hydrophobic-hydrophilic contact (destabilizing contact).+/-indicates presence/absence of a specific contacts.* indicates residues contacting ligand by their side chain (including C α atoms).

Table 2 .
Mulliken versus electrostatic potential fitted (EPF) a charges for fidarestat atoms Partial charges were determined by fitting the electrostatic potential derived from the ab initio calculation in the presence of special constraints.Since these problems, in general, degenerate, in the sense that two different charges distribution can a

Table 3 .
Free energy differences from the simulation.

Table 4 .
Average distances between the atoms of the fidarestat and the ALR2 active site a I is Hydrophilic atom type : N and O that can donate and accept hydrogen bonds, II is Acceptor : N or O that can only accept a hydrogen bond, III is Donor : N that can only donate a hydrogen bond.b Dist is distance (Å) between the fidarestat and ALR2 atoms.

Table 5 .
The order parameters S A and S B , as described in the text, calculated for normal and higher charged structures.