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.     

Far-field analysis of coupled bulk and boundary layer diffusion toward an ion channel entrance.

We present a far-field analysis of ion diffusion toward a channel embedded in a membrane with a fixed charge density. The Smoluchowski equation, which represents the 3D problem, is approximated by a system of coupled three- and two-dimensional diffusions. The 2D diffusion models the quasi-two-dimensional diffusion of ions in a boundary layer in which the electrical potential interaction with the membrane surface charge is important. The 3D diffusion models ion transport in the bulk region outside the boundary layer. Analytical expressions for concentration and flux are developed that are accurate far from the channel entrance. These provide boundary conditions for a numerical solution of the problem. Our results are used to calculate far-field ion flows corresponding to experiments of Bell and Miller (Biophys. J. 45:279, 1984).     

A Comparison of Bimolecular Reaction Models for Stochastic Reaction–Diffusion Systems

Stochastic reaction–diffusion models have become an important tool in studying how both noise in the chemical reaction process and the spatial movement of molecules influences the behavior of biological systems. There are two primary spatially-continuous models that have been used in recent studies: the diffusion limited reaction model of Smoluchowski, and a second approach popularized by Doi. Both models treat molecules as points undergoing Brownian motion. The former represents chemical reactions between two reactants through the use of reactive boundary conditions, with two molecules reacting instantly upon reaching a fixed separation (called the reaction-radius). The Doi model uses reaction potentials, whereby two molecules react with a fixed probability per unit time, λ, when separated by less than the reaction radius. In this work, we study the rigorous relationship between the two models. For the special case of a protein diffusing to a fixed DNA binding site, we prove that the solution to the Doi model converges to the solution of the Smoluchowski model as λ→∞, with a rigorous $O(\lambda^{-\frac{1}{2} + \epsilon})$ error bound (for any fixed ?>0). We investigate by numerical simulation, for biologically relevant parameter values, the difference between the solutions and associated reaction time statistics of the two models. As the reaction-radius is decreased, for sufficiently large but fixed values of λ, these differences are found to increase like the inverse of the binding radius.     

Diffusion of oxygen in water and hydrocarbons using an electron spin resonance spin-label technique.

The Smoluchowski equation for the bimolecular collision rate of dissolved oxygen molecules with spin labels yielded values for the diffusion constant of oxygen in water that are in agreement with the Stokes-Einstein equation (D infinity T/eta, where eta is the macroscopic viscosity) and with published values obtained by conventional methods. Heisenberg exchange at an interaction distance of 4.5 A occurs with a probability close to one for each encounter. In mixed hydrocarbons (olive oil, paraffin oils) and sec-butyl benzene, D infinity (T/eta)rho, where rho lies between 0.5 and 1. Oxygen diffuses in the hydrocarbons between 10 and 100 times more rapidly than predicted from the macroscopic viscosity. Similar results would be expected for diffusion of oxygen in model and biological membranes. Parallel measurements of rotational diffusion of the spin labels show little correlation with measurements of translational diffusion of oxygen. Dipolar interactions between spin labels and oxygen appear negligible except in the limit of highest viscosities.     

Local models of cell aggregation kinetics

Mathematical models of the phenomena of the sorting out of embryonic cells are derived which predict multiple clustering of like cell types into many small islands. The experiments of Trinkaus and Lentz, and other essentially two-dimensional cell sorting results, provide examples. Two continuous models are introduced, both based on the Smoluchowski theory of Brownian motion. A mass-interaction model predicts local clumping of the internally segregating cell type within a continuum of the other type. A diffusion model simulates the suspension-aggregate experiments of Roth and Weston. A discrete model of sorting out, the exchange model, is re-derived in a clearer manner than before, and certain recent criticisms are responded to. Results in recent two-dimensional cell sorting experiments are discussed and a stability criterion for patterns of multiple clusters is introduced. The site frequency model is re-examined and criticized in terms of its application to the Roth-Weston experiments.     

