Tautomerism in folic acid: Combined molecular modelling and NMR study

Accepted Manuscript Tautomerism in folic acid: Combined molecular modelling and NMR study

Gergana Gocheva, Nikolay Petkov, Andrea Garcia Luri, Stoyan Iliev, Nikoleta Ivanova, Jasmina Petrova, Yavor Mitrev, Galia Madjarova, Anela Ivanova PII: DOI: Article Number: Reference:

S0167-7322(19)33325-2 https://doi.org/10.1016/j.molliq.2019.111392 111392 MOLLIQ 111392

To appear in:

Journal of Molecular Liquids

Received date: Revised date: Accepted date:

13 June 2019 15 July 2019 18 July 2019

Please cite this article as: G. Gocheva, N. Petkov, A.G. Luri, et al., Tautomerism in folic acid: Combined molecular modelling and NMR study, Journal of Molecular Liquids, https://doi.org/10.1016/j.molliq.2019.111392

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Tautomerism in folic acid: combined molecular modelling and NMR study Gergana Gocheva,1 Nikolay Petkov,1,2 Andrea Garcia Luri,1,3 Stoyan Iliev,1 Nikoleta Ivanova,1,4 Jasmina Petrova,1 Yavor Mitrev,2 Galia Madjarova,1,*Anela Ivanova1,* 1

Laboratory of Quantum and Computational Chemistry, Faculty of Chemistry and Pharmacy, Sofia

Institute of Organic Chemistry with Centre of Phytochemistry, Bulgarian Academy of Sciences, Acad. G.

IP

2

T

University, 1 James Bourchier blvd., 1164 Sofia, Bulgaria

Bonchev str., bl. 9, 1113, Sofia, Bulgaria 3

CR

Faculty of Chemistry, University of Barcelona, 1-11 Martí i Franquès blvd., 08028 Barcelona, Spain

4

University of Chemical Technology and Metallurgy, 8 St. Kliment Ohridski blvd., 1756, Sofia, Bulgaria

US

*Corresponding authors:

E-mail: [email protected]; Tel.: ++35928161431; FAX ++35929625438

AN

E-mail: [email protected]; Tel.: ++35928161520; FAX ++35929625438

M

Abstract

ED

The determination of the predominant tautomeric form of folate in solution is addressed theoretically and experimentally. DFT and atomistic molecular dynamics simulations are used to compare the

PT

most stable tautomers of folate. B3LYP/6-311+G(d) model chemistry shows that the lactam form of folate FA-N3 is the most energetically favorable tautomer. The lactam FA-N1 differs by about 2.5

CE

kcal/mol. The conformation of the molecule does not influence markedly the energy ordering of the tautomers. Molecular dynamics results for the two lactams in saline outline similar behavior of the two molecules and flexible docking predicts comparable interactions with the folate receptor- (FR-

AC

). 1H NMR spectra of folic acid in water and DMSO solvents also identify the lactam tautomers as present in the solutions. Both 1H and 13C chemical shifts show fast exchange of protons between the N1 and N3 positions of pterin.

Keywords Folate tautomers, DFT calculations, atomistic molecular dynamics simulations, NMR spectra, molecular conformation

ACCEPTED MANUSCRIPT

CR

IP

T

Graphical abstract

Introduction

Folic acid (FA) is one of the vitamins B and its conjugates with drugs [1] are considered as

US

promising candidates for active targeting drug delivery systems for cancer treatment [2]. The molecule of FA consists of three subparts: pterin, 4-aminobenzoyl, and glutamic acid residues. The

AN

latter is usually deprotonated in physiological conditions to give folate (with charge -2), which is the bioactive form [3]. The chemical structure of folic acid implies substantial structural diversity. The

M

pterin fragment has a possibility for different lactam and lactim tautomers [4] (Figure 1), depending on the position of attachment of one of the polar hydrogen atoms. Their relative stability defines the

ED

ratio of existence in the human organism. The presence of different FA tautomers may influence its biological activity [5], for example, the highly specific binding to target proteins – the folate

PT

receptors or the dihydrofolate reductase (DHFR) [6]. Literature search of the structure of folic acid results in a puzzling outcome because in the various sources of information different tautomers are

CE

shown as representative. For example, in the 9th Edition of the European Pharmacopoeia the lactam tautomer FA-N1 is presented [7], whereas the NIST Chemistry WebBook [8] and the 17th Edition of

AC

the Japanese Pharmacopoeia Database [9] show the lactim form. There are several experimental reports on the existence of tautomers of pterin. The comparison of UV spectra and pKa for specifically methylated pterins with pterin itself outlines the N3 lactam tautomer (Figure 1) as the most stable one [10-12]. The feasibility of the various tautomers may be characterized theoretically as well. This problem has been tackled in several previous studies, which are limited mainly to the lactam-lactim tautomerism in pterin alone. The study of Gready [13] on pterin structures includes gas phase geometry optimization with HF/STO-3G, followed by HF/3-21G single point energy calculations of 30 pteridines. The set includes tautomers, protonated, deprotonated, and hydrogenated forms. The

ACCEPTED MANUSCRIPT results reproduce the experimental findings that the tautomer N3 is the most stable one. The lactamlactim tautomerism of 6-methyl-7,8-dihydropterin is studied theoretically as a model for the proton transfer mechanism of folate to DHFR [14]. The computational protocol is the same as in the work of Gready [13] but additional single point calculations are done with the basis set 6-31G. The work shows that the use of the split-valence basis sets 3-21G and 6-31G leads to correct reproduction of the energy ordering of the tautomers, the lactam form being more stable than the lactim one by 1.4

T

and 1.0 kcal/mol, respectively.

IP

The tautomerism of the pterin system is studied also in relation to the binding of folate to DHFR [15]. Hartree-Fock geometry optimization with basis sets STO-3G, 3-21G, and 6-31G(d) is

