modeling the electrostatic potential and 
 energies in systems comprised 
 of biological macromolecules is an essential step for each study aimed 
 at understanding the macromolecules function , stability , and 
 interactions .
however , this is not a trivial task , especially for 
 huge systems made of large biomolecules and their assemblages and 
 for modeling phenomena occurring in microseconds and longer timeframes . 

continuum electrostatics offers an advantage over explicit methods 
 in such cases since ( a ) the atomic details of the water phase are 
 reduced and ( b ) continuum electrostatics intrinsically provides equilibrium 
 solutions . typically , the macromolecule is considered to be a low 
 dielectric medium while the water phase is modeled as a homogeneous 
 medium with a dielectric constant of 80 .
while there is a consensus 
 in the community that a dielectric constant of about 80 is appropriate 
 for describing dielectric properties of bulk water in modeling equilibrated 
 systems , the optimal value of the protein ( macromolecular ) dielectric 
 constant is still an ongoing debate in the literature .
this inconsistency is indicated by the use of numerous 
  optimal  dielectric constant values in various studies . 

investigations modeling the macromolecule as a rigid object or using 
 snap - shots obtained from molecular dynamics ( md ) simulations to deliver 
 the energies via molecular mechanics poisson  boltzmann ( mmpb ) 
 or generalized born ( mmgb ) methods typically use a low dielectric 
 constant of  = 1 or  = 2 ( to account for electronic 
 polarizability ) , although larger values were reported 
 as well . in works devoted to modeling 
 protein stability ,
numerous dielectric constant values were used , 
 from as low as  = 1 or 2 to as 
 high as  = 40 , including multidielectric 
 regions .
similarly , in the field of modeling 
 macromolecular interactions , researchers were using various values 
 for the protein internal dielectric constant .
perhaps the most widespread of the  optimal  values 
 for the dielectric constant is seen in the continuum methods for pka calculations .
the 
 most commonly used dielectric constant value is  = 4 , which 
 is believed to account for electronic polarization and small backbone 
 fluctuations .
however , larger values , 
 such as  = 8 ,  = 10 ,  = 11 , and 
  = 20 , were also reported .
other 
 examples can be listed as well , but it is clear that there is no universal 
 value of the dielectric constant that is appropriate for all models 
 and methods .
perhaps the largest body of work on the dielectric 
 constant of 
 proteins is due to warshel and co - workers .
it was demonstrated that the dielectric  constants  
 in semimacroscopic models depend on the definition and the specific 
 treatment . using various models and the 
 discriminative benchmark , warshel and co - workers demonstrated that 
 the protein dielectric constant is not a universal constant but simply 
 a parameter
that also depends on the model being used . in terms of structural reorganization occurring 
 during the process
being modeled , they pointed out that semimicroscopic 
 models which do not model the structural relaxations are forced to 
 use a large value for the dielectric constant .
in addition , including 
 the reorganization of ionized residues and plausible water penetration 
 would require applying an even higher dielectric constant in the modeling 
 algorithm . altogether , the works of warshel 
 and many other researchers point out that the 
 specific treatment of the conformational reorganization ( both of the 
 solute and the surrounding water phase ) in the continuum electrostatics 
 is the major determinant for the  optimal  value of 
 the internal dielectric constant .
modeling reactions involving large 
 conformational changes which are not explicitly treated in the computational 
 protocol requires the usage of large values for the dielectric constant , 
 while other reactions that do not induce conformational changes can 
 be successfully modeled with low dielectric constant even when a rigid 
 structure is used .
the goal is to develop an automatic 
 procedure which can assign an appropriate dielectric constant value 
 by utilizing the 3d structure of the macromolecule alone .
as 
 mentioned above , the continuum electrostatics approach is an 
 alternative to the explicit models , and it is expected to perform 
 the best if it can capture as many details as possible from the explicit 
 model .
explicit models treat the system as a multitude of atoms , some 
 of which are connected via chemical bonds , other interacting via van 
 der waals ( vdw ) and electrostatic nonbonded interactions ( figure 1a ) .
replacing the atomic details of the explicit 
 models with a continuum should address three important questions : 
 ( a ) how to treat the water molecules within the immediate water shell 
 surrounding the protein , ( b ) how to model the cavities inside the 
 protein , and ( c ) how to model the inhomogeneous protein matrix ( figure 1b ) .
many investigations demonstrated 
 the importance of water molecules 
 in the first layer of water shell and 
 that their properties differ from those of bulk water .
the main reason is that a significant fraction 
 of these water molecules may be involved in direct interactions with 
 protein atoms , and thus their ability to move and reorient can be 
 severely restricted in comparison with those of the bulk water .
in addition , the surface exposed amino acid 
 side chains are quite flexible , and if these alternative conformations 
 are ensemble averaged , the resulting shell surrounding the protein 
 will be a mixture of protein side chains and water molecules .
obviously 
 the dielectric properties of this shell can not be modeled with the 
 bulk dielectric constant of 80 .
in addition , the dielectric constant should depend 
 on the topology of the protein surface , such that the immediate water 
 shell around convex surfaces should be assigned a lower dielectric 
 constant ( water and side chain reorganization is more restricted ) 
 than that of the concave protein surfaces . because of this , it is also important to be able to determine the 
 artificial boundary between the protein and water , which results in 
 generating the so - called molecular surface
