Reduced basis method for the nonlinear Poisson–Boltzmann equation regularized by the range-separated canonical tensor format

The Poisson–Boltzmann equation (PBE) is a nonlinear elliptic parametrized partial differential equation that arises in biomolecular modeling and is a fundamental tool for structural biology. It is used to calculate electrostatic potentials around an ensemble of fixed charges immersed in an ionic solution. Efficient numerical computation of the PBE yields a high number of degrees of freedom in the resultant algebraic system of equations, ranging from several hundred thousand to millions. Coupled with the fact that in most cases the PBE requires to be solved multiple times for a large number of system configurations, for example, in Brownian dynamics simulations or in the computation of similarity indices for protein interaction analysis, this poses great computational challenges to conventional numerical techniques. To accelerate such onerous computations, we suggest to apply the reduced basis method (RBM) and the (discrete) empirical interpolation method ((D)EIM) to the PBE with a special focus on simulations of complex biomolecular systems, which greatly reduces this computational complexity by constructing a reduced order model (ROM) of typically low dimension. In this study, we employ a simple version of the PBE for proof of concept and discretize the linearized PBE (LPBE) with a centered finite difference scheme. The resultant linear system is solved by the aggregation-based algebraic multigrid (AGMG) method at different samples of ionic strength on a three-dimensional Cartesian grid. The discretized LPBE, which we call the high-fidelity full order model (FOM), yields solution as accurate as other LPBE solvers. We then apply the RBM to the FOM. DEIM is applied to the Dirichlet boundary conditions which are nonaffine in the parameter (ionic strength), to reduce the complexity of the ROM. From the numerical results, we notice that the RBM reduces the model order from {mathcal {N}} = 2times 10^{6} to N = 6 at an accuracy of {mathcal {O}}(10^{-9}) and reduces the runtime by a factor of approximately 7600. DEIM, on the other hand, is also used in the offline-online phase of solving the ROM for different values of parameters which provides a speed-up of 20 for a single iteration of the greedy algorithm.

Conventional reduced order modeling techniques such as the reduced basis (RB) method (relying, e.g., on proper orthogonal decomposition (POD)) may incur in severe limitations when dealing with nonlinear time-dependent parametrized PDEs, as these are strongly anchored to the assumption of modal linear superimposition they are based on. For problems featuring coherent structures that propagate over time such as transport, wave, or convection-dominated phenomena, the RB method may yield inefficient reduced order models (ROMs) when very high levels of accuracy are required. To overcome this limitation, in this work, we propose a new nonlinear approach to set ROMs by exploiting deep learning (DL) algorithms. In the resulting nonlinear ROM, which we refer to as DL-ROM, both the nonlinear trial manifold (corresponding to the set of basis functions in a linear ROM) as well as the nonlinear reduced dynamics (corresponding to the projection stage in a linear ROM) are learned in a non-intrusive way by relying on DL algorithms; the latter are trained on a set of full order model (FOM) solutions obtained for different parameter values. We show how to construct a DL-ROM for both linear and nonlinear time-dependent parametrized PDEs. Moreover, we assess its accuracy and efficiency on different parametrized PDE problems. Numerical results indicate that DL-ROMs whose dimension is equal to the intrinsic dimensionality of the PDE solutions manifold are able to efficiently approximate the solution of parametrized PDEs, especially in cases for which a huge number of POD modes would have been necessary to achieve the same degree of accuracy.

The present study is concerned with unsteady electroosmotic flow (EOF) in a microchannel with the electric charge distribution described by the Poisson-Boltzmann (PB) equation. The nonlinear PB equation is solved by a systematic perturbation with respect to the parameter λ which measures the strength of the wall zeta potential relative to the thermal potential. In the small λ limits (λ<<1), we recover the linearized PB equation - the Debye-Hückel approximation. The solutions obtained by using only three terms in the perturbation series are shown to be accurate with errors <1% for λ up to 2. The accurate solution to the PB equation is then used to solve the electrokinetic fluid transport equation for two types of unsteady flow: transient flow driven by a suddenly applied voltage and oscillatory flow driven by a time-harmonic voltage. The solution for the transient flow has important implications on EOF as an effective means for transporting electrolytes in microchannels with various electrokinetic widths. On the other hand, the solution for the oscillatory flow is shown to have important physical implications on EOF in mixing electrolytes in terms of the amplitude and phase of the resulting time-harmonic EOF rate, which depends on the applied frequency and the electrokinetic width of the microchannel as well as on the parameter λ.

