The thermodynamic and kinetic aspects of midazolam ring closure from benzophenone to benzodiazepine form, its acid–base equilibria and aromaticity: a quantum-chemical study

The process of midazolam ring closure was studied from the thermodynamic and kinetic points of view by means of quantum-chemical methods. B3LYP/6–311 +  + G(d,p) model was employed for gas phase and water environment (polarizable continuum model) calculations. It was concluded that the reaction rate determining step is the first step—carbinolamine formation from amine and carbonyl ends of the opened benzodiazepine ring. The Gibbs free energy of activation was calculated as 35.1 kcal/mol for gas- phase and 33.9 kcal/mol for water environment. Intrinsic reaction coordinate calculations were performed to verify that the transition state really connects the substrate and product. Thermodynamically, this reaction is endoergic with ΔG = 9.6 kcal/mol for gas phase and 10.0 kcal/mol for water environment. However, the next step—carbinolamine protonation with immediate water molecule loss is expected to be fast and activation barrierless, which enables further progress of the ring closure, despite the positive ΔG of the fist step. Next, the protonated imine undergoes deprotonation to final closed ring, pharmacologically active molecule of midazolam. The whole chain of reaction is exoergic with ΔG equal to − 5.6 kcal/mol for gas phase and − 7.7 kcal/mol for water environment. In order to understand the role of other than benzodiazepine/imidazole molecular fragments on the ring closure process, a model was build which contains only benzodiazepine and imidazole rings. The activation barrier for the carbinolamine formation of the model is similar to midazolam in the gas phase but higher by about 10 kcal/mol for water environment. The most interesting difference is however that for the model, the carbinolamine formation step is exoergic with ΔG equal to − 2.2 kcal/mol for gas phase and − 1.8 kcal/mol for water environment. This difference can be connected to complicated conformational shape of the midazolam molecule, which during the ring closure undergoes unfavorable deformations with accompanying rise of the energy of the molecule. The protonation sites for both midazolam and the model were also studied. In the case of midazolam, the preferred protonation site is the imidazole ring nitrogen atom, but in the case of the model it is the benzodiazepine ring nitrogen atom. The aromaticity of the 5- and 7-membered rings were analyzed using two aromaticity—HOMA and pEDA. It follows that larger stability of the cation protonated at the benzodiazepine ring is accompanied with substantial increase of the 7-ring aromaticity in the model of midazolam. The complexes of midazolam and its model with a water molecule were analyzed because they are needed for evaluation of the energy of the whole process of ring closure. In the case of midazolam, the water molecule preferentially connects to imidazole ring, but in the case of the model, complexes with both imidazole and benzodiazepine rings have similar stability.


