The Implementation of the Colored Abstract Simplicial Complex and its Application to Mesh Generation.

We introduce CASC: a new, modern, and header-only C++ library which provides a data structure to represent arbitrary dimension abstract simplicial complexes (ASC) with user-defined classes stored directly on the simplices at each dimension. This is accomplished by using the latest C++ language features including variadic template parameters introduced in C++11 and automatic function return type deduction from C++14. Effectively CASC decouples the representation of the topology from the interactions of user data. We present the innovations and design principles of the data structure and related algorithms. This includes a metadata aware decimation algorithm which is general for collapsing simplices of any dimension. We also present an example application of this library to represent an orientable surface mesh.

Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit.

Triggered release of Ca2+ from an individual sarcoplasmic reticulum (SR) Ca(2+) release unit (CRU) is the fundamental event of cardiac excitation­contraction coupling, and spontaneous release events (sparks) are the major contributor to diastolic Ca(2+) leak in cardiomyocytes. Previous model studies have predicted that the duration and magnitude of the spark is determined by the local CRU geometry, as well as the localization and density of Ca(2+) handling proteins. We have created a detailed computational model of a CRU, and developed novel tools to generate the computational geometry from electron tomographic images. Ca(2+) diffusion was modelled within the SR and the cytosol to examine the effects of localization and density of the Na(+)/Ca(2+) exchanger, sarco/endoplasmic reticulum Ca(2+)-ATPase 2 (SERCA), and calsequestrin on spark dynamics. We reconcile previous model predictions of approximately 90% local Ca(2+) depletion in junctional SR, with experimental reports of about 40%. This analysis supports the hypothesis that dye kinetics and optical averaging effects can have a significant impact on measures of spark dynamics. Our model also predicts that distributing calsequestrin within non-junctional Z-disc SR compartments, in addition to the junctional compartment, prolongs spark release time as reported by Fluo5. By pumping Ca(2+) back into the SR during a release, SERCA is able to prolong a Ca(2+) spark, and this may contribute to SERCA-dependent changes in Ca(2+) wave speed. Finally, we show that including the Na(+)/Ca(2+) exchanger inside the dyadic cleft does not alter local [Ca(2+)] during a spark.

Numerical analysis of Ca2+ signaling in rat ventricular myocytes with realistic transverse-axial tubular geometry and inhibited sarcoplasmic reticulum.

