Theoretical study of mechanisms for the hydrolytic deamination of cytosine via steered molecular dynamic simulations.

Gibbs free energy profiles of the cytosine deamination assisted by a water molecule in a discrete aqueous medium were obtained by the application of Steered Molecular Dynamic (SMD) simulations. Two pathways were considered to explain the mechanism of this process, where the water molecule attacks the C-N bond to give an intermediate (an amino-hydroxy-oxo structure in the A-path, and a hydroxy-oxo in the B-path) as the determinant step of reaction. Stationary structures along both energy profiles were analyzed at molecular dynamics level, obtaining states with higher free energies than those from electronic calculations in the gas phase and in solution described as a continuous medium. From the results obtained, the more complex A-pathway, with five steps, was kinetically the most favorable (with an endergonic reaction energy of 7.41 kcal mol-1, a high barrier of 67.53 kcal mol-1, and a small velocity constant k 2 = 1.80 × 10-37 s-1), concluding that the uracil base can participate in a spontaneous genetic mutation since the uracil-ammonia complex has a long lifetime of 6.10 × 1027 s. This process turns out exergonic and faster when carried out in gas phase simulation or electronic calculation with a continuous medium, due to the disappearance of explicit water molecules that can compete with the assistant molecule. This journal is © The Royal Society of Chemistry.

Introduction

Theoretical structural and thermodynamic studies are essential to understand the behavior of a variety of biochemical processes.[1-8] These studies are especially important for DNA and RNA acids, as these molecules are responsible for encoding and transmitting the genetic information necessary for proper cellular function, and for converting this information into functional products such as proteins. Despite protection from the cellular medium, these nucleic acids undergo alterations to their structure from attacks by various agents, both physical (ionizing radiation) and chemical (free radicals), which can lead to genetic mutations.[9] It is important to note that the nitrogenous bases of these acids are exposed to the physiological medium formed mainly by water molecules. Therefore, the canonical or non-canonical forms of these bases will depend on the surrounding environment in which they are, this medium being capable of modifying the stability and equilibrium constants of the tautomeric processes.[10-17] It has been shown that other effects, such as the base-pair stacking in DNA and RNA acid strands, have relevance in proton-transfer reactions.[18,19]

While DNA is a reasonably regular structure formed by two nucleotide chains wrapped around an axis which are held together by hydrogen bonds formed between the base pairs guanine–cytosine (G–C) and adenine–thymine (A–T) (so-called Watson–Crick pairs), RNA is formed by a simple chain of nucleotides similar to that of DNA, except that thymine is replaced by uracil (U). However, modified nucleobases in DNA and RNA are also known to occur, forming rare bases or base pairs that can participate in genetic mutation.

Cytosine is the most unstable of the DNA bases, interacting with other nitrogenous bases via hydrogen bonds, mostly with guanine to give the guanine–cytosine base pair. This DNA base has been extensively studied experimentally and computationally, in both gas and aqueous phase. Its most stable (canonical) form has an amino–keto structure, although there may be other less frequently observed structures, such as the imino–enol structure resulting from the exchange of hydrogen between the amine and keto groups. Its interaction with water has also been the subject of numerous studies due to its conversion into a uracil base and its participation in the genetic mutation. The cytosine(keto)–water complex is more stable than the cytosine(enol)–water complex by 1 kcal mol−1,[20] with an interconversion energy barrier between 30 and 37 kcal mol−1.[21]

As shown in Fig. 1, the hydrolysis of cytosine (C) can lead to a deamination obtaining the uracil base (U) and ammonia (A), i.e. the loss of an amino functional group with conversion to a carbonyl functional group. This process is important in spontaneous genetic mutation since uracil is potentially mutagenic and changes the coding information during DNA replication and RNA transcription. Although this deamination is a rare process under normal physiological conditions, its rate can significantly change in the presence of some reactants (OH−, H+, H2O, NO, HNO2, HSO3− …) that act as catalysts. The reactants are able to reduce the energy barrier of deamination by up to 28 kcal mol−1, with rate constants between 10−10 and 10−13 s−1.[22,23]

Fig. 1

Cytosine–water (C–W) to uracil–ammonia (U–A) conversion.