Introduction
Benzodiazepines [1] are among the most frequently prescribed medicines, mainly as sedative agents, anxiolytics, and hypnotics. They almost completely surpassed barbiturates, because of wider therapeutic index and much lesser danger of overdose. The core of their structure is the 7-membered heterocyclic ring containing two nitrogen atoms (most commonly in positions 1 and 4), fused with the benzene ring-thus the name "benzodiazepines." Diazepam [2] was the first commonly used drug of this class, still in use to this day. An interesting subgroup of these medicines are ones containing additional azole ring fused with benzodiazepine ring. This introduces more lipophilicity which is connected with better penetration of the blood-brain barrier, they are characterized also by rapid onset and short half-life in the organism. The prototypic drug for this classmidazolam, contains imidazole ring [3], and others contain mostly the 1,2,4-triazole ring: alprazolam [4], estazolam [5], and triazolam [6].
In this paper, we will focus on midazolam (MDZ) as a typical and widely used drug for premedication, induction of anesthesia, and dealing with sleep disorders. In Scheme 1, there is depicted a schematic representation of diazepam and the topic of our study-midazolam.
The important difference of midazolam is that the imidazole ring makes this molecule not only more lipophilic but also more basic and susceptible to benzodiazepine ring opening. In acidic conditions, it undergoes protonation, at first to the cation, and in more acidic conditions to the dication (Scheme 2).
The dication is, however, unstable and undergoes ring opening as suggested first by Gerecke [7]. This process has several steps and transition products and the final product is the open form of midazolam (MDZo), which is called also the benzophenone form of midazolam, because of the characteristic benzophenone fragment (Scheme 3). The same author estimated the half-life of ring closure of MDZ to be about 10 min.
Because of forming of these various ionic forms, MDZ was the first soluble benzodiazepine [8]. This is especially important for intravenous formulations, where precipitation of the medicine is not acceptable [9,10]. The open form of midazolam does not possess pharmacological activity, as it lacks the essential part of the pharmacophore-the benzodiazepine ring. The existence of two different forms of midazolam (MDZ and MDZo) was confirmed by measuring UV spectra of midazolam in various pH conditions [11].
The process of ring opening is fully reversible and the closedring form of midazolam can be formed from the open-ring form, in a reversed chain of reactions. It begins by the amino group nucleophilic attack on carbonyl atom with the formation of carbinolamine. This is, however, very difficult in acidic conditions, when for the most of the molecules, the first-order amine nitrogen atom is protonated (Scheme 3) and lacks the free lone electron pair. Thus, to enable the ring closure process, the pH have to be basic, neutral, or only little acidic, for the first-order amino group must not be protonated. The dependence of solubility and MDZ/MDZo ratio on pH was measured by UV spectrometry [12]. It was established that the solubility of MDZ greatly increases under pH < 4. The pKa = 6 was estimated, which accounts for the fact that MDZ is a weak base. Regarding the ratio of the two forms-MDZ and MDZo, it was established that for pH > 4 the MDZ form dominates, at lower pH the amount of MDZo form gradually increases, and for pH = 2.5 the two forms exist in approximately equal amounts. In more acidic conditions, there is more MDZo form which eventually dominates for pH < 1.5.
In this paper, we will analyze the MDZ/MDZo equlibrium from the perspective of ring closure process, for which we assume the basic/neutral conditions, as the ring closure in vivo takes place in such conditions. Thus, we focus on two neutral (not cationic) forms of midazolam-open (MDZo) and closed (MDZ) which are depicted in Scheme 4.
The process of ring closure has several steps and is in the essence a condensation reaction between first-order amino group and carbonyl group, leading to cyclic imine, which is the final MDZ. The first step of the reaction involves nucleophilic attack of the amino nitrogen atom on the carbonyl carbon atom. In the idealized literature picture from a textbook [13], this is done in two sub-steps: 1. The formation of tetrahedral bipolar transition product, where the oxygen atom is negatively and amino group positively polarized (Scheme 5). 2. The proton transfer from NH 2 + group to O − with the formation of the carbinolamine.
However, our quantum-chemical calculations show that the transition product depicted in Scheme 5 is unstable and cannot be formed. Instead, when moving the NH 2 group closer and closer to the carbonyl carbon atom, at certain distance a proton transfer occurs simultaneously. Thus, the carbinolamine (6A in Scheme 6) is formed in one step instead of steps (1) and (2), thus the proton transfer is rather synchronous, not subsequent to bringing together the NH 2 and CO groups. We obtained the proper transition state for this stage, which will be shown and analyzed later. Following further the literature description of the ring closure, there should exist three postulated transition products depicted in Scheme 6.
First, the carbinolamine 6A is formed, second the protonated carbinolamine is formed-6B to facilitate removal of the water molecule, and third, the protonated imine-6C is formed, which is actually our closed form of midazolam protonated at the benzodiazepine ring. As the form protonated at the imidazole ring is slightly more stable, the proton is expected to be transferred to the imidazole ring, forming the 2A molecule MDZ + (Scheme 2). Eventually, in basic environment, it can be completely deprotonated to neutral MDZ.
The aims of this study are as follows: (1) to model each step of the process of ring closure by means of quantumchemical methods, establishing the transition products and transition states, whenever possible. It will be done for gas phase and water environment conditions; (2) to establish the role of other than benzodiazepine ring fragments of the molecule in the thermodynamics and kinetics of the process of ring closure. Do these fragments ease or make harder this process to occur? For this aim, a simplistic model of MDZ will be build where all fragments of the molecule apart from the benzodiazepine and imidazole rings will be cut out; (3) to study the acid-base equilibrium of midazolam cations and neutral species. To our knowledge, there are no papers dealing with these subjects from theoretical point of view.

