The Beginning of HCN Polymerization: Iminoacetonitrile Formation and Its Implications in Astrochemical Environments

Hydrogen cyanide (HCN) is known to react with complex organic materials and is a key reagent in the formation of various prebiotic building blocks, including amino acids and nucleobases. Here, we explore the possible first step in several such processes, the dimerization of HCN into iminoacetonitrile. Our study combines steered ab initio molecular dynamics and quantum chemistry to evaluate the kinetics and thermodynamics of base-catalyzed dimerization of HCN in the liquid state. Simulations predict a formation mechanism of iminoacetonitrile that is consistent with experimentally observed time scales for HCN polymerization, suggesting that HCN dimerization may be the rate-determining step in the assembly of more complex reaction products. The predicted kinetics permits for iminoacetonitrile formation in a host of astrochemical environments, including on the early Earth, on periodically heated subsurfaces of comets, and following heating events on colder bodies, such as Saturn’s moon Titan.


■ INTRODUCTION
In this work, we use steered ab initio molecular dynamics to unveil the reaction mechanism for base-catalyzed formation of iminoacetonitrile ( Figure 1, compound 2), a suspected key prebiotic reaction intermediate, in liquid hydrogen cyanide (HCN). 1,2 HCN is one of the most ubiquitous small molecules in the universe, having been observed in the interstellar medium, 3 in the coma of several comets, 4 in the atmosphere of the giant planets, 5,6 on Pluto, 7 as well as on Saturn's moon Titan. 8 Because of its high energy content and reactivity, HCN is a potential driving force for chemistry in various astrochemical environments. 9−11 For example, both the Voyager Flyby and the Cassini mission to Titan have detected clouds made of HCN-based aerosols that can be expected to contribute to concentrated deposits on the considerably colder surface. 10 Observations and modeling suggest that HCN can react to form various complex organic materials on Titan despite the low temperature. 12,13 While there are large uncertainties regarding the nature of early Earth's atmosphere, HCN is expected to have also been formed there as a product of reactions between nitrogen and methane following solar UV radiation, superflares, shockwaves, and discharge chemistry. 14 −16 A long-standing hypothesis first put forth by Oro, 1 Volker, 20 Ferris, and colleagues 2 posits that iminoacetonitrile, the most stable HCN dimer, 21 is the first intermediate in the basecatalyzed polymerization of HCN in the liquid phase ( Figure  2). The potential importance of this molecule has since been pointed out by many (e.g., refs 19, 22) and it has been proposed as an early intermediate in the HCN-based synthesis of various biologically relevant molecules, such as purines, pyrimidines, pterins, and amino acids. 1,23−25 HCN-derived polymers are also of interest as functional materials, e.g., in catalysis or as adhesives and coatings for biomedical applications. 13,22 Given the presumed key role of iminoacetonitrile in prebiotic chemistry, it is important to answer the questions: where can this molecule form and how?
The formation routes to iminoacetonitrile can be expected to depend greatly on the chemical and physical environment. Iminoacetonitrile has recently been detected in the interstellar molecular cloud Sagittarius B2, using the Green Bank Telescope, and in the molecular cloud G + 0.693, using the Institute for Radio Astronomy in the Millimeter Range (IRAM) 30 m telescope. 26,27 Several neutral, 28−30 radical, anionic, 31 and cationic 28 reaction mechanisms have been explored computationally, in efforts to describe the formation of iminoacetonitrile from HCN in the gas phase. However, such mechanisms appear to correspond to either prohibitively high activation barriers 28−31 or unfavorable thermodynamics 28 to allow for meaningful (thermal) reaction rates in the interstellar medium. Computational studies of HCN dimerization on ice grains have predicted similar unfavorable kinetics. 32 A reaction network study suggests that the cyanide radical and methylene imine might be important precursors to iminoacetonitrile in molecular clouds. 33 In this work, we focus on explaining the formation of iminoacetonitrile in neat liquid HCN. HCN polymerization reactions in polar solutions are known to be base-catalyzed and to proceed in a wide range of temperatures (195−373 K) 34,35 at concentrations ranging from 0.01 M to pure HCN (26.2 M). 17,18 However, while it has been hypothesized that iminoacetonitrile is required for the formation of many larger molecules, 1,19,36 the compound has never been directly observed during HCN polymerization experiments. Iminoacetonitrile has been generated by pyrolysis of cyanoformamide salts and has been well characterized in the gas phase and in argon matrices by vibrational spectroscopy, mass spectrometry, and UV photoelectron spectroscopy. 37−39 A 1 H and 13 C nuclear magnetic resonance spectroscopy study has also shown that iminoacetonitrile polymerizes rapidly above 233 K. 38 Ferris and colleagues have shown that N-alkyliminoacetonitriles in aqueous cyanide solutions react to form maleonitrile derivatives, which are, in turn, likely polymer intermediates. 36 To the best of our knowledge, there has only been one previous theoretical study of base-catalyzed iminoacetonitrile formation in polar solution by Kikuchi et al. 40 That study relied on an implicit solvent model for water and predicted a reaction barrier (28.7 kcal/mol), which is markedly lower than the corresponding uncatalyzed pathway (103 kcal/mol 40 ), but prohibitively large for reactions at ambient and lower temperatures.