. the treatment of 
 cavities or channels inside the biomolecules is 
 another crucial problem for continuum electrostatics methods .
the water molecules in small cavities are very few , typically their 
 mobility is highly restricted and should be modeled with a dielectric 
 constant lower than that of bulk water .
in contrast , in large cavities , the waters may form clusters and 
 act collectively .
in addition , some cavities 
 may be filled with water molecules with short or long residential 
 times ; thus cavities or channels containing 
 water molecules with low occupancy should be treated as regions with 
 a low dielectric constant and vice versa .
perhaps the most crucial 
 issue for continuum electrostatics is 
 the treatment of the protein matrix .
atoms in biomolecules have different 
 charges and flexibilities depending on the packing and their mutual 
 interactions . because of that , macromolecules dielectric properties 
 are not homogeneous , resulting in a position - dependent dielectric 
 constant ( figure 1b ) .
different flexibilities 
 are accounted when comparing the binding site involved in  lock - and - key  

versus  induced - fit  binding as well as 
 those in allosteric regulations . the proteins 
 involved in electron and proton transfer are shown to be inhomogeneous 
 in their dielectric response to the charge translocation as well .
in general , the hydrophobic core is much better 
 packed and contains significantly fewer charged atoms than the molecular 
 surface and therefore has a much lower capacity to respond to the 
 local electrostatic field , and thus , from a continuum electrostatics 
 point of view , should be modeled as a low dielectric medium as indicated 
 in ref ( 35 ) .
these 
 observations prompted the development of approaches to assign a specific 
 local dielectric constant to each amino acid type or even to each 
 amino acid .
some of these approaches include the intrinsic polarizabilities 
 of amino acids delivered from the results of md simulations , using residue - specific dielectric constants 
 for pka calculations and optimizing the amino - acid - specific dielectric constant 
 to predict the folding free energy changes due to mutations .
all of these approaches treat the protein matrix 
 as a kind of mosaic structure made of small blocks with different 
 dielectric properties . bypassing the question about the applicability 
 of a macroscopic definition of the dielectric constant for such small 
 structural segments ,
the main problem with such approaches ( including 
 our previous work ) is the nonrealistic nature of the dielectric boundary 
 between the segments with different dielectric constants .
the 
 smooth dielectric function should reflect the physical nature 
 of the macromolecule , and the most straightforward approach is to 
 be delivered from atomic densities of the corresponding atoms .
space 
 occupied by loosely packed and charged atoms is expected to have the 
 potential to reorganize and to respond to the local electrostatics 
 field and thus should be modeled with a high dielectric constant . 

such regions correspond to a protein surface which is loosely packed 
 and rich in charged and polar residues .
in contrast , space regions 
 made of tightly packed uncharged atoms , as in the hydrophobic core , 
 have little ability to respond to the local electrostatic field and 
 should be modeled with a low dielectric constant .
all of these features 
 can be captured with the so termed gaussian method of representing 
 atomic densities and then using it to deliver a smooth dielectric 
 function .
it should be clarified that we emphasize 
 using such an approach to deliver a smooth dielectric function for 
 the entire space but not to model the water  molecule boundary or derive the surface of molecules . in this work ,
we report a gaussian - based approach to deliver 
 a 
 smooth dielectric function for the entire space domain ( the macromolecule(s ) 
 and the water phase ) . to assess the effects of this development on 
 various biophysical quantities which can be computationally modeled
, 
 the method was implemented in the delphi program , and several tests were carried out .
the correctness of 
 the approach in reflecting the expected physical properties of macromolecules 
 is demonstrated on a large set of protein structures by showing that 
 the generated smooth dielectric function results in values of about 
 67 in the protein interior and values of 2030 at the 
 protein  water interface , which is consistent with previous 
 md - based work .
in addition , it is shown 
 that the space occupied by charged and polar groups is assigned larger 
 dielectric values as compared to the space occupied by hydrophobic 
 amino acids .
furthermore , in a benchmarking test of a large set of 
 experimentally determined pka s 
 in staphylococcal nuclease ( snase ) , the gaussian - based smooth dielectric 
 function ( with low reference dielectric constant value of 4 ) delivers 
 much better results ( rmsd = 1.77 ) as compared with the homogeneous 
 dielectric method ( rmsd = 2.43 ) , despite the fact that the homogeneous 
 dielectric method optimal results were obtained at a higher dielectric 
 constant value of 10 .
the results of homogeneous and gaussian - based 
 smooth dielectric methods for each individual titratable group are 
 shown in supporting information ( table 
 s1 ) .
given a molecule 
 in the water phase , we applied the gaussian equation and implemented 
 three steps as follows to calculate the dielectric distribution of 
 a protein from its density distribution as originally described by 
 nicholls et al .
given a macromolecule 
 with n atoms , the density of an atom i is represented by a gaussian distribution ( figure 2a , b):1where i(r ) is 
 the density at position r , ri is the distance 
 between the center of the atom i and position r , ri is the 
 vdw radius of atom i , and  is the relative 
 variance .