Computational details
All calculations were done by using the Gaussian 16 suite of programs [14] with the use of universal B3LYP density functional [15,16]. Preliminary geometry optimizations were performed with the use of Pople's 6-31G(d) orbital basis set [17]. Next, to obtain more reliable structures and energies, the more extended 6-311 + + G(d,p) basis set [18,19] was employed. All optimizations were performed for gas  [20] and followed by frequency calculations to confirm the nature of stationary points-if they are minima, transition states or higher order saddle points on the potential energy hypersurface. To verify that the transition state truly connects the substrate and product, the IRC (intrinsic reaction coordinate) calculations were performed at the B3LYP/6-31G(d) level of theory using transition state optimized at the same level of theory as the starting point. The maximum number of points was set to 50 and the force constant matrix was recalculated every 3 steps. Regarding the energetic parameters, total electronic energy, enthalpy, and Gibbs free energy of the reaction were analyzed. The latest is the most important as it directly corresponds to the reaction equilibrium constant. The starting structures for transition state search were obtained by the relaxed potential energy surface scans with the opt = modredundant option in Gaussian. Aromaticity of the rings of midazolam and midazolam model were analyzed using geometric HOMA [21,22] and electronic pEDA [23,24] aromaticity indices, calculated by using the free AromaTcl [25] software.

Midazolam
In the text, the energetic data for B3LYP/6-311 + + G(d,p) results, both for the gas phase and water environment modeled by PCM model are presented. The structure of midazolam (MDZ) optimized at the B3LYP/6-311 + + G(d,p)/ PCM level of theory is shown in Fig. 1. All further figures are done at this model chemistry.
From comparison of Scheme 1B and Fig. 1, it follows that the molecule is in reality highly nonplanar. The methylene group is mostly responsible for this, it also breaks the pielectron conjugation of the 7-membered ring and reduces significantly its aromaticity. It formally contains 7 pi-electrons (2 from the pyrrole-like nitrogen atom, common to both rings), which according to the Huckel rule makes it formally non-aromatic. The fused benzene ring as well as one connected by a single, rotatable bond are of course aromatic and planar, and so is the imidazole ring. To better understand the aromaticity of midazolam, two aromaticity indices were employed-geometric HOMA and electronic pEDA. However, calculating pEDA makes sense only for planar rings, thus for the benzodiazepine ring only HOMA is available. We will analyze here only 7-and 5-membered rings and will call accordingly HOMA for 7-membered ring-HOMA7, and HOMA for 5-membered ring-HOMA5. The pEDA index will be calculated only for the 5-membered ring and will be called-pEDA5. For midazolam HOMA7 = − 0.30 and HOMA5 = 0.79, also pEDA5 = − 0.07. It follows that the 7-membered ring is non-aromatic (negative HOMA) and 5-membered ring is aromatic with slight pi-deficiency (below the perfect electronic sextet). Later in the text, the changes of HOMA and pEDA also for cations will be analyzed.

Formation of carbinolamine
The first and essential step is the formation of carbinolamine via a transition state. Looking at the open form of midazolam, it is virtually impossible to propose an initial structure of the transition state, sufficiently close to the true transition state, to obtain it by means of direct optimization. The most reliable procedure in such situations is to perform a relaxed potential energy surface scan. A series of geometry optimizations is performed, but in each case a chosen interatomic distance is kept fixed. In this case, we found that the best procedure is to start not from the MDZo, but from carbinolamine. The atoms in consideration are the hydrogen from OH group and the nitrogen atom. We start from their distance in carbinolamine (2.6 Å) and gradually bring the hydrogen closer to the nitrogen atom, up to the distance of about 1.0 Å. The energy at first increases and when the structure is near the transition state, the energy is highest. Then, suddenly the ring opens (the proton now belong to the amino group), the energy lowers and the structure becomes the open-ring form-MDZo. The structure of the highest energy from the scan is a good starting point for full optimization of the transition state. The substrate-MDZo-in proper conformation for upcoming nucleophilic attack, the transition state and the product-carbinolamine are depicted in Fig. 2 below. All molecules are oriented spatially in the way to best show the essential bond lengths during the process of cyclization. It follows from Fig. 2 that the bond lengths in the transition state 2b is much closer to the product 2c that to the reactant 2a which is in accordance with the highly endoergic character of this reaction. The transition state much more resembles the product than the substrate which is in agreement with the Hammond's postulate. The energetic parameters of the reaction are summarized in Table 1. Only the relative energetic parameters are given here, the absolute values are available as Supplementary Information.
Regarding the thermodynamics, it follows from Table 1 that the first step of ring closure is an endothermic and endoergic process. The Gibbs free energy is even more positive than enthalpy because of the entropy loss during the process of ring closure-the number of degrees of freedom of the molecule is substantially reduced. With the carbinolamine less stable by almost 10 kcal/mol, it follows that it is created in very small amounts-the reaction equilibrium constant for this reaction is very low. In respect to kinetics, the calculated activation barrier is high enough-about 35 kcal/ mol which results in limited speed of this reaction. Interestingly, the difference between gas phase conditions and water environment is not large. In the later, the activation barrier is slightly reduced, but the product is slightly less stable thermodynamically.
To verify that the transition state really connects the substrate and product IRC (intrinsic reaction coordinate) calculations were performed at the B3LYP/6-31G(d) level of theory. The results are shown in Fig. 3 below.
In the middle of Fig. 3, there is a maximum which is the transition state-2b, if we go to the left we go towards the open form of midazolam-2a, and to the right we go towards the carbinolamine-2c. It follows that the energy decreases much faster towards 2c as the carbinolamine is a covalent structure and we quickly reach the minimum. The situation is different on the left side where we have an open ring weakly bonded by electrostatic interactions between the carbonyl carbon atom and the amino nitrogen atom. The energy decrease is much less steep as the carbon and nitrogen atoms slowly increase their distance minimizing the energy of the molecule. To reach the situation where the open form is lower energy than carbinolamine would require much more steps of IRC calculations and is not shown here. But finally, the IRC calculations proved that the reactants are truly 2a and 2c. The transition state optimized at B3LYP/6-31G(d) level and structures from beginning and ending of IRC are given in the Supplementary Information.