The C–W ⇆ U–A equilibrium produced by the cytosine deamination is the focus of this work. Specifically, the transference of hydrogen atoms from the water molecule to the nitrogen atoms of cytosine and the transference of the oxygen atom from water to the carbon atom linked to the amine group of cytosine, forming the corresponding uracil–ammonia (U–A) complex. It is known that interconversion between cytosine and uracil bases depends on the surrounding aqueous medium, hence the interest of studies that include the solvent effect. On one side, solvation methods[24-26] are the most practical way of taking into account solvent effects when calculating free energies in solution using electronic structure methods. Classical molecular dynamic simulations[27] are also a desirable alternative for predicting the evolution of the tautomeric process over time and for visualizing the solute–solvent system at the molecular level. The steered molecular dynamic (SMD) technique,[28,29] which applies external steering forces in a particular direction, allows us in a simpler way to determine the changes in Gibbs free energy over the course of any elementary process, and is useful for the study of reaction mechanisms.

There have been some theoretical studies on cytosine–water interactions.[20,21,30-43] Most of these studies focus on the gas phase using ab initio methods, although some authors employed a polarizable continuum model (PCM) or a mechanic molecular (MM) to describe the solvent. We can say that this is the first study in which the SMD technique has been applied to calculate the thermodynamic and kinetic properties of the deamination of cytosine.

The first computational study of cytosine deamination with the OH− radical was studied by Šponer et al.[39] using the PCM model to describe the solvent, obtaining an activation barrier of 51.0 kcal mol−1 that is considerably higher than the experimental value.

Almatarneh et al.[40] performed an ab initio study in the gas phase (at HF, MP2, and B3LYP calculation levels) using two mechanisms for the cytosine deamination process with H2O and OH− as reactants, optimizing the stationary points of the reaction profile and evaluating thermodynamic magnitudes associated with this process. For hydrolysis with H2O, they obtained high values for free energy barriers ranging from 53.83 kcal mol−1 (for the five-step pathway initiated by the protonation of the imino nitrogen) up to 63.89 kcal mol−1 (for the simpler two-step pathway initiated by the addition of water oxygen to the carbon of the C–N bond). This barrier was reduced considerably up to 37.30 kcal mol−1 when the reaction was catalyzed by the OH− ion to initially deprotonate the amino group of the cytosine and subsequently add water. In both mechanisms, the determinant step of the reaction was the attack of water oxygen at the carbon of the cytosine to give a tetrahedral intermediate, and the reaction energy of the process was −8.18 kcal mol−1.

More recently, Labet et al.[41] performed a study of this deamination in gas and in solution phase, but using two water molecules as reactants, one acting as a nucleophilic agent and the other as assistant to the first water molecule. They considered two mechanisms: (a) adding water to the protonated cytosine through the single C–NH2 bond and (b) adding water to the neutral cytosine through the double CN bond. They found that the first pathway to be more favorable with a free energy barrier of 33.5 kcal mol−1 for the determinant step produced in the addition of the water molecule (relatively close to the experimental value and lower than that obtained by Almatarneh[40]). The results show how the aqueous medium decreases the barrier and stabilizes the intermediates with respect to the reaction in the gas phase.

Labet et al. expanded the study of the gas phase deamination from the perspective of the reaction force.[42] The energy barrier results were ΔG‡ = 51.2 kcal mol−1 and ΔG‡ = 38.7 kcal mol−1 using two and one water molecule, respectively. They concluded that the reaction force profile involved a concerted mechanism composed of two asynchronous steps: the protonation of the imino nitrogen and the attack of the water molecule on the carbon atom linked to the amino group, favoring the assistance of two water molecules. The reaction energies were −10.1 kcal mol−1 and −13.4 kcal mol−1 depending on the mechanism used.

The most recent study was realized by Zhang et al.[43] using a combined QM(DFT)/MM molecular dynamics simulations. The cytosine to uracil conversion was catalyzed by yeast cytosine deaminase (yCD), where the cytosine reacts with the water molecule coordinate to zinc and the adjacent Glu64 specie serves as a general acid/base to proton transfer from water to cytosine. The overall reaction consists of several proton transfer processes and nucleophilic attacks, forming a tetrahedral intermediate adduct of cytosine and water binding to zinc. They obtained a free energy barrier for the rate determinant step of 18.0 kcal mol−1.