CR

carried out in the gas phase for the pterin tautomers N1 and N3 (Figure 1). The results show that the most stable molecule is N3. N1 is higher in energy by 5.50 kcal/mol with the basis set STO-3G,

US

which increases to 6.17 kcal/mol with 6-31G(d). This illustrates detectable influence of the basis set. A theoretical study of the relative energies of all tautomers of neutral pterin is made by Soniat

AN

and Martin [4]. The rest of FA is replaced by a methyl group at position 6 (Figure 1). The authors systematically study a set of 5 tautomers using the DFT functional B3LYP and comparing the basis

M

sets 3-21G(d), 6-31G(d), and 6-31G(d,p). In all calculations the aqueous environment is treated with the implicit PCM solvation model. For the different tautomers within this set, two conformers are

ED

obtained by changing the position and orientation of the hydrogen atom of the lactim hydroxyl group, which results in total of 10 simulated structures. The isomer with the lowest energy is 6-

PT

methylpterin with a proton bound to N3 (Figure 1). The other lactam, with the hydrogen bound to N1, is also considered significant in aqueous solution, as it has an energy difference of about 1

CE

kcal/mol. The ordering of the tautomers remains the same except for a change of one tautomer with the 3-21G(d) basis set. The rest of the isomers are less probable, since their relative energies vary

AC

from ca. 4 to ca. 16 kcal/mol above the lowest one. Pterin tautomers behave as weak acids in aqueous solution and Soniat and Martin extend their study to the tautomers of anionic pterin [16]. The negative charge is introduced by removing a proton and the calculations predict the most probable position of deprotonation. The structures are calculated

at

the

B3LYP/3-21G(d),

B3LYP/6-31G(d)//B3LYP/3-21G(d),

and

B3LYP/6-

31+G(d)//B3LYP/3-21G(d) levels of theory. The anionic 6-methylpterin structure with the lowest energy is the one with the negative charge on the oxygen of the carbonyl group both in the lactam and in the lactim form. It is also found that the particular tautomer is not important for the position of

ACCEPTED MANUSCRIPT the negative charge, as the excess electron density is mostly localized in the vicinity of the atom N(2a) (Figure 1). Later, Nekkanti and Martin report a DFT study of the tautomers of cationic pterin [17], which are generated by protonating 6-methylpterin at all possible positions. The cationic molecules are optimized at the B3LYP/6-31+G(d,p) level of theory and the six most stable structures are also optimized using MP2/6-31+G(d,p) for comparison. The MP2 results agree with the B3LYP data and

T

both coincide on the most stable structure. All the calculated relative energies using the two levels of

IP

theory do not differ much – just within 1.5 kcal/mol. The most stable pterin cation tautomer is the lactam N3 with protonated N1.

CR

Pterin neutral radicals, generated by eliminating an N-bonded hydrogen atom, are also studied [18]. They are assumed to be more stable than those obtained by removing a C-bonded

US

hydrogen atom. Geometry optimization with UB3LYP/6-31G(d) is done both in the gas phase and in implicit SMD aqueous environment. It is followed by single point energy calculations with

AN

UB3LYP/ 6-311+G(2d,p). The most stable structure in the gas phase is the radical N(2a), followed by N3 separated in energy by 4.8 kcal/mol. In water, the energy ordering of these two structures is

Recently,

M

reversed and the difference reduces to about 0.7 kcal/mol. a B3LYP/6-31+G(d) and M06-2X/6-31+G(d) study, including geometry

ED

optimization in the gas phase of different conformers of folic acid, is reported [19], followed by calculation of the valence photoelectron spectrum obtained with general-R-SAC-CI/D95(df,pd). The

PT

possibility for different tautomers within the studied conformers is not addressed. The above literature overview of the lactam-lactim tautomerism related to folic acid shows

CE

that in fact only the pterin fragment or its hydrogenated derivatives have been modelled so far, while the rest of the molecular structure of folic acid has not been included. Hence, it is not known whether

AC

and how it will influence the lactam-lactim energy ordering. The scope of this study is to characterize computationally and experimentally the lactam-lactim tautomerism in the entire molecule of folate and also to trace the influence of the molecular conformation on the lactam-lactim tautomerization by DFT calculations, molecular dynamics simulations, and NMR spectra. The experimental NMR data are compared to the computational predictions.

ACCEPTED MANUSCRIPT Models and Methods Molecular models Six different lactam and lactim tautomers may exist in the pterin fragment of the molecule of

CE

PT

ED

M

AN

US

CR

IP

T

FA, the chemical structures of which are presented in Figure 1.

Figure 1: Chemical formulae of the studied tautomers of folate with notations used throughout the

AC

text, starting with the two lactams on the first row; monitored dihedral angles are denoted on the structure of FA-N3; the numbering of the atoms is shown on the structure of FA-N1.

The potential energy surface of folic acid is, however, even more complex. In addition to the tautomers, the molecule of FA contains nine single bonds along the longest molecular backbone that allow rotation and, hence, each tautomer may have many different conformers. In the current study, the effect of both the type of tautomer and of the conformation is evaluated. The choice of initial structures for the DFT calculations of folic acid is based on the results from a conformational search on the nine rotatable bonds of the tautomers FA-N1 and FA-N3 made with the force filed AMBER99 [20] in explicit TIP3P [21] water. 5000 randomly generated and geometry optimized structures are

ACCEPTED MANUSCRIPT screened using the Metropolis method [22] at temperature 310 K to sample the conformational space. HyperChem 7.0 [23] is used for the conformational search. More details on the computational procedure are given in our recent publications [24, 25]. The same procedure is applied to generate the stable conformers of the two tautomers of neutral folic acid in vacuo. The representative structures from the conformational search in water are used as starting geometries for geometry optimization and calculation of NMR chemical shifts in water, while those from the gas-phase

