Background

Assessing the binding potency of millions of different compounds
against a common protein target presents many challenges in a wet
laboratory. Even with advanced lab robotics and automation platforms,
wet laboratory high-throughput screening (HTS) is still hindered by
time, resources, and finance.

As such, it is now commonplace to first consult structure-based drug
design (SBDD) software platforms in attempt to filter out any compounds
that are predicted to have low binding potency against the target of
interest. This, in turn, frees up time, resources and finance as
subsequent wet laboratory screening is only considering compounds we
believe to have most promise.

Flare™ is
a drug-design software platform which can be used to predict whether
compounds will bind to a given protein target, using techniques such as molecular docking and Electrostatic Complementarity™ (EC) analysis
which compares the electrostatics of the bound putative ligand with the
binding site of the target protein. The Flare GUI can seamlessly run
such calculations on a maximum of 10,000 ligands, but what happens when
we have millions of ligands? In such a case, it would not be possible to
store all these structures and run calculations within the GUI. This is
where using the pyflarecommand line can be very useful.

pyflareis a
dedicated Python binary that gives you access to Flare functionalities
from within a command line Python script. This is incredibly powerful as
it enables Flare calculations, such as docking and EC, to be executed
on a command line interface on your dedicated computing cluster without
the need to load and store molecules and proteins within the GUI. With
the GUI capacity limitation lifted, this leaves us free to perform Flare
calculations on much larger compound libraries.

This article will talk you through how we can write a pyflare  script to execute a docking and EC hit identification workflow. The real benefit of such a pyflarescript
is the transferability: the script can be executed on the different
drug discovery targets or compound libraries without being altered. To
launch this workflow, all that you need are three files: your protein
target stored as a PDB file, the associated crystallographic ligand
stored as an SDF file, and the compounds you wish to test in your
computational HTS also stored within a single SDF file.

So, let’s look at what such a hit identification workflow looks like when encoded into apyflarescript.

Initializing the project

The first thing we need to do is to create a new, blank, Flare project.

Figure 1. Creating a new Flare project using pyflare.

In the example shown in Figure 1, the Flare
software library is first imported, giving us access to Flare’s
calculation methods; a Flare project is then created and assigned to the
variable ‘project’. By assigning the Flare project to a variable, it
allows the user to perform operations such as saving your progress
(i.e., using Python to perform ‘File > Save As’), clearing the
current project data (i.e., using Python to perform ‘File > New
Project’) as well as reading in and exporting structure files.

Loading protein and crystallographic ligand data

Now we have created a new Flare project, let’s begin by loading our
prepared protein and crystallographic ligand into the project. These
structure files are located in the working directory where we are
executing thispyflarescript.

Figure 2. Loading protein and ligands into a Flare project.

The proteins.extend()and ligands.extend()
operations enable us to load our prepared protein target and
crystallographic ligand into the Flare project. Conceptually, we are
‘extending’ the current ligand and protein tables through the addition
of these new structures. As these new structures are stored within
separate files, we need to tell Flare to read these files using the ‘flare.read_file()’ command. The exact same procedure is also done for the SDF file encompassing our compound library.

Now that all the necessary structures for our hit identification
workflow have been loaded into the Flare project, we now assign the
target protein and library compounds to separate variables to be used in
subsequent calculations. In this example, we only have a single protein
loaded into our Flare project, which occupies the first data entry in
the proteins table. As such, we can assign it the variable target_protein by calling the zero-index (1st entry) of the proteins table. For docking and EC calculations in pyflare, the ligands to be used in the calculation must be stored within a list format. Therefore, by creating the ligand_designs
variable as an empty list, we can then use a for-loop to iterate
through all the ligands loaded into the Flare project, assigning each to
the ligand_designs list variable. The ligands are now stored in the correct format for calculation.

Preparing ligands for docking

