Application of Site-Identification by Ligand Competitive Saturation in Computer-Aided Drug Design

A. SILCS Simulations

SILCS simulations are initially performed to pre-compute an ensemble of the solutes, water and protein conformations, yielding the SILCS FragMaps, which are used for the subsequent SILCS calculations. The approach was developed to allow for competition between all the solutes and water for the entire region occupied by the target macromolecule. The simulation systems consist of a protein or other macromolecule and randomly distributed solute molecules at a concentration of 0.25 M immersed in an aqueous solution of 55 M water. The most common solute molecules include benzene, propane, methanol, formamide, dimethyl ether, methylammonium, imidazole, and acetate in the standard version of the SILCS methodology. The solutes represent apolar (aromatic and aliphatic) groups, hydrogen bond donors and acceptors, and charged functional groups. Moreover, alternate solute molecules such as fluoroethane, trifluoroethane, chloroethane, fluorobenzene, chlorobenzene, and bromobenzene may be considered for accommodating halogen atoms, along with the dimethyl ether and methanol, in the halogen “SILCS-X” version.⁶²

The SILCS simulation protocol includes a hybrid oscillating excess chemical potential (μ_ex) grand-canonical Monte Carlo and molecular dynamics (GCMC-MD) approach in order to improve the overall sampling of the solute molecules and water. The GCMC moves comprise molecular translations, molecular rotations, and insertion or deletion of the solutes and water. All these moves are accepted or rejected based on the Metropolis criteria. The use of the oscillating μ_ex involves altering μ_ex to attain adequate solute insertion acceptance rates, representing an enhanced solute sampling approach. Further details of the GCMC method can be accessed in our previous work.^{61, 62} The MD simulations provide additional conformational sampling of the protein, solutes and water in the NPT ensemble at a system temperature and pressure of 298 K and 1 atm, respectively. The backbone alpha carbon atoms of the protein are typically restrained with a force constant of 0.12 kcal mol⁻¹ Å^-2. The timestep for the MD simulation is 2 fs with periodic boundary conditions. The simulation box size is extended 15 Å beyond the proteins in all three dimensions. For transmembrane proteins, in the SILCS-Membrane tool, the principal axis of the protein is aligned with the bilayer normal (z-axis) and relocates the protein to the center of the bilayer. The GCMC-MD simulations are performed in a sequential manner for a total of 126 cycles, composed of 26 equilibration cycles in which the solutes and water are redistributed on the simulation system followed by 100 production cycles. Each cycle involves 200,000 GCMC steps, 5,000 steps of steepest descent minimization, 100-ps equilibration MD, and a 1-ns production MD simulation. The SILCS simulation protocol recommends ten independent simulations for each target in order to obtain adequately converged ensembles of the distributions of the solutes and water and the conformations of the protein. The force field parameters for the protein, solute molecules and ligands are obtained from the additive CHARMM36m and CHARMM general force field (CGenFF).^67–70 The GCMC code and additional workflows were initially developed in the MacKerell laboratory combined with various GROMACS⁷¹ utilities, which is also used for the MD simulations. In addition to studies of soluble proteins, the SILCS-Membrane module is available that allows for SILCS simulations of proteins, such as GPCRs and transporters, in an explicit membrane environment. The SILCS Software Suite is available from SilcsBio LLC and is free to not-for-profit researchers.

B. FragMap Construction

From the SILCS GCMC-MD simulations probability distributions of the solutes and water are calculated. The 100 1-ns production MD trajectories from each of the 10 independent simulations, yielding a total of 1 microsecond of MD, are used for producing FragMaps by binning selected atoms of the solute molecules into 1 Å x 1 Å x 1 Å cubic volume elements (voxels). The resulting occupancies are then normalized for the number of simulation snapshots and for the concentration of the solute atoms in the bulk solution. The bulk concentration is typically based on the average number of solute atoms relative to total number of water molecules at a 55 M concentration but may also be calculated based on an assumed concentration of 0.25 M or based on the total volume of the simulation system.⁷² The normalized occupancies are then converted to grid free energy (GFE) values by using the Boltzmann transformation. Therefore, GFE is an atom-based score representing free energy of the functional groups represented by the solutes relative to being in the bulk aqueous solution. Accordingly, the GFE FragMaps can assume both negative and positive free energies indicating favorable and unfavorable interactions with the protein, respectively, with zero kcal/mol indicating a free energy equivalent to being in aqueous solution. Along with the FragMaps, an “exclusion map” is generated based on regions not sampled by water or any solute molecules during the entire production run. This accounts for the conformational flexibility of the protein indicating the extent to which ligands binding to the protein can penetrate under the traditional solvent accessible surface. Thus, once the SILCS simulations are completed and the FragMaps and exclusion maps calculated, the protein structure itself is no longer used for subsequent SILCS calculations, though the structures used to initiate the simulations and from the SILCS MD simulations themselves are available for visualization and other purposes.

