understanding and exploiting 
 the chemical driving forces responsible 
 for tight and specific molecular interactions remains a fundamental 
 challenge in science and engineering . in practical molecular design 
 applications such as structure - based drug design , a lead compound 
 often serves as a starting point for optimization , in which the compound 
 ( or ligand ) is rationally modified in the context 
 of the target protein ( or receptor ) , to improve its 
 pharmacological parameters , including potency and selectivity .
geometrical 
 and physicochemical complementarity is often required for tight binding , 
 and the quality of this complementarity is generally reflected in 
 the van der waals and electrostatic contributions to binding
. however , 
 whereas the van der waals binding contribution is generally strongly 
 favorable , the net electrostatic contribution is often neutral or 
 even somewhat unfavorable .
ligand electrostatic interactions 
 are acquired through the loss of their individual interactions with 
 solvent in the unbound state , which results in a desolvation penalty . 

the task of ligand optimization often involves incremental improvement 
 of the shape or electrostatic complementarity , although other sources 
 can also be exploited .
charge optimization 
 theory has been developed in the context of 
 linear response theories for the study of electrostatic contributions 
 to binding affinity .
the fundamental observation is that the ligand 
 contribution to the binding affinity can be expressed in a simple 
 mathematical form comprising a desolvation penalty that goes as the 
 square of the ligand charge distribution and a generally favorable 
 screened intermolecular interaction term that is linear in the ligand 
 charge distribution .
the tradeoff of these terms and their different 
 dependencies on the ligand charge distribution result in an optimum 
 charge distribution corresponding to the most favorable electrostatic 
 contribution to binding .
initially , the geometry of both the free 
 ligand and the bound complex were treated as spherical , and the ligand 
 charge distribution was expressed as a set of multipoles in order 
 to solve the problem .
nonspherical geometries and alternative 
 basis sets for ligand charge distributions were then addressed , and the method was extended to exact molecular shapes . a framework for optimizing electrostatic specificity
calculation of the electrostatic potentials 
 is generally done by solving the linearized poisson  boltzmann 
 ( lpb ) equation , but other forms of linear response theory are applicable . 

a series of algorithmic advances has been made to accelerate the calculation and produce more - accurate lpb solutions .
charge 
 optimization theory has been broadly applied to probe or 
 design ligands in various molecular systems .
kangas and tidor applied charge optimization 
 theory to study the binding between the bacillus subtilis chorismate mutase and an endo - oxabicyclic transition - state analogue .
they found that , although the inhibitor showed 
 very good electrostatic complementarity to the enzyme active site , 
 a carboxylate group lost more in desolvation penalty than it gained 
 in interactions with the enzyme ; the calculations suggested that substitution 
 with a nitro group would improve the binding affinity .
mandal and 
 hilvert synthesized the new compound and found , experimentally , a 
 1.7-kcal / mol improvement in binding affinity in a context corresponding 
 to the calculational study , thus identifying the most potent known 
 inhibitor of this enzyme .
other applications 
 include e. coli glutaminyl - trna synthetase 
 binding to its cognate substrates , protein 
 inhibitors of hiv-1 cell entry , the interface 
 between protein kinases and their ligands , small - molecule influenza neuraminidase inhibitors , and the celecoxib ligand binding independently to cox2 
 and caii .
recently , charge optimization 
 and protein design together identified tighter binding peptides to 
 hiv-1 protease that were studied experimentally .
binding specificity optimization probed ligand binding 
 in the model system of hiv protease , other proteases , and their inhibitors .
the charge optimization approaches described 
 above are based on 
 considering the ligand to remain rigid as it binds to a receptor . 
 a comprehensive theory that includes ligand conformational change 
 on binding has been lacking in the field . indeed , as one of the authors 
 has shown , a direct replacement of the unbound ligand structure in 
 charge optimization produced charges that strongly stabilized the 
 bound ( preconformed ) over the unbound conformation , which is physically 
 implausible .
in this study , we explore 
 the electrostatic complementarity of an induced - fit ligand to its 
 receptor , and we generalize charge optimization theory to such cases . 
 for our purposes , an induced - fit ligand has one conformation in the 
 unbound state that may or may not be the same as the preconformed 
 structure in the bound state .
the origin of the previous unphysical 
 results when the rigid theory was applied to flexible ligands is analyzed , 
 and new methodology is developed and applied to the binding affinity 
 optimization problem for a designed hiv-1 protease inhibitor , mit-2-kb-98 .
it generates an optimum partial atomic charge 
 distribution that can be different from that with the rigid - ligand 
 assumption .
the key tenet of the new theory is the constraint that 
 the preconformed ligand can not be lower in free energy than the unbound 
 conformation ( otherwise , it would be the unbound 
 state ) .
( theory ) , we briefly introduce 
 rigid - ligand 
 charge optimization theory and then describe our treatment of the 
 induced - fit charge optimization problem , in which the ligand adopts 
 an unbound conformation potentially different from the bound . in our 
 formulation ,
the unbound ligand is either predefined or selected from 
 a set of candidate conformers .
the latter framing is useful for future 
 extensions of the theory in which the unbound state is treated as 
 an ensemble representing a population distribution .
whereas nonelectrostatic 
 binding contributions cancel for rigid charge optimization , they can 
 persist in the induced - fit theory developed here .
section 3 ( methods ) introduces the test system of hiv-1 protease 
 and its designed inhibitor , as well as methods for calculating continuum 
 electrostatic potentials and nonelectrostatic energies , and approaches 
 to solving the optimization problem . in section 4 
 ( results ) , we analyze the decomposed thermodynamic consequences 
 of the derived results and compare them to those with the rigid - ligand 
 charge optimization approach .
