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

doi:10.21203/rs.3.rs-2201755/v1

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. 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 unfavourable 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 analysed using two aromaticity indices – 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 analysed 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.

midazolam

open form

benzophenone

benzodiazepine

transition state

Gibbs free energy

B3LYP

aromaticity

HOMA

pEDA

acid

base

equilibrium

tautomerism

ring opening

ring closure

condensation

Benzodiazepines¹ are amongst 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² 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 class – midazolam, contains imidazole ring³ and others contain mostly the 1,2,4-triazole ring: alprazolam⁴, estazolam⁵ and triazolam⁶.

In this paper we will focus on midazolam (MDZ) as a typical and widely used drug for premedication, induction of anaesthesia 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⁷. 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 minutes.

Because of forming of these various ionic forms, MDZ was the first soluble benzodiazepine⁸. 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¹¹.

The process of ring opening is fully reversible and the closed-ring 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, were measured by UV spectrometry¹². 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 analyse 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¹³ 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₂⁺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₂ 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₂ and CO groups. We obtained the proper transition state for this stage, which will be shown and analysed 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 quantum-chemical 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.

All calculations were done by using the Gaussian 16 suite of programs¹⁴ with the use of universal B3LYP density functional¹⁵,¹⁶. Preliminary geometry optimizations were performed with the use of Pople’s 6-31G(d) orbital basis set¹⁷. Next, to obtain more reliable structures and energies, the more extended 6-311++G(d,p) basis set¹⁸,¹⁹ was employed. All optimizations were performed for gas-phase and water solution modelled via Polarizable Continuum Model (PCM)²⁰ 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. Regarding the energetic parameters, total electronic energy, enthalpy and Gibbs free energy of the reaction were analysed. 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 analysed using geometric HOMA²¹,²² and electronic pEDA²³,²⁴ aromaticity indices, calculated by using the free AromaTcl²⁵ software.

Midazolam

In the text the energetic data for B3LYP/6-311++G(d,p) results, both for the gas-phase and water environment modelled 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 Figure 1. All further Figures are done at this model chemistry.

From comparison of Scheme 1B and Figure 1 it follows that the molecule is in reality highly nonplanar. The methylene group is mostly responsible for this, it also breaks the pi-electron 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 analyse 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 analysed.

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 Figure 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 Figure 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.

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

	Gas-phase			Water environment (PCM)
	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]
2a	0	0	0	0	0	0
2b (TS)	33.84	31.38	35.08	32.64	30.22	33.93
2c	5.83	6.19	9.63	6.18	6.54	9.98

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.

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₃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 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 3a is shown in Figure 3.

The geometry of structure 3a 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 3b in Figure 3. 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 3a and 3b shows that the C-N bond now becomes clearly double, with its length equal to 1.3 Å. The relative energetic data for 3a and 3b are shown in Table 2.

Table 2. Relative values of total electronic energy, enthalpy and Gibbs free energy of protonated carbinolamine with fixed C-N bond length, optimized at B3LYP/6-311++G(d,p) in gas-phase and water environment (PCM)

	Gas-phase			Water environment (PCM)
	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]
3a	0	0	0	0	0	0
3b	-58.97	-59.22	-62.15	-48.95	-49.80	-52.90

In water environment the enthalpy of 3b is lower by 50 kcal/mol and the Gibbs free energy by 53 kcal/mol than 3a. 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 Figure 3b. In the Figure 4 there are shown two possible tautomers of protonated MDZ (without external water molecule)

The Gibbs free energy difference between 4a and 4b is 1.1 kcal/mol for Gas Phase and 3.4 kcal/mol for water environment in favour of 4b at the B3LYP/6-311++G(d,p) level of theory, thus it is expected that the initially created form 4a will quickly tautomerize to more stable form 4b, especially in water environment. Let us look at the aromaticity of 5- and 7-membered rings for 4a and 4b, the data is gathered in Table 3.

Table 3. HOMA for 7- and 5-membered rings and pEDA for 5-membered ring for midazolam (MDZ) protonated at 7- and 5-membered ring.

	HOMA7	HOMA5	pEDA5
