Finite element simulations of acetylcholine diffusion in neuromuscular junctions

Nanosensing of Acetylcholine Molecules: Influence of the Association Mechanism

A bullet-shaped nanopore surface modified by two polyelectrolyte (PE) layers, an inner polyethyleneimine (PEI) layer and an outer p-sulfonatocalix[4]-arene (SCX4) layer, is applied to sense trace levels of acetylcholine (Ach) molecules. We show that the higher the order of the association reaction of Ach with SCX4, the smaller the difference between the ionic current when Ach is present and that when it is absent, and so is the difference in the space charge density. In addition, the larger the binding constant K of that reaction, the lower the detection limit but narrower the detection range. Choosing pH 7 is most appropriate because if the pH is low, the concentration polarization of H⁺ is significant, and as it gets high, both PE layers become uncharged. At pH 7 and K = 2 × 10⁷ L/mol, the detection limit of the nanopore ranges from 1 to 10 nM, which is orders of magnitude lower than that of the other approaches.

Potential alternatives to current cholinesterase inhibitors: an in silico drug repurposing approach

Acetylcholinesterase/Butyrylcholinesterase inhibitors are considered an effective method for treating Alzheimer's disease (AD). In this current work, we have computationally analyzed 11 new small molecule drugs used in various neurological diseases and Donepezil, a known inhibitor of acetylcholinesterase, as a positive control. We investigated these drugs for possible fundamental interactions with acetylcholinesterase and butyrylcholinesterase as both are critical in the pathophysiology of Alzheimer's disease. We have selected FDA approved compounds for repurposing as possible inhibitors of these enzymes and novel therapeutic option for Alzheimer's disease. We selected the top two molecules for each protein for their binding energies, interactions, and Donepezil, the most commonly used drug for AD treatment. Molecular simulation and dynamics studies of the top 2 drugs in each case and free energy analysis helped us reach further conclusions about the best possible drugs for repurposing. Brexipirazole and Deutetrabenazine produce encouraging results as butyrylcholinesterase and acetylcholinesterase inhibitors, respectively.

Kinetics of Chemical Processes in the Human Brain. The Cholinergic Synapse-Mechanisms of Functioning and Control Methods

In the framework of the kinetic model, the functioning of the cholinergic synapse is considered. The results of mathematical modeling of changes in the level of acetylcholine, induced pH impulse, the influence of the frequency of impulse transmission and inhibition of acetylcholinesterase are presented. Physicochemical explanation for a number of important physiological phenomena, such as neuromuscular paralysis, the molecular mechanism of neurological memory, and actions of nerve poisons and toxins, is given.

Kinetics of Chemical Processes in the Human Brain. Proton Blockade of Acetylcholinesterase and pH-Impulse in the Mechanism of Functioning of the Cholinergic Synapse

A kinetic model describing the dynamics of synaptic "discharge" taking into account the kinetics of the injection of the neurotransmitter into the synaptic cleft, the pH-dependence of the catalytic activity of the enzyme, and diffusion withdrawal of protons is proposed. The model provides a physicochemical explanation for a number of important physiological phenomena, such as the neuromuscular paralysis, the molecular mechanism of neurological memory, and the effect of some neurotoxins and drugs.

Improvements to the APBS biomolecular solvation software suite

The Adaptive Poisson-Boltzmann Solver (APBS) software was developed to solve the equations of continuum electrostatics for large biomolecular assemblages that have provided impact in the study of a broad range of chemical, biological, and biomedical applications. APBS addresses the three key technology challenges for understanding solvation and electrostatics in biomedical applications: accurate and efficient models for biomolecular solvation and electrostatics, robust and scalable software for applying those theories to biomolecular systems, and mechanisms for sharing and analyzing biomolecular electrostatics data in the scientific community. To address new research applications and advancing computational capabilities, we have continually updated APBS and its suite of accompanying software since its release in 2001. In this article, we discuss the models and capabilities that have recently been implemented within the APBS software package including a Poisson-Boltzmann analytical and a semi-analytical solver, an optimized boundary element solver, a geometry-based geometric flow solvation model, a graph theory-based algorithm for determining pK_a values, and an improved web-based visualization tool for viewing electrostatics.

Hybrid finite element and Brownian dynamics method for charged particles