molecular binding in an aqueous 
 solvent can be usefully viewed 
 not as an association reaction , in which only new intermolecular interactions 
 are introduced between receptor and ligand , but rather as an exchange 
 reaction in which some receptor  solvent and ligand 
solvent 
 interactions present in the unbound state are lost to accommodate 
 the gain of receptor  ligand interactions in the bound complex . 

the net effect of the exchange nature of binding is potentially interesting 
 and nonintuitive for electrostatics , and leads to charge optimization 
 theory . here
, we first review the theory for the binding of a rigid 
 ligand to form a complex and then extend the theory to treat induced - fit 
 binding , in each case solving for the ligand charge distribution that 
 leads to the most favorable binding free energy .
consider the binding 
 of receptor ( r ) and ligand ( l ) to form complex ( c ) in an aqueous environment . 

the charge distribution  is modeled as a set of point charges qi at atomic centers ri .
the macromolecules are treated as 
 low dielectric volumes defined by their molecular shapes that are 
 surrounded by high dielectric solvent with ions at physiological concentration . 

the electrostatic binding free energy of the complex can be written 
 as1where gxes is the electrostatic 
 free energy of molecule x , which has a particularly 
 convenient form in any linear response theory . here , we illustrate 
 using solutions to the linearized poisson  boltzmann ( lpb ) equation , 
 without loss of generality:2equation 2 inter - relates 
 positional maps of electrostatic potential  , charge density 
 distribution  , and dielectric constant ( or permittivity )  , 
 where (r ) = [ 8ei(r)/(kt ) ] is the 
 inverse debye length with i the ionic strength , e the electron charge magnitude , k the 
 boltzmann constant , and t the absolute temperature .
the equation is generally solved for the electrostatic potential . 
 in linear response theory ,
the electrostatic free energy of a molecule x is given by3where the summation runs over the partial 
 atomic charges qj located 
 at positions rj .
in the lpb 
 model , individual point charges act independently and obey superposition , 
 so their contributions can be calculated separately and summed,4where i(rj ) ( the electrostatic 
 potential at position rj due 
 to a unit point charge at ri ) can be calculated by solving the lpb equation ( eq 2 ) .
we note that procedures are generally implemented to avoid 
 the coulombic singularity for i = j , but we will leave the singularity in the equations 
 here . by plugging this equation into eq 3 , one 
 can write5after introducing a matrix x for molecule x ( x(i , j ) = ( 1/2)i(rj ) ) and splitting it 
 into solvent 
 reaction field and coulombic terms , as well as expressing the partial 
 atomic charges of x as a vector qx = ( q1,q2,q3, ) , one can rewrite eq 5 as6combining eqs 1 and 6 , the electrostatic binding free energy of a complex 
 can be expressed as7denote the matrices csolv and ccoul as cy and consider 
 their structure .
each element cy(i , j ) = ( 1/2)i(rj ) is half the potential 
 ( solvent reaction field or coulombic , depending on y ) at point rj due to a unit 
 charge at ri .
the locations ri and rj can both be in the ligand , both be in the receptor , 
 or one in each . it is useful to indicate these three cases by representing 
 cy in block form as8 in this format , rcy is a square matrix of dimension nr ( the 
 number of receptor partial atomic charges ) giving half the potential 
 at partial atomic charge locations in the receptor due to the presence 
 of other individual receptor unit charges , and lcy is the same for the ligand .
ccy indicates a half - potential 
 at a ligand atom due to a unit charge on a receptor atom and is not 
 generally square .
because of reciprocity , cy , rcy , and lcy are all symmetric , and ccy and ( ccy ) are indeed transposes 
 of each other . for all cases ,
the subscript c indicates that the potentials 
 are taken in the protein complex . while this does not affect the coulombic 
 potentials ( taken in uniform low dielectric environment ) , it is important , 
 because it defines the molecular boundary for the reaction field potentials . 

the first two terms of eq 7representing 
 the electrostatic free energy of the complex  then can be rewritten 
 as9the factors of 2 in the l 
similarly , for the unbound receptor or 
 ligand , the matrices ry and ly can be represented 
 as rxy and lxy , respectively , which give the potential at a receptor 
 or ligand atom generated by atomic charges in state x ( x = u for the unbound state and x = c for the bound complex ) .
therefore , the electrostatic contribution 
 to the binding free energy in eq 1 can be rewritten 
 as10 when both the receptor 
 and the ligand are regarded as rigid molecules with no conformational 
 change upon binding ( u = c for both molecules ) , the intramolecular 
 coulombic terms cancel ( lccoul = lucoul and rccoul = rucoul ) . moreover
the matrices l and 
 r are expected to be positive semidefinite , because 
 desolvation should represent a penalty and symmetric due to reciprocity . in this rigid  docking
 
 mode with no conformational change , eq 10 can 
 be simplified to11where rt = 2qrtcc = 2qrt(ccsolv + cccoul ) is the screened electrostatic potential ( including both the solvent 
 reaction field and coulombic terms ) at the ligand atoms generated 
 by receptor atomic charges qr . also , we have 
 renamed gbindes to gdockes to describe the rigid case 
 for convenience later . of the three terms in eq 11 , ligand desolvation is quadratic in the ligand charge distribution ql,12the intermolecular screened electrostatic 
 interaction is linear in ql,13and receptor desolvation is a constant , independent 
 of ql:14thus , the electrostatic binding free energy 
 for rigid molecules ( gdockes ) is a paraboloid in the space 
 of ql .
moreover , 
 there is a unique global optimum qlopt = lcctqr that balances the 
 unfavorable ligand desolvation penalty and the favorable intermolecular 
 interactions , and gives the minimum electrostatic binding free energy . 

