Category Archives: NBO

Aurides Chemistry – New Paper in Organometallics


Compound 2 represents the first structural example of a 12 e− auride complex, with a pseudohalide/hydride nature in bonding. According to our NBO calculations, this electron deficient gold center is stabilized by weak intramolecular interactions between Au p orbitals and σC−C and σC−H bonds of adjacent aromatic rings together with a Ga−Au−Ga 3 centers−2 electrons bond (I like the term ‘banana bond‘, don’t you?).

Fig. 1 Crystal structure for Compound 2. Au in the center is effectively an auride.

I was invited to participate in this wonderful venture by my good friend and colleague Dr. José Oscar Carlos Jiménez-Halla, from the University of Guanajuato, Mexico, with whom we’re now working with Prof. Rong Shang at the Hiroshima University. Prof. Shang has synthesized this portentous Auride complex and over the last year, Leonardo “Leo” Lugo has worked with Oscar and I in calculating their electronic structure and bonding properties.

Gold catalysis is an active area of research but low valent Au compounds are electron deficient and therefore highly reactive and elusive; that’s why researchers prefer to synthesize these compounds in situ, to harness their catalytic properties before they’re lost. Power’s digalladeltacyclane was used as a ligand framework to bind to a Au(I) center, which became reduced after the addition and breaking of the Ga−Ga bond while the opposite face of the metallic center became blocked by the bulky aromatic groups on the main ligand. NBO calculations at the M05-2X/[LANL2TZ(f),6-311G(d,p)] and QTAIM BCP analysis show the main features of Au bonding in 2, noteworthy features are the 3c−2e bond (banana) and the σC−C and σC−H donations (See figure 2).

Fig.2 Natural Hybrid Composition for the Ga−Au−Ga ‘banana‘ bond (left). Bond Critical Points (BCPs) for Au in 2 (right).

One of the most interesting features of this compound is the fact that Au(PPh3)Cl reacts differently to the digallane ligand than it does to analogous B−B, Si−Si, or Sn−Sn bonds. The Au−Cl bond does not undergo metathesis as with B−B, nor does it undergo an oxidative addition, so to further understand the chemistry of−and leading to−compound 2, the reaction mechanism energy profile was calculated in a rather painstakingly effort (Kudos, Leo, and a big shoutout to my friend Dr. Jacinto Sandoval for his one on one assistance). Figure 3 shows the energy profile for the reaction mechanism for the formation of 2 from Power’s digallane reagent and Au(PPh3)Cl.

Fig. 3 Free Energy profile for the formation of 2. All values, kcal/mol

You can read more details about this research in Organometallics DOI:10.1021/acs.organomet.0c00557. Thanks again to Profs. Rong Shang and Óscar Jiménez-Halla for bringing me on board of this project and to Leo for his relentless work getting those NBO calculations done; this is certainly the beginning of a golden opportunity for us to collaborate on a remarkable field of chemistry, it has certainly made me go bananas over Aurides chemistry. OK I’ll see myself out.

Density Keyword in Excited State Calculations with Gaussian


I have written about extracting information from excited state calculations but an important consideration when analyzing the results is the proper use of the keyword density.

This keyword let’s Gaussian know which density is to be used in calculating some results. An important property to be calculated when dealing with excited states is the change in dipole moment between the ground state and any given state. The Transition Dipole Moment is an important quantity that allows us to predict whether any given electronic transition will be allowed or not. A change in the dipole moment (i.e. non-zero) of a molecule during an electronic transition helps us characterize said transition.

Say you perform a TD-DFT calculation without the density keyword, the default will provide results on the lowest excited state from all the requested states, which may or may not be the state of interest to the transition of interest; you may be interested in the dipole moment of all your excited states.

Three separate calculations would be required to calculate the change of dipole moment upon an electronic transition:

1) A regular DFT for the ground state as a reference
2) TD-DFT, to calculate the electronic transitions; request as many states as you need/want, analyze it and from there you can see which transition is the most important.
3) Request the density of the Nth state of interest to be recovered from the checkpoint file with the following route section:

# TD(Read,Root=N) LOT Density=Current Guess=Read Geom=AllCheck

replace N for the Nth state which caught your eye in step number 2) and LOT for the Level of Theory you’ve been using in the previous steps. That should give you the dipole moment for the structure of the Nth excited state and you can compare it with the one in the ground state calculated in 1). Again, if density=current is not used, only properties of N=1 will be printed.