The protonation of the carbinolamine and the water molecule loss
The next essential step in the formation of final midazolam molecule is the protonation of the hydroxy group of carbinolamine to facilitate the removal of the water molecule. In the idealized textbook image this is 6B transition product (Scheme 6). However, all attempts to obtain an energetic minimum of such a structure failed, because the water molecules immediately is expelled from the protonated carbinolamine molecule. We tried to perform a relaxed PES scan gradually approaching the H 3 O + ion to the carbinolamine but the energy during the entire scan lowers and at certain distance the proton suddenly jumps to carbinolamine with immediately expelling the water molecule with further Table 1 Relative values of total electronic energy, enthalpy, and Gibbs free energy of molecules involved in carbinolamine formation, optimized at B3LYP/6-311 + + G(d,p) in gas phase, and water environment (PCM) for the first step of ring closure of midazolam lowering of the energy. This shows that the process of protonation of carbinolamine is activation barrierless and there is no transition state involved. This reaction would undergo much faster in acidic environments, but as we must remember, the first step-the formation of carbinolamine demanded neutral of basic conditions. However, even in such conditions, some hydronic ions exist in the water environment, so we can expect that protonation of carbinolamine occurs also in these conditions to some degree, which should be sufficient, because after protonation-it turns out-the protonated carbinolamine immediately breaks up into protonated imine and water molecule. Thus, it is not possible to obtain a stable transition product-protonated carbinolamine, except when we fix the length of the weak C-O bond at the value observed in neutral carbinolamine which is 1.450 Å for water environment (and 1.454 Å for the gas phase). We decided to use the distance of 1.45 Å for fixed optimizations in both environments. The process of optimizing the geometrical parameters of a molecule is called constrained optimization. The result 4a is shown in Fig. 4. The geometry of structure 4a differs to some extend in regard to the structure 2c, for example, the C-N bond in the benzodiazepine ring is shortened from 1.44 to 1.42 Å. But in general, their geometrical shapes are very similar, and we can conclude that the protonation does not change much the geometry of the carbinolamine. If the C-O bond length is not fixed, there undergoes immediately the loss of water molecule from the protonated carbinolamine which results in protonated imine and water molecule complex 4b in Fig. 4. It is needless to say that the process of expelling the water molecule is activation barrierless and there is no transition state. Comparison of structures 4a and 4b shows that the C-N bond now becomes clearly double, with its length equal to 1.3 Å. The relative energetic data for 4a and 4b are shown in Table 2.  In water environment, the enthalpy of 4b is lower by 50 kcal/mol and the Gibbs free energy by 53 kcal/mol than 4a. In gas phase, these values are even more negative by almost 10 kcal/mol. This is an illustration of great instability of 3a (if the C-O bond is not fixed). This reaction is connected with the gain of entropy, thus the more negative Gibbs free energy than enthalpy.