Our main objective was the theoretical study of cytosine deamination from its canonical structure assisted by one water molecule. The calculation was performed in aqueous solution using SMD simulations to analyze the response of this methodology in biological process and to show the effect that a discrete medium has on atoms transfer. In addition, the conversion of C–W to U–A following two different mechanisms (similar to those studied by Almatarneh[40]) was considered, with the aim of finding the most favorable process, thermodynamically and kinetically. Simulations in the gas phase and electronic calculation in a continuous medium were also made to see the influence of the level of calculation in this process of deamination. Furthermore, a structural and energetic analysis of the stationary points in the different reaction mechanisms was carried out. Finally, an analysis of the energy profiles, rate constants, and lifetimes was realized, with the aim of predicting the possible participation of the uracil base in genetic mutation processes.

Computational details

Electronic structure calculations

Structure of the C–W complex was fully optimized with the M06-2X[44] functional and the 6-311++G(d,p) basis set, including vibrational analysis, the rest of stationary geometries were taken from the work of Almatarneh.[40] This particular meta-hybrid functional was chosen based on its adequate description of non-covalent intermolecular interactions in biological systems.[45-48] Solvent effects were taken into account by means of the IEF-PCM method[49] for geometry optimization and frequency calculations. Electronic structure calculations were performed using the Gaussian09 package.[50]

SMD simulations

Simulations using classical molecular dynamics allow the adequate representation of chemical processes in solution at a molecular level. In this work the solute–solvent interaction has been represented by a Lennard-Jones (12-6-1) functionwhere A and B are van der Waals parameters taken from the AMBER force field,[51]qs the net charge on each solute atom and qw the charge on the solvent atoms. Initially, qs was assigned using the Mulliken procedure[52] at the above mentioned level of calculation, but during the simulation process solute charges were recalculated at every step by the RESP method[53] implemented in the software used in simulations. On the other side, qw were pre-assigned as the TIP3P charges,[54-57]i.e., δ(O) = −0.83 and δ(H) = 0.415.

More than 2000 water molecules were included in the SMD simulations to represent an aqueous environment around the C–W complex. All simulations were performed using AMBER software[58] at 298.15 K. The initial C–W geometry was obtained from electronic structure calculations by fully optimizing the complete system, including the water molecule that assist the deamination process. Bond length, bond angle, torsional, and van der Waals parameters were taken form Kollman et al.[51]. Water molecules located initially at distances less than 1.6 Å from any solute atom were eliminated from the simulations. The long-range electrostatic interactions were treated by the Ewald method.[59] A cut-off of 7 Å was applied to the water–water and C–W–solvent interactions in order to simplify the calculations, and periodic boundary conditions were used to keep the number of solvent molecules constant in an NVT ensemble. A 300 K temperature and the Andersen temperature coupling scheme was employed.

The free-energy change along the reaction coordinate was represented by the potential of mean force (PMF). To relate Gibbs energy differences between two equilibrium states, the Jarzynski's equality[60] was used. So, this free energy change (ΔG) is related to an average over all possible works 〈W〉 of an external process that takes the system from the A equilibrium state to a new B state, as shown by eqn 2where W was calculated by integrating the force over the distance of the steered atom from its initial to its final position, according to eqn 3

SMD technique allows the calculation of Gibbs energy changes during the course of an elementary process in a simpler way, by applying external forces in a particular direction. SMD simulations were carried out using the QM/MM method with the semi-empirical Hamiltonian SCC-DFTB,[61,62] implemented in the AMBER12 software by Case et al.,[63,64] despite constraints and limited flexibility of this approximated quantum chemical method in its application to biological systems.[65] The system was partitioned during the simulations, applying quantum calculations to the solute (QM subsystem) while the aqueous solvent molecules were described in a classical way (MM subsystem). Specifically, the cytosine base and the water molecule that assists in the proton transfer process have been the only molecules included in the QM subsystem.