Orbital Contributions to Excited States


This is a guest post by our very own Gustavo “Gus” Mondragón whose work centers around the study of excited states chemistry of photosynthetic pigments.

When you’re calculating excited states (no matter the method you’re using, TD-DFT, CI-S(D), EOM-CCS(D)) the analysis of the orbital contributions to electronic transitions poses a challenge. In this post, I’m gonna guide you through the CI-singles excited states calculation and the analysis of the electronic transitions.

I’ll use adenine molecule for this post. After doing the corresponding geometry optimization by the method of your choice, you can do the excited states calculation. For this, I’ll use two methods: CI-Singles and TD-DFT.

The route section for the CI-Singles calculation looks as follows:

%chk=adenine.chk
%nprocshared=8
%mem=1Gb

#p CIS(NStates=10,singlets)/6-31G(d,p) geom=check guess=read scrf=(cpcm,solvent=water)

adenine excited states with CI-Singles method

0 1
--blank line--

I use the same geometry from the optimization step, and I request only for 10 singlet excited states. The CPCP implicit solvation model (solvent=water) is requested. If you want to do TD-DFT, the route section should look as follows:

%chk=adenine.chk
%nprocshared=8
%mem=1Gb

#p FUNCTIONAL/6-31G(d,p) TD(NStates=10,singlets) geom=check guess=read scrf=(cpcm,solvent=water)

adenine excited states with CI-Singles method

0 1
--blank line--

Where FUNCTIONAL is the DFT exchange-correlation functional of your choice. Here I strictly not recommend using B3LYP, but CAM-B3LYP is a noble choice to start.

Both calculations give to us the excited states information: excitation energy, oscillator strength (as f value), excitation wavelength and multiplicity:

Excitation energies and oscillator strengths:

 Excited State   1:      Singlet-A      6.3258 eV  196.00 nm  f=0.4830  <S**2>=0.000
      11 -> 39        -0.00130
      11 -> 42        -0.00129
      11 -> 43         0.00104
      11 -> 44        -0.00256
      11 -> 48         0.00129
      11 -> 49         0.00307
      11 -> 52        -0.00181
      11 -> 53         0.00100
      11 -> 57        -0.00167
      11 -> 59         0.00152
      11 -> 65         0.00177

The data below corresponds to all the electron transitions involved in this excited state. I have to cut all the electron transitions because there are a lot of them for all excited states. If you have done excited states calculations before, you realize that the HOMO-LUMO transition is always an important one, but not the only one to be considered. Here is when we calculate the Natural Transition Orbitals (NTO), by these orbitals we can analyze the electron transitions.

For the example, I’ll show you first the HOMO-LUMO transition in the first excited state of adenine. It appears in the long list as follows:

35 -> 36         0.65024

The 0.65024 value corresponds to the transition amplitude, but it doesn’t mean anything for excited state analysis. We must calculate the NTOs of an excited state from a new Gaussian input file, requesting from the checkpoint file we used to calculate excited states. The file looks as follows:

%Oldchk=adenine.chk
%chk=adNTO1.chk
%nproc=8
%mem=1Gb

#p SP geom=allcheck guess=(read,only) density=(Check,Transition=1) pop=(minimal,NTO,SaveNTO)

I want to say some important things right here for this last file. See that no level of theory is needed, all the calculation data is requested from the checkpoint file “adenine.chk”, and saved into the new checkpoint file “adNTO1.chk”, we must use the previous calculated density and specify the transition of interest, it means the excited state we want to analyze. As we don’t need to specify charge, multiplicity or even the comment line, this file finishes really fast.

After doing this last calculation, we use the new checkpoint file “adNTO1.chk” and we format it:

formchk -3 adNTO1.chk adNTO1.fchk

If we open this formatted checkpoint file with GaussView, chemcraft or the visualizer you want, we will see something interesting by watching he MOs diagram, as follows:

We can realize that frontier orbitals shows the same value of 0.88135, which means the real transition contribution to the first excited state. As these orbitals are contributing the most, we can plot them by using the cubegen routine:

cubegen 0 mo=homo adNTO1.fchk adHOMO.cub 0 h

This last command line is for plotting the equivalent as the HOMO orbital. If we want to plot he LUMO, just change the “homo” keyword for “lumo”, it doesn’t matter if it is written with capital letters or not.