kangas and tidor showed that this minimum value is upper - bounded by 
 zero , and thus the optimized ligand charge distribution leads to favorable 
 electrostatic binding free energy under certain conditions .
natural proteins seem to use this type of optimization 
 to achieve tight binding : the tight - binding inhibitor barstar was 
 found to be electrostatically optimized to its partner barnase , as 
 described by this charge optimization theory . in this extension 
 of charge
optimization theory , the unbound ligand is treated as a 
 set of individual conformational candidates , one of which is selected 
 by the optimization as a single unbound state conformation . in the 
 current work ,
of these 
 two simplifications , the former has no effect on ligand charge optimization , 
 because the unbound receptor contribution to the binding free energy 
 is a constant for all ligand charge distributions considered ; multiconformational 
 effects for the bound complex could be significant and will be treated 
 in future work .
also considered in future work are ensemble treatments 
 in which each state is represented as a boltzmann distribution of 
 conformations rather than a single conformation . in this framework
, 
 the binding free energy can be deconstructed as the sum of two terms 
 as shown in the blue row of the thermodynamic pathway shown in figure 1 . in the first step , the unbound ligand changes 
 conformation to the preconformed ligand structure ( i.e. , that adopted 
 by the ligand in the bound complex )
; the free energy change for this 
 step is termed the preconformation contribution gpre . in the second step ,
the preconformed ligand and receptor 
 dock rigidly to form the bound complex ( as modeled in the rigid - ligand 
 charge optimization problem above ) , and the free - energy change is 
 termed the docking contribution gdock .
thus , the total binding affinity is represented as the sum of two 
 terms corresponding to the two processes:15each term is expressed as a summation of electrostatic 
 and nonelectrostatic contributions.1617here , the nonelectrostatic terms include covalent 
 energy terms ( bonds , angles , torsions , etc . ) , intramolecular and intermolecular 
 van der waals interactions , and nonpolar solvation interactions .
the process 
 of induced - fit binding of a ligand to its receptor . 
 binding is decomposed into the two steps of preconformation and docking . 

the ligand has two conformations , unbound ( u ) and preconformed ( p ) , 
 which can be the same .
calculation of the binding free energy is indicated 
 schematically , using terms from eqs 20 and 21 .
the blue row indicates solvated molecules and 
 the upper row indicates calculations in uniform dielectric . in the preconformation step , the ligand is solvated 
 in an aqueous 
 environment alone and changes its conformation from the unbound ( u ) 
 to the preconformed ( p ) structure .
we stress that even for the preconformed 
 structure , the ligand is free in solution .
this is equivalent to a 
 special case of eq 10 , in which the receptor 
 is not present .
therefore , the electrostatic contribution to the change 
 in free energy ( gprees ) can be expressed as the sum of a solvent 
 reaction field term ( gpresolv ) and a coulombic term ( gprecoul ) , given by18the nonelectrostatic contribution to the change 
 in free energy is19where exnon - es denotes 
 the nonelectrostatic contribution to the free energy for the ligand 
 in a given conformation ( x ) . in the docking step
therefore,20and21equations 20 and 21 can be interpreted from the thermodynamic 
 pathway of figure 1 .
gpre is equal to the difference in free energy of the ligand 
 in two conformations (  u  and  p  ) . in 
 a given conformation x , the free energy of the ligand 
 includes the cost of assembling the charge distribution ql in a uniform low dielectric ( qltlxcoulql ) and that of moving it to high - dielectric 
 solvent ( qltlxsolvql ) , as well 
 as the nonpolar contributions to the chemical potential ( exnon - es ) . in the original formulation of charge 

optimization theory , the 
 binding reaction was treated as the rigid docking of preconformed 
 ligand and receptor ; with no ligand conformational 
 change ( gpre = 0 for that case ) . 
 here ,
, 
 we permit the unbound ligand to adopt one conformation that may be 
 the same or different from the bound ( preconformed ) ligand conformation . 

charge optimization in the current framing requires optimizing 
 the binding free energy with respect to the ligand charge distribution . 
 in the formulation here ,
a set of candidate unbound ligand conformations 
 is available , and optimization selects an unbound conformation and 
 charge distribution for the ligand that optimizes binding free energy . 

one physical constraint that we explore in this study is that the 
 conformation assigned as the ligand unbound state should correspond 
 to the lowest energy conformation in the set of candidate conformations 
 when evaluated with the optimized charge distribution .
said another 
 way , the unbound structure is constrained to be the ground unbound 
 state . to carry out the desired optimization
, we construct a 
 variational 
 binding free energy expression that includes all terms that depend 
 on the ligand charge distribution , either explicitly or implicitly . 

all terms that depend on the unbound ligand conformation must be kept 
 in the optimization , because they contribute to defining the unbound 
 ligand ground - state energy . here ,
those include all terms in gpre except epnon - es , and , for simplicity , 
 we keep all terms of gpre in our 
 optimization function .
the receptor electrostatic desolvation penalty 
 and the nonelectrostatic contribution to gdock are both independent of ql and , therefore , have been removed to construct the optimization 
 function gtotal , which can be written as22 when the unbound ligand 
 conformation is predefined , the optimization 
 task involves the following:23subject towhere qlk denotes the kth element of ql ( i.e. , the partial 
 atomic charge for atom k of the ligand ) .
the first 
 constraint limits the magnitude of any single charge to qmax ( generally 0.85 units of the electron charge magnitude ) , 
 and the second fixes the total charge of the ligand at qtotal ( usually an integer ) .
if the unbound and bound ligand conformations are the same , this 
 reduces to the original charge optimization approach ; if the two ligand conformations differ , this corresponds 
 to the treatment of gilson .