T

conformational analysis represent the initial set for the respective NMR calculations in DMSO.

IP

Two different stable conformers of FA-N1 and FA-N3 (Figure 2), found by the conformational search, are used as starting geometries for the DFT calculations. For each of these two initial

CR

structures, the six tautomers (Figure 1) are built and their geometry is optimized. Only one position of the hydrogen atom in the hydroxyl group of lactim and in the imine group of the guanidinyl forms

AC

CE

PT

ED

M

AN

US

is modeled. Overall, a set of 12 molecules is investigated.

Figure 2: The two different conformers Conf-1 (top) and Conf-3 (bottom) of the respective tautomers of FA used as initial structures for generation of all tautomers shown in Figure 1.

All calculations in water are done for the dianionic form of folic acid (Figure 1) and those in the gas phase and in DMSO – for the neutral one.

ACCEPTED MANUSCRIPT Computational protocol The geometry optimization and calculation of the vibrational frequencies are performed with B3LYP/6-311+G(d,p) in implicit water. The diffuse functions are indispensable since the modelled molecules are dianions. Natural bond orbital (NBO) analysis is performed, for which NBO 3.1 is used [26]. The NMR chemical shifts are calculated with B3LYP/6-311+G(2d,p) and with modified BLYP/6-311+G(2d,p) [27] in implicit water and DMSO. The proton NMR chemical shifts are

T

estimated with the modification for protons (WP04) and those of carbon – with the modified

IP

functional for carbons (WC04), as recommended in the paper. The chemical shifts are averaged over all representative conformers of each set, using their Boltzmann weights. The geometry optimization,

CR

frequency calculations, NBO analyses, and NMR spectra predictions are done with PCM [28] to describe the solvent and with the program package Gaussian 09 [29]. Two indices of delocalization

US

(MCDI and HOMA) are computed with Multiwfn [30] for the optimized geometries of FA-N1 and FA-N3. The NAO basis is employed for these estimates.

AN

The MD simulations are performed in NPT ensemble in 5 x 5 x 5 nm periodic box, where one molecule FA-N1 or FA-N3 is dissolved in saline to mimic physiological conditions. The negative

M

charges of FA are neutralized with Na+ cations. The force field for FA and for the inorganic ions is CGenFF [31] and that for water is TIP3P [21]. The average temperature is 310 K (V-rescale

ED

thermostat [32]) and the pressure is 1 bar (Berendsen barostat [33]). The full details of the MD simulation protocol can be found in [24, 25]. The production part of the MD trajectories is 100 ns

PT

and data for statistical analysis are saved every 0.2 ps. Gromacs 5.1.2 [34] is used to generate the trajectories and VMD 1.9.1. [35] – for visualization.

CE

The most populated MD geometries of the two most stable tautomers are docked into the Xray structure (PDB ID: 4LRH) of the folate receptor-α [36]. Flexible docking is done with AutoDock

AC

Vina 1.5.6 [37] with 300 initial randomized structures of each tautomer and grid size 2.0 x 1.5 x 2.7 nm, originating at the geometrical center of the active binding site of the receptor. NMR measurements All spectra are recorded on Bruker Avance II+ spectrometer (14.09 T magnet), operating at 600.11 MHz 1H frequencies, equipped with 5 mm BBO probe with z-gradient coil. The temperature is maintained at 310 K, using Bruker B-VT 3000 temperature unit with airflow of 535 L/h. All chemical shifts are reported in parts per million (ppm), referenced against tetramethylsylane (TMS, 0.00 ppm) for DMSO solutions and trimethylsilylpropanoic acid (TSPA, 0.00 ppm) for water

ACCEPTED MANUSCRIPT solutions. The aqueous sample was prepared by dissolving 11 mg folic acid in 1 ml distilled water and the pH was adjusted to 7 with K2CO3. The DMSO sample was prepared by dissolving the same quantity (11 mg/ml, 25 mM) of folic acid in the corresponding deuterated solvent. Folic acid was purchased from Sigma-Aldrich and recrystallized from hot water. A capillary with acetone-d6 was added in the NMR tube for lock and shimming in the water sample.

IP

T

Results and discussion

DFT calculations

CR

Due to the flexibility of the folic acid molecule, the geometry optimization of the tautomers formed from the two conformers ends in different local minima on the potential energy surface. The

AC

CE

PT

ED

M

AN

US

two optimized structures of each tautomer are superimposed in Figure 3.

Figure 3: DFT optimized geometry of folic acid tautomers started from two different conformations: the blue structures are built from Conf-1 and the black ones – from Conf-3.

ACCEPTED MANUSCRIPT The two initial conformers (Figure 2) differ not only in the shape of the molecule but also in the conformation of the amide bond. Conf-1 is characterized by a more compact form and trans amide bond, while Conf-3 has a more extended structure and cis amide bond. After the geometry optimization the tautomers built from Conf-1 open to a more extended conformation, similar to the one of Conf-3. In both cases the cis or trans position of the amide bond remains unchanged during the geometry optimization. Even though the overall molecular shape of each tautomer is similar,

T

there are some dissimilarities between the optimized structures obtained from the two different

IP

conformers. One of them is the position of the phenyl ring – it closes a slightly smaller (by about 10) angle (2, denoted on Figure 1) with pterin in the tautomers of Conf-3. There are also

CR

differences in the orientation of the carboxylate ends of the molecules. The rotation of the glutamate fragment is not the same in the tautomers formed from the two conformers, which arises from

US

different configuration of the amide bond: cis in the structures from Conf-3 and trans in those from Conf-1. This is in accordance with the relatively high rotation barrier (5.81 kcal/mol) between the

AN

two isomers of FA [24]. The last major difference in this terminal part is related to the rotation of the -carboxylate group of the tautomers. The Guanidinyl lactim structure has a more perpendicular

