home

Archive for the 'qc/mm/md progs' Category

The Binding Of Vitamin B12 To Transcobalamin(II); Structural Considerations For Bioconjugate Design – A Molecular Dynamics Study

Thursday, July 22nd, 2010

In press, in the journal Molecular Biosystems. A first official foray into molecular dynamics-only (MD-only) computational work and I am pleased to report that the computational results not only make sense with respect to the experimental results, they also indicate a possible new way to use vitamin B12 for the oral delivery of bio-active molecules more complicated than the binary bioconjugates considered to date.

The Interesting Result

The conclusion from the previous study was that the insulin B Chain (figure below) acts as a tether to separate the structured region of insulin (the region with the largest inflexible steric bulk, see below) from the region of the transcobalamin II (TCII) that bind vitamin B12. It was then determined that the approach employed for the B12-insulin bioconjugate, simply linking one biomolecule onto another with known binding and transport properties (this is a common theme in all bioconjugate design), worked because the last 10 residues in the insulin B Chain (B22 to B30) are flexible in solution (they, in fact, cover the insulin binding region in the crystal form, then uncover this region in the biologically active form).

As a general procedure for B12 bioconjugate design, one of the key requirements for a functional product is a tether length that provides sufficient separation between B12 and any molecular structure large enough to affect B12 binding within its transport proteins (makes sense, as a tethered structure that does not enable B12 binding in its transport proteins will find the B12 bioconjugate delivered to the gut where acids and digestive enzymes will hide the failed binding). This leads to the question, “How long must a tether be to meet this rather general criterion?” This is, partly, the correct question, as the retention of B12 binding within its transport proteins is a function of both proper tether length and [transport protein]-["other molecule"] interaction (in this first case, “other molecule” = insulin).

Saving the exhaustive analysis for the paper, this new study used this flexible region of human insulin (that is, B22 to B30, with the B12 linkage occurring on the B29 lysine side chain) as a proxy for any arbitrary tether, then used MD simulations to consider how the flexibility of this tether might lead to changes in B12 binding within its TCII pocket (the transport protein for which we have the best crystal structure). The result of these simulations was the identification of the side chain of lysine itself being just long enough to separate the B Chain tether region from the TCII protein surface. This does not mean that lysine will always serve as a perfect linkage. This means that, if the tether structure is effectively non-interacting with TCII (so not sterically demanding by itself), the lysine side chain is long enough to span the solvent-accessible hole produced by the encapsulation of B12 in (in this case) TCII.

The result is a design constraint when using lysine that is quite fortuitous! If the target peptide (insulin or whatnot) has a surface-accessible lysine side chain within a region that is flexible in solution, some simple amide chemistry may produce a viable B12 bioconjugate for delivering that peptide orally (thereby avoiding complete peptide degradation in the G.I. tract).

The More Interesting Result

Buried deep within the bottom of the Discussion section. If you watch the dynamics simulation of the TCII-[B12-tether] complex (shown below for a 300 K 50 ns simulation with 1.5 fs time steps in 14,000 waters (not shown)), you see that the binding of B12 within TCII and the geometry of the encapsulation complex are strongly linked. That is, TCII (and, presumably, its cohorts in the B12 transport pathway) can be thought of as two quite rigid fragments (Red and Blue in the animation) connected by a long tether (Green) that are separated in solution but brought into contact by the binding of vitamin B12 (Gold). The B12 is a glue that holds the fragments together, and a simple tabulation of hydrogen-bonding interactions in the crystal structure reveal that the B12 has more interactions to the A and B fragments of TCII individually than A and B have with each other (which is to say, the B12-A Segment interaction and B12-B Segment interaction are stronger than the A-B Segment interaction). From a biological perspective, this should make perfect sense. B12 is a large, extremely important biomolecule that, since we do not make it ourselves, is to be captured and transported as effectively as possible. The best way to bind this molecule is not to wait for it to burrow into a binding pocket, but rather to encapsulate it in a “clam shell” maneuver that provides “maximum embedding.” The tether between the A and B Segments technically would not have to be present if the A and B fragments were present in large quantities (although, as you might expect, the A-B tether does considerably reduce the time to complete encapsulation by forcing these fragments within close proximity).

According to the crystal structure, the B12 is entirely embedded within TCII, with only the solvent-accessible hole at the 5′-ribose position readily accessible for bioconjugate formation. If the overall structure were as rigid as a crystal structure might lead one to believe, functionalization at the cobalt position in the corrin ring would be out of the question.

As I just stated that such a binding mode would otherwise be unlikely, you can guess that there are B12 bioconjugates linked at the cobalt ring that are bio-active.

If you watch the dynamics simulation of the TCII-[B12-tether] complex, you see that the clam shell binding mode of TCII is one with a “loose hinge.” This loose hinge is really a result of the flexibility of the two protein fragments (typical protein motion) and flexibility in the short propionamide side chains of vitamin B12 that provide a bit of “spring” in the complete complex. In effect, the flexibility within the structure provides a means for cobalt to be coordinated to something without loss of B12 binding provided that the tether linking the cobalt and the “other” molecule is small enough that it does not require a large change in the A-B binding arrangement (that is, does not affect B12-A and B12-B binding).