when , 
 instead , the unbound ligand conformation is selected from 
 a set of candidates , the optimization problem is somewhat different . 

the optimized ligand charge distribution qlopt must be selected 
 such that ( 1 ) gtotal is minimized and ( 2 ) the chosen 
 unbound ligand conformation has the lowest free energy of all conformations 
 in the unbound set with qlopt applied .
the second requirement implements 
 the physical constraint that the conformation claimed for the unbound 
 state be the lowest free energy of those considered , which we show 
 avoids the pathological results observed when rigid theory was applied 
 directly to the nonrigid case . denoting the set of unbound candidates 
 by u and 
 indexing each member by i , we write the following 
 mathematical programming problem:24subject to note that z is not 
 smooth in ql ( i.e.
, z has 
 discontinuous derivatives ) , 
 because of the function min( ) .
the lack of smoothness presents 
 a fundamental challenge in solving the optimization problem and is 
 highly undesirable in practice .
thus , we reformulate the problem into 
 the following smooth one with the same optimum solution , using a standard 
 transformation technique that replaces the explicit min function to 
 define z with an implicit one,25subject to this formulation is general , in that 
 it includes the special cases 
 of rigid ligand ( u containing only the preconformed , 
 bound state ) and a two - state induced - fit 
 ligand ( u containing only an unbound state different 
 than the preconformed , bound state ) .
the 
 formulation is useful in that it can be extended to treat conformational 
 ensembles representing population distributions rather than unique 
 structures for individual states .
also , note that , although the presentation 
 here is in terms of a basis of point charges , the theory is readily 
 applicable to other bases .
the crystal structure of hiv-1 protease 
 alone was extracted from its bound complex with the inhibitor darunavir 
 ( pdb i d 1t3r ) and prepared as described 
 by altman et al . with both catalytic 
 residues ( asp 25 ) deprotonated .
three conserved water molecules were 
 retained , including the so - called  flap water molecule  

hydrogen 
 atoms were added to the complex with the hbuild module of the computer program package charmm .
parameters for mit-2-kb-98 were assigned from 
 the charmm22 set except for the partial atomic charges , which were 
 derived here .
fock ( rhf ) 3 - 21 g level 
 as implemented in the program gaussian 03 .
partial atomic charges were calculated for the new geometry by restrained 
 electrostatic potential ( resp ) fitting to rhf/6 - 31 g * 
 potentials .
a set of 1000 unbound mit-2-kb-98 conformations was computed from 
 a 30 dihedral grid using dead - end elimination ( dee ) and a * to select the 1000 lowest energy 
 structures different from the preconformed design structure .
this 
 unbound conformer set covered the lowest 1.5 kcal / mol in free energy 
 and over 92% of the boltzmann distribution of all the enumerated states . 

the energy function used was the gas - phase energy with a dielectric 
 constant equal to 4 and without nonbond cutoffs .
these structures 
 were indexed by rank in energy from 1 ( lowest ) to 1000 ( highest ) , 
 and , for convenience , the preconformed structure was given an index 
 of 0 .
solvation potential matrices 
 were computed with continuum electrostatic calculations implemented 
 in a locally modified version of the program delphi .
calculations were performed for the preconformed ligand in the free 
 state and complexed with the protein ( lpsolv  lu,0solv and l ) , as well as for each of the other 1000 members of the unbound 
 conformer set alone ( lu , isolv , i = 1 ,  ... , 1000 ) . details of the calculations 
 were as given in the work of altman et al .
partial atomic 
 charges for hiv-1 protease and atomic radii of both the protein and 
 the ligand were obtained from the parse parameter set .
the dielectric constant was set to 4 for the molecular 
 interior and 80 for the exterior , with the dielectric boundary represented 
 by the molecular surface computed with a 1.4- probe .
a 2.0- 
 ion - excluding stern layer surrounded all molecules , and a bulk ionic strength of 0.145 m was employed .
linearized 
 poisson  boltzmann ( lpb ) calculations were performed on a 129 
  129  129 grid with focusing boundary conditions ( successively 
 23% , 92% , and 184% fill ) .
the intramolecular coulombic electrostatic potential matrices lu , icoul were calculated directly from coulomb s 
 law ; 12 and 13 interactions were excluded and 14 
 interactions were scaled by a factor of 0.50 , which are consistent 
 with the charmm22 treatment .
charmm was used with the 
 charmm22 parameter set to evaluate the nonelectrostatic energies eu , inon - es for the 1001 members of the set 
 of unbound candidates , including internal energy terms ( bond , angle , 
 dihedral , improper , and urey  bradley ) and van der waals ( excluding 
 the 12 and 13 interactions , and with special treatment 
 of 14 interactions ; without nonbonded smoothing or truncation ) . 
 the nonpolar contribution to the hydration free energy
the induced - fit ligand charge optimization 
 problem formulated in eq 23 is not guaranteed 
 to be convex , which implies possible multiple local minima . in practice , 

there is no guarantee for a generic solver to determine the global 
 minimum of a nonconvex optimization problem . in this study , we used 
 multistart trajectories with a local minimizer .
we chose conopt3 as the local minimizer for its numerical efficiency , as available 
 through gams .
conopt3 is a nonlinear 
 optimization solver based on generic reduced gradient algorithms . for each problem ,
100 starting ligand charge distributions were 
 randomly generated by uniformly sampling in the feasible space ( the 
 intersection of a hypercube defined by |qlk|  qmax , k = 1 ,  ... , natom and a hyperplane 
 by k=1natomqlk = qtotal ) .
for this purpose , we 
 adopted an algorithm decomposing the space into different types of 
 simplexes and sampled within and across these simplexes uniformly .