The next stage in our hit identification workflow is to ensure all
our ligand designs have adopted a minimum energy 3D conformation prior
to docking. This can be done by performing a pop-to-3D calculation
followed by a minimization calculation (Figure 3).

Figure 3. Pop-to-3D and minimization of ligand structures to be docked.

By establishing a for-loop, we can loop through all the ligands designs we have stored in the ligand_designs
list variable and execute the pop-to-3D and minimize methods to obtain a
low energy 3D conformation for each ligand design. Note the use of ligand_designs[1:], which allows us to skip the first entry within the ligand_designs list variable. Recall, the crystallographic ligand is the first entry within the ligand_designs list variable and therefore does not need converting to 3D.

The docking experiment

Now the ligand designs have adopted a low energy 3D conformation, we
are at the right stage to carry out the first of our HTS experiments:
docking (Figure 4).

Figure 4. Setting up the docking experiment where our input ligands will be docked into the target protein.

We first initialize the docking operation through loading the docking class using flare.Docking(). We specify which ligands we wish to dock (e.g., our ligand_designs
defined above), and which protein we wish to dock into. When
configuring a docking calculation in Flare, we also need to be sure to
define the energy grid box. This should be defined to encompass the
target protein active site and, as such, can be designated as a region
of 3D space that hosts the crystallographic ligand. This was the first
ligand we loaded into the Flare project and, therefore, occupies the
zero-index (1st entry) within the ligands table. In our example, we are
only interested in storing the most stable binding pose for each ligand,
which is achieved by setting max_poses
to 1. We can customize the docking calculation quality by loading a
defined pre-set tolerance. In our case, we have chosen ‘Normal’
tolerances. We can now start the docking calculation and wait for the
result.

Calculating the EC for the docked poses

Once the docking is complete, the next stage is to calculate the EC score for each of the docked poses.

Figure 5. Calculating the Electrostatic Complementarity of the docked poses against the target.

After completion of the docking
calculation, each ligand within the Flare project is now being stored in
its most stable docked pose and is associated with the protein target.
Therefore, upon initializing the EC calculation, through loading the
Electrostatic Complementarity class using flare.ElectrostaticComplementarity(),
we simply need to direct Flare in the direction of the ligand_designs
variable which holds all the relevant information necessary to execute
an EC calculation. We can now launch the EC calculation and wait for the
result.

Triaging the results

The final step in the workflow is to determine which ligands show
promise as good binders to the protein target. This is achieved by
triaging the compound library by means of the docking and EC scores.
Explicitly, we are interested in triaging based on the LF Rank Score
(calculated during the docking experiment) and the EC Score (calculated
during the EC experiment). A molecule with a good LF Rank Score and EC
score could potentially be a strong binder to the target, and these are
the ligands we wish to isolate from our dataset.

Figure 6. Analyzing and extracting the best-scoring result.

Firstly, let’s create a new project role, aptly named Most stable pose,
to store the triaged ligands. We can first check whether such a role
has been created before, and if so, delete and create it again.
Otherwise, just create the role as normal using project.roles.append().
Now, for all the ligands in our project, we wish to identify the
ligands with good EC scores (> 0.25) and good LF Rank Scores (<
-9.00). Notice how we have used a nested if-statement here to ensure we
are capturing ligands that satisfy both triage criteria. For each ligand
we identify in our project that satisfies these criteria, we can store
it in the Most stable pose role created above.

Finally, we can export our triaged ligands to a separate SDF file using flare.write_file() and save the project in case we wish to come back to it later.

Conclusions

In this article, we have executed a hit identification workflow, scripted in Python using thepyflarebinary,
which has performed a computational HTS on a very large compound
library and identified which ligands show promise as good binders to a
target protein. This was all done without the need to load and store
structures in the GUI. This workflow is completely automated,
computationally fast, and does not suffer from GUI (data) storage
limitations and is all made possible by the Python extension to Flare.

Request a free evaluation of Flare today to further explore its full portfolio of molecular modelling capabilities.