And Then There Were Three…

The expectation/prediction/untested hypothesis is that vitamin B12 may be able to happily accommodate two additional molecules at the 5’-ribose and cobalt positions (properly designed) that then provide for the transport of two molecules and/or the delivery of three molecules (one being vitamin B12). This opens the door to a wealth of possibilities, from trinary delivery to combined drug delivery + radiopharma characterization. This is the possibility I’m most interested in pursuing in the next rounds of calculations, with the theory (presumably) providing a very good initial guess about the ideal tether designs to use with B12 for enabling delivery and bio-activity.

And Now For The Hard Work

Stepping back from the theoretical analysis for a moment, the most difficult obstacles to overcome in this study were the generation AND incorporation of force field parameters for vitamin B12 and a B12-Lysine mini-bioconjugate into GROMACS, a problem that I’ve addressed only in passing in several previous posts. What I won’t do in this post is explain the procedure (a single blog post will not do the procedure justice given the complexity of force field parameter generation). What I will do is provide the files for the topology for these systems and a short list of the modifications one needs to make in order to get these systems working. For additional reference, the same topology files are provided in the Supplemental Material for the paper (so, if you find yourself using these, obviously cite the paper and not my humble blog).

Files And Contents:

These are not files to be placed in a single directory, but are segments of file that are going to be placed directly into pre-existing topology files. This is not the best way to do it but is the procedure I began with and will not be changing without finding a very simple tutorial on how-to (which, if you have, I’d be happy to read).

The contents of the topology file (which I assume for you will be ffG53a6 but should work generally) are provided below:

ffG53a6_B12_BCN_LYB_LCB_topology.txt

The topology specifications for vitamin B12 (nothing bound to the cobalt in the corrin ring), cyanocobalamin (CN-B12, with a cyanide bound to the cobalt), B12 with a lysine residue attached to the 5’-ribose hydroxyl position (the tether linkage for the GROMACS prep programs), and CN-B12 with a lysine residue attached to the 5’-ribose hydroxyl position.

I am assuming that you’re using the ffG53a6 force field, meaning you add the topology sets to the bottom of the ffG53a6.rtp file.

GROMACS Modifications:

GROMACS force field and topology files must be modified slightly in order to read the topologies generated above and, depending on where you got the B12 structure, add/correct the hydrogen atoms in the B12 molecule.

In a typical UNIX/Linux installation (which I have provided compilation instructions for in a previous post), the files to be modified can be found in /usr/local/gromacs. And, if you’re using Ubuntu like I am, you’ll need to “sudo” these modifications.

1. aminoacids.dat

If you open this file, you see a list of three- and four-letter codes in the format:

50
ABU
ACE
...
VAL
PGLU

The “50” refers to the number of codes. As we’re going to be adding the codes B12, BCN, LYB, and LCB into GROMACS, we first change 50 to 54, then just list the four codes at the bottom of the file:

54
ABU
ACE
...
VAL
PGLU
B12
BCN
LYB
LCB

You’ll note that B12 and BCN aren’t like the others, LYB is not LYS, and LCB is also nowhere to be seen. The codes in this file are STANDARD and make sure you don’t inadvertently name your inserted structure one of the structures in the list.

2. ffG53a6.hdb

I specifically used the ffG53a6 force field for the TCII-B12 work, meaning I only made modifications to these force field files. The ffG53a6.hdb file is responsible for adding/correcting hydrogen atoms in your structure (just because the crystallographers do not see them does not mean they aren’t there) and contains hydrogen-beautification information for all of the three/four-letter codes recognized in aminoacids.dat. The content below is the hydrogen-correcting data for the B12, BCN, LYB, and LCB structures. Simply paste this into the bottom of the ffG53a6.hdb file.