4a	-0.098	0.774	-0.100
4b	-0.362	0.811	0.031

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 on OH^- ion is close to 3b, the proton will be removed and the water will be attached to the seven-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.

The Gibbs free energy difference between 5a and 5b 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 5b. Lower energy of 5b 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 5b (see Figure 5). 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.

Table 4. Relative values of total electronic energy, enthalpy and Gibbs free energy of closed form midazolam with water molecule attached to nitrogen atom of 7-membered and 5-membered ring in respect to open form, optimized at B3LYP/6-311++G(d,p) in gas-phase and in water environment (PCM)

	Gas-phase			Water environment (PCM)
	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]
MDZo	0	0	0	0	0	0
5a	-4.66	-5.12	-5.64	-7.18	-7.72	-7.67
5b	-5.96	-6.42	-6.94	-7.91	-8.55	-9.53

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 5a, 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 5b the is about 1 kcal/mol lower than . 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 Figure 6 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 Figure 7.

From a quick comparison of Figures 7 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.

Table 5. 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 the gas-phase and water environment (PCM) for midazolam model

	Gas-phase			Water environment (PCM)
	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]
7a	0	0	0	0	0	0
7b (TS)	32.44	30.27	33.26	40.48	38.62	41.59
7c	-6.15	-4.91	-2.20	-5.48	-4.35	-1.83

When comparing Table 4 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 unfavourable 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 unfavourable interactions.

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 Figures 8 and 3 that the C-O and C-N bonds are shorter in the model. Without constrains the structure 8a immediately loses the water molecule during optimization leading to 8b.

The energetic data for the water loss process are gathered in Table 6.

Table 6. Relative values of total electronic energy, enthalpy and Gibbs free energy of protonated carbinolamine with fixed C-N bond length, optimized at B3LYP/6-311++G(d,p) in the gas-phase and in water environment (PCM)

	Gas-phase			Water environment (PCM)
	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]	∆H [kcal/mol]	∆E [kcal/mol]	∆G [kcal/mol]
8a	0	0	0	0	0	0
8b	-50.03	-50.28	-54.38	-38.96	-39.65	-44.06

The difference between 8b and 8a is smaller than in the case of 3b and 3a by about 10 kcal/mol. It can be explained by the nature of 8a 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 Figure 9.

It turns out however that the order of stability of these form is reversed. The molecule protonated at the 7-membered ring nitrogen atom - 9a 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 3.

Table 3. HOMA for 7- and 5-membered rings and pEDA for 5-membered ring of midazolam model protonated at 7- and 5-membered ring.

	HOMA7	HOMA5	pEDA5
9a	0.242	0.645	-0.296
9b	-0.201	0.841	-0.005

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 9a 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 9a. 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 non-planar 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 Figure 10.

The energetic situation becomes now very interesting as both molecules 10a and 10b are very close in energy. In gas-phase the 10a has slightly lower energy and in water environment 10b 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 10a and 10b have very similar stability. It follows from Figure 10 that the hydrogen bond length (connected to its strength) is only a little shorter for 10a. On the contrary in Figure 5 this bond for 5b is much shorter than 5a which is connected with better stability of 5b. Thus, a comparison of the final product – midazolam model with water molecule to midazolam model open will be done only for 10a.

Table 7. Relative values of total electronic energy, enthalpy and Gibbs free energy of closed form midazolam with water molecule attached to nitrogen atom of 7-membered ring in respect to open form, optimized at B3LYP/6-311++G(d,p) in the gas-phase and in water environment (PCM)

	Gas-phase			Water environment (PCM)
	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]	∆E [kcal/mol]	∆H [kcal/mol]	∆G [kcal/mol]
MMDZo	0	0	0	0	0	0
10a	-8.67	-8.96	-10.54	-9.88	-10.33	-11.53

It follows from Table 7 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 unfavourable steric interactions.

The mechanism of midazolam ring closure was analysed from thermodynamic and kinetic point of view at the B3LYP/6-311 + + G(d,p) level of theory.