M

position of the pterin fragment with respect to the rest of the molecule compared to the other tautomers. There are no other significant differences between the optimized geometries of the

ED

various tautomers generated from the same conformer. The optimized structures have Z-shape for three of the tautomers (FA-N1, FA-N3, and

PT

Lactim), V-shape for one tautomer (Guanidinyl lactam), and U-shape for Guanidinyl lactim. Although small, this difference in shape may be important for the binding of the molecule into the

CE

active site of the folate receptor. It should be noted that the shape of Guanidinyl lactim resembles the one of bound FA from the crystallographic structure of its complex with the folate receptor- [36].

AC

The shape of the lactam tautomer FA-N3, however, coincides with its most stable Z-type geometry obtained at body temperature in saline [25]. This indicates that the conformation of FA is likely to change from the thermodynamically most stable one upon binding to the receptor but not so much upon temperature increase. Table 1 summarizes the relative free energies of all 10 studied tautomers calculated with DFT. The most stable tautomer is FA-N3, followed by the other lactam FA-N1. These results are in agreement with the previous studies on the isolated pterin subsystem [4-13, 15]. Comparing the tautomers that are built from the same conformer (Gset in Table 1), the lowest in energy is always the FA-N3 structure, no matter which conformer it originated from. The second in energy is FA-N1,

ACCEPTED MANUSCRIPT then the guanidinyl lactam, followed by the lactim, and the highest in energy is the guanidinyl lactim tautomer. The energy difference between one tautomer and the next is preserved for both conformers. The influence of the molecular conformation amounts to only 0.3 – 0.6 kcal/mol relative energy stabilization or destabilization of a given tautomer. This means that, regardless of the original conformer, the tautomers are always ordered in the same way in energy. This ordering coincides with the one found earlier for isolated pterin at the same computational level [4], which confirms the

T

minor influence of the molecular conformation of FA on the stability of the various tautomers.

IP

The third column in Table 1 compares the energy of all tautomers with respect to the obtained most stable structure in the entire series – tautomer FA-N3 optimized from the initial conformer with

CR

a more compact structure and trans amide bond. The different initial geometry and high flexibility of the structure strongly influences the final results in the sense that all tautomers obtained from Conf-3

US

are systematically shifted up in energy by 12.5 kcal/mol. One of the reasons for that may be the different configuration of the amide bond. Since the difference in the energy of cis and trans FA is of

AN

the order of 6 kcal/mol (for the most populated geometry in saline at body temperature) [24], this energy change can be attributed to a significant extent to the variance of the amide bond

M

configuration. The rest is due to contributions distributed along the other conformational degrees of

ED

freedom.

Table 1: DFT relative free energy of tautomers calculated with respect to the most stable structure in

PT

the entire series of studied molecules (Gall) and only in the respective set of tautomers (Gset). Tautomer

Conf-1

FA-N3 FA-N1 Guanidinyl Lactam Lactim Guanidinyl Lactim-N1 Guanidinyl Lactim-N3 FA-N3 FA-N1 Guanidinyl Lactam Lactim Guanidinyl Lactim-N1 Guanidinyl Lactim-N3

AC

CE

Conformer

Conf-3

Gall (kcal/mol) 0.00 2.44 6.84 7.13 17.61 32.66 12.49 15.22 18.69 20.3 30.15 32.89

Gset (kcal/mol) 0.00 2.44 6.84 7.13 17.61 32.66 0.00 2.73 6.20 7.81 17.66 45.38

ACCEPTED MANUSCRIPT From the results in Table 1 (fourth column) it can be seen that under certain conditions the second lactam form FA-N1 can also be realized since it is only 2.44 kcal/mol (from Conf-1) or 2.73 kcal/mol (from Conf-3) above the most stable FA-N3. However, the relative population of FA-N1 may not be significant since the energy difference exceeds the thermal energy at body temperature. The free energy difference between the two most stable lactam tautomers in the pterin molecule only is 1.98 kcal/mol, indicating that the flexible part of the folate molecule, i.e., 4-aminobenzoyl and

T

glutamic residues, adds to the stabilization of FA-N3.

IP

In an attempt to explain the relative energy stability of the tautomers, the NBO charge distribution is analyzed. The sum of the charges of the atoms from the pterin part, the 4-

CR

aminobenzoyl ring, the amide bond, and the glutamic part for all tautomers obtained from Conf-1 is shown in Table 2. The total charges of these subgroups in the studied tautomers do not differ much.

US

The pterin part of the molecules is slightly negatively charged, more substantial electron density is localized on the amide bond, and the most significant contribution to the total charge of -2 is

AN

expectedly concentrated in the glutamic part of the tautomers. Hence, the relative destabilization of

M

most of the tautomers is not electrostatically driven.

Table 2: Total charge (in a.u.) on pterin, aminobenzoyl, amide bond, and glutamic part of the

ED

structure of all tautomers obtained from Conf-1. Pterin

4-Aminobenzoyl

Amide bond

Glutamic acid

FA-N3

-0.043

0.056

-0.274

-1.739

-0.046

0.058

-0.273

-1.739

-0.046

0.058

-0.274

-1.739

-0.042

0.055

-0.273

-1.739

Guanidinyl Lactim-N1

-0.046

0.056

-0.271

-1.739

Guanidinyl Lactim-N3

-0.039

0.052

-0.273

-1.740

FA-N1

AC

Lactim

CE

Guanidinyl Lactam

PT

Tautomer

Spatial crowding is also not a likely reason because there are no close-lying atoms in any of the tautomers (Figure 3). Since electrostatic and structural parameters do not provide a clue for the relative stabilization of FA-N3 compared to FA-N1, the cause may be rooted in the degree of -electron density delocalization within the pterin ring. To check that, two aromaticity indices are computed for the optimized structures of FA-N1 and FA-N3 – the multicenter delocalization index (MCDI) and the

