A computational model to predict the Diels–Alder reactivity of aryl/alkyl-substituted tetrazines

Abstract The tetrazine ligation is one of the fastest bioorthogonal ligations and plays a pivotal role in time-critical in vitro and in vivo applications. However, prediction of the reactivity of tetrazines in inverse electron demand Diels–Alder-initiated ligation reactions is not straight-forward. Commonly used tools such as frontier molecular orbital theory only give qualitative and often even wrong results. Applying density functional theory, we have been able to develop a simple computational method for the prediction of the reactivity of aryl/alkyl-substituted tetrazines in inverse electron demand Diels–Alder reactions. Graphical Abstract Electronic supplementary material The online version of this article (10.1007/s00706-017-2110-x) contains supplementary material, which is available to authorized users.

Due to the high reaction rates, these ligations can be used in time-critical applications such as rapid radiolabeling and pretargeted PET imaging [10][11][12][13][14][15][16] and provide high yields within short reaction times even at low concentrations as usually encountered in radiochemistry and in vivo. Therefore, kinetics is one of the most important characteristics of bioorthogonal reactions. However, prediction of reactivities using the chemist's understanding of organic chemistry, especially of IEDDA reactions, might lead to wrong predictions [17] and only qualitative estimates. In addition, synthesis of tetrazines is often low yielding and involves handling or even requires the production of anhydrous hydrazine (not commercially available in Europe), which limits the feasibility of screening for high Diels-Alder reactivity. Hence, there is the need of reliable computational tools to predict the reactivity of various tetrazines in TLs.
Herein, we introduce a computational model for the prediction of the reactivity of aryl/alkyl-substituted Tz in the cycloaddition with trans-cyclooctene (TCO), thus eliminating the need for expensive and even dangerous synthetic work, finally enabling in silico screening for tetrazines with desired reactivity.

Results and discussion
Recently, we have investigated the reactivity of several 3-aryl-6-(3-fluoropropyl)-1,2,4,5-tetrazines 1-8 as chemical probes for rapid radiolabeling and pretargeted PET imaging (Fig. 2). While the alkyl substituent is the same for all eight tetrazines the aryl component shows considerable variation including electron-rich and electron-poor aryl groups. In addition, Tz 9 and 10 were included to investigate the influence of the alkyl group and an orthosubstituted aryl residue, respectively.
These experimental results were selected as a basis for the construction of a predictive computational tool. DFT was successfully used in the past by our group [16,22] and others [7,17,23,24] to predict or explain the reactivity of dienophiles and tetrazines in the tetrazine ligation. Therefore, the Minnesota density functional M06-2X in combination with the 6-311?G(d,p) basis set, was used as model chemistry. This density functional has been proven to produce accurate results for thermodynamics of cycloaddition reactions [25,26].
Diels-Alder reactions can be described by HOMO/ LUMO interactions using the frontier molecular orbital (FMO) theory. In case of the IEDDA cycloaddition the main orbital interaction is between a low-lying unoccupied orbital of the dienophile, usually being the LUMO?1 for aryl/alkyl tetrazines (Fig. 4a) [17], and the HOMO of the electron-rich dienophile (in this case TCO, Fig. 4b). According to the FMO theory, a smaller energy gap between the interacting orbitals facilitates the reaction. Thus, one might expect that a more electron withdrawing and thus LUMO?1-lowering substituent accelerates the reaction, while an electron-rich aryl substituent will decrease reactivity. HF/6-311?G(d,p)//M06-2X/6-311?G(d,p)-calculated orbital energies are shown in Fig. 4c. Tz 4 and 7 bearing an electron-withdrawing trifluoromethyl or sulfone group, respectively, show the lowest LUMO and LUMO?1 energies. However, the tetrazine with the highest reactivity, Tz 8, has one of the highest LUMO and a rather high LUMO?1 energy within the series. As shown in Fig. 5, there is no significant correlation between the LUMO?1-energy levels and the rate constants (R 2 = 0.07). This can be rationalized by the fact that FMO interactions are not the only major contributors to activation energies in such cycloadditions [17,[27][28][29]. Therefore, the FMO theory cannot be applied for the reliable prediction of the Diels-Alder reactivity of aryl/alkyl-substituted Tz.
However, the calculated energies of activation (DE à ) of the reactions between tetrazines 1-10 and TCO (11) show excellent linear correlation with the natural logarithm of the measured rate constants, as per transition state theory (Fig. 6  Correlation between natural logarithm of second-order rate constants and M06-2X/6-311?G(d,p) calculated energies of activation activation can, therefore, be used to reliably predict the reactivity of aryl/alkyl-Tz in IEDDA reactions with transcyclooctenes.
Using the linear correlation between ln(k) and calculated DE à , Eq. (1) can be constructed which allows the prediction of the rate constant of new aryl/alkyl-Tz for the reaction with TCO (11) in anhydrous 1,4-dioxane at 25°C.

Conclusion
Prediction of the reactivity of Tz in bioorthogonal ligation reactions is essential to reduce synthetic work, the resulting costs and associated risks. However, the chemist's prediction of these reactivities, mainly based on frontier molecular orbital theory, can only yield qualitative or even wrong results as shown in this work. We were able to develop a computational method for the prediction of reactivities of aryl/alkyl-substituted 1,2,4,5-tetrazines in the reaction with trans-cyclooctene based on M06-2Xcalculated energies of activation. This method is computationally cheap as it requires only the optimization of the tetrazine and the corresponding transition state with transcyclooctene at the M06-2X/6-311?G(d,p) level of theory, which can even be done on an average desktop PC within hours.
We are convinced that our method will aid the development of new aryl/alkyl tetrazines for bioorthogonal applications, and represents a step forward to the development of a universal computational tool being able to predict the reactivity of tetrazines with various substitution patterns.

Experimental Computational methods
Calculations were performed at the M06-2X/6-311?G(d,p) level of theory in the gas phase using the Gaussian 09 Rev. D.01 software package [30]. Vibrational analysis was performed to confirm stationary points are energetic minima or transition states, respectively. Orbital energies were calculated by performing a HF/6-311?G(d,p) single-point calculation on M06-2X/6-311?G(d,p) optimized structures. DE à was determined by calculating the difference in energies between transition state and reactants. Data analysis was preformed using Gaussview 5 and Chemcraft. XYZ coordinate files of all reactants and transition states are available as electronic supplementary material. Calculated energies for Tz 1-10 and the respective transition states (for the reaction with TCO) are shown in Table 1.

Kinetic measurements
Kinetic measurements were performed on a SX20 stoppedflow spectrophotometer (Applied Photophysics, UK) using a 535 nm LED light source. A 4 mM solution of 11 and approximately 0.1 mM solutions of Tz 1-10 were prepared in anhydrous 1,4-dioxane (note: for correct measurements it is of utmost importance to use anhydrous 1,4-dioxane,  since even small amounts of water will accelerate the reaction, leading to irreproducible data). These solutions were loaded into the driver syringes and equal volumes of TCO and Tz solution were mixed, resulting in concentrations of 2 and 0.05 mM, respectively. The reaction progress was followed by measuring the absorption around 535 nm. All measurements were performed in triplicates. Observed reaction constants (k obs ) were determined by non-linear regression (one-phase-decay) using Prism 6 (Graphpad) and second-order rate constants (k) were calculated from k obs . Table 2 lists k obs , c TCO , c Tz and k for all reactions.