Diffusional channeling in the sulfate-activating complex: combined continuum modeling and coarse-grained brownian dynamics studies

Enzymes required for sulfur metabolism have been suggested to gain efficiency by restricted diffusion (i.e., channeling) of an intermediate APS^2- between active sites. This article describes modeling of the whole channeling process by numerical solution of the Smoluchowski diffusion equation, as well as by coarse-grained Brownian dynamics. The results suggest that electrostatics plays an essential role in the APS^2- channeling. Furthermore, with coarse-grained Brownian dynamics, the substrate channeling process has been studied with reactions in multiple active sites. Our simulations provide a bridge for numerical modeling with Brownian dynamics to simulate the complicated reaction and diffusion and raise important questions relating to the electrostatically mediated substrate channeling in vitro, in situ, and in vivo.     

A new approximate whole boundary solution of the Lamm differential equation for the analysis of sedimentation velocity experiments

Sedimentation velocity is one of the best-suited physical methods for determining the size and shape of macromolecular substances or their complexes in the range from 1 to several thousand kDa. The moving boundary in sedimentation velocity runs can be described by the Lamm differential equation. Fitting of suitable model functions or solutions of the Lamm equation to the moving boundary is used to obtain directly sedimentation and diffusion coefficients, thus allowing quick determination of size, shape and other parameters of macromolecules. Here we present a new approximate whole boundary solution of the Lamm equation that simultaneously allows the specification of sedimentation and diffusion coefficients with deviations smaller than 1% from the expected values.     

Inverse diffusion methods for data peak separation

Previous methods for separation of overlapping data peaks include geometrical assessment and Fourier deconvolution. On the basis of inverse diffusion theory, we present new separation methods suitable for convenient programming and rapid calculation of Gaussian area contributions. Both continuum and discrete inverse diffusion models are described. Example computations are given for biological data: density gradient centrifugation, isoelectric focusing electrophoresis, and high-pressure liquid chromatography.     

Effects of Hydraulic Soil Properties on Vegetation Pattern Formation in Sloping Landscapes

Current models of vegetation pattern formation rely on a system of weakly nonlinear reaction–diffusion equations that are coupled by their source terms. While these equations, which are used to describe a spatiotemporal planar evolution of biomass and soil water, qualitatively capture the emergence of various types of vegetation patterns in arid environments, they are phenomenological and have a limited predictive power. We ameliorate these limitations by deriving the vertically averaged Richards’ equation to describe flow (as opposed to “diffusion”) of water in partially saturated soils. This establishes conditions under which this nonlinear equation reduces to its weakly nonlinear reaction–diffusion counterpart used in the previous models, thus relating their unphysical parameters (e.g., diffusion coefficient) to the measurable soil properties (e.g., hydraulic conductivity) used to parameterize the Richards equation. Our model is valid for both flat and sloping landscapes and can handle arbitrary topography and boundary conditions. The result is a model that relates the environmental conditions (e.g., precipitation rate, runoff and soil properties) to formation of multiple patterns observed in nature (such as stripes, labyrinth and spots).     

On the Protein Crystal Formation as an Interface-Controlled Process with Prototype Ion-Channeling Effect

A superdiffusive random-walk action in the depletion zone around a growing protein crystal is considered. It stands for a dynamic boundary condition of the growth process and competes steadily with a quasistatic, curvature-involving (thermodynamic) free boundary condition, both of them contributing to interpret the (mainly late-stage) growth process in terms of a prototype ion-channeling effect. An overall diffusion function contains quantitative signatures of both boundary conditions mentioned and indicates whether the new phase grows as an orderly phase or a converse scenario occurs. This situation can be treated in a quite versatile way both numerically and analytically, within a generalized Smoluchowski framework. This study can help in (1) elucidating some dynamic puzzles of a complex crystal formation vs biomolecular aggregation, also those concerning ion-channel formation, and (2) seeing how ion-channel-type dynamics of non-Markovian nature may set properly the pace of model (dis)ordered protein aggregation.     