An adaptive Cartesian grid (ACG) concept is presented for the fast and robust numerical solution of the 3D Poisson-Boltzmann Equation (PBE) governing the electrostatic interactions of large-scale biomolecules and highly charged multi-biomolecular assemblies such as ribosomes and viruses. The ACG offers numerous advantages over competing grid topologies such as regular 3D lattices and unstructured grids. For very large biological molecules and multi-biomolecule assemblies, the total number of grid-points is several orders of magnitude less than that required in a conventional lattice grid used in the current PBE solvers thus allowing the end user to obtain accurate and stable nonlinear PBE solutions on a desktop computer. Compared to tetrahedral-based unstructured grids, ACG offers a simpler hierarchical grid structure, which is naturally suited to multigrid, relieves indirect addressing requirements and uses fewer neighboring nodes in the finite difference stencils. Construction of the ACG and determination of the dielectric/ionic maps are straightforward, fast and require minimal user intervention. Charge singularities are eliminated by reformulating the problem to produce the reaction field potential in the molecular interior and the total electrostatic potential in the exterior ionic solvent region. This approach minimizes grid-dependency and alleviates the need for fine grid spacing near atomic charge sites. The technical portion of this paper contains three parts. First, the ACG and its construction for general biomolecular geometries are described. Next, a discrete approximation to the PBE upon this mesh is derived. Finally, the overall solution procedure and multigrid implementation are summarized. Results obtained with the ACG-based PBE solver are presented for: (i) a low dielectric spherical cavity, containing interior point charges, embedded in a high dielectric ionic solvent - analytical solutions are available for this case, thus allowing rigorous assessment of the solution accuracy; (ii) a pair of low dielectric charged spheres embedded in a ionic solvent to compute electrostatic interaction free energies as a function of the distance between sphere centers; (iii) surface potentials of proteins, nucleic acids and their larger-scale assemblies such as ribosomes; and (iv) electrostatic solvation free energies and their salt sensitivities - obtained with both linear and nonlinear Poisson-Boltzmann equation - for a large set of proteins. These latter results along with timings can serve as benchmarks for comparing the performance of different PBE solvers.

New solution decomposition and minimization schemes for Poisson–Boltzmann equation in calculation of biomolecular electrostatics

The Poisson-Boltzmann (PB) equation describes the electrostatic potential of a biomolecular system composed by a molecule in a solvent. The electrostatic potential is involved in biomolecular models which are used in molecular simulation. In consequence, finding an efficient method to simulate the numerical solution of PB equation is very useful. In this paper, we establish a probabilistic interpretation of the nonlinear PB equation with Backward Stochastic Differential Equations (BSDEs). This interpretation requires an adaptation of existing results on BSDEs.

Reduced basis methods are projection-based model order reduction techniques for reducing the computational complexity of solving parametrized partial differential equation problems. In this work we discuss the design of pyMOR, a freely available software library of model order reduction algorithms, in particular reduced basis methods, implemented with the Python programming language. As its main design feature, all reduction algorithms in pyMOR are implemented generically via operations on well-defined vector array, operator and discretization interface classes. This allows for an easy integration with existing open-source high-performance partial differential equation solvers without adding any model reduction specific code to these solvers. Besides an in-depth discussion of pyMOR's design philosophy and architecture, we present several benchmark results and numerical examples showing the feasibility of our approach.

In this paper, we present the a posteriori error analysis for the reduced basis method (RBM) applied to nonlinear variational problems that depend on a parameter in a nonaffine manner. To this end, we generalize the analysis by Veroy and Patera [Int. J. Numer. Methods Fluids, 47 (2005), pp. 773-788] to nonaffine parametrized partial differential equations. We use the empirical interpolation method (EIM) in order to approximate the nonaffine parameter dependencies by a linear combination of affine functions. We also investigate a standard dual problem formulation, in particular for the computation of a general output functional, also in combination with the EIM. First, we study the well-posedness in terms of the Brezzi-Rappaz-Raviart theory. Then, we develop a posteriori error estimates and investigate offline/online decompositions. The a posteriori error analysis allows us to introduce an adaptive sampling procedure for the choice of the modes. Numerical experiments for a convection-diffusion problem around a rotating propeller show the effectivity of the scheme.

Electrostatic Contributions to the Binding Free Energy of the λcI Repressor to DNA

Experimentally observed sequence-selective binding of metal ion to DNA oligonucleotides have been compared with variations of electrostatic potential (EP) along the helix. Calculations of EP have been performed for three atomic models of the oligonucleotide duplex [d(CGCGAATTCGCG)2] using several variants of EP calculations, including a solution of non-linear Poisson-Boltzmann equation (NPBE). N7 atom of guanine adjacent to adenine base was identified as a region with the most negative electrostatic potential in the major groove. The EP value for the Me ion binding site surpasses the value for N7 of other guanines by 10–26% depending on particular duplex conformation. Qualitatively, the sequence dependent variations of EP near guanine N7 atoms are in agreement with the sequence-selective behavior of Mn(II) and Zn(II) ions as revealed by NMR experiments. But the difference in EP between the two most negative regions near guanine N7 atoms does not exceed 1.25 kT/e. Simple model suggests that metal ions are capable to form ion-hydrate complexes with G-Pu steps of DNA duplex. These complexes are formed via one Me…G and five Me…water coordination bonds with water molecules hydrogen bonded to two adjacent purine bases in the same chain. We suppose that such a stereospecific structural possibility is the main factor which control the sequence-selectivity in the metal ion binding. A combination of both mechanisms allows to explain sequence specific Mn(II) and Zn(II) binding to a set of oligonucleotides.