ACCEPTED MANUSCRIPT harmonic oscillator model of aromaticity (HOMA) index [38]. The former accounts for stabilization via delocalization of -electrons over more than two atomic centers while the latter quantifies the degree of aromaticity by measuring the differences of bond lengths. These two indices have been shown [39] to provide the most reliable results for fused heterocycles such as those forming pterin. (MCDI)1/10 was taken to allow comparison to rings with different number of atoms. Higher values of the two indices signify more enhanced aromaticity.

T

The following values were obtained: for FA-N1 (MCDI)1/10 = 0.5372, HOMA = 0.6526 and for

IP

FA-N3 (MCDI)1/10 = 0.5447, HOMA = 0.7319. It is evident that both in terms of molecular structure and of electron density delocalization FA-N3 has more pronounced aromatic character than FA-N1.

CR

This can be outlined as the probable origin of the lower energy of the former tautomer. To assess the extent of aromatic stabilization in pterin, it is meaningful to compare these magnitudes to those of

US

benzene and naphthalene (calculated with the same computational protocol): (MCDI)1/6 = 0.6658, HOMA = 0.9848 and (MCDI)1/10 = 0.6203, HOMA = 0.8283, respectively. As could be expected, the

AN

-electrons are somewhat less delocalized in naphthalene than in benzene. In pterin, the presence of heteroatoms decreases the degree of delocalization even further. This is reflected both in HOMA and

M

in MCDI but the values of the two indices show that the electron density is still partially distributed

retains some aromatic character.

PT

Molecular dynamics simulations

ED

over multiple centers in the ring. Overall, the delocalization assessment implies that the pterin system

The above summary of the DFT calculations clearly outlines the two lactams FA-N1 and FA-

CE

N3 as the most energetically stable ones and it may be assumed that they are those, which could be observed experimentally. Therefore, they are studied further by atomistic molecular dynamics

AC

simulations to characterize their structure in saline. Figure 4 shows the evolution of the RMSD of the atomic coordinates of the tautomers during a 100 ns productive trajectory. The results reveal that both molecules rapidly change their geometry over time but the energetically most favorable tautomer FA-N3 has faster dynamics than FA-N1. In the second half of the simulation there are two substructures of FA-N1 that become pronounced and interchange more slowly between each other. In FA-N3, three substructures are probably populated, but they differ less in geometry than those of FA-N1. This is evidenced by the highly overlapping peaks in the RMSD distribution.

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

Figure 4: (left) Time evolution and (right) distribution of the RMSD of the atomic coordinates of the

AC

CE

PT

lactams FA-N1 and FA-N3 [24, 25] compared to the initial minimized structures.

Figure 5: Values of the dihedral angles θ1 (left) and θ2 (right) (denoted in Figure 1) extracted from the last 50 ns of the MD trajectories of FA-N1 and FA-N3 [25].

ACCEPTED MANUSCRIPT For deciphering the observed difference in the structures, Figure 5 presents the evolution of the dihedral angles between the pterin part and the rest of the molecule. The first dihedral angle (θ1) labels rotation around the single bond that connects the pterin and the methylene group and the second dihedral angle (θ2) represents rotation about the single bond connecting the methylene and the NH group (Figure 1). These two angles are very important due to the fact that they largely define the overall shape of the molecule [25] and, hence, may influence the interactions of FA in the active

T

site of folate receptor-, especially those of the pterin part [36]. θ1 of FA-N3 takes values in a broad

IP

range 0÷360° while θ1 of FA-N1 is localized in the range 120÷220°. The rotation about the second single bond is also distinguishable in the two lactams: FA-N1 is characterized with one peak with

CR

strong preference for the conformations with θ2∼280÷300° and FA-N3 has two maxima – a more populated one in the range 60÷110° and one spanning 260÷300°. These results indicate clear

US

dominance of a significant local bending in the molecular backbone, which is different for the two tautomers and explains their dissimilar shape (Figure 3). The distribution of the data also implies that

AN

the structural rearrangement of pterin in FA-N3 may take place more easily than in FA-N1. The differences in the structural dynamics of the two tautomers may be due to the formation of

and their evolution is given in Figure 6.

M

different intra- or intermolecular (with the water molecules) hydrogen bonds. They are calculated

ED

In general, the two isomers are strongly hydrogen-bonded but there are also some differences between them, mostly related to the intramolecular hydrogen bonds. FA-N1 forms longer living H-

PT

bonds within the molecule than FA-N3, which may explain the slower structural dynamics of the former. This assumption is confirmed by the fact that the only intramolecular H-bonds in FA-N1 are

CE

between the NH-group of the amide bond and the two oxygen atoms of the glutamate -carboxylate residue. One of these bonds is present throughout almost the entire trajectory and in some periods it

AC

is supplemented by a second H-bond between the two groups. This through-space coupling is very likely to restrict the conformational freedom of FA-N1 to some extent. FA-N3 forms the same intramolecular hydrogen bonds but they are more short-lived in total and only one of them is present most of the time. On the other hand, the preferred length of the intramolecular hydrogen bonds in FA-N3 (0.19 nm) is a bit shorter than that in FA-N1 (0.20 nm). This is one of the possible reasons for stabilizing this positional isomer against FA-N1. It should also be noted that a third intramolecular H-bond forms only in FA-N3. It is between the amide NH group and the carbonyl oxygen of pterin. The fact that it is present only in this tautomer is in line with its more expressed structural flexibility.

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

PT

Figure 6: Evolution of the number of intramolecular (top) and folate-water (bottom) hydrogen bonds

CE