Regarding classification of solute molecules and associated GFE FragMaps, there are four main atom classification schemes (ACS): Generic Apolar Scale (GAS21), Generic Apolar Scale with Halogens (GAX21), Specific Standard (SS21), and Specific Halogen (SX21). In the Generic ACS, several “specific” FragMap types are grouped together to yield “generic” FragMaps. The formamide nitrogen (FORN) and imidazole donor nitrogen (IMIH) voxel distributions were merged to yield a generic hydrogen bond donor (GEND) FragMap, while formamide oxygen (FORO), dimethyl ether oxygen (DMEO), and imidazole acceptor nitrogen (IMIN) distributions are combined into a generic hydrogen bond acceptor (GENA) FragMap. Generic FragMaps also include apolar FragMaps based on the carbons of benzene and propane, heterocycle carbons of imidazole (GEHC), and methanol oxygen (MEOO) distributions. The charged functional groups include a methylammonium nitrogen (MAMN) atom type for positively charged groups, and an acetate oxygens (ACEO) atom type for negatively charged groups. In the Specific ACS, all specific FragMaps are used, supplemented by generic maps. Additionally, both the Generic and Specific ACS can be supplemented with the halogen FragMaps obtained from SILCS-X simulations.⁶²

The convergence of the FragMaps is verified by calculating a metric named the overlap coefficient between the first five (1 to 5) and the remaining five (6 to 10) independent SILCS simulations. Higher values (ranges between 0 to 1) of OC confirm the convergence of the FragMaps. Additional details regarding the calculation of FragMaps, overlap coefficient, GFE, and ACS can be found in these references.^{55, 61, 62} The FragMaps obtained from the SILCS simulations can be used in a qualitative manner to assist ligand design by visualizing 3D regions where fragments can be placed and linked to make drug-like molecules or ligands can be modified by adding or removing functional groups for improving binding affinity and specificity. An example of the PTP1B FragMaps visualization is shown in Figure 1A where the FragMaps are overlaid on the crystal structure of the protein (PDB 2QBS). As shown, a variety of FragMaps occupy the active site of the protein (dotted oval). In addition, the distribution of FragMaps around the entire protein indicates the potential presence of additional putative binding sites on the protein; such sites may potentially act as allosteric binding sites. Figure 1B shows the FragMaps overlaid on the ligand from the PTP1B crystal structure at the binding site. The apolar FragMaps clearly replicate the presence of both aromatic rings. The representation of hydrogen bond donor and acceptor maps shows the possibility of such groups in that region. In addition, the carboxylate groups of the crystal ligand are accurately captured by the ACEO FragMaps as shown on the left side of the figure. The acetate group makes multisite interactions with the surrounding protein residues including lysine and arginine, along with waters. This example shows the capability of SILCS FragMaps in capturing prominent molecular interactions or binding hotspots that can be further utilized for several CADD approaches as shown in the later sections of this article. The FragMaps and exclusion maps can be readily visualized using Pymol⁷³ or visual molecular dynamics (VMD)⁷⁴ program as well as in other programs as the maps may be supplied in CCP4 format traditionally used for visualization of electron densities.

C. SILCS Toolset

The SILCS toolset comprises a number of applications (Figure 2) based on the FragMaps obtained from the SILCS simulations. Most of these applications are quantitative in nature, providing information that may be used for various aspects of drug development as well as macromolecular modelling. These include the ability to perform spatial and conformational sampling of ligands ranging from fragment molecules to drug-like molecules to larger molecules such as polypeptides, oligonucleotides and oligosaccharides, termed SILCS-MC. This technology is the foundation for additional tools including SILCS-Hotspots for allosteric site identification and SILCS-Biologics for use in the formulation of mAbs and other protein-based drugs. To facilitate database screening, SILCS may be used to create Pharmacophores, termed SILCS-Pharm. In addition, SILCS-PPI allows for the mapping of protein-protein interactions (PPI) including the full range of possible interactions between protein partners. Details of these methods and associated applications are presented below.

