Molecular Sciences Ab Initio Studies on the Preferred Site of Protonation in Cytisine in the Gas Phase and Water

Ab initio calculations (HF, MP2, DFT for isolated and PCM for solvated molecules) were performed for cytisine (1) and its model compounds: N-methyl-2-pyridone (2) and piperidine (3). Among three heteroatomic functions (carbonyl oxygen, pyridone and piperidine nitrogens) considered as the possible sites of protonation in 1, surprisingly the carbonyl oxygen takes preferentially the proton in the gas phase whereas in water the piperidine nitrogen is firstly protonated. For model compounds, the piperidine nitrogen in 3 is more basic than the carbonyl oxygen in 2 in both, the gas phase and water.


Introduction
Proton-transfer reactions play an important role in interactions of molecules with biological systems [1].Biologically active molecules are usually polyfunctional derivatives in which the sites of protonation or deprotonation strongly depend on acid-base properties of individual functional groups, their intra-and intermolecular interactions and environment [2,3].Examples are alkaloids, naturally synthesized by various plants, some of which have medicinal or poisonous properties [4].Although, they belong to the family of amines, most of them contain additionally other functions, which may compete with the amino group in the proton-transfer reactions.
Cytisine (1) -a toxic quinolizidine alkaloid (synthesized by various plants, e.g.Cytisus, Baptisia, Laburnum vulgare, Thermopsis rhombifolia, Baptisia australis) [4] -contains three heteroatoms, the pyridone nitrogen and oxygen linked in the amide function, and separately the piperidine nitrogen.Two of them, the pyridone oxygen and piperidine nitrogen can be considered as potential basic sites.Depending on environment one or other site may take the proton.From biological point of view, cytisine belongs to the family of agonists of nicotinic acetylcholine receptors (nAChRs) [5,6].It specifically binds to nAChRs localized within the central nervous system and the neuromuscular junctions [7].In large doses, cytisine can kill by respiratory failure [4].In small doses, cytisine gives positive results in treatment of schizophrenia, Parkinson's disease and Alzheimer's disease [8].
To understand interactions of cytisine with biological systems, it is very important to known structure of cytisine and its favoured site of protonation in different environments.Crystal structure of neutral cytisine has been determined by X-ray crystallography [9], and FT-IR and Raman spectroscopy [10].Two independent molecules of identical conformation -linked together by the intermolecular Hbond between the NH group of one molecule and the C=O group of the other molecule -have been identified.Due to low solubility of cytisine, such kind of H-bond has not been detected in CDCl 3 solutions using NMR technique [9a,11] as well as in eight solvents of different polarities, polarizabilities and acid-base properties, CCl 4 , CS 2 , CHCl 3 (and also CDCl 3 for comparison with the NMR results), CH 2 Cl 2 , MeOH, Et 2 O and Et 3 N using FT-IR spectroscopy [12].Only conformational isomers have been identified.The corresponding NMR and FT-IR signals have been attributed to two conformational isomers, which proportion slightly depends on the polarizability and H-bonding capacities of the solvent.Recently, we performed semiempirical calculations for all possible conformations of isolated cytisine to find its conformational preferences in the gas phase, and next, we used ab initio calculations for the most stable conformations obtained at the semiempirical level to estimate their physicochemical properties [10,13].The basicity parameters have not yet been investigated for cytisine in the gas phase.Only the pK a in water has been measured [14].In this paper, ab initio calculations were extended to investigate the protonation reaction in cytisine (1) and to indicate the favoured basic site.For calculations, two model compounds: N-methyl-2pyridone (2) and piperidine (3), which contain the same functional groups as 1, were also considered.Calculations were performed for the most stable conformations.We used three methods for the gas phase [15]: the RHF (restricted Hartree-Fock), MP2 (second-order Möller Plesset perturbation) and DFT (density functional theory) with the B3LYP functional (the combination of the Becke threeparameter hybrid exchange functional with nonlocal correlation functional of Lee, Yang and Parr), and the PCM model (polarizable continuum method) [16] with water selected as the solvent.The 6-31G* basis set [15a] was selected to calculate relative basicities [17].To verify the effect of the basis set on the relative basicities, calculations were also performed using the 6-311G** basis set.The results of calculations enable us to indicate the preferred site of protonation in the gas phase and in water for cytisine, and to compare the calculated basicity parameters of each potential basic site in 1 with those in model compounds 2 and 3.