The t-tubules of mammalian ventricular myocytes are invaginations of the cell membrane that occur at each Z-line. These invaginations branch within the cell to form a complex network that allows rapid propagation of the electrical signal, and hence synchronous rise of intracellular calcium (Ca(2+)). To investigate how the t-tubule microanatomy and the distribution of membrane Ca(2+) flux affect cardiac excitation-contraction coupling we developed a 3-D continuum model of Ca(2+) signaling, buffering and diffusion in rat ventricular myocytes. The transverse-axial t-tubule geometry was derived from light microscopy structural data. To solve the nonlinear reaction-diffusion system we extended SMOL software tool (http://mccammon.ucsd.edu/smol/). The analysis suggests that the quantitative understanding of the Ca(2+) signaling requires more accurate knowledge of the t-tubule ultra-structure and Ca(2+) flux distribution along the sarcolemma. The results reveal the important role for mobile and stationary Ca(2+) buffers, including the Ca(2+) indicator dye. In agreement with experiment, in the presence of fluorescence dye and inhibited sarcoplasmic reticulum, the lack of detectible differences in the depolarization-evoked Ca(2+) transients was found when the Ca(2+) flux was heterogeneously distributed along the sarcolemma. In the absence of fluorescence dye, strongly non-uniform Ca(2+) signals are predicted. Even at modest elevation of Ca(2+), reached during Ca(2+) influx, large and steep Ca(2+) gradients are found in the narrow sub-sarcolemmal space. The model predicts that the branched t-tubule structure and changes in the normal Ca(2+) flux density along the cell membrane support initiation and propagation of Ca(2+) waves in rat myocytes.

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.

Three-dimensional geometric modeling of membrane-bound organelles in ventricular myocytes: bridging the gap between microscopic imaging and mathematical simulation.

A general framework of image-based geometric processing is presented to bridge the gap between three-dimensional (3D) imaging that provides structural details of a biological system and mathematical simulation where high-quality surface or volumetric meshes are required. A 3D density map is processed in the order of image pre-processing (contrast enhancement and anisotropic filtering), feature extraction (boundary segmentation and skeletonization), and high-quality and realistic surface (triangular) and volumetric (tetrahedral) mesh generation. While the tool-chain described is applicable to general types of 3D imaging data, the performance is demonstrated specifically on membrane-bound organelles in ventricular myocytes that are imaged and reconstructed with electron microscopic (EM) tomography and two-photon microscopy (T-PM). Of particular interest in this study are two types of membrane-bound Ca(2+)-handling organelles, namely, transverse tubules (T-tubules) and junctional sarcoplasmic reticulum (jSR), both of which play an important role in regulating the excitation-contraction (E-C) coupling through dynamic Ca(2+) mobilization in cardiomyocytes.

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.

Feature-preserving adaptive mesh generation for molecular shape modeling and simulation.

We describe a chain of algorithms for molecular surface and volumetric mesh generation. We take as inputs the centers and radii of all atoms of a molecule and the toolchain outputs both triangular and tetrahedral meshes that can be used for molecular shape modeling and simulation. Experiments on a number of molecules are demonstrated, showing that our methods possess several desirable properties: feature-preservation, local adaptivity, high quality, and smoothness (for surface meshes). We also demonstrate an example of molecular simulation using the finite element method and the meshes generated by our method. The approaches presented and their implementations are also applicable to other types of inputs such as 3D scalar volumes and triangular surface meshes with low quality, and hence can be used for generation/improvement of meshes in a broad range of applications.

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 pKa values, and an improved web-based visualization tool for viewing electrostatics.

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.

High-order finite element methods for cardiac monodomain simulations.

Computational modeling of tissue-scale cardiac electrophysiology requires numerically converged solutions to avoid spurious artifacts. The steep gradients inherent to cardiac action potential propagation necessitate fine spatial scales and therefore a substantial computational burden. The use of high-order interpolation methods has previously been proposed for these simulations due to their theoretical convergence advantage. In this study, we compare the convergence behavior of linear Lagrange, cubic Hermite, and the newly proposed cubic Hermite-style serendipity interpolation methods for finite element simulations of the cardiac monodomain equation. The high-order methods reach converged solutions with fewer degrees of freedom and longer element edge lengths than traditional linear elements. Additionally, we propose a dimensionless number, the cell Thiele modulus, as a more useful metric for determining solution convergence than element size alone. Finally, we use the cell Thiele modulus to examine convergence criteria for obtaining clinically useful activation patterns for applications such as patient-specific modeling where the total activation time is known a priori.

Modeling effects of L-type ca(2+) current and na(+)-ca(2+) exchanger on ca(2+) trigger flux in rabbit myocytes with realistic T-tubule geometries.

The transverse tubular system of rabbit ventricular myocytes consists of cell membrane invaginations (t-tubules) that are essential for efficient cardiac excitation-contraction coupling. In this study, we investigate how t-tubule micro-anatomy, L-type Ca(2+) channel (LCC) clustering, and allosteric activation of Na(+)/Ca(2+) exchanger by L-type Ca(2+) current affects intracellular Ca(2+) dynamics. Our model includes a realistic 3D geometry of a single t-tubule and its surrounding half-sarcomeres for rabbit ventricular myocytes. The effects of spatially distributed membrane ion-transporters (LCC, Na(+)/Ca(2+) exchanger, sarcolemmal Ca(2+) pump, and sarcolemmal Ca(2+) leak), and stationary and mobile Ca(2+) buffers (troponin C, ATP, calmodulin, and Fluo-3) are also considered. We used a coupled reaction-diffusion system to describe the spatio-temporal concentration profiles of free and buffered intracellular Ca(2+). We obtained parameters from voltage-clamp protocols of L-type Ca(2+) current and line-scan recordings of Ca(2+) concentration profiles in rabbit cells, in which the sarcoplasmic reticulum is disabled. Our model results agree with experimental measurements of global Ca(2+) transient in myocytes loaded with 50 µM Fluo-3. We found that local Ca(2+) concentrations within the cytosol and sub-sarcolemma, as well as the local trigger fluxes of Ca(2+) crossing the cell membrane, are sensitive to details of t-tubule micro-structure and membrane Ca(2+) flux distribution. The model additionally predicts that local Ca(2+) trigger fluxes are at least threefold to eightfold higher than the whole-cell Ca(2+) trigger flux. We found also that the activation of allosteric Ca(2+)-binding sites on the Na(+)/Ca(2+) exchanger could provide a mechanism for regulating global and local Ca(2+) trigger fluxes in vivo. Our studies indicate that improved structural and functional models could improve our understanding of the contributions of L-type and Na(+)/Ca(2+) exchanger fluxes to intracellular Ca(2+) dynamics.

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.

Three-dimensional electron microscopy reveals new details of membrane systems for Ca2+ signaling in the heart.

In the current study, the three-dimensional (3D) topologies of dyadic clefts and associated membrane organelles were mapped in mouse ventricular myocardium using electron tomography. The morphological details and the distribution of membrane systems, including transverse tubules (T-tubules), junctional sarcoplasmic reticulum (SR) and vicinal mitochondria, were determined and presumed to be crucial for controlling cardiac Ca(2+) dynamics. The geometric complexity of T-tubules that varied in diameter with frequent branching was clarified. Dyadic clefts were intricately shaped and remarkably small (average 4.39x10(5) nm(3), median 2.81x10(5) nm(3)). Although a dyadic cleft of average size could hold maximum 43 ryanodine receptor (RyR) tetramers, more than one-third of clefts were smaller than the size that is able to package as many as 15 RyR tetramers. The dyadic clefts were also adjacent to one another (average end-to-end distance to the nearest dyadic cleft, 19.9 nm) and were distributed irregularly along T-tubule branches. Electron-dense structures that linked membrane organelles were frequently observed between mitochondrial outer membranes and SR or T-tubules. We, thus, propose that the topology of dyadic clefts and the neighboring cellular micro-architecture are the major determinants of the local control of Ca(2+) in the heart, including the establishment of the quantal nature of SR Ca(2+) releases (e.g. Ca(2+) sparks).

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.

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.

Tetrameric mouse acetylcholinesterase: continuum diffusion rate calculations by solving the steady-state Smoluchowski equation using finite element methods.

The tetramer is the most important form for acetylcholinesterase in physiological conditions, i.e., in the neuromuscular junction and the nervous system. It is important to study the diffusion of acetylcholine to the active sites of the tetrameric enzyme to understand the overall signal transduction process in these cellular components. Crystallographic studies revealed two different forms of tetramers, suggesting a flexible tetramer model for acetylcholinesterase. Using a recently developed finite element solver for the steady-state Smoluchowski equation, we have calculated the reaction rate for three mouse acetylcholinesterase tetramers using these two crystal structures and an intermediate structure as templates. Our results show that the reaction rates differ for different individual active sites in the compact tetramer crystal structure, and the rates are similar for different individual active sites in the other crystal structure and the intermediate structure. In the limit of zero salt, the reaction rates per active site for the tetramers are the same as that for the monomer, whereas at higher ionic strength, the rates per active site for the tetramers are approximately 67%-75% of the rate for the monomer. By analyzing the effect of electrostatic forces on ACh diffusion, we find that electrostatic forces play an even more important role for the tetramers than for the monomer. This study also shows that the finite element solver is well suited for solving the diffusion problem within complicated geometries.