D. Computational Requirements of the SILCS Method

Central to the utility of the SILCS toolset is the huge computational efficiency of the method. This is based on the SILCS approach being based on the use of a pre-computed ensemble of the solute, water, and protein conformations. This pre-computed ensemble, in the form of the SILCS FragMaps and exclusion map, is then used for all the calculations associated with the SILCS toolset. Thus, once the initial SILCS simulations are completed, subsequent calculations take on the order of seconds to hours. For example, SILCS-MC using an exhaustive docking approach may require on the order of minutes on a single core per ligand, whereas a local pose refinement calculation, analogous to FEP methods, requires only seconds. SILCS-PPI calculations are more demanding but still only require on the order of hours. This speed allows for large number of ligands to be evaluated against all possible binding sites on a protein in a matter of hours on a multicore small computer cluster. Concerning the initial SILCS simulations, their computational demands are based on the systems size. In the context of a cluster of 10 GPUs, (e.g. NVIDIA 2080TI GPUs on 8 core Ryzen CPUs or on the cloud such as AWS) for a small protein such as Mcl-1 or Bcl-xL (35–40K atoms) the calculation is completed within one day while a GPCR with full membrane representation (~140,000 atoms) requires approximately 3.5 days. Thus, while the initial SILCS simulations require a significant amount of time, once the SILCS FragMaps are available, they may be used for the full range of SILCS toolset in a rapid, and accurate manner.

A. SILCS-MC Docking

SILCS FragMaps can be utilized to perform rapid ligand docking and compute approximate ligand binding affinities. To calculate the binding free energies, each non-hydrogen atom of the ligand is either classified as one of the FragMap types or assigned as non-classified. All hydrogens are non-classified. The Cartesian coordinates of each classified atom are then used to determine its overlap with a FragMap voxel, with that atom being assigned the GFE value of the voxel. Atoms that are non-classified are assigned a GFE of zero. The Ligand Grid Free Energy (LGFE) is simply the sum of its individual atomic GFEs, which SILCS uses as an approximation of the ligand binding free energy. LGFEs are only approximations of binding free energy due to the lack of additional factors, such as the entropy loss of combining multiple smaller fragments into a larger ligand. The SILCS exclusion map, voxels for which no solute or water atoms were sampled during the simulations, acts as the protein surface and imposes an energetic penalty to prevent the ligand from penetrating into the protein. Figure 1C details the FragMap classifications for an example PTP1B ligand using all four atom classification schemes.

The SILCS-MC docking procedure involves a series of energy minimization, Markov chain Monte Carlo, and simulated annealing steps to determine the lowest LGFE pose of the ligand. The initial ligand pose used by SILCS-MC can either be generated randomly (“exhaustive” or “docking” protocol) or taken from a predetermined set of coordinates (“local” or “pose refinement” protocol). For the docking protocol, the ligand is randomly placed within a sphere of a given radius (ligand placement radius, r_LP) centered at a user-determined coordinate. Previous work has shown the docking protocol to generally perform better than the pose refinement protocol, though at more computational expense.⁶¹ For each docking or pose refinement, multiple SILCS-MC runs are used for each ligand to improve sampling and convergence of the LGFEs. There are several advantages to the SILCS-MC docking procedure for computing ligand binding affinities compared to current alchemical FEP^75–77 methods. First, while some knowledge of the location of the binding pocket is needed to seed to the docking simulations, the exhaustive protocol negates the need for detailed crystal structure coordinates of the parent ligand. Second, each ligand binding affinity and binding pose is independent of any other ligand, so there is no need to perform many relative binding affinity calculations as is required by FEP methods, especially with cycle-closure corrections.⁷⁸ Lastly, LGFE approximates the ligand’s absolute binding free energy, so no experimental data is necessary for computing such quantities, though that data can be used to improve SILCS absolute binding free energies as is done for FEP methods.⁶²