the solution was chosen as that with the lowest 
 value of the goal function from the multistart set .
results are presented for induced - fit 
 charge optimization on a 
 previously designed ligand binding to the enzyme hiv-1 protease .
the 
 ligand , mit-2-kb-98 , has an ( r)-(hydroxyethylamino)sulfonamide 
 scaffold functionalized to make favorable interactions in substrate 
 recognition subsites ( see figure 2b ) .
it has a total of 13 rotatable torsion angles 
 ( neglecting the two methyl groups and the amide bond ) , which were 
 varied to create the set of unbound ligand candidates ( see the methods section ) .
two different but related schemes 
 were used here . for each , a collection of 1001 different conformations 
 of the ligand mit-2-kb-98 was used to represent conformational candidates 
 for the unbound ligand .
the structure of the bound complex with the 
 receptor hiv-1 protease was taken as the design structure from the 
 study in which this ligand was designed , synthesized , and assayed .
of the 1001 unbound ligand conformers , conformer 
 number 0 corresponds exactly to the preconformed ( bound ) ligand structure , 
 whereas conformations 11000 correspond to ordered low - energy 
 rotameric variants of increasing energy with dihedral angles selected 
 from a 30 grid , identified with the a * algorithm . in one of 
 the approaches , termed  single - conformer unbound state  , 
 each of the 1001 unbound ligand conformers was considered in turn 
 to individually be the unbound state . with conformer 0
treated as 
 the unbound state , this corresponded to rigid - ligand charge optimization . with conformer 11000 as the unbound state ,
this corresponded 
 to one form of induced - fit charge optimization . in the second approach ,
termed  dominant - conformer 
 unbound state  , an additional constraint was placed on the 
 ligand .
namely , the conformation selected as unbound state was required 
 to be of the lowest free energy in the set of unbound structure candidates , 
 consistent with physical principles . that is , the unbound state is 
 the ground state of the candidate set . in the work presented here , 

the unbound candidate set was grown by progressively adding members 
 to examine size effects .
this second form of induced - fit optimization 
 is the main contribution of the current work .
( a ) designed inhibitor 
 mit-2-kb-98 ( stick representation ; green ) 
 in model complex with hiv-1 protease ( cartoon representation ; cyan ) , 
 with the three retained water molecules ( line representation ) .
( b ) 
 chemical structure of mit-2-kb-98 . to facilitate analysis of the results ,
the overall 
 binding free 
 energy gtotal is considered as the sum of two processes : 
 a preconformation step , in which the unbound ligand adopts the preconformed 
 structure in the absence of the receptor ( gpre ) , and a docking step , in which preconformed ligand 
 and receptor associate ( gdock ) .
the prime symbol 
 indicates that receptor desolvation and nonelectrostatic contributions 
 in the docking step , which are independent of ligand unbound conformation 
 and ligand charge distribution , and , thus , are effectively constant , 
 are neglected .
each member of the collection 
 of 1001 ligand conformations was treated individually as the unbound 
 state , and the rigid - ligand charge optimization approach was used 
 directly to optimize for chemically reasonable partial atomic charge 
 distributions , leading to the most favorable binding free energy . 

the nonconvex optimization problem was solved using 100 multiple starts 
 with conopt3 with the constraint qtotal = 0 ( neutral ligand ) .
all trajectories reached the 
 same solution , which indicates a wide basin of attraction for gtotal , in terms of the ligand partial atomic charges .
figure 3a shows the computed optimal binding affinity gtotal for each of the 1001 individual unbound ligand conformations ( black 
 symbols ) , as well as the preconformation ( gpre , blue symbols ) and docking ( gdock , red symbols ) contributions .
dramatically different results were 
 observed for rigid binding ( structure index i = 0 ; 
 the unbound and bound ligand conformation are identical ) and nonrigid 
 binding ( structure index i  { 1 , 2 , ... , 1000 } ; 
 the ligand changes conformation on binding ) .
the binding affinity 
 improved from roughly 20 kcal / mol for rigid binding to approximately 
 150 kcal / mol for nonrigid cases .
kcal / mol ) and a modestly favorable docking change 
 ( gdock  20 kcal / mol ) ; nonrigid 
 binding involved a surprisingly favorable preconformational gain ( gpre of roughly 300 to 200 kcal / mol ) 
 and an unusually large unfavorable docking change 
 ( gdock of approximately + 100 to + 150 kcal / mol ; 
 see figure 3a ) .
that is , rigid binding charge 
 optimization resulted in optimized charges that were truly complementary 
 to the binding site ; flexible charge optimization resulted in charges 
 that were not especially complementary to the binding site ( because 
 gdock was substantially unfavorable ) but 
 that pathologically recovered a large free - energy benefit from changing 
 the conformation from unbound to preconformed free in solution .
further 
 analysis shows the favorable energetics result from a dramatically 
 favorable coulombic gain ( gprecoul  500 kcal / mol ) 
 and a more moderate desolvation loss ( gpresolv  
 200 kcal / mol ; see figure 4a ) .
although 
 the calculation assumes a given ligand conformation as the unbound 
 state , optimization produces a charge distribution for which another 
 conformation , namely , the preconformed structure , is lower in free 
 energy . that is , the assumed unbound conformation is not the ground 
 unbound state , because , conceptually , the system has a choice of two 
 conformations ( p and u ) .
we reasoned that elimination of this unphysical 
 situation might lead to the charge optimization producing useful and 
 meaningful results for flexible binding , which motivated the dominant - conformer 
 unbound state method described in the next subsection .
contributions to gtotal from three types of 
 optimization : ( a ) the single - conformer unbound state method , ( b ) the 
 dominant - conformer unbound state method , and ( c ) the dominant - conformer 
 unbound state method , with the preconformed structure removed from 
 the unbound candidate set . filled symbols are used for nl = 0 , whereas open symbols are used 
 for all other values .