Computational Details
The choice of computational method, which could be applied to study the basicity of any compound, is very difficult [18].Smith and Radom [19] recommended the G2 method as one of the best approaches to calculate thermodynamic basicity parameters (such as proton affinity, PA).This method is based on the MP2 geometries.Since the MP2 optimization is rather time-consuming, the G2 method has only been applied to relatively small molecules [18][19][20].In the case of large molecules, single-point calculations {DFT(B3LYP) and/or MP2} have often been applied to structures optymized at the HF/6-31G* [21] or HF/6-31G** level [22].Taking into account the number of heavy atoms in the cytisine molecule, the literature propositions and our experience in the computed PAs [17][18][19][20][21][22], we chose the RHF method and the 6-31G* basis set for geometry optimization.To verify the effect of the basis set, calculations were also performed using the 6-311G** basis set.
Energies of protonation in the gas phase (E prot/gas ) for each potential basic site of isolated cytisine (1) and its model compounds (2 and 3) were calculated for reaction (1) according to eq (2).The energy of isolated proton is equal to zero, and thus it is neglected in eq 2. The effect of solutesolvent interactions in water was studied using the PCM model and the geometries optimized at the HF/6-31G* level.In the PCM calculations the van der Waals atomic radii was employed.According to the PCM method, solute was placed in a cavity formed by a union of spheres centered on each atom.The spheres defining the cavity were 1.2 times the van der Waals radii.E prot/aq s in water were estimated for reaction (1) using eq (3).In this model, the energy of the proton was calculated as the difference between the energy of H 3 O + (-76.4232408728a.u.) and water (-76.0209107822a.u.).

B + H + BH +
(1) E prot/gas = E(BH + ) gas -E(B) gas (2) Proton affinities in the gas phase {PA -the negative of the direct enthalpy change of the protonation reaction (1)} for the most basic site in the most stable structures of model compounds (2b and 3a) were estimated using eq (4) as previously described [3c,17b,18].Eq (4) includes the changes in electronic energy, in zero-point vibrational energy (ZPVE), in vibrational energy on going from 0 to 298.15 K, in rotational and translational energy, and the work term {∆(pV) = -RT = -0.592kcal mol -1 at 298.15 K}.For the proton, only the translational energy term is not equal to zero (3/2RT = 0.889 kcal mol -1 at 298.15 K).Gas-phase basicities {GB -the negative of the direct Gibbs free energy of the protonation reaction (1), GB = -∆G 298 (1) = G 298 (B) + G 298 (H + ) -G 298 (BH + )}, were estimated using eq (5).They differ from the PA by the corresponding entropy term {T∆S 298 (1) = TS 298 (B) + TS 298 (H + ) -TS 298 (BH + ), where S 298 is the sum of the rotational, vibrational and translational entropies}.For the proton, only the translational entropy is not equal to zero {S transl (H + ) = 26.040cal mol -1 K -1 }.For small simple molecules, the entropy terms of the neutral and ionic forms are almost the same (within 1-2 kcal mol -1 ).Therefore, the GBs differ from the PAs mainly by the translational entropy term of the proton (ca 7.8 kcal mol -1 ).All thermal corrections were computed at the HF/6-31G* and HF/6-311G** levels.
Results and Discussion