upper panels : cartoon presentation in the 2d grid plane , 
 where intensity of the colors reflects the value of either the density 
 or . lower panels : profiles of the density and resulting . after the density of each atom 
 within the macromolecule is generated , 
 the density in the overlapping areas occupied by multiple atoms is 
 calculated by2where the mol(r ) denotes the density at position r coming 
 from 
 multiple atoms , and i(r ) is the density of a single atom i , which is obtained 
 from eq 1 .
this function guarantees that the 
 density at the overlapping region is higher than the density generated 
 by any involved single atom , but the density mol ( r ) always stays between 0 and 1 ( figure 2c , d ) .
finally , the smooth dielectric function is 
 delivered from the density distribution using the linear function:3where  on the left denotes the dielectric 
 distribution function , in denotes the reference 
 dielectric value when the density is 1 , out denotes 
 the reference dielectric value for the water phase , and  is 
 the density obtained from eq 2 ( figure 2e , f ) . the above formalism ,
provided that standard 
 force field atomic 
 vdw radii are used , has only one adjustable parameter , the variance 
 in the gaussian distribution .
the functional forms of eqs 13 can also be considered 
 as  adjustable  since atomic density in principle can 
 be modeled with almost any symmetrical smooth function ( different 
 from gaussian ) , and the corresponding molecular density and dielectric 
 function can be delivered by other means  different from eqs 2 and 3 .
however , it seems to 
 us that the most important parameter is still the variance in the 
 gaussian distribution , since its variation will essentially mimic 
 the usage of different functions in eq 1 . in 
 the current work ,
the functional forms of eqs 2 and 3 are kept as originally suggested in 
 ref ( 40 ) , but alternatives 
 will be explored later .
thus , in order to proceed further with the 
 analysis , the value of the normalized variance in eq 1 must be selected . in this work
, we chose to select it based 
 on a benchmarking test of pka values , 
 because dielectric relaxation is the most profoundly associated with 
 ionization / deionization phenomena ( see discussion for more details ) . in order 
 to assess the impact of the smooth dielectric function on the accuracy 
 of pka calculations and to deliver the 
 optimal value of the normalized variance