are on the same scale . energetic and atomic - charge analysis of single - conformer 
 studies : 
 ( a ) contributions to gpre , ( b ) 
 contributions to gdock , ( c ) root - mean - square deviation 
 ( rmsd ) between single - conformer optimum and rigid optimum partial 
 atomic charge distributions , and ( d ) rmsd between single - conformer 
 optimum and nominal charge distributions . before presenting those results , we briefly comment 
 on the mathematical 
 source of the unphysical results .
the ligand desolvation matrix ( the 
 hessian of the electrostatic binding free energy ) for induced - fit 
 charge optimization is l + lpsolv + lpcoul  lusolv  lucoul ( eq 22 ) , which is positive semidefinite 
 for rigid binding ( where lpsolv = lusolv and lpcoul = lucoul and thus 
 only l remains ) but need not be for single - conformation 
 unbound state flexible binding .
eigenvectors corresponding to negative 
 eigenvalues present opportunities for improving ligand affinity without 
 limit but are fixed by constraints placed on the ligand charge distribution . 

these directions produced outsized gains in gpre at the expense of gdock. indeed , 
 examination of optimized charges indicates that many terminated at 
 the constraint and differed substantially from both the rigid optimum 
 charge distribution and the nominal one ( figures 4c and 4d ) .
interestingly , the solvent - screened 
 electrostatic contribution to gdock was 
 favorable ( about 50 kcal / mol ; blue symbols in figure 4b ) but overpowered by the ligand desolvation penalty 
 ( 100150 kcal / mol ; red symbols in figure 4b ) , which resulted in the unusual large unfavorable gdock. these mathematical sources of unphysical behavior correspond 
 to 
 explanations in terms of chemical interaction energetics .
eigenanalysis 
 revealed two mechanisms operating in tandem that were responsible 
 for unphysical behavior .
one mechanism was that atoms whose solvent 
 exposure increased in going from the actual unbound to preconformed 
 unbound state could improve electrostatic affinity by increasing charge 
 magnitude on these atoms and gaining solvation interactions upon preconformation . 

the second was that atom pairs whose interatomic distance changed 
 upon preconformation could gain affinity by accumulating like charge 
 on atoms that became more distant ( or opposite charge on atoms that 
 became closer ) upon preconformation .
we next enforced the 
 notion that the unbound state , if represented as a single conformation , 
 should be the ground state among conformations represented in the 
 set of unbound candidates .
the new procedure involved a different 
 formulation in which the optimization simultaneously chose a ligand 
 partial atomic charge distribution and a presumed unbound ligand conformation 
 from a set of candidates .
the choice was made self - consistently , such 
 that the presumed unbound ligand conformer did have the lowest free 
 energy in the set when computed with the optimum charge distribution , 
 and the optimum charge distribution was selected to minimize the total 
 binding free energy change for preconformation and docking summed . 

if multiple unbound conformers had the same unbound state free energy 
 ( to within 10 kcal / mol ) , they were considered 
 equal contributors to the unbound state , each with the same gtotal. in a post - processing step , the energetic cost to gpre of using degenerate conformers for the unbound 
 state was treated by adding an entropic penalty of rt ln  (  is the number of energetically equivalent conformers ) .
a series of individual optimizations was run using an unbound candidate 
 set ( u ) of increasing size .
i = 0 , the 
 set was of the form u = { 0 , 1 ,  ... , nl}. in successive runs , nl was incremented from 0 to 
 1000 . the continuous and smooth formulation of the optimization in 
 eq 25 was used .
the computed optimal binding affinity gtotal was 
 approximately 20 kcal / mol and did not vary systematically 
 with nl ; this affinity 
 was due to a zero preconformation energy gpre and dominated by gdock (  
 20 kcal / mol ) .
these results are dramatically different from 
 those obtained with the single - conformer unbound state method ( see 
 figure 3a ) .
the new formulation incorporating 
 the ground - state constraint both eliminated the large favorable gpre and maintained a favorable gdock. this latter point is consistent with intuition of electrostatic 
 complementarity .
the preconformation 
 free - energy change was zero , which suggests the preconformed structure 
 was selected either as the only unbound ligand conformation or as 
 one of several .
the coulombic , solvation , and nonelectrostatic contributions 
 to the preconformation free energy were not always individually zero 
 but summed to zero ( different from the single - conformer case ) , which 
 suggests degenerate unbound conformers were sometimes selected .
the 
 docking free - energy contribution was similar to that for the rigid 
 case .
the charge distribution was closer to the rigid optimum than 
 to the nominal distribution , but it could be distinct from either . 
 thus , interestingly ,
much of the energetics of rigid binding was reproduced 
 for this dominant - conformer case , although differences in energetics 
 and charge distributions did result .
the optimized ligand charge distribution qflex,2opt ( when nl = 1000 ) is 
 reported in table 1 , along with that derived 
 with a rigid ligand ( qrigidopt ) and the nominal charge distribution 
 ( q ) .
energetic and atomic - charge analysis of 
 dominant - conformer studies : 
 ( a ) contributions to gpre , ( b ) 
 contributions to gdock , ( c ) rmsd between dominant - conformer 
 optimum and rigid optimum partial atomic charge distributions , and 
 ( d ) rmsd between dominant - conformer optimum and nominal charge distributions .
( a ) indices of selected unbound conformations for incrementally 
 increasing unbound conformer candidate set for the dominant - conformer 
 studies .
( b  d ) some selected 
 conformations of unbound mit-2-kb-98 aligned to each other .