In a recent study, SILCS-MC docking was performed for all four classification schemes, as detailed in Figure 1C on a collection of 8 target proteins. Multiple ligand-placement radii were used for the exhaustive protocol: 1, 2, 5 and 10 Å such that a total of 16 models were run for each protein (four ligand-placement radii for each of the four classification schemes). Several performance metrics were used to compare the computed binding free energies to experimental values, including mean unsigned error (MUE), Pearson’s correlation coefficient (R), and percent correct (PC). The best performing model across multiple target proteins according to this metric was SX21 using the exhaustive protocol with r_LP = 5 Å,⁶² such that this was assigned as the default protocol for SILCS-MC docking if no experimental data is available. If some experimental data is available, the best performing model can be chosen for further testing and model optimization. For these purposes, PC (representing rank-ordering of ligands) was chosen as the desired performance metric. Performance of SILCS for PTP1B ligands⁷⁹ is shown in Figure 3 and Table 1 for a full set of 48 ligands as well as a smaller curated set of 23 ligands first used by Wang et al. with the FEP+ method⁷⁸ and subsequently used by other groups to test AMBER GPU-TI,⁸⁰ Flare⁸¹ FEP, and pmx⁸² (GROMACS utility). The SILCS results for PTP1B compare favorably to these published FEP based approaches that used the smaller ligand set, even when all 48 ligands are used. The Default model performs similarly to Flare FEP⁸¹ (PC = 0.78) and pmx⁸² (PC = 0.68 to 0.76), while the Best PC model approaches the AMBER GPU-TI⁸⁰ results (PC = 0.814). Among these methods, only FEP+⁷⁸ (PC = 0.844) surpasses the performance of SILCS-Best PC. When all 8 targets were considered, the SILCS-MC method is similarly competitive with the FEP methods, though at a computational savings of more than 300 fold.⁶² Moreover, ML optimization of the SILCS model proves to be a particularly powerful tool not currently available to FEP methods, producing an average percent correct score of 0.82, surpassing FEP+’s performance of 0.78 with the OPLS 2.1 force field.^{62, 78}

A machine learning (ML) procedure, Bayesian Markov State Monte Carlo Simulated Annealing (Bayesian-MCSA) allows for improvements in the predictability of the SILCS-MC method when experimental information is available. In the published study, the ML approach was used to further optimize the weighting of each FragMap as it contributes to the atomic GFE values.⁶¹ Using a set of ligand binding poses generated from non-optimized SILCS-MC docking simulations, FragMaps weights were adjusted via Bayesian-MCSA while ligand coordinates were fixed in order to optimize the PC metric. This optimization procedure was performed on five different training sets, each composed of 80% of the available 48 ligands to produce five different PC-optimized models for each ACS (20 total models). Each new model was then used for new SILCS-MC docking on all 48 ligands, with the best performing model (as determined by PC) chosen as the “ML-Optimized” model. While the PC-based ML-optimized model only marginally improves the performance for the full 48-ligand set, for the smaller, curated set of 23 ligands, the performance of SILCS surpasses that of FEP+. A particular advantage of the Bayesian ML SILCS method is the use of the FragMaps as priors such that the optimization may be undertaken with small training sets of ligands as previously reported.

While SILCS-MC can supply information on changes in ligand affinity upon chemical modification at an accuracy competitive with computationally demanding FEP methods, it offers additional information on the contributions of individual atoms in the ligand to the changes in affinity that are not accessible to FEP methods. Atomic GFE values can be used to guide lead optimization by detailing how different moieties in a ligand contribute to the change in binding affinity. Figure 4 details a modification between two PTP1B ligands where a cyclopentyl group is added in order to occupy a hydrophobic pocket. While the new group contributes favorably to binding affinity, its addition shifts the other moieties into less favorable positions, reducing their favorable contributions, thereby limiting the overall gain in affinity. This can be observed in the summed GFE score of −2.19 kcal/mol for the full cyclopentyl group, whereas the GFE score for other rings and nitrogen atom becomes less favorable. Overall, this modification results in a computed ΔΔG_SILCS score of −1.0 kcal/mol, matching well with the ΔΔG_Expt value of −0.87 kcal/ mol. The FragMaps in this pocket further suggest that acetate groups flanking both ends of the ligand may be a new direction for further optimization. Such detail could not be extracted from typical FEP-based methods. The results for PTP1B exemplify the information content of SILCS that is not accessible in other CADD binding affinity methods.

B. SILCS-Pharm

Pharmacophore modelling is a widely used CADD approach that is based on the types and spatial relationship of functional groups that are crucial for interactions of ligands with the protein target.⁸³ The resulting pharmacophore models can be used to quickly screen against a database to identify ligands with functional groups that match the required features and its performance is often superior to ligand docking based virtual screening (VS).⁸⁴ SILCS FragMaps at negative GFE values represent the favorable binding regions for the respective functional groups on a given target. When the desired binding site on the target is confirmed, SILCS-Pharm converts these favorable FragMaps into pharmacophore features using a clustering approach with the resulting target-based pharmacophore used for VS.^{58, 59}