Diffusion is often the rate-determining step in many biological processes. Currently, the two main computational methods for studying diffusion are stochastic methods, such as Brownian dynamics, and continuum methods, such as the finite element method. A previous study introduced a new hybrid diffusion method that couples the strengths of each of these two methods, but was limited by the lack of interactions among the particles; the force on each particle had to be from an external field. This study further develops the method to allow charged particles. The method is derived for a general multidimensional system and is presented using a basic test case for a one-dimensional linear system with one charged species and a radially symmetric system with three charged species.

Synaptic bouton properties are tuned to best fit the prevailing firing pattern

The morphology of presynaptic specializations can vary greatly ranging from classical single-release-site boutons in the central nervous system to boutons of various sizes harboring multiple vesicle release sites. Multi-release-site boutons can be found in several neural contexts, for example at the neuromuscular junction (NMJ) of body wall muscles of Drosophila larvae. These NMJs are built by two motor neurons forming two types of glutamatergic multi-release-site boutons with two typical diameters. However, it is unknown why these distinct nerve terminal configurations are used on the same postsynaptic muscle fiber. To systematically dissect the biophysical properties of these boutons we developed a full three-dimensional model of such boutons, their release sites and transmitter-harboring vesicles and analyzed the local vesicle dynamics of various configurations during stimulation. Here we show that the rate of transmission of a bouton is primarily limited by diffusion-based vesicle movements and that the probability of vesicle release and the size of a bouton affect bouton-performance in distinct temporal domains allowing for an optimal transmission of the neural signals at different time scales. A comparison of our in silico simulations with in vivo recordings of the natural motor pattern of both neurons revealed that the bouton properties resemble a well-tuned cooperation of the parameters release probability and bouton size, enabling a reliable transmission of the prevailing firing-pattern at diffusion-limited boutons. Our findings indicate that the prevailing firing-pattern of a neuron may determine the physiological and morphological parameters required for its synaptic terminals.

Hybrid finite element and Brownian dynamics method for diffusion-controlled reactions

Diffusion is often the rate determining step in many biological processes. Currently, the two main computational methods for studying diffusion are stochastic methods, such as Brownian dynamics, and continuum methods, such as the finite element method. This paper proposes a new hybrid diffusion method that couples the strengths of each of these two methods. The method is derived for a general multidimensional system, and is presented using a basic test case for 1D linear and radially symmetric diffusion systems.

A solvable model for the diffusion and reaction of neurotransmitters in a synaptic junction

BACKGROUND

The diffusion and reaction of the transmitter acetylcholine in neuromuscular junctions and the diffusion and binding of Ca2+ in the dyadic clefts of ventricular myocytes have been extensively modeled by Monte Carlo simulations and by finite-difference and finite-element solutions. However, an analytical solution that can serve as a benchmark for testing these numerical methods has been lacking.

RESULT

Here we present an analytical solution to a model for the diffusion and reaction of acetylcholine in a neuromuscular junction and for the diffusion and binding of Ca2+ in a dyadic cleft. Our model is similar to those previously solved numerically and our results are also qualitatively similar.

CONCLUSION

The analytical solution provides a unique benchmark for testing numerical methods and potentially provides a new avenue for modeling biochemical transport.

A synergic simulation-optimization approach for analyzing biomolecular dynamics in living organisms

A synergic duo simulation-optimization approach was developed and implemented to study protein-substrate dynamics and binding kinetics in living organisms. The forward problem is a system of several coupled nonlinear partial differential equations which, with a given set of kinetics and diffusion parameters, can provide not only the commonly used bleached area-averaged time series in fluorescence microscopy experiments but more informative full biomolecular/drug space-time series and can be successfully used to study dynamics of both Dirac and Gaussian fluorescence-labeled biomacromolecules in vivo. The incomplete Cholesky preconditioner was coupled with the finite difference discretization scheme and an adaptive time-stepping strategy to solve the forward problem. The proposed approach was validated with analytical as well as reference solutions and used to simulate dynamics of GFP-tagged glucocorticoid receptor (GFP-GR) in mouse cancer cell during a fluorescence recovery after photobleaching experiment. Model analysis indicates that the commonly practiced bleach spot-averaged time series is not an efficient approach to extract physiological information from the fluorescence microscopy protocols. It was recommended that experimental biophysicists should use full space-time series, resulting from experimental protocols, to study dynamics of biomacromolecules and drugs in living organisms. It was also concluded that in parameterization of biological mass transfer processes, setting the norm of the gradient of the penalty function at the solution to zero is not an efficient stopping rule to end the inverse algorithm. Theoreticians should use multi-criteria stopping rules to quantify model parameters by optimization.

Diffusion-convection effects on drug distribution at the cell membrane level in a patch-clamp setup