( b ) conformer 
 44 ( cyan ) and its replacement conformer 953 ( gray ) ; conformer 0 ( black 
 sticks ) was aligned to conformer 44 with all atoms .
( c ) conformer 
 907 ( cyan ) and 953 ( gray ) coexisted with conformer 0 as selected unbound 
 conformers .
( d ) conformer 907 and its replacement conformer 996 only 
 differed in the position of a hydrogen atom in the hydroxyl group 
 attached to the benzene ring .
q is the nominal charge distribution , qrigidopt is optimized using a rigid 
 ligand , and qflex,2opt is optimized using a flexible ligand treated 
 with the dominant - state unbound conformer approach using all 1001 
 candidates .
figure 6a indicates the indices 
 of conformations 
 selected as the unbound state for each optimization .
the preconformed 
 structure ( index i = 0 ) was selected in every optimization , 
 almost always with a small number ( up to 2 ) of other conformations 
 of equivalent free energy . as the size of the unbound candidate set 
 was progressively increased , certain structures were selected until 
 replaced by newer ones . after nl reached 50 , conformer 44 was always chosen , along with conformer 
 0 , until replaced by conformer 953 .
conformer 907 joined conformer 
 44 as one of the selected conformations until being replaced by conformer 
 996 .
conformers 44 and 953 ( which replaced conformer 44 ) only differ 
 in the orientation of the hydroxylated benzene ring ( see figure 6b ; upper right ) . in conformer 953 , it was almost 
 in the same plane as the neighboring scaffold atoms n15 , c16 , n17 , 
 and c4 ( labeled ) .
conformer 907 resembled conformer 953 , except for 
 a small torsional angle difference for the benzene ring and in the 
 position of the hydroxyl group hydrogen atom ( figure 6c ) .
these two conformers were selected , together with conformer 
 0 , until conformer 907 was replaced by conformer 996 , which only differs 
 in the position of the same hydroxyl hydrogen atom .
thus , the 
 dominant - conformer unbound state method removed the pathological 
 behavior observed in the single - conformer method through elimination 
 of the unphysical situation of the unbound conformer not being the 
 ground state structure .
interestingly , this produced binding energetics 
 extremely similar to rigid - binding charge optimization , and the preconformed 
 ligand structure was consistently chosen as the unbound ground state , 
 sometimes with one to two other degenerate structures .
indeed , an 
 analysis of the lagrange multipliers from the constrained optimization 
 shows that even when degenerate conformers populated the unbound state , 
 the cost for losing conformer 0 dominated the others substantially . 

the somewhat varying values for contributions to gpre are due to averaging different multiple members of 
 the unbound ligand state for each optimization . taken together , these 
 results suggest that there is little or no affinity benefit to changing 
 conformation upon binding , and that adoption of unbound structures 
 corresponding to complementary geometry in the bound state is sufficient 
 to produce optimal affinity . this suggestion is strengthened by similar 
 results obtained when constraining the total ligand charge to other 
 values ( qtotal  { 1 , + 1 } ) , 
 as well as running similar jobs using a crystal structure rather than the design structure as the bound 
 state and updating structure 0 accordingly ( data not shown ) .
as one final control , we reran the dominant - conformer method with 
 the preconformed structure ( index 0 ) removed from the unbound set , 
 so u = { 1 ,  ... , nl}. this was done to isolate the role of the preconformed 
 ligand structure in the unbound candidate set , because of the many 
 methodological differences between the single - conformer and dominant - conformer 
 methods .
the overall results and detailed breakdown ( figure 3c ) were generally similar to the single - conformer 
 case , with pathologically favorable affinity arising from unphysically 
 favorable gpre values .
these data 
 imply that the preconformed ligand structure should always be considered 
 as a candidate for the unbound conformation , because the alternative 
 not only produces physically unrealistic results but also is physically 
 untrue .
here , we 
 have extended charge optimization theory to treat potentially 
 flexible ligands . one additional concept was used to ensure physically 
 meaningful results  any conformational change from the unbound 
 state to the preconformed bound conformation must be uphill in free 
 energy ( or at least neutral ) .
if it were downhill , then the presumed 
 unbound state would not correspond to the ground state .
the addition 
 of this constraint produces physically reasonable optima and eliminates 
 pathologies that appear in its absence .
the flexible binding 
 case was extended from previous work on rigid 
 binding , in which electrostatics could be clearly separated from nonelectrostatic 
 contributions to the binding free energy .
the theoretical treatment 
 of the flexible case developed here intimately links nonelectrostatic 
 with electrostatic contributions in the preconformational free - energy 
 contribution in eq 22 .
the conformational free - energy 
 surface of the unbound ligand can , in principle , affect the optimization 
 of the unbound conformation and its charge distribution .
note that 
 the parameters for calculating nonelectrostatic contributions are 
 assumed independent of ligand partial atomic charges in the treatment 
 presented here . in practice , these parameters are often calibrated 
 to be compatible with these charges in force - field development , and 
 so may be somewhat interdependent . moreover , to create a perturbed 
 ligand with charge distribution closer to optimal , presumably one 
 would need to alter covalent chemistry , which could further change 
 molecular flexibility and the corresponding internal potential .
nevertheless , 
 an understanding of the charge optimum , in the presence of the current 
 covalent potential , could be quite valuable .
a feature of the 
 problem formulation adopted here is that a set 
 of candidate conformers ( including the preconformed ligand structure ) 
 is presented as input to the optimization procedure .
optimization 
 selects a conformer ( or a few degenerate conformers ) for the unbound 
 state and an optimized charge distribution that together minimize 
 the total binding free energy while satisfying all constraints . as 
 the set of unbound conformer candidates
