Probability of reaction pathways of amine with epoxides in the reagent ratio of 1:1 and 1:2

The algorithm for generating and estimating the probability of possible reaction pathways for multichannel bimolecular interactions was used to predict the reaction products in the reagent ratio of 1:1 and 1:2. Here, we have considered the possible reaction pathways of the reaction of amine (1S,2S,4S)-bicyclo[2.2.1]hept-5-en-2-ylmethanamine (1) with epoxides 2-((cyclohexyloxy)methyl)oxirane (2), 2-(phenoxymethyl)oxirane (3), N-(oxiran-2-ylmethyl)-N-phenylbenzenesulfonamide (8) in order to explain experimental observed data, which indicate differences in the reactivity of glycidyl ethers and glycidylsulfonamide with framework amines. Based on the proposed algorithm (Borysenko et al. 2021), we have investigated the reaction in the reagent ratio of 1:1 and 1:2. Calculated values of activation barriers indicate a low probability of formation of interaction products of amine (1) with epoxide (8) with a (1:2) reagent ratio due to steric hindrances in the reaction center.


Introduction
Mono-and diepoxide compounds have gained wide popularity due to the broad range of their applications as building blocks for the construction of biologically active systems [2][3][4][5][6], monomers for highly biocompatible and biodegradable polymers [7][8][9][10][11][12], bases for adhesives [13], and modifiers of polymer compositions [14]. Epoxy cycle opening reactions can be considered as model in the study of complex processes of crosslinking of epoxy monomers with hardeners or their polymerization. In this regard, the study of the mechanisms of opening of the epoxy cycle is undoubtedly an important and urgent task of modern chemistry.
For such compounds, the study of mechanism and the dependence of the ratio of products on the ratio of the initial reagents was already published [16]. The aim of this work was to study the possible reaction pathways of the aminolysis of epoxides (2,3,8) considering the conformational lability of the reagents in order to explain the experimentally observed formation of products (6,7) and the absence of product (10).
Over the last decades, several computer-assisted automatic methods for systematic searching reaction pathway and determining reaction mechanisms have been proposed. Among these methods, the following should be noted: anharmonic downward distortion following (ADDF) method [17][18][19][20]; artificial force-induced reaction (AFIR) method [21][22][23][24][25]; transition state search using chemical dynamics simulations (TSSCDS) method [26][27][28][29]; the ZStruct2 method developed by Zimmerman's research group [30][31][32]; and molecular dynamics and coordinate driving (MD/CD) method proposed by Yang et al. [33]. Since each method for searching reaction pathway aims to provide a balance between high accuracy, comprehensive potential energy surface search and low computational costs, the development of new, efficient algorithms for studying reaction paths remain an urgent issue for computational chemistry.

Computational details
The program PCModel v 9.2 [34] was used to perform conformational analysis at the MMX force field level [35]. Conformational search calculations were carried out using the GMMX technique.
Semi-empirical calculations by PM7 [36] method were carried out using MOPAC2016 [37] program, DFT calculations at the M062X/6-31G(d) [38] level of theory were performed using the Gaussian 16 program [39]. Barrier heights evaluated at semi-empirical level were calculated using PM7-TS method as single point calculations based on PM7 geometry. Vibrational frequencies were calculated for all stationary points to confirm whether the optimized geometry corresponds to a minimum or a TS. The intrinsic reaction coordinate [40,41] calculations were performed to approve the reaction path and the relationship of the given transition structure to the local minima of the reagent and product. IRC calculations were performed for the most favorable reaction pathways in the gas phase and in the solvent. The solvent effect was taken into account in two ways: explicitly two solvent molecules in vacuo and in the bulk solvent within the continuum solvation model SMD [42].

