Article Information

• PDF (896 kb)
• Full Text with Large Figures
• Cited by in Scopus (4)

PubMed

• Articles by Yuhui Cheng
• Articles by J. Andrew McCammon

Related Articles

• Finite Element Solution of the Steady-State Smoluchowski Equ...

  Finite Element Solution of the Steady-State Smoluchowski Equation for Rate Constant Calculations Biophysical Journal, Volume 86, Issue 4, 1 April 2004, Pages 2017-2029 Yuhua Song, Yongjie Zhang, Tongye Shen, Chandrajit L. Bajaj, J. Andrew McCammon and Nathan A. Baker Abstract This article describes the development and implementation of algorithms to study diffusion in biomolecular systems using continuum mechanics equations. Specifically, finite element methods have been developed to solve the steady-state Smoluchowski equation to calculate ligand binding rate constants for large biomolecules. The resulting software has been validated and applied to mouse acetylcholinesterase. Rates for inhibitor binding to mAChE were calculated at various ionic strengths with several different reaction criteria. The calculated rates were compared with experimental data and show very good agreement when the correct reaction criterion is used. Additionally, these finite element methods require significantly less computational resources than existing particle-based Brownian dynamics methods. Abstract | Full Text | PDF (587 kb)

• Continuum Diffusion Reaction Rate Calculations of Wild-Type ...

  Continuum Diffusion Reaction Rate Calculations of Wild-Type and Mutant Mouse Acetylcholinesterase: Adaptive Finite Element Analysis Biophysical Journal, Volume 87, Issue 3, 1 September 2004, Pages 1558-1566 Yuhua Song, Yongjie Zhang, Chandrajit L. Bajaj and Nathan A. Baker Abstract As described previously, continuum models, such as the Smoluchowski equation, offer a scalable framework for studying diffusion in biomolecular systems. This work presents new developments in the efficient solution of the continuum diffusion equation. Specifically, we present methods for adaptively refining finite element solutions of the Smoluchowski equation based on a posteriori error estimates. We also describe new, molecular-surface-based models, for diffusional reaction boundary criteria and compare results obtained from these models with the traditional spherical criteria. The new methods are validated by comparison of the calculated reaction rates with experimental values for wild-type and mutant forms of mouse acetylcholinesterase. The results show good agreement with experiment and help to define optimal reactive boundary conditions. Abstract | Full Text | PDF (306 kb)