You must realize that the Natural Transition Orbitals are quite different from Molecular Orbitals. For visual comparisson, I’ve printed also the molecular orbitals, given from the optimization and from excited states calculations, without calculating NTOs:

These are the molecular frontier orbitals, plotted with Chimera with 0.02 as the isovalue for both phase spaces:

The frontier NTOs look qualitatively the same, but that’s not necessarily always the case:

If we analyze these NTOs on a hole-electron model, the HOMO refers to the hole space and the LUMO refers to the electron space.

Maybe both orbitals look the same, but both frontier orbitals are quite different between them, and these last orbitals are the ones implied on first excited state of adenine. The electron transition will be reported as follows:

If I can do a graphic summary for this topic, it will be the next one:

NTOs analysis is useful no matter if you calculate excited states by using CIS(D), EOM-CCS(D), TD-DFT, CASSCF, or any of the excited states method of your election. These NTOs are useful for population analysis in excited states, but these calculations require another software, MultiWFN is an open-source code that allows you to do this analysis, and another one is called TheoDORE, which we’ll cover in a later post.

Natural Transition Orbitals (NTOs) Gaussian


The canonical molecular orbital depiction of an electronic transition is often a messy business in terms of a ‘chemical‘ interpretation of ‘which electrons‘ go from ‘which occupied orbitals‘ to ‘which virtual orbitals‘.

Natural Transition Orbitals provide a more intuitive picture of the orbitals, whether mixed or not, involved in any hole-particle excitation. This transformation is particularly useful when working with the excited states of molecules with extensively delocalized chromophores or multiple chromophoric sites. The elegance of the NTO method relies on its simplicity: separate unitary transformations are performed on the occupied and on the virtual set of orbitals in order to get a localized picture of the transition density matrix.

[1] R. L. Martin, J. Chem. Phys., 2003, DOI:10.1063/1.1558471.

In Gaussian09:
After running a TD-DFT calculation with the keyword TD(Nstates=n) (where n = number of states to be requested) we need to take that result and launch a new calculation for the NTOs but lets take it one step at a time. As an example here’s phenylalanine which was already optimized to a minimum at the B3LYP/6-31G(d,p) level of theory. If we take that geometry and launch a new calculation with the TD(Nstates=40) in the route section we obtain the UV-Vis spectra and the output looks like this (only the first three states are shown):

Excitation energies and oscillator strengths:

Excited State 1: Singlet-A 5.3875 eV 230.13 nm f=0.0015 <S**2>=0.000
42 -> 46 0.17123
42 -> 47 0.12277
43 -> 46 -0.40383
44 -> 45 0.50838
44 -> 47 0.11008
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -554.614073682
Copying the excited state density for this state as the 1-particle RhoCI density.

Excited State 2: Singlet-A 5.5137 eV 224.86 nm f=0.0138 <S**2>=0.000
41 -> 45 -0.20800
41 -> 47 0.24015
42 -> 45 0.32656
42 -> 46 0.10906
42 -> 47 -0.24401
43 -> 45 0.20598
43 -> 47 -0.14839
44 -> 45 -0.15344
44 -> 47 0.34182

Excited State 3: Singlet-A 5.9254 eV 209.24 nm f=0.0042 <S**2>=0.000
41 -> 45 0.11844
41 -> 47 -0.12539
42 -> 45 -0.10401
42 -> 47 0.16068
43 -> 45 -0.27532
43 -> 46 -0.11640
43 -> 47 0.16780
44 -> 45 -0.18555
44 -> 46 -0.29184
44 -> 47 0.43124

The oscillator strength is listed on each Excited State as “f” and it is a measure of the probability of that excitation to occur. If we look at the third one for this phenylalanine we see f=0.0042, a very low probability, but aside from that the following list shows what orbital transitions compose that excitation and with what energy, so the first line indicates a transition from orbital 41 (HOMO-3) to orbital 45 (LUMO); there are 10 such transitions composing that excitation, visualizing them all with canonical orbitals is not an intuitive picture, so lets try the NTO approach, we’re going to take excitation #10 for phenylalanine as an example just because it has a higher oscillation strength:

%chk=Excited State 10: Singlet-A 7.1048 eV 174.51 nm f=0.3651 <S**2>=0.000
41 -> 45 0.35347
41 -> 47 0.34685
42 -> 45 0.10215
42 -> 46 0.17248
42 -> 47 0.13523
43 -> 45 -0.26596
43 -> 47 -0.22995
44 -> 46 0.23277