in the two positional isomers of FA formed during the last 50 ns of the MD simulations. As for the H-bonds with the solvent, for both isomers they are much larger in number (28  3

AC

on average for both tautomers) and stronger than the intramolecular ones. This identifies them as a major stabilizing factor of folate in aqueous solution. The analysis of the crystal structure of folic acid [36] is consistent with these results, since it is strongly stabilized by intermolecular hydrogen bonds with amino acids from the receptor pocket. Docking into the folate receptor- The two most stable tautomers FA-N3 and FA-N1 are flexibly docked into the binding site of the folate receptor-. Free rotation of the backbone single bonds of the folate is allowed, except for the amide bond, while the receptor structure is frozen.

CR

IP

T

ACCEPTED MANUSCRIPT

US

Figure 7: Superposition of the most favorable binding poses of FA-N3 (thick cyan lines) and FA-N1 (thick turquoise lines) in the X-ray structure of FR-α [36]; only the protein restudies forming the

AN

binding pocket of the receptor are shown as thin lines.

M

As can be seen from Figure 7, the best poses of the two studied tautomers are practically identical and resemble the experimentally determined position, in which folic acid binds to the

ED

receptor [36]. The most favorable orientation is when the pterin part of the molecules is deeply buried into the pocket of the active site of the macromolecule. FA-N3 and FA-N1 have similar

PT

affinity (-11.1 and -11.3 kcal/mol) to the protein. It can be concluded from these results that the

CE

position of the proton will not be discriminative for the interactions with the receptor.

NMR characterization of folate

AC

Comparison between the calculated and the experimental NMR shifts of folate in water (Tables 3 and S1, Figures S1 and S2) and of folic acid in DMSO (Tables S2 and S3, Figures S3 and S4) shows very good correspondence in the obtained values. A key factor for that is the averaging over all representative conformers, which became evident in the process of data analysis. All used computational methods represent properly the order of the signals both in 1H and 13C NMR spectra. The calculations with B3LYP tend to somewhat overestimate the chemical shifts (Table S1 and S2). The closest to the experimental shifts are those calculated with WP04 for 1H NMR and WC04 – for 13

C. The DFT calculations for the methylene group hydrogens next to the asymmetric carbon in the

glutamate residue are also in good agreement with the experimentally observed ones, which is an

ACCEPTED MANUSCRIPT additional verification of the accuracy of the calculations. The experimental results from the 1H and 13

C NMR spectra are in agreement with previously described chemical shifts [40, 41].

Table 3: Experimental and WP04/WC04 calculated 1H and 13C NMR shifts of folic acid in aqueous solution. Atom

1

1

H, exp,

H, FA-N1,

1

13

H, FA-N3,

C, exp,

13

C, FA-N1,

ppm

ppm

ppm

ppm

2

-

-

-

159.4

4

-

-

-

169.3

4a

-

-

-

6

-

-

7

8.65

8a

ppm

T

ppm

C, FA-N3,

156.7±0.2

172.4±0.1

165.7±0.2

130.2

127.8±0.3

124.8±0.3

-

151.7

151.0±0.9

148.7±1.0

8.76±0.01

9.04±0.05

151.4

145.9±1.7

150.8±0.9

-

-

-

156.5

149.4±0.3

160.9±0.2

9

4.46

4.68±0.10

4.65±0.07

US

№

13

48.5

46.8±1.5

48.0±0.6

11

-

-

-

153.5

155.1±0.6

155.3±0.5

12

6.66

6.69±0.10

6.45±0.29

115.1

107.1±0.6

106.9±0.4

13

7.63

8.10±0.16

8.02±0.13

131.8

133.5±0.3

133.6±0.3

14

-

-

-

124.5

119.8±0.7

119.5±0.7

15

-

-

-

172.2

168.8±0.6

169.7±0.5

17

4.33

4.20±0.11

4.00±0.23

58.8

55.5±0.6

57.6±3.3

18

-

-

-

182.0

180.8±0.3

182.0±1.8

19a

2.17

19b

2.04

31.4

32.4±3.0

31.8±1.0

20 21

IP

CR

AN

M

ED

PT

CE

158.3±0.1

2.06±0.08

1.87±0.16

1.87±0.17

2.32

2.29±0.06

2.25±0.15

37.1

37.7±2.4

38.2±1.0

-

-

-

185.1

183.7±0.9

183.7±0.8

AC

1.99±0.05

When looking at the experimental spectra in the two solvents, it can also be seen that folate in water and folic acid in DMSO have analogous spectral characteristics. This signifies that protonation of the two carboxylic groups has mild influence on the 1H and

13

C chemical shifts. The observed

difference between the calculated shifts for H18 and H21 is a result of the formation of intramolecular hydrogen bonds between the carboxylic groups and between them and the carbonyl and NH group from the amide bond – when such a H-bond is formed, the signal is shifted to lower

ACCEPTED MANUSCRIPT field. Because carboxylate hydrogens do not undergo fast dissociation in DMSO, they can be seen in the 1H spectra in DMSO as a much broadened peak at 12.1 ppm, which is another evidence for the formation of H-bonds. The absence of a signal for hydroxyl protons in the 1H spectra in both solvents, along with the chemical shift of carbon C4, suggest that no lactim tautomers are present, at least in quantities observable by NMR. Additionally, no direct indications for the presence of an amide bond with cis-configuration were found in both solvents. The latter could also be indirectly

T

confirmed by the absence of measurable heteronuclear coupling constants between the amide proton

IP

and carbon C14. In addition, the experimentally determined value for heteronuclear coupling constant between H16 (NH) and C14 of ~0 Hz suggests cis-orientation for the H16-N-C15-C14

CR

dihedral angle and, correspondingly, preferred trans-configuration of the amide bond. This is in line

high calculated energy barrier for rotation [24].

US

with the predicted energy difference between the two configurational isomers and with the relatively

Unfortunately, it cannot be concluded unequivocally also from the experimental data which

AN