The standard Gibbs energies of activation (ΔG‡) and reaction (ΔG) can be determined by following the time evolution of the process, as previously shown in several works,[66-71]i.e., from the configuration of the reactant to the product (or vice versa) in elementary steps focusing on the reaction coordinate. In a previous work,[66] several calculations at different post-equilibrium simulation times (200, 400, 600, 800, and 1000 ps) were performed to verify the Gaussian distribution of W in eqn (2). Based on those results, long simulations have been used for each step of the deamination process in this work, verifying in each simulation that the variations of the chosen reaction coordinates do not deviate of the minimum of the harmonic restraint. The structures involved in the reaction mechanism can be visualized to represent the evolution of the system at each step of the process, especially intermediate and transition states.

The system was initially minimized with a 1 ps simulation time and with 1 fs steps. Subsequently, a 100 ps equilibration phase was performed to obtain the initial configuration to begin the propagation and data storage phase, verifying that the energy balance was reached. Finally, during the final 100 ps of the simulation, the position and velocity of each atom was stored every 100 fs. That is, the simulation were performed with a simulation speed that varies in between 10 to 30 Å ns−1, depending on the reaction coordinate considered. The force constant used for the reaction coordinates was 1000 kcal mol−1 Å−2 for distances and 10 kcal mol−1 rad−2 for the angles, throughout the simulation. The reaction coordinates chosen for describing the proton transfer processes are discussed in the next section.

Results and discussion

Process and mechanisms

The C–W ⇆ U–A equilibrium has been studied by applying SMD simulations in which the aqueous solvent was treated as a discrete medium. The structure of the C–W system and the denomination used to identify their atoms are shown in Fig. 2. The initial orientation and separation between the assistant water molecule and the cytosine base come from the optimal geometry obtained by the electronic structure calculation using a M06-2X/6-311++G(d,p) level and the IEF-PCM polarizable continuous model. Cartesian coordinates of the C–W geometry are reported in Table S1 of the ESI†

Fig. 2

Structure and annotation of atoms in the C–W complex.

All the reaction coordinates (RC) considered to carry out the SMD simulations for mechanisms studied can be expressed as:RCRCRCRCRCRCRCRCRCRCRCRCRC

Two mechanisms were considered in this work depending whether it starts with a protonation on the imino nitrogen N10 or with the attack of the water molecule on the C6–N7 bond of the cytosine base:

• A-pathway: initially the oxygen OW and hydrogen H2W atoms of the water molecule attack the hydrogen atom H9 and the nitrogen atom N10 of the cytosine, respectively, in a concerted way via RC1a1 and RC2a1 reaction coordinates, forming an imine–oxo intermediate (I1a). In a second step, the water molecule attacks again the cytosine base but now the oxygen atom OW on the carbon atom C6 and the hydrogen atom H9 on the nitrogen atom N7, by means of RC1a2 and RC2a2 reaction coordinates, forming an amino–hydroxy–oxo intermediate (I2a). In this second step, torsions of –OH and –NH2 groups adopt the optimal conformation of the intermediate to perform the protonation in the next step, via the reaction coordinates RC3a2 and RC4a2. Finally, in the third step, the intermediate I2a converts to the uracil–ammonia (U–A) complex as a result of the break in the C6–N7 bond and the transfer of the hydroxyl hydrogen H2W to the nitrogen N7 atom, through RC1a3 and RC2a3 reaction coordinates. All these steps and structures are schematized in Fig. 3.

Fig. 3

Stationary structures in both pathways.

• B-pathway: initially the assistant water molecule is displaced from its planar position (θO = 0°) to a more adequate perpendicular position (θO = 90°) to attack the cytosine base on the C–N bond. Then, the OW and H1W atoms of the water molecule concertedly attack carbon C6 and nitrogen N7 atoms of the cytosine base, respectively, following RC1b1 and RC2b1 reaction coordinates. Simultaneously, the torsion of the –NH2 group is realized and the C6–N7 bond starts breaking to give the hydroxy–oxo intermediate (I1b) and the ammonia molecule, via RC3b1 and RC4b1 reaction coordinates. In a second step, the hydroxyl hydrogen H2W atom is transferred intramolecularly to nitrogen N10 atom, through the RC1b2 coordinate, to form the uracil base, the ammonia molecule moves away from the uracil base. These two steps and structures are also schematized in Fig. 3.

In the simulations of both mechanisms, the free movement of all solvent molecules, the assistant water molecule, and some of the atoms of the cytosine base (specifically, C6, N7, H8, H9 and N10 atoms) were allowed. The remaining atoms were fixed during simulations to avoid distorting of the cytosine base. Cartesian coordinates of stationary structures in each mechanism from SMD simulations are collated in Table S2 of the ESI.†