■ RESULTS AND DISCUSSION
Our study of this condensed phase chemistry relies on molecular dynamics simulations. The use of dynamics simulations allows for the explicit consideration and sampling of the solvent environment and is important for three reasons: first, HCN is an exceptionally polar molecule. The dielectric constant of liquid HCN is 144.8 at 278 K, 41 a factor of 1.7 larger than for water at the same temperature (85.8). 42 Moreover, HCN is capable of donating and accepting strong directional hydrogen bonds. These properties of HCN share striking similarities with those of water and are suggestive of an unexplored potential for self-catalysis. Finally, the considered reaction mechanism involves charged species that interact both strongly and directionally with the environment.
Reaction Exploration in Liquid HCN. We here rely on density functional theory (DFT) metadynamics simulations to identify viable formation routes of iminoacetonitrile from HCN while minimizing bias with regard to mechanism or solvent role. Our metadynamics approach relies on pathcollective variables to define reaction coordinates. These variables are unbiased in the sense that they are only based on aspects of the coordination pattern of reactants and products. In our exploration, the reaction path is defined by an interpolation between the reactant and the product states. We use two path-collective variables: the s coordinate, which describes a structure's position along the reaction path, and the z coordinate, which describes how distant a given structure is from the path. 43 The simulation temperature of 278 K was chosen to be in the middle of the liquid range of HCN, which is close to that reported in recent studies on HCN polymerization. 18 The reaction pathway observed in our simulations ( Figure 3, panel a) resembles a concerted version of the mechanism initially proposed by Oro. 1 The reaction proceeds through the formation of a carbon−carbon bond (d CC ) between HCN and a cyanide anion that occurs concurrently with a proton transfer (d NH ) from a second HCN, which, in turn, regenerates the catalyzing base. Representative snapshots from our simulations displaying geometries of reactant I, transition state TS, and product II are shown in Figure 3a. The TS has been identified through committor analysis, as the structure from which multiple simulation trajectories reach reactant and product basins with equal probability (see Methods section).
Solvent participation remains distinct along the entire reaction pathway, and at minimum two hydrogen-bonded HCN molecules actively partake. This is in line with previous simulations of liquid HCN, which have predicted HCN, on average, to have seven molecules as closest neighbors. 44 The s coordinate allows us to assign the transition state as early (see the Methods section for details), in agreement with the Hammond postulate. The lowest-energy reaction profile for the observed formation of iminoacetonitrile was retrieved with umbrella sampling. Structures I, TS, and II are indicated in path-collective variable space in the resulting free-energy landscape shown in Figure 3b.
Kinetics and Thermodynamics of Iminoacetonitrile Formation. Our simulations predict a relative Gibbs energy of −7.1 ± 0.8 kcal/mol for the overall reaction, suggesting that the mechanism is thermodynamically allowed (see Supporting Information Figure 3). The uncertainty range here refers to the Figure 1. Selection of compounds proposed to form in HCN reaction mixtures. This work explores how HCN 1 can react to form iminoacetonitrile 2. Iminoacetonitrile is a suspected key intermediate in the formation of more complex molecules, such as the HCN tetramer diaminomaleonitrile 3, 17 biologically relevant molecules such as adenine 4 1 , and various polymers (e.g., 5, 6). 18,19 Iminoacetonitrile can also be referred to as C-cyanomethanimine.  Figure 3a) predict the formation of the E-form of iminoacetonitrile. The other possible conformer, Z-iminoacetonitrile, has previously been calculated to be lower in energy in vacuum. 38,45,46 However, in this study, the context is different: we consider iminoacetonitrile in a strongly interacting environment, and our best estimate puts the E-conformation at ∼1 kcal/mol below the Zconformation at 278 K (Supporting Information Figure 5). The Gibbs energy difference between the two structures is small and within the accuracy of the methods we use. Consequently, the specific conformation of iminoacetonitrile does not affect our conclusion on the reaction mechanism nor its kinetics.
The Gibbs energy activation barrier, which corresponds to the TS structure in Figure 3, computes as 15.5 ± 1.2 kcal/mol from our simulations. However, there is a caveat: because of limitations in accuracy of practically feasible electronic structure methods, chemical rate estimation from ab initio molecular dynamics is a well-known challenge (see, e.g., ref 47). The DFT method we use is a necessary compromise of quality and feasibility and one that is known to underestimate barrier heights. 48 To improve our estimate of the reaction profile, we have implemented a correction to the simulated data using calculations at higher levels of theory (see the Methods section). Our final best estimate for the activation barrier and reaction energy lies in a range of 21.8 ± 1.2 and −1.4 ± 0.8 kcal/mol, respectively. Thus, we predict the reaction to be thermodynamically allowed, but only marginally so.
Our barrier estimation is markedly lower than that by Kikuchi et al. (28.7 kcal/mol), who used an implicit solvation model of water to study the same nucleophilic attack but in the presence of ammonium and hydronium cations. 40 Of course, the preference for forming ion pairs in highly polar media will depend on the ion concentration. However, given that our modeling results predict a lower barrier and because the concentration of ammonia (or other bases) used in typical HCN polymerization experiments tends to be low (0.03−0.5 M), 18,49 we can conclude that (a) cation coordination is not necessary for iminoacetonitrile formation and that (b) active involvement of solvent molecules (and explicit consideration of these interactions in modeling) is necessary to explain the base catalysis of this reaction.
Implications of Our Results for the Polymerization of HCN. Polymerization of HCN typically results in complex reaction mixtures with many high-molecular-weight products that have poor solubility in common solvents. 50 Product compositions of such experiments also vary significantly with reaction conditions such as initial concentrations, 51 the presence of oxygen, and the degree of conversion. 52 These often ill-understood products may play important roles in astro-and prebiotic chemistry. As a consequence, definite identification and characterization of products of HCN polymerization are both important and exceedingly challenging. Here, we aim to better understand the potential role of iminoacetonitrile in the formation of HCN polymers. We do so by comparing our theoretical predictions of iminoacetonitrile formation kinetics with observed reaction rates in HCN polymerization experiments. Of course, the rate of HCN dimerization is not necessarily the same as that of its polymerization or processes leading to other reaction products. Implicit in our following comparison with the experiment are two assumptions: first, that the initial step in HCN's basecatalyzed polymerization proceeds through the formation of iminoacetonitrile, and, second, that this initial step is ratedetermining. Aside from its known rapid polymerization into a "red-brown-black material" above 233 K, 38 the precise reactivity and scope of iminoacetonitrile chemistry are largely unknown and will be the subject of future work.
In Figure 4, we have assumed pseudo-first-order kinetics for iminoacetonitrile formation and extrapolated the correspond- ACS Earth and Space Chemistry http://pubs.acs.org/journal/aesccq Article ing reaction time scales to temperatures that are both higher and lower than our simulation temperature of 278 K. We stress that the extrapolation of our results to temperatures other than 278 K comes with various approximations, the first one being negligible changes to the height of the Gibbs reaction barrier with temperature. Our extrapolation also implies that mass transfer limitations remain similar to those in the liquid state and that the reaction mechanism does not change. It is, consequently, expected that our estimates may disagree with observed time scales of reactions in low-temperature (<200 K) and high-temperature (>300 K) conditions. Shown together with our data in Figure 4 are experimental time scales for HCN polymerization reported for several different conditions. Close to a simulation temperature of 278 K, a barrier of 22 kcal/mol should proceed on the order of days (Figure 4). In other words, our rate estimate of condensed phase HCN dimerization is in reasonable agreement with the time scales of pure HCN polymerization observed by Mamajanov and Herzfeld at 278 K 18 and He et al. 19 at room temperature. The temperature dependence of the tetramerization rate measured at 283−313 K by Sanchez et al. 17 suggests an apparent activation barrier of 20.5 kcal/mol, also in close agreement with our barrier estimate. Our predicted reaction time scale also agrees well with aqueous HCN polymerization experiments performed by Miller and colleagues 34 as well as Sanchez et al. 2 near 250 K. We here remind that aqueous HCN solutions that are close to or below the eutectic point (250 K) form a eutectic phase in which HCN is concentrated (to 74.5 mole percent). 55 We can speculate that effects associated with such a phase may be one reason for the sometimes strikingly different experimental results near the eutectic temperature. For example, whereas Sanchez et al. have reported the formation of the HCN tetramer diaminomaleonitrile after days at 251 K, 2 Marin-Yaseli et al. were unable to observe meaningful reactivity after a month of reaction time at 253 K. 49 As expected, our approximate extrapolation in Figure 4 does not agree as well with experimental measurements performed at low temperatures. In particular, the one-of-a-kind experiment by Miller and colleagues, in which polymerization of HCN at 195 K was observed over ∼25 years, 34 implies a lower reaction barrier in the largely frozen solution compared to that predicted for the liquid (albeit due to a low ammonia−water eutectic temperature of ∼173 K, the aqueous NH 4 CN solution may not have been entirely frozen 34 ). The latter comparison implies that crystal surfaces and defects in frozen HCN solutions may be better polymerization catalysts compared to the liquid system studied herein. Relatedly, reaction kinetics studies performed between 323 and 363 K by Mas et al. suggest that different reaction mechanisms, and barrier heights, are at play at higher temperatures. 53 Possible Role of Iminoacetonitrile in Different Astrochemical Environments. The time scale required for the studied dimerization drops rapidly with temperature (c.f., Figure 4), with important consequences for the possible environments in which the reaction mechanism may be relevant. Our simulations correspond well to setups of some laboratory experiments, but one might question the direct relevance of these conditions in nature. Where might there exist, or have existed, concentrated solutions of HCN? The atmospheric production rate of HCN on the early Earth has been predicted to be 10 −8 −10 −6 cm 2 /s and much of what was produced is believed to have been deposited on Earth's surface (30 million tons/year). 14 However, despite a substantial production in the atmosphere, the concentration of HCN in Earth's early oceans is likely to have been much too low to allow for polymerization. 56 Eutectic freezing could have occurred in shallow pools or under very harsh glacier conditions and corresponds to a plausible route to concentrated aqueous HCN. 17,23 Our calculations predict that the 250 K eutectic temperature of a HCN−water mixture should allow for the formation of iminoacetonitrile, and possibly subsequent HCN polymerization, over a time scale of months.
Observations of cometary coma 4 have established comets as another environment rich in HCN, and the possibility of HCN polymerization in comets has been pointed out by many (see, e.g., refs 57, 58). The surface temperature of comets varies considerably with their proximity to a star (from ∼190 K at 4 AU to ∼340 K at 1 AU). 54 Heating-and-cooling cycles of comets may allow the formation of eutectic solutions of HCN with a high enough concentration to allow for iminoacetonitrile formation on relevant time scales. We stress that the formation of such solutions requires pressures of hundreds of Pa, an order of magnitude higher than that expected on the surface of comets. Whether subsurface liquids can form during heating of the cometary surface, or during cometary explosive eruptions, is not known and is the subject of ongoing research. 54 HCN is one of the main products of the atmospheric photochemistry on Saturn's moon Titan. 8 The temperature of Titan's atmosphere varies drastically with altitude but ranges from 200 K at 300 km 59 to ∼94 K at the surface. 60 The formed HCN can be expected to deposit on the surface and in the hydrocarbon seas of Titan in large amounts. 10,61 The abundance of HCN dissolved in the hydrocarbon seas is uncertain but predicted to be low 62,63 and orders of magnitudes smaller than the lowest reported molar concentration used in aqueous HCN polymerization experiments  19 and e Mas et al. 53 Colored vertical lines indicate the temperature in our simulations (orange), representative measures of cometary surfaces in the inner solar system (blue and red), 54 and the eutectic freezing temperature of HCN and water mixtures (green). 55 ACS Earth and Space Chemistry http://pubs.acs.org/journal/aesccq Article (0.02 mole %). 17 The cyanide anion, the catalyst in our study, is believed to be the most common anion in the chemical haze surrounding the large moon. 64 However, provided that we assume a similar reaction barrier in the hydrocarbon seas of Titan as in neat HCN, the low surface temperature (<94 K 65 ) would hinder the reaction from occurring within relevant time scales. Whereas solid-state equivalences to the reaction considered herein might prove operable over very long time scales, our current simulations describe HCN in the liquid state. In other words, for the iminoacetonitrile formation reaction considered in this work to be feasible on Titan, temporary heating of HCN surface deposits or melting of atmospheric aerosols would be required. Such heating might be generated by impact events or, possibly, coincide with runaway exothermic polymerization of acetylene, another main product of Titan's atmospheric chemistry. 66−68 ■ CONCLUSIONS HCN in the liquid state is complex, featuring an extremely large dielectric constant, strong directional intermolecular interactions, and acid−base equilibria that allows for its participation in catalytic cycles. Our simulation of basecatalyzed dimerization of HCN explains the formation of iminoacetonitrile in the liquid state. The time scales estimated for iminoacetonitrile formation are similar to those observed for HCN polymerization experiments, which suggests that the dimer formation may be the rate-limiting step for a host of subsequent reaction chemistry, including polymerization. The predicted rate of formation of iminoacetonitrile allows for its formation in several environments and supports a potential central role of the molecule in astro-and prebiotic chemistry. Barring detection by nuclear magnetic resonance and vibrational spectroscopy or successful trapping experiments, we suggest that the laser desorption ionization of frozen reaction mixtures coupled to time-of-flight mass spectroscopy may provide one route to experientially proving the role of iminoacetonitrile in HCN polymerization.  Table 1) and were omitted in the interest of reducing computational costs. The GTH pseudopotential 73 was used with a cutoff of 280 Rydberg. Liquid HCN was simulated in a cubic box with a side length of 15.94 Å subjected to periodic boundary conditions. The number of molecules were tailored to 64, including one cyanide anion, so as to accurately reproduce the density of HCN at 278 K (0.709 g/cm 3 ). 74 The system was propagated in an NVT ensemble at 278 K with a 0.5 fs time step. The steered simulations and the subsequent committor analysis simulations were performed with the canonical sampling through velocity rescaling thermostat 75 and a thermostat time constant of 50 fs.
Metadynamics. The s and z path-collective variables were constructed following the example of Branduardi et al. 43 (see Supporting Information Section 4). The path-collective variables used during the metadynamics 76 simulations were defined based on the bond topology of the carbon in the cyanide anion, through the distance function by Pietrucci and Saitta. 77 The PLUMED library version 2.5.0 78 was used to steer the dynamics simulations. Two reference structures were used: a cyanide dissolved in HCN and a solvated iminoacetonitrile, leaving the path between those two states open for exploration. Following a 5 ps equilibration, the reference structures were simulated for 20 ps prior to an analysis of their coordination patterns. The metadynamics simulation was initiated from the reactant state.
Committor Analysis and Umbrella Sampling. The transition state (TS) of the reaction observed in the metadynamics simulation was identified through a committor analysis. 79 Ten out of 20 trajectories starting from the TS structure ended up in the reactant basin and in the product basin, respectively. A free-energy profile of the reaction was then constructed using umbrella sampling 80 along the s coordinate using configurations from trajectories obtained in the committor analysis as starting points. The umbrella sampling simulations were run for 11 ps after 1.5 ps of equilibration. The force constants, k, used for the harmonic potentials in these windows (Supporting Information Table 2) were chosen so that k where σ 2 is the variance of the s coordinate in a typical simulation. The umbrella sampling windows were spaced at most 0.025 s coordinate units apart (corresponding to 3% of the distance between the s coordinates of reactants and products). Together with the chosen force constants, this separation provided a good overlap of the sampling in different windows (Supporting Information Figure 2). The weighted histogram analysis method (WHAM) 81 version 2.0.9 was used to retrieve the free-energy profile of the reaction. Uncertainties in the energy estimates were estimated using a block averaging method in WHAM. 82 Energy Correction to the Free-Energy Profile. As a necessary compromise between accuracy and computational cost, all simulations were done using the PBE-D3 functional. The largest contribution to the error in the free-energy profile can be attributed to inherent errors in the potential used to calculate the electronic (Born−Oppenheimer) energy ΔE. This conclusion was reached by comparing the thermal corrections ΔG − ΔE when calculated with PBE-D3 and the more accurate B3LYP 83,84 -D3 method (Supporting Information Section 1). Hybrid exchange−correlation functionals, such as B3LYP, represent the most accurate level of theory that can be practically implemented on our system. The B3LYP functional has reported average errors of 2 and 4 kcal/mol for association reactions and hydrogen transfer, respectively. 48 To reduce overall errors in the free energy, an electronic correction term, ΔE corr , was added to the free-energy profile ΔG PBE , according to ΔG final = ΔG PBE + ΔE corr , where ΔG final is our best estimate of the relative Gibbs energy (see Figure S4). The correction term, ΔE corr , was obtained by averaging 105 energy evaluations of configurations from the umbrella sampling windows corresponding to reactant (s value 1.09), transition state (s value 1.32), and product (s value 1.88) at two levels of theory. The configurations were taken at 0.1−0.2 ps intervals from the trajectories. ΔE corr was calculated as ΔE corr = ΔE B3LYP-D3,avg − ΔE PBE-D3,avg , where the latter two terms are average energy differences calculated with B3LYP-D3 and PBE-D3, respectively. The ΔE corr calculations were