• Finite Element Simulations of Acetylcholine Diffusion in Neu...

  Finite Element Simulations of Acetylcholine Diffusion in Neuromuscular Junctions Biophysical Journal, Volume 84, Issue 4, 1 April 2003, Pages 2234-2241 Kaihsu Tai, Stephen D. Bond, Hugh R. MacMillan, Nathan Andrew Baker, Michael Jay Holst and J. Andrew McCammon Abstract A robust infrastructure for solving time-dependent diffusion using the finite element package FEtk has been developed to simulate synaptic transmission in a neuromuscular junction with realistic postsynaptic folds. Simplified rectilinear synapse models serve as benchmarks in initial numerical studies of how variations in geometry and kinetics relate to endplate currents associated with fast-twitch, slow-twitch, and dystrophic muscles. The flexibility and scalability of FEtk affords increasingly realistic and complex models that can be formed in concert with expanding experimental understanding from electron microscopy. Ultimately, such models may provide useful insight on the functional implications of controlled changes in processes, suggesting therapies for neuromuscular diseases. Abstract | Full Text | PDF (451 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 10, 3397-3406, 15 May 2007

doi:10.1529/biophysj.106.102533

Biophysical Theory and Modeling

Finite Element Analysis of the Time-Dependent Smoluchowski Equation for Acetylcholinesterase Reaction Rate Calculations

Yuhui Cheng^*, †, ,  , Jason K. Suen^†, Deqiang Zhang^‡, Stephen D. Bond^§, Yongjie Zhang^¶, Yuhua Song^||, Nathan A. Baker^||, Chandrajit L. Bajaj^¶, **, Michael J. Holst^†† and J. Andrew McCammon^{*, †, ‡‡}

^* Howard Hughes Medical Institute, University of California at San Diego, La Jolla, California
^† Department of Chemistry and Biochemistry, and Center for Theoretical Biological Physics, University of California at San Diego, La Jolla, California
^‡ Accelrys, San Diego, California
^§ Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, Illinois
^¶ Institute of Computational Engineering and Sciences, Center for Computational Visualization, University of Texas at Austin, Texas
^|| Department of Biochemistry and Molecular Biophysics, Center for Computational Biology, Washington University in St. Louis, Missouri
^** Department of Computer Sciences, University of Texas at Austin, Texas
^†† Department of Mathematics, University of California at San Diego, La Jolla, California
^‡‡ Department of Pharmacology, University of California at San Diego, La Jolla, California

Address reprint requests to Y. Cheng, Tel.: 858-822-2771.

Diffusion plays an 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 Carlo 5,12,13,14, Brownian dynamics (BD) 15,16,17, and Langevin dynamics 18,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 need not deal with the individual Brownian particles and the computational cost can be substantially less than for the discrete methods.

In previous work with continuum methods, Song et al. 9,10 have presented finite element methods for solving the steady-state Smoluchowski equation (SSSE), which describes the steady-state behavior of diffusion-limited ligand binding. These methods have been shown to be significantly more efficient than traditional BD approaches for evaluating reaction rate constants for diffusion-limited binding of simple ligands. Recently, Zhang et al. 20 applied this approach for studies of several conformations of tetrameric mouse acetylcholinesterase (mAChE). However, the SSSE solution only provides the answer at the time-independent stage of diffusion. In other words, we only obtain the concentration distribution and rate constant when diffusion and reaction between the ligand and the enzyme reach the steady state. Physiological conditions, however, can be expected to include non-steady-state kinetics. One possible way to study the diffusion dynamics on biomolecular interface binding energy landscape is mean first-passage time, which was introduced recently by Wang et al. 21. The theory suggests a way of connecting the models/simulations with single molecule experiments by analyzing the kinetic trajectories. However, it is still an open question for the diffusional problem in a large spatial and timescale.

In the present work, we apply adaptive finite element methods to solve the time-dependent Smoluchowski equation (TDSE), using a posteriori error estimation to iteratively refine the finite element meshes. The binding of charged and noncharged ligands to mAChEs has been described at each timestep. The diffusion results have been compared with those from steady-state Smoluchowski diffusion studies and experimental results. AChE is a serine esterase that terminates the activity of acetylcholine (ACh) within cholinergic synapses by hydrolysis of the ACh ester bond to produce acetate and choline 22. Hydrolysis of ACh occurs in the active site of AChE, which lies at the base of a 20Å-deep gorge within the enzyme. The rate-limiting step of ACh hydrolysis by AChE is the diffusional encounter 23,24,25, making the system a popular target for both experimental 26,27,28 and computational diffusion studies 29,30.

Our time-dependent SMOL solver (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 ( being a particle’s position) under conditions of high friction, where the Smoluchowski equation applies.

The starting point for development of the time-dependent SMOL solver is the steady-state SMOL solver described by Song et al. 9,10. The original Smoluchowski equation has the form of a continuity equation,

(1)

(2)

Here is the distribution function of the ensemble of Brownian particles, is the diffusion coefficient, β=1/k_BT is the inverse Boltzmann energy, k_B is the Boltzmann constant, T is the temperature, and is the potential of mean force (PMF) for a diffusing particle due to solvent-mediated interactions with the target molecule. For simplicity, can be assumed 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, induced by the systematic forces, and friction quantified by the friction constant γ. The relation between diffusion coefficient D and friction constant γ is given by Stokes-Einstein equation: Dβγ=1.

The TDSE 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. 31,32,33, the application of the TDSE to this question involves the solution of Eq. (1) in a three-dimensional domain Ω, with the following boundary and initial conditions. A bulk Dirichlet condition is imposed on the outer boundary Γ_b⊂∂Ω,

(3)

(4)

(5)

For the nonreactive surface parts of the inner boundary Γ_r⊂∂Ω, a reflective Neumann condition is employed:

(6)

Finally, we set up the initial conditions as

(7)

Therefore, the diffusion-determined biomolecular reaction rate constant during the simulation time can be obtained from the flux by integration over the active site boundary, i.e.,

(8)

To numerically solve the TDSE, we employed the Galerkin finite element approximation to discretize the differential equation 34. The original TDSE (Eq. (1)) can be written as described below 10,35,36.

Let Ω ⊂³ be an open set, and let ∂Ω denote the boundary, which can be thought of as a set in ². Consider now the TDSE, a member of the class of elliptic equations

(9)

According to Holst et al. 37, the solution to the original problem also solves the problem

(10)

(11)

We have used the fact that a boundary integral vanishes because the test function v vanishes on the boundary.

For a discrete solution to Eq. (11), taking span{ϕ₁, ϕ₂, …, ϕ_N} ⊂ Eq. (11) reduces to a set of N nonlinear algebraic relations (implicitly defined) for N coefficients {α_j} in the expansion:

(12)

According to the Galerkin approximation, N equals the number of finite element nodes.

Therefore, the corresponding “weak form” of the TDSE is

(13)

To obtain an unconditionally stable solution, two implicit algorithms have been implemented in our codes: Crank-Nicolson and backward Euler’s methods.

Finally, the concentration distribution can be obtained by .

As described by Holst et al. 37, the adaptive mesh refinement procedure follows a solve-estimate-refine algorithm and has been implemented in the FEtk software (http://www.fetk.org/). Because of the inefficiency to “estimate” and “refine” in each time step, we only estimated and refined the mesh while solving the SSSE. With the refined mesh, TDSE diffusion studies were implemented. In the “estimate” step, we introduced the a posteriori error estimator η_s below holding for a Galerkin approximation u_h satisfying

(14)

(15)

(16)

Currently we provide two options to map the PMF to each finite element node in the time-dependent SMOL solver code. First, it can input the PMF obtained by boundary element methods 39,40. The PMF corresponds to the electrostatic potential obtained by solving the Poisson-Boltzmann equation. Second, APBS 0.4.0 (http://sourceforge.net/projects/apbs) is used to calculate the PMF, which is the potential field W(r) in Eq. (2)36. The partial charges and radii of each atom in the mAChE monomer and tetramer molecules have been assigned using the CHARMM22 force field, and the dielectric constant is set as 4.0 inside the protein and 78.0 for the solvent. The solvent probe radius is set as 1.4Å, and the ion exclusion layer is set as 2.0Å. Ionic strengths varying between 0 and 0.67M were used in the PMF calculations and following diffusion studies.

To allow the potential to approach zero at the outer boundary, a large space of 40 times the radius of the biomolecule is required. A series of nested potential grids is constructed in a multiresolution format where higher resolution meshes provide PMF values near the molecular surface while coarser meshes are used away from the molecule. The dimensions of the finest grid are given by the psize.py utility in the APBS software package, and the coarsest grid dimensions are set to cover the whole problem domain plus two grid spacings (to allow gradient calculation) in each dimension. The setup for the rest of the grid hierarchy is calculated using a geometric sequence for grid spacing. For mAChE monomer, the finest grid has dimensions of 86.3Å×76.4Å×101.4Å with 161, 129, and 193 grid points in each direction, respectively. This corresponds to a 0.5Å×0.6Å×0.5Å grid spacing setup. The coarsest grid has dimensions of 3400Å×3000Å×4300Å with 161 grid points in each direction. The corresponding grid spacing settings are 21.1Å×18.6Å×26.7Å.

For the mAChE monomer case, similarly with previous studies 9,10, we used a mouse AChE (mAChE) structure adapted from the crystal structure of the mAChE-fasciculin II complex 26 and perturbed by Tara and co-workers via molecular dynamics simulations with an ACh-like ligand in the active site gorge 30 to produce gorge conformations with wider widths than the original x-ray structure. The diffusing ligand was modeled as a sphere with an exclusion radius of 2.0Å and a diffusion constant of 7.8×10⁴Ų/μs. This perturbation was necessary for computational diffusion simulations with a fixed biomolecular structure. Reactive boundaries were defined using the biomolecular surface, which is the same as that used in Song et al. 10.

For the mAChE tetramer cases, we used three structures: a loose, pseudo-square planar tetramer with antiparallel alignment of the two four-helix bundles and a large space in the center (PDB: 1C2B); a compact, square nonplanar tetramer with parallel arrangement of the four-helix bundles that may expose all the four t peptide sequences on a single side (PDB: 1C2O); and in addition to the crystal structures, an intermediate structure (INT) was generated by morphing the two crystal structures using the morph script in visual molecular dynamics 41. Reactive boundary definitions are the same as the above mAChE monomer case.

The tetrahedral meshes were obtained and refined from the inflated van der Waals-based accessibility data for the mAChE monomer and tetramers using the level-set boundary interior exterior-mesher 42,43,44. Initially the region between the biomolecule and a slightly larger sphere centered about the molecular center of mass, was discretized by adaptive tetrahedral meshes. It generated very fine triangular elements near the active site gorge, while coarser elements everywhere else. The mesh is then extended to the entire diffusion domain and the inside of the biomolecule with spatial adaptivity in that the mesh element size increases with increasing distance from the biomolecule. The number of tetrahedral elements varies from 50,000 to 70,000 for different tetramer geometries.

Before applying the time-dependent SMOL program to a biomolecular system with complicated geometries, we first tested it with the classic spherical system 45 and compared the calculated result with the known analytical solution. For this test case, we chose a diffusing sphere with a 2Å radius and neutral charge. The entire problem domain is a sphere with a 400Å radius, which was discretized with 745,472 tetrahedral elements. A detailed view of the surface mesh for the stationary sphere is also shown in Fig. 1. The time-dependent Smoluchowski equation then becomes the Einstein diffusion equation. The diffusing particle’s dimensionless bulk concentration was set to 1. Ignoring hydrodynamic interactions, the diffusion constant D is calculated as 7.8×10⁴Ų/μs using the Stokes-Einstein equation with a hydrodynamic radius of 3.5Å, solvent viscosity of 0.891×10⁻³kg/(m · s), and 300K temperature.

Display large version of this figure

Figure 1

Illustration of the discretized problem domain for the spherical analytical test. The green represents the outer boundary, in which the ligand concentration is kept as a constant.

For a spherically symmetric system without external potential, the TDSE can be written as

This analytical form of the solution was expressed by the sum of zero-order spherical Bessel functions. Fig. 1 presents the concentration distributions during the simulation time with our TDSE solver, comparing with the above analytical solution.

According to Fig. 2, the performance of the SMOL program is good, with almost the same concentration distribution as in the analytical solution. It must be noted that the analytical solution for the time-dependent diffusion with the Columbic potential cannot be addressed with a simple formula; however, we have implemented our solver to test the same steady-state case addressed in Table 1 of Song et al. 10, and obtained very consistent results.

Display large version of this figure

Figure 2

Time evolution of the two-dimensional concentration distribution contour in the problem domain. (a) Our SMOL solution; (b) spherical analytical solution.

One of the major advantages of continuum methods such as the time-dependent SMOL solver is the ability to simulate the complete diffusion dynamics for large biological systems with complicated geometries with significantly lower computational cost than the Brownian dynamics approach. This section demonstrates the implementation of TDSE to study the ligand binding kinetic process of the mAChE monomer under various ionic strength conditions 46.

With the original mesh, we measured the diffusion-controlled reaction rates during the simulation time with the timestep at 50ps, as shown in Fig. 3. Separate calculations were performed at ionic strengths of 0.000, 0.050, 0.100, 0.150, 0.200, 0.250, 0.300, 0.450, 0.600, and 0.670M. At the zero ionic strength, the whole system reaches the steady state in over 15ns. The value of k_on at the end of the simulation is 9.535×10¹¹M⁻¹ · min⁻¹, which is very consistent with the experimental value at (9.80±0.60)×10¹¹M⁻¹ · min⁻¹27. Meanwhile, the k_on value for the neutral ligand at the steady state is 9.297×10¹⁰M⁻¹ · min⁻¹, which is consistent with the previous steady-state calculations 20. Table 1 (this article) listed the final k_on value derived from the TDSE calculations and the corresponding sets from SSSE calculations 9. When the ionic strength becomes higher, the time to reach the steady state decreases substantially. Obviously, we have obtained consistent results, comparing with the previous SSSE and BD calculations.

Display large version of this figure

Figure 3

k_on(t) values in time-dependent ACh diffusion under the various ionic strength conditions.

Furthermore, our TDSE solver can report all the ligand concentration distribution histories during the 20μs simulation. In this case, we recorded a concentration distribution every 100 timesteps (5ns) and a restart checkpoint every 1000 timesteps (50ns). Fig. 4 demonstrates the two-dimensional concentration distribution around the mAChE molecule at the end of 20μs simulation. The origin and the normal of the clip plane have been set at (0 16.6Å 0) and (1 0 1), respectively. The value k_on exhibits an ionic strength-dependence strongly indicative of electrostatic acceleration. The high ionic strength environment lessens the electrostatic interactions between the active site of mAChE and the ligand. Therefore, the relatively low ligand concentration area shrinks while ionic strength increases. Specifically, at 0.150M ionic strength, several dynamic snapshots have been plotted, as depicted in Fig. 2. These snapshots demonstrate the whole diffusion process of ACh-like ligands from the area far from the enzyme until they reach the active site and react and disappear.

Display large version of this figure

Figure 4

Two-dimensional ligand concentration distribution contour around mAChE at the steady state under various ionic strengths. The red color represents high concentration area, while the blue represents low concentration area.

In this section, we explore the use of the adaptive finite element methods to implement TDSE calculations on the mAChE monomer. The first step is to interactively solve the SSSE based on the a posteriori error estimation 10,37. The iterative error-based refinement of the initial 656,823-simplex mesh was performed until the global error is less than a value chosen to provide reaction rates that did not change appreciably upon further refinement. The refined mesh has 824,746 simplexes at 150mM ionic strength. Then, we implemented another TDSE calculation with the refined mesh. Fig. 5 shows the kinetic curves of both the original and refined mesh. Again, the two calculations are in good overall agreement but do show some differences at the final steady state. Specifically, the refinement of the adaptive meshes shows a little increase of the k_on value after reaching the steady state.

Display large version of this figure

Figure 5

The time-dependent ligand concentration distribution at 0.150M ionic strength for the mAChE monomer.

A previous SSSE study described the effect of electrostatic forces on ACh steady-state diffusion to the mouse acetylcholinesterase tetramer 20. Here, we extend the previous study using the same meshes and potential files with our time-dependent SMOL solver. The time step for the three tetramer models was set at 10ns, and concentration histories were recorded every 100 steps. Two crystal structures (1C2O and 1C2B) and an intermediate structure (INT) are all studied by solving the TDSE. The actual conformational dynamics of the mouse acetylcholinesterase tetramer has been neglected in this work.

Figure 6a shows the time-dependent rate constant per active site at 0.150M ionic strength for the above three mouse acetylcholinesterase tetramer structures. It takes >75μs for each active site to reach the steady state. For structure 1C2O, the entrances to two of the four active gorges (AS2 and AS4) are partially blocked by another subunit in the complementary dimer, while the other two gorges are completely accessible from outside (AS1 and AS3). As a result, the four kinetic curves in 1C2O can be classified into two subgroups: one subgroup corresponds to active site 1 (AS1) and active site 3 (AS3), in each of which the gorge is open, and at the end of the simulation, the reaction rates are 1.61×10¹¹M⁻¹min⁻¹ and 1.50×10¹¹M⁻¹min⁻¹, respectively. Another subgroup corresponds to active site 2 (AS2) and active site 4 (AS3), where the gorges are sterically shielded by nearby subunits, and the final reaction rates are 8.47×10¹⁰M⁻¹min⁻¹ and 9.62×10¹⁰M⁻¹min⁻¹, respectively, which is a little more than half of that for AS1 or AS3. Figure 6b demonstrates the k_on(t) values for the neutral ligand. Comparing with the +1.0e charged ligand, the neutral ligand still shows similar time-dependent curves for individual active sites, while the k_on(t) value is much less than the corresponding +1.0e charged case.

Display large version of this figure

Figure 6

The comparison of k_on values in time-dependent ACh diffusion between the original and refined meshes.

Similarly, we tested the final steady-state concentration distribution in the diffusion domain under various ionic strength conditions for the AChE tetramers. For example, Fig. 8 illustrates the different concentration profiles at 0.025, 0.050, 0.100, 0.300, 0.450, and 0.670M ionic strength solutions for structure 1C2O. Comparing with the monomer case, the ACh-like ligand concentration around the 1C2O molecule is much lower when the ionic strength is small, due to the stronger electrostatic attraction between the ligands and the tetramer molecule. While the ionic strength becomes higher, the electrostatic effect on the steady-state concentration distribution turns out to be weaker.

We also obtained the time-dependent rate constant per active site for the structure 1C2B, in which all the four gorges are nearly completely accessible to the solvent (Fig. 7). The profiles of k_on(t) of AS3 and AS4 are almost the same, while the value of the final steady state for AS1 or AS2 is a little smaller, but still above 1.10×10¹¹M⁻¹min⁻¹. The sum of the four active sites is 5.28×10¹¹M⁻¹min⁻¹, whereas the total steady-state k_on in structure 1C2O is 4.92×10¹¹M⁻¹min⁻¹. It must be noted that the steady-state k_on in the mAChE monomer is 1.97×10¹¹M⁻¹min⁻¹ at 0.150M ionic strength. Therefore, the average reaction rate per active site for the tetramer is ∼64% that of the monomer, which is close to the result of the previous SSSE studies 20. The four active sites show similar kinetic profiles, and reach the steady states at nearly the same time. Meanwhile, the time-dependent rate constant per active site in the structure INT appears more like that in the structure 1C2B, although the difference between AS1 and AS3 is still similar with that in the structure 1C2O. Additionally, our time-dependent SMOL solver can show the detailed diffusion process. For the 1C2O case, Fig. 8 describes the ligand concentration distribution in the problem domain. The red represents high concentration areas, while the blue represents the low concentration areas.

Display large version of this figure

Figure 7

The dependency of k_on(t) values on the simulation time for the structures 1C2O, INT, and 1C2B: (a) 0.150M ionic strength and +1.0e ligand; and (b) neutral ACh-like ligand.

Display large version of this figure

Figure 8

Steady-state ligand concentration distribution under six different ionic strength conditions for structure 1C2O.

Fig. 9In this study, we describe new continuum-based methods for studying diffusion in biomolecular systems. Specifically, we present the time-dependent SMOL software package, a finite element-based set of tools for solving the TDSE 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 can be addressed as below: first, the new solver has removed the drift term (Eq. (2)) which is discontinuous for ∇W, as well as the asymmetry 47. Theoretically, the new SMOL solver can utilize the conjugate gradient method, which is a direct method for symmetric and positive definite linear systems, while the old solver can only solve SSSE with the bi-conjugate gradient method. Compared with the new solver, the bi-conjugate gradient method is slower and harder to converge. Therefore, our new SMOL solver can solve both steady-state and time-dependent problems much more efficiently and stably than the old version.

Display large version of this figure

Figure 9

Ligand concentration distribution in the diffusion domain during the simulation time for structure 1C2O.

With the new code, we solved the time-dependent diffusion in the analytical case of a reactive sphere, mAChE monomer and tetramer cases. Comparing with previous steady-state studies, our research extends the study into the nonequilibrium diffusion dynamics and obtained very consistent results. Moreover, the calculated rates of the mAChE monomer were compared with experimental data 27 and show very good agreement with experiment while requiring substantially less computational effort than existing particle-based Brownian dynamics methods. Additionally, the value of k_on(t) seems to be underestimated with the coarser meshes, which is consistent with previous observations 10. Similarly, the k_on values in mAChE tetramers should increase if we refine the original mesh. In the previous study 20 and this one, we have found the activity of one subunit in a mAChE tetramer equals ∼60–70% that of a free monomer. With the appropriate meshes, we would expect to obtain an activity closer to that in the free monomer and the catalytic activity might not be too affected by subunit association as suggested in the experiment 48.

Additionally, we describe new adaptive meshing methods developed to discretize biomolecular systems into finite element meshes, which respect the geometry of the biomolecule. Although not presented in this study, it is important to note that the new meshing methods could be useful in a variety of biological simulations including computational studies of biomolecular electrophoresis 49, elasticity 42,43, and electrostatics 35,36,50,51.

Finally, this research lays the groundwork for the integration of molecular-scale information into simulations of cellular-scale systems such as the neuromuscular junction 6,11,52. In particular, this new finite element framework should facilitate the incorporation of other continuum mechanics phenomena into biomolecular simulations. The ultimate goal of this work is to develop scalable methods and theories that will allow researchers to begin to study biological dynamics in a cellular context efficiently and robustly.