lactam tautomer is predominant – FA-N1 or FA-N3. Both NMR experiments and quantum calculations confirm that in solution there is dynamic equilibrium between FA-N1 and FA-N3. The

M

fast exchange of protons attached to nitrogen at positions N1 and N3 is the reason for the very low intensity of some carbon signals in the pterin residue and for the same reason the exchangeable

ED

protons in positions N1 and N3 could not be seen in both solvents. The NMR shifts of almost all carbons in pterin obtained from the experiment are very close to the DFT calculated ones. The only

PT

exceptions are C4 and C8a, which are the most NMR-sensitive ones to the position of the proton in question. There, the experimental values have intermediate values between those calculated for FA-

AC

Conclusions

CE

N1 and FA-N3, which is another proof for the existing equilibrium.

DFT study of six tautomers of folate shows that the most stable tautomer is the FA-N3 lactam form. However, the energy difference between the two lactam forms FA-N3 and FA-N1 is small, FA-N3 being lower in energy by 2.44 kcal/mol. This energy ordering coincides with previous findings for the pterin molecule only. It signifies that the remaining more flexible part of folate does not change the energetically favored tautomer compared to the pterin molecule. Experimental 1H and 13C NMR spectra confirm that only lactam tautomers are present in water and DMSO solution. Atomistic MD simulations of the two lactams in saline do not outline striking differences. The

ACCEPTED MANUSCRIPT structure of FA-N3 is more flexible and this could be related to the ability of this lactam to interact with the active center of the folate receptor-. Docking of the two lactam tautomers in the pocket of the receptor reveals that the most preferred pose is not sensitive to the lactam form and both molecules feature similar interactions with the protein. All computational techniques and the NMR experiment used in this study suggest the existence of equilibrium between the two lactams FA-N3 and FA-N1 of folate in aqueous solution. Our results

T

clearly identify the lactam form as the most stable tautomer and this isomer is recommended as the

IP

focal structure for studying the mechanism of the highly specific interactions with the folate

CR

receptor-.

Acknowledgement

US

The research was funded by a project of the Bulgarian Scientific Research Fund, contract №

M

AN

DN09/14 from 16.12.2016.

Highlights

DFT study reveal the lactam forms of folate as the most stable tautomers.



Experimental 1H and

13

ED



PT

C NMR spectra confirm that only lactam tautomers are present in

water and DMSO solution.

Atomistic MD simulations in saline outline the folate tautomer with hydrogen attached to N3

CE



(FA-N3) as more flexible than FA-N1. Docking in the active site of the FR- reveals that the most preferred pose is not sensitive to

AC



the lactam form.

References

1

Caliceti P., Salmaso S., Semenzato A., Carofiglio T., Fornasier R., Fermeglia M., Pricl S.,

Synthesis and physicochemical characterization of folate-cyclodextrin bioconjugate for active drug delivery, Bioconjug. Chem. 2003, 14, 899-908.

ACCEPTED MANUSCRIPT 2

Yan, J. J., Liao, J. Z., Lin, J. S., He, X. X. Active radar guides missile to its target: receptor-based

targeted treatment of hepatocellular carcinoma by nanoparticulate systems. Tumor Biol. 2015, 36, 55-67. 3

Wu, Z., Li, X., Hou, C. Qian, Y. Solubility of Folic Acid in Water at pH Values between 0 and 7 at

Temperatures (298.15, 303.15, and 313.15) K, J. Chem. Eng. Data, 2010, 55, 3958-3961. 4

Soniat M., Martin C. B., Theoretical Study on the Relative Energies of Neutral Pterin Tautomers,

IP

Martin Yv., Let’s not forget tautomers, J Comput Aided Mol Des 2009, 23, 693–704.

6

Brzeziñska A., Wiñska P., Baliñska M., Cellular aspects of folate and antifolate membrane

transport, Acta Biochimica Polonica, 2000, 47, 735–749.

CR

5

T

Pteridines, 2008, 19, 120-124.

Ph. Eur. 9.0, 1742 (07/2010:0067, p.2524).

8

https://webbook.nist.gov/cgi/cbook.cgi?ID=C59303&Mask=80

9

http://jpdb.nihs.go.jp/jp/DetailList_en.aspx?submit=Detail(en)&keyword=Folic+Acid

AN

10

US

7

Pfleiderer W., Liedek E., Lohrmann R., Rukwied M., Pteridine X. Zur Struktur des Pterins, Chem.

11

M

Ber., 1960, 93, 2015-2024.

Brown D. J., Jacobsen N. W., Pteridine Studies. Part XIV. Methylation of 2-Amino- 4 -

12

ED

hydroxypteridine and Related Compounds, J. Chem. Soc. 1961, 4413-4420. Rokos H., Pfleiderer W., Pteridine, XL Zur Synthese von 3.N2. N2-Trimethyl-pterinen, Chem.

13

PT

Ber., 1971, 104, 739-747.

Gready J. E., Theoretical studies on Pteridines. 2. Geometries, Tautomer, Ionization and Reduction

CE

Energies of Substrate and Inhibitors of Dihydrofolate Reductase, Journal of Computational Chemistry, 1985, 6, 377-400.

Uchimaru T., Tsuzuki S., Tanabe K., Benkovic S. J., Furukawa K., Taira K., Computational

AC

14

studies on pterins and speculations on the mechanism of action of dihydrofolate reductase, Biochemical and Biophysical research communications, 1989, 161, 64-68. 15

Schwalbe C. H., Lowis D. R., Graham Richards W., Pterin 1H-3H Tautomerism and its Possible

Relevance to the Binding of Folate to Dihydrofolate Reductase, J. Chem. Soc., Chem. Commun., 1993, 1199-1200. 16

Soniat M., Martin C. B., Theoretical Study on the Relative Energies of Anionic Pterin Tautomers,

Pteridines, 2009, 20, 124-129. 17

