Optimal drug cocktail design: methods for targeting molecular ensembles and insights from theoretical model systems

Abstract

Drug resistance is a significant obstacle in the effective treatment of diseases with rapidly mutating targets, such as AIDS, malaria, and certain forms of cancer. Such targets are remarkably efficient at exploring the space of functional mutants and at evolving to evade drug binding while still maintaining their biological role. To overcome this challenge, drug regimens must be active against potential target variants. Such a goal may be accomplished by one drug molecule that recognizes multiple variants or by a drug "cocktail"--a small collection of drug molecules that collectively binds all desired variants. Ideally, one wants the smallest cocktail possible due to the potential for increased toxicity with each additional drug. Therefore, the task of designing a regimen for multiple target variants can be framed as an optimization problem--find the smallest collection of molecules that together "covers" the relevant target variants. In this work, we formulate and apply this optimization framework to theoretical model target ensembles. These results are analyzed to develop an understanding of how the physical properties of a target ensemble relate to the properties of the optimal cocktail. We focus on electrostatic variation within target ensembles, as it is one important mechanism by which drug resistance is achieved. Using integer programming, we systematically designed optimal cocktails to cover model target ensembles. We found that certain drug molecules covered much larger regions of target space than others, a phenomenon explained by theory grounded in continuum electrostatics. Molecules within optimal cocktails were often dissimilar, such that each drug was responsible for binding variants with a certain electrostatic property in common. On average, the number of molecules in the optimal cocktails correlated with the number of variants, the differences in the variants' electrostatic properties at the binding interface, and the level of binding affinity required. We also treated cases in which a subset of target variants was to be avoided, modeling the common challenge of closely related host molecules that may be implicated in drug toxicity. Such decoys generally increased the size of the required cocktail and more often resulted in infeasible optimizations. Taken together, this work provides practical optimization methods for the design of drug cocktails and a theoretical, physics-based framework through which useful insights can be achieved.

Figure 1

Representations of model molecules used in this study. Model drug molecules are shown in blue, bound to model targets, in red. Numbered spheres are ones that were allowed to bear charge in creating the target ensembles and drug libraries. (a) Model molecules used for the studies in Sections 3.1–3.3. The yellow dashed line indicates the axis of a rod, as described in the text. (b) Model molecules used for the studies in Sections 3.4–3.5. The target and drug shapes are also shown individually for clarity. In certain designs, target spheres labeled with “*” or “+” were constrained to lie within narrow charge ranges (see text).

Figure 2

Optimal cocktails designed toward six model target variants that differ in charge values at two spheres. Each variant is represented by an ‘x’ in 2-D charge space that represents the charge values at its two spheres. Variants that are covered by the same drug are grouped within the same rectangle. Each rectangle corresponds to a drug molecule in the cocktail that covers the drugs within it. The affinity threshold used was −3.3 kcal/mol. The charge distribution of the drug molecule(s) in each cocktail is listed below the plots. The order of the charge values for both the drug molecules and the variants corresponds to the numbering in the schematic above the plots. The quadrant-numbering convention used throughout this work is indicated with Roman numerals.(a) Optimal cocktail that covers all target variants. (b) Optimal coverage achievable from only one drug.

Figure 3

Optimal cocktails to cover a model mutation space of a target. The space is discretized with target variants whose charge distributions at each sphere are sampled at 0.2e intervals within the space (black x’s). Each variant’s location on the x-axis corresponds to its charge value at sphere 1, while its location on the y-axis is its charge value at sphere 2. (a) Optimal cocktail to cover the entire space. Four drugs were needed, and their charge distributions are indicated. Charges are listed in order (1– 4) according to the schematic shown in Fig. 2. (b) Optimal coverage achieved with one drug. The drug with all spheres uncharged was chosen as optimal. (c–d) optimal coverages achieved with two and three drugs, with their charges indicated.

Figure 4

Comparison of the coverages of two drugs when the mutation space is extended to include target variants with charge magnitudes up to 1.6e on either sphere. Now there are charged drugs, such as drug 2 (green) that cover more target variants than the completely hydrophobic drug (blue). Charges are listed in the order according to the schematic shown in Fig. 2.

Figure 5