Acid-base equilibrium between two possible sites of protonation of MDZ
After the formation of the protonated imine, the rest is actually the question of acid-base equilibrium of MDZ molecule. As was mentioned in the introduction, there exist two possible sites of protonation-the benzodiazepine ring nitrogen atom and the imidazole ring nitrogen atom. The first possibility was already shown in Fig. 4b. In Fig. 5, there are shown two possible tautomers of protonated MDZ (without external water molecule).
The Gibbs free energy difference between 5a and 5b is 1.1 kcal/mol for gas phase and 3.4 kcal/mol for water environment in favor of 5b at the B3LYP/6-311 + + G(d,p) level of theory; thus, it is expected that the initially created form 5a will quickly tautomerize to more stable form 5b, especially in water environment. Let us look at the aromaticity of 5-and 7-membered rings for 5a and 5b, the data is gathered in Table 3.
It follows from Table 3 that in midazolam cation protonated at 7-membered ring the aromaticity of this ring increases from HOMA7 = − 0.30 to − 0.10 but aromaticity of 5-membered ring is slightly reduced from HOMA5 = 0.79 to 0.77. Also, the pi-deficiency of the 5-membered ring is slightly increased from pEDA5 = − 0.07 to − 0.1. On the contrary, in case of protonation at the 5-membered ring HOMA7 = − 0.36, so it is slightly reduced, comparing to midazolam, but HOMA5 = 0.81 which means slight increase. Also, pEDA5 is now positive-the ring is little pi-excessive. These changes are all rather small. We can look at it also from a more visual perspective-the pattern of single and double bonds shows a better bond conjugation and equalization of the 5-membered ring in the case of MDZ protonated at imidazole ( Supplementary Information-Figure S1). Similar picture of better bond equalization cannot be seen in case of protonation of 7-membered ring. Possible sources of this include (1) the methylene group which breaks the conjugation of bonds and (2) the fused benzene ring interferes with bond equalization. Thus, the protonation at the imidazole nitrogen atom is much more beneficial.

Deprotonation of the protonated imine to MDZ
The final step in our chain of reactions is the deprotonation of the protonated imine to actual closed-ring midazolam MDZ. This process should be easiest of course in basic conditions where the OH − ions are available. To estimate if the whole chain of reactions is endothermic or exothermic, we have to guarantee that the number of atoms and electron is exactly the same for the molecules which we  compare. Otherwise, any comparisons of the energy of two molecules have no physical sense. Fortunately, we can achieve this by comparing the MDZo with MDZ complex with a water molecule. In this case, there is exactly the same number of atoms and electrons as the whole reaction is just an amine-carbonyl condensation. Also, the charge of the molecules is the same-both are neutral. Thus, we need the MDZ molecule with a water molecule attached to it. There are two possible comparisons: of water molecule attachment: (1) to benzodiazepine ring nitrogen atom and (2) to imidazole ring nitrogen atom. The first possibility is a "natural" one as we can expect that when OH − ion is close to 4a, the proton will be removed and the water will be attached to the 7-membered ring. However, because of the possible tautomeric equilibrium described in the previous chapter and for the sake of completeness, we should also make the comparison with MDZ with water molecule attached to imidazole ring. Both possibilities are presented in Fig. 6.
The Gibbs free energy difference between 6a and 6b is 1.3 kcal/mol for gas phase and 1.9 kcal/mol for water environment at the B3LYP/6-311 + + G(d,p) level of theory, and the more stable molecule is again 6b. Lower energy of 6b is also reflected in stronger hydrogen bond between water molecule and MDZ, a rough measure if this strength is the bond length, which is a little shorter in case of 6b (see Fig. 6). The comparison of energetic parameters of the first molecule-midazolam open (benzophenone) and final molecule-midazolam closed (+ water molecule) shows the thermodynamics of the whole chain of reactions. The data is summarized in Table 4.
It follows that the whole chain of reactions is exothermic and exoergic, with similar enthalpy and Gibbs free energy difference, no matter to which final product we make the comparison. There is one subtle difference-in water environment if we make the comparison to 6a, there is neither entropy gain or loss because the loss in the formation of carbinolamine was compensated by the gain in the process  of water secretion. In the case of comparison to 6b, the ΔG is about 1 kcal/mol lower than ΔH . In the gas phase, the Gibs free energy for both comparisons is lower by about 0.5 kcal/ mol. We can conclude that from the thermodynamic point of view, the reaction with progress completely to the closed form of midazolam-the reaction equilibrium constant for the whole chain of reactions is high. From the kinetic point of view, the reaction rate limiting stage is the first step-the formation of carbinolamine with its quite high activation barrier-over 30 kcal/mol. The rest of the process undergoes fast without activation barriers.