Nekkanti S., Martin C. B., Theoretical study on the relative energies of cationic pterin tautomers,

Pteridines, 2015, 26, 13-22.

ACCEPTED MANUSCRIPT 18

Reibnegger G., An ab initio and density functional theory study on neutral pterin radicals,

Pteridines, 2015, 26, 135-142. 19

Abyar F., Novak I., Conformational analysis and electronic structure of folic acid: A theoretical

study. Journal of Molecular Liquids 2019, 276, 819–825. 20

Jorgensen W. L., Maxwell D. S., Tirado-Rives J., Development and Testing of the OPLS All-

Atom Force Field on Conformational Energetics and Properties of Organic Liquids, J. Am. Chem.

IP

21

T

Soc., 1996, 118, 11225-11236.

Jorgensen W. L., Chandasekhar J., Madura J. D., Impey R. W., Klein M. L., Comparison of

22

CR

Simple Potential Functions for Simulating Liquid Water, J. Chem. Phys., 1983, 79, 926-935. Daan Frenkel and Berend Smit, Understanding Molecular Simulation. From Algorithms to

23

HyperChem 7.0; Hypercube: Gainesville, FL, 2001.

Iliev S., Gocheva G., Ivanova N., Atanasova B., Petrova J., Madjarova G., Ivanova A.,

AN

24

US

Applications, Ed. Academic Press, 2002

Identification and computational characterization of isomers with cis and trans amide bonds in folate 25

M

and its analogues, Phys. Chem. Chem. Phys., 2018, 20, 28818-28831 Petrova J., Gocheva G., Ivanova N., Iliev S., Atanasova B., Madjarova G., Ivanova A., Molecular

ED

simulation of the structure of folate and antifolates at physiological conditions, Journal of Molecular Graphics and Modelling 2019, 87, 172-184

Glendening E. D., Reed A. E., Carpenter J. E., Reed A. E., Curtiss L. A., Weinhold F.,

PT

26

Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint. Chem. Rev., 27

CE

1988, 88, 899–926.

Wiitala K. W., Hoye T.R., Cramer Ch., Hybrid Density Functional Methods Empirically

AC

Optimized for the Computation of 13C and 1H Chemical Shifts in Chloroform Solution, J. Chem. Theory Comput. 2006, 2, 1085-1092 28

(a) Tomasi J., Mennucci B., Cammi R., Quantum mechanical continuum solvation models, Chem.

Rev. 2005, 105, 2999-3094; (b) Scalmani G., Frisch M. J., Continuous surface charge polarizable continuum models of solvation. I. General formalism, J. Chem. Phys. 2010, 132, 114110-114125. 29

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., et al. Gaussian 09, revision D.01; Gaussian, Inc., Wallingford CT, 2009 (full citation is given in the Supporting Information) 30

Lu, T., Chen, F. Multiwfn: A multifunctional wavefunction analyzer. J. Comput. Chem. 2012, 33,

580-592.

ACCEPTED MANUSCRIPT 31

Vanommeslaeghe K., Hatcher E., Acharya C., Kundu S., Zhong S., Shim J., Darian E., Guvench

O., Lopes P., Vorobyov I., MacKerell Jr. A. D. CHARMM general force field: a force field for druglike molecules compatible with the CHARMM all-atom additive biological force field, J. Comput. Chem. 2010, 31, 671-690. 32

Bussi, G., Donadio, D., Parrinello, M. Canonical sampling through velocity rescaling. J. Chem.

Phys. 2007, 126, 014101.

T

Berendsen, H. J. C., Postma, J. P. M., DiNola, A., Haak, J. R. Molecular Dynamics with Coupling

IP

33

to an External Bath. J. Chem. Phys. 1984, 81, 3684-3691.

Abraham M. J., Murtola T., Schulz R., Páll S., Smith J. C., Hess B., Lindahl E., GROMACS: High

CR

34

performance molecular simulations through multi-level parallelism from laptops to supercomputers. 35

US

SoftwareX 2015, 1, 19-25.

Humphrey, W., Dalke, A., Schulten, K. VMD - Visual Molecular Dynamics. J. Molec. Graph.

36

AN

1996, 14, 33-38.

Chen C., Ke J., Zhou X. E., Yi W., Brunzelle J. S., Li J., Yong E. L., Xu H. E., Melcher K.,

M

Structural basis for molecular recognition of folic acid by folate receptors, Nature 2013, 500, 486489.

Trott O., Olson A. J., AutoDock Vina: improving the speed and accuracy of docking with a new

ED

37

scoring function, efficient optimization and multithreading, J. Comput. Chem. 2010, 31 455–461. (a) Feixas, F., Matito, E., Poater, J., Solà, M. Quantifying aromaticity with electron delocalisation

PT

38

measures. Chem. Soc. Rev. 2015, 44, 6434-6451; Krygowski, T. M., Cyrański, M. K. Structural 39

CE

Aspects of Aromaticity. Chem. Rev. 2001, 101, 1385-1419. Szczepanik, D. W., Andrzejak, M., Dominikowska, J., Pawełek, B., Krygowski, T. M.,

AC

Szatylowicz, H., Solà, M. The electron density of delocalized bonds (EDDB) applied for quantifying aromaticity. Phys.Chem.Chem.Phys., 2017, 19, 28970. 40

Bonechi C., Donati A., Lampariello R., Martini S., Picchi M. P., Ricci M., Rossi C., Solution

structure of folic acid Molecular mechanics and NMR investigation, Spectrochimica Acta Part A 2004, 60 1411–1419. 41

Cheung H. T. A., Birdsall B., Frenkie T.A., Chau D. D., Feeney J., 13C NMR Determination of

the Tautomeric and Ionization States of Folate in Its Complexes with Lactobacillus casei Dihydrofolate Reductase?, J Biochemistry 1993, 32, 6846-6854.