Structural and energies from electronic calculations

Analysis of the structural and energetic results for the cytosine deamination obtained from single point calculation using a PCM-M06-2X/6-311++G(d,p) level, taking the geometries of Almatarneh,[40] gives the following observations:

(a) Describing the A-pathway: in the TS1a state, the H2W atom is at 1.41 Å from the donor atom OW and at 1.13 Å from the acceptor atom N10 (with an angle of θOHN = 149.9°), whereas the H9 hydrogen is at 1.13 Å from the donor nitrogen and at 1.39 Å from the acceptor oxygen, forming a six-membered ring. The I1a intermediate has the hydrogen linked to the N10 nitrogen, dN = 1.0 Å, forming an imine–oxo tautomer of the cytosine. In the second transition state TS2a, OW atom is at 1.76 Å from C6 carbon and the hydrogen at 1.26 Å from N7 nitrogen, forming a four-membered ring in a perpendicular plane to the cytosine ring. The second intermediate I2a shows a tetrahedral structure on carbon C6, forming a bond with the hydroxyl group at a distance dC = 1.46 Å, however, the orientation of the hydroxyl group –OH is not adequate to transfer the hydrogen H1W to the amino group –NH2, and therefore both groups perform rotations to reach the third transition state TS3a with the hydrogen H1W at 1.2 Å distance from nitrogen N7, whereas C6–N7 bond was extended to 1.6 Å.

(b) For the B-pathway case, TS1b transition state shows a has similar structure to the TS2a described above (except that the N10 atom is not linked to any hydrogen) and distances are dH = 1.09 Å, dH = 1.46 Å, dO = 1.92 Å, and dC = 1.48 Å. The water oxygen is then linked to C6 carbon, dO = 1.32 Å, in the hydroxyl–oxo intermediate Ilb, resulting in a uracil tautomer. Finally, in the second transition state TS2b, the hydroxyl hydrogen is at 1.35 Å from the oxygen donor and at 1.47 Å from the acceptor nitrogen, forming a four-membered ring in the same plane as the uracil ring.

(c) The free energies, including thermal corrections, obtained from quantum calculations with a continuous description of the surrounding medium (see Fig. 4 and Table S3 of the ESI†) have, in the A-pathway, a barrier of ΔG‡(TS1a) = 15.42 kcal mol−1 for the imine nitrogen protonation, another of ΔG‡(TS2a) = 50.39 kcal mol−1 for the attack of water on the CN double bond, and the last of ΔG‡(TS3a) = 46.75 kcal mol−1 for the hydroxyl hydrogen transfer to amine nitrogen N7, respect to the C–W complex energy considered as zero energy. The barriers in the B-pathway have values of ΔG‡(TS1b) = 55.88 kcal mol−1 for the tetrahedral intermediate formation and of ΔG‡(TS2b) = 31.50 kcal mol−1 for the hydroxyl hydrogen transfer to imine nitrogen N10. Therefore, we can say that the highest barrier in both mechanisms is obtained when atom OW attacked carbon C6 to form the tetrahedral structures on C6 atom. The greater stability of the U–A complex with respect to the C–W complex indicates the slight exothermic character of this deamination (ΔG = −3.58 kcal mol−1), whereas the higher barrier in the B-pathway makes this mechanism slower.

Fig. 4

Electronic energy levels for A-pathway and B-pathway in solution phase.

Energies from SMD simulations

Activation (ΔG‡) and reaction (ΔG) free energies for A- and B-pathways from the SMD simulations in solution and gas phase, together with the structure electronic results in gas phase of Almatarneh et al.,[40] are presented in Table 1.

Activation and reaction energies from SMD simulations

Mechanisms	Phase	ΔG‡a	ΔGb
A-pathway (1-step)	Solutionc	32.28	10.73
	Gasc	29.62	4.50
	Gasd	19.04	2.27
A-pathway (2-step)	Solutionc	67.53	28.78
	Gasc	60.10	18.54
	Gasd	53.82	13.66
A-pathway (3-step)	Solutionc	62.74	7.41e
	Gasc	53.65	−3.64
	Gasd	36.18	−8.18e