Generation of pharmacophore features from GFE FragMaps by SILCS-Pharm involves four key steps. The first step is selection of voxels at a user assigned GFE cutoff for each type of FragMap within the user-defined binding region. All significant voxels which have more favorable GFE values than the cutoff value are selected within the binding pocket and then used in the second step for clustering. For each FragMap type, voxels that are adjacent to each other are put into the same cluster and form a FragMap feature. At the third step, FragMap features are converted into pharmacophore features. Since there are multiple FragMap types, while pharmacophore feature types need to be generic for VS purposes, FragMap features are then combined and generalized into pharmacophore features for practical use. After generation of the pharmacophore features in the third step, they are prioritized according to their feature GFE (FGFE) values. Since GFE represents the binding strength of a fragment atom at a given position, FGFE, which is the sum over GFEs of all voxels within a FragMap feature cluster, can be used to quantify the significance of a pharmacophore feature.⁵⁹ In addition to the normal pharmacophore features, volume constraints that define the shape of the binding region are constructed based on the SILCS exclusion map. Compared to the traditional volume restraints that are built from a rigid structure, use of exclusion map introduces target flexibility in the definition of the binding site shape in a pharmacophore model and thus can help to lower the number of false negatives in a VS.⁵⁹ At the last step, final pharmacophore hypotheses (i.e., models) are generated by combining selected pharmacophore features for use in VS. Typically, a given binding site will contain 5 to 10 pharmacophore features, numbers that are too large for VS as the number of ligands that match to type and location of large numbers of feature become vanishingly small. Accordingly, the user can select a small set of features, typically 4, to comprise an individual pharmacophore hypothesis for VS. The users can further select features to be included in the final hypotheses with the help of external information such as chemical structures of known binders. Models containing 4 features gave the best performances in general according to our test using the targets from the DUD database.⁸⁵ In practice, multiple pharmacophore hypotheses are selected for a given binding site and each used for individual VS. The top ranked compounds, for example, 10K per hypothesis, are then pooled and subjected to SILCS-MC pose refinement from which the LGFE scores are used for final ranking and selection of compounds for experimental testing.⁸⁶

In our previous study, SILCS pharmacophore models were tested in VS using the DUD data sets for eight protein targets including HIVPR, Fxa, DHFR, FGFr1, P38 MAP, ADA, ER and AmpC.⁵⁹ Enrichment results of SILCS-Pharm on the eight targets were compared with docking results using available docking programs including DOCK,⁸⁷ AutoDock,⁸⁸ and AutoDock Vina.⁸⁹ Better enrichment results were observed in general for SILCS-Pharm models when compared to the docking results along with significantly lower computational costs. One example is shown in Figure 5 for the protein target ER. As shown in Figure 5A, apolar (green meshes) and hydrogen bonding donor (blue meshes) and acceptor (red meshes) FragMaps are seen in the binding pocket of ER. SILCS-Pharm followed by user intervention produced a four-feature pharmacophore model as shown in Figure 5B composed of three hydrophobic features (green) and one joint hydrogen bonding donor and acceptor feature (cyan). This pharmacophore model is consistent with the crystal binding mode of a ligand from PDB entry 3ERT as shown in Figure 5B. When we tested this model against DUD data set for ER, an early stage enrichment factor (EF₁) of 14.3 and an area under the curve (AUC) of 0.81 were achieved, where these values were 6.0 and 0.64 for AutoDock 4 which yielded the best results among all tested docking methods.⁵⁹

C. SILCS-Hotspots

Two areas of particular interest to drug discovery are the identification of novel allosteric binding sites and the use of fragment-based drug design (FBDD). FBDD is a promising methodology to develop novel ligands that target macromolecules.^90–93 Instead of using large, drug-like compounds for screening, FBDD first identifies small fragments that bind to a set of binding “hotspots”. Once the hotspots and promising fragments are identified, chemists design larger, more drug-like molecules either by linking them or adding R-groups to a scaffold that can capture the fragment-protein interaction.^{94, 95} With respect to allosteric binding sites, while such sites have long been exploited for therapeutic agents, interest in identifying such sites has recently increased. Allosteric sites offer the possibility of improved specificity, especially against targets with common substrates, and to minimize side effects with targets that participate in biochemically essential pathways. In such targets, either positively or negatively modulating the activity may be preferable to completely inhibiting the protein activity. Towards identification of allosteric sites, it’s important to account for binding pockets that are not evident in experimental protein structures as well as being amenable to targeting by drug-like molecules.