To model soft-landing of peptide ions on surfaces, it is important to have accurate intermolecular potentials between these ions and surfaces. As part of this goal, ab initio calculations at the MP2/aug-cc-pVTZ level of theory, with basis set superposition error (BSSE) corrections, were performed to determine both the long-range attractive and short-range repulsive potentials for CH(4) interacting with the -NH(3)(+) group of CH(3)NH(3)(+). Potential energy curves for four different orientations between CH(4) and CH(3)NH(3)(+) were determined from the calculations to obtain accurate descriptions of the interactions between the atoms of CH(4) and those of -NH(3)(+). A universal analytic function was not found that could accurately represent both the long-range and short-range potentials for collision energies as high as those obtained in surface-induced-dissociation (SID) experiments. Instead, long-range and short-range analytic potentials were developed separately, by simultaneously fitting the four ab initio potential energy curves with a sum of two-body interactions between the atoms of CH(4) and -NH(3)(+), and then connecting these long-range and short-range two-body potentials with switching functions. Following a previous work [J. Am. Chem. Soc., 2002, 124, 1524], these two-body potentials may be used to describe the interactions of the N and H atoms of the -NH(3)(+) group of a protonated peptide ion with the H and C atoms of alkane-type surfaces such as alkyl thiol self-assembled monolayers and H-terminated diamond. Accurate short-range and long-range potentials are imperative to model protonated peptide ion soft-landing experiments. The former controls the collision energy transfer, whereas the latter describes the binding of the ion to the surface. A comparison of the ab initio potential energy curves for CH(3)NH(3)(+)/CH(4) with those for NH(4)(+)/CH(4) shows that they give nearly identical two-body interactions between the atoms of -NH(3)(+) and those of CH(4), showing that the smaller NH(4)(+)/CH(4) system may be used to obtain the two-body potentials. A comparison of the four ab initio potential energy curves reported here for CH(3)NH(3)(+)/CH(4), with those given by the AMBER and CHARMM molecular mechanical potentials, show that these latter potentials "roughly" approximate the long-range attractions, but are grossly in error for the short-range repulsions. The work reported here illustrates that high-level ab initio calculations of intermolecular potentials between small model molecules may be used to develop accurate analytical intermolecular potentials between peptide ions and surfaces.

The use of reduced‐order models (ROMs) in physics‐based modeling and simulation almost always involves the use of linear reduced basis (RB) methods such as the proper orthogonal decomposition (POD). For some nonlinear problems, linear RB methods perform poorly, failing to provide an efficient subspace for the solution space. The use of nonlinear manifolds for ROMs has gained traction in recent years, showing increased performance for certain nonlinear problems over linear methods. Deep learning has been popular to this end through the use of autoencoders for providing a nonlinear trial manifold for the solution space. In this work, we present a non‐intrusive ROM framework for steady‐state parameterized partial differential equations that uses convolutional autoencoders to provide a nonlinear solution manifold and is augmented by Gaussian process regression (GPR) to approximate the expansion coefficients of the reduced model. When applied to a numerical example involving the steady incompressible Navier–Stokes equations solving a lid‐driven cavity problem, it is shown that the proposed ROM offers greater performance in prediction of full‐order states when compared to a popular method employing POD and GPR over a number of ROM dimensions.

A microfluidic pumping flow model driven by electro-osmosis mechanisms is developed to analyze the flow characteristics of aqueous electrolytes. The pumping model is designed based on a single propagative rhythmic membrane contraction applied on the upper wall of a microchannel. The flow lubrication theory coupled with a nonlinear Poisson–Boltzmann equation is used to model the microchannel unsteady creeping flow and to describe the distribution of the electric potential across the electric double layer. A generic solution is obtained for the Poisson–Boltzmann equation without the Debye–Hückel linearization. The effects of zeta potential, Debye length, and electric field on the potential distribution, pressure distribution, velocity profiles, shear stress, and net flow rate are computed and interpreted in detail. The results have shown that this electrokinetic membrane pumping model can be used to understand microlevel transport phenomena in various physiological systems. The proposed model can also be integrated with other microfluidic devices for moving microvolume of liquids in artificial capillaries used in modern biomedical applications.

We apply the reduced basis method to solve Navier–Stokes equations in parametrized domains. Special attention is devoted to the treatment of the parametrized nonlinear transport term in the reduced basis framework, including the case of nonaffine parametric dependence that is treated by an empirical interpolation method. This method features (i) a rapid global convergence owing to the property of the Galerkin projection onto a space WN spanned by solutions of the governing partial differential equation at N (optimally) selected points in the parameter space, and (ii) the offline/online computational procedures that decouple the generation and projection stages of the approximation process.This method is well suited for the repeated and rapid evaluations required in the context of parameter estimation, design, optimization, and real‐time control. Our analysis focuses on: (i) the pressure treatment of incompressible Navier–Stokes problem; (ii) the fulfillment of an equivalent inf–sup condition to guarantee the stability of the reduced basis solutions. The applications that we consider involve parametrized geometries, like e.g. a channel with curved upper wall or an arterial bypass configuration. © 2007 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 23: 923–948, 2007

A variational solution of the nonlinear poisson-boltzmann equation inside a spherical cavity