(a) Optimal “cocktail” to cover the space of model mutant target variants when drug molecules were allowed to bind in any of their four shape-invariant conformations. Such orientational freedom was used as a model for conformational flexibility of a portion of a drug molecule, and also for multi-mode binding. The cocktail contains only one drug, whose charge distribution is indicated. Charges are listed in the order corresponding to the schematic shown in Fig. 2. (b) The tiles corresponding to each of the four individual orientations of the drug in (a). Note that some target variants are not covered by any individual orientation, but the entropic contribution due to multiple orientations allows for overall coverage.

Figure 6

Schematic representation of a drug molecule’s coverage in the continuum electrostatic framework. The coverage for a hypothetical drug with charge distribution $\overrightarrow{q_d}$ is shown. The paraboloid traces out the ΔG_bind,elec of this drug toward identically-shaped target variants as a function of target charge distribution. If a threshold is set for all variants at ΔG_bind = b, then only those variants contained within the shaded ellipse will be covered by the drug. For a different $\overrightarrow{q_d}$, this paraboloid would be translated both horizontally and vertically; therefore different drugs cover different regions and sizes of target charge space. ${\overrightarrow{q_d}}^T{M}_d\overrightarrow{q_d}$ determines the minimum value of the paraboloid for a drug with charge distribution $\overrightarrow{q_d}$, and so it therefore relates charge distribution to coverage.

Figure 7

Coverages of model drug molecules whose charge distributions are sampled along the eigenvectors of the M_d matrix for this system. Each of the four eigenvectors and corresponding eigenvalues are indicated in a–d. Charges in each eigenvector are listed in order (1 −4) according to the schematic shown in Fig. 2. Red, blue, black, green, cyan, magenta, and red tiles correspond to drug molecules sampled at −3, −2, −1, 0, 1, 2, and 3 times the eigenvector, respectively. If a color is absent from a plot, then the corresponding drug covered no target variants within the space (the affinity cutoff used was −3.3 kcal/mol). Note that drugs sampled along the eigenvector with negative eigenvalue (shown in (a)) covered many variants, and the tiles became bigger as the drugs’ charge magnitudes increased. The opposite was true for eigenvectors with positive eigenvalues.

Figure 8

General trends in cocktail design toward model ensembles as a function of target ensemble properties. (a) Average number of drugs needed in the optimal cocktail vs. the number of variants within the ensemble. Error bars indicate ±2 × (standard error). The barely visible error bars correspond to magnitudes of approximately < 0.1 drug. (b) Fraction of infeasible designs (out of 1000 attempts) vs. the number of variants in the ensemble. A design was infeasible if any ensemble member could not be bound by any potential drug better than the threshold affinity (−7.5 kcal/mol). (c) Average cocktail size vs. the target variant affinity threshold. 20 targets were used in all designs. Error bars indicate ±2 × (standard error). No feasible designs were sampled at a threshold of −9.5 kcal/mol. (d) Fraction of infeasible designs (out of 1000) vs. the affinity threshold. (e)/(f): Average size of optimal cocktail/fraction of infeasible designs (out of 1000) for ensembles with different constraints on their members’ charge distributions. “N” indicates no charge constraint (other than the maximum magnitude at any charged sphere being 0.5e). “B” indicates that the four highly buried, interfacial spheres (1–4 in Fig. 1b) must bear charges between −0.1e and −0.4e. “P” indicates that four partially buried interfacial spheres (5–8 in Fig. 1b) are constrained to charges between −0.1e and −0.4e. In all cases, 20 targets were in each ensemble and the affinity threshold was −7.5 kcal/mol.

Figure 9

General trends in cocktail design toward model ensembles when one of the twenty ensemble members is a decoy. (a) The average size of the cocktail is plotted against the allowed lower-bound free energy of binding toward the decoy. The required threshold affinity toward targets is −7.5 kcal/mol. 1000 designs were attempted, and the error bars indicate ±2 × (standard error). (b) Fraction of infeasible designs (out of 1000) as a function of the decoy lower-bound free energy threshold. Each infeasible design is grouped into one of two categories; “target infeasibility” means that at least one target variant within the ensemble cannot feasibly be bound with appropriate affinity by any potential cocktail member. “Decoy infeasibility” means that all targets in the ensemble can be feasibly covered, but there is no way to cover all targets while still avoiding the decoy.