SILCS-Hotspots is an extension of the SILCS-MC methodology designed to identify fragment binding hotspots that are spatially distributed across the global protein structure including both surface and interior binding sites.⁶³ It is performed by comprehensively partitioning the full 3D space of the protein into a collection of overlapping regions into which each fragment is docked a minimum of 1000 times using SILCS-MC. Once the fragments have been docked on the entire protein, two rounds of spatial clustering for the individual fragments and the overall fragment types are undertaken from which the final hotspots are identified. Accordingly, the method not only identifies sites but also a collection of fragments occupying those sites. This approach has been applied to a set of ~90 mono and bicyclic rings present in drug molecules⁹⁶ and the Astex MiniFrag⁹⁰ set, though user-supplied fragment libraries may be employed. In addition, the methodology can be applied to larger drug-like molecules to identify their binding sites. As SILCS-Hotspots use SILCS FragMaps, which are computed from GCMC-MD simulations with an enhanced sampling strategy, the method can identify a large number of sites on the target protein, including many sites not accessible in crystal structures. SILCS-Hotspots has successfully recapitulated the location of known drug-like molecules in both allosteric and orthosteric binding sites on diverse proteins targets including the nuclear receptor, kinases, enzyme, and GPCRs⁶³ as well as successfully leading to the identification of novel ligand binding and allosteric sites.^{86, 97}

An example is presented in Figure 6 for Androgen receptor (PDB 2AM9). The figure shows potential binding hotspots on the protein surface as van der Waals spheres and are colored based on the average LGFE scores of the fragments occupying each site binding pose. The SILCS-Hotspots method successfully identified the orthosteric binding site, which is occluded and is not accessible in a typical MD simulation. The identification of the internal orthosteric binding site was possible due to GCMC sampling performed during the SILCS simulation where solute molecules are removed/inserted anywhere in the simulation system. This and the use of the SILCS exclusion maps allows for the fragments to be docked into the interior orthosteric site. In addition, a known allosteric site was identified by the SILCS-Hotspots method. In the figure, flufenamic acid, an allosteric modulator, from the crystal structure PDB:2PIX is overlaid on the SILCS FragMaps and Hotspots calculated using PDB:2AM9. Notably, this allosteric site is not discoverable using other static pocket detection algorithms, such as FTMap⁵² or Fpocket,⁹⁸ because the pocket is closed in the crystal structure used for the SILCS simulations (PDB:2AM9) indicating the ability of the SILCS approach to identify deep and occluded pockets. This illustrates the importance of accounting for protein dynamics in finding allosteric pockets.

For both identification of allosteric sites and FBDD our analysis indicates the importance of considering all hotspots, versus just the most favorable sites. For FBDD to work successfully, it is essential to identify fragment binding sites that are spatially adjacent such that fragments occupying those sites may be linked to create drug-like ligands. Current experimental and computational approaches typically can identify only a limited number of fragment sites as well as use a limited number of fragment types. These limitations are particularly challenging as fragments that are of utility for FBDD alone may have binding affinities in the high mM range, making them inaccessible to experimental methods. Importantly, SILCS-Hotspots largely overcomes such limitations.

D. SILCS-PPI

Computational identification of protein-protein interactions (PPI) is challenging due to the large dimensionality of the docking process as well as the need to take into account protein flexibility. SILCS FragMaps, by accounting for protein flexibility and offering grid-based GFE distributions, were therefore investigated for use in PPI, leading to the SILCS-PPI protocol.⁶⁴ A variety of protein-docking algorithms have been developed⁹⁹ with a large portion adopting global searching algorithms that utilize fast Fourier transforms (FFT).¹⁰⁰ The FFT approach expedites the sampling process but requires a scoring function expressed on 3D grids. Since SILCS GFE FragMaps are naturally expressed on 3D grids, our SILCS-PPI protocol adapted the FFT-based algorithm to address the PPI prediction problem. Superior to classic FFT-based PPI prediction methods, the FragMap free energy grids contain intrinsic contributions from protein flexibility and desolvation effects. To take advantage of these grids for PPI docking, the probability distributions of sidechains were extracted from the SILCS simulations, yielding protein probability grids (PPG). The PPG are then assigned to the corresponding FragMap types, such that the spatial overlap of the FragMaps of the receptor protein with the PPG of the ligand protein acts as a rapid, grid-based estimation of the binding between the proteins. Notably, complementary information from the receptor and ligand proteins may be considered as part of the scoring function, from which the correlations between protein A and B are evaluated through FFT for both pairs. That is, FragMaps of protein A with PPGs of protein B and PPGs of protein A with FragMaps of protein B are all being considered. The docking score is summed over all FragMap types and the best score will be found for the PPI pose, which has maximized matches between FragMaps and PPGs of protein A and B.