The model of midazolam
In order to understand what is the impact of other fragments of the molecule on the thermodynamics and kinetics of benzodiazepine ring closure, we constructed a model where all other parts of the molecule except the imidazole ring were removed and replaced by hydrogen atoms. The model is shown in Fig. 7 and will be called MMDZ.
Regarding the aromaticity of the model, the indices are: HOMA7 = − 0.07, HOMA5 = 0.84, and pEDA5 = − 0.005. It follows that both rings are more aromatic than in midazolam itself. Thus, it follows that other parts of the molecule decrease the aromaticity of the core rings.

Formation of carbinolamine
The first step was to model the carbinolamine molecules formation for which reactant, transition state, and product are shown in Fig. 8.
From a quick comparison of Figs. 8 and 2, it follows that there are significant differences in geometries of MDZo and MMDZo and transition states. The carbinolamines show a very similar geometry. Regarding the open forms the C-N distance is much smaller in the model (2.9 Å comparing to 3.3 Å). It is probably connected to simpler structure of the model with less steric repulsions, allowing for more flexibility of the chain of atoms in the model. In the transition state, all essential bond lengths (except C-N) are longer in the model. The most interesting is however the energetic data which is gathered in Table 5.
When comparing Table 5 to Table 1, the most striking fact is that the process of carbinolamine formation is exothermic and exoergic in the case of the model of midazolam. A possible explanation of this fact is that for midazolam, the closure of the ring is connected with many unfavorable steric interactions. The open form of midazolam has very complicated spatial structure, its conformation is heavily bent in several directions to allow for the energy lowering. During the ring closure many bonds have to rotate and the closed carbinolamine puts many constrains on the  structure. On the other hand, the structure of the model is much simpler and its closure does not involve these unfavorable interactions. The IRC profile of the model is very similar to the midazolam itself, and thus was not shown here.

The protonation of carbinolamine and the water molecule loss
The protonated carbinolamine with C-O bond length fixed at 1.446 (the value for optimized model carbinolamine) is unstable, like in the case of MDZ. The important difference is that in the case of model, it was not possible to obtain a true minimum, but only the first-order saddle point. It follows from comparison of Figs. 9 and 4, that the C-O and C-N bonds are shorter in the model. Without constrains, the structure 9a immediately loses the water molecule during optimization leading to 9b.
The energetic data for the water loss process are gathered in Table 6.
The difference between 9b and 9a is smaller than in the case of 4b and 4a by about 10 kcal/mol. It can be explained by the nature of 9a which is a first-order saddle point, thus is richer in energy than a minimum.

Acid-base equilibrium between two possible sites of protonation of MMDZ
As in the case of MDZ, there exist two possible tautomeric forms of MMDZ cation and they are showed in Fig. 10.
It turns out, however, that the order of stability of these form is reversed. The molecule protonated at the 7-membered ring nitrogen atom-10a has lower Gibbs free energy. In gas phase, the difference is 1.1 kcal/mol and in water environment it is increased up to 2.8 kcal/mol. In order to better understand these data, let us look at the aromaticity indices gathered in Table 7.
Comparing the situation of protonation at 7-membered ring and 5-membered ring, we see that aromaticity of 7-membered ring significantly increases in the first case and significantly reduces in second case. For the 5-membered ring, the situation is reversed, but the changes of aromaticity are much smaller. The pEDA5 = − 0.30 for 10a shows some pi-electron deficiency of the 5-membered ring and means that the pi-electron density is shifted to 7-membered ring. If we look at the bond equalization scheme ( Figure S2 -Supplementary Information), we can see a conjugated, butadiene-like fragment in the 7-membered ring in case of 10a. The fused benzene ring is not present in the model; thus, it does not interfere with the bond equalization pattern and this increases the aromaticity of the 7-membered ring, thus stabilizes the cation protonated at the 7-membered ring. This ring is less nonplanar than analogous ring in the midazolam cation.