B12     19
1    2    HAO    N62    C61    O63
1    2    HAN    N62    C61    C60
1    2    HAM    N52    C50    O51
1    2    HAL    N52    C50    C49
1    2    HAK    N45    C43    O44
1    2    HAJ    N45    C43    C42
1    2    HAI    N40    C38    O39
1    2    HAH    N40    C38    C37
1    2    HAE    N29    C27    O28
1    2    HAD    N29    C27    C26
1    2    HAG    N33    C32    O34
1    2    HAF    N33    C32    C31
1    2    HAA    O7R    C2R    C1R
1    2    HAB    O8R    C5R    C4R
1    2    HAC    N59    C57    O58
1    1    H2B    C2B    N1B    N3B
1    1    H4B    C4B    C5B    C9B
1    1    H7B    C7B    C8B    C6B
1    1    H10    C10    C9     C11
LYB     20
1    1    H      N      -C     CA
1    4    HZ1    NZ     CE     CD
1    2    HAO    N62    C61    O63
1    2    HAN    N62    C61    C60
1    2    HAM    N52    C50    O51
1    2    HAL    N52    C50    C49
1    2    HAK    N45    C43    O44
1    2    HAJ    N45    C43    C42
1    2    HAI    N40    C38    O39
1    2    HAH    N40    C38    C37
1    2    HAE    N29    C27    O28
1    2    HAD    N29    C27    C26
1    2    HAG    N33    C32    O34
1    2    HAF    N33    C32    C31
1    2    HAA    O7R    C2R    C1R
1    2    HAC    N59    C57    O58
1    1    H2B    C2B    N1B    N3B
1    1    H4B    C4B    C5B    C9B
1    1    H7B    C7B    C8B    C6B
1    1    H10    C10    C9     C11
BCN     19
1    2    HAO    N62    C61    O63
1    2    HAN    N62    C61    C60
1    2    HAM    N52    C50    O51
1    2    HAL    N52    C50    C49
1    2    HAK    N45    C43    O44
1    2    HAJ    N45    C43    C42
1    2    HAI    N40    C38    O39
1    2    HAH    N40    C38    C37
1    2    HAE    N29    C27    O28
1    2    HAD    N29    C27    C26
1    2    HAG    N33    C32    O34
1    2    HAF    N33    C32    C31
1    2    HAA    O7R    C2R    C1R
1    2    HAB    O8R    C5R    C4R
1    2    HAC    N59    C57    O58
1    1    H2B    C2B    N1B    N3B
1    1    H4B    C4B    C5B    C9B
1    1    H7B    C7B    C8B    C6B
1    1    H10    C10    C9     C11
LCB     20
1    1    H      N      -C     CA
1    4    HZ1    NZ     CE     CD
1    2    HAO    N62    C61    O63
1    2    HAN    N62    C61    C60
1    2    HAM    N52    C50    O51
1    2    HAL    N52    C50    C49
1    2    HAK    N45    C43    O44
1    2    HAJ    N45    C43    C42
1    2    HAI    N40    C38    O39
1    2    HAH    N40    C38    C37
1    2    HAE    N29    C27    O28
1    2    HAD    N29    C27    C26
1    2    HAG    N33    C32    O34
1    2    HAF    N33    C32    C31
1    2    HAA    O7R    C2R    C1R
1    2    HAC    N59    C57    O58
1    1    H2B    C2B    N1B    N3B
1    1    H4B    C4B    C5B    C9B
1    1    H7B    C7B    C8B    C6B
1    1    H10    C10    C9     C11

As brief explanation, the three-letter code is followed by the number of Hydrogen atoms that are to be added. Each line can be read:

First Column – The number of hydrogen atoms added (so all of these entries on the far left mean “add ONE hydrogen”)

Second Column – The manner by which the hydrogen atom is to be added (this is listed in section 5.5 of the GROMACS 3.3 Manual (page 93))

Third Column – The name of the Hydrogen atom to be added

Fourth Column – The atom to which the H is going to be directly linked in the topology file

Fifth – Seventh Columns
– atoms that define how the Hydrogen is added with respect to (1) the code in Column 2 and (2) the atom to which the Hydrogen is added.

3. ffG53a6bon.itp

There are a few subtle tweaks to the force constants for a few bonds that I perform here right within the file and that proper MD people likely would scream at. I note that, when you do this, you are making changes to numbers that will affect the results if you somehow start doing heme MD simulations.

Change the gb_NN values to those provided below.

#define gb_34        0.198  0.6400e+06
; NR  -   FE    120
#define gb_4         0.1142  3.7000e+07
; C - O (CO in heme)  2220
#define gb_14       0.1340  1.1000e+07
; C  -  NR (heme)       1000
#define gb_30       0.1880  2.7200e+06
; FE  -  C (Heme)

You will note that I have not done anything to make cobalt appear in the topology or force field files. For the sake of running a simulation, Fe and Co are close enough that simply replacing CO for FE in the PDB file is sufficient. You can do the completely proper job of adding cobalt to the force field to get the mass right.

And that is the bare basics for getting a run to happen. A proper tutorial on how to generate force field parameters and topologies may be forthcoming, depending largely on interest and my ability to find time to do it.

Article citation: Damian G. Allis, Mol. BioSyst., 2010, DOI: 10.1039/c003476b

Damian G. Allis1, Timothy J. Fairchild2 and Robert P. Doyle1

1. Department of Chemistry, Syracuse University, Syracuse, NY 13244, USA
2. School of Chiropractic and Sports Science, Murdoch University, Murdoch, WA 6150, Australia

As part of ongoing research into the use of vitamin B12 (B12; cobalamin; Cbl)-based bioconjugate approaches for the oral delivery of peptides/proteins, a molecular dynamics (MD) study of the binding of a cyanocobalamin–insulin (CN–Cbl–insulin) conjugate to human transcobalamin(II) (TCII) was recently reported that provides a qualitative picture of how the human insulin protein in its open T-state geometry affects CN–Cbl binding to TCII. This initial analysis revealed that the B22–B30 segment of the insulin B-chain acts as a long tether that connects the larger combined insulin A/B region to CN–Cbl when this conjugation is performed at the CN–Cbl ribose 5-hydroxy position. The experimental support for this model of the binding interaction is provided by the consequences of the successful delivery of the CN–Cbl–insulin conjugate in the production of significantly decreased blood glucose levels in diabetic STZ-rat models. In efforts to provide a more detailed description of the (CN–Cbl)–TCII complex for modeling Cbl-based bioconjugate designs, the (CN–Cbl)–TCII system and a CN–Cbl conjugate incorporating a flexible tether composed of only the B22–B30 segment of human insulin have been examined by MD simulations. The implications of these simulations are discussed in terms of successful conjugate positioning on Cbl, especially when such sites are not apparent from the diffraction studies alone, and the possibilities, as yet not reported, for dual-tethered Cbl bioconjugates for multi-component drug delivery applications.