SILCS-PPI was tested using eight protein-protein complexes mainly from the protein docking benchmark 5.0¹⁰¹ covering a variety of protein types. Competitive performance by SILCS-PPI was observed when compared with three other commonly used PPI docking methods include ZDOCK,¹⁰² RosettaDock,¹⁰³ and ClusPro.¹⁰⁴ Moreover, as SILCS-PPI does global docking and samples the entire receptor surface, alternative PPI binding orientations can be also being considered, with the top 2000 orientations generated by default. The validity of alternative binding poses was confirmed by experimental crystal contacts in some of the test cases. While SILCS-PPI in the context of a single PPI calculation is computationally expensive due to the need to perform the SILCS simulations, once the FragMaps are available, they can be used for additional PPI calculations with other proteins as well as for use in the context of the range of SILCS tools.

E. SILCS-Biologics

Protein therapeutics are an increasingly common drug class that leverages the inherent specificity of proteins, including mAbs, for their targets.^105–108 However, protein therapeutic typically require a concentrated protein formulation, which can lead to aggregation, precipitation, and degradation due to self-association of the therapeutic molecule.^109–112 Minimizing these issues is important to allow for the drug to be delivered intravenously (IV) and to maintain the integrity of the drug during storage. Excipients are small molecules other than the active pharmaceutical ingredients that are added to the formulation of the therapeutic to improve the stability, bioavailability, and manufacturability.^113–115 It has been shown that certain excipients can increase stability and lower viscosity of the protein therapeutics significantly. However, the selection of the excipient often relies on experimental trial and error due to a lack of methodologies to rationally select excipients.

The need for improved approaches for protein formulations led to the development of SILCS-Biologics.^{65, 66} The approach combines SILCS-PPI⁶⁴ and SILCS-Hotspots⁶³ to discover excipients and rationally select excipients using a data-science and machine learning approach. The SILCS-Biologics workflow provides a set of features for each excipient from the SILCS-PPI and SILCS-Hotspots calculations and then identifies features that are most responsible for experimental physico-chemical data, such as the stability or viscosity of a protein therapeutic. SILCS-Biologics has been successfully applied to the National Institute of Standards and Technology mAb (NISTmAb)^{116, 117} along with the CNTO607¹¹⁸ mAb. CNTO607 is a model mAb to examine PPIs and NISTmAb was used to examine excipient−protein interactions taking advantage of the experimental data available on that system.¹¹⁸ Features from SILCS-Biologics based on the NISTmAb Fab for the different excipients include the distribution of binding sites and number of sites for each excipient, predicted affinity of excipients, competition of excipients, and the relation of the interactions to predicted PPI encompassing the entire protein surface. From the validation study, a range of SILCS-Biologics derived data was shown to be predictive against various experimental data. Specifically, the number of favorable excipient binding sites that are involved in predicted PPIs correlate with experimentally determined viscosity for a range of excipients including amino acids and sugars. In addition, a combination of the number of binding sites and the predicted binding affinity is indicated to be predictive of relative protein stability.

Figure 7 shows how chemically similar excipients may lead to differential interaction with a protein therapeutic. In the figure, the Fab domain of the NIST antibody structure is shown in surface representation. The surface of the protein is colored based on the likelihood of being involved in self-PPI (blue (low probability) to red (high probability)). We hypothesize that large amounts of PPI of the protein drug with itself (self-PPI) in a protein therapeutic lead to aggregation and higher viscosity. The figure also shows the binding hotspots of arginine (A) and lysine (B), two chemically similar excipients, on the Fab domain. Although having very similar chemical structures, they show differences in their interaction patterns. The detailed knowledge of how these excipients interact at the surface of the protein therapeutic combined with their self-PPI pattern provides a depth of knowledge to the formulation scientists to select the best excipients that can reduce aggregation in the sample. In addition to qualitative exploration, SILCS-Biologics is of utility for quantitative and predictive excipient design. Comparison of arginine, trehalose, and sucrose with the NISTmAb Fab, all of which give the highest viscosity, indicates that higher viscosities are associated with a low number of predicted binding sites, with lower binding affinity of arginine, leading to its anomalously high impact on viscosity. Such analysis indicates the potential for the SILCS-Biologics approach to be of utility in the rational design of excipients during biologics formulation.