We present a model-based method for estimating the effective concentration of the active drug applied by a pressure pulse to an individual cell in a patch-clamp setup, which could be of practical use in the analysis of ligand-induced whole-cell currents recorded in patch-clamp experiments. Our modelling results outline several important factors which may be involved in the high variability of the electric response of the cells, and indicate that with a pressure pulse duration of 1s and diameter of the perfusion tip of 600 μm, elevated amounts of drug can accumulate locally between the pipette tip and the cell. Hence, the effective agonist concentration at the cell membrane level can be consistently higher than the initial concentration inside the perfusion tubes. We performed finite-difference and finite-element simulations to investigate the diffusion/convection effects on the agonist distribution on the cell membrane. Our model can explain the delay between the commencement of acetylcholine application and the onset of the whole-cell current that we recorded on human rhabdomyosarcoma TE671 cells, and reproduce quantitatively the decrease of signal latency with the concentration of agonist in the pipette. Results also show that not only the geometry of the bath chamber and pipette tip, but also the transport parameters of the diffusive and convective phenomena in the bath solution are determinant for the amplitude and kinetics of the recorded currents and have to be accounted for when analyzing patch-clamp data.

Poisson-Nernst-Planck Equations for Simulating Biomolecular Diffusion-Reaction Processes I: Finite Element Solutions

In this paper we developed accurate finite element methods for solving 3-D Poisson-Nernst-Planck (PNP) equations with singular permanent charges for electrodiffusion in solvated biomolecular systems. The electrostatic Poisson equation was defined in the biomolecules and in the solvent, while the Nernst-Planck equation was defined only in the solvent. We applied a stable regularization scheme to remove the singular component of the electrostatic potential induced by the permanent charges inside biomolecules, and formulated regular, well-posed PNP equations. An inexact-Newton method was used to solve the coupled nonlinear elliptic equations for the steady problems; while an Adams-Bashforth-Crank-Nicolson method was devised for time integration for the unsteady electrodiffusion. We numerically investigated the conditioning of the stiffness matrices for the finite element approximations of the two formulations of the Nernst-Planck equation, and theoretically proved that the transformed formulation is always associated with an ill-conditioned stiffness matrix. We also studied the electroneutrality of the solution and its relation with the boundary conditions on the molecular surface, and concluded that a large net charge concentration is always present near the molecular surface due to the presence of multiple species of charged particles in the solution. The numerical methods are shown to be accurate and stable by various test problems, and are applicable to real large-scale biophysical electrodiffusion problems.

Finite element analysis of drug electrostatic diffusion: inhibition rate studies in N1 neuraminidase