Running (Only) A Single-Point Energy Calculation In Crystal06/09, Proper Input Format For Long-Range Dispersion Contributions In Crystal09, And Removing The MPICH2 Content From The Output File In Pcrystal

Saturday, May 29th, 2010

Now enjoying the benefits of dispersion-corrected solid-state density functional theory (and a proper MPICH2 implementation for infrared intensity calculations, although this now a problem for reasons to be addressed in an upcoming post) in Crystal09, three issues in recent calculations caused me to think hard enough about keyword formats and job runs that I have opted to post briefly about what to do in case google and bing are your preferred methods of manual searching.

1. How To Run Only A Single-Point Energy Calculation In Crystal06/Crystal09

This had never come up before and, by the time I needed to find an input file to see what do to, the first google search provided Civalleri’s Total Energy Calculation page that currently has broken links to .zip files. There is quite a bit about the different geometry optimization approaches in the manual, but a search for “single-point” provides no information about what to do for only single-point energy calculations.

The solution, it should be obvious after, is simply to not include the geometry optimization section in the input file. What would otherwise be the following (with arbitrary geometry optimization-like info between [COORDINATES] and [BASIS SETS]…

[COORDINATES]
OPTGEOM
TOLDEG
0.000005
TOLDEX
0.000020
END
END
[BASIS SETS]

becomes…

[COORDINATES]
[BASIS SETS]

One problem solved by simply not having any optimization parameters (again, makes sense and is now google-able).

2. Proper GRIMME Input Format For Long-Range Dispersion Contributions In Crystal09

This is another example where one’s first efforts in translating the manual into calculations may lead to considerable confusion until the proper format is finally identified (by which time you’ve run many pruned-down input tests).

GRIMME
1.05 20. 25.
1.05 20. 25. s6 (scaling factor) d (steepness) Rcut (cutoff radius)
5
1  0.14 1.001 Hydrogen Conventional Atomic number , C6 , Rvdw
6  1.75 1.452 Carbon Conventional Atomic number , C6 , Rvdw
7  1.23 1.397 Nitrogen Conventional Atomic number , C6 , Rvdw
8  0.70 1.342 Oxygen Conventional Atomic number , C6 , Rvdw
17 5.07 1.639 Chlorine Conventional Atomic number , C6 ,'Rvdw

I’m not even sure where the final ,’Rvdw comes from. Your .out file may terminate with the following error (or something similar)…

rank 7 in job 8  korterquad_51438   caused collective abort of all ranks
  exit status of rank 7: killed by signal 9

And the ERROR.peN file with any content will show the following, clearly pointing to a GRIMME-specific error…

 ERROR **** GRIMME_INPUT **** ELEMENT NOT DEFINED:           1

The problem is the additional content within the manual pages for the GRIMME keyword that require pruning (or, at least, some identifier to show what is and what is not needed). The proper GRIMME section above is properly provided in the INPUT file as…

GRIMME
1.05 20. 25.
5
1  0.14 1.001
6  1.75 1.452
7  1.23 1.397
8  0.70 1.342
17 5.07 1.639

Where (see page 88 of the Crystal09 manual)…

GRIMME <- keyword is called
1.05 20. 25. <- scaling factor, steepness, cutoff distance
5 <- number of elements in the list (not the total number of atoms)
1  0.14 1.001 <- atomic number, dispersion coefficient, van der Waals radius
...

When all is properly run, the bottom of your output file will look something like the following:

 CYC  43 ETOT(AU) -5.784662098123E+03 DETOT  1.18E-11 tst  8.17E-15 PX  6.73E-08

 == SCF ENDED - CONVERGENCE ON ENERGY      E(AU) -5.7846620981229E+03 CYCLES  43

 ENERGY EXPRESSION=HARTREE+FOCK EXCH*0.20000+(BECKE  EXCH)*0.80000+LYP    CORR

 TOTAL ENERGY(DFT)(AU)( 43) -5.7846620981229E+03 DE 1.2E-11 tester 8.2E-15
 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT EDFT        TELAPSE     4705.82 TCPU     4651.41

 *******************************************************************************

 GRIMME DISPERSION ENERGY CORRECTION

 SCALE FACTOR (S6):     1.0500

 GRIMME DISPERSION ENERGY (AU) -1.9723347118951E-01
 TOTAL ENERGY + DISP (AU) -5.7848593315941E+03

 *******************************************************************************

The Crystal09 manual refers you to Table 1 of the Stefan Grimme paper, “Semiempirical GGA-type density functional constructed with a long-range dispersion correction” (Journal of Computational Chemistry, Volume 27, Issue 15, Pages 1787 – 1799), which I’ve put together into the proper format below. Be sure to (1) delete the elements in parentheses ( -> get rid of the (H) <- ), (2) remove those atoms you do not need, (3) be sure to change the “number of elements” number for your structure, and (4) get and reference the Grimme paper so you have the proper C6 parameters and van der Waals radii accounted for (you’ll be the right nitwit if I mis-copied something and you ran with it (although I trust my input), and you should have the reference regardless).

( H)   1   0.14 1.001
(Li)   3   1.61 0.825
(Na)  11   5.71 1.144
( K)  19  10.80 1.485
(Rb)  37  24.67 1.628
(Be)   4   1.61 1.408
(Mg)  12   5.71 1.364
(Ca)  20  10.80 1.474
(Sr)  38  24.67 1.606
( B)   5   3.13 1.485
(Al)  13  10.79 1.639
(Ga)  31  16.99 1.650
(In)  49  37.32 1.672
( C)   6   1.75 1.452
(Si)  14   9.23 1.716
(Ge)  32  17.10 1.727
(Sn)  50  38.71 1.804
( N)   7   1.23 1.397
( P)  15   7.84 1.705
(As)  33  16.37 1.760
(Sb)  51  38.44 1.881
( O)   8   0.70 1.342
( S)  16   5.57 1.683
(Se)  34  12.64 1.771
(Te)  52  31.74 1.892
( F)   9   0.75 1.287
(Cl)  17   5.07 1.639
(Br)  35  12.47 1.749
( I)  53  31.50 1.892
(He)   2   0.08 1.012
(Ne)  10   0.63 1.243
(Ar)  18   4.61 1.595
(Kr)  36  12.01 1.727
(Xe)  54  29.99 1.881
Y-Cd      24.67 1.639
Sc-Zn     10.80 1.562

Note that the d-block is identical for each row (so no atom numbers provided).

3. Removing The MPICH2 Content From The Output File In Pcrystal(/09)

This final issue does not occur in Pcrystal(/06) but does in Pcrystal(/09), with the reason being (I assume) the new use of MPICH2 in Pcrystal(/09) instead of MPICH in Pcrystal(/06).  The problem comes from running the following set of commands at the terminal window in MPICH2:

mpiexec -machinefile machine -np N /path/to/Pcrystal &>FILENAME.out &

Embedded within the FILENAME.out file will be all flavors of MPI-specific output, perhaps such as the following (in this case errors, but it happens in proper output as well):

application called MPI_Abort(MPI_COMM_WORLD, 1) - process 4
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 7
rank 7 in job 9  korterquad_51438   caused collective abort of all ranks
 exit status of rank 7: return code 1
rank 4 in job 9  korterquad_51438   caused collective abort of all ranks
 exit status of rank 4: killed by signal 9

or…

mpiexec_machine (handle_stdin_input 1089): stdin problem; if pgm is run in background...
mpiexec_machine (handle_stdin_input 1090):     e.g.: mpiexec -n 4 a.out < /dev/null &

The solution is to break up the mpiexec output from the Pcrystal output, performed by directing the mpiexec-specific content to, in this case, /dev/null (because it is not necessary except for diagnostic purposes).

mpiexec -machinefile machine -np N /path/to/Pcrystal < /dev/null &>FILENAME.out &

Which removes all traces of mpi-specific output from FILENAME.out.

B12-Insulin Bioconjugate/Transcobalamin(II)/Insulin Receptor Cover Image For The April Issue Of Clinical Chemistry

Sunday, April 11th, 2010

A brief post about some free research press (and the new addition to the Cover Gallery). Having already been featured on the cover of the ChemMedChem March 2009 issue (see the New B12-Insulin-TCII-Insulin Receptor Cover Image For This Month’s ChemMedChem (March 2009) post) , the side-on view of the B12-Insulin/TCII/Insulin Receptor structure was chosen for this month’s cover of Clinical Chemistry. While the originating article itself is not included in the issue (I should have recommended citing the ChemMedChem article in the image caption), several diabetes-related articles are featured in this month’s issue.

ON THE COVER: Scientists are investigating ways to develop effective oral insulin therapies. One such model is a vitamin B12–insulin conjugate bound to transcobalamin II and is shown here docked in the insulin receptor. The discovery of easier ways to deliver insulin into the blood stream would improve the lives of the millions of individuals living with diabetes. This month’s issue of Clinical Chemistry contains 4 articles related to diabetes. The first 2 articles provide readers with a point/counterpoint discussion of the value of reporting estimated glucose along with Hb A1c. Next is an article on the association of apolipoprotein B with incident type 2 diabetes. Lastly, the development of the first radioimmunoassay for insulin led to a Nobel Prize and is chronicled in this month’s Citation Classic feature. (See pages 545, 547, 666, and 671.) Image reproduced with permission from Damian G. Allis and Robert P. Doyle, Department of Chemistry, Syracuse University.

As a brief explanation of the image, this “scene” is meant to show (without proper molecular dynamics simulations to show how well it would work) that the Transcobalamin(II) transport/protection protein for cobalamin/cyanocobalamin (vitamin B12) and the B12-insulin bioconjugate discussed in the ChemMedChem article is small enough to fit within the Insulin Receptor protein such that insulin may still be able to bind to its receptor. This is the final piece of the puzzle in the proposed mechanism (and experimentally demonstrated event) by which the B12-insulin bioconjugate retains all of the benefits of free B12 (transport from the digestive system to the bloodstream) and insulin (proper receptor binding and the subsequent induction of cellular glucose uptake).

The figure caption and April 2010 Table of Contents can be found in PDF format at the Clinical Chemistry website (with a local copy of the PDF also available HERE.

www.somewhereville.com/?page_id=985
www3.interscience.wiley.com/journal/122250806/issue
www.somewhereville.com/?p=511
www.clinchem.org
en.wikipedia.org/wiki/Diabetes
en.wikipedia.org/wiki/Molecular_dynamics
en.wikipedia.org/wiki/Cyanocobalamin
en.wikipedia.org/wiki/Vitamin_B12
en.wikipedia.org/wiki/Bioconjugate
en.wikipedia.org/wiki/Insulin_receptor
www.clinchem.org/content/vol56/issue4/

Terahertz Spectroscopic Investigation Of S-(+)-Ketamine Hydrochloride And Vibrational Assignment By Density Functional Theory, “Function Follows Functional Follows Formalism”

Sunday, February 21st, 2010

Accepted in the Journal of Physical Chemistry A, with my fingers crossed for pulling off the rare double-header in an upcoming print edition of the journal (having missed it by three intermediate articles with the Cs2B12H12 and HMX papers back in 2006 (you’d keep track, too). A fortuitous overlap of scheduled defense dates between P. Hakey, Ph.D. and M. Hudson, A.B.D.). A brief summary of interesting points from this study is provided below, including what I think is a useful point about how to most easily interpret AND represent solid-state vibrational spectra for publications.

1. AS USUAL, YOU CANNOT USE GAS-PHASE CALCULATIONS TO ASSIGN SOLID-STATE TERAHERTZ SPECTRA. It will take a phenomenal piece of data and one helluvan interpretation to convince me otherwise. As a more subtle point (for those attempting an even worse job of vibrational mode assignment), if the molecule exists in its protonated form in the solid-state, do not use the neutral form for your gas-phase calculation (this is a point that came up as part of an MDMA re-assignment published (and posted here) previously).

2. It is very difficult to find what I would consider to be “complete data sets” for molecules and solids being studied by spectroscopic and computational methods. For many molecular solids, the influences of thermal motion are not important to providing a proper vibrational analysis by solid-state density functional theory methods. Heating a crystal may make spectral lines broader, but phase changes and unusual spectral features do not often result when heating a sample from cryogenic (say, liquid nitrogen) to room temperature. Yes, there are thousands of cases where this is not true, but several fold more cases where it is. We are fortunate to live in a temperature regime where characterization is reasonably straightforward and yet we can modify a system to observe its subtle changes under standard laboratory conditions. The THz spectrum of S-(+)-Ketamine Hydrochloride gets a bit cleaner upon cooling, which makes the assignment easier. As the ultimate goal is to be able to characterize these systems in a person’s pocket instead of their liquid nitrogen thermos, the limited observed change to the spectrum upon cooling is important to note.

3. Crystal06 vs. DMol3 – This paper contains what is hoped to be a level, pragmatic discussion about the strengths and weaknesses of computational tools available to terahertz spectroscopists for use in their efforts to assign spectra. This type of discussion is, as a computational chemist using tools and not developing tools, a touchy subject to present on not because of the finger-pointing of limitations with software, but because the Crystal06 team and Accelrys (through Delley’s initial DMol3 code) clearly are doing things that the vast majority of their users (myself included) could in no way do by themselves. The analysis for the theory-minded terahertz spectroscopist is presented comparing two metrics – speed and functionality (specifically, infra-red intensity prediction). What is observed as the baseline is that both DMol3 and Crystal06 make available density functionals and basis sets that, when used at high levels of theory and rigorous convergence criteria, produce simulated terahertz spectra with vibrational mode energies that are in good (if not very good) agreement with each other. For the terahertz spectroscopist, Crystal06 provides as output (although this is system size- and basis set size-dependent) rigorous infrared intensity predictions for vibrational modes, inseparable from mode energy as “the most important” pieces of information for mode assignments. While DMol3 does not produce infrared intensities (the many previous terahertz papers I’ve worked on employed difference-dipole calculations that are, at best, a guesstimate), DMol3 produces very good mode energy predictions in 1/6th to (I’ve seen it happen) 1/10th the time of a comparable Crystal06 calculation. This is the reason DMol3 has been the go-to program for all of the neutron scattering spectroscopy papers cited on this blog (where intensity is determined by normal mode eigenvectors, which are provided by both (and any self-respecting quantum chemical code) programs).

Now, it should be noted that this difference in functionality has NOTHING to do with formalism. Both codes are excellent for what they are intended to do. To the general assignment-minded spectroscopist (the target audience of the Discussion in the paper), any major problem with Crystal06 likely originates with the time to run calculations (and, quite frankly, the time it takes to run a calculation is the worst possible reason for not running a calculation if you need that data. Don’t blame the theory, blame the deadline). In my past exchanges with George Fitzgerald of Accelrys, the issue of DMol3 infrared intensities came up as a feature request that would greatly improve the (this) user experience and Dr. Fitzgerald is very interested (of course) in making a great code that much better. Neither code will be disappearing from my toolbox anytime soon.

4. The Periodicity Of The Molecular Solid Doesn’t Care What The Space Group Is – One of the more significant problems facing the assignment-minded spectroscopist is the physical description of molecular motion in a vibrational mode. In the simplest motions involving the most weakly interacting molecules, translational and rotational motions are often quite easy to pick out and state as such. When the molecules are very weakly interacting, often the intramolecular vibrational modes are easy to identify as well, as they are largely unchanged from their gas-phase descriptions. In ionic solids or strongly hydrogen-bonded systems, it is often much harder to separate out individual molecular motions from “group modes” involving the in- and out-of-phase motions of multiple molecules. In the unit cells of molecular solids, it can be the case that these group modes appear, by inspection, to be extremely complicated, sometimes too involved to easily describe in the confines of a table in a journal article.

S-(+)-Ketamine Hydrochloride is one such example where a great simplification in vibrational mode description comes from thinking, well, “outside the box.” The image below shows two cells and the surrounding molecules of S-(+)-Ketamine Hydrochloride. As it is difficult to see why the mode descriptions are complex from just an image, assume that I am right in this statement of complexity. Part of this complexity comes from the fact that the two molecules in the unit cell are not strongly interacting, instead packed together by van der Waals and dispersion forces more than anything else. The key to a greatly simplified assignment comes from the realization that the most polar fragments of these molecules are aligned on the edges of the unit cell.

An alternate view of molecular vibrational motion comes from considering not the contents of the defined unit cell but the hydrogen-bonding and ionic bonding arrangement that exists between pairs of molecules between unit cells. The colorized image below shows two distinct chains (red and blue) that, when the predicted vibrational modes are animated, become trivial to characterize as the relative motions of a hydrogen/ionic-bonded chain. Rotational motions appear as spinning motions of the chains, translational motions as either chain sliding motions or chain breathing modes. It appears as a larger macromolecule undergoing very “molecular” vibrations. In optical vibrational spectroscopy, selection rules and the unit cell arrangement do not produce in- and out-of-phase motions of the red and blue chains, as only one “chain” exists in the periodicity of the unit cell. In neutron scattering spectroscopy, these relative motions between red and blue would appear in the phonon region. This same discussion was had, in part, in a previous post on the solid-state terahertz assignment of ephedrine (with a nicer picture).

So, look at the cell contents, then see if there’s more structure than crystal packing would indicate. It greatly simplifies the assignment (which, in turn. greatly simplifies the reader’s digestion of the vibrational motions).

Patrick M. Hakey, Damian G. Allis, Matthew R. Hudson, Wayne Ouellette, and Timothy M. Korter

Department of Chemistry, Syracuse University, Syracuse, New York 13244-4100

Abstract: The terahertz (THz) spectrum of (S)-(+)-ketamine hydrochloride has been investigated from 10 to 100 cm-1 (0.3-3.0 THz) at both liquid-nitrogen (78 K) and room (294 K) temperatures. Complete solid-state density functional theory structural analyses and normal-mode analyses are performed using a single hybrid density functional (B3LYP) and three generalized gradient approximation density functionals (BLYP, PBE, PW91). An assignment of the eight features present in the well-resolved cryogenic spectrum is provided based upon solid-state predictions at a PW91/6-31G(d,p) level of theory. The simulations predict that a total of 13 infrared- active vibrational modes contribute to the THz spectrum with 26.4% of the spectral intensity originating from external lattice vibrations.

pubs.acs.org/journal/jpcafh
www.somewhereville.com/?p=29
www.somewhereville.com/?p=26
www.somewhereville.com/?p=126
en.wikipedia.org/wiki/Density_functional_theory
en.wikipedia.org/wiki/Ketamine
www.crystal.unito.it
accelrys.com/products/materials-studio/quantum-and-catalysis-software.html
en.wikipedia.org/wiki/Time_domain_terahertz_spectroscopy
en.wikipedia.org/wiki/Computational_chemistry
accelrys.com
en.wikipedia.org/wiki/Inelastic_neutron_scattering
en.wikipedia.org/wiki/Vibrational_spectroscopy
www.somewhereville.com/?p=680

The Vibrational Spectrum Of Parabanic Acid By Inelastic Neutron Scattering Spectroscopy And Simulation By Solid-State DFT

Sunday, February 21st, 2010

Available as an ASAP in The Journal of Physical Chemistry A. As a general rule in computational chemistry, the smaller the molecule, the harder it is to get right. As a brief summary, parabanic acid has several interesting properties of significance to computational chemists as both a model for other systems containing similar sub-structures and as a complicated little molecule in its own right.

1. The solid-state spectrum requires solid-state modeling. This should be of no surprise (see the figure below for the difference in solid-state (top) and isolated-molecule (bottom)). This task was undertaken with both DMol3 and Crystal06, with DMol3 calculations responsible for the majority of the analysis of this system (as has always been the case in the neutron studies reported on this site).

2. The agreement in the hydrogen-bonded N-H…O vibrations is, starting from the crystal structure, in poor agreement with experiment. You’ll note the region between 750 and 900 cm-1 is a little too high (and for clarification, the simulated spectrum is in red below). According to the kitchen sink that Matt threw at the structure, the problem is not the same anharmonicity one would acknowledge by Dr. Walnut’s “catalytic handwaving” approach to spectrum assignment (Dr. Walnut does not engage in this behavior, rather endeavors to find it in others where it should not be).

3. The local geometry of the hydrogen-bonding network in this molecular solid leads to notable changes in parabanic acid structure that, in turn, leads to the different behavior of the N-H…O vibrational motions. There is one potentially inflammatory comment in the Conclusions section that results from this identification. The parabanic acid molecule is, at its sub-structure, a set of three constrained peptide linkages that under go subtle but vibrationally-observable changes to their geometry because of crystal packing and intermolecular hydrogen bond formation. This means that the isolated molecule and solid-state forms are different and that peptide groups are influenced by neighboring interactions.

So, why should one care? Suppose one is parameterizing a biomolecular force field (CHARMM, AMBER, GROMOS, etc.) using bond lengths, bond angles, etc., for the amino acid geometry and vibrational data for some aspect of the force constant analysis. The structural data for these force fields often originates with solid-state studies (diffraction results). This means, to those very concerned with structural accuracy, that a geometry we know to be influenced by solid-state interactions is being used as the basis for molecular dynamics calculations that will NOT be used in their solid-state forms. Coupled with the different spectral properties due to intermolecular interactions, the description being used as the basis for the biomolecular force field likely being used in solution (solvent box approaches) is based on data in a phase where the structure and dynamics are altered from their less conformationally-restricted counterpart (in this case, solid-state).

A subtle point, but that’s where applied theoreticians do some of their best work.

Matthew R. Hudson, Damian G. Allis, and Bruce S. Hudson

Department of Chemistry, 1-014 Center for Science and Technology, Syracuse University, Syracuse, New York 13244-4100

Abstract: The incoherent inelastic neutron scattering spectrum of parabanic acid was measured and simulated using solid-state density functional theory (DFT). This molecule was previously the subject of low-temperature X-ray and neutron diffraction studies. While the simulated spectra from several density functionals account for relative intensities and factor group splitting regardless of functional choice, the hydrogen-bending vibrational energies for the out-of-plane modes are poorly described by all methods. The disagreement between calculated and observed out-of-plane hydrogen bending mode energies is examined along with geometry optimization differences of bond lengths, bond angles, and hydrogen-bonding interactions for different functionals. Neutron diffraction suggests nearly symmetric hydrogen atom positions in the crystalline solid for both heavy-atom and N-H bond distances but different hydrogen-bonding angles. The spectroscopic results suggest a significant factor group splitting for the out-of-plane bending motions associated with the hydrogen atoms (N-H) for both the symmetric and asymmetric bending modes, as is also supported by DFT simulations. The differences between the quality of the crystallographic and spectroscopic simulations by isolated-molecule DFT, cluster-based DFT (that account for only the hydrogen-bonding interactions around a single molecule), and solid-state DFT are considered in detail, with parabanic acid serving as an excellent case study due to its small size and the availability of high-quality structure data. These calculations show that hydrogen bonding results in a change in the bond distances and bond angles of parabanic acid from the free molecule values.

pubs.acs.org/doi/abs/10.1021/jp9114095
pubs.acs.org/journal/jpcafh
en.wikipedia.org/wiki/Computational_chemistry
accelrys.com/products/materials-studio/quantum-and-catalysis-software.html
www.crystal.unito.it
en.wikipedia.org/wiki/Anharmonicity
chemistry.syr.edu/faculty/walnut.html
en.wikipedia.org/wiki/Hydrogen_bond
en.wikipedia.org/wiki/Peptide
en.wikipedia.org/wiki/Force_field_%28chemistry%29
www.charmm.org
ambermd.org
gromacs.org
en.wikipedia.org/wiki/Molecular_dynamics

Obligatory


T R P Nanosys


Syracuse Astro


Ubuntu 4 Nano


N-Fact. Collab.

Nano Gallery

  • nano gallery
  • nano gallery
  • nano gallery
  • nano gallery
  • nano gallery
  • nano gallery
  • nano gallery
  • nano gallery
  • nano gallery
  • nano gallery
  • more...
  • Aerial Photos

  • aerial photos
  • aerial photos
  • aerial photos
  • aerial photos
  • aerial photos
  • aerial photos
  • aerial photos
  • aerial photos
  • aerial photos
  • aerial photos
  • more...
  • Syracuse Scenes

  • syracuse scenes
  • syracuse scenes
  • syracuse scenes
  • syracuse scenes
  • syracuse scenes
  • syracuse scenes
  • syracuse scenes
  • syracuse scenes
  • syracuse scenes
  • syracuse scenes
  • more...