, we utilized the experimental 
 data obtained in garcia - moreno s lab ( http://pkacoop.org/wordpress/?p=28 ) and used in the pka cooperative ( http://pkacoop.org/wordpress/ ) .
the reason for selecting this 
 data set is not only because it is a representative set for the pka cooperative but also because it involves pka s both surface - exposed and buried amino 
 acids .
this data set is comprised of 89 pka s for staphylococcal nuclease ( snase ) . in 19 of the cases , 

the pka calculations were based on the 
 structure of wild type ( pdb i d : 1stn ) snase and 
 its hyperstable variant ( pdb i d : 3bdc ) , which 
 is called +phs . in 20 of the cases ,
pka calculations were based on the x - ray structures of snase 
 with mutations , and in 50 of the cases , the pka calculations were based on the in silico generated mutant 
 from the wild type snase ( pdb i d : 3bdc ; list of the amino acids , structure used , 
 and experimental pka s are available 
 from the pka cooperative web page , http://pkacoop.org/wordpress/ ) . since the predictions 
 of the pka s in our list do not 
 require modeling multiple titration sites , but the prediction of a 
 single pka per structure , the following 
 surface - free approach ( sfa ) was developed and applied . note that because 
 energies are delivered as grid energies of the corresponding finite - difference 
 algorithm , there is no need to define a molecular surface .
this provides 
 a significant advantage over previous works , since defining the molecular 
 surface would add additional uncertainty in the protocol .
in addition , 
 the sfa reflects best the motivation and the development of the smooth 
 dielectric function : no need to draw a sharp border between the protein 
 and the water phase . for each pka calculation , 
 four structures were generated ( figure 3 ) : 
 ( 1 ) the deprotonated state of the concerned residue in the protein , 
 ( 2 ) the protonated state of the concerned residue in the protein , 
 ( 3 ) the deprotonated state of the concerned residue in water , and 
 ( 4 ) the protonated state of the concerned residue in water .
the structures 
 of the protein and the residue for which pka is calculated were kept identical in the four states in order to 
 cancel the artificial 
the extra proton of protonated asp / glu was not modeled , but its charge 
 was distributed evenly over the carbonyl oxygens ( the argument being 
 to avoid artificial grid energy ) .
similarly , the missing proton of 
 deprotonated lys / arg was not removed from the side - chains , but the 
 charge was distributed over all polar hydrogens to result in zero 
 net charge for deprotonated lys / arg .
pka of a residue 
 calculated via four 
 step procedures representing four states : ( a ) the deprotonated state 
 of the residue in the protein , ( b ) the protonated state of the residue 
 in the protein , ( c ) the deprotonated state of the residue in water , 
 and ( d ) the protonated state of the residue in water . for each state ,
the grid energy was calculated , 
 keeping the grid 
 and position of the protein and the residue of interest identical 
 among the runs , and the following energies were obtained : ( 1 ) g ( depro , protein ) , the 
 protein electrostatic grid energy with the deprotonated residue ; ( 2 ) g ( proto , protein ) , the 
 protein electrostatic grid energy with the protonated residue ; ( 3 ) g ( depro , water ) , the 
 electrostatic grid energy of the deprotonated residue in the water 
 phase ; and ( 4 ) g ( proto , water ) , the electrostatic grid energy of the protonated 
 residue in the water phase
. then , the value of the pka shift was calculated as4 the parameters used in delphi were 
 scale = 2 points /  , perfil = 70 , and out = 
 80.0 .
the convergence criterion was 0.0001 [ kt / e ] , and 
 the linear poisson  boltzmann equation was solved .
the internal 
 reference dielectric constant and the normalized variance were considered 
 to be adjustable parameters .
the charges and radii were taken from 
 the amber force field . in order to ensure 
 the statistical significance of the analysis and to assess the general 
 trend of the smooth dielectric function ,
a large set of diverse proteins was taken from the pdb bank 
 ( http://www.rcsb.org/pdb/home/home.do ) , and several filtering 
 steps were performed .
first , only structures determined by x - ray experiments 
 with a resolution less than 1.5  were selected .
finally , 
 the structures with cofactors which are not made of regular residues 
 were also deleted from the data set .
the final data set was made of 
 91 proteins : 1arb , 1bgf , 1bkr , 1byi , 1c8c , 1es9 , 1ew4 , 1ezg , 1g8a , 1gpp , 1gqi , 1h1n , 1hzt , 1i1j , 1i2 t , 1idp , 1ijy , 1jl1 , 1k0 m , 1kmt , 1kng , 1l3k , 1lkk , 1llf , 1lu4 , 1lzl , 1m1f , 1mf7 , 1mn8 , 1my7 , 1n62 , 1ng6 , 1nkd , 1nko , 1o7i , 1o7j , 1o82 , 1o9i , 1oai , 1oi7 , 1r29 , 1r7j , 1sg4 , 1sjy , 1smo , 1t3y , 1tqg , 1tvn , 1ucs , 1ukf , 1ulr , 1usm , 1utg , 1uxz , 1uz3 , 1v05 , 1vh5 , 1w5r , 1w7c , 1x6z , 1xod , 1xsz , 1z21 , 1zzk , 2a26 , 2a6z , 2b0a , 2car , 2e3h , 2end , 2fcj , 2fhz , 2fq3 , 2gec , 2gom , 2h3l , 2hox , 2hqx , 2ivy , 2j5y , 2j6b , 2ohw , 2p5k , 2pmr , 2pth , 2ve8 , 2vn6 , 3bbb , 3bzp , 3cjs , 3ct6 .
the data 
 set was used to analyze the following plausible relations and distributions : 
 ( 1 ) the average dielectric constant (  ) and radius of gyration 
 of the protein and ( 2 ) the average  and residue type distribution .
( 1 ) average dielectric constant distribution : for each protein 
 in the data set , the center of the mass and gyration radius were calculated . 

then , the protein interior was mapped into different shells with a 
 different radius from the center of the protein . for each shell
, the 
 average  value obtained from all grid - midpoints inside this 
 shell was calculated .
thus , an -radius map was generated for 
 each of the proteins ( figure 4a ) .
( 2 ) 
 -residue type distribution : the average dielectric constant 
 per residue was calculated using only the side chain atoms ( backbone 
 was not included ) .
then , a sphere of radius 5  was drown around 
 each side - chain atom , and the dielectric constants of all midgrid 
 points within the sphere were summed and averaged .
further , these 
 average dielectric constants were summed over all atoms of the side - chain 
 and averaged again to obtain the average dielectric constant per side 
 chain .
finally , the average dielectric constant for each type of residue 
 was obtained from all residues with the 91 protein set .
average dielectric 
 constant against the radius of gyration ( rg ) 
 of a protein .
( a ) schematic of calculating the average  value 
 of a shell with radius r and thickness d ; the thickness of the shell is set as 0.1 rg of a protein .
( b ) average 
 -radius distribution of 91 proteins in the data set . left upper 
 corner shows the structure of elongated protein ( 1uz3.pdb ) for which 
 the spherical shape is not a good approximation , and on the right 
 lower corner is shown a protein ( 1brk.pdb ) which is almost spherical 
 in shape .
the 
 test on small molecule free energies of transfer from a vacuum to 
 water was done on a data set of 504 neutral organic small molecules taken from david mobley s group .
the 
 solvation energies of all of these molecules have been experimentally 
 determined , with the range from 11.95 to 3.16 kcal / mol .
solvation energy gsol has two components , 
 polar and nonpolar:5where gpolar indicates 
 the polar ( electrostatic ) term and gnonpolar denotes the nonpolar term .
the polar component of solvation 
 energy was calculated as the grid 
 energy difference of the system in water and in a vacuum:6 the above grid energies were calculated keeping the corresponding 
 small molecule at the same grid position to cancel the grid artifacts . 

specific considerations were made for the calculations in a vacuum 
 since one has to define the molecular surface in this case ( the border 
 between molecule and vacuum ) .
note that in our approach the dielectric 
 function is continuous and runs throughout the entire space and is 
 designed to describe dielectric properties of the molecule in water . 
 here , we assume that the properties of molecules are unchanged as 
 they are moved from water to a vacuum . thus ,
following the strategy 
 implemented in zap , the molecular surface 
 of molecules is defined by applying a specific cutoff for the dielectric 
 constant , cutoff .
the cutoff was varied in the protocol 
 to obtain the best fit against experimental data .
the nonpolar 
 term of solvation energy gnonpolar is 
 calculated via the accessible surface area method:7where  
 and b are constants 
 and sa denotes the solvent accessible surface area , which is calculated 
 using naccess2.1.1 ( http://www.bioinf.manchester.ac.uk/naccess/ ) .
the force field used in the calculations was am1-bcc , which is part of general amber force field ( gaff ) . in order to optimize the parameters ,
reference 
 in was varied from 0.1 to 4.0 , the value of normalized 
 variance i was varied from 0.80 
 to 1.40 , and the value of epsilon cutoff bnd was 
 varied from 4.0 to 60.0 . for each combination of the i and bnd values ,
the least - squares 
 method was used to obtain the optimized  and b constants .
the best parameters are shown in the results section ( note that because of the different nature 
 of the process , these values are not expected to be the same as those 
 obtained in pka calculations .
in this 
 section , the results of parametrizing and using the smooth 
 gaussian - based dielectric function are presented .
it should be restated 
 that this development is aimed to better capture the effects seen 
 in explicit electrostatic calculations within the framework of continuum 
 electrostatics for proteins .
the proteins are the primary target because 
 of their intrinsic conformational flexibility , the presence of small 
 or large cavities inside , and an irregularly shaped protein  water 
 interface . however , introducing a smooth gaussian - based dielectric 
 function should not deteriorate the performance of continuum electrostatics 
 even in the cases involving modeling molecules rigid by nature with 
 no cavities inside , such as small molecules and drugs .
thus , the features 
 of the smooth gaussian - based dielectric function are manifested modeling 
 several classes of problems : ( a ) pka s 
 of various titratable groups with the pka cooperative initiative , ( b ) conceptual considerations and linkage 
 to the physical properties of proteins and amino acids , ( c ) the role 
 of electrostatic potential on the electron transfer in the reaction 
 center protein , ( d ) free energy of transfer of small molecules , and 
 ( e ) reducing grid dependency in energy calculations .
analysis of the results 
 at the last two pka cooperative meetings showed that 
 groups using rigid structures for pka calculations 
 obtain the best results with respect to the experimental data if a 
 large value is assigned for the protein dielectric constant .
the works 
 of groups utilizing md techniques to predict pka s indicated that almost in each case the ionization 
 induces small or large conformational change in the proteins .
because of these observations , the pka cooperative data set was chosen for testing the gaussian - based smooth 
 dielectric function .
thus , the pka values 
 of 89 residues in staphylococcal nuclease protein ( snase ) were calculated 
 using the sfa method described in the methods section , and predictions were compared between homogeneous dielectric 
 and gaussian - based methods . in both cases ,
the value of the internal 
 dielectric constant was independently varied to obtain the smallest 
 rmsd with respect to experimental data . in order to determine the 
 best protein value for the pka calculation in a homogeneous protein dielectric ,
the protein value was varied from 1.0 to 20.0 with an increment 
 of 1.0 . for each protein value
, the pka values were calculated for all 89 cases and then compared 
 to the experimental results .
it was found that the best protein value for homogeneous dielectric pka calculations is 10.0 , which resulted in a rmsd between 
 calculated and experimental pka values 
 of 2.43 pk ( figure 5 ) . using 
 a gaussian - based method on the same data set ,
two parameters need 
 to be determined , the values of the reference dielectric value in and the normalized variance . here , we varied the 
 reference in from 1.0 to 10.0 with an increment 
 of 1.0 and  from 0.80 to 1.20 with an increment of 0.01 . for 
 each combination of these two parameters , the calculated pka values for all 89 cases were compared with 
 the experimental values . from this test , the best parameters for pka calculations obtained are in= 4.0 and  = 0.93 , resulting in rmsd = 1.77 ( figure 5 ) .
results of pka calculations 
 of 89 residues 
 in staphylococcal nuclease protein ( snase ) , calculated using the homogeneous 
 method and gaussian method in delphi .
the fact that the gaussian - based method achieves a much smaller 
 rmsd than the homogeneous protein dielectric method is very encouraging . 

in addition , the best results using the gaussian - based method were 
 obtained at a low reference dielectric constant of 4 , while the best 
 dielectric constant using the homogeneous method was 10 .
this indicates 
 that the gaussian - based method mimics the effects of the conformational 
 changes occurring in the titration better than the homogeneous high 
 dielectric does by distributing dielectric values within the protein 
 structure .
the first task is to check if the generated smooth 
 gaussian - based dielectric function addresses the questions mentioned 
 in the introduction about modeling the protein 
water 
 interface , cavities inside the protein , and the protein matrix . for 
 the testing ,
the optimal values for the reference internal dielectric 
 constant and normalized variance obtained in the pka section were used .
the reaction center protein ( pdb 
 i d : 1aij ) was 
 taken as a test case ( although the test was done on many other proteins 
 as well ) .
figure 6 shows the results of ( a ) 
 a slice of dielectric distribution and the entire protein structure 
 and ( b ) a slice of dielectric distribution with atoms close to this 
 slice .
several differences from the traditional homogeneous dielectric 
 distribution method can be seen , the smooth gaussian - based method 
 results in three improvements : ( 1 ) the dielectric distribution near 
 the protein  water boundary region no longer contains a discrete 
 dielectric jump .
instead , it now smoothly increases from a low dielectric 
 value inside the protein to water = 80 in bulk water . 

( 2 ) the small cavities inside the protein are assigned proper dielectric 
 values . for small cavities , which may contain only a limited number 
 of water molecules , which are not free to move and reorient , the dielectric 
 value for such cavities is neither water nor protein . instead , the gaussian - based method assigns dielectric 
 values between water and protein depending on the cavity size and shape .
( 3 ) the dielectric distribution 
 inside the entire protein is not simply a constant value ( protein ) .
the regions with tightly packed atoms are assigned lower 
 dielectric values , while the regions with loosely packed atoms get 
 higher dielectric values .
this distribution reflects the physical 
 nature of the dielectric response : the regions filled with loosely 
 packed atoms should be treated as polarizable media , since the atoms 
 in those areas can move to respond to the local electrostatic field . 

in contrast , tightly packed atoms in the hydrophobic core do not carry 
 much charge and can not respond to the local field and therefore should 
 be considered as regions with low polarizability .
( b ) a slice 
 of dielectric distribution with atoms close to the slice surface .
it is expected that the average dielectric constant in the core of 
 the protein is lower than that on the surface , as indicated in a previous 
 md study . here ,
the same question is 
 addressed on a set of 91 diverse proteins ( see the methods section ) .
each black line in figure 4b denotes the dielectric distribution of a particular protein . 

it can be seen that the distributions of the dielectric constant of 
 all the proteins in the data set have the same tendency : the dielectric 
 value in the core is close to protein .
as the shell 
 moves toward the surface of the protein , the dielectric value grows . 

when the radius of the shell is about two gyration radii , the average 
 dielectric value of the shell region reaches water .
however , despite the overall similarity , different proteins have 
 different dielectric distributions . when the radius of the shell is 
 equal to the gyration radius , the maximum and minimum average dielectric 
 values of the shell are 64.7 for protein pdb i d 1uz3 and 8.4 for the protein pdb i d 1bkr .
the reason for this difference is a 
 different shape of these two proteins ( see figure 4b ) .
the protein 1uz3 is a very elongated , nonspherically shaped protein , 
 and the average dielectric value is large at the shell with a radius 
 equal to the gyration radius because it includes water as well .
if 
 the protein has a spherical - like shape , the average dielectric value 
 of the shell with a radius equal to the gyration radius is smaller , 
 as for the protein 1uz3 .
the red line in figure 4b indicates the 
 average behavior of the dielectric distribution of all 91 proteins . 

this general tendency reveals that the inner parts of the proteins 
 have lower average dielectric values and the outer parts have higher 
 values .
another analysis was carried out to show the distribution of dielectric 
 values per residue type . for each residue , the average dielectric 
 value of all side chains was calculated as described in the methods section . the calculations were performed 
 on all 91 proteins in the data set , and then the results were averaged 
 per amino acid type ( figure 7 ) .
figure 7 shows the average dielectric constant value for 
 each type of residue , and it can be seen that the range is from 11.0 
 to 25.6 . charged amino acids ( lys , arg , glu , and asp ) are associated 
 with the highest average dielectric values .
they have the propensity 
 to be located on the surface of the protein and to be loosely packed , 
 leaving room for structural rearrangement .
the observation that the 
 gaussian - based dielectric function has the largest value assigned 
 for the space occupied by such residues correctly reflects the physics 
 of dielectric relaxation and supports the usefulness of the developed 
 procedure .
thus assigning high dielectric values for charged residues 
 is physically sound , which has been shown to improve the accuracy 
 of pka s and the electrostatic potential calculation . 

the space occupied by hydrophobic residues is assigned relatively 
 low average dielectric value ( figure 7 ) .
hydrophobic 
 residues are more likely to be found in the core of the protein , typically 
 tightly packed and not able to adopt alternative conformations . in 
 addition , their side chains are made of atoms carrying little charge . 
 because of that , their ability to alter the local electrostatic field 
 is very limited .
the fact that the gaussian - based method automatically 
 assigns a low dielectric constant for the hydrophobic residues is 
 very encouraging and demonstrates that the method captures the correct 
 physics . in the middle
their average dielectric 
 values ( figure 7 ) are higher than that of hydrophobic 
 residues but lower than that of charged residues .
thus , the model 
 reflects the correct physics for these amino acids as well .
the analysis is focused on 
 investigating the role of the electrostatic potential on the electron 
 transfer from quinone a ( qa ) to quinone b ( qb ) in the reaction center 
 protein .
previous works demonstrated 
 that this process is slow and involves conformational changes in the 
 protein .
since the goal of the gaussian - based smooth dielectric method 
 is to mimic the effects of conformational changes via properly assigned 
 dielectric constant , this particular reaction is a perfect test for 
 the method .
to illustrate the advantages of the gaussian - based smooth 
 dielectric method over the standard homogeneous dielectric approach , 
 the electrostatic potential distribution in the reaction center protein 
 ( pdb i d : 1aij ) was calculated using both the regular homogeneous method and the 
 gaussian - based method implemented in delphi .
the electron transfer is from qa to qb , the electrostatic 
 potential is expected to be less negative ( more positive ) at the electron 
 acceptor site ( qb ) as compared with the electron donor site ( qa ) . 

previous work and figure 8 demonstrate that the potential delivered with a homogeneous 
 protein dielectric constant opposes the electron transfer ; i.e. , the 
 potential at qb is more negative than at qa site . in our previous 
 investigation
, this problem was solved 
 by utilizing the results of ref ( 20 ) , showing that the vicinity of the qb site is 
 more flexible than that of the qa site , which made us assign a high 
 dielectric constant of 20 for the residues with a qb pocket while 
 modeling the rest of the protein with a dielectric constant of 4 . 

such an approach requires both prior knowledge about the flexibility 
 of different protein regions and manual assignment of an appropriate 
 dielectric constant .
in contrast , the gaussian - based smooth dielectric 
 method automatically assigns adequate dielectric constants based on 
 the 3d structure of the protein and , as seen in figure 8 , greatly reduces the potential barrier for the electron transfer . 

compared to the potential in figure 8a , the 
 potential at the qa site in figure 8b is less 
 positive and the qb site is less negative .
potentials of qa and qb 
 in figure 8b are at a similar level , which 
 agrees with the experimental observation of the electron transfer 
 process .
electrostatic potential in the reaction center protein ( pdb i d : 1aij ) , calculated by 
 ( a ) the homogeneous method in delphi and ( b ) the gaussian - based method 
 in delphi .
previous work on the same data 
 set of 504 small molecules utilized md 
 simulations with the explicit solvent 
 model and reported rmsd of 1.24 kcal / mol . here
, we use the same data set of small molecules to test the performance 
 of both the homogeneous and gaussian - based dielectric .
it should be 
 clarified that small molecules are much less flexible than proteins 
 and do not have internal cavities .
because of that , they are not expected 
 to be very polarizable and the gaussian - based method is not expected 
 to offer any advantage over the homogeneous dielectric . however ,
if 
 the development is correct , it is anticipated that the results obtained 
 with the gaussian - based and homogeneous dielectrics should be quite 
 similar . keeping in mind that small molecules are less polarizable 
 than proteins
, it is expected that different ( as compared with those 
 in proteins ) reference values for the internal dielectric constant 
 and variance will deliver the best rmsd . here , we optimized the reference 
 internal dielectric constant , the variance of the gaussian function , 
 the cutoff for the molecule  water interface , and the coefficients 
 of the nonpolar component of the solvation energy ( calculated with 
 the solvent assessable surface area , sasa ) .
the best results are achieved 
 when the parameters are set as follows : the reference in = 1.0 , reference variance i = 1.0 , cutoff for the molecular boundary is at bnd = 16 . for the nonpolar term calculation ,  = 0.0028 
 kcal/(mol  ) and b = 0.0948 kcal / mol . 

the corresponding rmsd is 1.59 kcal / mol , which is close to the results 
 obtained with explicit water simulations .
the results of the calculations 
 are shown in figure 9 against the experimentally 
 measured transfer energies .
for homogeneous dielectric calculations , 
 the best rmsd = 1.45 kcal / mol ( figure 9 ) , when 
 the parameters are set as follows : in = 1.0 ,  
 = 0.0094 kcal/(mol  ) , and b = 1.1579 
 kcal / mol .
solvation energies of small molecules calculated by gaussian - based 
 dielectric methods . calculating electrostatic energies by 
 using finite difference methods
if the concerned biomolecule is shifted 
 to a new position , the newly obtained solvation energy might be different 
 from those obtained at the original position .
since the gaussian - based 
 method uses smoothed dielectric distribution rather than two - dielectric 
 constants with a sharp jump at the protein 
water boundary , 
 it is expected that the energies calculated by the smoothed method 
 should be less grid - sensitive than traditional methods .
figure 10 shows the grid sensitivity of the electrostatic 
 solvation energy of a small molecule calculated with both a gaussian - based 
 smooth dielectric and with a homogeneous dielectric .
the first test involved translation 
 such that the small molecule is moved along the x direction in steps 
 of 0.005  and the energy calculated for each position ( figure 10a ) .
the molecule was rotated in steps of 10 , and for each rotation 
 the energy was calculated ( figure 10b ) .
figure 10 shows that energies calculated with the gaussian - based 
 smooth dielectric function are much less sensitive to the grid details , 
 which is another significant advantage of the proposed approach . grid 
 sensitivity of the gaussian method and homodeneous dielectric 
 delphi .
( b ) 
 a triethylamine molecule is rotated in 360. ( c ) energy sensitivity 
 of shift .
since the nature and the  optimal  
 value of the dielectric 
 constant of proteins ( and macromolecules in general ) was and is the 
 subject of many investigations and scientific publications , it is 
 worth summarizing the outcome of this work with regard to frequently 
 carried research tasks and scenarios .
one of the most common tasks 
 in computational research is to calculate the electrostatic component 
 of the solvation energy of a protein or a small molecule .
as mentioned 
 in the introduction , such calculations were 
 reported in the literature utilizing various values for the internal 
 dielectric constant .
which dielectric value is the best ? reiterating 
 once more the lessons outlined by warshel and co - workers , we would like to clarify the difference between calculations involving 
 small molecules and those involving proteins .
it is useful to envision 
 the case of a molecule for which atomic details are revealed with 
 absolute precision ( including hydrogens ) .
further idealizing the model 
 considering that the thermal motions are not present , the electrostatic 
 interactions among the atoms of the molecule should be calculated 
 with a dielectric constant of 1 ( or 2 to account for electronic polarizability ) . 

this idealized case represents , to a great extent , the real case scenario 
 of modeling the solvation energy of rigid molecules and , in general , 
 modeling processes which do not induce conformational ( and ionization ) 
 changes .
the above discussion does not represent the best test 
 case of our 
 approach , since the dielectric response is expected to be manifested 
 the best upon introducing / removal of a charge .
the natural process 
 in proteins involving charge creation / removal is ionization / deionization 
 of titratable groups upon ph changes or other factors . because of 

that , in our opinion , the most appropriate phenomenon to illustrate 
 the advantages of the smooth dielectric function is the pka s modeling .
calculating the pka shift requires one to find the energy difference ( most 
 electrostatic component of the energy ) between the corresponding protein 
 with ionized and neutral groups of interest .
however , typically only 
 the structure of the protein with the group being ionized or neutral 
 is available .
depending on various structural factors , the ionization / deionization 
 of the group could induce small or large conformational changes , which 
 are not taken into account in the protocols using rigid structures . 
 to account for these structural changes , which in turn will affect 
 the electrostatic interactions ,
we would like to point out that despite the use of a unified reference 
 internal dielectric constant (  = 4 ) for our protein modeling , 
 the corresponding smooth dielectric function is unique for each protein . 
 this is because it results from structural characteristics such as 
 packing , the presence of internal cavities , and their shapes .
because 
 of that , it is expected that the proposed approach and delivered parameters 
 can be used without adjustment in future studies involving different 
 sets of proteins .
the analysis of the smooth dielectric function 
 done on a large 
 and diverse set of proteins proved that the approach is physically 
 sound .
indeed , without explicit considerations about physical properties 
 of the amino acid ( charged , polar versus hydrophobic ) , the results 
 on the average dielectric constant value over amino acid types showed 
 that the highest dielectric values were automatically assigned to 
 the space occupied by charged groups , and the lowest dielectric values 
 were attributed to the hydrophobic core .
this is obviously related 
 to the packing and the fact that charged and polar residues prefer 
 to be at the protein surface , while hydrophobic groups are typically 
 buried .
however , engineered proteins may not follow such a trend and 
 may have a titratable group buried in the packed hydrophobic core , 
 as with many of the cases in the pka cooperative 
 data set .
a buried titratable group , charged or not , represents a 
 highly polar structural element , for which the contribution to the 
 protein s polarizability will depend on the ability to change 
 conformation or ionization states and may not be captured by our model . 
 to address
the possibility that the local dielectric constant associated 
 with the space taken by titratable groups ( especially in case of non - naturally 
 occurring groups ) may have to be not only a function of atomic packing 
 but to be further increased , we repeated the pka calculations with a modified gaussian - based approach : the 
 radii of titratable atoms of asp , glu , arg , and lys were artificially 
 lowered by a factor of 2 , which resulted in reduced packing and a 
 higher local dielectric constant .
the best rmsd ( rmsd = 1.75 pk ) was found to be almost the same as it was with the original 
 gaussian - based method ( rmsd = 1.77 pk ) .
this observation 
 indicates that the current model captures most of the effects relevant 
 to the dielectric response upon ionization , while it does not rule 
 out further improvement by specific treatment of the ionizable groups .
 implicit methods for modeling protein 
 electrostatics require dielectric 
 properties of the system to be known , in particular , the value of 
 the dielectric constant of protein . 
 while numerous values of the internal 
 protein dielectric constant were reported in the literature , still 
 there is no consensus of what the optimal value is . 
 perhaps this is 
 due to the fact that the protein dielectric constant is not a  constant  
 but is a complex function reflecting the properties of the protein s 
 structure and sequence . here , we report an implementation of a gaussian - based 
 approach to deliver the dielectric constant distribution throughout 
 the protein and surrounding water phase by utilizing the 3d structure 
 of the corresponding macromolecule . 
 in contrast to previous reports , 
 we construct a smooth dielectric function throughout the space of 
 the system to be modeled rather than just constructing a  gaussian 
 surface  or smoothing molecule  water boundary . 
 analysis 
 on a large set of proteins shows that ( a ) the average dielectric constant 
 inside the protein is relatively low , about 67 , and reaches 
 a value of about 2030 at the protein s surface , and 
 ( b ) high average local dielectric constant values are associated with 
 charged residues while low dielectric constant values are automatically 
 assigned to the regions occupied by hydrophobic residues . in terms 
 of energetics , 
 a benchmarking test was carried out against the experimental 
 pka s of 89 residues in staphylococcal 
 nuclease ( snase ) and showed that it results in a much better rmsd 
 (= 1.77 pk ) than the corresponding calculations done 
 with a homogeneous high dielectric constant with an optimal value 
 of 10 ( rmsd = 2.43 pk ) . 