Observations on rate theory for rugged energy landscapes

The potential energy profile for many complex reactions of proteins, such as folding or allosteric conformational change, involves many different scales of molecular motion along the reaction coordinate. Although it is natural to model the dynamics of motion along such rugged energy landscapes as diffusional (the Smoluchowski equation; SE), problems arise because the frictional forces generated by the molecular surround are typically not strong enough to justify the use of the SE. Here, we discuss the fundamental theory behind the SE and note that it may be justified through a master equation when reduced to its continuum limit. However, the SE cannot be used for rough energy landscapes, where the continuum limit is ill defined. Instead, we suggest that one should use a mean first passage time expression derived from a master equation, and show how this approach can be used to glean information about the underlying dynamics of barrier crossing. We note that the potential profile in the SE is that of the microbarriers between conformational substates, and that there is a temperature-dependent, effective friction associated with the long residence time in the microwells that populate the rough landscape. The number of recrossings of the overall barrier is temperature-dependent, governed by the microbarriers and not by the effective friction. We derive an explicit expression for the mean number of recrossings and its temperature dependence. Finally, we note that the mean first passage time can be used as a departure point for measuring the roughness of the landscape.     

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.     

Proton transport across charged membrane and pH oscillations.

Based on Eyring's multibarrier activation process, a mathematical model and equation is developed to account for proton diffusion through an immobilized protein and enzyme membrane perfused with an electrolyte, substrate, and a buffer. With this model we find that, in the presence of a buffer, our solution approaches the continuum case very rapidly. We apply our model to membranes composed of papain and bovine serum albumin and find that our theory closely stimulates the experimental observations on the effect of salt and buffer on proton diffusion. Our theory shows that the pH oscillations observed in the diffusion controlled papain-benzoyl-L-arginine ethyl ester (BAEE) reaction may be the result of CO2 dissolved in the bath at high pH. In our theory, under certain conditions and in agreement with experimental observation, the buffer penetration depth oscillates near the boundary of a papain membrane in a solution containing BAEE and borate. We also find that at low ionic strength small ions as well as a buffer are seen to oscillate if a membrane is highly charged.     

Rate theory models for ion transport through rigid pores. I. Time-dependent analysis in the case of vanishing interactions

In this first of a series of papers concerning the theoretical analysis of rate theory models for ion transport through rigid pores, the case of vanishing interactions is investigated. "Rigidity" means that ions crossing membranes through pores see a fixed structure of the pores, not changing in time. A single pore is considered to be a sequence of (n + 1) activation barriers separated by n energy minima. The explicit analytical treatment is restricted to pores with regular internal barrier structure, including the nonequilibrium situation of an applied electric field. In this case the connection with continuum diffusion models is demonstrated by performing in the limit n leads to infinity (n = number of binding sites within the pores) the transition to continuum. Thus, from diffusion equations describing a discrete number of jumps, the corresponding diffusion-like partial differential equations and boundary conditions are generated. For regular pores, from the time dependent solutions of the discrete equations, the corresponding solutions of the continuum equations are explicitly generated. The time-dependent relaxation behaviour of the discrete model is in good agreement with the continuum model if one assumes more than two binding sites in the pores.     

Application of the Smoluchowski equation with a source term to the model of lipid pore formation during a phase transition

The Smoluchowski equation, which describes pore diffusion in the radius space, with a source term, is used in modeling the process of the formation of a hydrophilic pore in a lipid bilayer during phase transition. The introduction of a hydrophobic-pore source term into the equation reflects the emergence of additional defects in a bilayer caused by the decrease in the molecule area under the transition from the liquid crystalline to the gel phase. The distribution of the time probability density calculated within the model that is required for the formation of a hydrophilic pore is in good agreement with the previously published experimental data.     

Diffusion-controlled sensitization of photocleavage reactions on surfaces

