Role of Acid in the Co-oligomerization of Formaldehyde and Pyrrole.

Building on previous work (J. Phys. Chem. A 2017, 121, 8154-8166) under neutral conditions, we examined the co-oligomerization of CH2O and pyrrole to form porphryinogen under acidic conditions using density functional theory (B3LYP//6-311G**). Thermodynamically, we found that azafulvene intermediates were significantly stabilized under highly acidic conditions. Kinetically, energy barriers were lowered for C-C bond formation, discriminating in favor of reactions that lead to porphyrinogen. However, it was challenging to satisfactorily combine our thermodynamic and kinetic profiles into a unified free-energy profile because of difficulties in optimizing transition states of cationic species involving proton hops. Instead, we used neutral carboxylic acids as a proxy to study how energy barriers changed. By combining data from both neutral and acidic conditions, we estimate a free-energy profile for the initial steps of oligomerization under milder acidic conditions more relevant to prebiotic chemistry.

Introduction

In a previous computational study, we constructed a free-energy map for the co-oligomerization of formaldehyde and pyrrole to form porphyrinogen,[1] as shown in Figure . We calculated the thermodynamics and kinetics of various elementary reaction steps in this system under baseline conditions at neutral pH in aqueous solution. In that work, we alluded to a follow-up study of how the thermodynamics and kinetics might change under acidic conditions. We hereby present the results of our follow-up study—things turned out to be more complicated and generating an overall energy profile was not as straightforward.

Figure 1

Overall reaction for porphyrinogen formation.

However, before diving into those details, here is why we are studying this system: porphyrinogen is an intermediate in a proposed prebiotic synthesis of porphyrin,[2] and our lab is interested in calculating free-energy profiles of reactions related to the origin of life—before extant life evolved to synthesize porphyrins enzymatically. Although porphyrins from unsubstituted pyrrole are less relevant to prebiotic chemistry, we do not consider these as we are comparing the present work to our previous study that only considered unsubstituted pyrrole and did not include oxidation reactions.

A brief summary of the relevant experimental work is as follows. Starting from the simple reactants paraformaldehyde and pyrrole, Hodgson and Baker reported one of the earliest syntheses of porphyrin;[3] their reaction used neat solutions of the reactants along with metal salts in sealed tubes heated to 84 °C. Fox and Strasdeit prepared pyrroles by heating amino acids, then added millimolar concentrations of formaldehyde, and were able to form both porphyrinogens and porphyrins.[4] Nitrite/nitrate oxidants added to the reaction mixture were likely responsible for the subsequent oxidation of porphryinogen to porphyrin.[5,6] Notably, these solutions were quite acidic at pH 2. Subsequently, the Lindsey group synthesized porphyrins at micromolar concentrations of formaldehyde and close-to-neutral conditions by adding micelles and vesicles to the reaction mixture.[7] For a comprehensive history of porphyrin synthesis, going back to the initial work by Rothemund,[8,9] we refer you to The Porphyrins.[10] Relevant to our present study, these experiments showed that acidic conditions played a role in increasing yields of porphyrin synthesis.

In the present study, we generated a thermodynamic profile for the co-oligomerization of formaldehyde and pyrrole to form porphryinogen under highly acidic conditions by adding a proton to all reaction species. Note that this represents extreme acidic conditions unlikely to be encountered in prebiotic chemistry, and we explain our strategy and its limitations in the Results and Discussion section. This approach allowed us to see that the major difference in the thermodynamic profile, compared to neutral conditions, is the stabilization of azafulvene intermediates in cationic species. A smaller effect observed was the thermodynamic stabilization of reactions involving C–C bond formation and destabilization of reactions involving C–N bond formation. These effects favor the formation of porphyrinogen over other side-products.

However, we encountered numerous challenges optimizing transition states to generate an analogous kinetic profile for cationic species, and after little success, we pivoted to an approach utilizing “neutral” organic acids as catalysts on neutral reactants (without the added proton). Using formic acid as a catalyst, we found that barriers were significantly lowered for C–C bond formation en route to porphyrinogen. Besides determining which reaction steps are acid-catalyzed, our study also provides a qualitative feel for the relationship between acid strength and reaction barriers; qualitative because we cannot easily combine the thermodynamic and kinetic contributions on a single map because the former represents extreme acidic conditions whereas the latter utilizes a neutral molecular acid to directly donate and accept protons. Furthermore, directly comparing the relative free energies in aqueous solution of cationic versus neutral species is challenging because of large magnitude differences in solvation free energies.

Although this is a less-satisfying result, we think it valuable to discuss the challenges in this study, and hope that we (or others) will find some way to better connect the two in the future. In the present study, we provide one possible sketch of the initial steps of formaldehyde and pyrrole reacting to form dipyrrlmethane under milder acidic conditions more relevant to prebiotic chemistry, which emphasizes the potential role of transient (cationic) azafulvene intermediate species, while acknowledging the limits of our approach.

Computational Methods

The protocol discussed below is similar to our previous work, and therefore much of the text in this section is reproduced from earlier work for clarity and reading ease.[1,11−13] As we directly compare present calculations under acidic conditions with previous work under neutral conditions,[1] it is essential to maintain the same protocol. All calculations were carried out using Jaguar 10.0[14] at the B3LYP[15−18] flavor of density functional theory (DFT) with a 6-311G** basis set. For stable molecular species (with no negative eigenvalues), conformers were generated with MacroModel 11.5 using the Merck Molecular Force Field.[19,20] Only structures within 5 kcal/mol of the lowest energy conformer were included in a restricted torsional sampling search. If more than 64 conformers were generated, we selected 64 low-energy conformers with a variety of starting structures to maximize the probability of finding global minima. Transition states were built individually both in extended forms and more compact structures with internal hydrogen bonds. Although we made good faith efforts to use a wide range of starting structures, it is possible to fail in locating the lowest energy transition state.

The Poisson–Boltzmann continuum approximation[21,22] was used to describe the effect of water as a solvent applying the Jaguar defaults: a dielectric constant of 80.4 and a probe radius of 1.40 Å. The forces on the quantum mechanical solute atoms because of the solvent can be calculated in the presence of the solvent. However, as in previous work, the solvation energy was calculated at the optimized gas-phase geometry because in most cases there is practically no change between the gas-phase and implicit solvent-optimized geometries. It is important to note that even though the solvation energy contribution is to some extent a free-energy correction, it certainly does not account for all of the free energy. A comparison of our chosen level of theory, basis set, and implicit solvent scheme, with other methods can be found in previous work.[23] Although water is a reactant in hydration reactions, concentration corrections are not included in this landscape, the advantages and disadvantages of which are also discussed in previous work.[24]

The analytical Hessian was calculated for each optimized structure, and the gas-phase energy corrected for zero-point vibrations. Negative eigenvalues in transition-state calculations were not included in the zero-point energy (ZPE). The temperature-dependent enthalpy correction term is straightforward to calculate from statistical mechanics where we assume that translational and rotational corrections are a constant times kT, that low-frequency vibrational modes will generally cancel out when calculating enthalpy differences, and that the vibrational frequencies do not change appreciably in solution. The vibrational scaling factor of 0.967 for B3LYP//6-311G** was not applied because when relative energies are calculated, the difference to the enthalpy correction becomes negligible within the computational error.

The corresponding free-energy corrections in solution are much less reliable.[25−27] Changes in free-energy terms for translation and rotation are poorly defined in solution, particularly as molecular size increases. Additional corrections to the free energy for concentration differentials among species (to obtain the chemical potential) can be significant, especially if the solubility varies among the different species in solution. Furthermore, as the reactions being studied are in solution, the free energy being accounted for comes from two different sources: thermal corrections and implicit solvent. Neither of these parameters is easily separable, nor do they constitute all the required parts of the free energy given our approximations. To estimate the free energy, we followed the approach of Lau and Deubel[28] who assigned the solvation entropy of each species as half of its gas-phase entropy. Wertz[29] and Abraham[30] had previously proposed that upon dissolving in water, molecules lose a constant fraction (∼0.5) of their entropy. Recent computational studies in other unrelated systems have come to the same conclusion.[31−33]

Although we calculated multiple conformers, only the most stable conformer (both stable minima and transition states) for each unique molecular species is reported in our free-energy map. The relative free energies in our free-energy profiles are calculated with respect to the reference molecules: formaldehyde, (neutral or protonated) pyrrole, and water. Assigning these molecules relative free energies of zero allows us to quickly and easily visualize a map of the energy landscape for the myriad reactions that can take place. Our free-energy values, designated G298, include the ZPE, enthalpic, and entropic corrections to 298 K for reactions in aqueous solution. The parsing of energies into their various contributions can be found in the Supporting Information along with the relative free energies shown in our thermodynamic profiles.

Our protocol was put to the test with a detailed comparison of our computational results with NMR measurements for the self-oligomerization of a 1 M solution of glycolaldehyde.[24] Our calculated equilibrium concentrations of the dominant species in solution (monomers and dimers) agreed very well with experiment. In addition, our calculations allowed us to successfully predict the concentrations of trimers in solution, 2 orders of magnitude lower than the monomers. Our protocol performs well in calculating the relative free energies of stable species, typically within 0.5 kcal/mol of experimental results, or an uncertainty of within a factor of 2.3 in terms of equilibrium constant ratios.[24,34]

To optimize transition states, additional water molecules or organic acids were explicitly added to the system to find the lowest energy barrier for proton transfers. For cationic species, our attempts to find unconstrained transition states were largely unsuccessful, as will be described in the Results and Discussion section. This led us to a different approach: utilizing “neutral” organic acid molecules as potential catalysts. These transition states were much easier to optimize, and had one large negative eigenvalue corresponding to the reaction coordinate involving bond breaking/forming and accompanying proton transfer. Whereas these calculated barriers can be compared to our previous work in the absence of acid,[1] they assume the thermodynamic profile of neutral zero-charge species, and cannot be mapped directly on the thermodynamic profile generated for protonated species.

Results and Discussion

Thermodynamic Profile: Neutral Conditions

The top half of Figure , a subset of the free-energy map from our previous work,[1] shows the main pathway leading to porphyrinogen formation under neutral conditions. The relative free energies (ΔGr) are with respect to H2O, CH2O, and pyrrole (all three assigned a relative free energy of zero). The five main reaction types are color-coded according to the figure key. Below are two examples of how the values in Figure are obtained. Detailed energy contributions for each species can be found in Table S1 of the Supporting Information.

Figure 2

Thermodynamic profiles of porphryinogen formation.

Example 1: For the addition of CH2O to the α-position of pyrrole (1) to form 3, the free-energy change is G298(3) – [G298(CH2O) + G298(1)] = −9.0 kcal/mol; this reaction is exergonic. As CH2O and pyrrole are both reference states, the relative free energy of 3, ΔGr(3), is −9.0 kcal/mol. Example 2: the dehydration of 3 to form the azafulvene 5 is endergonic; the change in free energy is [G298(5) + G298(H2O)] – G298(3) = +11.0 kcal/mol; but as water is a reference state and ΔGr(3) = −9.0 kcal/mol, ΔGr(5) = +11.0 – 9.0 = +2.0 kcal/mol, as shown in Figure .

A brief summary of the key results are as follows; for full details, see our previous work.[1] When CH2O is added to pyrrole, addition to the nitrogen (forming 2) is less thermodynamically favorable compared to addition to the carbon. The α-position (forming 3) is marginally favored over the β-position (forming 4). Adding a second pyrrole requires accessing the azafulvene 5. The dehydration of 3 to form 5 is an endergonic reaction (Example 2 above). The subsequent reaction of 1 and 5 to form 6 (red arrow) is exergonic by 21.8 kcal/mol. Alternatively, 3 and 5 can form 6 (purple arrow) while eliminating CH2O; this reaction is less exergonic (by 12.8 kcal/mol). Successive additions of CH2O and pyrrole eventually lead to porphyrinogen (15) formation.

Thermodynamic Profile: Highly Acidic Conditions

To mimic acidic conditions, we added a proton to all species—a strategy we previously used to study reactions of aldehydes with imidazole.[35] However, unlike protonated imidazole with a pKa of 6.95, protonated pyrrole has a pKa of −3.8 (based on 1H NMR spectra in 16 M sulfuric acid).[36] Thus, the cationic species in our calculations mimic very high or extreme acidic conditions, unlike the mild acidic conditions more relevant to prebiotic chemistry. To compare the relative thermodynamic stability of multiple species, we chose to maintain the same overall charge across the board because large differences in solvation free energies between species with different net charge can lead to spurious comparisons. Although the extreme acidic conditions are artificial, our approach provides a bird’s eye view of overall profile differences between cationic species versus neutral species.

The thermodynamic profile of the cationic species is shown in the bottom half of Figure . The method used to calculate relative free energies is similar to the neutral case, except that protonated pyrrole (1a) is included in the list of reference molecules under acidic conditions. The most stable form of protonated pyrrole is when the added proton is in the α-position (1a). The β-protonated pyrrole (1b) is 4.4 kcal/mol higher in free energy, whereas addition to the nitrogen (1c) is 15.7 kcal/mol higher, as shown in Figure . Detailed energy contributions for each species are provided in Table S2 of the Supporting Information.

Figure 3

Relative free energies of selected protonated species.

Experimentally, the pKa of 1a is −3.8, the pKa of 1b is −5.9, and the pKa of 1c has not been directly measured but estimated to be approximately −14.[36,37] Using ΔG = −RT ln Ka, the acid dissociation reaction 1a + H2O → H3O+ + 1 yields ΔG = −5.2 kcal/mol for Ka = −3.8. Similarly, the reaction 1b + H2O → H3O+ + 1 yields ΔG = −8.0 kcal/mol for Ka = −5.9; and the reaction 1c + H2O → H3O+ + 1 yields ΔG = −19.1 kcal/mol for Ka = −14. Based on these numbers, 1a is the most stable of the protonated pyrroles; 1b would be 2.8 kcal/mol higher in free energy, whereas 1c would be 13.9 kcal/mol higher. These are reasonably close to our calculated estimates of 4.4 and 15.7 kcal/mol respectively.

Returning to Figure , addition of CH2O to the pyrrole nitrogen (forming 2a) is less thermodynamically favored under highly acidic conditions (ΔG = −2.5 kcal/mol) compared to neutral conditions (ΔG = −4.7 kcal/mol). In contrast, addition to carbon is more thermodynamically favorable at both the α (protonated: ΔG = −12.4 kcal/mol; neutral: ΔG = −9.0 kcal/mol) and β (protonated: ΔG = −10.5 kcal/mol; neutral: ΔG = −8.5 kcal/mol) positions. The gap between 3a and 4a has widened to ∼2 kcal/mol under highly acidic conditions, compared to a gap of 0.5 kcal/mol between 3 and 4 under neutral conditions. As 3a leads to porphyrinogen, but 2a and 4a do not, we surmise that this is one thermodynamic contribution favoring product formation under highly acidic conditions.

In the hydroxymethylpyrroles, the proton in the α-position is still favored, as shown in Figure . Comparing 2a and 2b, adding a proton to the α-carbon is more stable than adding to the nitrogen by 15.6 kcal/mol, very similar to the difference between 1a and 1c. If hydroxymethyl occupies one of the α-positions, the proton prefers the unsubstituted position (3a over 3b) with a free-energy difference of 6.8 kcal/mol. (A similar situation is found comparing 6a and 6b; the free-energy difference is 7.0 kcal/mol.) If hydroxymethyl occupies the β-position, the proton prefers the α-position adjacent to the CH2OH moiety (4a over 4b). Our calculations are in agreement with experimental 1H NMR studies showing that “protonation of a 2-alkylpyrrole occurs at the 5-position, whereas 3-alkyl groups direct protonation into the neighboring 2-position, as observed in other electrophilic substitution reactions”.[37]

The most significant difference in the thermodynamic profile under highly acidic conditions compared to neutral conditions is the stabilization of protonated azafulvene relative to α-hydroxymethylpyrrole. In Figure , under neutral conditions, the dehydration of 3 to 5 is endergonic by 11.0 kcal/mol. This is not surprising because 1 (and by extension 3) is considered aromatic, with considerable resonance energy stabilization (75% of benzene),[37] but the neutral azafulvene 5 is not. In the protonated case, the analogous dehydration reaction (3a to 5a) is endergonic by only 2.7 kcal/mol. Unlike their neutral counterparts, 1a and 3a are not aromatic, and therefore there is no equivalent loss of aromaticity. One might even expect some stabilization as two additional π-electrons are added to the delocalized system, and in fact, dehydration is favorable for the dimer (7a to 8a), trimer (10a to 11a), and tetramer (13a to 14a). As dehydration is the only endergonic step on the path to porphyrinogen under neutral conditions, our calculations suggest that the effect of adding acid is to significantly favor this step, and therefore the overall reaction profile, thermodynamically.

The intermediates formed from subsequent additions of pyrrole and CH2O have slightly lower relative free energies under highly acidic conditions (compare 6a to 6 and 7a to 7, e.g.) and the overall profile is similar to neutral conditions. Exceptions are the azafulvene-like species as mentioned above, and the final product porphyrinogen. 15a has a marginally higher relative free energy than 15 because the added proton has no choice but to be located on an α-carbon thus reducing resonance stabilization.

To summarize, a bird’s eye view of the thermodynamic profiles under neutral versus highly acidic conditions show overall broad similarities, as seen by comparing the relative free energies in Figure . There are two notable differences: (1) the significant stabilization of azafulvene-like species when a proton is added to the system and (2) increased discrimination favoring CH2O addition in the α-position. Both these effects favor porphyrinogen formation. The caveat: the extreme acidic conditions modeled by adding a proton to all species is unrealistic. The thermodynamic profile under more relevant milder acidic conditions might be extrapolated to lie in between the two profiles we have generated in Figure . We will return to this matter after discussing the issues we encountered exploring the kinetics of this system.

Kinetics under Acidic Conditions

Our attempts to find optimal transition states of cationic species (overall +1 charge) under highly acidic conditions were largely unsuccessful. We were only able to cleanly optimize one unconstrained transition state structure (Figure S1)—the direct addition of 1 and 5a to form 6a (red arrow in Figure ). The forming C–C bond is 2.32 Å and the activation barrier is 8.1 kcal/mol. For comparison, the analogous reaction under neutral conditions is 14.7 kcal/mol from our previous work.[1] We optimized two unconstrained transition states (Figure S2) for the dehydration reaction forming 5a, but neither transition state had the appropriate eigenvectors along the reaction coordinate. In one case, the proton hops from one α-carbon to the other; the C–H bond is partially breaking, but the C–O bond remains fully intact. In the other case, the proton has completely hopped off the pyrrole although the breaking C–O bond distance is at a reasonable 1.89 Å. Energies for all three transition states are in Table S2.

We were unable to find reasonable unconstrained transition states for any of the CH2O addition reactions forming 2a, 3a, or 4a. In most cases, the proton hops off (through a chain of water molecules transfers) to the incoming CH2O, forming the CH2OH moiety which is essentially the product structure. In some cases, the CH2OH moiety and its “water cluster” separate from the pyrrole. We attempted several constrained optimizations moving stepwise along (what seemed to us) desirable reaction coordinates, but the barriers were unrealistically high and no smooth “curve” joining reactant to product along the reaction coordinate could be found. Releasing the constraints led to reactant or product complexes. We were similarly unsuccessful in our attempts to find the transition state for dimerization with CH2O elimination, that is, the reaction of 3a and 5a to form 6a.

Given these difficulties, we pivoted to a different strategy—using an organic acid as a proxy acid catalyst of the reaction types studied. We are not the first to employ this approach; for example, it was used to compare acetaldehyde hydration in the gas phase, in neutral aqueous solution, and in the presence of a carboxylic acid.[38] With this new strategy, optimizing unconstrained transition states was now straightforward, but introduced other complications—we discuss these below as we present the results.

Instead of calculating reaction barriers for every reaction in Figure , we chose to study representative reactions: addition of CH2O to pyrrole, dehydration to azafulvene, and addition of the second pyrrole. These reactions are shown in Figure ; detailed energy contributions of each transition state are found in Table S3 of the Supporting Information. We initially chose formic acid as the proxy acid catalyst because it is considered an origin-of-life building block,[39] and is a catalyst for studying prebiotic chemistry in hydrothermal vents[40] and the RNA world.[41] It is also a major product in early Miller–Urey spark discharge experiments.[42] In Figure , the energy barriers for formic acid are in green font. The suffix fa indicates that the catalyst is formic acid, and the prefix TSX designates the transition state forming product X. Parenthetical values in blue font are the energy barriers under neutral (aqueous) conditions from our previous work.[1]

Figure 4

Representative acid-catalyzed reactions. Energy barriers are in kcal/mol.

We considered addition of CH2O to three pyrrole positions (nitrogen, α-carbon, β-carbon) in the presence of formic acid. In Figure , comparing the energy barriers where formic acid (green) facilitates proton transfer rather than solely relying on water molecules (blue), we find the following. There is practically no difference for CH2O addition to the nitrogen (TS2fa). However, there is significant stabilization for CH2O addition to the α (TS3fa) or β (TS4fa) position by ∼10 kcal/mol. Notably, addition to the α-carbon has a lower barrier (+13.1 kcal/mol) compared to the nitrogen (+16.9 kcal/mol) in the presence of formic acid. When only water was present, the opposite was true: adding to the nitrogen was kinetically more favorable (+16.7 kcal/mol instead of +22.8 kcal/mol to the α-carbon).

The dehydration of 3 to 5 has a marginally higher barrier with formic acid (+15.1 kcal/mol) compared to water (+13.8 kcal/mol). In contrast, the C–C bond-forming dimerization reactions both have lower barriers with formic acid. The direct addition of 1 and 5 to form 6 has a barrier of +9.0 kcal/mol in formic acid compared to +14.7 kcal/mol in water. The addition of 3 and 5 to form 6, while simultaneously eliminating CH2O, has a barrier of +7.3 kcal/mol in formic acid compared to +10.9 kcal/mol in water. Interestingly, the barrier of 8.1 kcal/mol in the one case (addition of 1 and 5a) where we successfully optimized a protonated transition state is close to the 9.0 kcal/mol barrier with formic acid, but this may be coincidental.

From these results, we conclude that the role of acid in the kinetic profile is to lower the barriers for C–C bond formation when CH2O and pyrrole moieties are added along the co-oligomerization pathway to form porphyrinogen. On the other hand, acid does not aid the dehydration reaction or the formation of new C–N bonds.

Can we relate this to the thermodynamic profile? In the previous section, we saw that adding a proton to the system slightly stabilized C–C bond formation products by 2–3 kcal/mol, whereas the main effect was to significantly stabilize the azafulvene intermediates formed via dehydration (by 8–10 kcal/mol). However, adding a proton in the thermodynamic profile mimics extreme acid conditions, whereas it is less clear what conditions are mimicked in the kinetic profile by using molecular formic acid as a proxy for proton transfer in the transition state. To further examine how our calculations might simulate varying acidic conditions, we completed similar transition-state optimizations with acetic acid (aa), fluoroacetic acid (faa), and trifluoroacetic acid (tfaa). The pKa values for these three acids are 4.75, 2.59, and 0.23 respectively; whereas the pKa of formic acid is 3.75.

Unfortunately, the picture becomes more complicated. In Figure , if we consider only the energy barriers with acetic acid and its fluoro-substituents, we see that barriers decrease in all three cases with increasing acidity. Compared to aa, the faa and tfaa barriers are lower by 3–5 and 7–9 kcal/mol respectively. This seems to qualitatively fit well with the decrease in pKa values. However, the formic acid (fa) barriers are very similar to faa, and quite different from aa, even though the pKa of fa is almost midway between faa and aa. Is formic acid an outlier?

The transition states for several of the fa-catalyzed reactions are shown in Figure ; analogous structures for aa shown in the Supporting Information look very similar. Nothing seems amiss, and the bond-forming and bond-breaking transition-state distances are typical. For CH2O addition, the forming N–C or C–C bonds are ∼1.8 Å. For dehydration, the breaking C–O bond is ∼2.1 Å. For the dimerization reactions, the forming C–C bond is ∼2.1 Å with fa (and 2.4 Å with aa) in the direct addition of 1 and 5, similar to the direction addition of 1 and 5a discussed earlier. As expected, the acidic proton is further away from the carboxylate “catalyst” and closer to the formaldehyde–pyrrole oligomers. The analogous optimized transition-state structures for the other acids look very similar; and formic acid does not look like an outlier structurally.

Figure 5

Transition-state structures with formic acid catalyst.

Another puzzle is that the barriers for acetic acid in two cases are noticeably higher than for the water-only case. These two cases, CH2O addition to the pyrrole nitrogen (TS2aa) and dehydration (TS5aa), are the ones where our initial calculations with formic acid suggested no lowering of the energy barriers. It seems odd that acetic acid is “worse” at proton transfer than employing a chain of only water molecules. If we add an additional water molecule to the TS2aa transition state, it can be placed either in between acetic acid and pyrrole, or between acetic acid and formaldehyde (see Figure S3 in the Supporting Information). The former has a barrier of +21.6 kcal/mol, that is, adding the water makes no difference; but the latter has a barrier of +15.6 kcal/mol, much closer to the water-only barrier of +16.7 kcal/mol (although marginally lower and within the computational error). However, using the same strategy with the fluoro-substituted acids showed no decrease in the barrier. After attempting to optimize multiple transition states with additional interposed water molecules for several reactions, and finding no consistency, we decided to stick with having only the carboxylic acid as the proton shuttle—this at least allows us to systematically cross-compare the energy barriers using proxy acids, and these are the numbers shown in Figure .

A reviewer suggested trying newer DFT functionals to ascertain if any of the problems we were encountering was an artifact of B3LYP, but we ran into the same issues, and the results were less satisfactory; the results and accompanying discussion can be found in the Supporting Information.

Perhaps a lesson here is that one should be cautious when using a molecular acid as a proxy catalyst to mimic the environmental conditions with a target pH, and that using a series of acids to mimic different pH conditions may result in unexpected computational artifacts. Although a clear trend of lower barrier with lower pKa is observed comparing acetic acid, fluoroacetic acid, and trifluoroacetic acid, it is less clear if formic acid can be mapped on the same scale as the other three. There are many approximations leading to our free-energy values. The entropic approximation that worked well in our previous work when only water molecules were added to a transition state may introduce significant errors, as might the implicit solvent approximation, and in future work we plan to study this in a simpler system with clear experimental data over a range of pH values.

In the meantime, although the results do not fit into a tidy package, we can at least summarize our main findings, focusing on formic acid as the prebiotically plausible proxy for acid catalysis. We find that the most significant kinetic effect is lowering the barrier for CH2O addition to pyrrole to form a new C–C bond. Although addition to both the α and β positions is more facile, α remains kinetically favored over β by ∼4 kcal/mol. Pyrrole addition to the azafulvene intermediate, also a C–C bond-forming reaction, is also kinetically favored in the presence of formic acid. On the other hand, for CH2O addition to the pyrrole nitrogen (forming a C–N bond) and for the dehydration reaction, the barriers are similar for both formic acid and water. Thus, the key C–C bond-forming reactions en route to porphyrinogen formation are catalyzed by formic acid. The picture is less clear when comparing different carboxylic acids.

It is also unclear if the kinetic profile can be directly combined with the thermodynamic profile given the different approach used to mimic (different) acidic conditions. One possibility is to combine the barriers for the strongest molecular acid in our study, trifluoroacetic acid, with the “highly acidic” thermodynamic profile. However, given that 1a (pKa −3.8) is still much more acidic than tfaa, it is unclear if this is a reasonable match-up. Furthermore, it is weak acid conditions that are relevant to prebiotic chemistry and the motivation for this study. Thus, we will explore one possible way of extracting data from our seemingly disparate thermodynamic and kinetic studies, to shed some light on the role of acid under plausibly prebiotic conditions; to this we now turn.

Combining Thermodynamics and Kinetics: A Free-Energy Profile

One way to combine the data in both parts of Figure on the same scale is to peg the relative free energy of 1a to be +5.2 kcal/mol above 1, given the −3.8 pKa value of 1a as previously discussed. This is shown in Figure for the initial reaction steps leading to the formation of dipyrrlmethane (6). The relative free energies of the relevant neutral molecules (1, 3, 5, 6) are represented with thick black bars, whereas the thick red bars refer to cationic species: 1a is pegged at +5.2 kcal/mol, and the free energies of all other cations that might possibly be involved (from Figure ) can be referenced to 1a. Structures of all these species are also included in Figure for clarity.

Figure 6

Example of a combined energy profile with weak-acid catalysis.

Thin horizontal lines indicate transition states: blue for water, green for formic acid. (We chose formic acid because of its prebiotically plausibility as previously discussed.) The relative free energies (Table S3, Supporting Information) are already pegged to 1 as the transition states are overall neutral. Note that the barriers in Figure are for specific reactions, whereas the values in Figure are all relative to 1. Thus, for the first step these are identical (+22.8 kcal/mol and +13.1 kcal/mol for water and formic acid, respectively). However, for the dehydration step, as 3 has a relative free energy of −9.0 kcal/mol and the calculated barriers are +13.8 and +15.1 kcal/mol for water and formic acid, respectively (see Figure ), the relative free energies of TS5wat and TS5fa are +4.8 and +6.1 kcal/mol, respectively, in Figure .

Blue diagonal lines connect the water-only transition states to the neutral molecules. This represents the uncatalyzed pathway: 1 → TS3wat → 3 → TS5wat → 5 → TS6wat → 6. We expect that azafulvene 5 will not be easily observed as the reverse barrier of 5 hydrating back to 3 is tiny (+2.8 kcal/mol). This pathway has already been extensively discussed in our previous work under neutral conditions.[1]

Green diagonal lines connect the formic acid transition states with both the neutral molecules and potential protonated intermediates, most of which are likely to be transient and unobservable. The acid-catalyzed barrier (+13.1 kcal/mol) for CH2O addition is significantly lowered compared to water-only. The transition state (TS3fa in Figure ) is structurally close to a complex of (the Wheland intermediate) 3b and formate. As previously discussed, 3a is more stable than 3b, an easy proton-hop away; and in either case subsequent removal of the proton leads to 3; black dashed lines indicate the connection 3b → 3a → 3. It is highly unlikely that 3a or 3b can be isolated experimentally, although 3b may be a very transient intermediate. By repositioning the formate, we can optimize the transition state TSpt3fa (pt for proton transfer, structure shown in Supporting Information) to pull the proton off the α-carbon. TSpt3fa is slightly lower in energy than TS3fa; however, given the instability of 3b relative to 3, the proton can simply be removed by solvent.

In the dehydration step, TS5fa is calculated to be slightly higher in free energy than TS5wat. In Figure , the transition state of TS5fa resembles a complex of the structure 3d (in Figure ) and formate. However, in our calculations, adding a proton to the oxygen in 3 to form 3d automatically results in water moving away. We were able to optimize an unconstrained state with a C···O distance of 2.6 Å (at relative free energy −0.7 kcal/mol) indicated by the dashed C–O line in 3d, but there is practically no barrier to push the water molecule further away to form separate 5a and water moieties; hence, 3d is unlikely to be experimentally observed as an intermediate.

Azafulvene 5a is the only case where the corresponding protonated species is lower in relative free energy than its neutral counterpart (5). However, 5a is still thermodynamically less stable than 3 and the low hydration barrier suggests that, if present, 5a is transient and would be difficult to observe experimentally. As TS5wat is marginally lower in energy than TS5fa, it could operate as the transition state connecting 3 to 5a (blue dashed line) even in acidic conditions.

Adding pyrrole to azafulvene, TS6fa is 5.7 kcal/mol lower in free energy than TS6wat, suggesting that this step is acid-catalyzed. Similar to the CH2O addition step, TS6fa resembles a transition state involving 6b and formate (see Figure ). A proton hop might lead to the lower energy 6a (with an energy gap similar to 3b and 3a), but abstracting the proton (TS6tfa) is also a low-barrier reaction, and either pathway leads to the neutral dipyrrlmethane 6. We expect 6b and 6a to be difficult to observe even if they represent intermediates along the reaction pathway.

Throughout this section, our language has been tentative. Our scheme in Figure may not be the best way to combine the thermodynamics and kinetics; using molecular formic acid as the proxy catalyst representative of weakly acidic conditions is a significant assumption; and pegging the relative free energies of protonated species to their neutral counterparts essentially by a 5.2 kcal/mol “correction” factor might be erroneous. In weak acid, the free energies of local minima are likely to lie in between the neutral and protonated states, although closer to the neutral. Our scheme does highlight that thermodynamic stabilization of azafulvene (although transient as an intermediate) may be important in the acid-catalyzed reaction. We are embarking on future work to test this scheme in a more tractable system with good experimental data.

Conclusions

The present work explores the co-oligomerization of CH2O and pyrrole to form porphyrinogen under acidic conditions compared to neutral conditions in aqueous solution at 25 °C. We can summarize our results as follows. Azafulvene intermediates, less stable under neutral conditions, are significantly stabilized thermodynamically under acidic conditions. Kinetically, the presence of (formic) acid lowers barriers for C–C bond formation as successive CH2O and pyrrole monomers are added to the growing oligopyrrole. Acidic conditions also favor addition of CH2O to the α-carbon of pyrrole both thermodynamically and kinetically, thus favoring the formation of intermediates en route to porphyrinogen, compared to side reactions leading to other products.

It was challenging to connect our thermodynamic profile with our kinetic (barrier) estimates because of the different methodologies employed. Adding a proton mimics unrealistic extreme acidic conditions, but provides consistency when comparing the free energies of species within the thermodynamic profile of cationic species. The difficulty in locating transition states led us to use neutral carboxylic acids as proxy catalysts, although this brought up other complications. We generated one example of a unified free-energy map, to provide some sense of the effect of a weak acid on this system, but it rests on less secure assumptions. Future work will explore other approaches to the problem, and we hope to generate a scheme allowing finer discrimination of thermodynamic and kinetic effects as a function of pH in aqueous solution.