Applications with other co-solvent Methods

Similar to the SILCS applications, other co-solvent based methods have shown potential in the context of different CADD applications. Table 2 summarizes most of the available co-solvent methods and their published applications for CADD purposes. Details of the different methods may vary based on the number and type of co-solvents, sampling approach, force field parameters, and calculations of the distribution of co-solvents. However, the ultimate gain from using such approaches generally is to identify biologically relevant sites on the protein surfaces and using that information in rational drug design. In spite of the technical differences, most of the approaches successfully report identification of crucial binding positions including active and allosteric sites for a diverse collection of proteins indicating a general utility of co-solvent approaches. However, the major success of these methods is their utilization in computing protein-ligand relative binding affinities, hit identification, and lead optimization. Prominently, the methods MDmix, MixMD, and pMD have shown several such capabilities for drug design against a range of protein targets. All of these studies indicate the general utility of co-solvent methods and it is always encouraging to see how different labs leverage such methodologies in innovative applications.

Future Directions and Summary

While the SILCS toolset has a significant number of capabilities, ongoing efforts involve new concepts and applications based on the use of the pre-computed FragMaps. An interesting avenue is the impact of inclusion of the explicit treatment of electronic polarization in the calculation of the FragMaps using the Drude polarizable force field. The first study moving in this direction involved a diverse set of targets and ligands using the SILCS framework with comparison to the results from the additive force field.¹²² Usage of the Drude polarizable force field produces better hydrogen bond donor or acceptor FragMaps including halogen atoms and results in better docking scores as compared to scores of the additive force field. An extension of the SILCS-MC docking approach is in progress to include water-mediated interactions between the protein and ligands for more accurately predicting relative binding affinities and better understanding the contributions of individual water molecules to ligand binding. Another ongoing development is focused on improving the sampling efficiency of the GCMC simulations by enhancing insertion and deletion probabilities via two main approaches: (i) implementing biases into the GCMC procedure, and (ii) implementing a new, more adaptive oscillating-μ_ex algorithm. The biases are built on the work of Woo et al.¹²³ by implementing a cavity-biased search algorithm and a conformational bias. These two biases select the optimal insertion points and conformations and orientations for each fragment. These efforts are being performed in conjunction with an extended oscillating-μ_ex algorithm that responds more strongly and directly to changes in solute concentrations, leading to greater sampling of the solutes in the GCMC simulation and ultimately a more dynamic physicochemical environment around the target protein. While the introduction of biasing methods in GCMC led to increased computational costs, efforts to port the GCMC code to GPUs show promise to significantly improve the speed of these aspects of the SILCS simulations. These advances are also showing significant progress in the application of SILCS to RNA molecules, an emerging class of drug targets. In the end, our laboratory is committed to implementing and validating useful computational tools to facilitate multiple aspects of the drug development process.

Among different co-solvent based methods, SILCS developments have made tremendous progress in the last decade, broadening the types of applications and applying them to a variety of protein targets. The primary advantage over other co-solvent methods is the utilization of multiple solutes in the SILCS simulations, allowing for the calculations of FragMaps for a broad range of chemical functional groups. Once available, the FragMaps may be used in a variety of contexts as described above in the overview of the SILCS toolset. Examples of interesting application of the methods have included the identification of host-based antivirals with the potential to treat a wide range of RNA viruses including influenza, coronaviruses and filoviruses,⁸⁶ the development of a model to predict hERG blockade,¹²⁴ an important problem in drug development, through SILCS simulations of the hERG channel itself, design of lipopolysaccharide analogs targeting TLR4,¹²⁵ and the modelling of interaction of IgG antibodies with the protein Endoglycosidase S2.¹²⁶ In an interesting variation of the SILCS approach, it was utilized to study the permeability of small drug molecules across phospholipid bilayers and their free energy profiles.⁷² Notable application by other laboratories include using SILCS to facilitate the design of covalent inhibitors¹²⁷ and the optimization of kinase inhibitors.¹²⁸ Overall, these various applications further indicate the utility of the SILCS approach indicating it to be a platform approach of wide utility beyond just being used for site identification by ligand competitive saturation.