This article describes a numerical solution of the steady-state Poisson-Boltzmann-Smoluchowski (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.

The spatiotemporal pattern of Src activation at lipid rafts revealed by diffusion-corrected FRET imaging

Genetically encoded biosensors based on fluorescence resonance energy transfer (FRET) have been widely applied to visualize the molecular activity in live cells with high spatiotemporal resolution. However, the rapid diffusion of biosensor proteins hinders a precise reconstruction of the actual molecular activation map. Based on fluorescence recovery after photobleaching (FRAP) experiments, we have developed a finite element (FE) method to analyze, simulate, and subtract the diffusion effect of mobile biosensors. This method has been applied to analyze the mobility of Src FRET biosensors engineered to reside at different subcompartments in live cells. The results indicate that the Src biosensor located in the cytoplasm moves 4–8 folds faster (0.93±0.06 µm²/sec) than those anchored on different compartments in plasma membrane (at lipid raft: 0.11±0.01 µm²/sec and outside: 0.18±0.02 µm²/sec). The mobility of biosensor at lipid rafts is slower than that outside of lipid rafts and is dominated by two-dimensional diffusion. When this diffusion effect was subtracted from the FRET ratio images, high Src activity at lipid rafts was observed at clustered regions proximal to the cell periphery, which remained relatively stationary upon epidermal growth factor (EGF) stimulation. This result suggests that EGF induced a Src activation at lipid rafts with well-coordinated spatiotemporal patterns. Our FE-based method also provides an integrated platform of image analysis for studying molecular mobility and reconstructing the spatiotemporal activation maps of signaling molecules in live cells.

Fluorescence biosensors have been widely used to report the spatial and temporal activity of target molecules in live cells. However, biosensors can move independently of the target molecule and carry its signal to other subcellular locations. Therefore, the observed images appear to be the combination of the target molecular activity and the artifacts introduced by the movement of the biosensors (mainly due to diffusion). The intriguing question is how to estimate and exclude the movement effect of biosensors from the observed fluorescent images and to reconstruct the real activity map of the target molecules. The Src molecule plays important roles in cell adhesion, migration, and cancer invasion. In this paper, we developed a novel computational method to analyze and simulate the movement of the Src biosensor, which was then subtracted from the original fluorescent images. With this computational method, we observed discrete clusters of high Src activity at relatively stationary locations on the plasma membrane. Therefore, our results highlight the coordination of molecular activities in space and time. In addition to Src, our computational method can be used to reconstruct the activity map of other signaling molecules.

Tortuosity and anomalous diffusion in the neuromuscular junction

The signal transfer from nerve to muscle occurs by diffusion across the neuromuscular junction. The continuum level analysis of diffusion processes is based on the diffusion equation, which in one dimension is partial differential c/partial differential t=D(partial differential(2)c/partial differential x(2)) , where c is the molecular concentration and D is the diffusivity. However, in confined systems such as the neuromuscular junction, the diffusion equation may not be valid, and even if valid the value of D may be altered by the confinement. In this paper, Monte Carlo simulations are used to probe diffusion at the molecular level in a realistic model of a neuromuscular junction. The results show that diffusion is anomalous (i.e., not described by the diffusion equation) for time scales less than approximately 0.01 s, which is the time scale relevant for signaling processes in the synapse. At longer time scales, the diffusion is normal (i.e., described by the diffusion equation), but with a value of D that is reduced by a factor of approximately 5 times compared to the value for diffusion in open space. As the width of the synaptic cleft decreases, these effects become even more pronounced. The physical basis of these results is described in terms of the structure of the neuromuscular junction.

Continuum simulations of acetylcholine consumption by acetylcholinesterase: a Poisson-Nernst-Planck approach

The Poisson-Nernst-Planck (PNP) equation provides a continuum description of electrostatic-driven diffusion and is used here to model the diffusion and reaction of acetylcholine (ACh) with acetylcholinesterase (AChE) enzymes. This study focuses on the effects of ion and substrate concentrations on the reaction rate and rate coefficient. To this end, the PNP equations are numerically solved with a hybrid finite element and boundary element method at a wide range of ion and substrate concentrations, and the results are compared with the partially coupled Smoluchowski-Poisson-Boltzmann model. The reaction rate is found to depend strongly on the concentrations of both the substrate and ions; this is explained by the competition between the intersubstrate repulsion and the ionic screening effects. The reaction rate coefficient is independent of the substrate concentration only at very high ion concentrations, whereas at low ion concentrations the behavior of the rate depends strongly on the substrate concentration. Moreover, at physiological ion concentrations, variations in substrate concentration significantly affect the transient behavior of the reaction. Our results offer a reliable estimate of reaction rates at various conditions and imply that the concentrations of charged substrates must be coupled with the electrostatic computation to provide a more realistic description of neurotransmission and other electrodiffusion and reaction processes.

Electrodiffusion: a continuum modeling framework for biomolecular systems with realistic spatiotemporal resolution

A computational framework is presented for the continuum modeling of cellular biomolecular diffusion influenced by electrostatic driving forces. This framework is developed from a combination of state-of-the-art numerical methods, geometric meshing, and computer visualization tools. In particular, a hybrid of (adaptive) finite element and boundary element methods is adopted to solve the Smoluchowski equation (SE), the Poisson equation (PE), and the Poisson-Nernst-Planck equation (PNPE) in order to describe electrodiffusion processes. The finite element method is used because of its flexibility in modeling irregular geometries and complex boundary conditions. The boundary element method is used due to the convenience of treating the singularities in the source charge distribution and its accurate solution to electrostatic problems on molecular boundaries. Nonsteady-state diffusion can be studied using this framework, with the electric field computed using the densities of charged small molecules and mobile ions in the solvent. A solution for mesh generation for biomolecular systems is supplied, which is an essential component for the finite element and boundary element computations. The uncoupled Smoluchowski equation and Poisson-Boltzmann equation are considered as special cases of the PNPE in the numerical algorithm, and therefore can be solved in this framework as well. Two types of computations are reported in the results: stationary PNPE and time-dependent SE or Nernst-Planck equations solutions. A biological application of the first type is the ionic density distribution around a fragment of DNA determined by the equilibrium PNPE. The stationary PNPE with nonzero flux is also studied for a simple model system, and leads to an observation that the interference on electrostatic field of the substrate charges strongly affects the reaction rate coefficient. The second is a time-dependent diffusion process: the consumption of the neurotransmitter acetylcholine by acetylcholinesterase, determined by the SE and a single uncoupled solution of the Poisson-Boltzmann equation. The electrostatic effects, counterion compensation, spatiotemporal distribution, and diffusion-controlled reaction kinetics are analyzed and different methods are compared.

The effects of tertiapin-Q on responses of the sinoatrial pacemaker of the guinea-pig heart to vagal nerve stimulation and muscarinic agonists

Using Langendorff preparations of the guinea-pig heart, we have examined the participation of the acetylcholine (ACh)-activated potassium channel, IK,ACh, in the bradycardia produced by electrical stimulation of the vagus (parasympathetic) nerve and muscarinic agonists (ACh and bethanecol, bolus i.a.). Hearts from young animals (160-250 g) were perfused with Krebs-Henseleit solution, and pacemaker frequency was determined from the P wave of an ECG. Tertiapin-Q was used to block IK,ACh. Vagal stimulation (10 s trains at 2, 5 and 10 Hz) produced graded reductions in atrial rate that were substantially attenuated, and to a similar extent, by 300 nm and 1 microm tertiapin-Q (to 0.42 +/- 0.12, mean +/- s.d., of the control values; P < 0.001). Acetylcholine (3 nmol) produced brief graded bradycardias that were also attenuated by tertiapin-Q (0.24 +/- 0.24; P = 0.006). Similar results were obtained when experiments were repeated in 2 mm Cs+ (to block the hyperpolarization-activated pacemaker current). Bethanecol (30, 50 and 70 nmol), a muscarinic agonist with no appreciable nicotinic activity, produced sustained bradycardias that were attenuated by 300 nm tertiapin-Q (0.36 +/- 0.21; P < 0.0001). The responses to vagal stimulation and ACh developed more slowly in tertiapin-Q, indicating that a rapidly acting mechanism had been blocked. Responses to vagal stimulation were faster in 2 mm Cs+. Together, these observations show that ACh released from parasympathetic nerve varicosities exerts a considerable part of its effect on the pacemaker by activating IK,ACh and acts in a manner not readily distinguishable from that of directly applied muscarinic agonists.

Finite element analysis of the time-dependent Smoluchowski equation for acetylcholinesterase reaction rate calculations

This article describes the numerical solution of the time-dependent Smoluchowski equation to study diffusion in biomolecular systems. Specifically, finite element methods have been developed to calculate ligand binding rate constants for large biomolecules. The resulting software has been validated and applied to the mouse acetylcholinesterase (mAChE) monomer and several tetramers. Rates for inhibitor binding to mAChE were calculated at various ionic strengths with several different time steps. Calculated rates show very good agreement with experimental and theoretical steady-state 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 rate accelerations of the monomeric and tetrameric mAChE that result from electrostatic steering are preserved under the non-steady-state conditions that are expected to occur in physiological circumstances.

Glutamate, water and ion transport through a charged nanosize pore

The transport of transmitter, ions and water through a positively-charged nanopore was investigated through computer simulations. The physics of the problem is described by a coupled set of Poisson-Nernst-Planck and Navier-Stokes equations in a computational domain consisting a cylindrical pore, whose radius ranged from 1 to 8 nm and which was flanked by two compartments representing the vesicular interior and extra-cellular space. The concentration of co-ions is suppressed and of counter-ions enhanced, especially near the pore wall owing to electrostatic interactions. Glutamate (i.e. the transmitter considered) is negatively charged and is simulated as a counter-ion. The electro-kinetically induced pressure due to the movement of ions is negative and very pronounced near the pore wall where the concentration and flux of counter-ions is very high. The water velocity peaks in the pore center, diminishes to zero at the pore wall, but is constant along the pore axis. The mean velocity of the water/fluid is proportional to the vesicular pressure and pore cross-sectional area. Interestingly it is inversely related to the vesicular glutamate concentration. The factors determining the glutamate flux are complex. The diffusive flux generally predominates for narrow pore, and convective flux may dominate for wide pore if the vesicular pressure is high. Surprisingly at low vesicular pressure the mean total glutamate flux per unit cross-sectional pore area is higher for narrow pores. Higher flux is probably due to the rise of glutamate concentration in the nanopore, which is much more pronounced for narrow nanopores, due to the maintenance of approximate neutrality of charges in the pore and on the pore wall. In conclusion intra-vesicular pressure helps 'flushing-out' the transmitter, but the induced pressure 'drags-out' the water into the extra-cellular space.

Continuum simulations of acetylcholine diffusion with reaction-determined boundaries in neuromuscular junction models

The reaction-diffusion system of the neuromuscular junction has been modeled in 3D using the finite element package FEtk. The numerical solution of the dynamics of acetylcholine with the detailed reaction processes of acetylcholinesterases and nicotinic acetylcholine receptors has been discussed with the reaction-determined boundary conditions. The simulation results describe the detailed acetylcholine hydrolysis process, and reveal the time-dependent interconversion of the closed and open states of the acetylcholine receptors as well as the percentages of unliganded/monoliganded/diliganded states during the neuro-transmission. The finite element method has demonstrated its flexibility and robustness in modeling large biological systems.

Finite concentration effects on diffusion-controlled reactions

The algorithm by Northrup, Allison, and McCammon [J. Chem. Phys. 80, 1517 (1984)] has been used for two decades for calculating the diffusion-influenced rate-constants of enzymatic reactions. Although many interesting results have been obtained, the algorithm is based on the assumption that substrate-substrate interactions can be neglected. This approximation may not be valid when the concentration of the ligand is high. In this work, we constructed a simulation model that can take substrate-substrate interactions into account. We first validated the model by carrying out simulations in ways that could be compared to analytical theories. We then carried out simulations to examine the possible effects of substrate-substrate interactions on diffusion-controlled reaction rates. For a substrate concentration of 0.1 mM, we found that the diffusion-controlled reaction rates were not sensitive to whether substrate-substrate interactions were included. On the other hand, we observed significant influence of substrate-substrate interactions on calculated reaction rates at a substrate concentration of 0.1M. Therefore, a simulation model that takes substrate-substrate interactions into account is essential for reliably predicting diffusion-controlled reaction rates at high substrate concentrations, and one such simulation model is presented here.

Spatial modeling of dimerization reaction dynamics in the plasma membrane: Monte Carlo vs. continuum differential equations

Bimolecular reactions in the plasma membrane, such as receptor dimerization, are a key signaling step for many signaling systems. For receptors to dimerize, they must first diffuse until a collision happens, upon which a dimerization reaction may occur. Therefore, study of the dynamics of cell signaling on the membrane may require the use of a spatial modeling framework. Despite the availability of spatial simulation methods, e.g., stochastic spatial Monte Carlo (MC) simulation and partial differential equation (PDE) based approaches, many biological models invoke well-mixed assumptions without completely evaluating the importance of spatial organization. Whether one is to utilize a spatial or non-spatial simulation framework is therefore an important decision. In order to evaluate the importance of spatial effects a priori, i.e., without performing simulations, we have assessed the applicability of a dimensionless number, known as second Damköhler number (Da), defined here as the ratio of time scales of collision and reaction, for 2-dimensional bimolecular reactions. Our study shows that dimerization reactions in the plasma membrane with Da approximately >0.1 (tested in the receptor density range of 10(2)-10(5)/microm(2)) require spatial modeling. We also evaluated the effective reaction rate constants of MC and simple deterministic PDEs. Our simulations show that the effective reaction rate constant decreases with time due to time dependent changes in the spatial distribution of receptors. As a result, the effective reaction rate constant of simple PDEs can differ from that of MC by up to two orders of magnitude. Furthermore, we show that the fluctuations in the number of copies of signaling proteins (noise) may also depend on the diffusion properties of the system. Finally, we used the spatial MC model to explore the effect of plasma membrane heterogeneities, such as receptor localization and reduced diffusivity, on the dimerization rate. Interestingly, our simulations show that localization of epidermal growth factor receptor (EGFR) can cause the diffusion limited dimerization rate to be up to two orders of magnitude higher at higher average receptor densities reported for cancer cells, as compared to a normal cell.

Continuum diffusion reaction rate calculations of wild-type and mutant mouse acetylcholinesterase: adaptive finite element analysis

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.

Finite element solution of the steady-state Smoluchowski equation for rate constant calculations

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.

Interfacing Brownian dynamics simulations

Starting from the flux of particles in a Brownian dynamics simulation we derive boundary conditions, which allow us (i) to couple a Brownian dynamics calculation to a reservoir of particles of a given density, i.e., setting up constant density boundary conditions, and (ii) to build an interface between Brownian dynamics and a diffusional treatment of adjacent simulation volumes. With these algorithms it is sometimes possible to dramatically reduce the system size--and therefore the necessary resources--of multiparticle Brownian dynamics calculations. In this paper we give one-dimensional examples which illustrate potential applications and savings.