Results and discussion
In order to study the features of aminolysis of epoxides (2, 3, 8), we applied the algorithm for generation and assessment of probability of possible reaction pathways for multiple channel bimolecular interactions proposed in our previous study [1], which was devoted to the interaction of epoxide (2) with amine (1) to obtain glycidyl ether (4).
This algorithm includes conformational search for reaction intermediate using Molecular Mechanics (MMX) approach, based on the obtained conformation construction of structures of transition states and pre-reaction complexes, and calculation of activation energies to further determine the probable reaction pathways.
In this work, initially, we studied the reaction of amine (1) with epoxides (3,8) in the ratio of reactants 1:1. Then, the investigation of the reaction in the ratio of reagents 1:2 was carried out according to the analogous procedure.
The reaction we have considered is exergonic and proceeds through S N 2-like mechanism forming bipolar ion (5i, 9i) as an intermediate on the rate-limiting endorgonic stage (Scheme 2).
Located conformations of bipolar ions (5i, 9i) could be straightforwardly applied for construction of starting geometry for optimization of transition states of the corresponding reactions. This statement follows from Hammond's postulate which states that in the case of endorgonic stage, the geometry of the intermediate should be close to the geometry of the transition state.
The proposed strategy consists of the following steps: (IV) Locating TSs, pre-reaction complexes at M062X/6-31G(d) level of theory. As starting geometries, we have used TS geometries that had been obtained in step (III). Similarly, to step (III), the repeating structures were excluded from the sample. (V) The reaction paths obtained by M062X/6-31G(d), whose contributions to the overall rate constant were the largest, have been selected for calculation of reaction paths considering the influence of the solvent.
Let us now apply the strategy described in the steps I-V in more detail for reaction of amine (1) with epoxides (3,8).
The assumed ratio of reactants is 1:1.

Conformational analysis
Calculations were carried out using MMX force field within GMMX technique as implemented in PCModel v 9.2 program.
Since the conformational transitions in the intermediates (5i, 9i) are associated with rotation around the bonds, in order to generate the initial structures by the GMMX method, the algorithm that randomly selects a subset of the bonds intended for rotation was chosen. The N-C-C-O torsion angle value was fixed at 180° because this angle value corresponds to the trans-opening state of the epoxy ring. All other torsion angles were used to create conformations by the GMMX method. Based on the previous study [1] for

Locating of TSs and pre-reaction complexes at semi-empirical (PM7) level
Geometry of conformations from previous step has been modified by setting up length of forming N-C and breaking O-C bond lengths equal to 1.782 Å and 1.999 Å, respectively. This corresponds to the transition state geometry parameters of model reaction (interaction of methylamine with oxirane) (see [1]).
After localization of TS structures and exclusion of repeating ones the pre-reaction complexes and energy barriers were calculated. Table 1 shows obtained activation barrier values for the first step of aminolysis reaction and contributions of each routes to the total reaction rate constant. The overall contributions to the total reaction rate constant (k i ) were calculated using Eq. (1).
According to PM7 method, the largest contributions from 14.9 to 15.7% to the total reaction rate are made by TSs that corresponds to the reaction path numbers 3, 4, 5 in the case of intermediate (5i). For intermediate (9i), the largest contribution to the total reaction rate equal to 67.3% is made by TS that corresponds to the reaction path 5.
For all these transition states, the presence of a hydrogen bond (NH•••O) is observed, which leads to the stabilization of the structure (Fig. 1).

Study of reaction in vacuo at M062X/6-31G(d) level of theory
The starting geometries for the localization of transition states were structures obtained by the PM7 method.
All of the 18 possible reaction paths obtained for intermediates (5i, 9i) by the PM7 method were studied at the M062X/6-31G(d) level of theory. After excluding TSs with the same geometries, 7 TS conformations were obtained for intermediate (5i) and 9 for intermediate (9i).
Starting geometry for optimization of pre-reaction complexes was generated by displacement of atoms in the TS structures along imaginary normal vibrational mode. The activation barriers were calculated for each reaction channel (Table 1). The largest contribution of 73.8% to the total reaction rate is made by TS of intermediate (5i) that corresponds to path 8. In accordance with the results for intermediate (9i), the largest contributions to the reaction rate within 21.1-26.2 are from the TS conformations that correspond to paths 2, 4, 7, 9. The structures of these TS conformations are shown in Fig. 2.
The most energetically favorable reaction paths were also examined considering solvent effects.

Study of reaction with a (1:2) reagent ratio in vacuo
The procedure for studying the aminolysis reaction with a (1:2) reagent ratio was similar to that described above for a (1:1) reagent ratio. Let us discuss the obtained results.
At the second step of the study, the transition states and pre-reaction complexes for each found conformer of all three intermediates were localized at PM7 and M062X/6-31G(d) levels of theory. Subsequently, the activation barriers of possible reaction channels and the contribution of each TS to the total reaction rate were calculated. As can be seen from Table 2, the results of the semi-empirical and the DFT methods are slightly different. The largest contribution to the overall reaction rate for intermediate 6i is from the TS number 2 (67.8%) obtained by PM7 method. Interestingly, at the M062X/6-31G(d) level, in addition to TS corresponding to path 2, which has the largest input of 18.2%, several TSs   6i, 7i, 10i) for possible reaction channels and the contribution of located pathways to overall rate constant of the reactions of glycidyl ethers (4, 5) with epoxide (2, 3) and glycidylsulfonamide (9) with epoxide (8) Fig. 3.