Preferred Structures for Isolated and Solvated Neutral Molecules
Cytisine (1) contains two asymmetric carbon atoms (C 7 and C 9 ) and one amino nitrogen atom in the piperidine cycle, which determine the number of isomers.Semiempirical calculations performed for possible isomers of cytisine indicated that structures 1a and 1b (Fig. 1) are the most stable [13a].Isomers 1a and 1b possess the same configuration (RS) on the asymmetric carbons, and different orientation of the lone electron pair and the hydrogen at the piperidine nitrogen.Structure 1a was found for the solid cytisine by X-ray, FT-IR and Raman measurements [9,10], and both isomers 1a and 1b were identified for the solvated cytisine in FT-IR and NMR spectra [10][11][12].For cytisine model compounds: N-methyl-2-pyridone (2) and piperidine (3), two stable conformations for each derivative (2a, 2b, 3a and 3b given in Fig. 1) were found at the semiempirical level.
Extended ab initio calculations {(RHF, MP2, DFT(B3LYP) and PCM} performed for the stable structures of cytisine (1a and 1b) and its model compounds (2a, 2b, 3a and 3b) indicate that the total electronic energies of investigated compounds in conformation a are close to those of b (Table 1).Enlarging the basis set to 6-311G** has a small effect on the computed relative energies.For cytisine, zero-point vibrational energies (ZPVE) and thermal corrections to the Gibbs energy calculated at the HF/6-31G*//HF/6-31G* and HF/6-311G**//HF/6-311G** levels are almost the same for both conformers (a and b), and thus, the Gibbs energies differ by less than 1.5 kcal mol -1 (1 cal = 4.184 J).This confirms previous results [9a, [10][11][12][13]] that cytisine may exist as a mixture of two conformers in the gas phase or apolar environment.The use of PCM model with water as the solvent does not change this behaviour.Similar situation is observed for the model compounds.

Geometrical Parameters
As shown previously, geometrical parameters (bond length in Å and angles in degrees) computed for both isomers of cytisine (1a and 1b) at the HF/6-31G* level are close to those reported in the literature for crystal cytisine [13a].Enlarging the basis set to 6-311G** has no important effect on general behaviour, which can be summarized as follows.
Change of the hydrogen orientation from the equatorial in 1a to the axial in 1b has significant influence only on particular geometrical parameters.The largest differences are found in the piperidine cycle of 1: (i) the N 12 -C 11 -C 9 and N 12 -C 13 -C 7 angles are larger in 1b than 1a by 4.8 and 4.3°, respectively, whereas differences in other angles are smaller than 0.5°, and (ii) the C 7 -C 13 and C 9 -C 11 bonds are longer in 1b than in 1a by 0.012 and 0.007 Å, respectively, whereas differences in other bond lengths do not exceed 0.005 Å.These exceptions indicate that 1b may be stabilized by an intramolecular interaction of the axial hydrogen at the piperidine nitrogen with the pyridone moiety.The distance between the pyridone nitrogen and the axial hydrogen at the piperidine nitrogen in 1b is equal to 2.821 Å, which is typical for a weak intramolecular N-H⋅⋅⋅N interaction.
The distance between the potential basic sites, the pyridone oxygen (sp 2 ) and piperidine nitrogen (sp 3 ) in cytisine (equal to 4.881 and 4.957 Å for 1a and 1b in vacuo and to 4.9±0.1 Å in the crystal) is close to that between nitrogen atoms in nicotine (4.8±0.1 Å) [25].This indicates that the interatomic distance O(sp 2 )⋅⋅⋅N(sp 3 ) in cytisine may play the same role in binding with the nAChRs sites as the N(sp 2 )⋅⋅⋅N(sp 3 ) in nicotine [25].
Although three dimensional cytisine-nAChRs and nicotine-nAChRs structures are not yet known, except some model structures [26], one can suppose that both agonists bind with the receptor not only by the two heteroatoms (one in sp 2 and the other in sp 3 hybridization) placed in particular distance in the molecule, but also by the aromatic ring present in both derivatives.Aromatic character of πelectron system can quantitatively be measured by the HOMA index based on experimental or computed bond lengths [27].Defining the HOMA index, it has been assumed that bonds of the system with a fully aromatic character have their optimal lengths, and the HOMA is equal to unity (e.g. for benzene HOMA is close to 1).For non-aromatic system with localized π electrons, the HOMA is equal to zero (e.g. for 1,3,5-cyclohexatriene -a Kekulé structure for benzene -HOMA is close to zero).
The HOMA indices estimated for the RHF computed (0.57 for 1a and 0.55 for 1b) and experimental (0.77, 0.65 for X-ray 1a and 1a') cytisine [13a] can be compared with those for nicotine isomers (4a and 4b, which differ only by the pyridine nitrogen orientation, ∆E=0.53 kcal mol -1 ), calculated on the basis of geometries optimized at the HF/6-31G* (0.99 for both isomers), and structure found by gas electron diffraction (0.98 for 4a) [28].Comparison shows that the RHF method predicts lower HOMA value for the pyridone moiety in cytisine than for the pyridine ring in nicotine.For nicotine, the pyridine ring has an exceptional aromatic character (HOMA index close to unity).The protonation at the piperidine nitrogen strongly decreases the HOMA index for the pyridone ring in 1 (HOMA = 0.32 for geometries computed at the HF/6-31G* level).This suggests that the pyridone ring losses its aromaticity due to an engagement of a part of π electrons of the pyridone ring in interactions with the protonated piperidine nitrogen.Opposite effect induces the protonation at the pyridone oxygen in 1.It strongly increases the HOMA index of the pyridone ring (HOMA = 0.97 for geometries computed at the HF/6-31G* level).On the other hand, protonation of the N-aza in the pyridine ring and the N-amino in the pyrrolidine cycle in 4a only slightly reduces aromaticity of the pyridine ring (HOMA, computed for the HF/6-31G* geometries, varies from 0.99 for neutral nicotine to 0.97 and 0.98 for the N-aza and N-amino protonated forms, respectively).The differences in aromatic character of cytisine and nicotine and in changes of aromaticity on protonation can explain some differences in physiological and behavioural effects of both agonists.

Possible Basic Sites
Since cytisine (1) contains three heteroatoms (one carbonyl oxygen, and two nitrogens in the pyridone and piperidine cycle) all of them can be considered as possible basic sites (Fig. 2).The carbonyl oxygen has two lone electron pairs, one synperiplanar to the pyridone nitrogen and the other antiperiplanar.Each of these electron pairs can be used in a formation of the covalent bond with the proton.Therefore, two monocations protonated at the carbonyl oxygen (O sp and O ap ) are possible.The same situation takes place in model N-methyl-2-pyridone (2).The nitrogens has one lone electron pair and thus only one monocation can be formed in protonation reaction of the pyridone (N pyr ) and piperidine nitrogen (N pyr ) in 1, as well as in model compounds 2 and 3.
The pyridone nitrogen is linked with the carbonyl oxygen in the amide function.This direct binding changes drastically the basic properties of both sites, due to a possible n-π resonance conjugation (O=C−N< ↔ -O−C=N< + ).The conjugation effect reduces strongly basicity of the pyridone nitrogen and augments basicity of the oxygen.Therefore, the oxygen atom is usually the favoured site of protonation in compounds with the amide function.The same situation may take place in 1 and 2. This suggests that in cytisine two heteroatoms, the pyridone oxygen (sp 2 hybridized site) and the piperidine nitrogen (sp 3 hybridized site) may be the potential basic sites similarly as in nicotine, for which two heteroatoms, the pyridine nitrogen (sp 2 hybridized site) and the pyrrolidine nitrogen (sp 3 hybridized site) are the basic sites.

Favoured Site of Protonation in the Gas Phase
Ab initio calculations {RHF, MP2 and DFT(B3LYP)} were performed for all possible monoprotonated forms of cytisine, and energies of protonation calculated according to eq (2).In calculations, both stable conformations of cytisine (1a and 1b) have been considered.Table 2 summarizes the E prot/gas values found for cytisine in the gas phase for all possible basic sites: O sp , O ap , N pyr and N pip .Enlarging the basis set from 6-31G* to 6-311G** does not change general behaviour on the preferred site of protonation.It slightly increases the computed relative protonation energies (by 1-3 kcal mol -1 ).Independently on the calculation methods, the lowest total electronic energy values are found for the O ap -protonated forms in both conformations of cytisine.This indicates that the O ap site is the most basic in the gas phase, and it has the largest E prot/gas values.The E prot/gas values for the O sp site are lower by ca 4-5 kcal mol -1 .The N pyr site is the less basic one.Its E prot/gas values are lower than those of the O ap by more than 40 kcal mol -1 , and thus the N pyr can be neglected in the gas phase.The N pip site is less basic than the O ap by ca 6-10 and 1-4 kcal mol -1 in 1a and 1b, respectively.This means that the order of basic site preference for cytisine in the gas phase is as follows: O ap >N pip >>N pyr .
The order of basic site preference in 1 does not follow that in model compounds, N-methyl-2pyridone (2) and piperidine (3).In calculations, the two stable conformations (a and b) were considered for both model compounds, and the N pyr site was neglected in 2. The calculated E prot/gas values are listed in Table 3. Independently on the method of calculations, the N pip in piperidine is more basic site in the gas phase than the O ap and O sp in pyridone.The difference between the E prot/gas values of the N pip and O ap is ca 5-8 kcal mol -1 .Larger basis set (6-311G**) led to slightly lower ∆E prot/gas values (by 1-3 kcal mol -1 ).According to definitions, the ∆E prot/gas values, however, cannot be directly compared to that found on the basis of experimental proton affinities and gas-phase basicities of 2 and 3.For this reason, the PA and GB values were calculated for model compounds.[29] proton affinities (PA) a and gas-phase basicities (GB) a of model compounds (2 and 3) with those calculated for the most basic sites in the most stable isomers of 2b and 3a in the gas phase using eqs ( 4) and ( 5) The comparison confirms that (i) the nitrogen atom in piperidine (3) is more basic site in the gas phase than the oxygen atom in pyridone (2), (ii) enlarging the basis set from 6-31G* to 6-311G** only slightly decreases the computed relative basicity parameters (∆PA and ∆GB), and (iii) the 6-31G* basis set seems to be sufficient for calculation of the relative basicity parameters.The differences between the PA (5-6 kcal mol -1 ) and GB values (5-7 kcal mol -1 ) calculated for the N pip and O ap sites in model compounds using the 6-31G* basis set are in good agreement with the experimental relative basicity parameters (∆PA = 6.7 kcal mol -1 and ∆GB = 6.2 kcal mol -1 [29]).
Strong basicity of the carbonyl oxygen of 1 in the gas phase may result from the perpendicular orientation of the pyridone ring to the piperidine cycle.This orientation favours various effects.First, it increases the polarizability effect of the aliphatic bicyclic system linked to the pyridone moiety.Therefore, the carbonyl oxygen in 1 has stronger basic properties than in 2. Next, it augments the electron density in 1a due to possible interaction of the lone electron pair of the piperidine nitrogen with the π electrons of the pyridone ring.This effect explains higher basicity (by ca 5-6 kcal mol -1 ) of the carbonyl oxygen in 1a than 1b.In 1b, the pyridone nitrogen and/or π electrons can interact with the hydrogen atom of the piperidine nitrogen.Finally, the field/inductive effect of the pyridone ring decreases basicity (by ca 5-6 kcal mol -1 ) of the piperidine nitrogen in 1 as compared to 3.
Similar behaviour has been observed for nicotine [30].For model compounds, 3-alkylpyridines are weaker bases than N-methyl-2-phenylpyrrolidine, whereas for nicotine the pyridine nitrogen (sp 2 hybridized site) seems to be more basic than the pyrrolidine one (sp 3 hybridized site).Several reasons have been proposed to explain this phenomenon.The most important are the polarizability effect of the pyrrolidine cycle, which increases the basicity of the pyridine nitrogen, and the field/inductive effect of the pyridine ring, which decreases the basicity of the pyrrolidine nitrogen.

Basic Site Preference in Water
The PCM model used for geometries optimized at the HF/6-31G* level was chosen to study the preferred site of protonation for cytisine (1) in water and also for its model compounds (2 and 3).Calculations were performed for neutral and all possible monoprotonated forms in water, and energies of protonation estimated according to eq (3).In calculations, both stable conformations (a and b) have been considered.Table 5 summarizes the E prot/aq values found for three possible basic sites: O sp , O ap and N pip .The N pyr site as the less basic one has been neglected.In 1, its E prot/aq value is lower than that of the N pip by more than 45 kcal mol -1 .The lowest total energy value is found for the N pip -protonated form of cytisine in both conformations.This means that contrary to the gas phase the N pip (sp 3 hybridized site) is the most basic in water, and it has the largest E prot/aq value.The E prot/aq values for the O ap or O sp (sp 2 hybridized site) are lower by ca 10-12 and 14-16 kcal mol -1 , respectively.This order of basic site preference in 1 follows that in model compounds, N-methyl-2-pyridone (2) and piperidine (3).The N pip in piperidine is more basic site than the O ap and O sp in pyridone.The difference between the calculated E prot/aq values of the N pip in 3 and O ap in 2 (19 kcal mol -1 ) is slightly larger than for cytisine.Indeed, in water solution, pyridone is a very weak base (pK a <1) in comparison to piperidine (pK a =11.22) [31].The difference between the experimental pK a values (> 14 kcal mol -1 in energy units) is close to that between the PCM calculated E prot/aq .The experimental pK a value (7.92 [14]) for cytisine indicates that the N pip is the favoured site of protonation in water.
Similar behaviour has been observed for nicotine [30].For both, nicotine and model compounds, the pyrrolidine nitrogen (sp 3 hybridized site) is more basic than the pyridine nitrogen (sp 2 hybridized site).This similarity between cytisine and nicotine in the distance between the basic sites, their preference in the protonation reaction, and its change when going from the gas phase to water can explain similarities in physiological and behavioural effects of both agonists of nAChRs.

Table 1 .
Differences in total electronic energies {∆E=E(a)-E(b)} a and Gibbs energies {∆G=G(a)-G(b)} a between two stable conformations a and b for cytisine (1), and its model compounds: N-methyl-2-pyridone (2) and piperidine (3) in the gas phase and water

Table 2 .
Energies of protonation (E prot/gas ) a calculated for heteroatoms in cytisine isomers (1a and 1b) in the gas phase using eq (2)

Table 3 .
Energies of protonation (E prot/gas ) a calculated for heteroatoms in model

Table 3 .
Cont. ∆E prot/gas = E prot/gas (N pip ) -E prot/gas (O ap ).c Geometries optymized at the HF/6-31G* level.dGeometriesoptymized at the HF/6-311G** level.Table4summarizes the PA and GB values calculated for the most basic sites (O ap and N pip ) in the most stable isomers of model compounds (2b and 3a).For comparison, their experimental PA and GB values are also given in this Table.

Table 4 .
Comparison of the experimental