the direct effect is that there are more 
 candidates from which to select the unbound conformation .
the indirect 
 effect is that if a new conformation is not selected as the unbound 
 ground state , its presence is an additional constraint that limits 
 the charge space available to the eventual ground state , because , 
 to be the ground state , it must be lower in unbound free energy than 
 all other conformers considered .
for example , if the previous ground 
 state is chosen as the ground state again , with the optimal charges 
 applied it has to be lower in free energy when unbound than all other 
 conformers including the new conformer added . for the test case used 
 here ,
the preconformed ligand structure was always selected as the 
 ( possibly degenerate ) ground unbound state by the optimization , even 
 when the number of candidates grew to 1001 . for increasing numbers 
 of candidates , the optimized charge distribution differed from the 
 rigid optimum due to the indirect effect , which diminished the optimum 
 binding affinity by over 0.5 kcal / mol for candidate sets with more 
 than 20 members . for larger sets than those studied here , or for other 
 cases , one imagines that a different conformation might be chosen 
 for the ground unbound state , but this would be a result of needing 
 to satisfy constraints rather than an adaptation to improve affinity . for the cases studied here , permitting flexibility
never improved 
 optimal affinity beyond that achieved by rigid binding optimization , 
 and , frequently , the induced - fit case incurred an affinity penalty , 
 relative to the rigid case .
one interpretation of these results is 
 that induced - fit binding , which may have a special role in molecular 
 recognition , may not be an optimization to enhance affinity .
this 
 interpretation is based on the theoretical framework developed here , 
 which includes partial atomic charge distributions that remain fixed 
 upon molecular conformational change .
it is formally possible that 
 conformational change is designed to create appropriate charge distribution 
 changes that do improve affinity , and we can not currently rule out 
 this possibility .
what we can say , however , is that conformational 
 change upon binding does not appear to enhance affinity by improving 
 the tradeoff between desolvation and interaction energetics , unless 
 it also can somehow simultaneously lead to affinity - enhancing charge 
 distribution changes .
these arguments can be generalized beyond 
 the specific molecular 
 example studied here by noting that induced - fit binding can be conceptualized 
 as a preconformation step , followed by a docking step , as illustrated 
 by the blue row of figure 1 .
constraining the 
 selected unbound conformer to the ground unbound state means that 
 the preconformation step has a free - energy change larger than or equal 
 to zero , with a lower bound of zero .
a lower bound on the docking 
 step is the rigid - ligand electrostatic optimum , which is also a lower 
 bound on gtotal. this lower bound will be an 
 optimum solution for the induced - fit case unless it violates the ground - state 
 constraint .
thus , the dominant - conformer optimum can not have higher 
 affinity than the rigid optimum , and any optimum producing an actual 
 ligand conformational change upon binding will have affinity either 
 equivalent to or weaker than the rigid optimization , and be a consequence 
 of having to satisfy the ground - state constraint . the practical 
 implication for molecular design of high - affinity 
 ligands
is that the construction of ligands preconformed for binding 
 and with rigid - ligand optimized partial atomic charges leads to highest 
 affinity in the theory developed here .
if the ligand does still change 
 conformation on binding , in many cases , the rigid - ligand optimized 
 partial atomic charges will still be excellent .
for the cases studied 
 here , the largest loss in electrostatic affinity from this strategy 
 would be 4.00 kcal / mol ( from 8.35 kcal / mol , which is calculated to 
 be the optimal value ) , and the largest loss in total affinity would 
 be 8.11 kcal / mol ( from 66.78 kcal / mol as the optimal value ) .
the current work adds molecular flexibility to charge optimization 
 theory but still treats each state as a single conformer .
molecular 
 systems actually exist as boltzmann - weighted distributions of conformations , 
 and the current theory is being extended to treat such ensemble distributions , 
 which could lead to new insights into the sources of conformational 
 changes that accompany binding .
 the design of ligands with high affinity and specificity 
 remains 
 a fundamental challenge in understanding molecular recognition and 
 developing therapeutic interventions . 
 charge optimization theory addresses 
 this problem by determining ligand charge distributions that produce 
 the most favorable electrostatic contribution to the binding free 
 energy . 
 the theory has been applied to the design of binding specificity 
 as well . 
 however , the formulations described only treat a rigid ligand  
 one 
 that does not change conformation upon binding . here , 
 we extend the 
 theory to treat induced - fit ligands for which the unbound ligand conformation 
 may differ from the bound conformation . 
 we develop a thermodynamic 
 pathway analysis for binding contributions relevant to the theory , 
 and we illustrate application of the theory using hiv-1 protease with 
 our previously designed and validated subnanomolar inhibitor . direct 
 application of rigid charge optimization approaches to nonrigid cases 
 leads to very favorable intramolecular electrostatic interactions 
 that are physically unreasonable , and analysis shows the ligand charge 
 distribution massively stabilizes the preconformed ( bound ) conformation 
 over the unbound . after analyzing this case , 
 we provide a treatment 
 for the induced - fit ligand charge optimization problem that produces 
 physically realistic results . 
 the key factor is introducing the constraint 
 that the free energy of the unbound ligand conformation be lower or 
 equal to that of the preconformed ligand structure , which corresponds 
 to the notion that the unbound structure is the ground unbound state . 
 
 results not only demonstrate the applicability of this methodology 
 to discovering optimized charge distributions in an induced - fit model , 
 but also provide some insights into the energetic consequences of 
 ligand conformational change on binding . specifically , the results 
 show that , from an electrostatic perspective , induced - fit binding 
 is not an adaptation designed to enhance binding affinity ; at best , 
 it can only achieve the same affinity as optimized rigid binding . 
Introduction
Theory
Methods
Results
Discussion and Conclusion