The reaction rate constant is determined by the Gibbs free energy of activation of the first step – carbinolamine formation. It amounts to \({\Delta }{G}^{\#}=35.1\) kcal/mol for gas-phase conditions and \({\Delta }{G}^{\#}=33.9\) kcal/mol for water environment

The process of carbinolamine formation for midazolam is endoergic with \({\Delta }G=9.6\) kcal/mol for gas-phase and \({\Delta }G=10.0\) kcal/mol for water environment.

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.

The whole series of reactions leading to ring closure is exoergic with \({\Delta }G=-5.6\) kcal/mol for gas-phase and \({\Delta }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.

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.

The main difference between midazolam and its simplistic model was that for model the fist step – carbinolamine formation – is exoergic with \({\Delta }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.

The \({\Delta }{G}^{\#}\) for carbinolamine formation for the model is lower for gas-phase (\({\Delta }{G}^{\#}=33.3\) kcal/mol) but higher for water environment (\({\Delta }{G}^{\#}=41.6\) kcal/mol) comparing to midazolam.

The whole series of reactions for the model is more exoergic than for midazolam with \({\Delta }G=-10.5\) kcal/mol for the gas-phase and \({\Delta }G=-11.5\) kcal/mol for water environment.

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.

Protonation of the midazolam occurs mainly at the imidazole nitrogen atom with \({\Delta }G\) difference of 1.1 kcal/mol for gas-phase and 3.4 kcal/mol in water environment

On the contrary the midazolam model prefers protonation at the benzodiazepine nitrogen atom with \({\Delta }G\) difference of 1.1 kcal/mol and 2.8 kcal/mol in water environment

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.

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 \({\Delta }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.

ACKNOWLEDGEMENTS

Computational Grant G36-9 from the Interdisciplinary Centre for Mathematical and Computational Modelling at Warsaw University (ICM UW) is gratefully acknowledged

Funding

Not applicable

Conflicts of interest/Competing interests

The authors declare that they have no conflict of interest

Availability of data and material

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

Code availability

Not applicable

Authors' contributions

Not applicable

Sellers, EM (1978) Clinical pharmacology and therapeutics of benzodiazepines. Can Med Assoc J 118(12):1533–1538
Calcaterra NE, Barrow JC (2014) Classics in chemical neuroscience: diazepam (valium). ACS Chemical Neuroscience 5 (4): 253–60
Olkkola KT, Ahonen J (2008) Midazolam and other benzodiazepines, Handb Exp Pharmacol 182: 335–360
Chouinard G, Annable L, Fontaine R, Solyom L (1982) Alprazolam in the treatment of generalized anxiety and panic disorders: a double-blind placebo-controlled study. Psychopharmacology (Berl.), 77(3): 229–233
Roehrs T, Zorick F, Lord N, Koshorek GL, Roth T (1983). Dose-related effects of estazolam on sleep of patients with insomnia. J Clin Psychopharmacol 3(3): 152–156
Pakes GE, Brogden RN, Heel RC, Speight TM, Avery GS (1981) Triazolam: a review of its pharmacological properties and therapeutic efficacy in patients with insomnia. Drugs 22(2): 81–110
Gerecke M (1983) Chemical structure and properties of midazolam compared with other benzodiazepines. Br J Clin Pharmacol 16(1): 11S-16S
Kanto JH (1985) Midazolam: the first water-soluble benzodiazepine. Pharmacology, pharmacokinetics and efficacy in insomnia and anesthesia. Pharmacotherapy 5(3): 138–155
Ota Y, Hamada A, Nakano M, Saito H (2003) Evaluation of Percutaneous Absorption of Midazolam by Terpenes. Drug Metab Pharmacokinet 18(4): 261-266
Kanto, J (1983) Drugs Today 19: 221-227
Orive MM, Gallo B, Alonso R, Vicente F, Vire JC, Patriarce GJ (1989) Spectrophotometric Study of the Acid-Base Equilibria of an Imidazobenzodiazepine, Midazolam, and Its Determination in Pharmaceutical Formulations. Mikrochimica Acta 1: 181-190
Andersin R (1991) Solubility and acid-base behaviour of midazolam in media of different pH, studied by ultraviolet spectrophotometry with multicomponent software. J Pharm Biomed Anal 9(6):451-455
McMurry J, Organic Chemistry 9th edition (2016), Cengage Learning
Gaussian 16 Revision B.01 Frisch MJ, Trucks GW, Schlegel HB,Scuseria GE, Robb MA, Cheeseman JR, Scalmani G, Barone V,Petersson GA, Nakatsuji H , Li X, Caricato M, Marenich AV, Bloino J, Janesko BG, Gomperts R, Mennucci B, Hratchian HP,Ortiz JV, Izmaylov AF, Sonnenberg JL , Williams-Young D, Ding F, Lipparini F, Egidi F, Goings J, Peng B, Petrone A, Henderson T, Ranasinghe D, Zakrzewski VG, Gao J, Rega N, Zheng G, Liang W, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Vreven T, Throssell K, Montgomery Jr JA, Peralta JE, Ogliaro F, Bearpark M J, Heyd JJ, Brothers EN, Kudin KN, Staroverov VN, Keith TA, Kobayashi R, Normand J, Raghavachari K, Rendell AP, Burant JC, Iyengar SS, Tomasi J, Cossi M, Millam JM, Klene M, Adamo C, Cammi R, Ochterski JW, Martin RL, Morokuma K, Farkas O, Foresman JB, Fox DJ. Gaussian Inc. Wallingford CT 2016
Becke AD (1988) Density-functional exchange-energy approximation with correct asymptotic behavior. Phys Rev A 38(6):3098
Lee C, Yang W, Parr RG (1988) Development of the Colle-Salvetticorrelation-energy formula into a functional of the electron density. Phys Rev B 37(2):785
Petersson GA, Bennett A, Tensfeldt TG, Al-Laham MA, Shirley WA, Mantzaris J (1988) A complete basis set model chemistry. I. The total energies of closed-shell atoms and hydrides of the first row elements. J Chem Phys 89(4):2193
McLean AD, Chandler GS (1980) Contracted Gaussian basis sets for molecular calculations. I. Second row atoms, Z=11–18. J Chem Phys 72(10):5639–5648
Clark T, Chandrasekhar J, Spitznagel GW, Schleyer PR (1983) Efficient diffuse function-augmented basis sets for anion calculations. III. The 3–21+G basis set for first-row elements. Li–F. J Comput Chem 4(3):294–301
Scalmani G, Frisch MJ (2010) Continuous surface charge polarizable continuum models of solvation. I General formalism J Chem Phys 132:114110–114115
Kruszewski J, Krygowski TM (1972) Definition of aromaticity basing on the harmonic oscillator model. Tetrahed Lett 13:3839
Krygowski TM (1993) Crystallographic studies of inter- and intramolecular interactions reflected in aromatic character of .pi.-electron systems. J Chem Inf Comput Sci 33:70
Oziminski WP, Dobrowolski JC (2009) σ- and π-electron contributions to the substituent effect: natural population analysis. J Phys Org Chem 22:769
Oziminski WP, Krygowski TM (2011) Effect of aromatization of the ring on intramolecular H-bond in 3-hydroxy-4-formylo derivatives of fulvene. Chem Phys Lett 510: 53-56
Oziminski WP (2020) AromaTcl: a program for calculation of sEDA, pEDA and HOMA indices. Internet:https://www.researchgate.net/project/AromaTcl-a-program-for-calculation-of-sEDA-pEDA-and-HOMA-indices

Scheme 1 to 6 are available in the Supplementary Files section.

No competing interests reported.

SIMidazolamRingClosingFINAL.docx
Scheme1.png
Scheme 1. Schematic structures of: (A) Diazepam and (B) Midazolam
Scheme2.png
Scheme 2. Schematic structures of midazolam ions: (A) cation (MDZ+) and (B) dication (MDZ++)
Scheme3.png
Scheme 3. Fully protonated open (benzophenone) form of midazolam in highly acidic medium
Scheme4.png
Scheme 4. Two neutral forms of midazolam: (A) open - MDZo and (B) closed - MDZ
Scheme5.png
Scheme 5. An idealized (and non-existent) bipolar transition product of the first step of ring closure
Scheme6.png
Scheme 6. Three postulated transition products in the process of ring closure: (A) Carbinolamine (B) protonated carbinolamine (C) protonated imine

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

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Computational Details

Results And Discussion

Conclusions

Declarations

References

Scheme

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1