Each set of NTOs for each transition must be calculated separately. First, copy you filename.chk file from the TD-DFT result to a new one and name it after the Nth state of interest as shown below (state 10 in this case). NOTE: In the route section, replace N with the number of the excitation of interest according to the results in filename.log. Run separately for each transition your interested in:

#chk=state10.chk

#p B3LYP/6-31G(d,p) Geom=AllCheck Guess=(Read,Only) Density=(Check,Transition=N) Pop=(Minimal,NTO,SaveNTO)

0 1
--blank line--

By requesting SaveNTO, the canonical orbitals in the state10.chk file are replaced with the NTOs for the 10th excitation, this makes it easier to plot since most visualizers just plot whatever set of orbitals they read in the chk file but if they find the canonical MOs then one would need to do some re-processing of them. This is much more straightforward.

Now we format our chk files into fchk with the formchk utility:

formchk -3 filename.chk filename.fchk
formchk -3 state10.chk state10.fchk

If we open filename.fchk (the file where the original TD-DFT calculation is located) with GaussView we can plot all orbitals involved in excited state number ten, those would be seven orbitals from 41 (HOMO-3) to 47 (LUMO+2) as shown in figure 1.

Figure 1. Canonical orbitals involved in the 10th excited state according to the TD-DFT calculation

If we now open state10.fchk we see that the numbers at the side of the orbitals are not their energy but their occupation number particular to this state of interest, so we only need to plot those with highest occupations, in our example those are orbitals 44 and 45 (HOMO and LUMO) which have occupations = 0.81186; you may include 43 and 46 (HOMO-1 and LUMO+1, respectively) for a much more complete description (occupations = 0.18223) but we’re still dealing with 4 orbitals instead of 7.

Figure 2. Natural Transition Orbitals for Phenylalanine. Orbital 44 (particle) and Orbital 45 (hole) exhibit the largest occupations for Excited State No. 10

The NTO transition 44 -> 45 is far easier to conceptualize than all the 10 combinations given in the canonical basis from the direct TD-DFT calculation. TD-DFT provides us with the correct transitions, NTOs just paint us a picture more readily available to the chemist mindset.

NOTE: for G09 revC and above, the %OldChk option is available, I haven’t personally tried it but using it to specify where the excitations are located and then write the NTOs of interest into a new chk file in the following way, thus eliminating the need of copying the original chk file for each state:

%OldChk=filename.chk
%chk=stateN.chk

NTOs are based on the Natural Hybrid orbitals vision by Löwdin and others, and it is said to be so straightforward that it has been re-discovered from time to time. Be that as it may, the NTO visualization provides a much clearer vision of the excitations occurring during a TD calculation.

Thanks for reading, stay home and stay safe during these harsh days everyone. Please share, rate and comment this and other posts.

A new paper on the Weak Link Approach


Chemically actuating a molecule is a very cool thing to do and the Weak Link Approach (WLA) allows us to do precisely that through the reversible coordination of one or various organometallic centers to a longer ligand that opens or closes a macrocyclic cavity. All this leads to an allosteric effect so important in biological instances available in inorganic molecules. Once again, the Mirkin group at Nortwestern University in Evanston, Illinois, has given me the opportunity to contribute with the calculations to the energetic properties of these actuators as well as their electronic properties for their use as molecular scavengers or selective capsules for various purposes such as drug delivery agents.

As in the previous WLA work (full paper), the NBODel procedure was used at the B97D/LANL2DZ level of theory, only this time the macrocycle consisted of two organometallic centers and for the first time the asymmetric opening of the cavity was achieved, as observed by NMR. With the given fragments, all possibilities shown in scheme 1 were obtained. The calculated bond energies for the Pt – S bonds are around 60 – 70 kcal/mol whereas for the Pt – Cl bonds the values are closer to 90 kcal/mol. This allows for a selective opening of the cavity which can then be closed by removing the chlorine atoms with the help of silver salts.

wla

For the case of complex mixture 4a, 4b, and 4c, the thermochemistry calculations show they are all basically isoenergetic with differences in the thousandths of kcal/mol. The possibilities for the groups in the weakly bonded ligands are enormous; currently, there is work being done about substituting those phenyl rings for calix[4]arenes in order to have a macrucyclic capsule made by macrocylic capusules.