B-pathway (1-step)	Solutionc	87.78	20.28
	Gasc	73.32	7.85
	Gasd	63.89	2.25
B-pathway (2-step)	Solutionc	54.85	4.53
	Gasc	46.24	4.11
	Gasd	30.07	−8.18e

Activation energy (kcal mol−1) evaluated as ΔG‡ = G(transition state) − G(reactant).

Reaction energy (kcal mol−1) evaluated as ΔG = G(product, intermediate) − G(reactant).

In this work.

Electronic energies of Almatarneh et al.[40] using a G3MP2 level.

This electronic energy does not correspond to the U–A conformation of our work.

In general, simulations in an aqueous medium for both cytosine deamination pathways describe energy profiles that are in agreement with those obtained from our structure electronic calculations (compare Fig. 4 and 5), although higher energies are observed in all stationary states making this process slower. The small energy difference for the U–A product according to the mechanism followed is due to the final reorientation of the ammonia molecule with respect to the uracil molecule and the different solvent structure around the solute in each case. That is, energy profiles with more energetic transition states and intermediates are obtained using SMD simulations in solution with respect to our simulations in the gas phase and with those electronic calculations realized by Almatarneh et al.[40]

Fig. 5

Free energy profiles from the SMD simulations.

At this point, it is important to analyze why the activation barriers of the different mechanisms obtained by simulation are superior to those derived from our quantum studies. There are different causes that justify the observed discrepancies. First, the different levels of calculation used to treat the complex structures in the electronic calculations (PCM-M06-2X/6-311++G(d,p)) and in the QM-MM simulations (SCC-DFTB), which presents a different representation of the electron density. Another reason comes from the description at the molecular level of the solvent that can cause some molecule of the surrounding environment to compete with the assistant water molecule, making more difficult the transference and therefore a higher the activation barrier compared to our quantum studies, where a continuous description of the solvent is employed. Finally, the profile represents the energy at each step of the chosen reaction coordinate and showing the adaptation of the system to the current position of the fixed atoms of cysteine (partially restrained flexibility) and to the parameters used in the simulations, unlike the electronic calculation in which each stationary point is optimized.

It is also important to highlight why the activation barriers are higher in solution that in the gas phase using SMD simulations. The main cause that justifies this behavior is in water molecules that surrounded the C–W system, which can compete with the assistant water molecule and hinder the atoms transfer from one molecule to another. A second water molecule assisting the process avoid this effect and decrease the energies of the reaction profiles as we have verified in recent studies,[70,71] where similar processes were assisted by several water molecules, mainly in those representative structures of transition states where it is possible that six-membered rings were formed which makes more stable stationary states.

A-pathway simulations

In A-pathway, the free energy barrier obtained for the first step is ΔGTS1a‡ = 32.28 kcal mol−1, increasing to ΔGTS2a‡ = 67.53 kcal mol−1 in the second step, and to ΔGTS3a‡ = 62.74 kcal mol−1 in the third one, with respect to the C–W complex as cero energy (see Fig. 5). We can observe that this energy profile is in concordance with the one obtained from quantum calculations, where the higher barrier belongs to the attack of water molecule on carbon C6 to form the tetrahedral intermediate. The U–A final product is also at a higher energy than in the electronic case given the difference in calculation level and the description of solvent, as well as the different final configuration in each case.

In spite of the endergonic reaction energy, ΔG = 7.41 kcal mol−1, the high energy barriers show the difficulty in reverting U–A complex to the I2a intermediate and to the initial C–W complex. The activation energy for the determinant step is significantly different from the experimental value, showing that this reaction rarely occurs in a neutral medium and requires another type of catalyst or the presence of a second assistant molecule to decrease the barrier. However, A-pathway adequately describes the deamination mechanism of cytosine.

It is interesting to note that in A-pathway, the intermediate I1a represents a stable tautomer of the canonical cytosine and between the second and third step there are several conformational changes (torsions of the –OH and –NH2 groups) with indistinguishable barriers in SMD calculations (see energy profile between TS2a and TS3a). The final product U–A has a long lifetime of τU–A = 6.10 × 1027 s (calculated as τ = 1/k−3, being k−3 the reverse rate constant of the equilibrium I2a ⇆ U–A), indicating its possible participation in spontaneous mutation process. Besides, the deamination process was slow with a small rate constant of the determinant step of k2 = 1.80 × 10−3 s−1, evaluated as