The kinetic rate equation for the photosensitized cleavage reaction of surface-bound photolabile chromophores with free diffusion of sensitizer molecules from the bulk of a solution to the surface is derived by determining the stationary solution of a diffusion equation with suitable boundary conditions. The relation between the phenomenological rate constant for the photosensitized reaction at the surface and in the bulk is established. Applying the result to the analysis of an experimental example, the origin of the quasi zeroth-order kinetics of the sensitized reaction is revealed. A theoretical comparison of intramolecular sensitization in photocleavable protecting groups with a molecular antenna and sensitization with the freely diffusing sensitizer shows that in a typical case sensitization with free diffusion is more effective than intramolecular sensitization for sensitizer concentrations higher than 5 mM.     

Exclusion processes on a growing domain

A discrete model provides a useful framework for experimentalists to understand the interactions between growing tissues and other biological mechanisms. A cellular automata (CA) model with domain growth, cell motility and cell proliferation, based on cellular exclusion processes, is developed here. Average densities can be defined from the CA model and a continuum representation can be determined. The domain growth mechanism in the CA model gives rise to a Fokker-Planck equation in the corresponding continuum model, with a diffusive and a convective term. Deterministic continuum models derived from conservation laws do not include this diffusive term. The new diffusive term arises because of the stochasticity inherited from the CA mechanism for domain growth. We extend the models to multiple species and investigate the influence of the flux terms arising from the exclusion processes. The averaged CA agent densities are well approximated by the solution of nonlinear advection-diffusion equations, provided that the relative size of the proliferation processes to the diffusion processes is sufficiently small. This dual approach provides an understanding of the microscopic and macroscopic scales in a developmental process.     

Analysis of synaptic transmission in the neuromuscular junction using a continuum finite element model.

There is a steadily growing body of experimental data describing the diffusion of acetylcholine in the neuromuscular junction and the subsequent miniature endplate currents produced at the postsynaptic membrane. To gain further insights into the structural features governing synaptic transmission, we have performed calculations using a simplified finite element model of the neuromuscular junction. The diffusing acetylcholine molecules are modeled as a continuum, whose spatial and temporal distribution is governed by the force-free diffusion equation. The finite element method was adopted because of its flexibility in modeling irregular geometries and complex boundary conditions. The resulting simulations are shown to be in accord with experiment and other simulations.     

Analysis of cell growth in a polymer scaffold using a moving boundary approach

Two mathematical models of chondrocyte generation and nutrient consumption are developed to analyze the behavior of cell growth in a biodegradable polymer matrix. Substrate reaction and diffusion are analyzed in two regions: one consisting of cells and nutrients and the other consisting of only nutrients. A pseudo-steady state approximation for the transport of nutrients in these two regions is utilized. The rate of growth is determined by a moving boundary equation that equates the rate at which the interfacial region between the cells and the void space moves to a substrate dependent growth reaction. The change in the location of this interfacial region with time therefore depicts the rate at which the cells propagate. The two limiting cases discussed in this article represent extremes in how the cells will grow in the polymer matrix; one case assumes that cells grow inward from the external boundary, and the other case assumes that cells grow parallel to the external boundary. The results of both models are compared to experimental data found in the literature. It is found through these comparisons that the model parameters, including the unit cell spacing parameter L, the metabolic rate constant k, the growth rate constant k(G), and external mass transfer coefficient, K, may vary as the thickness of the polymer matrix is changed, however, unrealistic and large changes in the diffusion coefficients were required to account for the full range of experimental data. Furthermore, these results suggest modification of the functional form of the growth kinetics to include substrate or product inhibition, or death terms. Based upon diffusion/reaction concepts, these models for cell growth in a biodegradable polymer give bounds for the upper and lower limits of the cellular growth rate and nutrient consumption in a polymer matrix and will aid in the development of more extensive models. (c) 1997 John Wiley & Sons, Inc. Biotechnol Bioeng 56: 422-432, 1997.