Research‎ > ‎

Research Blog

Free Energy Calculations of Xchem Fragments Targeting the SARS Cov-2 Main Protease

posted Jul 22, 2020, 3:45 PM by Rajat Kumar Pal   [ updated Jul 27, 2020, 6:23 AM ]

We have obtained computational binding free energy estimates of a series of molecular fragments found to bind the main protease enzyme of the SARS Coronavirus-2. The fragments have been discovered by crystallographic screening as part of the COVID19-Moonshot project ( The community is now designing and testing potential viral inhibitors based on these fragments.

The crystallographic screen did not provide the binding affinities of the fragments. Yet these could be valuable to prioritize designs which include the highest affinity fragments. Indeed we predict here that the fragments display a wide range of affinities. This data is provided in the hope that it is useful in the discovery of SARS Cov-2 inhibitors. 

A brief initial report on this study is available for download.

Testing and Comparing approaches for conformational sampling in Alchemical Binding Free Energy calculations

posted Jun 22, 2020, 11:03 AM by Sheenam   [ updated Jun 22, 2020, 11:07 AM ]

Finding the true binding free energy of a protein-ligand complex remains a crucial challenge in the field of computer-aided drug design. This is the case even for simple systems. One of the key factors to obtaining reliable binding free energy estimates is an adequate sampling of the protein side chains, which if not accounted for, can give an error in binding free energy values by several kcal/mol. To address this problem, here we test a heating & cooling conformational sampling strategy where we heat the system to enable side chains motions and then cool it before production at room temperature. We incorporate this conformational sampling strategy in our alchemical binding free energy estimation protocol and we apply it to a well-known benchmark system: the complex of the L99A mutant of T4 lysozyme with small organic molecules which display varying configurations of protein side chains. We show that unlike the standard calculations at a single temperature, the heating & cooling methodology explores multiple conformations and returns converged binding free energy estimates independent of the starting conformation.

Fig. 3: Orientation of p-xylene with Val111 sidechain in the modeled structure of  T4L99A observed at two main conformations A. trans B. gauche.

Running SDM calculations on XSEDE Comet

posted Apr 5, 2020, 10:29 AM by Rajat Kumar Pal   [ updated May 5, 2020, 1:58 PM by Emilio Gallicchio ]

File transfer from local machine to XSEDE Comet

Transferring of files to Comet from local machine and vice versa can be done using "rsync" or using the "Globus Transfer" utilities

Using rsync:

To rsync files and directories, do the following:

rsync -avz <somefolder> <USER><USER>/

Using Globus:

For more information about how to use Globus to set up a local "Endpoint" for file transfers to Comet, please visit :

Storage Resource on Comet:

The storage/simulation volume is located at
All user directories are located here under rut147/. The storage space for the home directory  (/home/<USER>) is significantly smaller. It is not suitable to keep simulation directories here. Software, applications, scripts etc. can be stored in the home directory. All user simulations files and directories can be kept here and periodically transferred to local machine using rsync or Globus transfer.

Running minimization/thermalization on Comet

Update: the procedure below which is based on one GPU on a shared node, might fail depending on which GPU is assigned to the job. An alternative is posted on our slack channel and will be described here soon.

Minimization and thermalization require only one GPU and cannot be run under the head node (where your files and software are located). For this, a single GPU needs to be requested at the Comet queuing system as follows:

srun --partition=gpu-shared --gres=gpu:1 --nodes=1 -t 01:00:00 -A TG-MCB150001 --pty /bin/bash

This should give you eventually a prompt on one of the gpu nodes with 1 gpu for 1 hour. Navigate to the simulation directory of the complex you want to do minimization/thermalization and do for example:


The runopenmm script that is required to run the minimization/thermalization and as well as SDM is given here:

export LD_LIBRARY_PATH=${openmm_dir}/lib:${openmm_dir}/lib/plugins:$LD_LIBRARY_PATH
${pythondir}/bin/python "$@"

Running SDM Calculations on Comet

Comet uses a queuing system to launch any jobs on its GPU nodes. The queuing system is based on "SLURM" and uses .slurm script to launch jobs.
To run jobs on Comet, each complex directory needs a .slurm script. An example .slurm script is given here:

#SBATCH -J 6vww-1

#SBATCH -p gpu
#SBATCH --gres=gpu:4
#SBATCH -N 1 # This is nodes, not cores
#SBATCH -n 1 # one process per node so we get one entry per node
#SBATCH -t 18:00:00 # Max time allotted for job

echo "Number of nodes: $SLURM_NNODES"
echo "Nodelist: $SLURM_NODELIST"
echo "Number of tasks: $SLURM_NTASKS"
echo "Tasks per node: $SLURM_TASKS_PER_NODE"


scontrol show hostname $SLURM_NODELIST > .slurm_nodes

awk '{ for(i=0;i<4;i++)print $1 ","0":"i",1,centos-OpenCL,,/tmp"}' < .slurm_nodes > nodefile
cd ..

rsync -av --exclude='slurm*.out' ${jobname} /scratch/$USER/$SLURM_JOBID/

cd /scratch/$USER/$SLURM_JOBID/${jobname}
for i in `seq 1 10` ; do
  ./runopenmm ~/software/async_re-openmm/ ${jobname}_asyncre.cntl
  rsync -av * $SLURM_SUBMIT_DIR/

Change the jobname in the slurm file to suit your job. If I have many ligands I put a placeholder for the -J and jobname entries:

#SBATCH -J <job>

and then use sed -i "/<job>/${complex}/g" ${complex}/${complex}.slurm in a for loop running over the directories in the complex/ directory

To launch the slurm script do for example: sbatch <jobname>.slurm

This is how the .slurm script works:

The slurm script runs the simulation in 10 batches.
1. It firsts rsync's the simulation directory to the /scratch/ disk of the node,
2. It then runs the batch,
3. rsync's back to the working directory, and
4. runs the next batch and so on.

So, for example, if you want the asyncre job to run for 960 minutes, the asyncre.cntl file should specify 96 minutes as the WALL_TIME.

Free Energy Analysis on Comet

Prepare your R environment (first time only):

mkdir ~/R/x86_64-pc-linux-gnu-library
cp -r ~u15684/R/x86_64-pc-linux-gnu-library/3.6 ~/R/x86_64-pc-linux-gnu-library/

Load the R package:

module load R

Navigate to the scripts/ directory of the SDM project. For example:

cd /oasis/projects/nsf/rut147/<username>/6vww-1/scripts/

Edit the
file to point it to the top-level SDM project directory (first entry, in this example it would be 
/oasis/projects/nsf/rut147/<username>/6vww-1) and the number of samples to discard for equilibration (last entry), and then run the analyze script:

bash ./

A CentOS Docker Image for OpenMM/SDM Development

posted Mar 29, 2020, 5:11 PM by Emilio Gallicchio   [ updated Apr 3, 2020, 11:25 AM ]

It is difficult to maintain an OpenMM build environment. There are a number of strict requirements for tools, libraries, and compilers that are not easy to satisfy for each operating system. An obvious workaround is to maintain a development build box. However, the build box consumes space and resources and, as during the ongoing emergency, it might not be always available over the network. In fact, for unknown reasons, we lost access to our lab internal network at Brooklyn College and no one can go there to fix it.

A better alternative is to implement the software development environment as Docker image. Docker is a framework to create and run virtual machines, called "images" in Docker's parlance. I summarize here the steps we took to create one such image, and how we use it to build OpenMM and related utilities and prepare molecular systems for the Single Decoupling Method for protein-ligand binding free energy estimation.

The latest image is called egallicchio/centos610-anvil-7. To load it do:

$ docker pull egallicchio/centos610-anvil-7

How to Use the Docker Image to Run the SDM workflow

Run the image:

$ docker run -it egallicchio/centos610-anvil-7

Follow the usual procedures to run the SDM workflow. Skip the minimization/thermalization in the script as there is no GPU in the docker image to do it. Also, for some reason, mae2dms needs the conda libs:

> export LD_LIBRARY_PATH=/opt/conda/lib/:$LD_LIBRARY_PATH
> bash ./

Finally, rsync the working directory to a computational server/cluster to do the minimization/thermalizations and perform the ASyncRE calculations.

Building and Maintaining the Docker Image

We have based our image on the one used by conda-forge, which is based on CentOS 6.10:

$ docker pull condaforge/linux-anvil
$ docker run -it condaforge/linux-anvil

you will be dropped into a
bash shell as user conda. A python 3.7 conda environment stored in /opt/conda is automatically activated. This user is authorized to install packages via sudo and yum. The Red Hat devtoolset-2 is also installed. The sudo command in this toolset conflicts with the system sudo. To run sudo use /usr/bin/sudo explicitly. To gain root do (while the image is running):

$ docker ps

to get the id of the running image, 8d0865c61538, say. Then:

docker exec -it 8d0865c61538 bash

To compile OpenMM we need gcc from devtoolset-6, also some reasonable text editors and such:

> /usr/bin/sudo install devtoolset-6
> scl enable devtoolset-6 'bash'
/usr/bin/sudo install nano emacs
> /usr/bin/sudo install wget rsync

This is how we built msys:

> conda install -c conda-forge boost
> conda install -c conda-forge scons
> mkdir src && cd src
git clone

I edited the SConscript file to add /opt/conda/include to the include path, also added the needed libraries in LIBS. The key section looks like:

if True:
    for p in env['CPPPATH']:
        if p.startswith('/proj') or p.startswith('/gdn'):
            flg.append('-I%s' % p)
    env.Append(CFLAGS=flg, CXXFLAGS=flg)


scons -j4
> scons -j4 PYTHONVER=37
scons -j4 PYTHONVER=37  install PREFIX=$HOME/local

The msys python tools such as dms-info fail complaining about some kind of python/C++ function argument mismatch. Here is a post that explains how to fix the problem if we need to. We mostly care about mae2dms, which works.

The next steps are to gather the tools necessary to compile OpenMM.

After a number of failed attempts related to cmake incompatibilities, I ended up installing cmake 3.6.3 from sources. As root:

# cd ~/src/
# wget
# tar xzvf v3.6.3.tar.gz
# yum install ncurses-devel
# cd CMake-3.6.3
./bootstrap && make && make install

I installed CUDA as root as well. I skipped the installation of the NVIDIA driver since the image does not have an NVIDIA GPU card:

# wget
# rm

The doxygen app packaged with conda-forge is broken. It does not scan directories with header files. It took a while to debug this. The centos version of doxygen is fine:

> /usr/bin/sudo yum install doxygen

The next steps are for building OpenMM. As the conda user:

> cd ~/src
> wget
> tar zxvf 7.3.1.tar.gz
> cd openmm-7.3.1

I modified CMakeLists.txt in src/openmm-7.3.1/wrappers/python/ to do the python installation under local/openmm-7.3.1/lib. Here is an excerpt:

#set(PYTHON_SETUP_COMMAND "install --root=\$ENV{DESTDIR}/")
set(PYTHON_SETUP_COMMAND "install --prefix=/home/conda/local/openmm-7.3.1/")

Next, do the actual build. For whatever reason CUDA_CUDA_LIBRARY needs to be specified explicitly:

> mkdir -p ~/local/openmm-7.3.1 && mkdir -p ~/devel/build_openmm_7.3.1
> cd ~/devel/build_openmm_7.3.1
ccmake -i ../../src/openmm-7.3.1/ -DCUDA_CUDA_LIBRARY=/usr/local/cuda/lib64/stubs/

In the ccmake interface I turned off the C and Fortran wrappers and pointed the installation directory to /home/conda/local/openmm-7.3.1. Then did the usual:

> make install && make PythonInstall

Now OpenMM is installed under ~/local/openmm-7.3.1. I followed similar steps to install the AGBNP and SDM plugins (see README).

> cd ~/src
> git clone
> git clone
> git clone

Etc. The SDM workflow does not need building. For the AGBNP and SDM plugins, the python wrapper CMakeLists.txt was modified to do the installation of python libraries under ~/local/openmm-7.3.1/lib as above.

The OpenMM installation is now ready for shipment for deployment on computational servers:

> cd ~/local
> tar zcvf openmm-7.3.1.tgz openmm-7.3.1
> scp openmm-7.3.1.tgz me@myfavoriteserver:~/software/

On the server, untar the distribution. Then use a launch script (runopenmm) such as:

export OPENMM_PLUGIN_DIR=${openmm_dir}/lib/plugins
export LD_LIBRARY_PATH=${openmm_dir}/lib:${openmm_dir}/lib/plugins:$LD_LIBRARY_PATH
export PYTHONPATH=${openmm_dir}/lib/python3.7/site-packages:$PYTHONPATH
${pythondir}/bin/python "$@"

For example:

 $ ~/software/bin/runopenmm

To finish up the development image, I installed a version of the academic Desmond-Maestro needed by the SDM workflow:

> mkdir -p ~/schrodinger/installers && cd ~/schrodinger/installers
> export SCHRODINGER=~/schrodinger/Desmond_Maestro_2018.4
> scp me@myfavoriteserver:~/software/Desmond_Maestro_2018.4.tar .
> tar xf Desmond_Maestro_2018.4.tar

and then proceed as usual with the Maestro installation.

Finally, exit from the docker image and commit the changes:

> exit
$ docker commit -m "build box" -a "Emilio Gallicchio" <image id> egallicchio/centos610-anvil-7

where <image id> is the id of the docker container one gets from docker ps -a.

For subsequent updates commit with a new tag:

$ docker commit -m "update make" -a "Emilio Gallicchio" <image id> egallicchio/centos610-anvil-7:version2

Protonation of an α-Hydroxytropolone HIV RNase H Inhibitor through QM/MM Methods

posted Feb 1, 2020, 1:42 PM by Judy Wei   [ updated Feb 1, 2020, 1:47 PM ]

α-Hydroxytropolones are potential anti-HIV drugs that inhibit the processing of the DNA/RNA hybrid by the RNase H enzyme. It is known that these molecules bind to the Mg2+ metal ion cofactors of the enzyme. However, the specific protonation state of the bound inhibitor is unclear. The molecule we intend to study is the bound structure of β-thujaplicinol to HIV-1 RNase H. In order to figure out the protonation state of the oxygen substituents, we proposed four different states to find out which one is most consistent with the known crystal structure by the method of QM/MM geometry optimization.
Figure 3. the structure comparison of the calculated molecule to the original structure.
Figure 3. the structure comparison of the calculated molecule to the original structure 

Application of Heating and Cooling Strategy on Trp-cage Folding

posted Jan 31, 2020, 10:42 AM by VJay Molino   [ updated Jan 31, 2020, 10:45 AM ]

In our body, the proteins that are synthesized fold to their native structure at physiological temperature. Heating a protein breaks interaction between the amino acid residues and extreme temperature causes a protein to unfold and be denatured. What if during protein folding the system is heated? How will it affect the folding process? Will this strategy provide better conformational search during folding that will eventually lead to the native structure of the protein?

Figure 2. Line representation (top) and cartoon representation (bottom) extended structure of the peptide Trp-cage from the sequence of PDB ID 1RIJ.


Calculating Free Energy Change by Displacing Water Molecules

posted Jan 6, 2020, 12:10 PM by Joe Z Wu   [ updated Jan 6, 2020, 1:55 PM ]

While it is most common to consider direct ligand-protein interactions, recent research tries to consider also the solvation effect of water displacement in the binding region. Our protein of choice for this project is the bromodomain of “Pleckstrin homology domain interacting protein” (PHIP)3, a small protein module that recognizes acetylated lysines on histones and serves as an important role in the study of regulation of gene expression.We have chosen to study the mediation effects of the specific W1 water molecule located within the bromodomain binding region. By displacing this W1 water molecule from neat water and the bromodomain active site with a bubble potential, we were able to calculate the free energy penalty, and as such, understand the amount of net work necessary to introduce a ligand into the site. The resulting free energy change of displacing W1 in neat water was 1.3 kcal mol-1, while for the bromodomain binding pocket, this value was 0.384 kcal mol-1To understand more extensively the energetic effects of water molecules on pocket stability and ligand binding affinity, further studies using additional water molecules can be considered.

Figure 1. Spherical bubble (translucent) in bromodomain binding pocket surrounded by water.

Maximum Likelihood Inference of the Symmetric Double-Well Potential

posted Nov 10, 2019, 4:19 PM by Solmaz Azimi

The double-well potential serves as an optimal, one-dimensional model for exploring physical phenomena. This project aimed to estimate the probability distribution of the double-well potential fitted to a Gaussian mixture model by maximum likelihood inference. Hypothetical datasets of the quadratic function were generated using smart-darting Monte-Carlo simulations under the Metropolis acceptance criterion. A major challenge in studying free energy landscapes is sampling efficiency, therefore a larger displacement was implemented in the Monte-Carlo simulations that are generally done at smaller displacements. Although a conventional Monte-Carlo approach can accomplish displacement from one free energy minimum to the next, this is done at the expense of accuracy, such that the acceptance rate of the simulation deviates from an optimal 50%. The generated datasets demonstrated close resemblance to normalized Boltzmann Distribution functions at two temperatures. TensorFlow was utilized to obtain the most optimal parameters for observing the generated datasets based on the Gaussian Mixture Model. Our results demonstrated that TensorFlow approximates a set of optimal parameters for the Gaussian Mixture Model that accurately resemble the Boltzmann Distribution at 300 K, however the method is unable to do so with similar accuracy at 2000 K. 

The Standard State in Binding Free Energy Calculations

posted Jul 30, 2019, 4:10 PM by Emilio Gallicchio   [ updated Jul 30, 2019, 7:01 PM ]

Phys. Chem. Intro

Q: What is the standard free energy of binding?

It is the free energy of binding between a ligand and a receptor when they are in an ideal solution at the standard concentration (C=1 M). The standard free energy of binding is usually denoted by ΔGb

Q: how is ΔGmeasured in practice?

There are many ways to measure a standard binding free energy, the most straightforward is to measure the binding constant Kb, or equivalently the dissociation constant Kd= 1/Kb

Q: But wait, the relation above does not make sense because it takes the log of a quantity in concentration units.

In Physical Chemistry equilibrium reaction constants are dimensionless. The dimensionless nature of the equilibrium constant for binding is clear from its definition:

where R is the receptor, L is the ligand and RL is the complex. Changing the units of concentration does not change the binding constant.

Q: yeah ... I am not sure, my friend works in this lab and they measured the Kd of a drug and he said it's 1.8 nM

In fields such as Biochemistry, Medicinal Chemistry, and Biology, equilibrium constants are very confusingly reported with units, sometimes very strange units. Unless they are doing something really odd, such as redefining the standard state, read the above as Kd = 1.8 x 10-9, no units, and you'll be fine.

Q: really, I can change the standard state?

If you want, sure, but maybe don't call it "standard". A standard it's just a reference that people have agreed on, a bit like the origin of a coordinate system. You can certainly define your own reference solution state, I don't know ... 34,512 molecules per cubic millimeter. (As an example, for acid/base reactions, biologists like to use 10-7 M as the standard concentration of the hydronium ion, causing endless confusion.) Anyway, if you do use your own standard do not expect to get the same "standard" free energy of binding and equilibrium constant as the rest of the world, which uses molar units, for the same reason that you can't expect to get the same distance from the origin if you change the origin.

Let's stick to the 1 M standard state, please ...

The Standard State in Binding Free Energy Calculations

Q: Why do I need to report the standard free energy of binding?

The binding free energy depends on the concentrations of receptor and ligand. For example, the binding free energy is zero at equilibrium concentrations. It is meaningless to report a free energy of binding without specifying the concentrations it corresponds to. By stating that a free energy of binding is "standard", we tell people that it corresponds to the standard state at 1 M concentrations.

Q: Do I need to do a simulation at the standard state? That is, do I need to place enough ligand and receptor molecules in the simulation box so that they are at 1 M concentration?

No. You can do a calculation at any concentration you want. As the simulation progresses, monitor the number of complexes, RL, and the number of free R and L to measure their equilibrium concentrations. That will give you the equilibrium constant and standard free energy of binding using the relation above.

However, the binding constant is rarely computed by "counting" in this way. We like to do computational alchemy. It is faster and more reliable.

Q: Okay, I did an alchemical calculation. I decoupled the ligand from the solution then I coupled it to the receptor. I fed the data into MBAR etc. which spit out a free energy. Can I call it a standard free energy of binding?

Generally, no. What you have calculated is the excess component of the free energy of binding. To turn it into a standard free energy of binding, add the ideal term:

where Vsite is the volume of the binding "region". 

Q: Okay, that sounds easy. Where does the quantity ΔGideal come from?

 ΔGideal  is the reversible work for transferring a ligand from an ideal solution at concentration C to the binding site of volume Vsite . Read these papers:
    1. MK Gilson, JA Given, BL Bush, JA McCammon. The statistical thermodynamic basis for computation of binding affinities: a critical review, Biophysical Journal, 72, 1047 (1997).
    2. M Mihailescu, MK Gilson. On the theory of noncovalent binding. Biophysical Journal 87 (1), 23-36 (2004).
    3. E Gallicchio and RM Levy, Recent Theoretical and Computational Advances for Modeling Protein-Ligand Binding Affinities. Advances in Protein Chemistry and Structural Biology, 85, 27-80, (2011). 

Q: But how do I get  Vsite?

Getting Vsite is easy. You have set it yourself. It is whatever volume the ligand is allowed to explore during the process of alchemically coupling the ligand to the receptor.

Q: What do you mean? I did not restrain the ligand.

Then Vsite is the volume of the simulation box. Did you run the simulation long enough so that the ligand visited the whole box? Probably not, eh!? Then the simulation is not converged. Even more troubling is that you also have implicitly defined the complex RL as any conformation in which R and L are in a region of solution of the same volume as the simulation box. This is probably not what you wanted. You probably wanted to measure the standard free energy of binding of the ligand to a specific region of the receptor, not for the whole receptor and even including the solvent. You should have restrained the ligand to that region. Vsite  is the volume of whatever that region is that you meant to consider.

Q: Right, right ... I forgot about it. I did apply a restraint potential to avoid the "wandering ligand" problem when the ligand and receptor are decoupled.

Ok, great. Than you have a Vsite after all. Is the restraint potential based on the CM-CM distance between ligand and receptor atoms? If so, simply integrate the Boltzmann factor of the restraint potential over the three coordinates of space to get Vsite. Or is the restraining potential something more complicated?

Q: Yeah ... more complicated. The orientation of the ligand was also restrained.

No problem. Integrate over the orientational angles as well and get the angular binding site volume Ωsite and add -kB T ln Ωsite⁄8 π 2 to the ideal standard state factor. Read the papers I told you about.

Q: No, no. I restrained the orientation of the ligand by tethering two atoms of the ligand to two receptor atoms.

That is not allowed. It would perturb the intramolecular conformational distributions of receptor and ligand when they are not coupled. Read this paper:

Boresch S, Tettinger F, Leitgeb M, Karplus M. Absolute binding free energies: A quantitative approach for their calculation. J. Phys. Chem. B. 107, 9535–9551 (2003)

Q: I see. But after coupling the ligand to the receptor I turned off the restraints. Does that change  Vsite?

The binding restraints are not meant to be turned off. They are on during the whole coupling alchemical leg and they stay on.

Q: But then I am not simulating the real thing! The real complex does not have restraints!

I never said that you are simulating the real thing. We are doing computational alchemy! The binding restraints define the complex. With the binding site restraints, you are essentially specifying the set of conformations of the receptor-ligand system that form what you call the "complexed" state. You can't change the definition of the complex depending on the alchemical state of the system.

Q: Okay, now I know you finally went off the deep end. You cannot be right because then the standard free energy of binding would depend on how I set the restraints!

Precisely. The free energy of the complex depends on how you defined it. This is not different from, say, having to define the ranges of φ and ψ backbone angles of the alpha helix state in order to compute the alpha-helix population of a peptide. Change the range of the angles, and the population, as well as the corresponding free energy, changes. You need to define the complexed state before you can measure its free energy!

Q: But then I do compare the calculated free energy to the experiment? The experimental system does not have restraints. What does the experiment measure?

Good question and a very interesting topic. Read paper no. 2 above.

Q: Okay, but now I am worried that I have used harmonic restraints centered on a specific pose of the ligand which I got from docking and I am not sure it is the correct bound pose. Is it a problem?

It could be. It is safer to use flat-bottom harmonic restraints for the binding site volume. That way you have more tolerance and a greater chance of including the "correct" pose.

Q: I am starting to get this. But I heard that imposing restraints and then releasing them can improve convergence. Is that wrong?

It is absolutely fine to impose and then release additional restraining potentials along the alchemical path to improve convergence. These do not need to obey the same strict rules for the binding site restraints discussed above. Just make sure that the binding site restraints stay on throughout the alchemical simulation.

Q: Now I understand why many people prefer to compute relative binding free energies (FEP) rather than the standard binding free energies! They do not need to worry about standard states, the binding site volume, and all this stuff!

Actually, all of this applies to relative binding free energies as well. The relative binding free energy is the difference between two standard binding free energies. Because these depend on the definition of the binding site, so does the relative binding free energy. In principle, without restraints, a ligand undergoing an FEP transformation will eventually dissociate and wander around the simulation box adversely affecting the relative free energy estimate. The longer the simulation the higher the chance that this will happen. You want your estimate to improve, not worsen, as you make the simulation longer. Binding site restraints should be used in FEP as well ... but people rarely do it. 

Don't ask. I really don't know why.

Binding Free Energy of Posaconazole in Two Protonation States to a Model of Captisol

posted May 17, 2019, 6:15 PM by Sheenam   [ updated May 19, 2019, 10:50 AM ]

Posaconazole is a triazole antifungal drug that is used to treat invasive infections by Candida species and Aspergillus species in severely immunocompromised patients. However, its property as a weakly basic and poorly aqueous soluble drug results in poor bioavailability and variable absorption. Hence, for a better intravenous administration of the drug under aqueous conditions, the drug comes in a composition which consists of a solubilizing agent, a modified β-cyclodextrin such as Captisol which supposedly solubilizes with the drug in acidic medium. To test this hypothesis, we calculate the free binding energy of the drug Posaconazole to Captisol in different protonated states to understand the pH dependence of the observed solubility of posaconazole in the presence of Captisol via computational methods. The results of the computational experiments conducted here confirm the hypothesis that this effect is due to the higher affinity of the protonated form of Posaconazole relative to the unprotonated one.

Fig 1: Captisol-Posaconazole Complex (Posaconazole in Green)

1-10 of 35