Thanks to Andrea D’Aquino for taking me into her project, for all the stimulating discussions and her great ideas for expanding WLA into new avenues; I’m sure she’ll succeed in surprising us with more possibilities for these allosteric macrocycles.

The full paper is published in Inorganic Chemistry from the ACS (DOI: 10.1021/acs.inorgchem.7b02745). Thanks for reading and -if you made it this far into the post- happy new year!

Python scripts for calculating Fukui Indexes


One of the most popular posts in this blog has to do with calculating Fukui indexes, however, when dealing with a large number of molecules, our described methodology can become cumbersome since it requires to manually extract the population analysis from two or three different output files and then performing the arithmetic on them separately with a spreadsheet or something.

Our new team member Ricardo Loaiza has written a python script that takes the three aforementioned files and yields a .csv file with the calculated Fukui indexes, and it even points out which of the atoms exhibit the largest values so if you have a large molecule you don’t have to manually check for them. We have also a batch version which takes all the files in any given directory and performs the Fukui calculations for each, provided it can find file triads with the naming requirements described below.

Output files must be named filename.log (the N electrons reference state), filename_plus.log (the state with N+1 electrons) and filename_minus.log (the N-1 electrons state). Another restriction is that so far these scripts only work with NBO population analysis as provided by the NBO3.1 program available in the various versions of Gaussian. I imagine the listing is similar in NBO5.x and NBO6.x and so it should work if you do the population analysis with them.

The syntax for the single molecule version is:

python fukui.py filename.log filename_minus.log filename_plus.log

For the batch version is:

./fukuiPorLote.sh

(Por Lote means In Batch in Spanish.)

These scripts are available via GitHub. We hope you find them useful, and you do please let us know whether here at the comments section or at our GitHub site.

Physical Chemistry Meeting and CCIQS Symposium


Materials Research Institute at Morelia Michoacan (southern Mexico)

Materials Research Institute at Morelia Michoacan (southern Mexico)

So many events going on and so little time to blog about.
Two weeks ago, four members of this group traveled to Morelia in southern Mexico to present their research at the XIII Mexican Physical Chemistry Meeting. The next week after that, they all brought their posters back to Toluca for the internal symposium at CCIQS, where a masters student, María Eugenia, gave a small talk about her research project concerning photosynthesis in bacteria. Below, a short description of their projects is presented in order of seniority.

María Eugenia “Maru” Sandoval
Maru is working in photosynthesis of green sulfur bacteria. Her research deals with the excited states calculations at the Time Dependent DFT level for describing the first stages of photon interaction in antennae complexes of the photosystem II, namely the Fenna-Matthews-Olsen (FMO) complex, which was selected due to its relative structural simplicity over that of more evolved organisms. Maru also gave a talk at the internal Symposium back in Toluca the very next week where she got a positive feedback which will be used in her project.

Howard Diaz
One of the many strategies out there for treatment of HIV-1 infections is to block those proteins used to anchor the virus to a healthy cell. Sort of getting the virus’ hands busy so they can’t attach to a host. 60+ new compounds derived from thiourea were screened and assessed in their interactions with protein GP120, the protein to which the attachment is made, through docking and DFT calculations. Lead compounds are reported. It must be stressed that Howard got an award at CCIQS for having one of the best posters out of 70 in the entire symposium. Kudos and thanks to you, Howard! We now have some MD simulations in order.

Howard Diaz -  GP120 blockers for HIV-1

Howard Diaz – GP120 blockers for HIV-1

Guillermo “Memo” Caballero
His project has some nice philosophical implications if you ask me. Memo started as an experimental chemist and when he ran into a wall trying to obtain a pyridine from the non-aromatic analogue (glutarimide), he came to our group to run some calculations and find out how to force the aromatization process, or at least rationalize if it could be performed at all. Two mechanisms were proposed and now we know that even when the reaction should be quite exothermic, the reaction barriers are too high to be overcome by conventional methods. We now need to find a way to decrease those barriers (cue transition metal simulations). So in a way we are dealing here with the mechanism of a reaction that never happens (at least in an intramolecular way), leading to a reverse reductio ad absurdum reasoning – we assumed the reaction(s) did happen and we found out why is it impossible for them to happen.

No pic. available as of yet