Study of reaction with explicit consideration of solvent at M062X/6-31G(d) level of theory
Since the experimental reaction was carried out in the presence of the protic solvent 2-propanol [16], the specific solvation with solvent was modeled in this study by adding molecules of alcohol to the investigated molecular system. The purpose of considering solvent molecules into the system was to simulate the "relay" proton transfer, which can be carried out at the expense of several molecules of 2-propanol. We performed a simulation of the effect of the explicit solvent in vacuo and in solution of 2-propanol with a dielectric constant of ε = 19.264. The initial position of the solvent molecules in the transition state was chosen so that simultaneous activation of epoxide and amine occurred through the formation of hydrogen bonds between the molecules alcohol and epoxide, alcohol and alcohol, alcohol and amine. The number of solvent molecules sufficient for "relay" proton transfer and the formation of the reaction product for all cases except epoxide (3) was equal to two molecules. To simulate the reaction of epoxide (3) with amine (1) in a solvent, taking into account the peculiarities of the conformations of the transition states, namely the distance of the amino group from the epoxide oxygen, we selected three solvent molecules.
Only those reaction pathways of epoxide aminolysis in which all of these hydrogen bonds were present were investigated. The TS structures of aminoalcohols (5,6,7,9,10) in the presence of the solvent are given in Fig. 4, and their Cartesian coordinates are included in the Supporting Information.
Simultaneous activation of epoxide and amine by solvent molecules led to a decrease activation barrier of the reaction. As can be seen from Table 3, all reactions with a (1:1) reagent ratio are characterized by barriers from 136.9 to 148.4 kJ/mol in vacuo, and from 71.7 to 85.6 kJ/mol within explicit consideration of solvent molecules.
For intermediate 10i with a (1:2) reagent ratio, the TS conformation of the most probable reaction path number 7, found in vacuo ( Table 2, Fig. 3f), is characterized by the formation of N-H•••O hydrogen bond, which, on the one hand, stabilizes this transition state, and on the other hand, is the cause of steric inaccessibility for interaction with solvent molecules. In other words, for the formation of the aminoalcohol product, conformational changes must take place first, so that the "relay" proton transfer from the amino group to the epoxide oxygen becomes possible. In contrast, TS of the less energy efficient path number 3 has an open amino group to interact with solvent molecules. In this case, the simulation of the reaction with the participation of the solvent gave the values of free Gibbs energy of activation larger, by 25.1-27.2 kJ/mol, than for the reaction with a (1:1) reagent ratio. The obtained results are consistent with the experimental data, according to which, in this case, only the product of the interaction with a (1:1) reagent ratio is formed.
The plausible reason for the observed increase in the energy barrier is the large steric hindrance for the interaction of bulk reagents and solvent molecules in a process characterized by a (1:2) reagent ratio. The obtained activation barriers are consistent with experimental data, which indicate differences in the reactivity of glycidyl ethers and glycidylsulfonamide with framework amines.

Conclusions
Conformational analysis of primary products of reactions could be considered as a general procedure for generation of transition state structures in the case of multiple channel reactions. Method of generation of transition state structures for multiple channel reactions that considers the high conformational mobility of the reactants and the possibility of varying their orientation during the reaction was tested in this study. Conformational analysis of investigated epoxy esters derivatives, conducted using method of molecular mechanics MMX, revealed structures with the lowest energies -the most stable conformers of primary products of epoxy compounds aminolysis. Transition state structures for the reaction of aminolysis of glycidyl ethers and glycidylsulfonamides with frame amines were localized using quantum-chemical methods PM7 and M062X/6-31G(d). Formation of 1:1 and 1:2 products of glycidyl ethers aminolysis is determined by steric effects of conformational flexible reagents and bulk solvent molecules and formation of N-H•••O hydrogen bond in transition state. Calculated values of activation barriers indicate a low probability of formation of interaction products of amine with epoxyglycidylsulfonamides with a (1:2) reagent ratio than with a (1:1) reagent ratio, which is consistent with and explains experimental data.