Deprotonation of the protonated imine to MDZ
Finally, let us compare the open-ring form of midazolam model to its closed form with a water molecules attached, again to guarantee the equal number of atoms for energy comparisons. The neutral form MMDZ with a water molecule attached are showed in Fig. 11. The energetic situation becomes now very interesting as both molecules 11a and 11b are very close in energy. In gas phase, the 11a has slightly lower energy and in water environment 11b but the difference is less that 0.5 kcal/mol, so it is below the 1 kcal/mol error and is not meaningful. We can assume that 11a and 11b have very similar stability. It follows from Fig. 11 that the hydrogen bond length (connected to its strength) is only a little shorter for 11a. On the contrary, in Fig. 6, this bond for 6b is much shorter than 6a which is connected with better stability of 6b. Thus, a comparison of the final product-midazolam model with water molecule to midazolam model open will be done only for 11a.   It follows from Table 8 that the whole ring closure process for the model is even more exothermic and exoergic than for the midazolam itself. It can be attributed to less congested structure of the model of the midazolam for which the process of ring closure not connected to any unfavorable steric interactions.

Conclusion
1. The mechanism of midazolam ring closure was analyzed from thermodynamic and kinetic point of view at the B3LYP/6-311 + + G(d,p) level of theory. 2. The reaction rate constant is determined by the Gibbs free energy of activation of the first step-carbinolamine formation. It amounts to ΔG # = 35.1 kcal/mol for gas phase conditions and ΔG # = 33.9 kcal/mol for water environment. IRC calculations proved that the transition state really connects the expected substrate and product. 3. The process of carbinolamine formation for midazolam is endoergic with ΔG = 9.6 kcal/mol for gas phase and ΔG = 10.0 kcal/mol for water environment. 4. Further reactions are carbinolamine protonation and immediate water loss with creation of the imine cation (which is the closed ring for of midazolam in protonated form). These reactions undergo without activations barriers and do not affect the kinetics of the process. 5. The whole series of reactions leading to ring closure is exoergic with ΔG = −5.6 kcal/mol for gas phase and ΔG = −7.7 kcal/mol for water environment. This fact and point 4 are the reasons why the ring closure is possible despite the endoergic character of the first step. 6. To understand the influence of the other than benzodiazepine and imidazole rings to ring closure process a model was built where these parts of the molecule were replaced by hydrogen atoms. 7. The main difference between midazolam and its simplistic model was that for model the fist step-carbinolamine formation-is exoergic with ΔG = −2.2 kcal/mol for gas phase and − 1.8 kcal/mol for water environment, despite large entropy loss connected with the ring closure. The reason for this difference is that the midazolam molecule has complicated and very nonplanar shape which during the ring closure is the source of unfavourable interactions and rise of the energy of the molecule. 8. The ΔG # for carbinolamine formation for the model is lower for gas phase ( ΔG # = 33.3 kcal/mol) but higher for water environment ( ΔG # = 41.6 kcal/mol) comparing to midazolam. 9. The whole series of reactions for the model is more exoergic than for midazolam with ΔG = −10.5 kcal/ mol for the gas phase and ΔG = −11.5 kcal/mol for water environment. 10. The fragments of midazolam which differs it from its simplistic model have the influence of slowing down the process of ring closure mainly by the endoergic character of the first step. 11. Protonation of the midazolam occurs mainly at the imidazole nitrogen atom with ΔG difference of 1.1 kcal/mol for gas phase and 3.4 kcal/mol in water environment. 12. On the contrary, the midazolam model prefers protonation at the benzodiazepine nitrogen atom with ΔG difference of 1.1 kcal/mol and 2.8 kcal/mol in water environment. 13. Points 11 and 12 can be explained taking into consideration the aromaticity of the rings. In the case of protonation at the benzodiazepine ring, better bond equalization scheme in the 7-membered ring can be observed for midazolam model a than for midazolam. 14. Complexes with a water molecule behave differently for midazolam and its model. In the case of midazolam, the water molecule is preferentially connected to imidazole ring and the ΔG difference is 1.3 kcal/mol for gas phase and 1.9 kcal/mol for water environment. For midazolam model, the two complexes have similar stability.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s11224-022-02108-6. Data availability Additional data (total electronic energies, enthalpies, Gibbs free energies, and geometric parameters of optimized molecules for both gas phase and water environment, GaussView pictures of selected molecules for single/double bonds patters scheme) are available as Supplementary Information.

Competing interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.