September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 FINITE ELEMENT ANALYSIS OF DRUG ELECTROSTATIC DIFFUSION: INHIBITION RATE STUDIES IN N1 NEURAMINIDASE YUHUI CHENG, MICHAEL J. HOLST AND J. A. MCCAMMON Uinversity of California, San Diego 9500 Gilman Dr. MC 0365 La Jolla, CA 92037, USA E-mail: ycheng@mccammon.ucsd.edu This article describes the numerical solution of the steady-state Poisson-BoltzmannSmoluchowski (PBS) and Poisson-Nernst-Planck (PNP) equations to study diffusion in biomolecular systems. Specifically, finite element methods have been developed to calculate electrostatic interactions and ligand binding rate constants for large biomolecules. The resulting software has been validated and applied to the wild-type and several mutated avian influenza neurominidase crystal structures. The calculated rates show very good agreement with recent experimental studies. Furthermore, these finite element methods require significantly fewer computational resources than existing particle-based Brownian dynamics methods and are robust for complicated geometries. The key finding of biological importance is that the electrostatic steering plays the important role in the drug binding process of the neurominidase. 1. Introduction Electrostatics and Diffusion play the important role in a variety of biomolecular processes, which have been studied extensively using various biophysical, biochemical and computational methods. Computational models of diffusion have been widely studied using both discrete 1,2,3,4,5 and continuum methods 6,7,8,9,10,11 . The discrete methods concentrate on the stochastic processes based on individual particles, which include Monte Carlo12,13,14,5, Brownian dynamics (BD) 15,16,17 and Langevin dynamics18,19 simulations. Continuum modeling describes the diffusional processes via concentration distribution probability in lieu of stochastic dynamics of individual particles. Comparing with the discrete methods, continuum approaches needn't deal with the individual Brownian particles This work is supported by nih gm31749, nsf mcb-0506593 and mca93s013 (to jam). additional support from the howard hughes medical institute, the nsf supercomputer centers, the san diego supercomputing center, accelrys, inc., the w.m. keck foundation, the national biomedical computational resource and the center for theoretical biological physics is gratefully acknowledged. 1 September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 2 and the computionational cost can be substantially less than for the discrete methods. In the present work, we apply adaptive finite element methods to solve the PBS equation and PNP equation in several neurominidase structures. The diffusion results have been compared with those from recent experimental kinetic studies. The H5N1 avian influenza neurominidase is a highly pathogenic virus that might acquire the ability to pass readily among humans and cause a pandemic 20,21. Our continuum modeling demonstrates that it is efficient but accurate enough to address the drug binding. To take account of the stochastic dynamics of diffusion particles, we are developing a state-in-art SMOL package by integrating the continuum modeling and Brownian dynamics together to study the drug diffusion in the neurominidase. 2. Theory and Modeling Details Our SMOL package (http://mccammon.ucsd.edu/smol/index.html) models the diffusion of ligands relative to a target molecule, subject to a potential obtained by solving the Poisson-Boltzmann equation. It is perhaps most easily explained by initially considering motion of an ensemble of Brownian particles in a prescribed external potential (R)(R being a particle's position) under conditions of high friction, where the Smoluchowski equation applies. 2.1. Boundaries and Initialization of the Poisson-Boltzmann and Smoluchowski Equations The original Poisson-Boltzmann and Smoluchowski equations have the form of continuity equations (see Fig. 1): · ((R)) + 2sinh((R)) = qi (R - Ri ) + q j p j (R; t ) i j (1) (2) p j (R; t ) = - · J j (R; t ) t where the particle flux is defined as: J j (R; t ) = D j (R)[ p j (R; t ) + (R) p j (R; t )] = D j (R)e-q j ( ) eq j ( ) p(R; t ) R R (3) Here p j (R; t ) is the distribution function of the ensemble of Brownian particles, including the drug molecule, coions and counterions. D j (R) is the diffusion coefficient of the diffusion particle, = 1/kB T is the inverse Boltzmann energy, kB is the Boltzmann constant, T is the temperature, and (R) is the potential of mean September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 3 F igure 1. The illustration of the problem domains for Poisson-Boltzmann and Smoluchowski Equations. The Poisson-Boltzmann equation is solved on , while the Smoluchowski equation is solved only on s . represents the molecular surface, while a and r correspond the reactive and nonreactive boundaries. force (PMF) for a diffusing particle due to solvent mediated interactions with the target molecule. For simplicity, D j (R) can be assumed to be constant. The two terms contributing to the flux have clear physical meanings. The first is due to free diffusional processes, as quantified by Fick's first law. The second contribution is due to the drift velocity - (R) induced by the systematic forces - (R) and friction quantified by the friction constant . The relation between diffusion coefficient2D and friction constant is given by Stokes-Einstein equation: D = 1. = n0 e2 /(0 kB T ), 0 is the vacuum permittivity and n0 is the bulk ionic concentration. To accurately solve the Poisson-Boltzmann equation, the potential ( R) has been decomposed into singularized (s ) and regularized (r ) components 22 . The singularized component s can be further decomposed into the harmonic (h ) and the singular green function G: s = h + G. The h can be obtained via solving the September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 4 below equation: · (h (R)) = 0 on m h = -G on (4) Where the green function G is simply the electrostatic potential induced by all the singular charges inside the biomolecule, i.e., · (m G(R)) = qi (R - Ri ) on m i (5) qi and can be directly given by G(R) = i |R-R | m i Eq. 4 can be solved using the finite element method or boundary element method. For r , we need to solve the below equation: · (r (R)) + (R)2 sinh(r + s ) = 0 on r = 0 on s r = -m on n n (6) Since (R) = 0 on m and s = 0 on s , therefore Eq. 6 can be further simplified as · (r (R)) + (R)2 sinh(r ) = [s ] Where [s ] is the jump on the molecular interface. Finally, the solution of the PB equation cab be given by: = r + h + G (8) (7) The Smoluchowski equation (Eq. 2) can be solved to determine biomolecular diffusional encounter rates before steady state is established. Following the work of Song et al.9,10 and Zhou et al.23,24,25 , the application of the Smoluchowski equation to this question involves the solution of Eq. 2 in a three-dimensional domain , with the following boundary and initial conditions. A bulk Dirichlet condition is imposed on the outer boundary b , p j (Rl ; t ) = pbul k , f or Rl b , j (9) where pbul k denotes the bulk concentration at the outer boundary. A reactive j Robin condition is implemented on the active site boundary a , n(R0 ) · p(R0 ; t ) = (R0 ) p(R0 ), f or R0 a , (10) September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 5 providing an intrinsic reaction rate (R0 ). Here, n(R0 ) is the surface normal. For the diffusion-limited reaction process, such as inhibitor binding in N1, the concentration of the inhibitor at the binding site is approximately zero. Therefore, the reactive Robin condition on the inner boundary can be simplified as: p j (R0 ; t ) = 0, f or R0 a , (11) Where p j in the inhibitor concentration. It must be noted that for coions and counterions, there are no reactive boundaries. For the non-reactive surface parts of the inner boundary r , a reflective Neumann condition is employed. n(R0 ) · j p(R0 ; t ) = 0, Finally, we set up the initial conditions as p p j (R; 0) = 0 bul k j (12) |R| = l |R| < l (13) The solution of the PBS equation can be obtained by sequentially solving the Eq. 1 and Eq. 2; To solve the PNP equation, we need to solve the Eq. 1 and Eq. 2 iteratively, until the potential value converges to the initial threshold set by the user. Therefore, the diffusion-determined biomolecular reaction rate constant during the simulation time can be obtained from the flux j(R; t ) by integration over the active site boundary, i.e. 1 k(t ) = p-bul k j; Z a n(R) · j(R; t )d S (14) Finally, the corresponding concentration distribution can be obtained by R p j (R; t ) = e-( ) u(R; t ). 2.2. Adaptive Finite Element Mesh Generation Before generating the tetrahedral meshes for the external and internal space for the biomolecule, we first generate the surface triangular meshes for the molecule using the MSMS code (http://www.scripps.edu/ sanner/html/msms home.html) 26 . Certainly, Although MSMS is very efficient for surface mesh generation, the quality of the surface mesh still need to improve. Fortunately, the "Adventure Project" provides the excellent surface mesh smoothing for the MSMS surface mesh (http://adventure.q.t.u-tokyo.ac.jp/). With the molecular surface mesh, Tetgen was implemented to generate adaptive finite element meshes for biomolecules (http://tetgen.berlios.de/). September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 6 In the present study, molecular dynamics studies on the wild-type N1 structure (PDB code: 2HU4) have obtained "open" and "closed" conformations by Amaro et al. (unpublished results). Adaptive meshes for the "open" and "closed" conformations have been generated respectively. To compare the difference between the wild type and the mutants, we also prepared the meshes for the His274Tyr mutant (PDB code: 3CL0 and 3CKZ) and the Asn294Ser mutant (PDB code: 3CL2). The active site of the anti-neuraminidase inhibitors has been labeled as reactive boundaries. Reactive boundaries were defined following the typical BD methods: a spherical reactive surface is defined at an arbitrary radius from the biomolecular active site 27 . Generally, the user can label the active site with several spheres and write into the input file. 3. Results and Discussion 3.1. Validation of the PBS and PNP solver with A Spherical Test Case Before applying our SMOL code to a biomolecular system with complicated geometries, we first tested it with the classic spherical system 28 and compared the calculated result with the known analytical solution. For this test case, we chose a ° diffusing sphere with a 2 A radius and +1e charge. The receptor molecule was rep° resented with a sphere with a 10 A radius and -1e. The diffusion domain is another ° radius, which was discretized with 141,792 tetraconcentric sphere with a 400 A hedral elements. A detailed view of the surface mesh for the stationary sphere is also shown in Fig. 2(a). The diffusing particle's dimensionless bulk concentration was set to 1 mM . The concentration of coion and counterion was set at 150 mM . The dielectric constant is 2 in the receptor and 80 in the diffusion domain. Ignoring hydrodynamic interactions, the diffusion constant D is calculated °2 as 7.8 × 104A /µs using the Stokes-Einstein equation with a hydrodynamic radius ° of 3.5 A, solvent viscosity of 0.891 × 10-3kg/(m · s), and 300 K temperature. For a unit sphere with radius r0 and charge q in the center of the sphere, the linearized Poisson-Boltzmann equation have the below analytical solution a : a = e-r qer0 · , rR (1 + r0) r (20) Where R is the outer boundary. Fig. 3 shows the calculated potential comparing with the analytical values. The relative error is less than %3. Generally, under physiological condition, the ionic strength is approximately 0.15 M, here we solve the PBS and PNP equations with 150 mM N a+ and Cl - and 1 mM +e charged inhibitor. If the active site of the receptor is as shown in September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 7 F igure 2. (a) the external diffusion domain; (b) the non reactive molecular surface (cyan), the active site (red) and the outer sphere (pink). Fig. 2(b), we can obtain the potential value ((R)) and concentration distribution ( p j (R; t )) on each node in the diffusion domain. The total free energy for a given distribution p j (R; t ), i.e., the quantity which develops towards a minimum during the diffusion process, is a functional defined through G[ p j (R; t )] = Z g( p j (R; t ))d r (21) where g( p j (R; t )) is the free energy density connected with p j (R; t ). g( p j (R; t )) = p j (R; t )((R) + -1l n p j (R; t ) ) pbul k j (22) Similarly, we define the local entropy density s j (R; t ) as s j (R; t ) = -kB p j (R; t )l n p j (R; t ) pbul k j (23) Solving the PBS equation, we obtained kon = 9.967 × 1010M -1 · min-1 when the steady state comes. The distribution of N a+ , Cl - and +1e inhibitor is shown September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 8 Figure 3. (a) the calculated and analytical potential in the diffusion domain; (b) the relative error for calculated potential in the diffusion domain. in Fig.4. Furthermore, the equilibrium has not been reached due to the influx of all the ions. To accurately describe the equilibrium state, we need to iteratively solve the PB and Smoluchowski equations until the potential and diffusion rate constant reach the stable value, which is the "so-called" PNP equation. The final kon is 6.517 × 1010M -1 · min-1 . It seems that coions and counterions play the role in compensating the fast consumption of the +1e charged inhibitor. Comparing to the PBS model in Fig. 4, it clearly shows that there are more N a+ ions are pumped out of the diffusion domain, while more Cl - ions distribute near the receptor in our PNP model. The local free energy density and entropy density are much more tight and close to the receptor surface in the PNP model (Fig. 5). 3.2. Application of the PNP solver on different neuraminidase models The active sites of all the neurominidase models are represented by two spheres, ° one is a 3 A radius with the center on the carboxyl oxygen of Thr-225 and the other ° is in the gorge and around 5 A away from the carboxyl-O of Thr-225. Finally, it ° 2 for the area of each active site. Similarly, we suppose turns out to be around 200 A September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 9 F igure 4. The concentration distribution calculated for N a+ , Cl - and +1e charged drug via the PBS and PNP models, respectively. that there is 150 mM N a+ and 150 mM Cl - , 1 mM -e inhibitor in the bulk. It must be noted that the potential drug inhibitor carries -e charge. First, the "open" and "closed" conformations of the wild-type neuraminidase have been modeled. When the steady state comes, the association reaction rate constant kon is 2.76 × 108 M -1 · min-1 for the "open" conformation, while 3.66 × 108M -1 · min-1 for the "closed" conformation. It turns out that the "closed" conformation have stronger binding affinity to anti-neuraminidase drugs. Experimental studies also observed the "tighter interaction with ligand" 20 . Meanwhile, we simulated three mutants 3CL0, 3CKZ and 3CL2. The calculated steady-state rate constants have been listed in Table.1. Comparing with the experimental data, our results are very consistent. Additionally, it must be noted that all the calculated rate constant is larger than experimental data. Actually, the active site of the neuraminidase is very flexible to change according to both the experimental finding and unpublished molecular dynamics simulations by Amaro et al. The size of the possible drugs, such as sialic acid, oseltmivir and zanamivir, is relatively big (at lease 30 atoms). Therefore, different from our previous acetylcholinesterase model 29 , the non-electrostatic interaction, such as VdW and steric effect, and the September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 10 F igure 5. The local free energy and entropy densities in the PBS and PNP models. wild-type open wild-type closed 3CL0 3CKZ 3CL2 rate(PBS) 18.1 23.1 11.9 12.4 16.5 rate(PNP) 14.6 19.1 7.32 8.17 11.8 experimental 2.52 2.52 0.24 0.35 1.1 steric effect, must be taken account into the simulation for accurately depicting the drug diffusion. Currently, we are trying to hybrid the Brownian dynamics and finite element method to obtain the collision reaction probability. Combining with the above diffusion coefficient, it can predict the diffusion association constant more accurately. September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 11 4. Conclusion In this study, we describe continuum-based methods for studying electrostatic diffusion in biomolecular systems. Specifically, we present the SMOL software package, a finite element-based set of tools for solving the electrostatics and diffusion to calculate ligand binding rate constants for large biomolecules under pre-steady-state and steady-state conditions. The main improvement from the previous SMOL solver 9,10,29 can be addressed as below: first, the new regularized Poisson-Boltzmann algorithm instead of the previous APBS. It makes the solver read the potential value from the same diffusion node and avoid the previous data mapping error; second, the iteratively solving Poisson-Boltzmann and Smoluchowski equations results in a stable PNP solver; third, local free energy and entropy densities have been calculated during the diffusion process. finally, a simple mesh generator has been included for preparing tetrahedral meshes for the diffusion domain and molecular domain. Although our above solution for the neuraminidase is qualitatively consistent with experiments, the continuum method on studying the diffusion still has two main limitations: the size of the diffusion ligand and non-electrostatic interaction are not included. Currently, we are combining the continuum modeling method and Brownian dynamics together for simulating the drug diffusion in the neuraminidase. According to Fig. 4, it must be noted that only a quite limited area close to the active site has substantial concentration/free energy density gradient. Our new method singles out the active site to implement Brownian dynamics while continuum method on other areas. The ultimate goal of this work is to develop scalable methods and theories that will allow researchers to begin to study biological macromolecules in a cellular context. References 1. 2. 3. 4. 5. D. L. Ermak and J. A. McCammon, J. Chem. Phys. 69, 1352 (1978). S. H. Northrup, S. A. Allison and J. A. McCammon, J. Chem. Phys. 80, 1517 (1984). N. Agmon and A. L. Edelstein, Biophys. J. 72, 1582 (1997). R. R. Gabdoulline and R. C. Wade, Methods 14, 329 (1998). J. R. Stiles and T. M. Bartol, Monte Carlo methods for simulating realistic synaptic microphysiology using MCell, in Computational Neuroscience: Realistic Modeling for Experimentalists, ed. E. D. Schutter (CRC Press, Inc., New York, 2000) pp. 87­ 127. J. L. Smart and J. A. McCammon, Biophys. J. 75, 1679 (1998). M. G. Kurnikova, R. D. Coalson, P. Graf and A. Nitzan, Biophys. J. 76, 642 (1999). Z. Schuss, B. Nadler and R. S. Eisenberg, Phys. Rev. E 6403 (2001). Y. H. Song, Y. J. Zhang, C. L. Bajaj and N. A. Baker, Biophys. J. 87, 1558 (2004). Y. H. Song, Y. J. Zhang, T. Y. Shen, C. L. Bajaj, A. McCammon and N. A. Baker, 6. 7. 8. 9. 10. September 21, 2008 21:30 Proceedings Trim Size: 9in x 6in chengpsb20092 12 Biophys. J. 86, 2017 (2004). 11. K. S. Tai, S. D. Bond, H. R. Macmillan, N. A. Baker, M. J. Holst and J. A. McCammon, Biophys. J. 84, 2234 (2003). 12. H. Berry, Biophys. J. 83, 1891 (2002). 13. D. Genest, Biopolymers 28, 1903 (1989). 14. M. J. Saxton, Biophys. J. 61, 119 (1992). 15. J. A. McCammon, Science 238, 486 (1987). 16. S. H. Northrup, J. O. Boles and J. C. L. Reynolds, Science 241, 67 (1988). 17. R. C. Wade, M. E. Davis, B. A. Luty, J. D. Madura and J. A. Mccammon, Biophys. J. 64, 9 (1993). 18. P. Eastman and S. Doniach, Proteins 30, 215 (1998). 19. L. Yeomans-Reyna and M. Medina-Noyola, Phys Rev E Stat Nonlin Soft Matter Phys. 64, p. 066114 (2001). 20. R. J. Russell, L. F. Haire, D. J. Stevens, P. J. Collins, Y. P. Lin, G. M. Blackburn, A. J. Hay, S. J. Gamblin and J. J. Skehel, Nature 443, 45 (2006). 21. P. J. Collins, L. F. Haire, Y. P. Lin, J. F. Liu, R. J. Russell, P. A. Walker, J. J. Skehel, S. R. Martin, A. J. Hay and S. J. Gamblin, Nature 453, 1258 (2008). 22. I. Chern, J. Liu and W. Wang, Methods Appl. Anal. 10, 309 (2003). 23. H. X. Zhou, J. Chem. Phys. 92, 3092 (1990). 24. H. X. Zhou, S. T. Wlodek and J. A. McCammon, Proc. Natl. Acad. Sci. U. S. A. 95, 9280 (1998). 25. H. X. Zhou, J. M. Briggs, S. Tara and J. A. McCammon, Biopolymers 45, 355 (1998). 26. M. F. Sanner, A. J. Olson and J. Spehner, Proc. 11th ACM Symp. Comp. Geom. C6-C7 (1995). 27. S. Tara, A. H. Elcock, J. M. Briggs, Z. Radic, P. Taylor and J. A. McCammon, Biopoly´ mers 46, 465 (1998). 28. E. B. Krissinel and N. Agmon, J. Comput. Chem. 17, 1085 (1996). 29. Y. H. Cheng, J. K. Suen, Z. Radic, S. D. Bond, M. J. Holst and J. A. McCmmon, ´ Biophys. Chem. 127, 129 (2007).