If we compare the energy profile with the one obtained in the gas phase (see values in Table 1), we see a similar behavior although all the structures are more energetic. The presence of discrete molecules of water can lead to some of them entering in competition with the assistant molecule, hindering the transfer processes. In turn, in profile gas phase, highlights the exergonic character (ΔG = −3.64 kcal mol−1) in better agreement with the one presented by Almatarneh et al.[40] in a similar study. It should be also noted that the U–A product in solution has a configuration different from that of Almatarneh's work because in our simulation the ammonia molecule is away from the uracil base along the direction of the C6–N7 bond without seeking an optimum position that could be found in the O12–H13 zone by forming hydrogen bonds.

B-pathway simulations

In B-pathway, the energy barrier for the first-step is higher ΔG‡TS1b = 87.78 kcal mol−1, decreasing to ΔG‡TS2b = 54.85 kcal mol−1 in the second step. Again, the activation energy for the determinant step is far from the experimental value. The first step of the mechanism requires an initial energetic cost to move the water molecule from its original position to the perpendicular plane of the cytosine base (first zone of energy profile until 10 kcal mol−1) to begin the attack of OW on carbon C6, this energy barrier is higher than the equivalent step in A-pathway (about 60 kcal mol−1 from I1a intermediate up to TS2a). Apart from this displacement, other movements are also considered in this first step, such as the nucleophilic attack of the OW on the C6 atom, the proton transfer to the amino group, the break of the C–N bond, and the torsion of the O–H bond, what justify the height of this barrier. The second step, corresponding with the proton transfer to the amino nitrogen N7, has a similar energetic cost to that of proton transfer to N10 atom in A-pathway (about 34 kcal mol−1 in both cases).

The overall reaction energy of ΔG = 4.53 kcal mol−1 indicates an endergonic process with difficulty in reverting to other previous structures since an initial barrier of 50.32 kcal mol−1 has to be overcome. The final U–A product has a long lifetime of 1.30 × 1024 s, indicating its possible participation in spontaneous mutation processes. Furthermore, the deamination process is slower in this pathway with a smaller rate constant of k1 = 2.50 × 10−52 s−1 for the determinant step.

Again the structures in solution turn out to be more energetic than those in the gas phase, highlighting the similarity of energy for the final specie U–A in both simulations.

Comparison of pathways

Comparing both mechanisms (Fig. 5), the barrier that describes the determinant step of the attack of OW on carbon C6, is less energetic when realized from tautomer I1a (A-pathway) than from the cytosine canonical structure (B-pathway), the difference lays on in the higher charge density on the CN double bond in the I1a structure. The step that describes the proton transfer from the –OH to –NH2 group in the A-pathway is of similar energy as when the hydrogen is transferred to the –NH group in B-pathway. Likewise, the tautomeric intermediate of uracil I1b is more stable than the cytosine tautomer I1a since reverting to the initial C–W complex required more energy. Finally, the simulation of both pathways allows us to see the solvent structure around the C–W and U–A complexes, obtaining values of the global reaction energy that differ by less than 3 kcal mol−1.

Conclusions

The cytosine deamination process in aqueous solution has been studied by applying SMD simulations. In between different mechanisms used, the simulation of A-pathway, where a proton transfer is made from the water molecule in the first step, is the most kinetically favorable. All transition state and intermediate structures have energies higher than those obtained from electronic calculations using a continuous solvent model and from the gas phase simulations. Better results can be obtained when it is considered more water molecules assisting the deamination process. The formation of the tetrahedral intermediate on C6 is the most energetic step in both mechanisms, while the transfer of the hydrogen of the hydroxyl group to the amino nitrogen N7 or imino nitrogen N10 are of similar energy. The stability of the U–A complex, with long lifetimes, indicates that uracil bases formed in this deamination reaction could participate in genetic mutation. Finally, we would like to emphasize the importance of performing SMD simulations to describe reaction mechanisms in nitrogenous bases at a molecular level, from which the evolution of the C–W ⇆ U–A process along the reaction coordinate can be studied and visualized.

Conflicts of interest

There are no conflicts to declare.