Luis Enrique “Kike” Aguilar
Luis continues to work with calix(n)arenes, this lab’s first love, in drug delivery systems. He is working with two drugs at once: Bosutinib and Sorafenib, second generation drugs for the treatment of Chronic Myeloid Leukemia in cases were resistance to Imatinib has been developed. One of his main goals is to find a calixarene system which is able to discriminate between Bosutinib and pseudo-bosutinib, a commercial isomer which has incorrectly been available for a few years now.

"Kike" Aguilar -  Calix[n]arenes as drug carriers for Bosutinib and Sorafenib

“Kike” Aguilar – Calix[n]arenes as drug carriers for Bosutinib and Sorafenib

Thanks a lot for your efforts; they are paying good results to your ca

The boys from the lab

The boys from the lab

reers and the advancement of our research group. Now back to work, guys!

 

New paper in JACS


Well, I only contributed with the theoretical section by doing electronic structure calculations, so it isn’t really a paper we can ascribe to this particular lab, however it is really nice to see my name in JACS along such a prominent researcher as Prof. Chad Mirkin from Northwestern University, in a work closely related to my area of research interest as macrocyclic recognition agents.

In this manuscript, a calix[4]arene is allosterically opened and closed reversibly by coordinating different kinds of ligands to a platinum center linked to the macrocycle. (This approach has been referred to as the weak link approach.) I recently visited Northwestern and had a great time with José Mendez-Arroyo, the first author, who showed me around and opened the possibility for further work between our research groups.

(Ligands: Green = Chloride; Blue = Cyanide)

Closed, semi-open and fully open conformations; selectivity is modulated through cavity size. (Ligands: Green = Chloride; Blue = Cyanide)

Here at UNAM we calculated the interaction energies for the two guests that were successfully inserted into the cavity: N-methyl-pyridinium (Eint = 57.4 kcal/mol) and Pyridine-N-oxide (Eint = +200.0 kcal/mol). Below you can see the electrostatic potential mapped onto the electron density isosurface for one of the adducts. Relative orientation of the hosts within the cavity follows the expected (anti-) alignment of mutual dipole moments. At this level of theory, we could easily be inclined to assert that the most stable interaction is indeed the one from the semi-open compound and that this in turn is due to the fact that host and guest are packed closer together but there is also an orbital issue: Pyridine Oxide is a better electron acceptor than N-Me-pyridinium and when we take a closer look to the (Natural Bonding) orbitals interacting it becomes evident that a closer location does not necessarily yields a stronger interaction when the electron accepting power of the ligand is weaker (which is, in my opinion, both logic and at the same time a bit counterintuitive, yet fascinating, nonetheless).

Electrostatic potential mapped onto the electron density surface of one of the aducts under study

Electrostatic potential mapped onto the electron density surface of one of the adducts under study

All calculations were performed at the B97D/LANL2DZ level of theory with the use of Gaussian09 and NBO3.1 as provided within the former. Computing time at UNAM’s supercomputer known as ‘Miztli‘ is fully acknowledged.

The full citation follows:

A Multi-State, Allosterically-Regulated Molecular Receptor With Switchable Selectivity
Jose Mendez-Arroyo Joaquín Barroso-Flores §,Alejo M. Lifschitz Amy A. Sarjeant Charlotte L. Stern , and Chad A. Mirkin *

J. Am. Chem. Soc., Article ASAP
DOI: 10.1021/ja503506a
Publication Date (Web): July 9, 2014

 Thanks to José Mendez-Arroyo for contacting me and giving me the opportunity to collaborate with his research; I’m sure this is the first of many joint projects that will mutually benefit our groups. 

 

New paper in Computational and Theoretical Chemistry


I always get very happy to have a new paper out there! I find it exciting but most of all liberating since it makes you feel like your work is going somewhere but most of all that it is making its way ‘out there’; there is a strong feeling of validation, I guess.

Two very different families of calix[n]arenes (Fig 1) were tested as drug carriers for a very small molecule with a huge potential as a chemotherapeutic agent against Leukemia, namely, 3-phenyl-1H-[1]benzofuro[3,2-c]pyrazole a.k.a. GTP which has proven to be an effective in vitro Tyrosine Kinase III inhibitor. Having such a low molecular weight it is expected to have a very high excretion rate therefore the use of a carrier could increase its retention time and hence its activity. This time we considered n = 4, 5, 6 and 8 for the size of the cavities and R = -SO3H and -OEt as functional groups on the upper rim as to evaluate only a polar coordinating group and a non-polar non-coordinating one since GTP offers two H-bond acceptor sites and one H-bond donor one along the π electron density that could form π – π stacking interactions between the aromatic groups on GTP and the walls of the calixarene.

Fig 1. Calixarenes under study and their complexes with GTP

Fig 1. Calixarenes under study and their complexes with GTP

Once again calculations were carried out at the B97D/6-31G(d,p) level of theory along with molecular dynamics simulations for over 100 ns of production runs. NBO Deletion interaction energies were computed in order to discern which hosts formed the most stable complexes.

NBO Del interaction energies B97D/6-31G(d,p)

NBO Del interaction energies B97D/6-31G(d,p)

You may find a link to the ScienceDirect website for downloading the paper from this link. Last, but certainly not least, I’d like to thank all coauthors for their contributions and patience in getting this study published: Dr. Rodrigo Galindo-Murillo; Alberto Olmedo-Romero; Eduardo Cruz-Flores; Dr. Petronela M. Petrar and Prof. Dr. Kunsági-Máté Sándor. Thanks a lot for everything!

fig8

Donor and acceptor H-bond sites increases the probability of keeping the drug in place for a higher retention rate

Donor and acceptor H-bond sites increases the probability of keeping the drug in place for a higher retention rate

New paper in Journal of Chemical Theory and Computation


Happy new year to all my readers!

Having a new paper published is always a matter of happiness for this computational chemist but this time I’m excedingly excited about anouncing the publishing of a paper in the Journal of Chemical Theory and Computation, which is my highest ranked publication so far! It also establishes the consolidation of our research group at CCIQS as a solid and competitive group within the field of theoretical and computational chemistry. The title of our paper is “In Silico design of monomolecular drug carriers for the tyrosine kinase inhibitor drug Imatinib based on calix- and thiacalix[n]arene host molecules. A DFT and Molecular Dynamics study“.

In this article we aimed towards finding a suitable (thia-) calix[n]arene based drug delivery agent for the drug Imatinib (Gleevec by Novartis), which is a broadly used powerful Tyrosine Kinase III inhibitor used in the treatment of Chronic Myeloid Leukaemia and, to a lesser extent, Gastrointestinal Stromal Tumors; although Imatinib (IMB) exhibits a bioavailability close to 90% most of it is excreted, becomes bound to serum proteins or gets accumulated in other tissues such as the heart causing several undesired side effects which ultimately limit its use. By using a molecular capsule we can increase the molecular weight of the drug thus increasing its retention, and at the same time we can prevent Imatinib to bind, in its active form, to undesired proteins.

We suggested 36 different calix and thia-calix[n]arenes (CX) as possible candidates; IMB-CX complexes were manually docked and then optimized at the B97D/6-31G(d,p) level of theory; Stephan Grimme’s B97D functional was selected for its inclusion of dispersion terms, so important in describing π-π interactions. Intermolecular interaction energies were calculated under the Natural Bond Order approximation; a stable complex was needed but a too stable complex would never deliver its drug payload! This brings us to the next part of the study. A monomolecular drug delivery agent must be able to form a stable complex with the drug but it must also be able to release it. Molecular Dynamics simulations (+100 ns) and umbrella sampling methods were used to analyse the release of the drug into the aqueous media.

Optimized geometries for all complexes under study (B97D/6-31G*)

Optimized geometries for the 20 most stable complexes under study (B97D/6-31G*)

Potential Mean Force profiles for the four most stable complexes for position N1 and  N2 from the QM simulations are shown below (Red, complexes in the N1 position, blue, N2 position). These plots, derived from the MD simulations  give us an idea of the final destination of the drug respect of the calixarene carrier. In the next image, the three preferred structures (rotaxane-like; inside; released) for the final outcome of the delivery process are shown. The stability of the complexes was also assessed by calculating the values of ΔG binding through the use of the Poisson equations.

PMF for the most stable compounds

PMF for the most stable compounds

General MD simulation final structures

General MD simulation final structures

Thanks to my co-authors Maria Eugenia Sandoval-Salinas and Dr. Rodrigo Galindo-Murillo for their enormous contributions to this work; without their hard work and commitment to the project this paper wouldn’t have been possible.

%d bloggers like this: